THE MASS–METALLICITY RELATION FOR GIANT PLANETS

Giant planets do not directly take on the composition of their parent stars. If some flavor of core-accretion formation is correct (Pollack et al. 1996), then a seed core of ∼5–10 M_⊕ of ice/rock must build up first, which begins accreting nebular gas. The gas need not share the exact composition of the parent star, due to condensation and migration of solids within the disk (Lodders 2009; Öberg et al. 2011). The growing giant planet accretes gas but is also bombarded by planetesimals, which may add to the core mass or dissolve into the growing H/He envelope. The amounts and variety of heavy elements (metals) accreted by the growing planet will depend on its formation location, formation time, disk environments sampled, and whether it forms near neighboring planets, among many other factors. A planet's present-day composition is our indirect window into its formation process.

The observed atmospheric composition of fully convective giant planets reflect the mixing ratios of heavy elements within their whole H/He envelope (with some caveats—see Section 3.1). Spectroscopy of the solar system's giants points to enhancements in the mixing ratio of carbon (as seen in CH₄) of ∼4, 10, 80, and 80, in Jupiter, Saturn, Uranus, and Neptune, respectively (Wong et al. 2004; Fletcher et al. 2009). The Galileo Entry Probe found that Jupiter's atmosphere is enhanced in volatiles like C, N, S, and noble gases by factors of two to five compared to the Sun (Wong et al. 2004).

Remarkably, one can make inferences about the bulk metallicity, or heavy element enhancement, Z_planet/Z_Sun, of a solar system giant planet by measuring only its mass and radius. These results are consistent with and refined by measurements of their gravitational moments. By comparing to structure models, it is straightforward to infer from mass and radius alone that Jupiter and Saturn are both smaller and denser than they would be if they were composed of solar composition material (Zapolsky & Salpeter 1969; Fortney et al. 2007). Therefore, we can infer that they are enhanced in heavy elements. As we describe below, this leads to an immediate connection with transiting giant planets, where mass and radius can be accurately measured.

Knowledge of the heavy element enrichment of our solar system's giant planets has led to dramatic advances in our understanding of planet formation. Models of the core-accretion model of planet formation have advanced to the point where they can match the heavy element enrichment of each of the solar system's giant planets (Alibert et al. 2005; Lissauer et al. 2009; Helled & Bodenheimer 2014). However, we are only beginning to attain similar data for exoplanets, which will provide a critical check for planet formation models over a tremendously larger phase space. Such constraints are particularly important for comparison with planetary population synthesis models that aim to understand the processes of core formation, H/He envelope accretion, and planetary migration, in diverse protoplanetary environments (Ida & Lin 2004; Mordasini et al. 2014). Constraints on planet formation from exoplanetary systems have been almost entirely driven by data on the frequency of planets in the mass/semimajor axis plane. A promising avenue is to move these comparisons into a complementary plane, that of planetary composition, either in bulk composition (Guillot et al. 2006; Burrows et al. 2007; Miller & Fortney 2011; Mordasini et al. 2014), or atmospheric composition (Madhusudhan et al. 2011; Fortney et al. 2013; Konopacky et al. 2013; Kreidberg et al. 2014; Barman et al. 2015) as shown in a schematic way in Figure 1.

Zoom In Zoom Out Reset image size

The dramatic rise in the number of observed transiting exoplanets provides a unique opportunity. With radii derived from transit observations and masses derived from radial-velocity or transit-timing variation measurements, we get especially detailed information about these objects. This gives us a measured density, and therefore some rough information about their bulk composition. A more advanced analysis uses models of planet structural evolution (e.g., Fortney et al. 2007) to constrain the quantity of heavy elements.

Most of the giant planets we have observed are strongly irradiated hot Jupiters, whose radii are inflated beyond what models predict. Much effort has been put into understanding this discrepancy. A thorough discussion is outside the scope of this article, but the various proposed inflation mechanisms are extensively reviewed in Fortney & Nettelmann (2010), Baraffe et al. (2014), and Weiss et al. (2013). Unfortunately, without a definite understanding of the inflation process, this acts as a free parameter in modeling: the inflationary effect enlarges a planet and added heavy elements shrink it, resulting in a degeneracy that inhibits our ability to obtain useful composition constraints. Still, work has been done to use models to address a the star–planet composition connection, using plausible assumptions about the effect, as a heat source (Guillot et al. 2006) or as a slowed-cooling effect (Burrows et al. 2007). Both studies saw an increase in planet heavy element mass with stellar metallicity.

A promising avenue of investigation are the sample of transiting exoplanets, which are relatively cool. Planets that receive an incident flux below around 2 × 10⁸ erg s⁻¹ cm⁻² (T_eq ≲ 1000 K) appear to be non-inflated (Demory & Seager 2011; Miller & Fortney 2011), obviating the need for assumptions about that effect. Figure 2 shows this threshold. Miller & Fortney (2011; hereafter referred to as MF2011) studied these planets, finding correlations in the heavy element mass with planetary mass and stellar metallicity. In particular, they noted a strong connection between the relative enrichment of planets relative to their parent stars (Z_planet/Z_star) and the planet mass. However, that study was limited by a small sample size of 14 planets.

Zoom In Zoom Out Reset image size

In our work that follows, we consider the set of cool transiting giant planets, now numbering 47, and compare them to a new grid of evolution models to estimate their heavy element masses, and we include a more sophisticated treatment of the uncertainties on our derived planetary metal-enrichments. We then examine the connections between their mass, metal content, and parent star metallicity.

Our data was downloaded from the Extrasolar Planets Encyclopedia (exoplanets.eu, Schneider et al. 2011) and the NASA Exoplanet Archive (Akeson et al. 2013). These data were combined, filtered (see below), and checked against sources for accuracy. Some corrections were needed, mostly in resolving differing values between sources; we aimed to consistently use the most complete and up-to-date sources. We tried to include all known planets who met the selection (even if some of their data were found in the literature and not the websites). Critically, we use data from the original sources, rather than the websites (see Table 1).

Table 1. 

