Thorium in solar twins: implications for habitability in rocky planets

ABSTRACT

We have investigated the thorium (Th) abundance in a sample of 53 thin disc solar twins covering a wide range of ages. These data provide constrains on the mantle energy budget of terrestrial planets that can be formed over the evolution of the Galaxy’s thin disc. We have estimated Th abundances with an average precision of 0.025 dex (in both [Th/H] and [Th/Fe]) through comprehensive spectral synthesis of a Th ii line present at 4019.1290 Å, using very high resolution (R  =  115 000) high quality HARPS spectra obtained at the ESO La Silla Observatory. We have confirmed that there is a large energy budget from Th decay for maintaining mantle convection inside potential rocky planets around solar twins, from the Galactic thin disc formation until now, because the pristine [Th/H]_ZAMS is super-solar on average under a uniform dispersion of 0.056 dex (varying from |$+$|0.037 up to |$+$|0.138 dex based on linear fits against isochrone stellar age). Comparing to neodymium (Nd) and europium (Eu), two others neutron-capture elements, the stellar pristine abundance of Th follows Eu along the Galactic thin disc evolution, but it does not follow Nd, probably because neodymium has a significant contribution from the s-process (about 60 per cent).

stars: abundances, stars: fundamental parameters, stars: solar-type, Galaxy: disc, solar neighbourhood; planets and satellites: terrestrial planets

1 INTRODUCTION

The structure of a terrestrial planet is driven by the formation process itself based on the collision and accretion of planetesimals with differential mineralogical settling and also by the internal heat budget. As the prototype geologically dynamic planet, Earth has a partially crystallized metallic dynamo core, and a silicate liquid mantle, whose convection maintains plate tectonics of a rocky crust. The major internal heat sources of Earth come from the secular cooling of core and mantle, and the radioactive decay of Th, U and K isotopes in the mantle (Huang et al. 2013). The mantle thickness and convection inside telluric planets are directly linked to the abundances of these isotopes in the mantle (basically ²³²Th, ²³⁵U and ²³⁸U) (McDonough 2003). The volcanic activity generated by the plate tectonics recycles gases such as carbon dioxide between the atmosphere and the mantle, contributing to habitability by keeping the planet surface at moderate temperatures (Walker, Hays & Kasting 1981). A geologically active planet could even play a role for the origin of life and its long-term maintenance (e.g. Misra et al. 2015, and references therein).

Unterborn, Johnson & Panero (2015) were the first to speculate about the potential energy budget of terrestrial planets directly connected to the stellar photospheric abundance of thorium (Th). They measured log ε(Th)¹ varying from 59 up to 251 per cent of the solar value in a sample of 13 solar twins and analogues, indicating that possible terrestrial planets formed around these stars have sufficient internal energy budget to sustain a convective mantle in a dynamic planet and, consequently, increases the probability of having a habitable surface. They also found that the Sun is depleted in Th in comparison with 13 solar twins/analogues. They argued that a large variation in Th abundance could be explained by differences in stellar age within their sample, besides other main possible factors such as the continuum normalization of the analysed blend Fe-Ni-Mn-Th-Co-CN-Ce-CH at 4019 Å (covering λλ4018.9–4019.2 Å) and the difference in spectral resolution between two data sets (R  =  48 000 and R  =  115 000). Specifically, they did not report any deep analysis of the Th abundance as a function of stellar age (taken from Baumann et al. 2010). They concluded that the observed large Th abundance variation is likely due to primordial chemical inhomogeneity among the stars.

Thorium is an r-process element synthesized in the Galaxy through several astrophysical processes such as core-collapse supernovae, neutron star mergers, and/or neutron star – black hole mergers (r-process elements are those built up by rapid capture of neutrons into nuclei, i.e. over energetic processes in which the neutron flux is high enough to become negligible the effect of the free neutrons decay). The half time of its most abundant radioactive isotope ²³²Th (14.05 Gyr) is as high as the age of the Universe (Cayrel et al. 2001). Previous works have found a large variation in the abundance of r-process elements (Th included) covering both metal-poor old halo stars and the younger more metal-rich disc stars (Honda et al. 2004; del Peloso, da Silva & Arany-Prado 2005a; del Peloso et al. 2005b; Johnson & Bolte 2001). Butcher (1987) was the pioneer in estimating the age of the Galaxy by measuring Th abundance in G-dwarfs, taking the advantage of its sizeable decay during the lifetime of the Galaxy. How the abundance of Th has evolved during the Galactic disc’s life is still an open question – more specifically in the thin disc at the Sun’s Galactocentric distance. Nearby solar twins spanning a range of ages can provide clues about this question because their photospheric abundances reflect the composition of the interstellar medium at the time they formed. As planets and stars are likely formed from the same natal cloud, potential terrestrial planets should have a Th abundance related to the Th abundance of its host star.

We have homogeneously measured Th abundances in a sample of nearby solar twins (including two planet-host stars) of the Galactic thin disc, covering a broad range of stellar ages, in order to provide constrains about the mantle energy budget of terrestrial planets during the evolution of the thin disc. A Th ii line present in a multispecies blend at 4019 Å has been carefully analysed for this purpose, through comprehensive spectral synthesis differentially to the Sun. Section 2 presents the sample of solar twins, their main parameters, and the spectroscopic data. In Section 3, we describe the spectral synthesis approach employed for measuring Th abundance in a multispecies blend. Section 4 covers the analysis of resulting Th abundance ratios (current and primordial) against metallicity and stellar age. Finally, we show the concluding remarks in Section 5.

2 SAMPLE OF SOLAR TWINS AND OBSERVATIONAL DATA

