Composition of Early Planetary Atmospheres II: Coupled Dust and Chemical Evolution in Protoplanetary Disks
A.J. Cridland, Ralph E. Pudritz, Tilman Birnstiel, L. Ilsedore Cleeves, Edwin A. Bergin
aa r X i v : . [ a s t r o - ph . E P ] M a y Mon. Not. R. Astron. Soc. , 1–19 (2015) Printed 30 April 2019 (MN L A TEX style file v2.2)
Composition of Early Planetary Atmospheres II: CoupledDust and Chemical Evolution in Protoplanetary Disks
A.J. Cridland ⋆ , Ralph E. Pudritz , , , † , Tilman Birnstiel ,L. Ilsedore Cleeves & Edwin A. Bergin Department of Physics and Astronomy, McMaster University, Hamilton, Ontario, Canada, L8S 4E8 Origins Institute, McMaster University, Hamilton, Ontario, Canada, L8S 438 Zentrum f¨ur Astronomie der Universit¨at Heidelberg, Institut f¨ur Theoretische Astrophysik, Albert-Ueberle-Str. 2, 69120 Heidelberg, Germany Max Planck Institute for Astronomy K¨onigstuhl 17, D-69117 Heidelberg, Germany University Observatory, Faculty of Physics, Ludwig-Maximilians-Universit¨at M¨unchen, Scheinerstr. 1, 81679 Munich, Germany Harvard-Smithsonian Cneter for Astrophysics, 60 Garden Street, Cambridge, MA 02138, USA Department of Astronomy, University of Michigan, 1085 S. University Avenue, Ann Arbor, MI 48109, USA
30 April 2019
ABSTRACT
We present the next step in a series of papers devoted to connecting the composi-tion of the atmospheres of forming planets with the chemistry of their natal evolvingprotoplanetary disks. The model presented here computes the coupled chemical anddust evolution of the disk and the formation of three planets per disk model. Our threecanonical planet traps produce a Jupiter near 1 AU, a Hot Jupiter and a Super-Earth.We study the dependency of the final orbital radius, mass, and atmospheric chemistryof planets forming in disk models with initial disk masses that vary by 0.02 M ⊙ aboveand below our fiducial model ( M disk, = 0 . M ⊙ ). We compute C/O and C/N forthe atmospheres formed in our 3 models and find that C/O planet ∼ C/O disk , whichdoes not vary strongly between different planets formed in our model. The nitrogencontent of atmospheres can vary in planets that grow in different disk models. Thesedifferences are related to the formation history of the planet, the time and locationthat the planet accretes its atmosphere, and are encoded in the bulk abundance ofNH . These results suggest that future observations of atmospheric NH and an esti-mation of the planetary C/O and C/N can inform the formation history of particularplanetary systems. Key words: protoplanetary discs, interplanetary medium
What physical processes set the observable chemical abun-dances in exoplanetary atmospheres? The bulk abun-dance of the atmosphere will ultimately depend on thechemical structure of the gas that is present in theplanet’s natal protoplanetary disk ( ¨Oberg & Bergin 2016;Madhusudhan et al. 2014a, 2016; Mordasini et al. 2016).The disk astrochemistry depends on the physical state ofthe disk, which includes the temperature, density, and ion-ization of the gas, and radiative flux of high energy pho-tons (Fogel et al. 2011; Helling et al. 2014). In this and pastworks (Alessi et al. 2016; Cridland et al. 2016, 2017) wepresent a model of an evolving astrochemical disk which ⋆ E-mail: [email protected] † E-mail: [email protected] when combined with the core accretion model for planetaryformation (Ida & Lin 2004) can predict the bulk chemicalabundance of planetary atmospheres.An ‘early’ planetary atmosphere is a young ( t age . O, c (cid:13) A.J. Cridland et al. CO , CO, CH , K, and Na for about a dozen planets. Earlyconnections between theoretical work and observations iden-tified the carbon-to-oxygen elemental ratio (C/O) as an ex-cellent candidate for linking the planet formation processand the resulting atmospheric chemistry ( ¨Oberg et al. 2011;Madhusudhan et al. 2014b; Thiabaud et al. 2015). This ra-tio varies with disk radius as volatiles are either frozen outor sublimated into the gas phase. Hence C/O of a planet’satmosphere will depend on where that planet accreted itsgas, relative to the ice line of volatiles like H O, CO , andCO ( ¨Oberg et al. 2011).A recent synthetic model that calculated the links be-tween the formation and migration of hot Jupters in evolvingdisks, to the construction and composition of their atmo-spheres used simplified models of the composition of plan-etesimals (Mordasini et al. 2016). The resulting C/O ratiosin this analysis were subsolar because of the O rich natureof the accreting solid materials.In this paper we introduce a new method of computingthe coupled chemical and dust evolution of the disk, in or-der to give a first combined treatment of these two criticalaspects of disk evolution and early planetary formation. Weuse this approach to then address key problems in linkingplanet formation with the resulting composition of planetaryatmospheres. One application of this approach is to demon-strate that the famous C/O and the carbon-to-nitrogen ratio(C/N) may only weakly dependent on the formation his-tory of planets. And that the particular abundance of ni-trogen carrying molecules trace details in the planet’s for-mation history like the timing of the planet’s gas accretion(Cridland et al. 2016).Of particular interest to our work is whether thebulk chemical compositions of these ‘early’ planetary atmo-spheres can model the bulk compositions that are inferredby emission and transmission spectra of observed planets.In doing so we might learn from where the diversity of at-mospheric chemistry is inherited.The paper is structured as follows: in Section § § § § § § § Our model involves the combination of accretion diskphysics, dust physics (coagulation, fragmentation, radialdrift and settling), radiative transfer, astrochemsitry (pho-tochemistry, gas phase, and grain surface) and planet for-mation theory (core accretion and migration). A particu-larly important aspect of the planet formation model is the use of ‘planet traps’ to set the location of the protoplanetand limit its rate of migration. Planet traps are inhomo-geneities in disk properties where the net torque vanishes,so that migration occurs at the much slower rate of diskevolution. Introduced by Masset et al. (2006) to investigateplanetary dynamics at the inner edge of a disk dead zone,planet traps have been generalized and used in populationsynthesis models to replicate the number of planets observedin different regions of the mass-semi-major axis diagram (forex. Hasegawa & Pudritz (2011, 2012, 2013)).A planet trap works because of the relative strength ofthe Lindblad and co-rotation torques (Masset et al. 2006).The Lindblad torque is caused by the gravitational pull ofmaterial that is perturbed from Lindblad resonances nearthe planet back upon the planet (Ward 1997). The two near-est resonances are at radial positions with orbital frequen-cies that are twice the planet’s (outer resonance) and halfthe planet’s (inner resonance) orbital frequency. Generallythe outer resonance is stronger than the inner resonance andthe planet migrates to smaller orbital radii.The co-rotation torque is caused by gas with orbitalfrequencies that are close the planet’s. The gas enters intohorseshoe orbits where they oscillate between orbital radiithat are smaller than the planet’s to orbital radii that arelarger. At the turn of each horseshoe orbit the gas exchangesangular momentum with the planet (Paardekooper et al.2010). If each turn is symmetrical the co-rotation torquedoes not exert a net torque on the planet - this is calledsaturation. With enough turbulence the gas mixes with gason the opposite side of the horseshoe orbit which causes anasymmetry in the two turns and a net torque which generallycauses outward migration (Seager 2010; Paardekooper et al.2010). Outside of planet traps, the co-rotation torque is gen-erally weaker than the Lindblad torque.In our work we are interested in three planet traps: thewater ‘ice line’, the heat transition and the dead zone. Thewater ice line depends on the chemical state of the gas andis located where the water content of the disk transitionsfrom being primarily in the vapour phase to being in theice phase. The heat transition traces the location where thedisk’s heating mechanism transitions from being primarilydue to viscous stresses to direct irradiation from the hoststar. Finally the dead zone is located at the transition pointbetween a low ionization state where the disk turbulence isweak (‘dead’), to a high ionization state where the turbu-lence is active (Matsumura et al. 2009).Our fiducial model (CPA16) includes all of the aboveeffects, but assumed a constant dust to gas ratio at all timesthroughout the disk. In that work we show that the chem-ical abundance of planetary atmospheres depends on both where and when the planet accretes its gas. We observe thatC/O and C/N ratio encode information regarding the for-mation history of the planet. Atmosphere C/O ratios thatwere nearly the same as the disk’s initial C/O ratio werefound for planets that formed within the water condensa-tion front (ice line) of the disk. We expect that higher C/Oratios (C/O ∼
1, compared to solar: 0.54 and our initial disk:0.288) will be found in planets that accrete their atmosphereoutside the ice line, however this was not demonstrated inthat previous paper. The sub-solar C/O of our initial diskwas due to our choice of the initial molecular species that areinherited by the disk. Our initial abundances are based on c (cid:13) , 1–19 omposition of Early Planetary Atmospheres II the work of Fogel et al. (2011) and Aikawa & Herbst (1999)and represents the chemical state of a dense molecular cloud.The C/N ratio is encoded when the planet accretes itsatmosphere. Atmospheres that were accreted early ( t . t > r . and HCN which leads to higher C/N than later whenthe disk is cooler and nitrogen is found primarily as N .We incorporated the effects of dust coagulation, frag-mentation, and radial drift to our model in CPB17 anddemonstrated that the gas disk became more ionized whencompared to the simpler dust model of CPA16 because theradial drift and growth of the grains reduced the disk’s opac-ity to high energy radiation. The ionization structure of thedisk is important to our work because it sets the locationof one of the planet traps: the dead zone. We assume thatthe source of turbulence in the disk is driven by the magne-torotational instability (MRI) which is caused by a couplingof the disk’s magnetic field to free electrons in the gas (forex. Balbus & Hawley (1991)). The growth of turbulence dueto the MRI is counteracted by the diffusion of the magneticfield. We assume that the primary source of diffusion onthe midplane of the disk is Ohmic resistivity (Bai & Stone2013; Gressel et al. 2015). This assumption ignores othernon-ideal effects like ambipolar diffusion (Bai & Stone 2013;Gressel et al. 2015) and the Hall effect (Bai & Stone 2016).The impact of these other sources of diffusion are more com-plicated to estimate and are part of an ongoing discussionregarding the turbulent structure of protoplanetary disks.In this paper, we compute the formation of the threeplanets and the chemical structure of their atmospheres inthe resulting disk model from CPB17. We will compare theresults of our work to the results of CPA16 where we usethe same disk parameters, but change the way that the dustphysics is handled. Importantly, the more complicated dustphysics shows that the disk’s dead zone evolved radiallymuch faster than in the disk model with a constant dust togas ratio without grain coagulation or fragmentation. Thefaster evolution causes the planet forming in the dead zonetrap to have a different formation history, and potentially adifferent atmospheric composition. In CPA16 the dead zoneplanet failed to accrete an atmosphere in the 4.1 Myr disklifetime.With the need for turbulence to maintain an unsatu-rated co-rotation torque, the relative location of the deadzone edge and other planet traps dictate whether the planetremains trapped over its formation lifetime. Within the deadzone, we assume that the turbulent α drops by two ordersof magnitude when compare to the active region. This re-duction leads to a lengthened eddy turnover time, and asaturated corotation torque if the eddy turnover time ex-ceeds the libration time of the horseshoe orbit (see for ex-ample CPA16). In CPA16 we found that the planet formingin the heat transition trap saturated 1.1 Myr into forma-tion, causing the planet to quickly migrate and become aHot Jupiter. In CPB17 we showed that the dead zone waslocated at smaller radii than the heat transition trap by 1.3Myr. If the dead zone crosses the heat transition trap’s loca- tion the co-rotation will remain unsaturated and the planetwill remain trapped for its entire formation history. This willproduce a different planet than the one reported in CPA16.Finally it was pointed out in CPA16 that to fully under-stand the connection between the atmospheric abundancepredictions and our model, we would require a statisticallysignificant sampling of possible disk parameters like initialmass and lifetime. In Alessi et al. (2016) it was shown thatboth of these disk parameters impact the formation his-tory of planets in our formation model. In particular theinitial mass of the disk changes the location and evolu-tion of the planet traps, hence the location of the form-ing planet. Ideally, a full population synthesis model (ie.Ida & Lin (2004); Hasegawa & Pudritz (2013); Alessi & Pu-dritz (in prep.)), with values drawn from distributions of afew parameters, would be applied in our evolving astrochem-ical disk. In preparation for a future complete populationsynthesis study, we report results from two additional diskmodels where we have varied the initial disk mass by 0.02 M ⊙ above and below our fiducial disk of 0.1 M ⊙ . Our evolving astrochemical model as well as our method ofidentifying the planet trap locations is outlined in detail inCridland, Pudritz & Alessi (2016). Our dust model is out-lined in Cridland, Pudritz, Birnstiel (2017) and was basedon the model of Birnstiel et al. (2012). Here we outline theimportant aspect of each of these models.Generally we compute the temperature and surface den-sity radial distribution of the gas with an analytic solutionof the gas diffusion equation (see Chambers (2009), CPA16).This self-similar model for disk structure and evolution as-sumes that the heating of the gas is dominated by eitherviscous evolution or direct irradiation of the host star. Fur-thermore, as a computation constraint we assume that thedisk’s opacity and α -parameter (ie. in Shakura & Sunyaev(1973), due to the presence of disk winds and/or turbulence)are constant in space and time, and that the mass accre-tion rate is constant in space but evolves with time (seeAlessi et al. (2016) and CPA16 for details).This gas model is combined with a model for the evolu-tion of the dust surface density (Birnstiel et al. 2012). Theparticle size distribution is then reconstructed with the semianalytical treatment of Birnstiel et al. (2015). The dust andgas distributions are computed over many ‘snapshots’ overa disk lifetime of 4.1 Myr.In each snapshot we compute the flux of UV and X-rayphotons with RADMC3D, a Monte Carlo radiative transferscheme (Dullemond 2012). The wavelength dependent diskopacity is computed by the on-the-fly method in RADMC3Dand depends on the dust surface density and size distribu-tion using optical constants from Draine (2003). We com-pute the flux of 4 UV and 3 X-ray test wavelengths, thenextrapolate their results for the wavelengths: 930 - 2000angstroms for UV and 1 - 20 keV for X-ray using sampleT Tauri UV (Fogel et al. 2011; Bethell & Bergin 2011) andX-ray (Kastner et al. 1999, 2002; Cleeves et al. 2015) spec-tra as guides.We note that relying solely on dust as our source ofdisk opacity does ignore the X-ray cross section of the gas. c (cid:13) , 1–19 A.J. Cridland et al.
Over our range of wavelengths dust does dominate the X-ray cross section (Bethell & Bergin 2011), however we willlikely underestimate the disk’s opacity late in it’s evolutionas the dust surface density drops. At these late times, thegas will begin to dominate the cross section and set a floorto the disk’s opacity.Next we compute the astrochemistry using a chemicalkinetic code (ie. Fogel et al. (2011); Cleeves et al. (2014),CPA16) in each snapshot using the gas, dust, and radiationas inputs. This process produces an evolving astrochemicaldisk in which we computed the formation of a planet usingthe core accretion model and the migration of the formingplanet with a planet trap model. Some specific features ofthis theoretical framework are discussed below.
An important feature in our astrochemical model is the factthat the temperature and density profiles evolve with timeas the disk ages. The chemistry is handled with a non-equilibrium code from Fogel et al. (2011) and Cleeves et al.(2014). This astrochemical code includes, among others, ion-driven gas phase reactions, the freeze out of volatiles, photo-desorption, and the production of water through grain sur-face reactions. In the chemical model we compute an averagegrain size based on the results of our dust model (see below),and allow the average grain size to vary at different disk radiiand evolve in time.We combined this chemical model with an analytic diskmodel from Chambers (2009) which models a disk thatevolves due to mass accretion caused by viscous stresses.This analytic model produces temperature and surface den-sity profiles that scale as: T ∝ R β ; Σ ∝ R s , (1)where the exponents β and s are determined by the sourceof heating. Viscous heating and direct irradiation are theprimary heating sources. Viscous heating is caused by grav-itational energy being released as material accretes throughthe disk. This heating scales with the mass accretion rate ˙ M which is reduced as time passes. Because of this dependence,the region of the disk that is viscously heated will cool asthe disk ages. At larger radii the densities are lower and theheating is dominated by the direct irradiation of the diskfrom the host star. We assume that the stellar properties ofthe host star do not change, so the temperature profile ofthe disk is also static. This assumption ignores the effect ofexcess high energy photons due to the accretion shock onthe heating of the disk. In doing so we have likely underes-timated the temperature of the disk in the irradiated regionat the earliest times. However, as the disk ages and the ac-cretion rate drops the heating from the accretion shock isdrastically reduced. In a future work we will explore the im-pact of both the temperature structure and chemistry of anevolving stellar emission spectrum.Viscous heating dominates at smaller radii where thedensities are the highest, but as the disk accretes its densityand mass accretion rate both drop. As the mass accretionrate drops, the region of the disk where viscous heating dom-inates shrinks and the point where the temperature from vis-cous heating is the same as the temperature from irradiation(the heat transition trap) moves to smaller radii. The evolving viscously heated region of the disk also hasimportant implications to the chemistry in the disk. Featureslike the water ice line, which closely depends on the temper-ature of the gas, will evolve to lower radii as the disk cools.Additionally the main molecular carrier of some elements,like nitrogen, can change depending on the temperature ofthe gas. In hotter disks the nitrogen will generally be foundin NH and HCN, while in cooler disks it is mostly foundin N (CPA16). This change in elemental carrier will beimprinted in the bulk chemical abundance of the planetaryatmosphere based on where and when the planet accretes itsgas. Just as the temperature profile evolves, the surface den-sity profile of both the gas and dust also changes with time.The gas density drops as mass is accreted onto the hoststar due to turbulent viscous stresses and magnetocentrifu-gal winds. Likewise the dust also accretes onto the host starthrough viscous stresses, however the it is also susceptible toa faster source of evolution: radial drift. Radial drift involvesthe sub-Keplarian orbit of dust grains caused by a gas dragwhich scales with the surface area of the grain. As a resultthe largest grains are most affected by radial drift, whichwhen combined with coagulation causes a faster reductionin surface density than viscous evolution alone. The depen-dency of the radial drift rate with grain surface area extendsup to a point where the grain’s inertia becomes too great forgas drag to have a large impact on the grain’s orbit. Our sizedistribution extends up to a few tens of cm, and so we do notsimulate the dust evolution up to this cut off. The growthand drift of dust grains lead to an overall reduction in thedisk’s opacity to high energy photons, and an increase in theionization of the gas (CPB17). While we do evolve the ra-dial drift of dust grains (see below), we ignore the transportof frozen volatiles on the grains as they migrate inwards. Inprinciple, grains will bring volatiles from the outer disk tosmaller radii as they drift which could enhance the abun-dance of volatiles in the planet-forming regions of the disk(r .
10 AU). We leave an investigation into this effect tofuture work.Ions drive chemical reactions because of their lower ac-tivation barrier when compared to neutral-neutral gas re-actions. As mentioned above they also cause turbulencethrough the MRI when coupled to the disk’s magnetic field.Because the chemical and physical state of the gas is sensi-tive to the flux of ionizing photons, an accurate descriptionof disk opacity and evolution of the dust is an importantfeature of our model.We use the same initial chemical abundancesas in CPA16, derived from Fogel et al. (2011) and(Aikawa & Herbst 1999), see Table 1.
Our evolving dust model was presented in CPB17 andis based on the Two-population model presented inBirnstiel et al. (2012). This model involves estimating thesize and density distribution of the dust grains by comput-ing the evolution of two sample dust populations, one repre-senting the smallest grains (monomers) and one represent-ing the large population. The size distribution is assumed tobe a power law with a monomer size of 0.1 microns and a c (cid:13) , 1–19 omposition of Early Planetary Atmospheres II Table 1.
Initial abundances relative to the number of H atoms.Included is the initial ratio of carbon atoms to oxygen atoms(C/O) and the initial ratio of carbon atmos to nitrogen (C/N).Species Abundance Species AbundanceH . O 2 . × − He 0 .
14 N 2 . × − CN 6 . × − H +3 . × − CS 4 . × − SO 5 . × − Si + . × − S + . × − Mg + . × − Fe + . × − C + . × − GRAIN 6 . × − CO 1 . × − N . × − C 7 . × − NH . × − HCN 2 . × − HCO + . × − C H 8 . × − C/O 0 . . maximum size set by either fragmentation or radial drift de-pending on which physical process has the lower timescale.This simplified dust model has been tested againstfull coagulation simulations, reproducing the general trendsfrom those more complicated models (Birnstiel et al. 2012).A benefit of the Two-population model is how quickly theevolution can be computed. Depending on the treatmentof the fragmentation, standard dust evolution codes scale O ( N ) where N is the number of grain sizes used in thesimulation. So the Two-population model, which computesthe coagulation, fragmentation and radial drift of only twopopulations can be run much faster than a standard coagu-lation simulations that is resolved with at least 100 differentsizes.Our dust model includes the effect of a changing frag-mentation threshold speed across the water ice line, asfirst discussed in Birnstiel et al. (2010). The fragmentationthreshold speed is the minimum relative collision speed thatresults in the fragmentation of the grains. This thresholdspeed depends on the chemical nature of the dust grain, itsporosity, and the amount of ice coverage. For simplicity weassume a threshold speed of 10 m/s outside the water iceline because the grain is fully covered, and hence strength-ened by a layer of ice (Wada et al. 2009; Gundlach & Blum2015). Within the ice line the layer of ice is gone and thethreshold speed is reduced to 1 m/s. We used the chemicalresults of our fiducial model (CPA16) to fit the location andevolution of the water ice line as a function of time. Usingthis fit we allowed the location of the water ice line to evolvethroughout the evolution of the dust.The impact of this changing fragmentation thresholdspeed leads to a reduction in the size of the largest grainby about two orders of magnitude within the water ice line(CPB17, or see Figure 7 below). These smaller grains radi-ally drift slower than larger ones, with a drift time scale ofthe form (Birnstiel et al. 2010): τ drift = rV k Stc s γ − (2)at the radial position r . The Kepler speed ( V k ) and gas soundspeed ( c s ) are evaluated at the position r . The Stokes num-ber scales linearly with the size of the grain and encodesthe coupling between the grain and the gas. Longer drifttime scales within the ice line lead to an enhancement of the dust surface density of an order of magnitude within thewater ice line. This enhancement has two important affectson planet formation: first it shields the inner disk from highenergy photons, reducing the level of ionization early on inthe disk lifetime. Secondly, the enhancement will increasethe growth rate of the forming planet because there is morematerial available during oligarchic growth. Our planet formation model relies on the Core Accretionmodel as seen in Ida & Lin (2004) to compute the rate ofmass accretion onto the planet. This core accretion modelis used in conjunction with the planet trap model that wasoutlined above. The combination of these two models hasbeen used in population synthesis models to reproduce thepopulation statistics that have been observed in the exoplan-etary data (Hasegawa & Pudritz (2013); Alessi & Pudritz(in prep.)).
The Core Accretion model contains three important phasesof accretion. The first is known as oligarchic growth, wherea single core accretes solid material from a sea of 10-100 kmsized planetesimals. We do not model the growth of theseplanetesimals from micron sized dust grains up to their as-sumed sizes, instead we assume that all of the available solidmaterial is in planetesimals. In doing this, we neglect thereduction in the disk’s opacity caused by the rapid agglom-eration of dust particles caused by the streaming instability(see for ex. Sch¨afer et al. (2017)). The growing planet willimpact the surface density of the surrounding gas and dust,and as a result will locally change the temperature and ion-ization of the gas. These effects are beyond the scope of thiswork.The rate of accretion is fast, with timescales O (10 )yr. The solid accretion rate has the form (Ida & Lin (2004),their equations 5 and 6): dM c dt = M c τ c,acc ≃ M c . × (cid:18) Σ d − (cid:19) (cid:16) a (cid:17) − / × (cid:18) M c M ⊕ (cid:19) − / (cid:18) M s M ⊕ (cid:19) / × "(cid:18) Σ g . × gcm − (cid:19) − / (cid:16) a (cid:17) / (cid:18) m g (cid:19) / − g yr − . (3)The growing core quickly accretes all of the solid ma-terial in its immediate vicinity. At this point it has reachedits ‘isolation mass’ (Ida & Lin 2004), which represents themaximum available material to the planet if the planet doesnot migrate. Once it has reached its isolation mass the rateof planetesimal accretion drops, allowing for gas accretionto take place.We assume that the gas accretion rate is limited only c (cid:13) , 1–19 A.J. Cridland et al. by the Kelvin-Helmholtz timescale which we model with τ KH ≃ c (cid:18) M plnt M ⊕ (cid:19) − d yr , (4)so that the mass of the planet grows as dM plnt dt = M plnt τ KH , (5)for further details see Alessi et al. (2016). The form of equa-tion 4 has been demonstrated by numerical simulations, andthe constants c and d have been shown to have values be-tween 8-10 and 2-4 respectively (Ida & Lin 2004). We followIda & Lin (2004) in choosing c = 9 and d = 3. This choicewas also used in Hasegawa & Pudritz (2013) to recreate thedistribution of planets on the mass semi-major axis diagramderived from observations.While in practice the gas surface density is locally re-duced by the accreting planet (see for example Ormel et al.(2015)), we do not update the disk gas density and temper-ature while the planet draws down its atmosphere. As it forms, the protoplanet and its natal accretion diskexchange angular momentum through the Lindblad and co-rotation torques (Paardekooper 2014) which generally re-sults in a loss of angular momentum for the protoplanet.This loss causes the planet to ‘migrate’ to smaller orbitalradii in a process called Type-I migration.As was previously stated, Type-I migration is limitedby discontinuities in disk properties known as ‘planet traps’(Masset et al. 2006). Here we assume that a trapped planetwill exactly follow the radial location of the its planet trap.If the Ice line and Heat transition traps are found withinthe Dead zone, we check at every timestep whether the co-rotation torque saturates by comparing the eddy turnovertime of the gas to the gas libration time. If the libration timeis shorter than the eddy turnover time then the co-rotationtorque saturates.When the co-rotation torque saturates, the planetbreaks free of the trap, and we compute its radial evolutiondue to the Lindblad torques along (Seager 2010),Γ = C Γ Σ p Ω p a (cid:18) M p M s (cid:19) (cid:18) a | ∆ r | (cid:19) , (6)where Ω p , a , and M p are the optical frequency, radius, andmass of the planet, M s is the mass of the host star, and Σ p isthe gas surface density at the planet’s location. The constant C Γ = − (2 . − . s ) (Paardekooper et al. 2010) has beendetermined through numerical simulations. The minimumlength scale | ∆ r | below which the gas does not contributeto Lindblad torques depends on the relative impact of theplanet’s gravity to the gas pressure on the dynamics of thegas. Below the Hill radius ( R H ) gas particles become boundto the planet, while below the gas scale height ( H p ) thegravitational effects of the gas are smoothed out. Thereforewe set | ∆ r | = max ( H p , R H ).We evolve the orbital radius of the planet through ˙ a/a =Γ /J p where J p = M p a Ω p is the orbital angular momentumof the planet. As such the planet’s orbital radius evolves as,1 a dadt = C T ΣΩ p a M p (cid:18) M p M s (cid:19) (cid:18) amax ( H p , R H ) (cid:19) (7) when the co-rotation torque saturates.In what follows, when discussing the properties of ourforming planets we will refer to the planet by the trap whereit originated. The ability of a planet to accrete its atmosphere is limitedby the cooling rate of the gas as it shreds its gravitationalenergy. This cooling rate depends on the gas opacity whichmeasures the cross section of photon absorption or scatterper unit mass of gas.As the planet accretes gas it forms an envelope aroundthe planet core that then slowly contracts onto the planet.This process is modeled by two physical properties: the gasaccretion rate which depends on the planet mass and Kelvin-Helmholtz timescale, and the critical planet mass abovewhich the gas envelope can contract onto the core.In previous population synthesis models (ie.Hasegawa & Pudritz (2013)) and in our previous work(CPA16) the Kelvin-Helmholtz timescale and critical massare individually parameterized so that they have the formof Equation 4 and M crit = 10 · f crit (cid:18) ˙ M core − M ⊕ yr − (cid:19) q M ⊕ , (8)respectively. Ikoma et al. (2000) showed that the parame-ters c and f crit depend on the opacity, while the parameters d and q depend on the choice of opacity table (Ida & Lin2004).The equivalent form of Equations 4 and 8 fromIkoma et al. (2000) have the form: τ env = 3 × (cid:18) M core M ⊕ (cid:19) − /q (cid:18) κ cm g − (cid:19) s/q yr , (9)where s/q ≃ κ is the dust opacity, and M crit = 7 (cid:18) ˙ M core − M ⊕ yr − (cid:19) q (cid:18) κ cm g − (cid:19) s M ⊕ . (10)The parameters s and q are within the range of 0.2-0.3(Ikoma et al. 2000), and based on our choice of parametersfor Equation 4 imply that s = q = 0 .
25. In our normalizationof Equation 4 we assume that κ ∼ .
33 cm g − which isconsistent with the Rosseland mean opacity for a sub-mmgrain computed by our opacity table.With our assumed opacity we find that f crit ≃ . f crit = 0 . f crit is higher, gas accretion begins later in the planet’s for-mation history which changes the radial location and timingof the atmosphere formation.A smaller f crit assumes that the envelope opacity islower, which would imply some grain growth within the en-velope. Because the dust physics is not well constrained inthe envelope of accreting planets we will show the results for f crit = 0 .
2, which matches the parameters used in CPA16.In the Appendix we include the results for f crit = 1 . g − . c (cid:13) , 1–19 omposition of Early Planetary Atmospheres II The main purpose of this work is to assess the effect of dif-ferent dust models on disk astrochemistry and the resultingatmospheric abundances. To that end we first compare thechemical results of the CPA16 disk model to the disk chem-istry that is computed here, and based on the dust modelpresented in CPB17.
In Figures 1 and 2 we present the radial distribution of themidplane chemical abundances for a few chemical speciesin the CPA16 chemical model (solid line) and the chemicalresults from the disk model presented in CPB17 (dashedline). In Figure 1 we show the important chemical speciesto the observations of exoplanetary atmospheres, while inFigure 2 we show other molecules that have been observedor used as tracers in protoplanetary disks.
In the top 4 panels of Figure 1 we compare the primary car-bon and oxygen carriers in this disk midplane. An importantdifference between the chemical distributions that resultedfrom the CPA16 (solid line) and CPB17 (dashed line) dustmodels is the abundance of volatiles like H O and CO atlarger radii. For the results of this work, based on the dustmodel in CPB17, there is a higher abundance of volatilesfarther out in the disk, where one would expect the ice com-ponent of the above species would dominate. The retentionin the CPB17 dust model is connected to the fact that thelargest dust grains are more rapidly depleted than in theCPA16 dust model. As a result the amount of dust availablefor volatile freeze out is reduced in the CPB17 dust model.In addition to the reduction of available surface onto whichvolatiles can freeze, a reduction of the disk opacity to UVphotons can also lead to an enhancement of gas abundancethrough the photodesorption of frozen species.At smaller radii than the water ice line there is an en-hancement of dust surface density fueled by the inward ra-dial drift of dust grains. This enhancement shifts the loca-tion of the water ice line inward by ∼ . are also shifted from their position in the CPA16model. These changes are linked to the increase of photodes-orption at larger radii than the water ice line that resultsfrom the reduced dust opacity in the outer disk, and higherradiative flux along the midplane. In the bottom 2 panels of Figure 1 we compare the abun-dances of N and NH . As mentioned before (or see Fig-ure 7 below), there is a higher dust surface density withinthe water ice line of the CPB17 model when compared tothe CPA16 model. This enhancement leads to higher abun-dance of nitrogen carriers within the water ice line because ofthe electron capture and subsequent dissociation of ions likeN H + and NH +4 (Herbst & Klemperer 1973) are catalyzedby the presence of dust. Similarly with a higher dust surface density at smaller radii, the gas is shielded from destructivehigh energy radiation. This dust enhancement is not seen inthe CPA16 dust model because of the simpler treatment ofdust physics and as a result we find that NH is destroyedas the disk ages.An important caveat to our derived NH abundances isthat we currently ignore the formation of NH ice throughsuccessive capture of atomic hydrogen by atomic nitrogen onthe surface of grains (see for ex. Watson & Salpeter (1972)).This process is analogous to the formation of water ice ondust grains, which is included in our chemical model. Theextra NH ice could then contribute to its gas abundancethrough photodesorption and thermal desorption. This pos-sible extra source of NH formation implies that our pre-sented abundances here may represent a lower bound of pos-sible observed abundances. In the middle panels of Figure 2 we show the enhancementof chemical tracers of disk ionization (HCO + and N H + )relative to the chemical model presented in CPA16. This en-hancement confirms the high ionization rate that is presentin the CPB17 disk model. In the bottom right panel of Fig-ure 2 we show the abundance of HCN while in top right panelof Figure 2 we show the abundance of CN. Observationally,the relative abundance of these species has often been usedas a proxy for UV flux. In this work we see a reduction inthe abundance of HCN relative to the results of CPA16, andthe corresponding enhancement of CN which suggests thatUV photolysis is more prevalent in the CPB17 disk modelthan in our previous model. This prevalence is simply dueto the reduced disk opacity and increased UV flux along themidplane. We will elaborate on the link between the flux ofhigh energy and the abundances of HCN and CN in moredetail in the following subsection.At late times in the CPB17 disk model (red dotted line)tracers that are formed from photochemical processes (CN,HCO + , N H + , and H CO) follow similar radial distribu-tions: low abundances within the water ice line and highabundance outside. At late times photochemistry dominatesthe productions of these chemical species, and hence theirabundance is directly linked to the flux of high energy pho-tons on the midplane. This fact can most easily be seen inHCO + which shows a dramatic increase in its abundanceoutside the water ice line, where high energy photons canfirst penetrate down to the midplane. N H + is largely keptunder-abundant by destructive reactions with gaseous CO(Qi et al. 2013, 2015). At the largest radii, outside of theCO ice line, N H + is recovered in all models when the COgas freezes out. As an example of the subtle link between dust physics, ra-diative transfer, and chemistry we discuss the changes inthe abundance of HCN and CN at different disk radii. Therelative abundance of HCN and CN is often used as an ob-servational tracer of UV flux in protoplanetary disks (forex. Kastner et al. (1997)) because of the chemical pathway c (cid:13) , 1–19 A.J. Cridland et al. (R/AU)357911 l og ( n X / c m − ) CO t = 0.1 Myr t = 1.3 Myr t = 3.7 Myr (R/AU)-4-2024681012 l og ( n X / c m − ) H O0 1 2log (R/AU)-8-4048 l og ( n X / c m − ) CO (R/AU)-3-113579 l og ( n X / c m − ) CH (R/AU)-8-4048 l og ( n X / c m − ) NH (R/AU)-11357911 l og ( n X / c m − ) N Figure 1.
Radial distribution of the chemical species that are important to the observations of exoplanetary atmospheres for the CPA16disk model (solid lines) and for the results of the CPB17 dust model presented in this work (dashed lines). The coloured arrows denotethe approximate location of the water ice line in the CPB17 model, defined as the point where the water vapour abundance has droppedby half. c (cid:13) , 1–19 omposition of Early Planetary Atmospheres II (R/AU)-17-13-9-5-13 l og ( n X / c m − ) CH OH t = 0.1 Myr t = 1.3 Myr t = 3.7 Myr (R/AU)-6-4-202468 l og ( n X / c m − ) CN0 1 2log (R/AU)-7-5-3-113 l og ( n X / c m − ) HCO+ 0 1 2log (R/AU)-9-7-5-3-1 l og ( n X / c m − ) N H+0 1 2log (R/AU)-8-6-4-2024 l og ( n X / c m − ) H CO 0 1 2log (R/AU)-4-202468 l og ( n X / c m − ) HCN
Figure 2.
Same as in Figure 1 but for other chemical species that have been observed in protoplanetary disks, and are used as tracersof photochemical processes.c (cid:13) , 1–19 A.J. Cridland et al.
Figure 3.
The radiative flux of the highest energy bands of UV(top panel) and X-ray (middle panel) along with the midplaneabundances of HCN and CN at log t/yr = 5 .
0. On the middlepanel, we plot an estimate of the height where an incoming X-rayphoton reaches an optical depth of one. On the bottom panel,the dotted lines show the midplane location of the 10 − contoursof the UV and X-ray flux. The location of the water ice line isdenoted with the black arrow on the top and middle panels. Onthe contour plots, regions of white space denote zero flux. that produces CN through the photodissocation of HCN(Willacy & Langer 2000).In Figure 3 we relate the midplane abundance of HCNand CN to the flux of high energy photons through the disk.Within the water ice line (denoted by the black arrows) thelarge abundance of dust grains shields the gas from high en-ergy photons allowing ion chemistry, catalyzed by the pres-ence of dust, to dominate the production and destructionof HCN and CN. The ions responsible for these reactionsare the residual ions from our initial chemical state (Table1). Crossing the water ice line results in a strong reductionof dust surface density, and a subsequent reduction in theabundances of both HCN and CN.The radiative flux of both UV and X-ray photons showsimilar features in our disk models. There is a layer of highly ionized (electron fraction ∼ − ) gas far above the mid-plane of the disk. This layer is denoted by the brightest con-tours in Figure 3. In the middle panel of Figure 3 we plot anestimate of the height where photons reach an optical depthof one. Lower than this height the flux of photons comingdirectly from the source (ie. ignoring scattering) is attenu-ated, reaching zero in the white regions of the figure. Thehigh optical depth of the disk’s dust within the water iceline creates a shadowing outside the ice line which is even-tually filled in through scattering. As the disk ages and thedust surface density is reduced, the shadowing shrinks andits edge moves to lower radii. X-ray photons scatter most ef-ficiently, so their flux recovers very quickly just outside thelocation of the water ice line. By contrast, the lower energyUV photons take longer to recover along the midplane. Wenote the difference in the radial location of the 10 − contouras dotted lines on the bottom panel of Figure 3.Moving outward from the water ice line, the abundanceof HCN and CN initially drops because of the reduced sur-face density of dust, but HCN begins to recover when theflux of X-rays along the midplane suddenly increases. Theabundance of CN remains low until the flux of UV pho-tons begins to rise, at which point it is produced by thephotodissociation of HCN. While HCN is destroyed by pho-todissociation, it is efficiently produced through ion chem-istry, maintaining its abundance until the gas is sufficientlycold to allow HCN to freeze out. Finally at the largest radiiwe reach the HCN ice line, where it begins to freeze outonto dust grains. This freeze out also strongly reduces theabundance of CN since there is less HCN in the gas phaseto photodissociate. • We find minimal changes to the location of the waterand CO ice line between the two disk models from CPA16and CPB17 • The CO ice line shows larger variation between diskmodels, caused by the increased flux of UV at the midplane • Volatiles are more abundant in the outer disk of theCPB17 disk model, caused by a radial drift induced reduc-tion of dust surface density in the outer disk and less efficientfreeze out • An enhancement of NH and N is found within thewater ice line of the CPB17 disk model because of the higherdust surface density within the water ice line • The abundance of photochemical tracers (HCO + , HCN,and CN) are indicative of higher UV and X-ray flux on themidplane of the CPB17 disk model, outside of the water iceline • Strong variation in chemical abundances can be linkedto gradients in dust surface density, and radiative flux ofhigh energy photons
In the CPB17 disk model, the evolution of the chemical stateof the disk changes the formation history of dynamicallytrapped planets. Here we compare the planetary atmosphereresults from CPA16 model and the chemical model presentedin this work. c (cid:13) , 1–19 omposition of Early Planetary Atmospheres II -1 0 1 2log (R/AU)-2-10123 l og ( M p l n t / M ⊕ )
11 Myr 22 Myr 33 Myr 44 Myr1 1 Myr2 2 Myr3 3 Myr4 4 Myr 11 Myr2 2 Myr3 3 Myr4 4 Myr
Planet Formation Tracks
Dead zoneIce lineHeat trans (a) Planetary formation tracks from CPA16. -1 0 1 2log (R/AU)-2-10123 l og ( M p l n t / M ⊕ )
11 Myr22 Myr3 3 Myr4 4 Myr1 1 Myr2 2 Myr3 3 Myr44 Myr 11 Myr22 Myr33 Myr44 Myr
Planet Formation Tracks
Dead zoneIce lineHeat trans (b) Planetary formation tracks for this work based on the CPB17 dustmodel.
Figure 4.
Planet formation tracks for the two different disk models presented in CPA16 and CPB17. Each model had the same diskmass, size, gas evolution and planet formation model. The only difference is in the choice of dust model. The grey dashed lines showsthe location of the water ice line for the two disk models at 1,2,3 and 4 Myr. The annotated numbers shows where on the mass andsemi-major axis diagram the planet is at 1,2,3 and 4 Myr. In the CPB17 dust model, the chemical abundance of the disk was computedat a higher resolution than in CPA16. Hence why the water ice line appears to evolve to smaller radii in the right panel.
Table 2.
Final planetary properties at t = t life from the CPA16and CPB17 dust models CPA16 M plnt ( M ⊕ ) a final (AU)Dead zone 0 .
95 3 . . . . . CPB17 M plnt ( M ⊕ ) a final (AU)Dead zone 197 . . . . .
24 0 . In Figures 4(a) and 4(b) we show the evolution of the form-ing planets through the mass and semi-major axis diagramfor the model presented in CPA16 (4(a)) and the model pre-sented in CPB17 (4(b)). We find a significant difference inthe planet formation histories between these two disk mod-els. At the water ice line, the dust density is enhanced in thenew dust model. This enhancement can increase the dust-to-gas ratio up to 1-to-4 from the standard 1-to-100 that isused in CPA16, which drastically increases the rate of solidaccretion in the new dust model. Because of the enhance-ment the planet forming at the water ice line builds a largercore and rapidly draws down an atmosphere, faster thanin the CPA16 model. Formation happening quicker implies that the planet will not migrate inwards as far in the CPB17dust model which results in the larger final semi-major axisfor the planet in the water ice line.The planet forming at the dead zone edge in the CPA16model did not accrete an atmosphere because it spent itswhole formation history far out in the disk. This is in con-trast to the CPB17 model where the dead zone edge rapidlyevolves inward, bringing its forming planet with it. In theCPB17 model the dead zone planet results in a mass of ap-proximately 2 Saturn masses at ∼ .
12 AU.The heat transition planet in the CPB17 model resultsin a super Earth sized planet very close to its host star. Themass of this planet is significantly smaller than its counter-part in the CPA16 model planet which is approximately 1.4Jupiter masses. The difference between these two models isdue to the evolution of the dead zone edge which crosses theposition of the heat transition trap very early in the CPB17model. In this model, the dead zone spends most of its timeat smaller radii than the orbital radius of the heat transi-tion trap. As a result the planet remains trapped at the heattransition trap location for its entire formation history, stay-ing at larger radii than its counterpart in the CPA16 model.At larger radii, the surface density of gas and dust is lowerand hence the planet has less material on which to feed.In Figure 5(a) we show the resulting chemical abun-dances of all minor gases (most abundant gases other thanH and He) produced in the CPA16 model. Comparativelyin Figure 5(b) we show the resulting chemical abundancesfor the planets formed in the CPB17 model. Since the planetat the dead zone edge did not accrete an atmosphere in theCPA16 model, we instead plot the abundance of ice speciesthat were accreted and assumed 100% out-gassing. We donot believe this is an appropriate assumption for the pro-duction of an atmosphere around this planet and hence will c (cid:13) , 1–19 A.J. Cridland et al.
Mass Percent
Dead zoneIce lineHeat trans H OCOHCO NH HCNCH HNCN Mass Percent
Dead zoneIce lineHeat trans
Minor Gas Abundances
Dead zoneIce lineHeat transDead zoneIce lineHeat trans (a) Atmospheric results presented in CPA16.
Mass PercentDead zoneIce lineHeat trans H OCOHCO HCNHNCNH NH CNCH O N Mass PercentDead zoneIce lineHeat trans
Minor Gas Abundances (b) Atmospheric results of this work, based on the dust model pre-sented in CPB17.
Figure 5.
Atmospheric abundances for most abundant gases other than hydrogen and helium. The bottom panels of both plots are themost abundant gases after water and CO are also removed. Note that in the bottom panel of Figure 5(b) the scale is a factor of 5 largerthan in the bottom panel of Figure 5(a). limit our discussion to planets that directly accreted theiratmosphere from the disk gas.In our new dust model the water and CO abundancesare similar because all planets accreted their atmosphericgas near or within the water ice line. As the disk ages andcools the three planet traps converge at small radii and sincethe heat transition and dead zone planets accrete their gaslater in the disk lifetime, all planets accrete gas with similarchemical abundances.The less dominant species show some variation betweenthe two dust models. The heat transition planet accretes itsgas at a very late time at a colder region of the disk. Thisresults in a negligible amounts of methane and elemental hy-drogen, but a small quantity of oxygen gas being accreted.The ice line planet accreted its atmosphere in less than 1Myr in both models, however it accreted its solid core fasterin the CPB17 model, and did not migrate as close to the hoststar as it did in the CPA16 model. Because of the slowermigration the CPB17 ice line planet accreted its gas in aslightly cooler position of the disk compared to its CPA16counterpart which resulted in less methane, but more nitro-gen gas being accreted.
Table 3.
Carbon-to-Oxygen (C/O) and Carbon-to-Nitrogen(C/N) ratios for the planets from CPA16 and this work basedon the dust model from CPB17.
CPA16
C/O C/NDead zone N/A N/AIce line 0.23 47Heat transition 0.28 10
THIS PAPER
C/O C/NDead zone 0.29 4.12Ice line 0.29 4.08Heat transition 0.29 4.17
Two important tracers of planet formation history are shownin Table 3. The elemental ratios of carbon, nitrogen and oxy-gen trace when and where the planetary atmosphere was ac-creted. The carbon-to-oxygen ratio (C/O) indicates whethera planet accretes its atmosphere at smaller orbital radii thanthe water ice line, or a larger radii. Both the CPA16 and c (cid:13)000
Two important tracers of planet formation history are shownin Table 3. The elemental ratios of carbon, nitrogen and oxy-gen trace when and where the planetary atmosphere was ac-creted. The carbon-to-oxygen ratio (C/O) indicates whethera planet accretes its atmosphere at smaller orbital radii thanthe water ice line, or a larger radii. Both the CPA16 and c (cid:13)000 , 1–19 omposition of Early Planetary Atmospheres II Table 4.
Final planetary parameters for the models of higher andlower initial disk masses. M disk, = 0 . M ⊕ M plnt ( M ⊕ ) a final (AU)Dead zone 197 . . . . . . M disk, = 0 . M ⊕ M plnt ( M ⊕ ) a final (AU)Dead zone 203 . . . . . . CPB17 models result in planets that accrete their atmo-sphere close to the water ice line in the disk, and hencethey all have similar C/O which is similar to the disk’s C/O(0.288). This sub-solar C/O results from our choice of initialchemical abundances which we chose to replicate the chem-ical state of molecular clouds (ie. Aikawa & Herbst (1999)).In doing so we assume that the chemical abundance of thegas is not reset during the disk formation process - an effectthat can significantly alter the chemical state of the disk (seefor ex. Eistrup et al. (2016)).The carbon-to-nitrogen ratio is also sensitive to the lo-cation of the planet when it accreted its atmosphere relativeto the ice lines of nitrogen volatiles. While the abundance ofnitrogen carriers (NH , HCN, N ) is dependent on when theplanet accretes its gas, C/N is more dependent on our choiceof disk model. In the CPB17 disk model, more nitrogen (bymass) is accreted which reduces C/N by a factor of between2.5 and 10 relative to the results of the CPA16 model. As a prelude to a full population synthesis model, we presentthe results of two additional disk models with a higher (0.12M ⊙ ) a lower (0.08 M ⊙ ) initial disk mass than our fiducialmodel (0.1 M ⊙ ).The initial disk mass sets the normalization for the ini-tial gas surface density profile for a given initial disk size.Because of its dependence on gas surface density, the initialmass accretion rate is also set by the choice of initial diskmass. In our disk model, the temperature of the inner (r .
20 AU) disk is dependent on viscous heating, and henceis also dependent on the initial disk mass. All of these de-pendencies conspire to change the initial location of each ofthe planet traps, and the formation history for the planetsforming within the traps.In Figures 6(a) and 6(b) we show the planet formationtracks for the disk models with initial mass 0 . M ⊕ and0 . M ⊕ respectively. In the lower mass disk, the heat tran-sition and dead zone planets begin at smaller orbital radiithan the planets in the high mass disk. As a result the plan-ets in the lower mass disk accrete their material in a regionof the disk with higher gas and dust surface density, build-ing the core of the planet faster. In the case of the planettrapped at the heat transition trap, because its core is builtfaster in the lower mass disk it has time to accrete an at- mosphere. While in the higher mass disk model the heattransition planet did not. The ice line planet also accretesits core faster in the low mass disk model, the reason for thisis discussed below.The final masses and semi-major axes are shown in Ta-ble 4. The planet forming at the dead zone trap accretesa similar mass and ends its evolution at a similar semi-major axis in all three models. This feature is related tothe strength of the dust surface density enhancement thatoccurs at the ice line, and the tendency of the dead zoneto evolve towards this enhancement. This tendency stronglyties the evolution of both traps at the latest times in thedisk, and since the planet in the dead zone trap never ac-cretes enough material to open a gap, it is bound to theradial evolution of its trap.At late times the viscous regime shrinks and eventu-ally disappears completely, resulting in a disk that is heatedsolely by direct irradiation from its host star. The resultingtemperature structure does not evolve in time and hencethe location of the ice line and dead zone evolve to nearlythe same radius in every disk model before ceasing its radialevolution. In Figure 7 we show the dust surface density for the threedisk models shown at 0.1 Myr and 1.3 Myr. In the figureswe denote the location of the water ice line for the fiducialdisk model. The water ice line for the two other models isinitially located at smaller orbital radii for the lower massdisk and a larger radii in the higher mass disk. In the outerdisk ( r > r iceline ) the dust surface density scales with thetotal initial disk mass with higher initial disk masses leadingto higher dust surface densities, while in the inner disk ( r 11 Myr22 Myr3 3 Myr4 4 Myr1 1 Myr2 2 Myr33 Myr44 Myr 11 Myr22 Myr33 Myr4 4 Myr Planet Formation Tracks Dead zoneIce lineHeat trans (a) Planetary formation tracks for M disk, = 0 . M ⊕ -1 0 1 2log (R/AU)-2-10123 l og ( M p l n t / M ⊕ ) 11 Myr22 Myr33 Myr4 4 Myr1 1 Myr2 2 Myr3 3 Myr4 4 Myr 11 Myr22 Myr33 Myr44 Myr Planet Formation Tracks Dead zoneIce lineHeat trans (b) Planetary formation tracks for M disk, = 0 . M ⊕ Figure 6. Planet formation tracks for the two additional disk models (see § -1 0 1 2 log Disk radius (AU) -4-3-2-10123 l og S u r f a c eden s i t y ( g / c m ) Water ice line t = 0 . Myr t = 1 . Myr Dust surface density of the three dust models Total M disk = 0 . M ⊕ M disk = 0 . M ⊕ M disk = 0 . M ⊕ (a) -1 0 1 2 log Disk radius (AU) -4-3-2-10123 l og S u r f a c eden s i t y ( g / c m ) t = 0 . Myr t = 1 . Myr Dust surface density of the three dust models a < µ ma < µ m (b) -1 0 1 2 log Disk radius (AU) -4-3-2-10123 l og S u r f a c eden s i t y ( g / c m ) t = 0 . Myr t = 1 . Myr Dust surface density of the three dust models µ m < a < mm µ m < a < mm (c) -1 0 1 2 log Disk radius (AU) -4-3-2-10123 l og S u r f a c eden s i t y ( g / c m ) t = 0 . Myr t = 1 . Myr Dust surface density of the three dust models a > mma > mm (d) Figure 7. Dust surface density for the three presented models. Along with the total surface density (in black) the densities are binned intoseparate sizes: sub-micron, sub-millimeter and sizes greater than a millimeter. These dust surface densities are shown for two snapshotsof our disk model, with 0.1 Myr represented by the top line in all figures, and 1.3 Myr represented by the bottom lines. Dust is clearedbetween the two snapshots because of radial drift and viscous evolution. This clearing is particularly effective on the largest grains (Figure7(d)), which shows a truncation in the disk at larger radii. c (cid:13) , 1–19 omposition of Early Planetary Atmospheres II Table 5. Elemental ratios at t = t life from the low and highinitial mass dust models M disk, = 0 . C/O C/N Dead zone 0 . 288 4 . . 288 4 . . 289 4 . M disk, = 0 . C/O C/N Dead zone 0 . 288 4 . . 288 4 . shown in Figure 6(a). The planet starts its gas accretionwith a heavier core and hence undergoes a shorter period ofslow gas accretion. Because this slow gas accretion phase isresponsible for the majority of semi-major axis evolution forforming planets, the resulting gas giant does not evolve farfrom its initial position.Also shown in the figure is the temporal evolution ofthe dust surface density distributions from the earliest time(0.1 Myr) to an intermediate time (1.3 Myr) through thelifetime of the disk. As the disk ages, radial drift and viscousevolution causes the surface density of the dust to drop atall radii. Interestingly, the outer edge of the large dust graindisk (see Figure 7(d)) shrinks, while the outer edge of theintermediate sized grains ( µm < a < mm , see Figure 7(c))grows. In Figures 8(a) and 8(b) we show the resulting bulk chemicalcomposition of the minor gases in each planet’s atmospheres,and in Table 5 we show the elemental ratios. The water andCO abundances (by mass) of each of the planets are nearlyequal, attributed to the fact that each planet accreted itsgas near the ice line, or at smaller radii. Within the radiusof the water ice line the relative abundance of water vapourand CO gas is constant.The abundances of the nitrogen carriers do vary be-tween the planets, in a similar manor as in the fiducialmodel. As in the fiducial model, in the high mass disk modelthe ice line planet accretes its solid core rapidly, accretingan atmosphere before 1 Myr. Early in the disk’s evolutionit is warmer and the nitrogen content is generally in NH .In the lower mass disk model the initial mass accretion rateis lower and hence the gas begins cooler. Because of its ini-tially cooler state, nitrogen is primarily in molecular gas(N ) rather than NH . This difference is then reflected inthe bulk chemical composition of the water ice line planet.Additionally in the lower mass disk model the water ice lineplanet does not migrate very far before accreting an atmo-sphere, and hence it accretes its gas in a cooler region of thedisk than the water ice line planets in the two other diskmodels.In Figure 9 which is one of the key results of this pa-per, we show the mixing ratio of H O, CO , CO, CH for our predicted atmospheres. Additionally, we show the range of mixing ratios that are inferred by retrieval studies(Madhusudhan & Seager 2009, 2010; Miguel & Kaltenegger2014; Lee et al. 2012; Line et al. 2011, 2014) of atmosphericemission spectra. We see that the H O and CO content inour predicted atmospheres does not vary across our threedisk models. The reason for this similarity is because eachdisk model is initialized by the same chemical abundance,and each planet accreted its atmosphere within the water iceline of their respective disks. These results suggest that asan observing tool, if the atmospheric mixing ratios of H Oan CO are similar to the host star then this would suggestthat the planet accreted its gas within the water ice line.Within the water ice line, most of the carbon and oxygencarriers are found in the gas phase ( ¨Oberg et al. 2011). Sincethere is little to no freeze out of carbon and oxygen carryingmolecules, the chemical composition of the gas is similar tothe initial conditions of the disk. This result could change ifthe majority of carbon and oxygen come from the late ac-cretion of planetesimals or pebbles (for eg. Mordasini et al.(2016)). This late stage of accretion occurs after the evapo-ration of the gas disk and is not explored in this work.Even with the same initial chemical abundances in eachdisk model we find small variations in the mixing ratio ofCO , and large variation in the mixing ratios of CH . Thesevariations are directly linked to the different formation histo-ries, and are sensitive to when and where the planet accretesits gas. These variations do depend both on the initial massof the protoplanetary disk and when the planet accretedits gas, which complicates future interpretation of observedchemical abundances in exoplanetary atmospheres. The dust component of a protoplanetary disk has impor-tant observational and theoretical effects on the structure ofthe disk. We have seen that the dust component of our diskmodels impacts the ionization structure and hence the for-mation history and size of planets forming in planet traps.In the context of exoplanetary atmospheres, different for-mation histories lead to different bulk compositions in theatmospheres of exoplanets. Understanding how differencesin these bulk compositions arise could lead to a method ofbroadly interpreting the formation history of an exoplanetbased on the observation of its atmospheric chemical con-tent.These interpretations are complicated by the chemicaland physical evolution of the atmosphere as the planet ages,and it remains to be seen by how much these global mix-ing ratios are changed when the atmosphere is allowed toevolve. Of particular importance could be the interactionbetween the atmosphere and the core. Heavy elements inthe atmosphere could be injected by out-gassing or throughthe erosion of the core (for ex. Ali-Dib (2017)). We leavethese calculations to future work.The fact that we form Jupiter sized planets very early( t < . c (cid:13) , 1–19 A.J. Cridland et al. Mass PercentDead zoneIce lineHeat trans H OCOHCO HCNHNCNH NH CNCH O N Mass PercentDead zoneIce lineHeat trans Minor Gas Abundances (a) Final bulk composition for M disk, = 0 . M ⊕ Mass PercentDead zoneIce lineHeat trans H OCOHCO HCNHNCNH NH CNCH O N Mass PercentDead zoneIce lineHeat trans Minor Gas Abundances (b) Final bulk composition for M disk, = 0 . M ⊕ Figure 8. Resulting bulk compositions of minor gases for the low and high mass disk models. The atmospheres of each of the planetsformed in these models had similar abundances (by mass) of water and CO because they accreted their atmospheres at similar positionsrelative to the water ice line. The nitrogen content shows some variation based on the timing of atmosphere accretion and the globaltemperature structure of their natal disks. Figure 9. A comparison between the mixing ratios of molecules inferred by observations (red bars) with our theoretically derivedplanetary atmospheres. The three disk models with varying initial disk mass are shown, with the left three points represent the low diskmass, the middle three representing the fiducial mass disk and the right two points showing the high mass disk results. The colour of thepoints denote the planets formed in the dead zone (black), the water ice line (blue), and heat transition (orange) traps. The observationaldata comes from Madhusudhan & Seager (2009, 2010); Miguel & Kaltenegger (2014); Lee et al. (2012); Line et al. (2011, 2014).c (cid:13)000 A comparison between the mixing ratios of molecules inferred by observations (red bars) with our theoretically derivedplanetary atmospheres. The three disk models with varying initial disk mass are shown, with the left three points represent the low diskmass, the middle three representing the fiducial mass disk and the right two points showing the high mass disk results. The colour of thepoints denote the planets formed in the dead zone (black), the water ice line (blue), and heat transition (orange) traps. The observationaldata comes from Madhusudhan & Seager (2009, 2010); Miguel & Kaltenegger (2014); Lee et al. (2012); Line et al. (2011, 2014).c (cid:13)000 , 1–19 omposition of Early Planetary Atmospheres II (eg. Zhang et al. (2015); Okuzumi et al. (2016)). While thecondensation front interpretation of disk gaps like in HL Taumay be more consistent with reality than a planet gap, wenote that the concept of a Jupiter massed planet cannot, atthis stage, be disregarded as it appears to be easily producedat a water ice line that is fed a substantial amount of dustby radial drift. O and CO mixing ratios? For our three models presented here we find that each of theplanets show similar ratios of H O and CO regardless of thedisk model, and their C/O are similar to the ratio in theinitial chemical conditions of the disk. This appears to becaused by the fact that each of the planets presented hereeither accreted their gas close to the location of the waterice line or at smaller radius than the location of the waterice line. The convergence of the planet traps is a featureof the particular disk models presented here, however it isnot clear whether this result will stand until an explorationthrough parameter space can be made.These results contrast with some observations that sug-gest that planets have C/O higher than their host starsBrewer et al. (2017). If we assume that the stellar photo-sphere represents the chemical state of the protoplanetarydisk that the planets formed from, then these observationssuggest that most Hot Jupiters accreted their gas beyond thewater ice line. Our work suggests that the opposite is true forthe case of planet formation in the presence of planet traps.Since all of our planets either accreted their gas within orvery near to the water ice line, thereby accreting ‘pristine’gas - with the same C/O as their disk and host star. We notethat a recent retrieval survey for eight hot Jupiters reportsupper limits of C/O < . > r iceline ).This latter feature of our model is inconsistent with sub-millimeter observations of protoplanetary disks, which haveobserved dust out to hundreds of AU. Such an issue as raisedin CPB17 where we find that radial drift clears out a sub-stantial amount of dust quickly in the outer disk. A potentialfix for such a rapid clearing of material is dust trapping at pressure maxima (Pinilla et al. 2012). Indeed in our planetformation model we assume that planet traps act as thesedust traps, and we enhance the density of solids accordingly.However in our model of dust evolution these traps are notimplemented. A possible outcome of these traps that moresolid material is maintained in the outer disk, thereby in-creasing the initial formation rate of planetary cores for theplanets trapped at the dead zone and heat transition.Additionally, the location and evolution of the deadzone may have an important impact on the evolution of thedust size and density distributions. As was shown in CPB17the dead zone suppresses fragmentation and enhanced dustsettling, increasing the density of large grains and reducingthe density of sub-micron grains along the midplane of thedisk. In a future work we will explore a model of co-evolvinggas and dust to incorporate the effects of an evolving deadzone over then entire disk lifetime. as a tracer of formationhistory? The carbon-carrying molecule that showed the greatest levelof variation between each of our modeled planets was CH .Generally, we find that the planets which accreted their at-mosphere in the colder parts of the disk showed the smallestmixing ratios for this molecule. This implies that an observa-tion of a small CH mixing ratio could imply that the planetformed near or outside the water ice line. We note, howeverthat CH tells you very little about when a planet accretesits atmosphere. As an example in the disk with initial massof 0.08 M ⊕ the ice line and dead zone planets accreted atvery different times, but resulted with nearly identical CH mixing ratios.These conclusions will become complicated by the equi-librium chemistry that occurs in exoplanetary atmospheres.In equilibrium chemistry, there is a sharp transition at atemperature of ∼ 750 K above which CH is converted intoCO (see for ex. Pignatale et al. (2011)). This transition willgreatly impact our ability to observe CH in Hot Jupitersbecause their equilibrium temperatures can exceed 1000 K(Wakeford et al. 2017). With the James Webb Space Telescope (JWST) launch in2018 we will have a new tracer of formation history, NH .The camera JWST-MIRI will study the mid-infrared whichgives us the first chance to directly detect features caused bythe presence of NH in the emission spectra of atmospheres.The first possible detection of nitrogen chemistry has beenrecently reported by MacDonald & Madhusudhan (2017).In that work, the authors report a new retrieval technique(known as POSEIDON) which they apply to the spectra ofHD 209458b. MacDonald & Madhusudhan (2017) report arange of mixing ratios for NH between 0.01 - 2.7 ppm rela-tive to the abundance of hydrogen atoms, depending on theirchoice of atmospheric model. Across all of our disk modelspresented here, we find mixing ratios as low as 2 × − ppm,and as high as 52 ppm. These results results suggest thatwith current HST observations and advanced retrieval algor-thims (eg. Lavie et al. (2016); MacDonald & Madhusudhan c (cid:13) , 1–19 A.J. Cridland et al. (2017)) we can begin to constrain the formation histories ofobserved planets with high abundances of NH .As was pointed out in CPA16, the detection of NH in an atmosphere might be indicative of a planet that ac-creted its gas early ( t < does not necessarily indicate a planetthat accreted gas later, because some atmospheric evolu-tion could have removed the molecule from the upper atmo-sphere. In this work we have demonstrated that complex dust phys-ical models can drastically change the formation history ofplanets forming in planet traps. With the changing forma-tion history comes a change in when and where the planetaccretes its atmosphere and hence a change in the bulk at-mospheric abundances of the gas.Changing the surface density of dust grains throughoutthe disk impacts the freezing efficiency of volatiles and theirformation pathways which are dependent on the presence ofdust. This has lead to: • A retention of volatile H O and CO when comparedto the CPA16 (constant dust-to-gas ratio) chemical modelat r > r iceline • An enhancement of nitrogen carriers N and NH arefound at r < r iceline due to formation reactions that arecatalyzed by the presence of dustThe radial drift of dust results in a rapid depletion ofdust in the outer parts of the disk, and a higher ionizationfraction at larger radii ( r > r iceline ). These higher ionizationfractions change the formation history of planets forming inour model producing: • Hot Saturns that form at the dead zone edge • Earth and Super-Earth sized planets at the heat tran-sition trap • Hot Jupiters and 1 AU Jupiters from the water ice linetrap We have begun to explore the available parameter spaceof disk initial masses in order to understand what range ofmolecular abundances could exist in an ‘early’ exoplanetaryatmosphere. This early study has suggested that: • The mixing ratio of H O and CO appear to be con-stant across planets that formed within the water ice line,and produce C/O that are the same as the ratio of the pro-toplanetary disk • CH shows a large variation between different formationhistories and is tied to the planet’s proximity to the disk’swater ice line • NH could act as a measure of when the planet accretesit’s gas– Earlier on, the disk is warmer and NH dominatesthe nitrogen carriers– Later the disk cools and the nitrogen is found pre-dominately in molecular nitrogen gasUnderstanding the differences in the abundance of ele-mental carriers is important as next generation of telescopes come online and begin to observe the emission spectra of ex-oplanetary atmospheres. As the library of detectable chemi-cal species grow our models can evolve to understand wherethe observable diversity of chemical species arise. ACKNOWLEDGEMENTS REFERENCES Aikawa Y., Herbst E., 1999, A&A, 351, 233Alessi M., Pudritz R. E., Cridland A. J., 2016, MNRASAli-Dib M., 2017, MNRAS, 464, 4282Bai X.-N., Stone J. M., 2013, ApJ, 769, 76Bai X.-N., Stone J. M., 2016, ArXiv e-printsBalbus S. A., Hawley J. F., 1991, ApJ, 376, 214Benneke B., 2015, ArXiv e-printsBethell T. J., Bergin E. A., 2011, ApJ, 739, 78Birnstiel T., Andrews S. M., Pinilla P., Kama M., 2015,ApJL, 813, L14Birnstiel T., Dullemond C. P., Brauer F., 2010, A&A, 513,A79Birnstiel T., Klahr H., Ercolano B., 2012, A&A, 539, A148Brewer J. M., Fischer D. A., Madhusudhan N., 2017, AJ,153, 83Chambers J. E., 2009, ApJ, 705, 1206Cleeves L. I., Bergin E. A., Adams F. C., 2014, ApJ, 794,123Cleeves L. I., Bergin E. A., Qi C., Adams F. C., ¨ObergK. I., 2015, ApJ, 799, 204Cridland A. J., Pudritz R. E., Alessi M., 2016, MNRAS,461, 3274Cridland A. J., Pudritz R. E., Birnstiel T., 2017, MNRAS,465, 3865Draine B. T., 2003, ApJ, 598, 1026 c (cid:13) , 1–19 omposition of Early Planetary Atmospheres II Dullemond C. P., 2012, RADMC-3D: A multi-purpose ra-diative transfer tool. Astrophysics Source Code LibraryEistrup C., Walsh C., van Dishoeck E. F., 2016, A&A, 595,A83Fogel J. K. J., Bethell T. J., Bergin E. A., Calvet N., Se-menov D., 2011, ApJ, 726, 29Gressel O., Turner N. J., Nelson R. P., McNally C. P., 2015,ApJ, 801, 84Gundlach B., Blum J., 2015, ApJ, 798, 34Hasegawa Y., Pudritz R. E., 2011, MNRAS, 417, 1236Hasegawa Y., Pudritz R. E., 2012, ApJ, 760, 117Hasegawa Y., Pudritz R. E., 2013, ApJ, 778, 78Helling C., Woitke P., Rimmer P. B., Kamp I., Thi W.-F.,Meijerink R., 2014, Life, 4Herbst E., Klemperer W., 1973, ApJ, 185, 505Ida S., Lin D. N. C., 2004, ApJ, 604, 388Ikoma M., Nakazawa K., Emori H., 2000, ApJ, 537, 1013Jin M., Lee J.-E., Kim K.-T., Evans, II N. J., 2016, ApJS,225, 21Kastner J. H., Huenemoerder D. P., Schulz N. S., CanizaresC. R., Weintraub D. A., 2002, ApJ, 567, 434Kastner J. H., Huenemoerder D. P., Schulz N. S., Wein-traub D. A., 1999, ApJ, 525, 837Kastner J. H., Zuckerman B., Weintraub D. A., ForveilleT., 1997, Science, 277, 67Lavie B. et al., 2016, ArXiv e-printsLee J.-M., Fletcher L. N., Irwin P. G. J., 2012, MNRAS,420, 170Line M. R., Knutson H., Wolf A. S., Yung Y. L., 2014, ApJ,783, 70Line M. R., Vasisht G., Chen P., Angerhausen D., YungY. L., 2011, ApJ, 738, 32MacDonald R. J., Madhusudhan N., 2017, ArXiv e-printsMadhusudhan N., Ag´undez M., Moses J. I., Hu Y., 2016,Space Sci. Rev.Madhusudhan N., Amin M. A., Kennedy G. M., 2014a,ApJL, 794, L12Madhusudhan N., Amin M. A., Kennedy G. M., 2014b,ApJL, 794, L12Madhusudhan N., Seager S., 2009, ApJ, 707, 24Madhusudhan N., Seager S., 2010, ApJ, 725, 261Masset F. S., Morbidelli A., Crida A., Ferreira J., 2006,ApJ, 642, 478Matsumura S., Pudritz R. E., Thommes E. W., 2009, ApJ,691, 1764Miguel Y., Kaltenegger L., 2014, ApJ, 780, 166Mordasini C., Molli`ere P., Dittkrist K.-M., Jin S., AlibertY., 2015, International Journal of Astrobiology, 14, 201Mordasini C., van Boekel R., Molli`ere P., Henning T., Ben-neke B., 2016, ArXiv e-printsMoses J. I., Madhusudhan N., Visscher C., Freedman R. S.,2013, ApJ, 763, 25¨Oberg K. I., Bergin E. A., 2016, ArXiv e-prints¨Oberg K. I., Murray-Clay R., Bergin E. A., 2011, ApJL,743, L16Okuzumi S., Momose M., Sirono S.-i., Kobayashi H.,Tanaka H., 2016, ApJ, 821, 82Ormel C. W., Kuiper R., Shi J.-M., 2015, MNRAS, 446,1026Paardekooper S.-J., 2014, MNRAS, 444, 2031Paardekooper S.-J., Baruteau C., Crida A., Kley W., 2010,MNRAS, 401, 1950 Pignatale F. C., Maddison S. T., Taquet V., Brooks G.,Liffman K., 2011, MNRAS, 414, 2386Pinilla P., Birnstiel T., Ricci L., Dullemond C. P., UribeA. L., Testi L., Natta A., 2012, A&A, 538, A114Qi C., ¨Oberg K. I., Andrews S. M., Wilner D. J., BerginE. A., Hughes A. M., Hogherheijde M., D’Alessio P., 2015,ApJ, 813, 128Qi C. et al., 2013, Science, 341, 630Ranjan S., Charbonneau D., D´esert J.-M., MadhusudhanN., Deming D., Wilkins A., Mandell A. M., 2014, ApJ,785, 148Sch¨afer U., Yang C.-C., Johansen A., 2017, A&A, 597, A69Seager S., 2010, ExoplanetsShakura N. I., Sunyaev R. A., 1973, A&A, 24, 337Sing D. K. et al., 2015, MNRAS, 446, 2428Thiabaud A., Marboeuf U., Alibert Y., Leya I., Mezger K.,2015, A&A, 574, A138Wada K., Tanaka H., Suyama T., Kimura H., YamamotoT., 2009, ApJ, 702, 1490Wakeford H. R., Sing D. K., 2015, A&A, 573, A122Wakeford H. R., Visscher C., Lewis N. K., Kataria T., Mar-ley M. S., Fortney J. J., Mandell A. M., 2017, MNRAS,464, 4247Ward W. R., 1997, ApJL, 482, L211Watson W. D., Salpeter E. E., 1972, ApJ, 174, 321Willacy K., Langer W. D., 2000, ApJ, 544, 903Zhang K., Blake G. A., Bergin E. A., 2015, ApJL, 806, L7This paper has been typeset from a TEX/ L A TEX file preparedby the author. APPENDIXAPPENDIX A: CONSISTENT TREATMENTOF ENVELOPE OPACITY As an aside we parameterize the critical mass above whichgas accretion begins, as well as gas accretion rate so that thetwo physical processes are assuming the same dust opacityof 3.33 cm /g in the planetary envelope. This parameteriza-tion differs from our previous work because historically theparameters are treated separately due to the fact that theopacity scaling of the critical mass and gas accretion ratechange for different opacity tables. We noted in the textthat the choice of parameter greatly change the formationhistory, and hence will have an impact on the resulting bulkchemical composition of the atmosphere.In Figures 1(a) and 1(b) we show the formation trackof the fiducial version of the CPB17 disk model and thesame disk model with the parameter f crit = 1 . 6. Here we seethat the formation history changes drastically for the largestplanets because the timing of their gas accretion changes.Importantly the core mass has to grow to a larger massbefore it can begin to bring down an atmosphere.Because of its importance on the formation history ofthe planet, we will incorporate variations in the envelopeopacity in our future population synthesis models. c (cid:13) , 1–19 A.J. Cridland et al. -1 0 1 2log (R/AU)-2-10123 l og ( M p l n t / M ⊕ ) 11 Myr22 Myr33 Myr4 4 Myr1 1 Myr2 2 Myr3 3 Myr44 Myr 11 Myr22 Myr33 Myr44 Myr Planet Formation Tracks Dead zoneIce lineHeat trans (a) Planetary formation tracks for the fiducial CPB17 model -1 0 1 2log (R/AU)-2-10123 l og ( M p l n t / M ⊕ ) 11 Myr22 Myr33 Myr44 Myr 1 1 Myr22 Myr33 Myr 44 Myr 11 Myr22 Myr33 Myr44 Myr Planet Formation Tracks Dead zoneIce lineHeat trans (b) Planetary formation tracks for CPB17 model, with f crit = 1 . Figure A1. Comparison of the formation tracks for the fiducialCPB17 disk model with the CPB17 disk model with the planetformation parameter f crit = 1 . 6. c (cid:13)000