#	Name	Mass (MJ)	Radius (RJ)	Age (Gyr)	Flux ($\tfrac{\mathrm{erg}}{{\mathrm{cm}}^{2}\cdot {\rm{s}}}$)	[Fe/H] (Dex)	Metal (M⊕)	Zp	Zp/Zs
1	CoRoT-8 b	0.22 ± 0.03	0.57 ± 0.02	0.1–3.0	1.22 × 108	0.30 ± 0.10	${45.91}_{-6.83}^{+6.99}$	0.66 ± 0.03	${24.18}_{-6.67}^{+5.00}$
2†	CoRoT-9 b	0.84 ± 0.07	1.05 ± 0.04	1.0–8.0	6.59 × 106	-0.01 ± 0.06	${18.87}_{-18.77}^{+13.30}$	${0.07}_{-0.07}^{+0.04}$	${5.21}_{-5.21}^{+2.81}$
3	CoRoT-10 b	2.75 ± 0.16	0.97 ± 0.07	0.1–3.0	5.38 × 107	0.26 ± 0.07	${202.98}_{-78.08}^{+84.09}$	0.23 ± 0.09	${9.26}_{-4.44}^{+3.63}$
4	CoRoT-24 c	0.09 ± 0.04	0.44 ± 0.04	8.0–14.0	5.67 × 107	0.30 ± 0.15	${21.64}_{-8.17}^{+9.38}$	0.77 ± 0.06	${29.18}_{-13.28}^{+8.55}$
5	GJ 436 b	0.07 ± 0.01	0.37 ± 0.02	0.5–10.0	4.63 × 107	-0.32 ± 0.12	${20.07}_{-1.62}^{+1.68}$	0.85 ± 0.03	${132.02}_{-43.68}^{+31.66}$
6	HAT-P-11 b	0.08 ± 0.01	0.42 ± 0.01	2.4–12.4	1.31 × 108	0.31 ± 0.05	${20.31}_{-2.36}^{+2.45}$	0.79 ± 0.02	${27.76}_{-3.56}^{+3.09}$
7	HAT-P-12 b	0.21 ± 0.01	${0.96}_{-0.02}^{+0.03}$	0.5–4.5	1.91 × 108	-0.29 ± 0.05	${11.83}_{-3.23}^{+3.17}$	0.18 ± 0.05	${24.75}_{-7.43}^{+6.92}$
8†	HAT-P-15 b	1.95 ± 0.07	1.07 ± 0.04	5.2–9.3	1.51 × 108	0.22 ± 0.08	${42.45}_{-42.35}^{+30.67}$	${0.07}_{-0.07}^{+0.04}$	${3.00}_{-3.00}^{+1.65}$
9	HAT-P-17 b	0.53 ± 0.02	1.01 ± 0.03	4.5–11.1	8.97 × 107	0.00 ± 0.08	${14.08}_{-6.24}^{+6.56}$	0.08 ± 0.04	${6.02}_{-3.25}^{+2.72}$
10	HAT-P-18 b	0.18 ± 0.03	0.95 ± 0.04	6.0–16.8	1.18 × 108	0.10 ± 0.08	${6.08}_{-3.29}^{+4.35}$	${0.10}_{-0.06}^{+0.05}$	${5.88}_{-4.06}^{+3.11}$
11	HAT-P-20 b	7.25 ± 0.19	0.87 ± 0.03	2.9–12.4	2.00 × 108	0.35 ± 0.08	${662.48}_{-109.11}^{+111.85}$	0.29 ± 0.05	${9.35}_{-2.66}^{+2.06}$
12	HAT-P-54 b	0.76 ± 0.03	0.94 ± 0.03	1.8–8.2	1.04 × 108	-0.13 ± 0.08	${47.33}_{-9.22}^{+9.60}$	0.20 ± 0.04	${19.06}_{-5.80}^{+4.54}$
13	HATS-6 b	0.32 ± 0.07	1.00 ± 0.02	0.1–13.7	5.84 × 107	0.20 ± 0.09	${8.67}_{-4.70}^{+7.17}$	${0.08}_{-0.06}^{+0.04}$	${3.80}_{-3.24}^{+1.98}$
14	HATS-17 b	1.34 ± 0.07	0.78 ± 0.06	0.8–3.4	9.97 × 107	0.30 ± 0.03	${196.42}_{-34.17}^{+35.50}$	0.46 ± 0.08	${16.56}_{-3.12}^{+2.93}$
15	HD 17156 b	3.19 ± 0.03	1.09 ± 0.01	2.4–3.8	1.98 × 108	0.24 ± 0.05	${58.84}_{-8.75}^{+8.81}$	0.06 ± 0.01	${2.40}_{-0.49}^{+0.43}$
16	HD 80606 b	3.94 ± 0.11	0.98 ± 0.03	1.7–7.6	1.65 × 107	0.43 ± 0.03	${215.69}_{-52.23}^{+52.63}$	0.17 ± 0.04	${4.58}_{-1.19}^{+1.13}$
17	K2-19 b	${0.19}_{-0.04}^{+0.02}$	0.67 ± 0.07	8.0–14.0	1.49 × 108	0.19 ± 0.12	${28.72}_{-7.56}^{+9.21}$	0.48 ± 0.09	${22.99}_{-9.31}^{+6.59}$
18	K2-24 b	0.06 ± 0.01	0.52 ± 0.05	1.3–8.9	8.24 × 107	0.42 ± 0.04	${13.38}_{-2.86}^{+3.30}$	0.68 ± 0.06	${18.44}_{-2.56}^{+2.27}$
19	K2-24 c	0.08 ± 0.02	0.72 ± 0.07	1.3–8.9	3.20 × 107	0.42 ± 0.04	${10.81}_{-3.32}^{+3.55}$	${0.42}_{-0.09}^{+0.10}$	11.33 ± 2.73
20	K2-27 b	0.09 ± 0.02	0.40 ± 0.03	0.5–10.0	1.58 × 108	0.14 ± 0.07	${24.40}_{-6.42}^{+6.86}$	0.84 ± 0.04	${43.88}_{-8.19}^{+6.65}$
21	Kepler-9 b	0.25 ± 0.01	0.84 ± 0.07	2.0–4.0	8.53 × 107	0.12 ± 0.04	${23.63}_{-6.80}^{+7.15}$	${0.30}_{-0.09}^{+0.08}$	${16.06}_{-5.11}^{+4.69}$
22	Kepler-9 c	0.17 ± 0.01	0.82 ± 0.07	2.0–4.0	3.29 × 107	0.12 ± 0.04	${16.14}_{-4.58}^{+4.81}$	0.30 ± 0.08	${16.17}_{-4.95}^{+4.61}$
23*	Kepler-16 b	0.33 ± 0.02	0.75 ± 0.00	0.5–10.0	4.19 × 105	-0.30 ± 0.20	${41.49}_{-2.78}^{+3.11}$	${0.39}_{-0.02}^{+0.01}$	${62.05}_{-40.47}^{+22.31}$
24†	Kepler-30 c	2.01 ± 0.16	1.10 ± 0.04	0.2–3.8	1.12 × 107	0.18 ± 0.27	${42.54}_{-42.44}^{+30.62}$	${0.07}_{-0.07}^{+0.04}$	${3.81}_{-3.81}^{+2.30}$
25	Kepler-30 d	0.07 ± 0.01	0.79 ± 0.04	0.2–3.8	4.05 × 106	0.18 ± 0.27	${8.70}_{-2.02}^{+1.95}$	${0.38}_{-0.07}^{+0.08}$	${21.39}_{-21.53}^{+10.03}$
26*	Kepler-34 b	0.22 ± 0.01	0.76 ± 0.01	5.0–6.0	3.21 × 106	-0.07 ± 0.15	${24.42}_{-1.90}^{+1.97}$	0.35 ± 0.02	${31.24}_{-14.14}^{+8.98}$
27*	Kepler-35 b	0.13 ± 0.02	0.73 ± 0.01	8.0–12.	5.06 × 106	-0.34 ± 0.20	${14.56}_{-2.40}^{+2.62}$	0.36 ± 0.02	${62.22}_{-38.71}^{+21.89}$
28	Kepler-45 b	0.51 ± 0.09	0.96 ± 0.11	0.4–1.5	8.83 × 107	0.28 ± 0.14	${37.08}_{-17.71}^{+23.05}$	${0.23}_{-0.12}^{+0.11}$	${9.12}_{-7.11}^{+4.46}$
29†	Kepler-75 b	10.10 ± 0.40	1.05 ± 0.03	3.4–9.7	1.29 × 108	0.30 ± 0.12	${0.00}_{-0.00}^{+133.99}$	${0.00}_{-0.00}^{+0.04}$	${0.00}_{-0.00}^{+10.84}$
30	Kepler-89 d	0.16 ± 0.02	0.98 ± 0.01	3.7–4.2	1.57 × 108	0.01 ± 0.04	${5.25}_{-1.10}^{+1.23}$	0.10 ± 0.01	${7.02}_{-1.21}^{+1.10}$
31†	Kepler-117 c	1.84 ± 0.18	1.10 ± 0.04	3.9–6.7	5.79 × 107	-0.04 ± 0.10	${25.14}_{-25.04}^{+22.10}$	${0.04}_{-0.04}^{+0.02}$	${3.44}_{-3.44}^{+2.02}$
32	Kepler-145 c	0.25 ± 0.05	0.39 ± 0.01	2.1–3.1	8.03 × 107	0.13 ± 0.10	${73.86}_{-16.42}^{+17.21}$	0.92 ± 0.03	${50.39}_{-13.43}^{+10.20}$
33*	Kepler-413 b	0.21 ± 0.07	0.39 ± 0.01	0.5–10.0	3.20 × 106	-0.20 ± 0.20	${60.44}_{-20.85}^{+21.91}$	${0.89}_{-0.03}^{+0.04}$	${112.32}_{-70.60}^{+40.09}$
34	Kepler-419 b	2.50 ± 0.30	0.96 ± 0.12	1.6–4.1	4.81 × 107	0.09 ± 0.15	${198.16}_{-100.72}^{+126.80}$	${0.25}_{-0.15}^{+0.13}$	${15.41}_{-13.45}^{+8.16}$
35	Kepler-420 b	1.45 ± 0.35	0.94 ± 0.12	6.3–12.3	1.57 × 107	0.27 ± 0.09	${110.51}_{-59.61}^{+88.89}$	${0.24}_{-0.16}^{+0.12}$	${9.25}_{-7.07}^{+4.90}$
36†	Kepler-432 b	5.84 ± 0.05	1.10 ± 0.03	2.6–4.2	1.73 × 108	-0.02 ± 0.06	${67.39}_{-67.29}^{+60.95}$	${0.04}_{-0.04}^{+0.02}$	${2.74}_{-2.74}^{+1.61}$
37	Kepler-539 b	0.97 ± 0.29	0.75 ± 0.02	0.4–1.2	5.54 × 106	-0.01 ± 0.07	${154.94}_{-48.76}^{+51.90}$	0.50 ± 0.03	${36.74}_{-7.07}^{+5.75}$
38*	Kepler-1647 b	1.52 ± 0.65	1.06 ± 0.01	3.9–4.9	7.91 × 105	-0.14 ± 0.05	${29.22}_{-15.14}^{+19.76}$	0.05 ± 0.02	${5.34}_{-2.11}^{+2.10}$
39	WASP-8 b	${2.24}_{-0.09}^{+0.08}$	${1.04}_{-0.05}^{+0.01}$	3.0–5.0	1.75 × 108	0.17 ± 0.07	${84.34}_{-38.00}^{+43.33}$	${0.12}_{-0.06}^{+0.05}$	${5.79}_{-3.27}^{+2.65}$
40	WASP-29 b	0.24 ± 0.02	${0.79}_{-0.04}^{+0.06}$	7.0–13.0	2.05 × 108	0.11 ± 0.14	${25.65}_{-6.09}^{+6.63}$	${0.33}_{-0.08}^{+0.07}$	${19.28}_{-9.52}^{+6.18}$
41	WASP-59 b	0.86 ± 0.04	0.78 ± 0.07	0.1–1.2	4.41 × 107	-0.15 ± 0.11	${128.37}_{-25.06}^{+26.73}$	0.47 ± 0.09	${48.78}_{-18.26}^{+13.21}$
42†	WASP-69 b	0.26 ± 0.02	1.06 ± 0.05	0.5–4.0	1.94 × 108	0.14 ± 0.08	${7.23}_{-7.13}^{+4.96}$	${0.09}_{-0.09}^{+0.05}$	${4.56}_{-4.56}^{+2.56}$
43	WASP-80 b	0.54 ± 0.04	1.00 ± 0.03	0.5–10.0	1.06 × 108	$-{0.13}_{-0.17}^{+0.15}$	${20.31}_{-7.74}^{+8.47}$	${0.12}_{-0.05}^{+0.04}$	${12.28}_{-8.82}^{+5.37}$
44	WASP-84 b	0.69 ± 0.03	0.94 ± 0.02	0.4–1.4	9.26 × 107	0.00 ± 0.10	${53.82}_{-6.48}^{+6.79}$	0.24 ± 0.03	${17.88}_{-5.38}^{+4.00}$
45	WASP-130 b	1.23 ± 0.04	0.89 ± 0.03	2.0–14.0	1.09 × 108	0.26 ± 0.10	${108.83}_{-17.35}^{+17.58}$	0.28 ± 0.04	${11.22}_{-3.63}^{+2.74}$
46	WASP-132 b	0.41 ± 0.03	0.87 ± 0.03	0.5–14.0	7.62 × 107	0.22 ± 0.13	${33.70}_{-6.40}^{+6.87}$	0.26 ± 0.05	${11.63}_{-4.92}^{+3.38}$
47	WASP-139 b	0.12 ± 0.02	0.80 ± 0.05	0.3–1.0	1.61 × 108	0.20 ± 0.09	${15.64}_{-2.48}^{+2.93}$	0.42 ± 0.05	${19.40}_{-5.31}^{+4.07}$