The sample is composed of 67 solar twins (stars with effective temperature, surface gravity, and metallicity around the solar values within ± 100 K in T_eff and within ± 0.1 dex in log g and [Fe/H]), which were recently analysed by Spina et al. (2018) and Bedell et al. (2018). Their spectroscopic data and main parameters adopted are described in this section.

Spina et al. (2018) derived their photospheric parameters (and abundance of 12 neutron-capture, n-capture, elements) by applying a line-by-line differential spectroscopic analysis relative to the Sun through equivalent width (EW) measurements of Fe i and Fe ii lines. The estimated typical errors in T_eff, log g, [Fe/H], and ξ (micro-turbulence velocity) are, respectively, 4 K, 0.012, 0.004 dex and 0.011 km s⁻¹. We have also adopted the stellar masses m and isochrone ages t as derived by Spina et al. (2018), which are represented by the more probable values (i.e. that corresponds to the peak of a probability distribution) by applying the isochrone method with the q² code of Ramírez, Meléndez & Asplund (2014a) and Ramírez et al. (2014b). Only for the two very young stars HIP 3203 and HIP 4909, whose more probable age values are null, the ages are better represented by the average values, as also adopted by our previous works (Spina et al. 2018). The typical age uncertainty is 0.4 Gyr, which is derived from the probability distribution itself. Besides the 12 n-capture elements (Sr, Y, Zr, Ba, La, Ce, Pr, Nd, Sm, Eu, Gd, and Dy), we have also considered the abundances of 17 additional elements analysed by Bedell et al. (2018) (C, O, Na, Mg, Al, Si, S, Ca, Sc, Ti, V, Cr, Mn, Co, Ni, Cu and Zn). The aforementioned works measured elemental abundances by a strict line-by-line EW analysis relative to the Sun.

High Accuracy Radial velocity Planet Searcher (HARPS) spectra are used in this work to extract Th abundances. HARPS is an ultra-stable echelle spectrograph installed on the 3.6m telescope of the European Southern Observatory (ESO) at La Silla Observatory in Chile (Mayor et al. 2003). The HARPS spectra covers λλ3780–6910 Å under a resolving power R  =  115 000. The spectra adopted in the current work are those analysed by Spina et al. (2018) and Bedell et al. (2018). Each one-dimensional spectrum comes from more than 50 stacked spectra that are previously Doppler-corrected and carefully continuum normalized. The majority of spectra were collected for a large ESO observing program on HARPS, which aims at hunting for planets around solar twins (Bedell et al. 2015; Meléndez et al. 2017, 2015), complemented with HARPS spectra from other programs extracted from the online ESO Science Archive Facility. The spectral signal-to-noise ratio (SNR), measured at 6000 Å, has an average value of about 800 per pixel, varying from 300 up to 1800 per pixel. The solar reference spectrum is a combined spectrum of the reflected light by the asteroid Vesta. Its flux is continuum normalized in the same way as the spectra of the solar twins, and it has SNR  =  1300 per pixel at 6000 Å. Further details about the analysed HARPS spectra are given in Spina et al. (2018). Table 1 shows the main parameters of sample solar twins, covering 58 stars as described at the end of the next section.

Table 1.

Stellar parameters of the sample stars collected from previous published works (covering 58 solar twins): photospheric parameters and isochrone age (Spina et al. 2018), macro-turbulence and rotation velocities (dos Santos et al. 2016), and elemental abundances (Bedell et al. 2018). Full table online.

Star ID	Teff	Log g	[Fe/H]	ξ	Vmacro	V.sin(i)	Age	[Si/H]	[Nd/H]	[Eu/H]
	(K)		(dex)	(km s−1)	(km s−1)	(km s−1)	(Gyr)	(dex)	(dex)	(dex)
Sun	5777	4.440	0.000	1.00	3.20	2.04	4.56	0.000	0.000	0.000
HIP 003203	5868 ± 9	4.540 ± 0.016	|$-$|0.050 ± 0.007	1.16 ± 0.02	3.27	3.82	0.50 ± 0.30	|$-$|0.112 ± 0.009	0.160 ± 0.017	0.057 ± 0.014
HIP 004909	5861 ± 7	4.500 ± 0.016	0.048 ± 0.006	1.11 ± 0.01	3.33	4.01	0.60 ± 0.40	|$-$|0.047 ± 0.012	0.219 ± 0.012	0.099 ± 0.013
HIP 007585	5822 ± 3	4.445 ± 0.008	0.083 ± 0.003	1.01 ± 0.01	3.37	1.90	3.50 ± 0.40	0.058 ± 0.008	0.138 ± 0.009	0.124 ± 0.015
HIP 008507	5717 ± 3	4.460 ± 0.011	|$-$|0.099 ± 0.003	0.96 ± 0.01	2.88	0.77	4.90 ± 0.45	|$-$|0.127 ± 0.028	0.044 ± 0.012	|$-$|0.019 ± 0.022
–	–	–	–	–	–	–	–	–	–	–
HIP 118115	5798 ± 4	4.275 ± 0.011	|$-$|0.036 ± 0.003	1.10 ± 0.01	3.55	0.89	8.00 ± 0.30	|$-$|0.029 ± 0.012	0.023 ± 0.006	0.069 ± 0.007

3 HOMOGENEOUS THORIUM ABUNDANCES

