Chemical Diversity in Protoplanetary Disks and Its Impact on the Formation History of Giant Planets

The protoplanetary disk considered in Paper I and in this work is modeled over the observed disk HD 163296 (Isella et al. 2016; Turrini et al. 2019). The gas surface density and temperature profiles of HD 163296 were rescaled to match the total disk mass and the estimated temperature profile of the Minimum Mass Solar Nebula (MMSN; Hayashi 1981). The physical parameters of the resulting disk are summarized in Table 1. Our disk model assumes radial profiles of gas surface density and temperature on the midplane that are constant in time and parameterized as:

$$\begin{eqnarray}&&{{\rm{\Sigma }}}_{\mathrm{gas}}(r)={{\rm{\Sigma }}}_{0}{\left(\displaystyle \frac{r}{165\,\mathrm{au}}\right)}^{-0.8}\exp \left[-{\left(\displaystyle \frac{r}{165\,\mathrm{au}}\right)}^{1.2}\right]\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&T(r)={T}_{0}{\left(\frac{r}{1\,{\rm{a}}{\rm{u}}}\right)}^{-1/2},\end{eqnarray} \tag{ 2 }$$

where Σ₀ = 3.3835 g cm⁻² and T₀ = 280 K.

Regarding the solid component of the disk, we assumed that it is distributed between 1 and 150 au with a surface density profile:

$$\begin{eqnarray}{{\rm{\Sigma }}}_{\mathrm{solids}}(r)=\left\{\begin{array}{ll}2\,{Z}_{i}(r){{\rm{\Sigma }}}_{\mathrm{gas}}(r) & \quad x\leqslant 150\,\mathrm{au}\\ 0 & \quad x\gt 150\,\mathrm{au},\end{array}\right.\end{eqnarray} \tag{ 3 }$$

where Z_i (r) is the mass fraction of condensed material. A factor of 2 was introduced as a concentration factor to account for the inward drift of dust and pebbles as a consequence of their dynamical coupling with the gas (see Paper I and references therein for further discussion).

Following the approach of Paper I, we assumed that the bulk of dust in the disk is rapidly converted into planetesimals over a timescale of 1 Myr, as suggested by comparisons between the masses of exoplanetary systems and of the dust and gas in protoplanetary disks (e.g., Manara et al. 2018; Mulders et al. 2021). The assumption is also consistent with recent observations of dust temporal evolution in disks (see Testi et al. 2022 and Bernabò et al. (2022) for further discussion) and with meteoritic data from the solar system (e.g., Scott 2007; Lichtenberg et al. 2022, and references therein).

The conversion of dust into planetesimals reduces the efficiency of the gas–grain chemistry, hence slowing down the chemical evolution of the disk itself. When the disk is severely depleted in dust with respect to the initial stages of its evolution, its chemical composition can be reasonably approximated as fixed. We then assumed that the composition of our disk evolves until 1 Myr and remains fixed thereafter (see Section 2.4 for details on the compositional model). The implications of physically and chemically evolving disks will be explored in future works.

The planetesimals in our model are all characterized by a fixed radius of r_p = 50 km (Klahr & Schreiber 2016; Johansen & Lambrechts 2017) and are divided into two populations, depending on whether they are located within or beyond the water snowline. Planetesimals inside the water snowline are rock dominated and characterized by a density of ρ_rock = 2.4 g cm⁻³, while those beyond the water snowline are enriched in ices and characterized by a density of ρ_ice = 1 g cm⁻³. See the Appendix for further discussion.

The radius and density values are used to compute the effects of gas on the planetesimal dynamics. All the planetesimals evolve dynamically under the influence of both the forming giant planet and the disk itself. Specifically, the planetesimals interact gravitationally with the forming giant planet, whose collisional cross section determines whether planetesimals are accreted or scattered by planetary encounters. Moreover, at each location in the disk the planetesimals are subject to two competing forces: the dynamical excitation due to the disk self-gravity and the damping effect of the gas drag. The implementation of the disk self-gravity is based on the analytical treatment for thin disks by Ward (1981), following the approach of Marzari (2018) and Nagasawa et al. (2019). The effect of the gas drag on the dynamical evolution of the planetesimals is modeled following the treatment by Brasser et al. (2007), with the updated drag coefficients from Nagasawa et al. (2019). See Paper I and references therein for more details.

The planetesimal disk is simulated by means of a set of dynamical tracers distributed randomly between 1 and 150 au with a uniform probability distribution and a spatial density of 2000 tracers/au. Each tracer represents a swarm of planetesimals. The mass of the swarm is computed by dividing the total solid mass in an annular region of the disk (computed by integrating Equation (3)) by the number of dynamical tracers it contains. One can then associate the flux of impacting tracers recorded by the simulations with a mass flux of accreted planetesimals on the giant planet (see Paper I for details).

In Figure 1, the first four panels from the top are snapshots of the N-body simulations showing the dynamical evolution of the planetesimals at different stages of the planet formation and migration history. As the planet forms and migrates, planetesimals from different compositional regions of the disk are dynamically excited and end up being scattered or accreted onto the planet. An example of the normalized flux of accreted planetesimals is shown by the blue curve in the bottom-right panel of Figure 1. See Paper I and the Appendix for more details.

The N-body simulations were performed with Mercury-Ar χ es (Turrini et al. 2019, 2021), a high-performance implementation of the hybrid symplectic algorithm of the Mercury 6 software from Chambers (1999). Besides the improvements in numerical stability and computational efficiency with respect to Mercury, the algorithm of Mercury-Ar χ es allows for simulating the mass growth, the radius evolution, and the orbital migration at each stage of giant planet formation, as well as the effects of the disk self-gravity and of the gas drag. We refer the reader to Paper I and the Appendix for more details on the theoretical treatment of these processes and on their numerical implementation in the N-body simulations. In the following, we provide a brief overview of what is included in the planet formation and migration model.

We simulated a giant planet that forms and migrates on the disk midplane (i.e., on an orbit with inclination i = 0) over a timescale of 3 Myr, growing from a planetary embryo to a Jupiter-like giant planet. We adopted a two-phase approach to model both the mass growth and the evolution of the physical radius of the planet, following the growth tracks from Lissauer et al. (2009); Bitsch et al. (2015), and D'Angelo et al. (2021) using the parametric approach from Turrini et al. (2011, 2019). In the first 2 Myr, the planet accretes its core and extended atmosphere and grows from an initial mass of M₀ = 0.1 M_⊕ (mass of the planetary embryo) to a critical mass of M_c = 30 M_⊕, equally shared between the core and the atmosphere. In this phase, the physical radius of the planet grows following the approach described by Fortier et al. (2013), which is based on the hydrodynamical simulations by Lissauer et al. (2009). Over the last 1 Myr, the planet undergoes its runaway gas accretion phase and its mass grows until the final value of M_F = 317.8 M_⊕ (equal to 1 Jovian mass) is reached. At the onset of the runaway gas accretion, the physical radius of the planet starts to shrink due to the gravitational infall of the gas, reaching a final value of R_I = 1.15 × 10⁵ km (equal to 1.6 Jovian radii). An example of the mass growth track normalized to its final value is shown in the bottom-right panel of Figure 1 (orange curve) for the case of a giant planet that starts forming at 19 au.

To model the migration of the giant planet we adopted the realistic nonisothermal migration tracks from the population synthesis models by Mordasini et al. (2015) following a piecewise approach based on the analytical treatments of Hahn & Malhotra (2005) and Walsh et al. (2011). The simulations assume that the protoplanet initially undergoes a damped Type I migration. During the growth of the planetary core from M₀ to the critical value M_C , the protoplanet proceeds through a linear regime of slow Type I migration that accounts for 40% of the total radial displacement of the giant planet. When the critical mass is reached and the runaway gas accretion phase begins, the protoplanet enters a faster regime of Type I migration. Once the protoplanet becomes massive enough to open a gap in the disk, it transitions to the slower Type II migration regime. These two last phases take the form of a power-law migration regime that accounts for the remaining 60% of the total radial displacement of the giant planet. The green curve in the bottom-right panel of Figure 1 shows an example of the migration track normalized to its initial value for the case of a giant planet that starts forming at 19 au.

The simulations model a total of six migration scenarios, with the protoplanet starting at 5, 12, 19, 50, 100, and 130 au from the star and ending at 0.4 au. The initial positions span the ranges of the observed architectures of giant planets in the solar system (5–10 au), in exoplanetary systems (from fractions of au to 20 au), and in circumstellar disks (100–150 au). The choice of a wide range of distances is meant to investigate the compositional implications of such diverse formation regions.

As part of our exploration of the parameter space, we investigated how the final composition of the planet changes as we vary the level of chemical activity in the disk. To this aim, we took advantage of the results by Eistrup et al. (2016), who modeled the ionization environment of the protoplanetary disk as set by two key sources. Specifically, they considered the ionization from the decay of short-lived radionuclides (SLRs) and from cosmic rays (CRs). For both disk chemical setups they analyzed (inheritance and reset, see Section 2.4), they explored a case of low ionization level, in which SLRs are the only source of ionization, and a case of high ionization level, in which an additional contribution from CRs of external origin is also taken into account. This results in a total of four disk chemical scenarios that we used to derive realistic planetary compositions from the outcomes of the simulations of Paper I. Future works will address the dependency of our results on a more detailed treatment of the interaction between the disk and the sources of ionization (Padovani et al. 2016, 2018; Rodgers-Lee et al. 2020).

In the low ionization scenario, the dominant contribution to ionization comes from the decay products of ²⁶Al, ³⁶Cl, and ⁶⁰Fe, which have half-lives t_half of 0.74, 0.30, and 2.6 Myr, respectively (Cleeves et al. 2014). Eistrup et al. (2016) adopted a simplified version of the analytical prescription given in Equation (30) by Cleeves et al. (2013b) for the ionization rate per H₂ molecule at the disk midplane:

$$\begin{eqnarray}&&{\zeta }_{\mathrm{SLR}}(r)=(1.25\times {10}^{-19}\,{{\rm{s}}}^{-1}){\left(\displaystyle \frac{{\rm{\Sigma }}(r)}{{\rm{g}}\,{\mathrm{cm}}^{-2}}\right)}^{0.27},\end{eqnarray} \tag{ 4 }$$

where Σ(r) is the surface density of the disk as a function of the radius r. Equation (4) does not account for the time decay of the ionization rate. Therefore, in the SLRs-dominated environment by Eistrup et al. (2016), the ionization rate is constant in time and higher in the inner and denser regions of the midplane.

For the high ionization scenario, in addition to SLRs, Eistrup et al. (2016) included a contribution of cosmic rays originating from outside the system. They parametrically modeled the CR ionization rate per H₂ molecule as:

$$\begin{eqnarray}&&{\zeta }_{\mathrm{CR}}(r)=(1\times {10}^{-17}\,{{\rm{s}}}^{-1})\cdot \exp \left(\displaystyle \frac{-{\rm{\Sigma }}(r)}{96\,{\rm{g}}\,{\mathrm{cm}}^{-2}}\right),\end{eqnarray} \tag{ 5 }$$

where the characteristic value of ζ_CR ∼ 10⁻¹⁷ s⁻¹ is the interstellar rate that is typically assumed in models of disk chemistry (Webber 1998; Cleeves et al. 2014). The addition of CRs allows for higher ionization in the outer disk, where the surface density is lower and CRs can more efficiently penetrate into the midplane.

The aim of considering the two scenarios is to investigate the impact of ionization-driven chemical activity on the chemical composition of the disk midplane. Both the prescriptions given by Equations (4) and (5) should be interpreted as conservative estimates of the rates that regulate the ionization environment of a typical midplane. A recent study by Padovani et al. (2018) suggests that the CR ionization rate in high-density environments is higher than the typically assumed value of 10⁻¹⁷ s⁻¹. Moreover, the model by Eistrup et al. (2016) includes only attenuation due to the disk's surface density. Works by Cleeves et al. (2013a, 2013b, 2014) revealed that the CR flux at the disk's surface may be strongly attenuated by winds and/ or magnetic fields. In particular, modulation by stellar winds may reduce the CR flux reaching the midplane by many orders of magnitude (ζ_CR ≲ 10⁻²⁰ s⁻¹), leaving SLRs as the dominant midplane ionization source (Cleeves et al. 2013a). Nevertheless, as mentioned before, the abundances of SLRs do evolve with time. As pointed out by Cleeves et al. (2013b), the SLRs ionization rate scales with the inverse of t_half and long-lived radionuclides (e.g., ⁴⁰K) only produce ionization rates of the order of ∼10⁻²² s⁻¹ or even less. Moreover, although there is evidence of an enhanced abundance of SLRs in the early solar system, little is still known about their actual contribution to the total ionization rate in the midplane of a typical disk (Cleeves et al. 2013b).

Regarding other sources of ionization, Cleeves et al. (2013a) showed that photoionization from stellar and interstellar UV acts largely on C-bearing species in the upper atmosphere of the disk, leaving the midplane essentially unaffected. On the other hand, X-ray photons from the star penetrate deeper and dominate the intermediate layers of the disk. Only the scattered component of stellar X-rays is able to reach the midplane (Ercolano & Glassgold 2013). However, for a stellar X-ray luminosity of L_X = 10^29.5 erg s⁻¹, the scattered X-rays are expected to produce ionization rates falling in the range ζ_CR ∼ (1–10) × 10⁻²¹ s⁻¹, far below the rates associated with SLRs and CRs, with minimal impact on the midplane chemistry (see Cleeves et al. 2014, and Figure 2 and 3 therein). Concerning X-ray background fields, as in the case of a protoplanetary disk embedded in a cluster, recent works revealed that these would have relatively little impact on disk chemistry (Meijerink et al. 2012; Rab et al. 2018). In particular, for an ordinary disk in a typical low-mass star-forming region, the ionization rate produced by an X-ray background flux of 2 × 10⁻⁴ erg cm⁻² s⁻¹ would equal the interstellar ζ_CR only in the outer disk, at r ∼ 200 au from the central star (see Rab et al. 2018, and Figure 10 therein).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To characterize the chemical environment in which the giant planet forms and migrates, we used an updated version of the compositional model presented in Paper I. The model quantifies the composition of the disk over its key components: rocks, organics, ices, and gas. In this work, we introduce a more realistic abundance profile for the refractory organic carbon, and we revise the treatment of the volatile component to allow for different evolutionary scenarios and ionization levels. Here, we outline the key features of the model, and we provide details on its refinements. We refer the reader to Paper I for a complete discussion.

We adopted the protosolar composition for both the star and its hosting protoplanetary disk. The protosolar elemental abundances, which characterize the original mixture of gas, were taken from Asplund et al. (2009) and Scott et al. (2015a, 2015b). Overall, they result in mass fractions X = 0.7148, Y = 0.2711, and Z = 0.0141, of H, He, and heavy elements, respectively. The gas composition sets the initial conditions for the condensation sequence across the disk, which regulates the distribution of the elements among the different phases (gas and solid) and carriers (rocks, organics, gas, and ices). As in Paper I, we focus on four tracing elements: C, O, N, and S. Their partitioning between the three solid components and between the solid and gas phases is derived from solar system data and recent results from astrochemical models, as discussed below.

2.4.1. Refractory Elements and Rocks

2.4.2. Volatiles

2.4.3. Refractory Organics

The final composition of giant planets is generally expected to reflect the chemical structure of the protoplanetary disk, especially that of the region of the disk where the planet was born (e.g., Öberg et al. 2011; Lothringer et al. 2021). We can assess the degree to which this assumption is reasonable by comparing the elemental ratios in the two environments. The results of such analysis are illustrated in Figures 4–7. Each figure examines one of the four chemical scenarios. The top two panels show the normalized elemental ratios in the atmosphere of gas-dominated (on the left) and solid-enriched (on the right) giant planets at the end of the six simulations. The lines connecting the different scenarios are provided as a visual aid only. The bottom four panels show the radial profiles of the normalized elemental ratios in the disk, as computed from Figure 3. Such profiles can be seen as snapshots of the midplane chemical structure, where solid and dashed lines characterize the solid and gas phase, respectively.

Overall, the four figures show that the atmosphere of the giant planet inherits only part of the chemical properties of the disk. Specifically, elemental ratios that are sub- or superstellar in the disk will generally be so in the planet as well. However, the results also show that this disk-planet link cannot be used to make quantitative predictions. Specifically, the plots reveal that for some migration scenarios the normalized elemental ratios in the planet and in the disk do not immediately relate. For instance, although the S/N* is markedly superstellar throughout the disk and increases toward the star, its value in the solid-enriched giant planet that starts forming at 5 au drops to the substellar level. The final composition of the planet envelope depends crucially on how efficiently the gas and the solids are accreted from each region of the disk along the migration pathway. For gas-dominated giant planets, the accretion of gas reaches its maximum during the runaway gas accretion phase (see Section 2.2). The final atmosphere of the planet will therefore be strongly influenced by the chemical properties of the region of the disk traversed by the planet during this phase. In solid-enriched giant planets the accretion of gas and solids are coupled, but their relative contribution to the planetary metallicity depends on the extent of the migration. In particular, the larger the migration the larger the number of planetesimals encountered by the planet and the larger the contribution of planetesimal accretion to the final metallicity. As the planet migrates for shorter distances, it encounters and accretes fewer planetesimals, and the contributions to its metallicity by gas and solids become comparable. As a consequence, the elemental ratios of solid-enriched giant planets follow the global trends of the disk solid component only in large-scale migration scenarios, while diverging from them for limited migrations.

The decreasing trend of the planetary S/N* ratio from large- to short-scale migrations can be easily understood considering how S and N are distributed in the disk and how they are accreted on the planet. In our disk, S is entirely locked in the solid phase, while N is almost entirely in the gas phase as N₂, slightly deviating from the stellar value beyond the NH₃ snowline (see Figure 3). As a direct consequence of this partitioning, in the planet envelope the total S scales linearly with the accreted planetesimals, while the total N is accreted in stellar proportion almost independently on the radial migration. Consequently, even though the S/N* ratio of solids in the inner disk is about 40 times higher than the stellar value, for short-scale migrations the contribution of the accreted N becomes large enough to bring the S/N* ratio in the planet to substellar value. Similar arguments can be used to explain why the N/O*, C/O*, and C/N* ratios of solid-enriched giant planets approach stellar values for short-scale migrations.

Our results emphasize how the composition of giant planets and their atmospheres can be connected to that of their native disks only by coupling the disk chemical structure with the growth and migration tracks that characterize planet formation. With this in mind, we can proceed with the discussion of our results and their implications for planet formation.

Figure 8 compares the normalized elemental ratios of solid- and gas-dominated giant planets in the inheritance (blue points) and reset (red points) scenarios. The results are shown for each of the six simulations, identified by the location of the planet at the onset of its formation. The C/O* ratio shows almost flat trends, revealing a rather limited diagnostic power.

For short-scale migrations, the C/O* ratios in all scenarios vary between 1 and 2 times the stellar value. The limited variations between the different disk chemistry and planet migration scenarios are too small to be unequivocally resolved with the current accuracy of retrieval methods, which is of the order of 20% (see Barstow et al. 2020 and Paper I for further discussion). Therefore, giant planets that formed close to the star in a disk that inherited its composition from the prestellar core cannot be observationally distinguished from those that formed in a disk that experienced a complete reset of the chemistry, nor can we resolve whether they are gas dominated or solid enriched.

The degeneracy is partially broken for large-scale migrations, where the C/O* ratio of solid-enriched and gas-dominated giant planets follows distinct trends. Specifically, solid-enriched giant planets are characterized by stellar values of the C/O* ratio, independent of how far the planet started its migration and the chemical structure of the birth environment. On the contrary, gas-dominated giant planets are characterized by a C/O* ratio significantly deviating from the stellar value. In this case, the C/O* ratio provides constraints on both the migration and the disk chemical scenario. Detailed discussion on each of these individual aspects is provided in Sections 3.3, 3.4, and 3.5.

As anticipated in Section 3.2, the C/O* ratio provides information on the accretion history of giant planets only for large-scale migration scenarios. In this case, stellar values are associated with solid-enriched giant planets, whereas marked sub- and superstellar values provide an indication of accretion dominated by gas. Specifically, substellar C/O* ratios are associated with gas-dominated giant planets in the reset scenarios, whereas superstellar values characterize gas-dominated giant planets in the inheritance scenarios.

To understand the origin of such trends we refer to the disk composition, as described by the molecular abundance profiles in Figure 2. Beyond 10 au, the gas phase in the reset scenarios is characterized by a lower abundance of CH₄ and higher abundances of CO and O₂ with respect to the inheritance scenarios. Because both CO and O₂ typically condense at very low temperatures, they act as reservoirs of gaseous O in most of the disk extension. Consequently, the abundance of total O* in the gas phase is higher in the reset scenarios than in the inheritance ones. In particular, beyond 10 au, the abundance of O* exceeds that of C*. Such behavior is shown in Figure 3 and is the reason why the C/O* ratio of the disk gas phase drops below the stellar value in the reset scenarios (see Figures 6 and 7). As higher abundances in the gas correspond to lower abundances in the solids, the opposite trend is observed for the C/O* ratios of the disk solid phase. Moving from the disk to the planet, we recall that the final elemental ratios in the planet envelope depend on the mass accretion rate (see the discussion in Section 3.1). Such rate is a nonlinear function of the orbital distance, and it is highly influenced by the position on the disk where the runaway gas accretion occurs. As discussed in Section 2.2, in our model the runaway gas accretion phase begins once the first 40% of the total radial displacement of the planet has been covered. For migrations starting beyond 20 au, in our simulations, this means that the region traversed by the planet during the runaway gas accretion phase falls between 10 au and 80 au. Therefore, gas-dominated giant planets in the reset scenarios accrete most of their gas from the disk region in which the C/O* ratio is substellar, which is why their final C/O* ratio is also substellar. When it comes to solid-enriched giant planets, the contribution of planetesimal accretion to the envelope metallicity increases with the extent of migration. For large-scale migrations, the metallicity is dominated by the planetesimals accreted from the disk region in which the C/O* ratio of the solids is superstellar, which is why the final C/O* ratio in the planet envelope is also superstellar.

We highlight that in the inheritance scenario with a high ionization level, the C/O* ratio of gas-dominated giant planets slightly decreases with increasing length of migration. For the accuracy of the current retrieval methods (see Section 3.2), this introduces a degeneracy of the planet's accretion history. Specifically, for very large-scale migration scenarios, gas-dominated giant planets in the inheritance high scenario would be indistinguishable from solid-enriched ones. Such a decreasing trend in the inheritance high scenario originates from the enhanced abundance of O* between 50 au and 100 au in our disk model (see Figure 3), which in turn results from the peak abundance of O₂ in the outer disk (see Figure 2). In this region, the abundance of O* exceeds that of C*, hence bringing the disk C/O* ratio to substellar values (see Figure 5). Gas-dominated giant planets that accrete most of their gas from this region will therefore be characterized by lower C/O* ratios than in the inheritance low scenario.

The information provided by the C/O* ratio on the planet accretion history can find independent confirmation in the total metallicity of the planet envelope. Gas-dominated giant planets are characterized by substellar envelope metallicities that decrease with migration, as highlighted also in Paper I. Solid-enriched giant planets are instead characterized by the opposite trend, i.e., superstellar envelope metallicities that increase with migration. Note that here the relevant quantity is the planetary metallicity normalized to that of the host star. The choice of the right stellar reference value is therefore key to the correct interpretation of observational data. The C/O ratio and metallicity of the Sun are widely used, either explicitly or implicitly, as referenced in interpreting the composition of exoplanetary atmospheres (see, e.g., Madhusudhan 2019 for a recent review). However, both the C/O ratio and the metallicity of stars in planetary systems can significantly deviate from their respective values in the Sun (e.g., Delgado Mena et al. 2010; Mulders 2018 and references therein; Magrini et al. 2022). Therefore, to avoid introducing biases, one has to look at the planetary elemental ratios and metallicity in the reference frame of the host star. For instance, a planetary C/O ratio estimated to be 0.55 (equal to the solar C/O ratio) would be stellar only if the planet formed around a solar-type star. If, instead, the planet formed around a star characterized by a subsolar C/O ratio (e.g., 0.45), the planetary C/O ratio should be correctly interpreted as superstellar. Our methodology has been recently successfully applied by Carleo et al. (2022) and Guilluy et al. (2022) to the interpretation of observational data from the GAPS 2 program.

It is important to notice that when the C/O* ratio can only be estimated as being either sub- or superstellar, it is not possible to constrain the planet's accretion history, unless independent measurements of the envelope metallicity are available.

Our findings reveal that using the C/O* ratio alone to probe the formation pathways of giant planets could lead to wrong conclusions. The limited diagnostic power of the C/O* ratio is essentially due to the high volatility of both C and O, which makes the ratio unable to unequivocally trace the accretion of gas and solids. The inclusion of N, one of the most volatile elements, in the set of elemental tracers allows the computing of the C/N* and N/O* ratios and the breaking of the degeneracy on the accretion history. Specifically, Figure 8 shows that for both the C/N* and N/O* ratios, the trends for gas-dominated and solid-enriched giant planets are well separated and always characterized by opposite behaviors. Such behaviors result from the lower volatility of C and O with respect to N. Specifically, while the bulk of N in disks remains in the gas phase as N₂, the abundances of gaseous C and O decrease toward the outer disk (see Figure 3). Variations of the C/N* and N/O* ratios in disks are therefore driven by variations of C* and O*. As a consequence, the C/N* (N/O*) ratio of gas-dominated giant planets is globally substellar (superstellar) and decreases (increases) with the radial migration. The behavior of solid-enriched giant planets is the opposite, with C/N* (N/O*) globally superstellar (substellar) and increasing (decreasing) with migration.

It is worth noticing that the N/O* ratio of gas-dominated giant planets is always markedly superstellar and never drops below 2. Such a trend is due to the fact that in our model about 50% of O is early sequestered into rocks, hence subtracted from the gas phase. As a consequence, in short-scale migration scenarios, even if the other elemental ratios approach the stellar value, the N/O* ratio does not.

What emerges from our study is that in order to get a complete picture of planet formation, one has to look at multiple elemental ratios. One of the advantages of using normalized ratios is that their trends can be more readily and intuitively compared. The comparison of the top two panels in Figures 4–7 reveals the existence of fixed relations between the ratios of gas-dominated and solid-enriched giant planets. In particular, regardless of the chemical scenario, the N/O* ratio of gas-dominated giant planets is always larger than the C/O* ratio, which is in turn larger than the C/N* ratio. The opposite relation holds for solid-enriched giant planets for moderate and large-scale migrations. Such relations result from the differences in relative volatility between C, O, and N. More volatile elements are characterized by lower condensation temperatures. As such, they condense in solid form farther away from the star and remain in gas form over wider disk regions. As discussed in Section 2.4, in our compositional model, N is almost entirely present in the gas phase as N₂ in all four scenarios. Therefore, as a result of its high volatility, N will be proportionally more abundant than C and O in giant planets that accreted only gas, hence leading to N/O* > C/O* > C/N*. For solid-enriched planets, which accreted from the solid phase that is progressively depleted in highly volatile elements, the opposite is true, and C/N* > C/O* > N/O*. The emerging picture is that a joint evaluation of normalized elemental ratios allows us to immediately and unequivocally discriminate between planets that derive their metallicity from the accretion of gas and planets that experienced planetesimal enrichment. These results expand those of Paper I, showing how the relative values of the normalized ratios follow the same patterns independently on the underlying chemical scenario. These behaviors have been recently confirmed in the studies of Biazzo et al. (2022), Carleo et al. (2022), and Guilluy et al. (2022), based on Paper I .

Among all the ratios, C/O* appears to be the least robust as a diagnostic tool for the migration history of the planet. In particular, Figure 8 shows that the C/O* ratio exhibits limited changes among the migration pathways in almost all chemical scenarios. Its overall flatness is essentially due to the limited range of variation of the C/O* ratio in the disk with respect to the other ratios, as it is shown in the bottom four panels of Figures 4–7 (note the different scales on the y-axes). The strongest effects on the planetary C/O* ratio are observed in the two reset scenarios for gas-dominated giant planets, where the values decrease by a factor of 3 (i.e., from 1.6 to 0.5), for formation regions increasingly farther away from the star. In particular, the trend decreases steadily up to 50 au and progressively flattens further out. We discussed the origin of such a trend in Section 3.3. On the contrary, for both gas-dominated and solid-enriched giant planets, the C/N* and N/O* ratios are significantly more sensitive to the radial migration than the C/O* ratio. Specifically, they show monotonic trends with deviations from the stellar value that increase with direct dependency on the migration, leading to separations up to one order of magnitude between short- and large-scale migrations.

The S/N* ratio provides additional information for solid-enriched giant planets. As it is shown in the top-right panel of Figures 4–7, and highlighted in Figure 9, large-scale migrations of solid-enriched giant planets result in S/N* ratios systematically above C/N* ratios. Although the difference between the two may appear marginal, the trend results from the lower volatility of S with respect to C, which causes a larger fraction of S to be locked in the solid phase than C. Given that the total budget of S and C is dominated by the accretion of solids for large-scale migrations, the abundance of S grows faster than that of C, resulting in S/N* > C/N*. On the contrary, short-scale migrations of solid-enriched giant planets result in S/N* ratios markedly lower than C/N* ratios. Such a trend is a consequence of the increasingly dominant contribution of the gas accretion (see the discussion in Section 3.1) combined with the fact that the gas is the major carrier of C within the snowline of refractory organic carbon. Our results confirm that the relation between the S/N* and C/N* ratios holds for all the disk chemical scenarios in which the bulk of N remains in gaseous form while the bulk of S is early trapped in refractory material, so that the planetary S/N* ratio becomes a direct tracer of the planetesimal accretion.

Zoom In Zoom Out Reset image size

Given its intrinsic properties, the S/N* ratio is an extremely versatile diagnostic tool. In particular, the joint evaluation of the S/N* ratio and the envelope metallicity of the planet allows discrimination between two scenarios that are not directly modeled by our simulations. The first one describes a giant planet that formed in situ and close to the star, in a hot environment with T > 700 K. At such a high temperature, all the volatile elements, including S, are in the gas phase. In cosmochemistry, S is indeed the most abundant among the moderately volatile elements (Palme et al. 2014), with condensation temperature T_cond = 672 K (Wood et al. 2019). Consequently, the planet would be characterized by stellar values of both the envelope metallicity and the S/N* ratio. The second scenario is that of a solid-enriched hot Jupiter that accreted most of its heavy-element budget beyond the N₂ snowline, where even the highly volatile molecules (CO and N₂) are condensed in the solid phase (Öberg & Wordsworth 2019). At the end of its migration, when the planet reaches its close-in orbit, its envelope would be characterized by a markedly superstellar metallicity and a stellar S/N* ratio. Therefore, the joint evaluation of the metallicity and the S/N* ratio allows distinguishing hot Jupiters that formed in situ and close to the star from those that underwent extreme migrations.

The combined use of metallicity, C/N*, and S/N* ratios allows the investigation of another alternative scenario that was not directly modeled by our simulations. This is the case of a giant planet that formed in a disk in which the inward drift of dust and pebbles causes the ices of C, O, and N to sublimate as they cross their respective snowlines (Booth & Ilee 2019; Schneider & Bitsch 2021), thus directly enriching the disk gas phase in heavy elements. In such scenario, gas accretion may then contribute more than solid accretion to the total metallicity of the planet, hence mimicking the effect of short-scale migration in our simulations. Due to the higher volatility of C with respect to S, sublimation would be more efficient in enriching the gas in C rather than in S. Therefore, hot Jupiters that form in such scenario, will be characterized by superstellar metallicity and C/N* > S/N*.

As discussed in Sections 3.3 and 3.4, our new analysis confirms and expands the general conclusion of Paper I about the limited predictive power of the C/O* ratio. Nevertheless, the C/O* ratio does retain some diagnostic power in our model. Specifically, the C/O* ratio of gas-dominated giant planets that underwent moderate and large-scale migrations (beyond 10 au in our model) is always superstellar in the inheritance scenarios and substellar in the reset ones (see the top-left panel in Figure 8). Therefore, provided that it is possible to independently unveil the origin of the accreted material, e.g., by looking at the other elemental ratios or through density and metallicity estimates, the C/O* ratio of gas-dominated giant planets can be used to constrain the chemical structure of the birth environment. On the other hand, little can be said about solid-enriched giant planets, whose C/O* ratio does not deviate appreciably from the stellar value in both scenarios.

For gas-dominated giant planets, the N/O* ratio in the atmosphere provides additional insights into the chemical characterization of the disk. As it is shown in the bottom-left panel of Figure 8, the results for the inheritance and the reset scenarios follow two distinct trends that diverge with respect to each other. The separation between the two increases with the radial migration, resulting in a N/O* ratio systematically higher in the inheritance scenarios than in the reset ones. Specifically, for large-scale migrations, the N/O* ratio in the inheritance scenarios can be more than a factor of three larger than in the reset ones. Such a trend is a direct consequence of the higher abundance of gaseous O in the outer disk in the reset scenarios with respect to the inheritance ones (see Figure 3). We discussed this feature in Section 3.3 to explain why the C/O* ratio of gas-dominated giant planets is substellar in the reset scenarios.

To provide a more immediate visualization of the characteristics of planetary atmospheres and to facilitate the interpretation of observational data, we compared pairs of normalized elemental ratios through binary plots. Figure 10 illustrates this comparison for three pairs of ratios, C/N* versus C/O*, N/O* versus C/O*, and N/O* versus C/N* for both gas-dominated and solid-enriched giant planets in the four chemical scenarios. The gray dashed lines, indicating stellar values, delimit four quadrants in which the ratios can assume sub- or superstellar values. One of the advantages of this representation is that the results for different combinations of the parameters of the system (i.e., the extent of migration, the accretion, and chemical scenarios) occupy distinct quadrants. One can then immediately identify the conditions that favor one scenario over the others.

Zoom In Zoom Out Reset image size

For instance, a point in the third quadrant (bottom left) of the C/N* versus C/O* plot, i.e., substellar values for both the C/N* and C/O* ratios, identifies a planet that started its migration far from the star (beyond the CO₂ snowline) and accreted only gas in a reset scenario. Conversely, a superstellar C/O* ratio coupled with a C/N* ratio an order of magnitude lower than the stellar value indicates that the same formation process occurred in an inheritance scenario.

Moreover, the comparison between the N/O* and C/N* ratios allows for constraining the accretion history and the migration of giant planets. Specifically, the deviations of the ratios from their respective stellar value are directly proportional to the radial migration and grow in opposite directions for gas-dominated and solid-enriched giant planets. In the inheritance scenarios (see the two bottom-left panels), the N/O* and C/N* ratios of gas-dominated giant planets can deviate up to one order of magnitude above and below the stellar value, respectively. For the same planets, moderate deviations of the N/O* ratio with respect to the C/N* ratio provide a strong indication that the evolution occurred in a reset scenario (see the two bottom-right panels). When considering solid-enriched giant planets, both the N/O* and C/N* ratios are less sensitive to migration than in gas-dominated planets and follow the same trend independently on the chemical scenario.

Alongside favored regions of the parameter space for the different planet formation histories, our results suggest forbidden (as indicated by the yellow shaded areas in Figure 10) or marginally permitted quadrants in pairwise elemental ratios analysis. For instance, in no case among those we investigated in our modeling can planets populate the quadrant where both the N/O* and the C/N* ratios are substellar (see the bottom four panels of Figure 10). The opposite case, i.e., N/O* and C/N* ratios simultaneously superstellar, is very marginally compatible only with solid-enriched giant planets that experienced very short-scale migrations, starting within the H₂O snowline.

One of the main findings of Paper I, confirmed by this study, is the untapped potential of the S/N* ratio as a diagnostic tool for planet formation. As discussed in Section 3.4, this is essentially due to the marked volatility contrast between S and N, which makes their ratio an effective tracer of the accretion of solids. It is important to point out, however, that any element X more refractory than O can be used in place of S to perform the same analysis. From an observational perspective, this means that our set of diagnostic tracers can be extended to elements that may be more easily detectable than others in giant planet atmospheres. For instance, for increasing equilibrium temperatures of the planets, the atmospheres get enriched in elements with increasingly lower volatility. For warm and hot Jupiters, moderately volatile elements are all good alternatives to S. Such category includes elements that condense at temperatures between 1250 K and 250 K, such as Na, K, Cl, and F (Palme et al. 2014; Wood et al. 2019). For ultrahot Jupiters the list extends to refractory (e.g., Al, Ca, and Ti) and rock-forming (e.g., Mg, Si, Ni, and Fe) elements, which have condensation temperatures in the range 1850–1250 K (Palme et al. 2014; Wood et al. 2019).

In Section 3.4 we stressed that the planetary elemental ratios strongly depend on the phase of the main carrier of the two involved elements (e.g., the gas phase for N₂, which is the main carrier of N and the solid phase for the rocks that are the main reservoir of S) and on how efficiently the gas and the solids contribute to the final envelope metallicity. In particular, we examined the case of the C/N* and S/N* ratios in solid-enriched giant planets. Which of the two ratios is larger than the other depends on whether the planet's formation is dominated by the accretion of gas or solids. More generally, the comparison between two elemental ratios X/N* and Y/N*, where the elements X and Y have a marked volatility contrast, allows for discrimination between scenarios in which one of the two elements is predominantly accreted by different phases.

In Section 3.4 we also highlighted that the S/N* ratio allows to resolve scenarios of in situ formation very close to the star, in a hotter environment than the condensation temperature of S. Specifically, giant planets that formed by accreting gas from such a region will be characterized by stellar values of both metallicity and S/N* ratio. It is important to notice that the replacement of S with a more refractory element would push this limit toward higher temperatures, allowing to resolve formation scenarios even closer to the star.

In Section 3.4 we also discussed how the comparison between the metallicity and the C/N* and S/N* ratios allows the identification of giant planets that formed by accreting gas enriched in heavy elements due to the sublimation of inwardly drifting dust and pebbles, contrary to planetesimal accretion. Sublimation is more efficient in enriching giant planets in highly volatile elements rather than in refractory ones. In particular, any refractory element X would produce less enrichment in the gas than in C. Therefore, planets that formed by accreting high-metallicity gas would be characterized by superstellar metallicity and C/N* > X/N*.

In the previous sections we discussed how the atmospheric elemental ratios can be used to trace the formation history of giant planets. From an observational perspective, such elemental ratios need to be retrieved from the abundances of the molecular carriers of C, O, N, and S in exoplanetary atmospheres.

The observable chemistry in giant planet atmospheres is shaped by multiple and diverse processes, such as thermochemical reactions at equilibrium, mixing processes, photochemistry, and chemical diffusion. Each of them dominates in different atmospheric layers (e.g., Madhusudhan et al. 2016 and references therein).

Deep in the atmosphere, the chemistry is at equilibrium and the chemical composition is a function only of the elemental abundances, the temperature, and the pressure. For a solar composition, the main molecular carriers of C and N are CO, CH₄, N₂, and NH₃. Two regimes in the pressure–temperature space can be identified. At high temperatures and low pressures, C and N are essentially locked in CO and N₂, while at low temperatures and high pressures, the main reservoirs of C and N are CH₄ and NH₃. At the pressure of 1 bar, the transition between the two regimes occurs at temperatures of 1200 K and 700 K for the C- and N-bearing molecules, respectively. O is largely present as water vapor, while S is mainly locked into H₂S (e.g., Madhusudhan et al. 2016; Fortney et al. 2021 and references therein). The chemical equilibrium in atmospheres with supersolar metallicity generally favors the production of molecules that contain two or more heavy elements. In such atmospheres, CO₂ is therefore favored over CO and N₂, which in turn are favored over CH₄ and NH₃ (e.g., Madhusudhan et al. 2016 and references therein).

In the upper layers of the atmosphere, photochemical reactions induced by stellar UV radiation are proved to alter the chemical composition. Among the direct and indirect products of photochemistry there are H, HCN, C₂H₂, OH, O₂, and NO. Regarding S, photochemistry converts H₂S in S, HS, S₂, SO, and SO₂ (e.g., Madhusudhan et al. 2016; Fortney et al. 2021).

Chemistry in observed exoplanetary atmospheres may be more complex than discussed above, being also influenced by the cooling history of the planet (Fortney et al. 2020). At intermediate layers, volatile species may be sequestered into clouds, although these are not expected to form at the high atmospheric temperatures of hot Jupiters (e.g., Madhusudhan et al. 2016; Madhusudhan 2019 and references therein).

Ground-based and high-resolution spectral observations have recently succeeded in providing the first characterizations of exoplanetary atmospheres. For instance, observations by Giacobbe et al. (2021) of the hot Jupiter HD 209458b, at the equilibrium temperature of 1500 K, revealed the presence of H₂O, CO, HCN, CH₄, NH₃, and C₂H₂ in the atmosphere. Carleo et al. (2022) and Guilluy et al. (2022) have recently extended this result to giant planets with equilibrium temperature T ≤ 1000 K. Future observations with JWST (Greene et al. 2016), Twinkle (Edwards et al. 2019b), and Ariel (Tinetti et al. 2018), as well as advances in atmospheric retrieval techniques, will soon revolutionize the field. The broad spectral range and the high spectral resolution of JWST will allow for precise determination of the abundances of a number of species of interest besides H₂O, like CO, CO₂, CH₄, and NH₃ (Greene et al. 2016). An even greater number of molecules will be accessible with Ariel, by the time the telescope starts its operations in 2029 (Tinetti et al. 2018).

Recent studies have validated our methodology on real data in preparation of its use for interpreting the observations of future ground-based and space-based facilities. Alongside Carleo et al. (2022) and Guilluy et al. (2022), Biazzo et al. (2022) put constraints on possible formation scenarios for the giant planets in four planetary systems, around the stars HAT-P-26, WASP-39, HAT-P-12, and WASP-10. Specifically, they analyzed the trends of the normalized planetary metallicity and elemental ratios based on the atmospheric composition obtained by Kawashima & Min (2021) through spectral disequilibrium retrieval models. A similar study was performed independently by Kolecki & Wang (2022) on an additional set of four giant planets. As discussed before, Carleo et al. (2022) and Guilluy et al. (2022) jointly used the planetary metallicity and the C/O* ratio to constrain the formation histories of WASP-69b and WASP-80b, as well as the characteristics of their native circumstellar disks (see also Section 3.3). More systematic studies of the observational implications of our results will be the focus of dedicated future works.