Note. A * indicates circumbinary planets. A † indicates those planets with results adjusted to reflect the fact that a certain portion of their samples could not be modeled (see Section 4 for discussion). Sources—1: Bordé et al. (2010), 2: Deeg et al. (2010), 3: Bonomo et al. (2010), 4: Bordé et al. (2010), 5: Bean et al. (2006); Southworth (2010), 6: Bakos et al. (2010), 7: Hartman et al. (2009), 8: Kovács et al. (2010), 9: Howard et al. (2012), 10: (Hartman et al. (2011)), 11: Bakos et al. (2011), 12: Bakos et al. (2015), 13: Hartman et al. (2015), 14: Brahm et al. (2015), 15: Nutzman et al. (2011), 16: Pont et al. (2009) & Naef et al. (2001) & Saffe et al. (2005), 17: Nespral et al. (2016), 18: Petigura et al. (2015), 19: Petigura et al. (2015), 20: Van Eylen et al. (2016), 21: Holman et al. (2010), 22: Holman et al. (2010), 23: Doyle et al. (2011), 24: Sanchis-Ojeda et al. (2012), 25: Sanchis-Ojeda et al. (2012), 26: Welsh et al. (2012), 27: Welsh et al. (2012), 28: Johnson et al. (2012), 29: Bonomo et al. (2015), 30: Masuda et al. (2013) & Hirano et al. (2012), 31: Bruno et al. (2015), 32: Xie (2014), 33: Kostov et al. (2014), 34: Dawson et al. (2014) & Dawson et al. (2012), 35: Santerne et al. (2014), 36: Ortiz et al. (2015), 37: Mancini et al. (2016), 38: Kostov et al. (2015), 39: Queloz et al. (2010), 40: Hellier et al. (2010), 41: Hébrard et al. (2013), 42: Anderson et al. (2014), 43: Triaud et al. (2015), 44: Anderson et al. (2014), 45: Hellier et al. (2016), 46: Hellier et al. (2016), 47: Hellier et al. (2016).