The abundance of Th has been measured from a complex blend at 4019 Å, which has been quite employed for this purpose, both in the study of Galactic disc stars (e.g. Morell, Kallander & Butcher 1992; del Peloso et al. 2005b; Unterborn et al. 2015) and Galactic halo stars (e.g. Cowan et al. 1997; Johnson & Bolte 2001; Roederer et al. 2009). The analysis of this multispecies blend requires a careful spectral synthesis at very high resolution and high SNR. The Th ii 4019.1290-Å line together with the Co i 4019.1260-Å, CH 4019.1390-Å, and CN 4019.2060-Å lines make the blend asymmetric at the red wing. The very weak CH line is from the (3,2) vibrational band of the A-X electronic system (¹²C¹H), and the relatively weak CN line is from the (4,5) vibrational band of the B-X electronic system (¹²C¹⁴N). In the case of solar-type stars, the major atomic contributor of this blend is the Fe i 4019.0420-Å line followed by the Ni i 4019.0667-Å and Mn i 4019.0658-Å lines. A minor contribution comes from the Ce ii 4019.0569-Å line. We refer the reader to Fig. 1 to see the different contributions to the blend, based on our calibration through spectral synthesis of the solar spectrum, as detailed in this section.

Figure 1.

Spectral synthesis calibration of the Fe-Ni-Mn-Th-Co-CN-Ce-CH blend at 4019 Å to the solar spectrum after the procedure of gf-values adjustment as explained in Section 3: multispecies atomic and molecular line contributions at the left top panel, blended global profile reproduction at the right top panel, comparisons between the observed spectrum and seven synthetic spectra at the left bottom panel, and χ² minimization graph at the right bottom panel. The Th ii 4019.1290-Å line is present in the red wing.

There is no measurable line of uranium in the HARPS spectral coverage. However, there is another Th ii line in the HARPS range (at 4086.5210 Å), but it is 4.5 times weaker than the Th ii 4019.1290-Å line. Uranium is another r-process radioactive element that also contributes for the energy budget of terrestrial planets (mantle convection/thickness).

The spectral synthesis has been performed in the region λλ4000–4030 Å using the 2014 version of the moog code (Sneden 1973) that solves the photospheric radiative transfer with line formation in FGK stars under the LTE (local thermodynamic equilibrium) approximation. The spectral synthesis based on moog has been handled through our own automated procedure developed under the python language (version 2.7), which covers several steps from the solar gf calibration to the derivation of the Th abundance with its uncertainty, as explained in this section. We have also adopted the same stellar parameters and model atmospheres as our previous analyses of the sample stars (Bedell et al. 2018; Spina et al. 2018) to be self-consistent with them, i.e. interpolated models in the grid of Castelli & Kurucz (2004) (ATLAS9) assuming solar abundance ratios. Throughout our differential analysis, we have also adopted the standard parameters for the Sun: T_eff  =  5777 K, log g = 4.44, [Fe/H] = 0.00 dex, and ξ  =  1.00 km s⁻¹ (e.g. Cox 2000). The solar chemical composition by Asplund et al. (2009), and by Grevesse et al. (2015) specifically for the heavy elements Cu to Th, have been adopted as reference, such that the adopted Sun’s Th abundance is logε(Th) = 0.03 dex. The wavelength step of synthetic spectra has been fixed to 0.01 Å.

The atomic line list has been compiled from the VALD database (Ryabchikova et al. 2015), and we have added the following molecular lines from the Kurucz database (Kurucz 2017): ¹²C¹H and ¹³C¹H lines of the A-X system, ¹²C¹H and ¹³C¹H lines of the B-X system, ¹²C¹²C lines of the d-a system, ¹²C¹⁴N lines of the B-X system, and ¹⁴N¹H lines of the A-X system. In the case of solar-like stars, only lines of the A-X system of ¹²CH, B-X system of ¹²CH, d-a system of ¹²C¹²C and B-X system of ¹²C¹⁴N, actually produce detectable absorptions in the studied region at our SNR (∼ 800). The other molecular lines (the main isotopes of carbon and nitrogen) are included in the line list to make more realistic the reproduction of the stellar continuum, because they are supposed to be within the spectral noise for SNR ≤ 1000 (¹³CH B-X system), or even 10 000 (¹³CH A-X system, ¹⁴NH A-X system). We have fixed solar isotopic ratios for ¹²C/¹³C and ¹⁴N/¹⁵N that are, respectively, 89.4 and 435 (Asplund et al. 2009). We have also confirmed that the lines of CH A-X and CN B-X in the investigated 30 Å-wide region are all insensitive to variation of C and N isotopic ratios. There are very few lines of CH B-X that are somehow sensitive to the ¹²C/¹³C ratio, but none affects the analysed Fe-Ni-Mn-Th-Co-CN-Ce-CH blended feature at λλ4018.9–4019.2 Å. The adopted dissociative energies of CH, C₂, CN, and NH are those by Barklem & Collet (2006). We tested somewhat different values of the dissociation constants, and found no changes in the results within about 0.01 dex. Table 2 presents the list of atomic and molecular lines covering 2 Å around the Th ii feature, whose gf values have been calibrated to the cited solar spectrum, as explained below.

Table 2.

Line list after the gf calibration to the solar spectrum: wavelength, species code, excitation potential of transition inferior level, log gf and species identification. The list covers 2 Å around the Fe-Ni-Mn-Th-Co-CN-Ce-CH blend at 4019 Å. The species code is the moog standard notation, i.e. atomic number(s) before the decimal point (listed in crescent order for molecules) followed by the ionization level immediately after the decimal point (0: neutral, 1: first ionized, and so on) and also by mass numbers for molecules in the case of isotopic species discriminated (listed in crescent order). Full table online.

wavelength	species code	χe	gf	species
(Å)		(eV)
4018.5080	106.00113	1.391	0.276E-05	CH
4018.5090	607.0	3.754	0.154E |$+$| 00	CN
4018.5360	607.0	3.193	0.424E-04	CN
4018.5620	107.0	2.440	0.433E-02	NH
4018.5703	25.0	5.087	0.195E |$+$| 00	Mn i
–	–	–	–	–
4020.4940	106.00113	0.595	0.162E-04	CH