A machine-readable version of the table is available.

Download table as:  DataTypeset image

For our sample, we selected the cool giant planets that had well-determined properties. Typically, this means they were the subject of both transit and radial velocity studies. The mass and radius uncertainties were of particular importance. Planets needed to have masses between 20 M_⊕ and 20 M_J, and relative uncertainties thereof below 50%. Our sample's relative mass uncertainties were typically well below that cutoff, distributed as ${10}_{-5.7}^{+12.8} \%$. We also constrained relative radius uncertainty to less than 50%, but again saw values much lower (${5.0}_{-2.5}^{+4.6} \%$). Both uncertainty cuts were made to eliminate planets with only rough estimates or upper limits for either value.

As discussed in Section 1, we used a 2 × 10⁸ erg s⁻¹ cm⁻² upper flux cutoff to filter out potentially inflated planets. Consequently, candidates needed to have enough information to compute the time-averaged flux: stellar radius and effective temperature (for both stars, if binary), semimajor axis, and eccentricity. In addition, we needed measured values for the stellar metallicities in the form of the iron abundance [Fe/H]. These tended to have fairly high uncertainties, and were a major source of error in our determination of Z_planet/Z_star.

Determining the age of a planet is necessary to use evolution models to constrain its metal content; unfortunately, this is often a difficult value to obtain. Our ages come from stellar ages listed in the literature, typically derived from a mixture of gyrochronology and stellar evolution models. These methods can produce values with sizable uncertainties, typically given as plausible value ranges. We treat these as flat probability distributions (a conservative choice), and convert values given as Gaussian distributions to 95% confidence intervals.

Because planets do most of their cooling, and therefore contraction, early in their lives (see Fortney et al. 2007 for a discussion), large uncertainties in age are not a major obstacle in modeling older planets. For planets who may be younger than a few gigayears, we cannot rule out that they are very metal-rich, but young and puffy, as the two effects would cancel out. We account for this in our analysis; planets that may be young consequently exhibit higher upper bound uncertainties in heavy element mass. For Kepler-16 (AB)-b, Kepler-413 (AB)-b, and WASP-80 b, no age was given, so we used a range of 0.5–10 Gyr. This is reasonable because it represents the possibilities of the planets being either young or old and because the age was a second-order effect on our heavy element assessment (after mass and radius) as seen in Figure 3.

Zoom In Zoom Out Reset image size

We created one-dimensional planetary models consisting of an inert core composed of a 50/50 rock–ice mixture, a homogenous convective envelope made of a H/He–rock–ice mixture, and a radiative atmosphere as the upper boundary condition. The models are made to satisfy the equations of hydrostatic equilibrium, mass conservation, and the conservation of energy:

$$\begin{eqnarray}&&\displaystyle \frac{\partial P}{\partial m}=-\displaystyle \frac{{Gm}}{4\pi {r}^{4}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial r}{\partial m}=\displaystyle \frac{1}{4\pi {r}^{2}\rho },\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial L}{\partial m}=-T\displaystyle \frac{\partial S}{\partial T}.\end{eqnarray} \tag{ 3 }$$

Equation (3) is not used in the core, where luminosity is neglected. Structures are initially guessed, then iteratively improved until these conditions are met to within sufficiently small error. This computation was done in a new Python code created for this study. Its advantages are its relative speed (a six million planet grid takes <1 hr to create) and its ability to easily try different compositional structures (e.g., heavy elements in the core versus envelope) and/or equations of state. As a simple visual diagnostic, Figure 3 shows the output of our models.

For our atmosphere models, we interpolate on the solar metallicity grids from Fortney et al. (2007), as was done in several other works (including Miller et al. 2009 and Lopez & Fortney 2014). With these models, we compute the intrinsic luminosity (L) of a planet model (see Equation (3)), which describes the rate that energy escapes from the interior, at a given surface gravity and isentropic interior profile. These grids are then used to determine the rate of entropy change in the envelope, and the contraction history of a given model planet. The initial entropy is not important; reasonable initial values typically all converge within a few hundred million years (Marley et al. 2007). As such, we do not need to consider if planets form in a hot or cold start scenario. Following Miller et al. (2009), we include the small extension in radius due to the finite thickness of the radiative atmosphere.

A fully self-consistent treatment of the atmosphere would use a range of metal-enriched atmosphere grids to be interpolated to yield consistency between atmospheric metallicity and the H/He envelope metal mass fraction. Metal-enriched grids would tend to slow cooling and contraction (Burrows et al. 2007) and yield larger heavy element masses than we present here. However, given uncertainties in the upper boundary condition in strongly irradiated giant planets (see, e.g., Guillot 2010 for an analytic analysis, Spiegel & Burrows 2013 for a review of the chemical and physical processes at play, and Guillot & Showman 2002 for an application of a 2D boundary condition) and our uncertainty in where the heavy elements are within the planets (core versus envelope), it is not clear if such an expanded study is yet warranted.