The first step for the differential chemical analysis relatively to the Sun using spectral synthesis is the calibration of gf-values for reproducing the line profiles, adopting a solar spectrum obtained with the same instrument and resolution. We have performed the calibration of gf-values of all lines in the region λλ4000–4030 Å through an automatic procedure under a line-by-line basis, using the flux at the core of each individual line. Specifically, for a very narrow region around the 4019-Å blend (about 1Å wide), we have refined the gf calibration focusing on the contribution lines of the blend itself paying a special attention to the continuum adjustment between the observed and synthetic spectra. Due to the very close proximity between the cores of the Th ii and Co i lines under the HARPS spectral sampling (4019.1290 and 4019.1260 Å, respectively), we have fixed the gf-values for both lines, adopting accurate known experimental values (0.592 ± 0.018 for the Th ii line by Nilsson et al. 2002 and 0.537 ± 0.049E-02 for the Co i line by Lawler, Whaling & Grevesse 1990). Whilst the gf-value of the Co i line has been actually fixed, the gf-value of the Th ii line has been slightly refined to reproduce the spectral blend against the solar spectrum, ensuring a confident differential analysis to the Sun (the gf fine tuning of the Th ii line has been just 1.9 per cent greater than the experimental measurement that has a relative error of 3 per cent). We have got [Th/H] = [Th/Fe]  =  0.005 dex over this calibration step relatively to the Sun. Fig. 1 illustrates the spectral synthesis calibration of the Fe-Ni-Mn-Th-Co-CN-Ce-CH blend at 4019 Å to the solar spectrum. We have been also able to confirm the relative intensity between the Th ii and Co i lines as semitheoretically predicted by Grevesse et al. (2015). The ratio EW(Th ii 4019.1290-Å line)/EW(Co i 4019.1260-Å line) agrees within 1σ deviation taking into account the 20 per cent estimated uncertainty for the EW(Th ii 4019.1290-Å line) made by them. To finely reproduce the line broadening, we have adopted the line-of-sight rotational velocity V.sin(i) and macro-turbulence velocity V_macro that were, respectively, measured and computed by dos Santos et al. (2016) for our sample solar twins (the later one parametrically as a function of T_eff and log g). A limb darkening linear law for solar-type stars has been employed (limb darkening parameter u  =  0.60 adopted, to be consistent with dos Santos et al. 2016). Based on our calibration to the HARPS solar spectrum we estimate that the CH 4019.1390-Å and CN 4019.2060-Å lines are 81 and 28 per cent weaker than the Th ii line, respectively. We have computed their EW in each individual synthetic spectrum. Note that the cores of CH and CN lines are, respectively, 0.010 and 0.077 Å far away from the core of Th ii line (see Fig. 1).

An automatic procedure based on the χ² minimization has been used for performing the spectral synthesis fit in order to extract the Th abundance from the blended feature (χ²  =  |$\sum _{i{\, =\, }1,5}$| (O_i − S_i)²/σ(O_i)², where O_i is the flux of the observed spectrum, S_i is the flux of the synthetic spectrum, σ(O_i) is the error in the observed flux, and i represents the wavelength point). The wavelength window to compute χ² is centred at the core of Th ii line, as the spectral sampling (i.e. 4019.13 Å), and it covers a region in which the line absorption flux is stronger than or comparable to the other lines taking the solar spectrum as reference (see Fig. 1). The χ² window is λλ4019.11–4019.15 Å, ranging over 5 pixels. The error in observed flux is estimated as a function of the continuum SNR, such that σ(O_i)  = O_i/SNR_continuum. The continuum SNR in the blend region has been measured as a simple average over five small continuum windows (centred at 4002.24 Å, 4009.35 Å, 4019.67 Å, 4020.67 Å, and 4021.47 Å). The local SNR ranges from 147 (HIP 8507) up to 993 (HIP 54582), having an average value of 379 (σ(SNR)  =  152). For the Sun’s spectrum, we have SNR  =  439. The normalization of flux scale is fine tuned for the absorption line blend by applying a multiplicative correction factor taking as reference the continuum points at 4002.24 Å and/or 4021.47 Å.

The χ² procedure is particularly useful to directly estimate the error in Th abundance due to the spectral synthesis fitting itself, such that the error in [Th/Fe] is derived when χ² increases in ν unities of degree of freedom from the minimum χ² value along the polynomial fitting curve χ² versus [Th/Fe]. In this case, the degree of freedom ν is equal to 4, because the χ² window covers 5 pixels and there is a single free parameter that is the Th abundance itself (i.e. ν  =  n_pixels − n_parameters). The thorium abundance uncertainties due to the propagation of the errors of photospheric parameters have been added in quadrature to the error of the spectral synthesis. Fig. 2 illustrates the spectral synthesis procedure of the Th ii line in the case of HIP 85042.

Figure 2.

HIP 85042: Three spectral synthesis diagnostic plots (global profile, continuum windows and spectral comparisons) and χ² minimization graph (right bottom panel) for the Fe-Ni-Mn-Th-Co-CN-Ce-CH blend at 4019 Å, whose red wing is sensitive to the Th abundance variation (Th ii 4019.1290 Å-line). Seven theoretical spectra are computed assuming a uniform step of 0.10 dex in [Th/Fe]. A polynomial fit is applied for χ² as a function of [Th/Fe] to search for the minimum value and the resulting Th abundance ratio.

We have also verified the influence of errors in the abundances of carbon and cobalt in the Th abundance error due to the close proximity of those lines of Co i, CH, and CN to the core of the Th ii line. Taking the case of HIP 85042 as representative for the sample stars, the uncertainties in [Co/H] and [C/H] (Bedell et al. 2018) make changes in [Th/Fe] correspondent, respectively, to 0.003 and 0.002 dex, which are many times smaller than the other involved errors. Therefore, the impact of Co and C abundance errors in the Th abundance analysis is negligible, and the final error in [Th/Fe] is basically fixed by the spectral synthesis error itself. The average error in [Th/Fe] is 0.025 dex considering 53 stars in total, after removing objects whose derived Th abundances are not acceptable, as justified in the following paragraph. Finally, the error in [Th/H] is computed in quadrature taking into account the error in [Fe/H], and it coincides with that one in [Th/Fe].

The final step for measuring the Th abundance is to check the quality of spectral synthesis of the whole analysed absorption blend on a star-by-star basis (sample of 67 stars). First of all, the fit is very unsatisfactory in the cases of HIP 65708 and HIP 83276 (the later certainly due to an unidentified problem in the data reduction). Results are not acceptable due to the low spectral SNR at the studied 30 Å-wide region for HIP 6407, HIP 81746, and HIP 89650 (roughly for SNR ≲ 130 and/or σ[Th/Fe] ≳ 0.070 dex). The blend has not been reproduced well for other four stars: HIP 10303, HIP 14501, HIP 30037, and HIP 115577 (two of them are alpha-rich old solar twins and there is nothing special related to the two others). Three out of these nine eliminated stars (HIP 14501, HIP 65708, and HIP 115577) were also excluded from the fits made in the chemical evolution study of Bedell et al. (2018) because they were supposed to belong to the thick disc (ages above 8 Gyr and a visible enhancement in α-elements).

For 15 cases (HIP 3203, HIP 7585, HIP 18844, HIP 22263, HIP 25670, HIP 28066, HIP 30476, HIP 34511, HIP 40133, HIP 42333, HIP 73241, HIP 74432, HIP 79672, HIP 101905, and HIP 104045), we have improved the automated spectral synthesis fit by slightly changing the continuum level. We have also revised the spectral synthesis fit of the Th ii blended feature for a few very young stars of our sample (four stars in total) by modifying the rotational velocity in order to better globally reproduce the feature (decrease by about 8–20 per cent in V.sin(i)). They are HIP 3203, HIP 4909, HIP 38072, and HIP 101905, having, respectively, a decrease of 8, 10, 20, and 13 per cent in V.sin(i). Despite of few cases, for which the rotational broadening had to be modified, we can state that the Th abundance is more sensitive to the continuum level than the rotational broadening.

We have investigated for several cases the impact of changing the abundance of Fe in order to get a better spectral synthesis fit since the blended feature is dominated by an Fe line, specially at the blue wing. We conclude that the derived Th abundance is just slightly modified within its error, even when the log ε(Fe) increases or decreases up to 5–10 times its error. The same conclusion is reached by changing the abundances of Ni, Mn, and/or Ce, whose lines shape the blended red wing. Consequently, their derived Th abundances are not updated upon the effects of changing the Fe, Ni, and Mn individual abundances.

In order to be consistent with the work of Bedell et al. (2018), we have decided to exclude from our analysis HIP 28066, HIP 30476, HIP 73241, and HIP 74432 because they are α-rich old solar twins, as well as HIP 64150 because it exhibits anomalous s-process element abundances likely due to a past contamination from a close binary companion ( dos Santos et al. 2017; Bedell et al. 2018; Spina et al. 2018).

In total, we have measured and make available Th abundance for 58 solar twins. However, just 53 thin disc solar twins are adopted in our analysis of correlations of the Th abundance as function of [Fe/H] and isochrone stellar age (one chemically anomalous in s-elements and four α-rich old stars are excluded). The Sun is not included in the analysis too. Table 3 shows the derived Th abundances in 58 solar twin stars.

Table 3.

Thorium abundance measured in this work relatively to the Sun for 58 solar twins (Sun included as well). Full table online.

Star ID	[Th/H]	Error	[Th/Fe]	Error	SNR	|$\chi ^{2}_{\rm min}/\nu$|
	(dex)	(dex)	(dex)	(dex)
Sun	|$+$|0.005	0.021	|$+$|0.005	0.021	439	6.144
HIP 003203	|$+$|0.120	0.034	|$+$|0.170	0.033	319	0.041
HIP 004909	|$+$|0.180	0.023	|$+$|0.132	0.022	534	10.978
HIP 006407	|$+$|0.004	0.081	|$+$|0.062	0.081	121	1.901
HIP 007585	|$+$|0.193	0.018	|$+$|0.110	0.018	454	3.796
–	–	–	–	–	–	–
HIP 118115	|$+$|0.009	0.026	|$+$|0.045	0.026	295	0.071

4 THORIUM ABUNDANCE RATIOS VERSUS METALLICITY AND STELLAR AGE

We have analysed the derived Th abundances for 53 thin disc solar twins as a function of both metallicity and isochrone age. In this section, we compile our results.

We have verified that the current (observed) Th abundance [Th/H] in solar-twin stars follows the Fe abundance over the very restricted metallicity interval of them. There is no (anti)correlation between the observed values of [Th/Fe] and [Fe/H]. Fig. 3 shows the observed [Th/H] and [Th/Fe] as a function of [Fe/H], in which linear fitting curves are drawn for illustration purposes only. The fitting procedure is further described in this section.

Figure 3.

[Th/H]_observed and [Th/Fe]_observed as a function of [Fe/H]: A statistically robust linear fit is shown in each plot (blue solid line). The fitting equation with coefficients and their errors inside parenthesis are also displayed. The fit corresponds to 53 thin disc solar twins only (black filled circles). Data of excluded stars are plotted together for illustration purposes only (red empty triangles: four α-rich old stars, and magenta empty square: the chemically anomalous in s-elements HIP 64150). Sun’s data are also plotted as reference at the coordinates origin (green solar standard symbol), although are consistent with the fits.