In creating our models, we had several important factors to consider: where the heavy elements are within the planet, what the composition of these heavy elements is, and how to treat the thermal properties of the core. These questions do not have any clearly superior or well-established solutions. The uncertainties from them, however, were overshadowed by the large uncertainties from available values for mass, radius, and (to a lesser extent) age. Therefore, in designing our models, we chose plausible solutions instead of attempting to incorporate all modeling uncertainties into our results.

3.1. Heavy Element Distribution

3.2. Equations of State

To apply our models to an observed planet, we take draws from the probability distributions implied by the measurement uncertainties in mass, radius, and age. Each draw therefore consists of a mass, radius, and age for one planet. The probability distributions are normal, except for the age, which is conservatively modeled as a flat range (see Section 2). For each draw, we compute the inferred heavy-element mass. By making many (10,000) draws, we have an estimate of the range of heavy-element masses consistent with observations. This procedure is done for each planet in the sample.

The resulting distributions were single peaked and roughly Gaussian overall (see Figure 4). For some of our planets, the uncertainty was dominated by the mass or radius, as evidenced by a correlation among the draws for each individual planet. Mass error more commonly dominated the uncertainty at high Z_planet, and radius error at low Z_planet. Our reported values show the marginal mean of each distribution, with upper and lower uncertainties computed as the rms deviation from the mean from draws above the mean and below the mean, respectively. These represent the data reasonably well, but care should be taken not to overlook the correlation with input variables. For our part, we do all computations directly on the samples (see below) and report the uncertainties in the resulting distributions.

Zoom In Zoom Out Reset image size

Some planets whose Z_planet values were clustered near zero (pure H/He) or one (pure ice/rock) generated draws that could not be recreated in our models. This occurred if (for example) the draw was randomly assigned a low value for mass and a high value for its radius, such that it was larger than an analogous pure H/He object. These draws were discarded, but we noted how often this occurred. For the six planets where this occurred 15.8%–50% of the time (where at least a 1σ tail is outside the valid range), the error-bar on that side was adjusted to a Z_planet of exactly zero or one, as appropriate. These are marked with a † in Table 1. For the massive H/He dominated planet Kepler-75b, this occurred more than half the time, so its heavy element derived values are listed as upper limits.

To conduct the regressions described in our results (Section 5), we use a linear Bayesian regression on the log of the values of interest, using uninformative priors on the slope and intercept. The likelihood for the regression was

$$\begin{eqnarray}&&P({\boldsymbol{y}}| {\boldsymbol{X}},{\boldsymbol{\beta }},{\sigma }^{2})=\mathrm{Normal}({\boldsymbol{X}}{\boldsymbol{\beta }},{\sigma }^{2}{\boldsymbol{I}})\end{eqnarray} \tag{ 4 }$$

where ${\boldsymbol{y}}$ is the vector of y points, $X=[{\boldsymbol{1}},{\boldsymbol{x}}]$ is the matrix of covariates, ${\boldsymbol{\beta }}$ is the coefficients vector (essentially [b,m] from $y={mx}+b$), and I is the identity matrix. Using the standard noninformative prior $P(\beta ,{\sigma }^{2}| {\boldsymbol{X}})\propto 1/{\sigma }^{2}$ and the distribution of ${\boldsymbol{y}}$ and the covariate ${\boldsymbol{x}}$ as $P({\boldsymbol{y}},{\boldsymbol{x}})={\prod }_{i}{P}_{i}({y}_{i},{x}_{i})$, we derived the full conditional distributions and implemented a Gibbs sampler. For each fit, we initialized with the classical fit, and had a burn in of 1000 steps and thinned the results by keeping only every tenth step to ensure the results were well-mixed. Our fits were performed in logspace so that they were effectively power-law fits.

We considered the possibility that mass and radii observations are not independent. If they are correlated, then sampling more directly from observational data (or posteriors thereof) could improve the uncertainties in our heavy element masses. Such an operation would need to be careful not to extract more information than the data actually provide. Southworth et al. (2007) describes a method of computing surface gravity that is more precise than directly using a planet's derived mass and radius. They instead use orbital parameters, the planet's radius, and the stellar reflex velocity K. These are closer to the observed quantities, and so avoid unnecessarily compounding uncertainties. We used Southworth's formula as a proxy to estimate how much improvement in uncertainty we might see from such an approach. We found that most of our planets exhibited only modestly better gravity estimates (though a handful were substantially improved, such as HAT-P-18 b). As such, we determined that trying to work more directly from observation data would not be worthwhile for our set of planets. Still, we suggest that studies examining individual planets should consider this approach.

4.1. Modeling Uncertainty

In the preceding sections, we listed a number of possible sources of modeling uncertainty. We will now argue that these uncertainties, while present, do not significantly affect our results, especially the fits described in Section 5. To the extent that they are affected, the error should be concentrated in the coefficient of our fits (not the power) because to first order these uncertainties would affect all planets equally.

First, we consider the effect of the heavy-element distribution on the structure and evolution of the planets. Considering the two extreme cases, where the metal is either entirely confined to a core or entirely dissolved into a homogeneous envelope, we compared the resulting models of four representative giant planets against our preferred models. As can be seen in Figure 5, these different choices have a clear effect on the inferred metal abundance, but the effect is small compared to observational uncertainties.

Zoom In Zoom Out Reset image size

To evaluate the effect of EOS uncertainty, we considered the difference between the Saumon-Chabrier H/He EOS that we used and the Militzer & Hubbard (2013) EOS. This EOS is computed from DFT-MD simulations, which may be more accurate than the semi-analytic SCvH EOS. We were unable to use it for this work because it only covers densities up to those found in roughly Jupiter-mass planets. Figure 12 from Militzer & Hubbard (2013) shows that for envelope entropies typical of older planets, the deviation in the resulting radius is about 10%. To match this, we might require as much as 15% less metal. In practice, the amount is somewhat lower due to next-order effects: e.g., smaller planets evolve slower (less surface area to emit from) and the metal EOS of the planets is unaffected by a H/He EOS change.

To quantify this, we derived inferred metal-masses of our planets using the Militzer–Hubbard EOS where possible. The results for four planets are shown in Figure 6. This sample of planets is not representative in mass because the Militzer–Hubbard EOS does not extend to high enough pressures to model super-Jupiters. Most of the results were fairly similar (bottom row), but a few, generally young planets, exhibited more significant differences. The choice of EOS matters, but is usually a next-order effect after observational uncertainties.