The observed strong anticorrelations [Th/H]–age and [Th/Fe]–age (current [Th/H] and [Th/Fe] ratios) are partially explained by the radioactive decay of Th, because after correcting the Th abundances by this effect, the anticorrelation [Th/H]–age keeps existing rather than being erased. Fig. 4 shows the observed [Th/H] and [Th/Fe] as a function of the isochrone stellar age, in which linear fitting curves are drawn for illustration purposes only.

Figure 4.

[Th/H]_observed and [Th/Fe]_observed as a function of the isochrone stellar age: A statistically robust linear fit is shown in each plot (blue solid line). The fitting equation with coefficients and their errors inside parenthesis are also displayed. The fits correspond to 53 thin disc solar twins only (black filled circles). Data of excluded stars are plotted together for illustration purposes only (red empty triangles: four α-rich old stars, and magenta empty square: the chemically anomalous in s-elements HIP 64150). Sun’s data are also plotted as reference (green solar standard symbol), although are consistent with the fits.

The correction we have applied to transform the current values of [Th/H], [Th/Fe], [Th/Nd], [Th/Eu], and [Th/Si] into the pristine values on the zero-age main sequence (ZAMS) are those given in Frebel et al. (2007) and Cayrel et al. (2001) adopting 14.05 Gyr as half lifetime of ²³²Th (τ_1/2(Th)), i.e. [Th/X]_ZAMS  = [Th/X]_now  |$+$| (ln(2).log(e)/τ_1/2(Th)).(t_star - t_Sun), where [Th/X]_ZAMS is the corrected abundance ratio relative to the Sun, [Th/X]_now means the observed current [Th/X], X is a stable element, t is the star isochrone age in Gyr, and t_Sun  =  4.56 Gyr. We have estimated the uncertainties in [Th/X]_ZAMS through the propagation of the errors in [Th/X]_now and isochrone age as well.

The pristine (corrected) [Th/H]_ZAMS correlates well with [Fe/H] along the short metallicity scale of solar twins, making [Th/Fe]_ZAMS not correlated with [Fe/H] (see Fig. 5, in which we plot both [Th/H]_ZAMS and [Th/Fe]_ZAMS as a function of [Fe/H]).

Figure 5.

[Th/H]_ZAMS versus [Fe/H] and [Th/Fe]_ZAMS versus [Fe/H]: A statistically robust linear fit is shown in each plot (blue solid line). The fitting equation with coefficients and their errors inside parenthesis are also displayed. The fit corresponds to 53 thin disc solar twins only (black filled circles). Data of excluded stars are plotted together for illustration purposes only (red empty triangles: four α-rich old stars, and magenta empty square: the chemically anomalous in s-elements HIP 64150). Sun’s data are also plotted as reference at the coordinates origin (green solar standard symbol), although are consistent with the fits.

We have fitted [Th/X] as a function of [Fe/H] (observed and ZAMS values for H and Fe) and also [Th/X] as a function of isochrone stellar age (observed and ZAMS values for H and Fe, and ZAMS values only for Nd, Eu, and Si) through a linear fitting approach that minimizes the orthogonal distance of the data points to the fitting curve that accounts for the uncertainties in both variables (it is an orthogonal distance linear regression with variance weighting in both x and y). All fits to the data have been performed using the Kapteyn kmpfit package.²

Fig. 6 shows plots of [Th/X]_ZAMS versus isochrone age (X: H, Fe, Nd, Eu, and Si). Table 4 presents the results of the statistically robust linear fits of [Th/X]_observed and [Th/X]_ZAMS as a function of both [Fe/H] (where X: H and Fe only) and isochrone age t (i.e. [Th/X]_ZAMS  =  a.[Fe/H]  |$+$|  b and [Th/X]_ZAMS  =  a.t  |$+$|  b).

Figure 6.

[Th/X]_ZAMS as a function of the isochrone stellar age (where X: H, Fe, Nd, Eu, and Si): a statistically robust linear fit is shown in each plot (blue solid line). The fitting equation with coefficients and their errors inside parenthesis are also displayed. The fits correspond to 53 thin disc solar twins only (black filled circles). Specifically for the plots [Th/Fe]_ZAMS–age and [Th/Eu]_ZAMS–age, an additional fit is also shown (red solid line) by excluding one borderline 3σ outlier (HIP 101905, a 1.2 Gyr old star relatively poor in Th), but the results do not change significantly. Data of excluded stars are plotted together for illustration purposes only (red empty triangles: four α-rich old stars, and magenta empty square: the chemically anomalous in s-elements HIP 64150). Sun’s data are also plotted as reference (green solar standard symbol), although are consistent with the fits.

Table 4.

Results of the linear fits of [Th/X] against [Fe/H] and isochrone stellar age t: [Th/X]  = a(± σ_a).[Fe/H]  |$+$| b(± σ_b) (observed and ZAMS [Th/X] values for X: H and Fe only), and [Th/X]  = a(± σ_a).t  |$+$| b(± σ_b) (observed and ZAMS [Th/X] values for X: H and Fe, and ZAMS [Th/X] values for X: Nd, Eu and Si). The slope divided by its error is shown in the last column. Specifically for [Th/Fe]_ZAMS–age and [Th/Eu]_ZAMS–age, results are also presented by excluding one borderline 3σ outlier (HIP 101905, a young star relatively poor in Th). These two extra fits are marked with an asterisk (*). The fitting approach is explained in Section 4.

	a	σa	b	σb	rms	χ2/ν	|a|/σa
[Th/X] versus [Fe/H]
			(dex)	(dex)	(dex)