Zoom In Zoom Out Reset image size

As a reference, we applied our model to Jupiter and Saturn. Since these planets have well-determined properties that include some gravitational moments, we can use them as a test of our model's validity. Our inferred heavy-element mass should resemble estimates from better-constrained models that make use of these gravitational moments (e.g., Guillot 1999) but the same H/He EOS. A state-of-the-art model for Jupiter in Hubbard & Militzer (2016) favors metal masses around 22 M_⊕ (but note that it uses a different EOS).

As we see from Table 2, our inferred metal masses fall within a plausible range for Jupiter and Saturn. Furthermore, we show the errors resulting from uncertainties in mass and radius (10% each) to demonstrate that these are the dominant sources of uncertainty in our study. Also, we see that the error from age uncertainty for these somewhat old planets is not very significant.

Table 2.  Inferred Total Heavy-element Mass for Jupiter and Saturn, from Guillot (1999) and This Work

Source	Jupiter	Saturn
Guillot (1999)	10–40	20–30
This Work	37	27
±10% M,R	±20	±5.5
±2 Gigayears	±1.1	±.8

Note. For reference, we also show the uncertainties that would result if we had 10% uncertainties in mass and radius, and separately for a 2 Gyr uncertainty in age. Note that the central values lie within the estimate from Guillot.

Download table as:  ASCIITypeset image

We examined our results for connections between three quantities connected through the core-accretion model: the planetary mass M, its heavy element mass M_z, and the stellar metallicity [Fe/H]. We also considered the metal mass fraction of the planet Z_planet and its ratio to that of the parent star Z_planet/Z_star.

5.1. Relation to Planet Mass

There exists a clear correlation between planet mass and heavy element content. We can see that as we move toward more massive planets, the total mass of heavy elements increases, but the bulk metallicity decreases (Figure 7). Using Kendall'(s) Tau with the mean values of metal mass and total planet mass, we measure a correlation of 0.4787 and a p-value of 2.07 × 10⁻⁶, strongly supporting a correlation. This is consistent with a formation model where an initial heavy element core accretes predominantly, but not exclusively, H/He gas (e.g., Pollack et al. 1996). Indeed, it appears likely that all of our sample planets have more than a few M_⊕ of heavy-elements and usually far more (though we cannot determine a minimum exactly), consistent with both MF2011 and theoretical core-formation models (Klahr & Bodenheimer 2006). A fit to the log of the data gives ${M}_{z}=(57.9\pm 7.03){M}^{(.61\pm .08)}$, or roughly ${M}_{z}\propto \sqrt{M}$ and Z_planet $\propto \,1/\sqrt{M}$. Our parameter uncertainties exclude a flat line by a wide margin, but the distribution has a fair amount of spread around our fit. The intrinsic spread was the factor 10^σ = 1.82 ± .09 (because σ was calculated on the log of the variables), which means that 1σ of the data is within a factor of 1.82 of the mean line. While some of this may be from observational uncertainty, it seems likely that other effects, such as the planet's migration history and the stochastic nature of planet formation, also play a role. With this in mind, using planet mass alone to estimate the total heavy element mass appears accurate to a factor of a few.

Zoom In Zoom Out Reset image size

5.2. Effect of Stellar Metallicity

The metallicity of a star directly impacts the metal content of its protoplanetary disk, increasing the speed and magnitude of heavy element accretion. We examined our data for evidence of this connection. MF2011 observed a correlation for high metallicity parent stars between [Fe/H] and the heavy element masses of their planets (see also Guillot et al. 2006 and Burrows et al. 2007 for similar results from inflated planets). If we constrain ourselves to the 14 planets in MF2011, we see the same result. However, the relation becomes somewhat murky for our full set of planets (see Figure 8). Applying Kendall's Tau to the most likely values of metal mass and [Fe/H], we measure a correlation of 0.08845 and a p-value of 0.3805, which indicates no correlation. Some of the reason for this may lie with the high observational uncertainty in our values for stellar metallicity, but it is still difficult to believe that there is a direct power-law relationship.

Zoom In Zoom Out Reset image size

Transit surveys should not be biased in stellar metallicity, so we can instead consider how the distribution of planet mass and metallicity vary as a function of stellar metallicity. Most of the planets with heavy element masses above 100 M_⊕ orbit metal-rich stars; there is no clear pattern for planets with lower metal masses. Considering the connection between planet mass and heavy element mass, we note that planets more massive than two to three Jupiters are found far less often around low-metallicity stars (see Figure 9). Presumably, these trends are connected. This is similar to one of the findings of Fischer & Valenti (2005) in which the number of giant planets and the total detected planetary mass are correlated with stellar metallicity. The population synthesis models in Mordasini et al. (2012) also observe and discuss an absence of very massive planets around metal-poor stars. In the future, a more thorough look at this connection should take into account stellar metal abundances other than iron and put an emphasis on handling the high uncertainties in measurements of stellar metallicity.

Zoom In Zoom Out Reset image size

5.3. Metal Enrichment

A negative correlation between a planet's metal enrichment relative to its parent star was suggested in MF2011 and found in subsequent population synthesis models (Mordasini et al. 2014), so we revisited the pattern with our larger sample. We see a good relation in our data as well (Figure 11), and using Kendall's tau as before, we find a correlation of 0.4398 with a p-value of 1.3 × 10⁻⁵. The exponent for the fit shown in Figure 11 (−0.45 ± 0.09) differs somewhat from the planet formation models of Mordasini et al. (2014; between −0.68 and −0.88). The pattern appears to be stronger than if we considered only the planetary metal fraction Z_planet alone (shown in Figure 10). This supports the notion that stellar metallicity still has some connection to planetary metallicity, even though we do not observe a power-law-type of relation. Jupiter and Saturn, shown in blue, fit nicely in the distribution. Our results show that even fairly massive planets are enriched relative to their parent stars. This is intriguing, because it suggests that (since their cores are probably not especially massive) the envelopes of these planets are strongly metal-enriched, a result that can be further tested through spectroscopy. Note that we calculate our values of Z_star by assuming that stellar metal scales with the measured iron metallicity [Fe/H], using the simple approximation ${Z}_{\mathrm{star}}=0.014\times {10}^{[\mathrm{Fe}/{\rm{H}}]}$ (given our [Fe/H] uncertainties, a more advanced treatment would not be worthwhile). Considering other measurements of stellar metals in the future would be illuminating. For instance, since oxygen in a dominant component of both water and rock, perhaps there exists a tighter correspondence between Z_planet and the abundance of stellar oxygen, rather than with stellar iron.

Zoom In Zoom Out Reset image size

We also considered the possibility that orbital properties might relate to the metal content, perhaps as a proxy for the migration history. We plot the residual from our mass versus Z_planet/Z_star fit against the semimajor axis, period, eccentricity, and parent star mass in Figure 12, and the same residuals against our mass vs metallicity fit in Figure 13. No pattern is evident for any of these. Given the number of planets and the size of our error bars in our sample, we cannot rule out that any such patterns exist, but we do not observe them here. We also considered the residual against the stellar flux. Any giant planets with radius inflation that made it into our sample would appear as strong outliers below the fit, since we would have mistaken inflation for lower heavy element masses. Therefore, the lack of a pattern here suggests our flux cut is eliminating the inflated hot Jupiters as intended.

5.4. Heavy Element Masses in Massive Planets

One might expect that core accretion produces giant planets with total metal masses of approximately M_z = M_core + Z_* M_env, where the core mass M_core ∼ 10 M_⊕ depends on atmospheric opacity but is not known to depend strongly on final planet mass, a large fraction of disk solids are assumed to remain entrained with the envelope of mass M_env as it is accreted by the core, and the mass fraction of the disk in metals is assumed to match that of the star, Z_*. Such a model would predict that a planet of total mass M_p = M_core + M_env has a metallicity ${Z}_{{pl}}\equiv {M}_{z}/{M}_{{\rm{p}}}$ that is related to that of the star by ${Z}_{{pl}}/{Z}_{* }=1+({M}_{\mathrm{core}}/{M}_{{\rm{p}}})(1-{Z}_{* })/{Z}_{* }$. This expression falls off substantially more rapidly with planet mass than the ${Z}_{{pl}}/{Z}_{* }\approx 10{({M}_{{\rm{p}}}/{M}_{{\rm{J}}})}^{-0.5}$ fit in Figure 11.

Zoom In Zoom Out Reset image size

This lack of good agreement is not surprising given that the above model does not predict the high metallicities of solar system giants. These have long been interpreted as coming from late-stage accretion of additional planetesimal debris (e.g., Mousis et al. 2009). In keeping with this solar system intuition, we instead propose that the metallicity of a giant planet is determined by the isolation zone from which the planet can accrete solid material. We assume that a majority of solids—which we treat interchangeably with metals in this initial investigation—eventually decouple from the disk gas and can then be accreted from the full gravitational zone of influence of the planet. An object of mass M, accreting disk material with surface density Σ_a at distance r from a star of mass M_*, can accumulate a mass

$$\begin{eqnarray}&&{M}_{a}=2\pi r(2{f}_{{\rm{H}}}{R}_{{\rm{H}}}){{\rm{\Sigma }}}_{a}=4\pi {f}_{{\rm{H}}}{\left(\displaystyle \frac{M}{3{M}_{* }}\right)}^{1/3}{{\rm{\Sigma }}}_{a}{r}^{2},\end{eqnarray} \tag{ 5 }$$

where f_H ∼ 3.5 is the approximate number of Hill radii ${R}_{{\rm{H}}}=r{(M/3{M}_{* })}^{1/3}$ from which accretion is possible as long as the orbital eccentricity of accreted material is initially less than ${(M/3{M}_{* })}^{1/3}$ (Lissauer 1993). Thus, under our assumption that solids have decoupled from the gas, a planet of mass M_p can accrete a total solid mass of

$$\begin{eqnarray}&&{M}_{z}\approx 4\pi {f}_{{\rm{H}}}{f}_{e}{Z}_{* }{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{3{M}_{* }}\right)}^{1/3}{\rm{\Sigma }}{r}^{2}\,,\end{eqnarray} \tag{ 6 }$$

where Σ is the total surface density of the disk and ${f}_{e}{Z}_{* }{\rm{\Sigma }}$ is the surface density in solids. The parameter f_e allows for an enhancement in the metal mass fraction of the disk compared to the solar value, Z_*, for example due to radial drift of solid planetesimals through the gas nebula. No enhancement corresponds to f_e = 1. We note that Equation (6) applies independent of the solid mass fraction of the planet because M_p is taken to be the final observed planet mass, including any accreted solids.

For comparison, the standard isolation mass of a planet forming in a disk with total surface density Σ is

$$\begin{eqnarray}&&{M}_{\mathrm{iso}}={[4\pi {f}_{{\rm{H}}}{(3{M}_{* })}^{-1/3}{\rm{\Sigma }}{r}^{2}]}^{3/2}\end{eqnarray} \tag{ 7 }$$

which may be calculated using Equation (5) with ${M}_{a}=M={M}_{\mathrm{iso}}$ and Σ_a = Σ. Recalling that ${Z}_{{pl}}={M}_{z}/{M}_{{\rm{p}}}$, Equations (6) and (7) combine to yield

$$\begin{eqnarray}&&\displaystyle \frac{{Z}_{{pl}}}{{Z}_{* }}={f}_{e}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{\mathrm{iso}}}\right)}^{-2/3}\end{eqnarray} \tag{ 8 }$$

For f_e = 1, when M_p = M_iso, the total mass of the planet equals the total disk mass in its isolation region and ${Z}_{{pl}}/{Z}_{* }=1$, as, for example, would happen if an isolation-mass planet formed by accumulating all material within its isolation zone. For M_p < M_iso, the planet's metallicity exceeds the metallicity of the star because the planet has not been able to accrete all isolation-zone gas but we assume that it is able to accrete all of the region's solids.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Equation (8) encapsulates the physics of our model, but to compare with observed planets, we re-express this result in terms more easily related to the expected population of protoplanetary disks from which planets form—Toomre's Q parameter (Safronov 1960; Toomre 1964) and the disk aspect ratio H/r. In a Keplerian disk, $Q={c}_{s}{\rm{\Omega }}/(\pi G{\rm{\Sigma }})=(H/r){(\pi {r}^{2}{\rm{\Sigma }}/{M}_{* })}^{-1}$, where c_s is the isothermal sound speed, Ω is the orbital angular velocity, H = c_s/Ω is the disk scale height, and G is the gravitational constant. Thus,

$$\begin{eqnarray}&&\displaystyle \frac{{Z}_{{pl}}}{{Z}_{* }}\approx 3{f}_{{\rm{H}}}{f}_{e}\displaystyle \frac{H}{r}{Q}^{-1}{\left(\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{* }}\right)}^{-2/3}\end{eqnarray} \tag{ 9 }$$