[Th/H]obs	|$+$|0.904	0.044	|$+$|0.080	0.003	0.071	8.8	20.5
[Th/Fe]obs	−0.114	0.043	|$+$|0.080	0.003	0.071	9.2	2.6
[Th/H]ZAMS	|$+$|0.671	0.049	|$+$|0.081	0.003	0.040	2.6	13.7
[Th/Fe]ZAMS	−0.338	0.049	|$+$|0.082	0.003	0.040	2.6	6.9
[Th/X] versus t
	(dex.Gyr−1)	(dex.Gyr−1)	(dex)	(dex)	(dex)
[Th/H]obs	−0.0353	0.0016	|$+$|0.246	0.009	0.058	4.5	22.1
[Th/Fe]obs	−0.0247	0.0015	|$+$|0.195	0.008	0.048	3.5	16.5
[Th/H]ZAMS	−0.0117	0.0015	|$+$|0.138	0.008	0.056	5.2	7.8
[Th/Fe]ZAMS	−0.0009	0.0015	|$+$|0.086	0.008	0.047	3.6	0.6
*[Th/Fe]ZAMS	−0.0025	0.0015	|$+$|0.096	0.008	0.043	3.2	1.7
[Th/Nd]ZAMS	|$+$|0.0211	0.0017	−0.077	0.009	0.046	2.4	12.4
[Th/Eu]ZAMS	|$+$|0.0042	0.0016	−0.004	0.009	0.043	2.4	2.6
*[Th/Eu]ZAMS	|$+$|0.0027	0.0017	−0.005	0.009	0.039	2.1	1.6
[Th/Si]ZAMS	−0.0121	0.0018	|$+$|0.149	0.010	0.053	3.1	6.7

The four main results about [Th/H]_ZAMS and [Th/Fe]_ZAMS in solar twins as a function of [Fe/H] and isochrone age are as follows.

• Th/H_ZAMS correlates very well with [Fe/H] and the slope of this correlation is smaller than the unity (⁠|$+$| 0.671 ± 0.049). It is not a surprise that [Th/Fe]_ZAMS decreases with [Fe/H] over the restrict metallicity range of solar twins. Note that the relative uncertainty of the slope a makes the linear correlations [Th/H]_ZAMS–[Fe/H] and [Th/Fe]_ZAMS–[Fe/H] statistically significant (about 14 and 7 σ of significance, respectively). See Fig. 6 and Table 4. This suggests that the abundance of Th in the interstellar medium keeps growing in some way during the Galactic thin disc evolution, despite the natural radioactive decay of Th and the iron production driven by the contribution from SNIa relatively to SN-II. We refer the reader to the temporal increase of abundance ratios relative to iron of stable r-process elements such as Nd, Sm, Eu, Gd, and Dy, as reported by Spina et al. (2018) (see fig. 5 in their paper).

• The pristine Th abundance in the sample solar twins seems to have been super-solar on average during the evolution of Galactic thin disc and it has certainly become super-solar since the epoch of the Sun formation. There are no stars with [Th/H]_ZAMS < 0 dex with ages smaller than the solar one. We have obtained for the entire sample an average of [Th/H]_ZAMS  =  |$+$|0.080 dex with a standard deviation of 0.058 dex. It ranges from −0.117 up to |$+$|0.257 dex (from 76 to 181 per cent of the ZAMS solar value), showing a linear temporal increase from |$+$|0.037 dex 8.6 Gyr ago up to  |$+$|0.138 dex now (8.6 Gyr is the isochrone age of the oldest star in our analysed sample). Note that the uncertainty of the slope a of the relation [Th/H]_ZAMS–age corresponds to only 13 per cent of its value (about 8 σ of significance). See Fig. 6 and Table 4. The Sun at ZAMS seems to be somehow deficient in Th when compared against twin stars, specially for stars younger than itself. It is also interesting to note that all four α-rich stars are richer in Th on average than the other Galactic thin disc stars of similar ages; exactly as found for Nd, Sm and Eu in Spina et al. (2018).

• Th/Fe_ZAMS shows no correlation with the isochrone age. Note that the relative error of the slope is greater than 100 per cent. [Th/Fe]_ZAMS has likely kept super-solar during the Galactic thin disc evolution, showing an average of  |$+$|0.086 ± 0.008 dex (rms  =  0.047 dex). Super-solar Th/Fe ratios for youngest stars could be due to an increased production of Th relatively to Fe, perhaps due to neutron star mergers. Spina et al. (2018) also found that the slope of [X/Fe]–age for stable n-capture elements (in dex Gyr⁻¹ unity) seems to be anticorrelated with the solar s-process contribution percentage taken from Bisterzo et al. (2014) (see fig. 6 in Spina et al. 2018), such that the r-process elements Eu, Gd, Dy, and Sm (showing s-process contributions smaller than about 35 per cent) have the smallest slopes in modulus like Th, which exhibits a statistically null slope even after removing one borderline 3σ threshold outlier (HIP 101905 that is a young star relatively poor in Th for its metallicity). The slope of [Th/Fe]–age becomes −0.0025 ± 0.0015 without HIP 101905, which is not significantly different from zero (presenting a confidence level of about 1.7σ). See Fig. 6 and Table 4.

• There are uniform real dispersions of [Th/H]_ZAMS and [Th/Fe]_ZAMS over the whole isochrone age scale on average (rms close to 0.056 and 0.047 dex, respectively).

By comparing Th against neodymium (Nd) and europium (Eu), two other n-capture elements that are otherwise stable nuclei indeed, we have noticed that [Th/Nd]_ZAMS decreases with the stellar age, but [Th/Eu]_ZAMS has kept constant around solar ratio over the Galactic thin disc evolution. Whilst [Th/Nd]_ZAMS has a linear decrease from  |$+$|0.104 dex 8.6 Gyr ago up to −0.077 dex now on average (fitting rms  =  0.046 dex), the average of [Th/Eu]_ZAMS oscillates around  +0.014 dex with a standard deviation of 0.045 dex. The slope of the relation [Th/Nd]_ZAMS–age has a 12.4σ of significance, and the relation [Th/Eu]_ZAMS–age has just 2.4σ of significance. In case of removing one borderline outlier slightly beyond 3σ (HIP 101905 that is a young star relatively poor in Th for its metallicity), the slope of [Th/Eu]–age is even more compatible with zero (1.6σ of significance). See Fig. 6 and Table 4. These results confirm that Th follows Eu during the Galaxy’s thin disc evolution. However, it does not follow Nd. In fact, the contribution of s-process for the abundance Nd in the Solar system is estimated to be 57.5 ± 4.1 per cent, while for Eu it is 6.0 ± 0.4 per cent (Bisterzo et al. 2014; Spina et al. 2018). The decrease of [Th/Nd]_ZAMS with time may be explained by an increasing contribution from the s-process due to low-mass AGB stars in the Galactic thin disc as discussed by Spina et al. (2018). Therefore, our observations are in line with Nd not being a pure r-process element, but likely having an important s-process contribution, as already suggested by Bisterzo et al. (2014).

On the other hand, [Th/Si]_ZAMS has likely had a temporal increasing evolution from +0.045 dex 8.6 Gyr ago up to  +0.149 dex now on average (over the super-solar regime). The anticorrelation [Th/Si]_ZAMS–age has a slope with a relative uncertainty of 15 per cent only or about 7σ of significance (under a rms  =  0.053 dex; see Fig. 6 and Table 4). Silicon (Si) is one of the most abundant refractory elements in a terrestrial planet and its planetary abundance is closely related to the mantle mass (McDonough 2003). Thus, the abundance ratio Th/Si provides comparison of the Th abundance among telluric planets with different silicate mantle thickness (Unterborn et al. 2015). A higher [Th/Si]_ZAMS can indirectly indicate a higher probability of having a convective mantle in a terrestrial planet due to a higher radioactive energy budget from the Th decay.

Thus, our results suggest that most solar twins in the Galactic thin disc have the same capability as the Sun to have rocky planets with convective mantles, and perhaps suitable for life.

5 CONCLUSIONS

We have confirmed that there is a large energy budget from the Th decay for keeping the mantle convection and thickness of potential rocky planets around Sun-like stars, since the Galactic thin disc formation until now, because we have measured [Th/H]_ZAMS as super-solar on average during the Galactic thin disc life. Whilst the observed [Th/H] varies from −0.117 up to  +0.257 dex (from 76 up to 181 per cent of the current observed solar value), [Th/H]_ZAMS seems to have linearly changed on average from +0.037 dex 8.6 Gyr ago up to  +0.138 dex now (from 109 up to 137 per cent of the ZAMS solar value), showing a uniform dispersion of about 0.056 dex (linear fitting rms), such that [Th/H]_ZAMS seems to have certainly become super-solar for all analysed solar twin stars since the epoch of Solar system formation. [Th/Fe]_ZAMS has kept nearly constant and super-solar as well during the Galactic thin disc evolution (around  +0.086 ± 0.008 dex under a linear fitting rms equal to 0.047 dex). The Sun at ZAMS actually seems to be deficient in Th when specifically compared against younger solar twins. Whilst the observed [Th/Fe]_obs is linearly well anticorrelated with isochrone age, [Th/Fe]_ZAMS shows no relation with age. They both actually seem to have been super-solar on average during practically the entire Galactic thin disc evolution. [Th/Nd]_ZAMS has linearly decreased from  +0.104 dex 8.6 Gyr ago up to −0.077 dex now on average. [Th/Eu]_ZAMS has likely kept constant around solar ratio during the Galactic thin disc evolution, exhibiting a moderate dispersion of about 0.04 dex. This implies that Th follows Eu during the Galaxy’s thin disc evolution, but it does not follow Nd. The decrease of [Th/Nd]_ZAMS with time may be explained by a crescent contribution from the s-process in low-mass AGB stars, indicating that Nd is not exactly a pure r-process element. In fact, it has an estimated s-process production close to 60 per cent in the Solar system (Bisterzo et al. 2014). [Th/Si]_ZAMS has probably had a real increase from +0.045 dex 8.6 Gyr ago up to  +0.149 dex now on average.

Our results suggest that solar twin stars in the Galactic thin disc are as probable as the Sun to host rocky planets with convective mantles, and perhaps with suitable geological conditions for habitability.

SUPPORTING INFORMATION

Table 1. Stellar parameters of the sample stars collected from previous published works (covering 58 solar twins): photospheric parameters and isochrone age (Spina et al. 2018), macro-turbulence and rotation velocities (dos Santos et al. 2016), and elemental abundances (Bedell et al. 2018).

Table 2. Line list after the gf calibration to the solar spectrum: wavelength, species code, excitation potential of transition inferior level, log gf, and species identification.

Table 3. Thorium abundance measured in this work relatively to the Sun for 58 solar twins (Sun included as well).

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

ACKNOWLEDGEMENTS

This work has made use of the VALD database, operated at Uppsala University, the Institute of Astronomy RAS in Moscow, and the University of Vienna. RB acknowledges the PhD grant from Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES) (093.875.006-24). AM acknowledges the research grant from Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq) (309562/2015-5). JM thanks Fundação de Amparo à Pesquisa do Estado de São Paulo(FAPESP) (2012/24392-2). The research by MA has been supported by the Australian Research Council (grants FL110100012 and DP150100250). We all are grateful for the anonymous referee by his valuable revision to improve the paper.

Footnotes

1

Log ε(X)  =  log (n(X)/n(H))  |$+$|  12, where n is number density.

2

https://www.astro.rug.nl/software/kapteyn/index.html