Figure 14 replots the values and best-fit line of Figure 11, overlaying Z_pl/Z_* calculated using Equation (9) with f_e = 1 for each planet's M_p and M_*. We use a fiducial value of H/r = 0.04 (corresponding to, e.g., 2 au at 200 K) and plot curves for Q = 1, 5, and 20. The Q = 5 curve provides a good match for the best-fit line, while Q = 1 and Q = 20 bound the remaining data points. We omit data points corresponding to planet masses in excess of the total local disk mass, ${M}_{{\rm{p}}}\gt \pi {r}^{2}{\rm{\Sigma }}=(H/r){Q}^{-1}{M}_{* }$, which removes the portion of the Q = 20 curve corresponding to the highest planet masses. For f_e > 1, the values of Q plotted in Figure 14 should be multiplied by f_e. Figure 15 provides an example for f_e = 2, demonstrating that for modest enhancements in the solid to gas ratio in the disk, limits on the total disk mass could preclude formation of the most massive planets in the sample in the highest Q disks, potentially explaining the lack of points in the bottom-right portion of the plot. Values of Q < 1 imply gravitational instability and cannot be maintained for extended periods in a disk, so the fact that our modeled Z_pl/Z_* does not require Q < 1 for any value of f_e is encouraging. For reference, at Jupiter's location r = 5 au in the minimum mass solar nebula Σ = 2 × 10³g cm⁻² ${(a/\mathrm{au})}^{-3/2}$ (Hayashi 1981), Q ≈ 25.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Why might a planet remain smaller than its isolation mass? For our model, the reason is irrelevant. The planet could have accreted its envelope at a different location in the disk at an earlier time when solids were unable to decouple from the gas. It could even be a remnant fragment formed by gravitational instability. Whatever its formation process, our model only asserts that it accretes most of its metals at some point concurrent with or after it has accumulated most of the gas it will accrete from the nebula. Hence, its metallicity is determined by the mass in solids contained within its gravitational feeding zone. However, for a plausible scenario, we may again appeal to studies of Jupiter formation, which suggest that the planet's final mass is determined by the mass at which the planet truncates gas accretion by opening a large gap in the disk (e.g., Lissauer et al. 2009).

Gap opening and the planetary starvation that results remains a topic of continuing research (e.g., Crida et al. 2006; Fung et al. 2014; Duffell & Dong 2015). Here we simply note that the planet masses observed in our sample are consistent with a classic theory of gap starvation. Tidal torques opposed by viscous accretion yield a gap width

$$\begin{eqnarray}&&{\rm{\Delta }}={\left(\displaystyle \frac{{f}_{g}}{\pi \alpha }\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{* }}\displaystyle \frac{{r}^{2}}{{H}^{2}}\right)}^{1/3}{R}_{{\rm{H}}},\end{eqnarray} \tag{ 10 }$$

where f_g ≈ 0.23 is a geometric factor (Lin & Papaloizou 1993). Assuming that gas accretion through the gap becomes inefficient for Δ/R_H = f_S ∼ 5 (Lissauer et al. 2009; Kratter et al. 2010), a planet mass M_p set by gap truncation implies the disk viscosity parameter (Shakura & Sunyaev 1973; Armitage 2011)

$$\begin{eqnarray}&&\alpha =\displaystyle \frac{{f}_{g}}{\pi {f}_{S}^{3}}\displaystyle \frac{{M}_{{\rm{p}}}}{{M}_{* }}{\left(\displaystyle \frac{H}{r}\right)}^{-2}\,.\end{eqnarray} \tag{ 11 }$$

These values—displayed in Figure 16 for the same planets plotted in Figure 14—span a reasonable theoretical range for protoplanetary disks, particularly in dead zones, where most giant planets are thought to form (e.g., Turner et al. 2014; Bai 2016).

Zoom In Zoom Out Reset image size

We note that Equations (9) and (11) do not depend directly on the planet's distance from its star, r. Instead, they are functions of H/r, which depends weakly on r for typical disks, and of Q, which may be thought of as a parametrization of the disk mass. Hence, our model applies regardless of whether the observed planets have migrated from their formation locations. For relatively low-mass disks, such migration is in fact required. For the minimum mass solar nebula used above, for example, ${M}_{\mathrm{iso}}=0.6{M}_{{\rm{J}}}{(r/\mathrm{au})}^{3/4}$. Giant planets separated by less than 1 au from their host stars most likely did migrate from more distant formation locations (e.g., Dawson & Murray-Clay 2013).

There is a strong connection between planet mass and metal content. From a sample of 47 transiting gas giants, we find that the heavy-element mass increases as $\sqrt{M}$, so Z_planet decreases as $1/\sqrt{M}$. We also see that our planets are consistently enriched relative to their parent stars, and that they likely all have more than a few M_⊕ of heavy-elements. These results all support the core-accretion model of planet formation and the previous results from MF2011 that metal-enrichment is a defining characteristic of giant planets. We have also shown that our results for Z_planet/Z_star are comfortably consistent with a simple planet formation model using plausible values for disk parameters. Our results were not consistent with a more naive model of formation in which a fixed-mass core of heavy-elements directly accretes parent-star composition nebular material.

This work suggests that spectroscopy of the atmospheres of gas giants should also yield metal-enrichments compared to parent star abundances (Fortney et al. 2008), as is seen in the solar system. We suggest that this growing group of <1000 K planets should be a sample of great interest for atmospheric observations with the James Webb Space Telescope. Since the bulk metallicity of the cool planets can be determined, atmospheric studies to determine metal-enrichments can be validated. For most of our planets more massive than Saturn, the heavy elements' values are high enough that most metals are in the envelope, rather than the core (See Figure 7), so our Z_planet values are upper limits on Z_atmosphere. Atmospheric observations to retrieve the mixing ratios of abundant molecules like H₂O, CO, CO₂, and CH₄ could show if most heavy elements are found within giant planet envelopes or within cores. In comparison, for hot Jupiters >1000 K, we cannot estimate with confidence the planetary bulk metal enrichment since the radius inflation power is unknown.

Some connection between planetary heavy element mass and stellar metallicity is suggested by our work, but the correlation is not strong. A fruitful area of future investigation will be in analyzing stellar abundances other than iron. Learning which stellar metals most strongly correlate with planetary bulk metallicity (for instance, Fe, Mg, Si, Ni, O, and C) would hint at the composition of planetary heavy elements and could provide insights into the planet formation process. With the continued success of ground-based transit surveys, along with K2 (Howell et al. 2014), and the 2017 launch of TESS, we expect to see many more "cooler" gas giant planets amenable to this type of analysis, which will continue to provide an excellent opportunity for further exploring these relations. With this continuing work on bulk metal-enrichment of planets, and the spectroscopy of planetary atmospheres, the move toward understanding planet formation in the mass/semimajor axes/composition planes should be extremely fruitful.

The authors thank Johanna Teske and Natalie Hinkel for insightful discussions during the scope of this work. J. J. F. acknowledges the support of NASA XRP grant NNX16AB49G.