Planet-forming material in a protoplanetary disc: the interplay between chemical evolution and pebble drift
MMNRAS , 1–14 (2019) Preprint 31 May 2019 Compiled using MNRAS L A TEX style file v3.0
Planet-forming material in a protoplanetary disc: theinterplay between chemical evolution and pebble drift
R. A. Booth (cid:63) and J. D. Ilee † Institute of Astronomy, Madingley Road, Cambridge CB3 0HA, UK School of Physics & Astronomy, University of Leeds, Leeds LS2 9JT, UK
Accepted 2019 May 27. Received 2019 May 17; in original form 2019 April 10
ABSTRACT
The composition of gas and solids in protoplanetary discs sets the composition of plan-ets that form out of them. Recent chemical models have shown that the compositionof gas and dust in discs evolves on Myr time-scales, with volatile species disappearingfrom the gas phase. However, discs evolve due to gas accretion and radial drift of duston time-scales similar to these chemical time-scales. Here we present the first modelcoupling the chemical evolution in the disc mid-planes with the evolution of discs dueto accretion and radial drift of dust. Our models show that transport will always over-come the depletion of CO from the gas phase, and can also overcome the depletion ofCO and CH unless both transport is slow (viscous α (cid:46) − ) and the ionization rateis high ( ζ ≈ − ). Including radial drift further enhances the abundances of volatilespecies because they are carried in on the surface of grains before evaporating left attheir ice lines. Due to large differences in the abundances within 10 au for models withand without efficient radial drift, we argue that composition can be used to constrainmodels of planet formation via pebble accretion. Key words: astrochemistry — protoplanetary discs — planets and satellites: forma-tion — planets and satellites: composition
The structure and composition of protoplanetary discs arefundamental pieces in the puzzle of how planets form andevolve. In the most direct sense, planets form via the ac-cretion of gas and solids – in the form of planetesimals,pebbles, or dust, which together with the gas are the maincomponents of protoplanetary discs. This provides the op-portunity to learn about planet formation by connecting thecomposition of planets to the discs in which they form. Achallenge for the protoplanetary disc community is thus todetermine the composition of the planetary building blocks,a task complicated by uncertainties in the processes govern-ing disc evolution and the associated difficulties with the in-terpretation observational constraints. Here, we explore oneaspect of this problem – namely how differences in the trans-port of gas and solids through the disc compete with theirchemical evolution to control the composition of discs.We frame our investigation in terms of a recent idea,which is to use the atmospheric C/O ratio in hot Jupitersto determine their building blocks. The observational utility (cid:63)
E-mail: [email protected] (RAB) † E-mail: [email protected] (JDI) of the C/O ratio comes from the relative ease with whichit can be constrained in hot Jupiter atmospheres. This isbecause the abundance of the observed molecules can varyby orders of magnitude for only small changes in the C/Oratio (see Lodders 2010; Madhusudhan 2012). ¨Oberg et al.(2011) realised that the C/O ratio could be used to con-strain the radial location within a disc from which a planetformed. This is possible because the C/O ratio of the gasand ices changes across the ice lines (places where molecularspecies transition from being predominantly in the gas phaseto ices in the solid phase) of the dominant carbon and oxy-gen bearing species, such as CO, CO and H O. Numerousstudies have since explored how planet formation and migra-tion affect the C/O ratio of the planet (Madhusudhan et al.2014; Mordasini et al. 2016; Cridland et al. 2016, 2017; Ali-Dib 2017; Booth et al. 2017; Madhusudhan et al. 2017). Al-though there are differences in the composition of the modelplanets in these studies, there is general agreement that thepartitioning of the carbon and oxygen between the gas andsolid phases, along with the amount these species accreted,is important in determining the planet’s composition.However, the composition of discs is not static in time,which needs to be taken into account in planet formationmodels. Although Cridland et al. (2016, 2017) recently cou- © a r X i v : . [ a s t r o - ph . E P ] M a y Booth & Ilee pled a planet formation model to a disc model includingchemical kinetics, they focused on planets that accreted theirenvelopes inside the water ice line, where the majority ofvolatile carbon and oxygen bearing species are in the gasphase. Further out in the disc, in regions where the tem-peratures are lower and molecular species freeze out, chem-ical kinetics can have a greater impact on the gas and solidphase composition by exchanging carbon and oxygen be-tween species that are in ices on grains or in the gas phase.Using an extensive gas-grain network Eistrup et al. (2016)showed that CH and CO can be removed from the gasphase, with the net effect of lowering the gas-phase C/Oratio inside of around 10 au. However, this process is slow,taking several million years (Eistrup et al. 2018; Bosmanet al. 2018b), by which time gas giant must be well under-way (Greaves & Rice 2010; Najita & Kenyon 2014; Manaraet al. 2018).By comparison, it is becoming clear that the transportof molecular species is important on time-scales comparableto chemical processes. Bosman et al. (2018a) showed thateven conservatively slow models of transport were enoughto raise the CO abundance above levels that has previouslybeen assumed. Furthermore, transport is an essential com-ponent of a popular new model for planet formation – peb-ble accretion (Ormel & Klahr 2010; Lambrechts & Johansen2012) – which relies on dust grains large enough that theydecouple from the gas (e.g. Lambrechts & Johansen 2014;Morbidelli et al. 2015; Bitsch et al. 2015. Johansen et al.2018 invoke smaller pebbles, although they are still largeenough to migrate). These pebbles thus migrate rapidly to-wards the star (Weidenschilling 1977) carrying the volatilemolecular species frozen out onto their surfaces with them asthey migrate. These volatile species then enter the gas phasewhen the pebbles cross their respective ice lines. Under theconditions often assumed in pebble accretion models, radialdrift enhances the gas phase abundances by a factor of a fewwithin a Myr (Booth et al. 2017; Krijt et al. 2018).Previous studies investigating the interplay betweenchemical kinetics and transport typically assume that dustand gas are coupled (e.g. Aikawa et al. 1999; Tscharnuter& Gail 2007; Nomura et al. 2009; Semenov & Wiebe 2011;Heinzeller et al. 2011; Walsh et al. 2014; Gail & Trieloff 2017;see Henning & Semenov 2013, for a review). These stud-ies suggest that the transport has two predominant effects:1) diffusion weakens concentration gradients and 2) trans-port limits the amount of time molecular species spend at agiven location, reducing the influence of chemical processesthat act on time scales longer than the transport time-scale.While a number of studies have treated the transport of gasand solids separately (e.g. in the context of CO sequestra-tion in disc mid-planes; Piso et al. 2015, 2016; Bergin et al.2016; Kama et al. 2016; Booth et al. 2017; Krijt et al. 2018;or the destruction of carbon grains in the terrestrial planet-forming region; Klarmann et al. 2018), to our knowledge thisis the first study to treat gas and dust transport indepen-dently and gas-phase kinetics in detail. We emphasise thatdifferences in the transport of gas and dust is the only waythat the bulk atomic abundances (i.e. gas+dust abundanceof carbon, oxygen, and nitrogen, etc) can vary in the disc;the partitioning of different molecular species between the gas and ice phase changes the atomic abundance of the eachphase separately, but not the total abundance.In this work we investigate the competition betweenchemical kinetics and transport in the mid-plane of proto-planetary discs. To this end, we couple a chemical kineticsnetwork (Ilee et al. 2011) with a 1D disc evolution model.Our disc evolution model includes gas evolution (which istreated as viscous), grain growth and radial drift (Boothet al. 2017) . We allow for differences in the transport ofgas and ice phase molecular species, tracking their motionwith the gas and dust, respectively. We examine how trans-port and chemical kinetics compete, focusing on the classicalgiant planet forming region (1–10 au). We consider the physical evolution of gas and dust follow-ing Booth et al. (2017). The gas evolution is treated as-suming a viscous disc, with a constant α = − or − .The justification of these choices is given in subsection 2.2.Grain growth and radial drift are treated based upon Birn-stiel et al. (2012). The mid-plane temperature is computedtaking into account viscous heating and irradiation from thecentral star, along with external irradiation with a temper-ature of
10 K .We self-consistently take into account the motion ofeach of the chemical species included in the model, treat-ing the motion of gas and ice phase species independently.The gas phase species are advected at the gas velocity, whichis determined by viscous evolution. Ice phase species are ad-vected along with the dust. The dust radial velocity is setby a combination of both the gas velocity due to viscosityand radial drift (e.g. Takeuchi & Lin 2002).The grain growth model assumes that there are twopopulations of grains, small and large, here we assume thatthe ice species are partitioned across the two populations inthe same way as the dust. This results in the advection ofice phase species being dominated by the large grain popu-lation, in good agreement with calculations treating the fulldistribution of sizes (Stammler et al. 2017). In addition toadvection, we include turbulent diffusion. The diffusion co-efficient D of the gas is given by D = ν / Sc , where ν is theviscosity and the Sc the Schmidt number, for which we take Sc = . For the dust and ice phases the diffusion coefficientis scaled to the gas coefficient following Youdin & Lithwick(2007) and assuming the eddy turnover time-scale equalsthe dynamical time-scale, Ω − . For the dust sizes consideredhere, the diffusion coefficients of the dust and gas are thesame to within one per cent.The opacity used in the mid-plane temperature calcu-lation has been updated. In Booth et al. (2017) we used theRosseland mean opacity tables similar to Bell & Lin (1994),as computed by Zhu et al. (2012). Since these opacitiesassume a grain size distribution appropriate for the inter-stellar medium they can considerably overestimate the dustopacity once the grains have grown and begin to migrate. https://github.com/rbooth200/DiscEvolution MNRAS , 1–14 (2019) hemistry & transport in discs Here we instead use Rosseland mean opacity computed self-consistently for grain size distributions with number density n ( a ) ∝ a − . with the maximum grain size taken from thegrain growth model. The dust properties follow Tazzari et al.(2016), based on the composition of Pollack et al. (1994) andassuming a porosity of 30 per cent.In addition the models that include the effects of trans-port processes: viscous evolution, radial drift and diffusion;we also consider models that switch off some or all of theseprocesses. This allows us to separate out the effects of phys-ical and chemical processes on the composition. We explore two models for the evolution of the disc, withthe same initial surface density profile, but choices of theturbulent α parameter – and correspondingly different ac-cretion rates. Unless otherwise stated, in each model we setthe initial disc mass to . M (cid:12) and stellar mass to M (cid:12) .The surface density profile is initialised to a Lynden-Bell &Pringle (1974) profile, Σ ( R ) = M d π R d exp (cid:18) − RR d (cid:19) . (1)We assume that the initial ratio of refractory dust-to-gassurface density is 0.01 everywhere. The total dust-to-gas ra-tio is higher than this by approximately factor of 2 due tothe presence of ices frozen onto grain surfaces. The initialgrain size is also taken to be 0.1 µ m.For the turbulent α parameter we use α = − and α = − . While we note that recent observations suggestturbulence in the outer regions of discs is weaker than this(Pinte et al. 2016; Flaherty et al. 2017, 2018; Teague et al.2018), there is evidence that the accretion rate through discsis more compatible with higher values of α . In the absenceof direct measurements of the radial flow of gas throughdiscs, the accretion rate onto the star combined with the discmass provide the most reasonable estimates. For example,Boneberg et al. (2016) found that an effective α ∼ . wasrequired to match the accretion rate of HD 163296, while theturbulent mixing in the vertical direction was lower. Further-more, based on the surface density profile and accretion rateClarke et al. (2018) derive α ∼ . for CI Tau. Even TWHya, which is an old system and has a low accretion rate,has an accretion rate compatible with α ∼ − (Calvet et al.2002; Bergin et al. 2013; Ercolano et al. 2017). Based on dustdisc sizes and ages, Andrews & Williams (2007) and Guil-loteau et al. (2011) suggest α in the range − – − . Sim-ilarly, for most discs Rafikov (2017) inferred α in the range − – − based on the disc radii and masses inferred fromdust emission in Lupus and their associated accretion rates.High dust-to-gas ratios in Lupus, as suggested by Miotelloet al. (2017), would result in higher estimates of α for Lupus;however, Manara et al. (2016) suggest that the dust mass iswell matched to the accretion rate (but see Mulders et al.2017 and Lodato et al. 2017 for a discussion). Differences be-tween the accretion rate and inferred turbulence levels maypoint to MHD winds as the source of angular momentumtransport instead of viscosity (Suzuki & Inutsuka 2009; Fro-mang et al. 2013; Bai & Stone 2013). However, since we aremainly concerned with the speed of mass transport through the disc, this distinction is not critical. Thus we model thediscs as viscous, considering α in the range − – − to betypical for nearby, low mass protoplanetary discs.For the dust evolution, we use the same parameters asBooth et al. (2017), for which the standard parameters arebased on fits to more detailed models (Birnstiel et al. 2012),except that we use the ice fraction to specify the locationwhere the fragmentation velocity transitions from water ice-like (
10 m s − ) to silicate-like ( − ). Rather than varyingthese parameters, we choose instead to simply run modelswith radial drift included or neglected (in which case thedust moves with the gas). Combined with the two choicesof α , this spans the commonly studied parameter space. Forexample α = − is commonly considered (e.g. Birnstielet al. 2012), for which dust grows large enough to be inthe radial drift dominated regime at large radii, leading torapid evolution of the dust-to-gas ratio. This model is alsocomparable to the parameters assumed in pebble accretionmodels of Lambrechts & Johansen (2014); Morbidelli et al.(2015) and Bitsch et al. (2015). In the α = − model, dustgrowth is initially limited by fragmentation everywhere, re-sulting in smaller grains and slower radial drift. Slower radialdrift combined with a faster gas radial velocity results in thedust-to-gas ratio decreasing more slowly. Such a model maybe in better agreement with observed disc properties, wherelow dust-to-gas ratios are typically not inferred (Boneberget al. 2016; Ansdell et al. 2016; although note that massesdetermined from CO observations are uncertain), and is alsocomparable to the parameters assumed in the pebble accre-tion model of Johansen et al. (2018).The evolution of the surface density, dust-to-gas ratio,temperature, and maximum grain size for our two canonicalmodels are shown in Figure 1 ( α = − ) and Figure 2 ( α = − ). The evolution of combined dust-gas models has beendescribed at length in Birnstiel et al. (2012), and Boothet al. (2017) in the case where adsorption and desorptionat ice lines is included. Thus we only describe the salientdifferences here and refer the readers to the above referencesfor further details. The main features of note are: the slowerevolution of dust-to-gas ratio at α = − , in which the dust-to-gas ratio is enhanced by a factor ∼ inside 10 au, and onlydecreases outside the scaling radius (100 au). In comparisonthe dust-to-gas ratio falls to − on Myr time-scales in the α = − model, a consequence of larger grain maximumgrain size. Note also that gas surface density decreases in the α = − model, but barely changes after 1 Myr at α = − .Another consequence of the higher accretion rate at α = − is the warmer temperature in the disc inside of ∼ , dueto both the larger contribution of viscous heating and alsoa higher optical depth because of smaller grains.The time-scale for the evolution of the gas and dust sur-face densities are given by the time-scale for the accretion ofgas and dust onto the star. These time-scales will be usefulin the following sections for interpreting the competition be-tween chemical evolution and transport. In the case of thegas, the accretion time-scale is controlled by the viscous ve-locity, v ν , from which we can estimate a viscous time-scale, t ν = R / v ν , t ν ≈ × (cid:18) α − (cid:19) − (cid:18) R au (cid:19) yr . (2) MNRAS000
10 m s − ) to silicate-like ( − ). Rather than varyingthese parameters, we choose instead to simply run modelswith radial drift included or neglected (in which case thedust moves with the gas). Combined with the two choicesof α , this spans the commonly studied parameter space. Forexample α = − is commonly considered (e.g. Birnstielet al. 2012), for which dust grows large enough to be inthe radial drift dominated regime at large radii, leading torapid evolution of the dust-to-gas ratio. This model is alsocomparable to the parameters assumed in pebble accretionmodels of Lambrechts & Johansen (2014); Morbidelli et al.(2015) and Bitsch et al. (2015). In the α = − model, dustgrowth is initially limited by fragmentation everywhere, re-sulting in smaller grains and slower radial drift. Slower radialdrift combined with a faster gas radial velocity results in thedust-to-gas ratio decreasing more slowly. Such a model maybe in better agreement with observed disc properties, wherelow dust-to-gas ratios are typically not inferred (Boneberget al. 2016; Ansdell et al. 2016; although note that massesdetermined from CO observations are uncertain), and is alsocomparable to the parameters assumed in the pebble accre-tion model of Johansen et al. (2018).The evolution of the surface density, dust-to-gas ratio,temperature, and maximum grain size for our two canonicalmodels are shown in Figure 1 ( α = − ) and Figure 2 ( α = − ). The evolution of combined dust-gas models has beendescribed at length in Birnstiel et al. (2012), and Boothet al. (2017) in the case where adsorption and desorptionat ice lines is included. Thus we only describe the salientdifferences here and refer the readers to the above referencesfor further details. The main features of note are: the slowerevolution of dust-to-gas ratio at α = − , in which the dust-to-gas ratio is enhanced by a factor ∼ inside 10 au, and onlydecreases outside the scaling radius (100 au). In comparisonthe dust-to-gas ratio falls to − on Myr time-scales in the α = − model, a consequence of larger grain maximumgrain size. Note also that gas surface density decreases in the α = − model, but barely changes after 1 Myr at α = − .Another consequence of the higher accretion rate at α = − is the warmer temperature in the disc inside of ∼ , dueto both the larger contribution of viscous heating and alsoa higher optical depth because of smaller grains.The time-scale for the evolution of the gas and dust sur-face densities are given by the time-scale for the accretion ofgas and dust onto the star. These time-scales will be usefulin the following sections for interpreting the competition be-tween chemical evolution and transport. In the case of thegas, the accretion time-scale is controlled by the viscous ve-locity, v ν , from which we can estimate a viscous time-scale, t ν = R / v ν , t ν ≈ × (cid:18) α − (cid:19) − (cid:18) R au (cid:19) yr . (2) MNRAS000 , 1–14 (2019)
Booth & Ilee -3 -2 -1 § G ; D -4 -3 -2 -1 § D = § G -1 R [au] T [ K ] -1 R [au] -5 -4 -3 -2 -1 a [ c m ] Figure 1.
Evolution of the physical parameters for the model with α = − . Top left: gas (solid) and dust (dashed) surface density. Topright: dust-to-gas-ratio. Bottom left: Temperature. Bottom Right: maximum grain size. -3 -2 -1 § G ; D -4 -3 -2 -1 § D = § G -1 R [au] T [ K ] -1 R [au] -5 -4 -3 -2 -1 a [ c m ] Figure 2.
Same as Figure 1, but for α = − . The long viscous time-scale at 100 au explains the slowchange in gas surface density in α = − model. However,at 1 au the viscous transport time is always much less than1 Myr, and for α = − the viscous transport time is lessthan 1 Myr everywhere inside of approximately 50 au.In addition to transport with the gas, the dust surfacedensity also evolves due to radial drift, with the associatedradial drift time-scale, t drift = R / v drift . In this case, the evolu-tion of the dust surface density is always on the shorter of ra-dial drift and viscous time-scales. The radial drift time-scaledepends on both the dust grain size and gas surface densityprofile. Taking the gas surface density to be Σ = Σ ( R / au ) − ,we estimate the radial drift time-scale as t drift ≈ × (cid:18) Σ g cm − (cid:19) (cid:16) a mm (cid:17) − yr . (3)Typical grain sizes are 0.1–1 mm for α = − and 1–10 mm for α = − . Hence the dust evolves on approxi-mately the viscous time-scale for α = − , which explainsthe slow evolution of dust-to-gas ratio. However, for α = − the dust evolves faster than the gas, and the dust-to-gas ra-tio decreases. Within the physical model discussed above, we embed atime-dependent gas-grain chemical evolution model thatself-consistently calculates the chemical evolution of eachgrid cell. For this, we utilise the krome chemical evolutionpackage (Grassi et al. 2014). For the chemical reaction net-work, we adopt the network previously used in Ilee et al.2011; Evans et al. 2015 and Ilee et al. 2017, with modifica-tions described here to provide consistency with more recentdisc evolution models. The network consists of 1485 reac-tions involving 136 species containing the elements H, He,C, N, O, Na and S. These reactions were originally selectedfrom a subset of the UMIST Rate 95 database (Millar et al. MNRAS , 1–14 (2019) hemistry & transport in discs , were used to update some of the rates and rate co-efficients. In addition to the chemical reactions, adsorption,and desorption processes; we self-consistently follow the ad-vection and diffusion of each of the gas and ice phase speciesfollowing Booth et al. (2017). We overview the key aspectsof our chemical model below.The fundamental variables integrated in the disc evolu-tion and transport code are the total surface density, dustfraction, and mass fraction of each chemical species. Theseare used to compute the mid-plane number density, n ( i ) , anddust-to-gas ratio, (cid:15) , assuming a Gaussian vertical structureusing the scale-height computed from the mid-plane temper-ature. The fractional abundance relative to the total numberof hydrogen nuclei, n , is X ( i ) ≡ n ( i )/ n . The rate equation forthe gas-phase fractional abundance of species i is therefore ddt X ( i ) ≡ S ( i ) (4) = (cid:213) j , l , m k ( j ) X ( l ) X ( m ) n − (cid:213) j (cid:48) , m k ( j (cid:48) ) X ( i ) X ( m ) n − (cid:213) j (cid:48)(cid:48) k ( j (cid:48)(cid:48) ) X ( i ) n − S rad ( i ) X ( i ) + S ( i ) + S d ( i ) − S a ( i ) , where k ( j ) is the rate coefficient of the j th reaction, and thesummations are restricted so that only reactions involved inthe formation or removal of the i th species are included. Wedenote this total rate S ( i ) . These rate coefficients are of thestandard Arrhenius form k ( j ) = α ( j ) (cid:18) T (cid:19) β ( j ) exp (cid:18) − γ ( j ) T (cid:19) , (5)where α ( j ) is the room temperature rate coefficient of thereaction (at K), β ( j ) describes the temperature depen-dence and γ ( j ) is the activation energy of the reaction. Theterms S rad ( i ) X ( i ) , S , S a and S d represent radiative destruc-tion processes, three-body gas phase reactions, adsorptionon to dust grains, and desorption from dust grains, respec-tively.For three body reactions, we assume that the only thirdbody of importance is H , and therefore S ( i ) = (cid:213) j , l , m k ( j ) X ( l ) X ( m ) X ( H ) n (6) − (cid:213) j (cid:48) , m k ( j (cid:48) ) X ( i ) X ( m ) X ( H ) n − (cid:213) j (cid:48)(cid:48) k ( j (cid:48)(cid:48) ) X ( i ) X ( H ) n . The rate of destruction of species i due to radiative pro-cesses is given by S rad = α ( i ) exp (− γ A V ) (7)for photoreactions (where α is a proportionality constantcontaining the dust grain albedo, which we take to be 0.5). http://kida.obs.u-bordeaux1.fr For cosmic ray ionisation, we take S rad = ζ = − s − every-where, unless otherwise specified. Though photoabsorptionof radiation from external sources will affect the chemistry inthe upper layers of a disc, we focus on the bulk of the planetforming material toward the disc mid-plane. We therefore as-sume that all material is well shielded from external sourcesof photons, setting A V = mag everywhere. However, we doinclude photoreactions due to secondary photons generatedby cosmic rays, using the rates from Heays et al. (2017).For the adsorption of species from the gas phase ontothe surfaces of icy grain mantles, we assume S a = X dust S ( i ) (cid:115) kT π m ( i ) (cid:104) π a (cid:105) (8)where X dust is the dust number density, S ( i ) is the stickingcoefficient (taken to be 0.3), m ( i ) is the atomic mass of theadsorbed species, and (cid:104) π a (cid:105) is the average dust grain area.The average grain radius is computed using the maximumgrain size from the dust evolution code, assuming a grainsize distribution n ( a ) d a ∝ a − . d a with the minimum grain-size equal to . µ m. However, we note that our results areonly very weakly dependent on this choice (as shown in sec-tion A).For unsaturated C, N and O species on the surfacesof dust grains, we assume hydrogenation occurs via theEley-Rideal mechanism, at the rate of adsorption of H fromthe gas multiplied by the probability of encounter with thespecies on the grain surface, e.g., gC → gCH → . . . → gCH (where g denotes a species on a grain surface, see Visseret al. 2011).For the thermal desorption of species from dust grainsurfaces into the gas phase, we assume S d ( i ) = . × − σ ν ( i ) exp (cid:18) − E b ( i ) kT (cid:19) , (9)where σ = . × cm − is the surface density of bind-ing sites, E b ( i ) is the binding energy of the i th species onthe surface of the dust grain, and T is the dust temperature(which we assume to be in equilibrium with the gas temper-ature), and ν is the characteristic vibrational frequency ofthe species attached to the grain.Integration of the rate equations in each cell is per-formed independently, using the dlsodes package (Hind-marsh 1983) within the krome package (Grassi et al. 2014)to yield the time-dependent evolution of fractional abun-dance for each species throughout the disc. The coupling to the disc evolution model is achieved in a afirst-order operator-split fashion, which we briefly describehere. The full equation for the surface density of a givenspecies, Σ ( i ) , may be written as ∂ Σ ( i ) ∂ t + R ∂∂ R [ R Σ ( i ) v r ( i )] = Σ ( i ) S ( i ) , (10)where v r ( i ) is the radial velocity of the species (either the gasvelocity or dust velocity, depending on whether the speciesis in the gas or ice phase), and S ( i ) is the source terms due tochemical reactions (Equation 4). We discretize this equation MNRAS000
Same as Figure 1, but for α = − . The long viscous time-scale at 100 au explains the slowchange in gas surface density in α = − model. However,at 1 au the viscous transport time is always much less than1 Myr, and for α = − the viscous transport time is lessthan 1 Myr everywhere inside of approximately 50 au.In addition to transport with the gas, the dust surfacedensity also evolves due to radial drift, with the associatedradial drift time-scale, t drift = R / v drift . In this case, the evolu-tion of the dust surface density is always on the shorter of ra-dial drift and viscous time-scales. The radial drift time-scaledepends on both the dust grain size and gas surface densityprofile. Taking the gas surface density to be Σ = Σ ( R / au ) − ,we estimate the radial drift time-scale as t drift ≈ × (cid:18) Σ g cm − (cid:19) (cid:16) a mm (cid:17) − yr . (3)Typical grain sizes are 0.1–1 mm for α = − and 1–10 mm for α = − . Hence the dust evolves on approxi-mately the viscous time-scale for α = − , which explainsthe slow evolution of dust-to-gas ratio. However, for α = − the dust evolves faster than the gas, and the dust-to-gas ra-tio decreases. Within the physical model discussed above, we embed atime-dependent gas-grain chemical evolution model thatself-consistently calculates the chemical evolution of eachgrid cell. For this, we utilise the krome chemical evolutionpackage (Grassi et al. 2014). For the chemical reaction net-work, we adopt the network previously used in Ilee et al.2011; Evans et al. 2015 and Ilee et al. 2017, with modifica-tions described here to provide consistency with more recentdisc evolution models. The network consists of 1485 reac-tions involving 136 species containing the elements H, He,C, N, O, Na and S. These reactions were originally selectedfrom a subset of the UMIST Rate 95 database (Millar et al. MNRAS , 1–14 (2019) hemistry & transport in discs , were used to update some of the rates and rate co-efficients. In addition to the chemical reactions, adsorption,and desorption processes; we self-consistently follow the ad-vection and diffusion of each of the gas and ice phase speciesfollowing Booth et al. (2017). We overview the key aspectsof our chemical model below.The fundamental variables integrated in the disc evolu-tion and transport code are the total surface density, dustfraction, and mass fraction of each chemical species. Theseare used to compute the mid-plane number density, n ( i ) , anddust-to-gas ratio, (cid:15) , assuming a Gaussian vertical structureusing the scale-height computed from the mid-plane temper-ature. The fractional abundance relative to the total numberof hydrogen nuclei, n , is X ( i ) ≡ n ( i )/ n . The rate equation forthe gas-phase fractional abundance of species i is therefore ddt X ( i ) ≡ S ( i ) (4) = (cid:213) j , l , m k ( j ) X ( l ) X ( m ) n − (cid:213) j (cid:48) , m k ( j (cid:48) ) X ( i ) X ( m ) n − (cid:213) j (cid:48)(cid:48) k ( j (cid:48)(cid:48) ) X ( i ) n − S rad ( i ) X ( i ) + S ( i ) + S d ( i ) − S a ( i ) , where k ( j ) is the rate coefficient of the j th reaction, and thesummations are restricted so that only reactions involved inthe formation or removal of the i th species are included. Wedenote this total rate S ( i ) . These rate coefficients are of thestandard Arrhenius form k ( j ) = α ( j ) (cid:18) T (cid:19) β ( j ) exp (cid:18) − γ ( j ) T (cid:19) , (5)where α ( j ) is the room temperature rate coefficient of thereaction (at K), β ( j ) describes the temperature depen-dence and γ ( j ) is the activation energy of the reaction. Theterms S rad ( i ) X ( i ) , S , S a and S d represent radiative destruc-tion processes, three-body gas phase reactions, adsorptionon to dust grains, and desorption from dust grains, respec-tively.For three body reactions, we assume that the only thirdbody of importance is H , and therefore S ( i ) = (cid:213) j , l , m k ( j ) X ( l ) X ( m ) X ( H ) n (6) − (cid:213) j (cid:48) , m k ( j (cid:48) ) X ( i ) X ( m ) X ( H ) n − (cid:213) j (cid:48)(cid:48) k ( j (cid:48)(cid:48) ) X ( i ) X ( H ) n . The rate of destruction of species i due to radiative pro-cesses is given by S rad = α ( i ) exp (− γ A V ) (7)for photoreactions (where α is a proportionality constantcontaining the dust grain albedo, which we take to be 0.5). http://kida.obs.u-bordeaux1.fr For cosmic ray ionisation, we take S rad = ζ = − s − every-where, unless otherwise specified. Though photoabsorptionof radiation from external sources will affect the chemistry inthe upper layers of a disc, we focus on the bulk of the planetforming material toward the disc mid-plane. We therefore as-sume that all material is well shielded from external sourcesof photons, setting A V = mag everywhere. However, we doinclude photoreactions due to secondary photons generatedby cosmic rays, using the rates from Heays et al. (2017).For the adsorption of species from the gas phase ontothe surfaces of icy grain mantles, we assume S a = X dust S ( i ) (cid:115) kT π m ( i ) (cid:104) π a (cid:105) (8)where X dust is the dust number density, S ( i ) is the stickingcoefficient (taken to be 0.3), m ( i ) is the atomic mass of theadsorbed species, and (cid:104) π a (cid:105) is the average dust grain area.The average grain radius is computed using the maximumgrain size from the dust evolution code, assuming a grainsize distribution n ( a ) d a ∝ a − . d a with the minimum grain-size equal to . µ m. However, we note that our results areonly very weakly dependent on this choice (as shown in sec-tion A).For unsaturated C, N and O species on the surfacesof dust grains, we assume hydrogenation occurs via theEley-Rideal mechanism, at the rate of adsorption of H fromthe gas multiplied by the probability of encounter with thespecies on the grain surface, e.g., gC → gCH → . . . → gCH (where g denotes a species on a grain surface, see Visseret al. 2011).For the thermal desorption of species from dust grainsurfaces into the gas phase, we assume S d ( i ) = . × − σ ν ( i ) exp (cid:18) − E b ( i ) kT (cid:19) , (9)where σ = . × cm − is the surface density of bind-ing sites, E b ( i ) is the binding energy of the i th species onthe surface of the dust grain, and T is the dust temperature(which we assume to be in equilibrium with the gas temper-ature), and ν is the characteristic vibrational frequency ofthe species attached to the grain.Integration of the rate equations in each cell is per-formed independently, using the dlsodes package (Hind-marsh 1983) within the krome package (Grassi et al. 2014)to yield the time-dependent evolution of fractional abun-dance for each species throughout the disc. The coupling to the disc evolution model is achieved in a afirst-order operator-split fashion, which we briefly describehere. The full equation for the surface density of a givenspecies, Σ ( i ) , may be written as ∂ Σ ( i ) ∂ t + R ∂∂ R [ R Σ ( i ) v r ( i )] = Σ ( i ) S ( i ) , (10)where v r ( i ) is the radial velocity of the species (either the gasvelocity or dust velocity, depending on whether the speciesis in the gas or ice phase), and S ( i ) is the source terms due tochemical reactions (Equation 4). We discretize this equation MNRAS000 , 1–14 (2019)
Booth & Ilee
Table 1.
Initial abundances relative to the number of hydrogennuclei, X ( i ) ≡ n ( i )/ n . Based on Eistrup et al. (2016).Species X ( i ) Species X ( i ) H . × − CH . × − H . N . × − He . NH . × − H O . × − H S . × − CO . × − Na . × − CO . × − on a fixed Eulerian grid, as described in Booth et al. (2017).We then proceed by splitting these equations into two steps,first a transport step without the chemical reactions: ∂ Σ ( i ) ∂ t + R ∂∂ R [ R Σ ( i ) v r ( i )] = , (11)where the amount of each species being transported betweenthe different radial cells between times t and t + d t is com-puted as described in Booth et al. (2017). This produces anew estimate for the surface density in each cell, Σ ∗ ( i ) , attime t + dt .The second step in the computation is to integrate ∂ Σ ( i ) ∂ t = Σ ( i ) S ( i ) (12)over a time d t starting from Σ ∗ ( i ) . To do this, we firstcompute mid-plane density of each chemical species via ρ ( i ) = Σ ∗ ( i )/(√ π H ) , where H is the disc-scale height. Wethen pass n ( i ) = ρ ( i )/ m ( i ) , where m ( i ) is the mass of thespecies, into the chemical solver, krome , which integratesthe number density from t to t + d t . The new n ( i ) can thenbe converted back to Σ ( i ) . Finally, the grain size is updatedas described in Booth et al. (2017) and the temperature inthe disc updated. For the initial abundances, we assume the gas is molecu-lar with the abundance of key molecular species given inTable 1, which follow Eistrup et al. (2016), which are com-parable to the abundances in comets (Mumma & Charn-ley 2011). These abundances correspond to C/H, N/H, andO/H ratios of 0.47, 0.89, and 0.85 times solar, respectively(Asplund et al. 2009). The C/O ratio is 0.29, which is be-low the protosolar value of approximately 0.54. These lowerabundances are consistent with the remaining species beinglocked up in refractory solids (Pollack et al. 1994). Whererelevant, we thus assume the remaining C, N, and O atomsare locked up in the cores of dust grains, which for the pur-pose of this study are considered inert.
We first begin with a simple test case in which all trans-port terms are turned off in the disc evolution model, i.e.the chemical reaction network for each cell is integrated in-dependently with the density fixed at the initial value. This allows a direct comparison with recent works exploring theimpact of chemistry on the composition of planets (Eistrupet al. 2016; Cridland et al. 2016; Yu et al. 2016; Cridlandet al. 2017; Eistrup et al. 2018), with which we can bench-mark the smaller chemical network used in this study. Forthis test we present results from the standard disc modelwith a mass of . M (cid:12) and α = − . The evolution in thesemodels is almost entirely driven by ionization due to cosmicrays (see Appendix A or Eistrup et al. 2016). The abun-dances of the dominant carbon and oxygen carriers at 0 and yr are shown in Figure 3, which can be compared withFigures 2 and 3 of Eistrup et al. (2016). The time evolutioncan also be seen in Figure A2.Our model without transport is in good qualitativeagreement with the results from the more extensive networksused by Eistrup et al. (2016, 2018), and Yu et al. (2016).The most significant effect is the depletion of CH inside of10 au, where CH is in the gas phase. The primary path-way for the removal of CH are reactions with C + , which inturn is produced via reactions between CO and He + that isgenerated by cosmic ray ionization. The reactions betweenCH and C + produce C H + and C H + , which can recombinewith electrons to form C H and C H (Walsh et al. 2015;Yu et al. 2016). C H then freezes out, acting as the mainreservoir in our network. In reality, C H is likely furtherhydrogenated on the grain surfaces – Eistrup et al. (2018)find that the main reservoir in this region is the saturatedhydrocarbon C H , which is also frozen out. We thereforeargue that the conversion of CH to larger molecules is ad-equately captured in our model because the time-scale fordepletion (a few 1 yr ) is captured properly. Furthermore,since the binding energies of these two species are similar(Collings et al. 2004; ¨Oberg et al. 2009, see also Penteadoet al. 2017) they will be in the ice phase for similar rangesof parameter space, resulting in similar transport efficiencytoo. We thus argue that the differences are not significant forour goal, i.e. to asses how the carbon and oxygen abundanceof the gas and solids during planet-formation is affected bythe competition between transport and chemistry.The most significant differences come from the fact thatwe do not include detailed grain surface chemistry. However,besides the aforementioned hydrogenation of small hydro-carbons, these differences are most significant on time-scalesof a few Myr. E.g. Eistrup et al. (2018) show that this leadsto the conversion of CO into CO on time-scales of a few Myrin the outer disc where CO is frozen out (see also Bosmanet al. 2018b). Without grain surface chemistry CO can onlybe destroyed in the gas phase, which takes several to
10 Myr .For this reason we focus on the first , where the differ-ences in the CO abundance are within a factor of 2.
We now turn our attention to models involving transport,first considering models with viscous evolution but no radialdrift. In this model, the dust and gas move together as theymove towards the star, similar to the models by Aikawa et al.(1999). The differences in the chemical evolution betweenthe models with and without transport are thus due to thegas at a given location having come from larger radii in thedisc, where the temperatures and densities are lower. Since
MNRAS , 1–14 (2019) hemistry & transport in discs -1 R [au]10 -8 -7 -6 -5 -4 -3 a b un d a n c e no transport 10 R [au] ® = 10 ¡ R [au] ® = 10 ¡ COCO CH H O C H Figure 3.
Abundances relative to hydrogen, X ( i ) , of the major carbon and oxygen carriers after yr in models without radial drift.In the left panel transport is neglected entirely, while the other two panels show models in which the dust moves at the same speed asthe gas. Solid lines denote gas phase species and dashed lines denote ice phase species and the dotted lines show the initial abundance ofeach species. The vertical bars show the approximate ice line locations. Note that C H should be considered a proxy for hydrocarbonsmore generally. -1 R [au]10 -8 -7 -6 -5 -4 -3 a b un d a n c e no chem 10 R [au] ® = 10 ¡ R [au] ® = 10 ¡ COCO CH H O C H Figure 4.
Same as Figure 3, but for models including radial drift. The left hand panel shows the α = − case in including onlyadsorption and desorption. each parcel of gas now spends a limited amount of time inthe regions where a given chemical process is happening (e.g.the depletion of CH ), this naturally leads to competitionbetween the chemical time-scales and transport time-scales(which are given in subsection 2.2).Transport thus reduces the impact of chemical reactionson the composition in the inner regions. Even in our low ac-cretion rate model with α = − , viscous transport carriesgas from 2 au onto the star more rapidly than CO depletes( (cid:38) yr ), as discussed in detail by Bosman et al. (2018a).However, at α = − , viscous transport is not able to over-come the depletion of CH , due to a combination of a shorterchemical time-scale (a few yr ) and the longer time thatmaterial spends at 10 au than at 2 au.At α = − , viscous transport becomes efficient enoughthat gas-phase chemistry no longer has a significant effect on the abundance of the major carbon and oxygen carriers in-side 10 au. While there is still some conversion of CH tolarger hydrocarbons, this only reduces the abundance of CH by a factor of ∼ . Thus the composition of the mid-planeinside the CO ice line ( ∼
30 au ) is largely set by the compo-sition of the material being carried in from larger radii.
The primary effect of radial drift is to enhance the transportof molecular ices from the outer disc (Booth et al. 2017). Themolecular species frozen out enter the gas phase when thegrains pass each species’ respective ice lines, enhancing thelocal volatile abundance inside the ice line. The strength ofthis enhancement is controlled by the relative time-scaleson which the gas and dust evolve: rapid radial drift leads
MNRAS000
MNRAS000 , 1–14 (2019)
Booth & Ilee to a rapid increase of volatile abundance inside the ice line.In the opposite limit, when radial drift is slow, there is noenhancement of the abundances because the gas and dustmove together. Thus the volatile species are transportedaway from the ice lines in the gas phase at the same ratethey are brought there in ices on the grains.The influence of radial drift can be most clearly seen inthe left hand panel of Figure 4, which shows the α = − model including radial drift and adsorption/desorption, butneglecting other chemical processes. The rapid dust trans-port in this model means that 90 per cent of the volatilesare delivered into the inner disc within 0.5 Myr. Combinedwith the low gas accretion rate, this results in high gas phasemolecular abundances in the inner disc, along with a reduceddust-to-gas ratio (Figure 2).For CO, which enters the gas phase at around 40 auwhere viscous transport is slow, this results in a spike in thegas phase CO abundance. Inside 5 au the CO abundance isdominated by gas that originated inside the ice line and hasbeen carried inwards by the gas; thus, the CO abundancereturns to its initial value. For the species that have theirice lines closer to the star, such as CH , CO and H O, theabundance inside the ice lines has become close to constantbecause the transport time-scale is less than . Out-side their ice lines, the gas phase abundance of all molecularspecies remains low because any molecules that diffuse be-yond the ice line freeze out and are carried back on the grainsurfaces. This results in the ice-to-dust abundance of indi-vidual species being enhanced by as much as a factor of ∼ near their ice lines.The interplay between gas-phase chemistry, radial driftand viscous evolution can be seen in the middle panel of Fig-ure 4, for the model with α = − . Here we see again thatthe conversion of CH to larger hydrocarbons is efficient in-side of 10 au, removing the spike in CH at the ice line.However, radial drift is efficient enough that the the CH abundance is still high after 1 Myr. The high CH abundanceat ∼
10 au also translates into a higher rate of hydrocarbonformation, which can be seen through the higher abundanceof C H ice relative to water ice in the middle panel of Fig-ure 4 compared to Figure 3. However, efficient radial driftprevents the newly formed ice reaching high abundances,with the ice instead being transported in to its own ice line,as suggested in Booth et al. (2017). A similar process is re-sponsible for the higher CO abundance in the inner region.In this case it is the formation of H CO from CO near theice line that is responsible. Once formed, the H CO freezesout and is transported on the grains to around 5 au, whereit enters the gas phase and is converted back to CO. Theenhancement of CO in the inner regions thus relies on theformation of H CO. We note that when grain surface reac-tions are included Eistrup et al. (2016) find CO is formedinstead. While CO would freeze out also and be transportedinwards, the conversion of CO into CO is slow compared tothe transport. Thus the effect of CO conversion to CO andtransport on the grain surfaces would be to further enhancethe CO abundance in the inner regions, rather than CO.The evolution of species such as CO, CH , CO andH O after would be characterized with an inside-outdepletion of the species in order of their snow line locations,because species closest to the star accreted onto the star first (Booth et al. 2017). Including chemical reactions onlyslightly modifies this picture, with the depletion of CO andCH being slightly accelerated with respect to the resultswithout chemical reactions through their processing to CO and larger hydrocarbon ices, respectively.In the α = − model the evolution is similar bothwith and without radial drift. This is due to the viscoustime-scale ( t ν ) and radial drift time-scale ( t drift ) being com-parable, so radial drift only modestly affects the transportof dust. Furthermore, both of the transport time-scales areshorter than the chemical time-scales for the depletion of themajor carbon and oxygen carriers. At α = − the impactof radial drift is to increase the abundances by a factor ∼ due to the radial drift of ices. The enhancement in molecu-lar abundances is only slightly higher than the enhancementin dust-to-gas ratio because the dust drifts inwards at onlyslightly higher speeds than the gas. However, suppressingthe influence of radial drift entirely would require t ν (cid:28) t drift .Similar to the α = − models without radial drift, we againsee that gas-phase reactions only have a modest effect on theabundances, although the abundance of hydrocarbons pro-duced from CH still reaches 20 per cent of CH abundance.In Figure 5 we show how transport influences the evo-lution of the main nitrogen reservoirs, N and NH . In theouter disc we find that the evolution is slow(time-scales (cid:29) ), even in the absence of transport. We find thatN converts to NH , as found by Schwarz & Bergin (2014)in models with a high initial N abundance. Conversely,Eistrup et al. (2016) report the conversion of NH to N in the outer disc, again on time-scales much longer than . The difference is due to the destruction of NH iceby cosmic-ray induced photons, which is neglected in thiswork. Nevertheless, the evolution of the nitrogen abundancein the outer disc is expected to be slow.In the inner disc NH is depleted and the N abundanceincreases. The remaining NH abundance inside the NH iceline is sensitive to the rate of NH production from N . Wefind that this is fastest via the pathway N γ CR → N + H → NH + H → NH + → NH + O → NH + − → NH . (13)Outside the water ice line where water is not present inthe gas phase the NH + → NH + step occurs more slowlyvia reactions with H , resulting in lower abundances. In thecolder regions further out NH formation can also occur ongrain surfaces.As in the case of carbon and oxygen bearing species, wesee that even low levels of transport prevent the depletion ofNH , largely because the depletion time-scale, ∼ . M y r , islong compared to the transport time-scales in the inner disc.We note that high NH abundances inside the NH ice lineis in conflict with recent measurements of the NH abun-dance in the inner regions of discs ( < − ; Pontoppidanet al. 2019), even for α = − . However, it may be possiblethat the NH abundance observed in the upper layers is notrepresentative of the bulk abundance on the disc. A simi-lar conflict between observed CO abundances and modelsincluding transport was reported by Bosman et al. (2018a).Neither the conversion of N to NH nor the reverse pro-cess change the total gas phase nitrogen abundance signifi-cantly – in the inner disc both species are in the gas phase, MNRAS , 1–14 (2019) hemistry & transport in discs -1 R [au]10 -8 -7 -6 -5 -4 -3 a b un d a n c e no transport 10 R [au] ® = 10 ¡ R [au] ® = 10 ¡ N NH -1 R [au]10 -8 -7 -6 -5 -4 -3 a b un d a n c e no chem 10 R [au] ® = 10 ¡ R [au] ® = 10 ¡ N NH Figure 5.
Abundances relative to hydrogen of the key nitrogen bearing molecules. The models are the same as Figure 3 and Figure 4.Solid lines denote gas phase species, dashed lines denote ice phase species. Top: models with viscous transport only. Bottom: modelsincluding radial drift. while in the outer disc the destruction of N is slow. How-ever, due to the long viscous time in the outer disc (1 Myrat 50 au for α = − ), transport only has a modest effect onthe evolution of nitrogen in the outer disc. Conversely, in theinner disc where the viscous time is much shorter, transportreduces the conversion of NH to N .Models with radial drift also have a similar effect onthe nitrogen reservoirs as the carbon and oxygen reservoirs.Radial drift causes enhancements in the gas phase N abun-dance inside the N ice line by a factor 2 – 3. The NH abundance increases by factors of 3–10 due to radial driftand is further enhanced by the formation of extra NH inthe outer disc, which freezes out and is carried in by thedrifting grains.We briefly consider how sensitive the results presentedhere are to assumptions about the disc model. Notably, theradial drift and viscous time-scales are not very sensitive tomodel assumptions, apart from the dependence on α . Forthis reason, we investigated how the chemical time-scalesvary in models without transport, varying the disc mass,cosmic-ray ionization rate and average grain size used in thechemical models. The results are presented in detail in Ap- pendix A, which shows that the time-scale for the depletionof CH and CO are only sensitive to the cosmic-ray ioniza-tion rate. Since we adopt the canonical value of cosmic-rayionization rate ( − s − ), which assumes no shielding ofthe interstellar cosmic ray flux via stellar winds or magneticfields (Cleeves et al. 2013a), our results likely represent alower limit of the importance of transport relative to mid-plane chemical kinetics. The bulk abundances of the main atomic species, e.g. car-bon, oxygen, and nitrogen, are of particular interest in thecontext of planet formation since the composition of exo-planet atmospheres must be connected to the competition ofgas and solids in discs (e.g. ¨Oberg et al. 2011; Madhusudhanet al. 2014). We first explore the evolution of the C/O ratioin models without radial drift (Figure 6). In these models,the total abundance (gas and ice phase) of carbon and oxy-gen bearing species does not change, only the partitioningof these species between the gas and ice phases.The changes in composition driven by chemical pro-
MNRAS000
MNRAS000 , 1–14 (2019) Booth & Ilee -1 R [au] [ C = O ] solarno transport 10 -1 R [au] solar ® = 10 ¡ -1 R [au] solar ® = 10 ¡ Figure 6.
Evolution of the C/O ratio for models with viscous transport, but without radial drift or diffusion. For comparison the lefthand panel shows the evolution without transport. The solid lines denote the gas phase abundances, while the dotted lines show the iceabundances. -1 R [au] -1 R e l a t i v e c a r b o n a b un d a n c e no chem 10 -1 R [au] -1 ® = 10 ¡ -1 R [au] -1 ® = 10 ¡ -1 R [au] [ C = O ] solarno chem 10 -1 R [au] ® = 10 ¡ -1 R [au] ® = 10 ¡ Figure 7.
C/O ratio and carbon abundance evolution for models with radial drift. The carbon abundance shown is relative to the initialtotal carbon abundance ( . × − , . × solar ). Note the different scales for the C/O ratio between the left and right panels.MNRAS , 1–14 (2019) hemistry & transport in discs cesses are best seen in the model without transport (and α = − ). Inside the water ice line (0.5 au), all of the majorcarbon and oxygen are in the gas phase and thus chemicalprocessing does not change the gas phase carbon or oxygenabundance. Between the water and CO ice lines (1.5 au) thegas phase carbon and oxygen abundance is dominated by COand CO , with a small contribution from hydrocarbons (CH and C H – a proxy for C H ; see subsection 3.1). Here, theC/O ratio increases as CO is destroyed, producing CO andO, with O eventually forming H O, which freezes out. Be-tween the CO and CH ice lines (10 au), CO and CH areinitially the main gas phase carbon and oxygen carriers. CH gets converted to larger hydrocarbons (C H ), which freezeout causing the C/O ratio to reduce to 1. The spike in C/Oratio near the CO ice line is due to the slightly lower bind-ing energy of C H than CO (note that C H also has alower binding energy than CO , so Eistrup et al. 2018 see asimilar spike in C/O). Outside of the CH ice line CO is al-ways the dominant the gas phase carbon and oxygen carrier,so the C/O ratio does not evolve. The nitrogen abundanceevolution can be derived from Figure 5 and is is slow, thusthe C/N of N/O ratio evolution will be driven primarly bythe evolution of carbon and oxygen, respectively.As the speed of the viscous transport increases, thechanges in composition driven by chemical reactions becomerestricted to regions of the disc at larger and larger radii. Al-ready at α = − the effect of transport is first noticeablebetween the H O and CO ice lines, since CO no longer de-pletes (see also Bosman et al. 2018a), while at α = − thecarbon and oxygen abundance hardly evolve. At α = − the inward movement of the ice lines as the disc cools gen-erates the only significant changes in abundance ratio.The influence of radial drift on the bulk abundancesis sensitive to whether radial drift or viscous evolution isfaster. When radial drift is included we see spikes in thegas phase C/H ratio at α = − (Figure 7), as reportedby Booth et al. (2017). These spikes are driven by volatilespecies being left at their snow lines when grains drift pastthem. We also see a peak in the C/O ratio inside of 10 au,which is driven by CH entering the gas phase. Since radialdrift depletes the disc of dust (and ice phase species) withina few 0.1 Myr the C/O ratio at 10 au begins to decline. Thisdecline is due to a combination of the conversion of CH toC H and the viscous transport of CH inwards from the iceline. After ∼ . , conversion of CH to C H generatesa second spike in C/O ratio at a few au (where C H entersthe gas phase). At α = − the viscous and radial drifttimes are comparable, reducing the influence of radial drift.Nevertheless, we still find an enhancement in the gas phaseC/H ratio by a factor 2 – 3, since radial drift is still able tobring extra material into the inner regions. We again findC/O > between the CH and CO snow lines since CH is transported into this region faster than it is converted toC H . In this paper we have investigated how the abundances ofthe main carbon, oxygen, and nitrogen carriers in the mid- plane of protoplanetary discs evolve under the competitionbetween chemical kinetics and radial transport. The effectsof transport differ between gas dominated transport (viscousevolution, and the transport of small grains at the viscousspeed) and radial drift dominated transport. The main effectof gas dominated transport is to reduce the amount of timethat a parcel of gas spends in regions where the processingoccurs. Thus gas transport weakens the depletion of CO ,CH and CO from the gas phase that otherwise occurs insidetheir ice lines. Radial drift acts differently, enhancing thetransport of molecular species inwards on grain surfaces. Themolecules are then left at in the gas phase as they sublimateat their ice lines. Thus the main effect of radial drift is toenhance the abundances of molecular species inside their icelines. Their subsequent evolution is controlled by transportfurther inward in the gas phase, along with processing bychemical reactions. A secondary effect of radial drift is thatnew species produced by chemical processing (e.g. C H orH CO) can freeze out and be carried further in before re-entering the gas phase.The results presented here are relatively insensitive tothe parameters in the model, except those that govern therate of transport. In Appendix A we explore how the chemi-cal depletion time-scales depend on model parameters in theabsence of transport, showing that the results are only sen-sitive to the cosmic-ray ionization rate, ζ , which facilitatesthe removal of CH and CO from the gas phase. If ζ is lowerthan the canonical value of − (as suggested by Cleeveset al. 2013a), then transport is expected to dominate for α (cid:38) − , which is likely the case for most discs. Our resultsare otherwise insensitive to the disc mass or assumptionsabout the grain size distribution. This means that α and ζ essentially control the relative importance of transport andchemical reactions. There will also be a weak dependence onthe temperature of the disc as the ice lines move to largerradius where the viscous time is longer.The models presented here support the ideas put for-ward in Booth et al. (2017), using similar assumptionsfor the transport of chemical species but simpler assump-tions about their abundances. I.e. in models with efficientgrain growth and radial drift the abundances inside ice linesare mostly controlled by the inward transport of molecu-lar species from larger distances on grain surfaces, togetherwith transport in the gas phase inside the ice line. Reduc-ing the grain sizes weakens radial drift and results in thecompositions being set largely by the composition at largerradii. The exception to this is when transport times are long( (cid:29) at
10 au ) and ionization rates are high ( ζ (cid:38) − ),in which case chemical processing in the mid-plane may beable to remove CH (and eventually CO) from the gas phase.We also verify the suggestion that, by processing moleculesto less volatile species (such as CH to C H ), chemical re-actions act to aid the transport of volatile species to smallerradii (Figure 7).The models presented here assume a turbulent viscos-ity, but it is possible that accretion in protoplanetary discs isprimarily laminar, as discussed in subsection 2.2. However, alow level of radial diffusion associated with weak turbulencewould not affect the results presented here because the ef-fects of radial diffusion are already localised to regions closeto the snow lines. The primary effect of weak turbulence on MNRAS , 1–14 (2019) Booth & Ilee the models is thus on grain growth and radial drift, whichis limited by turbulent fragmentation in our α = . model.Lowering the strength of turbulence in that model by a fac-tor 10 would result in efficient grain growth and radial drift,as in the α = − model. Reducing the strength of turbu-lence further would have little effect because grain growthbecomes limited by radial drift instead.A potentially more important effect associated withlaminar accretion in discs is the possibility that the accretionoccurs in a narrow layer far away from the mid-plane (e.g.Gressel et al. 2015). In this case, vertical mixing is neededto couple the effects of accretion to the composition of themid-plane. Two-dimensional models are needed to explorethis in detail, but it is unlikely that the accretion flow canbe completely decoupled from the mid-plane because verti-cal mixing should occur a factor of ( R / H ) ∼ – timesfaster than the radial mixing for the same α . Thus an ex-traordinarily low α would be required to prevent verticalmixing altogether.Dust traps, if widely present in protoplanetary discs,could have a significant impact on the abundance evolution.Such traps have typically been invoked to explain why dustgrains can survive in discs for several Myr (e.g. Pinilla et al.2012). The impact of traps on the composition will dependon the evolution of the traps themselves. If the traps were tomove with the gas, then the net effect would be similar to ourmodels neglecting radial drift. However, if the traps remainstationary, preventing the inflow of dust entirely, the chem-ical evolution could be different as volatiles remain trappedat their original location. This would lead to the depletionof volatile species from the gas phase as the gas in the innerdisc is accreted onto the star. This effect has been invoked toexplain the depletion of refractory material accreting ontoHerbig stars with transition discs (Kama et al. 2015).Although the high CO abundance in our models ap-pears in conflict with the observed depletion of CO in TWHydra (Favre et al. 2013; Zhang et al. 2017), this can poten-tially be reconciled due to TW Hydra’s age (8 Myr) and lowaccretion rate (Calvet et al. 2002; Donaldson et al. 2016),as long as the ionization rate is not too low ( ζ ∼ − ).In this case the viscous time-scale near the CO ice line is ∼ , long enough that grain surface processes could de-plete the CO abundance (requiring ∼ , Bosman et al.2018b). The low CO line fluxes observed for a range of discs(e.g. those in Lupus; Miotello et al. 2017) point to a moregeneral depletion of CO within a few Myr. However, for sys-tems requiring a higher α (e.g. HD 163296; Boneberg et al.2016) or low ionization rate (Cleeves et al. 2013a), grain sur-face processes are unlikely to be able to deplete CO rapidlyenough. The large difference in the gas-phase abundances in mod-els with and without efficient radial drift (e.g. models with α = − and − , respectively) suggests that exoplanetcomposition will be a powerful way to constrain planet for-mation by pebble accretion. Most models of giant planetformation by via pebble accretion assume large pebble sizesthat drift efficiently (Lambrechts & Johansen 2014; Mor- bidelli et al. 2015; Bitsch et al. 2015, 2019; Ida et al. 2016).Planets forming in these models will have super-solar abun-dances unless they accrete their envelopes beyond the COice line (Booth et al. 2017). Our results suggest that theseplanets will likely also have C / O > and possibly C / O (cid:29) if they accrete their envelopes inside the CH ice line (whichis the case in the models of Bitsch et al. 2019, for exam-ple). Recently, Johansen et al. (2018) suggested that giantplanet formation via pebble accretion can occur for moremodest pebble sizes, as long as the pebbles have settled tothe mid-plane. Such models will be associated with moremodest enhancements in the abundances, more similar toour α = − model (including radial drift). We have explored the competition between chemical reac-tions in the mid-plane and transport due to viscous evolu-tion and radial drift on the composition of protoplanetarydiscs. Our models show that the radial transport associatedwith accretion is able to transport material faster than itcan be depleted due to chemical processing if the viscous α parameter is high ( α ∼ − , Figure 3) or the ionization rateis low ( ζ ∼ − , Figure A2). For α = − and ζ = − ,transport will still overcome the depletion of CO from thegas phase inside its ice line, but CO and CH will remaindepleted.Including radial drift enhances the gas phase abun-dances of the dominant molecular species inside their icelines (Figure 4). The strength of enhancement depends onhow much faster radial drift is than transport in the gas.For commonly assumed parameters in models of dust evolu-tion and planet formation pebble accretion (Birnstiel et al.2012; Lambrechts & Johansen 2014; Bitsch et al. 2015, e.g.),this leads to factor ∼ enhancements in the abundances.In models with more modest radial drift (grain sizes in therange 0.1–1 mm), radial drift can still lead to a factor 2–3enhancement in the abundances inside 10 au.In the coming years, the joint operation of the JamesWebb Space Telescope (JWST) and the Atacama LargeMillimetre/submillimetre Array (ALMA) will enable robustmeasurements of the elemental composition of hot Jupiteratmospheres (such as C/O) alongside measurements of thechemical composition of planet forming discs. Our resultsdemonstrate the need to consider the interplay between bothphysical and chemical processes during protoplanetary discevolution when interpretting such observations, with impli-cations for the connection of planet composition to the com-position of the disc from which they formed. ACKNOWLEDGEMENTS
We would like to thank Catherine Walsh, Cathie Clarke,and Mihkel Kama for useful discussions and careful read-ing of the manuscript. We additionally thank the referee,Alexander Cridland, for helpful suggestions which have im-proved this work. RAB and JDI gratefully acknowledge sup-port from the DISCSIM project, grant agreement 341137under ERC-2013-ADG. JDI also acknowledges support fromthe STFC under ST/R000549/1. This project has made use
MNRAS , 1–14 (2019) hemistry & transport in discs of the SciPy stack (Jones et al. 2001), including NumPy(Oliphant 2007) and Matplotlib (Hunter 2007). REFERENCES
Aikawa Y., Umebayashi T., Nakano T., Miyama S. M., 1999, ApJ,519, 705Ali-Dib M., 2017, MNRAS, 464, 4282Andrews S. M., Williams J. P., 2007, ApJ, 659, 705Ansdell M., et al., 2016, ApJ, 828, 46Asplund M., Grevesse N., Sauval A. J., Scott P., 2009, AnnualReview of Astronomy and Astrophysics, 47, 481Bai X.-N., Stone J. M., 2013, ApJ, 769, 76Bell K. R., Lin D. N. C., 1994, ApJ, 427, 987Bergin E. A., et al., 2013, Nature, 493, 644Bergin E. A., Du F., Cleeves L. I., Blake G. A., Schwarz K., VisserR., Zhang K., 2016, ApJ, 831, 101Birnstiel T., Klahr H., Ercolano B., 2012, A&A, 539, A148Bitsch B., Lambrechts M., Johansen A., 2015, A&A, 582, A112Bitsch B., Izidoro A., Johansen A., Raymond S. N., MorbidelliA., Lambrechts M., Jacobson S. A., 2019, A&A, 623, A88Boneberg D. M., Pani´c O., Haworth T. J., Clarke C. J., Min M.,2016, MNRAS, 461, 385Booth R. A., Clarke C. J., Madhusudhan N., Ilee J. D., 2017,MNRAS, 469, 3994Bosman A. D., Tielens A. G. G. M., van Dishoeck E. F., 2018a,A&A, 611, A80Bosman A. D., Walsh C., van Dishoeck E. F., 2018b, A&A, 618,A182Calvet N., D’Alessio P., Hartmann L., Wilner D., Walsh A., SitkoM., 2002, ApJ, 568, 1008Clarke C. J., et al., 2018, ApJ, 866, L6Cleeves L. I., Adams F. C., Bergin E. A., 2013a, ApJ, 772, 5Cleeves L. I., Adams F. C., Bergin E. A., Visser R., 2013b, ApJ,777, 28Collings M. P., Anderson M. A., Chen R., Dever J. W., Viti S.,Williams D. A., McCoustra M. R. S., 2004, MNRAS, 354,1133Cridland A. J., Pudritz R. E., Alessi M., 2016, MNRAS, 461, 3274Cridland A. J., Pudritz R. E., Birnstiel T., Cleeves L. I., BerginE. A., 2017, MNRAS, 469, 3910Donaldson J. K., Weinberger A. J., Gagn´e J., Faherty J. K., BossA. P., Keiser S. A., 2016, ApJ, 833, 95Eistrup C., Walsh C., van Dishoeck E. F., 2016, A&A, 595, A83Eistrup C., Walsh C., van Dishoeck E. F., 2018, A&A, 613, A14Ercolano B., Rosotti G. P., Picogna G., Testi L., 2017, MNRAS,464, L95Evans M. G., Ilee J. D., Boley A. C., Caselli P., Durisen R. H.,Hartquist T. W., Rawlings J. M. C., 2015, MNRAS, 453, 1147Favre C., Cleeves L. I., Bergin E. A., Qi C., Blake G. A., 2013,ApJ, 776, L38Flaherty K. M., et al., 2017, ApJ, 843, 150Flaherty K. M., Hughes A. M., Teague R., Simon J. B., AndrewsS. M., Wilner D. J., 2018, ApJ, 856, 117Fromang S., Latter H., Lesur G., Ogilvie G. I., 2013, A&A, 552,A71Gail H.-P., Trieloff M., 2017, A&A, 606, A16Grassi T., Bovino S., Schleicher D. R. G., Prieto J., Seifried D.,Simoncini E., Gianturco F. A., 2014, MNRAS, 439, 2386Greaves J. S., Rice W. K. M., 2010, MNRAS, 407, 1981Gressel O., Turner N. J., Nelson R. P., McNally C. P., 2015, ApJ,801, 84Guilloteau S., Dutrey A., Pi´etu V., Boehler Y., 2011, A&A, 529,A105Heays A. N., Bosman A. D., van Dishoeck E. F., 2017, A&A, 602,A105 Heinzeller D., Nomura H., Walsh C., Millar T. J., 2011, ApJ, 731,115Henning T., Semenov D., 2013, Chemical Reviews, 113, 9016Hindmarsh A. C., 1983, IMACS Transactions on Scientific Com-putation, 1, 55Hunter J. D., 2007, Computing in Science Engineering, 9, 90Ida S., Guillot T., Morbidelli A., 2016, A&A, 591, A72Ilee J. D., Boley A. C., Caselli P., Durisen R. H., Hartquist T. W.,Rawlings J. M. C., 2011, MNRAS, 417, 2950Ilee J. D., et al., 2017, MNRAS, 472, 189Johansen A., Ida S., Brasser R., 2018, arXiv e-prints, p.arXiv:1811.00523Jones E., Oliphant T., Peterson P., et al., 2001, SciPy: Opensource scientific tools for Python,
Kama M., Folsom C. P., Pinilla P., 2015, A&A, 582, L10Kama M., et al., 2016, A&A, 592, A83Klarmann L., Ormel C. W., Dominik C., 2018, A&A, 618, L1Krijt S., Schwarz K. R., Bergin E. A., Ciesla F. J., 2018, ApJ,864, 78Lambrechts M., Johansen A., 2012, A&A, 544, A32Lambrechts M., Johansen A., 2014, A&A, 572, A107Lodato G., Scardoni C. E., Manara C. F., Testi L., 2017, MNRAS,472, 4700Lodders K., 2010, Exoplanet Chemistry. Wiley, p. 157,doi:10.1002/9783527629763.ch8Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603Madhusudhan N., 2012, ApJ, 758, 36Madhusudhan N., Amin M. A., Kennedy G. M., 2014, ApJ, 794,L12Madhusudhan N., Bitsch B., Johansen A., Eriksson L., 2017, MN-RAS, 469, 4102Manara C. F., et al., 2016, A&A, 591, L3Manara C. F., Morbidelli A., Guillot T., 2018, A&A, 618, L3Millar T. J., Farquhar P. R. A., Willacy K., 1997, A&AS, 121,139Miotello A., et al., 2017, A&A, 599, A113Morbidelli A., Lambrechts M., Jacobson S., Bitsch B., 2015,Icarus, 258, 418Mordasini C., van Boekel R., Molli`ere P., Henning T., BennekeB., 2016, ApJ, 832, 41Mulders G. D., Pascucci I., Manara C. F., Testi L., Herczeg G. J.,Henning T., Mohanty S., Lodato G., 2017, ApJ, 847, 31Mumma M. J., Charnley S. B., 2011, Annual Review of Astron-omy and Astrophysics, 49, 471Najita J. R., Kenyon S. J., 2014, MNRAS, 445, 3315Nomura H., Aikawa Y., Nakagawa Y., Millar T. J., 2009, A&A,495, 183¨Oberg K. I., Garrod R. T., van Dishoeck E. F., Linnartz H., 2009,A&A, 504, 891¨Oberg K. I., Murray-Clay R., Bergin E. A., 2011, ApJ, 743, L16Oliphant T. E., 2007, Computing in Science Engineering, 9, 10Ormel C. W., Klahr H. H., 2010, A&A, 520, A43Penteado E. M., Walsh C., Cuppen H. M., 2017, ApJ, 844, 71Pinilla P., Birnstiel T., Ricci L., Dullemond C. P., Uribe A. L.,Testi L., Natta A., 2012, A&A, 538, A114Pinte C., Dent W. R. F., M´enard F., Hales A., Hill T., Cortes P.,de Gregorio-Monsalvo I., 2016, ApJ, 816, 25Piso A.-M. A., ¨Oberg K. I., Birnstiel T., Murray-Clay R. A., 2015,ApJ, 815, 109Piso A.-M. A., Pegues J., ¨Oberg K. I., 2016, ApJ, 833, 203Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., RoushT., Fong W., 1994, ApJ, 421, 615Pontoppidan K. M., Salyk C., Banzatti A., Blake G. A., WalshC., Lacy J. H., Richter M. J., 2019, arXiv e-prints, p.arXiv:1902.03647Rafikov R. R., 2017, ApJ, 837, 163Schwarz K. R., Bergin E. A., 2014, ApJ, 797, 113MNRAS000
Kama M., Folsom C. P., Pinilla P., 2015, A&A, 582, L10Kama M., et al., 2016, A&A, 592, A83Klarmann L., Ormel C. W., Dominik C., 2018, A&A, 618, L1Krijt S., Schwarz K. R., Bergin E. A., Ciesla F. J., 2018, ApJ,864, 78Lambrechts M., Johansen A., 2012, A&A, 544, A32Lambrechts M., Johansen A., 2014, A&A, 572, A107Lodato G., Scardoni C. E., Manara C. F., Testi L., 2017, MNRAS,472, 4700Lodders K., 2010, Exoplanet Chemistry. Wiley, p. 157,doi:10.1002/9783527629763.ch8Lynden-Bell D., Pringle J. E., 1974, MNRAS, 168, 603Madhusudhan N., 2012, ApJ, 758, 36Madhusudhan N., Amin M. A., Kennedy G. M., 2014, ApJ, 794,L12Madhusudhan N., Bitsch B., Johansen A., Eriksson L., 2017, MN-RAS, 469, 4102Manara C. F., et al., 2016, A&A, 591, L3Manara C. F., Morbidelli A., Guillot T., 2018, A&A, 618, L3Millar T. J., Farquhar P. R. A., Willacy K., 1997, A&AS, 121,139Miotello A., et al., 2017, A&A, 599, A113Morbidelli A., Lambrechts M., Jacobson S., Bitsch B., 2015,Icarus, 258, 418Mordasini C., van Boekel R., Molli`ere P., Henning T., BennekeB., 2016, ApJ, 832, 41Mulders G. D., Pascucci I., Manara C. F., Testi L., Herczeg G. J.,Henning T., Mohanty S., Lodato G., 2017, ApJ, 847, 31Mumma M. J., Charnley S. B., 2011, Annual Review of Astron-omy and Astrophysics, 49, 471Najita J. R., Kenyon S. J., 2014, MNRAS, 445, 3315Nomura H., Aikawa Y., Nakagawa Y., Millar T. J., 2009, A&A,495, 183¨Oberg K. I., Garrod R. T., van Dishoeck E. F., Linnartz H., 2009,A&A, 504, 891¨Oberg K. I., Murray-Clay R., Bergin E. A., 2011, ApJ, 743, L16Oliphant T. E., 2007, Computing in Science Engineering, 9, 10Ormel C. W., Klahr H. H., 2010, A&A, 520, A43Penteado E. M., Walsh C., Cuppen H. M., 2017, ApJ, 844, 71Pinilla P., Birnstiel T., Ricci L., Dullemond C. P., Uribe A. L.,Testi L., Natta A., 2012, A&A, 538, A114Pinte C., Dent W. R. F., M´enard F., Hales A., Hill T., Cortes P.,de Gregorio-Monsalvo I., 2016, ApJ, 816, 25Piso A.-M. A., ¨Oberg K. I., Birnstiel T., Murray-Clay R. A., 2015,ApJ, 815, 109Piso A.-M. A., Pegues J., ¨Oberg K. I., 2016, ApJ, 833, 203Pollack J. B., Hollenbach D., Beckwith S., Simonelli D. P., RoushT., Fong W., 1994, ApJ, 421, 615Pontoppidan K. M., Salyk C., Banzatti A., Blake G. A., WalshC., Lacy J. H., Richter M. J., 2019, arXiv e-prints, p.arXiv:1902.03647Rafikov R. R., 2017, ApJ, 837, 163Schwarz K. R., Bergin E. A., 2014, ApJ, 797, 113MNRAS000 , 1–14 (2019) Booth & Ilee
Semenov D., Wiebe D., 2011, The Astrophysical Journal Supple-ment Series, 196, 25Stammler S. M., Birnstiel T., Pani´c O., Dullemond C. P., DominikC., 2017, A&A, 600, A140Suzuki T. K., Inutsuka S.-i., 2009, ApJ, 691, L49Takeuchi T., Lin D. N. C., 2002, ApJ, 581, 1344Tazzari M., et al., 2016, A&A, 588, A53Teague R., et al., 2018, ApJ, 864, 133Tscharnuter W. M., Gail H. P., 2007, A&A, 463, 369Visser R., Doty S. D., van Dishoeck E. F., 2011, A&A, 534, A132Walsh C., Herbst E., Nomura H., Millar T. J., Weaver S. W.,2014, Faraday Discussions, 168, 389Walsh C., Nomura H., van Dishoeck E., 2015, A&A, 582, A88Weidenschilling S. J., 1977, MNRAS, 180, 57Youdin A. N., Lithwick Y., 2007, Icarus, 192, 588Yu M., Willacy K., Dodson-Robinson S. E., Turner N. J., EvansNeal J. I., 2016, ApJ, 822, 53Zhang K., Bergin E. A., Blake G. A., Cleeves L. I., Schwarz K. R.,2017, Nature Astronomy, 1, 0130Zhu Z., Hartmann L., Nelson R. P., Gammie C. F., 2012, ApJ,746, 110
APPENDIX A: DEPENDENCE ON MODELPARAMETERS
Here we investigate how the conclusions of our study mightdepend on the assumed parameters of the model. In partic-ular we vary the cosmic-ray ionization rate, the disc mass,and the average grain size used in the chemical models.For brevity, we report only the results of models neglect-ing transport.The cosmic-ray ionization rate, ζ , is the most impor-tant parameter. Figure A1 shows that the depletion of CH and CO from the gas phase are sensitive to this parameter.Reducing the cosmic-ray ionization rate to ζ = − s − issufficient to significantly reduce the depletion of CH . Nev-ertheless, CH is still depleted by a factor ∼ within aMyr, which is fast enough to compete with viscous trans-port at α = − , but not α = − . Thus our suggestionthat mid-plane chemical kinetics will have a significant im-pact at low accretions remains true. Even lower ionizationrates would slow the depletion of CH further. However,short-lived radioactive nuclei are expected to prevent ion-ization rates falling much below − s − (e.g. Cleeves et al.2013b).Higher cosmic-ray ionization rates would lead to qual-itatively different conclusions. At α = − and below, thedepletion of CO from the gas phase would now become sig-nificant within 1 Myr. Furthermore, the time-scale for thedepletion of CH reduces to ∼ . (Figure A2), whichis now comparable with the viscous transport time-scale at α = − . Under such conditions transport is unlikely tostrongly suppress the depletion of CO and CH . However,we note that such high ionization rates are not typicallyexpected for T Tauri star discs, instead the cosmic-ray ion-ization rates are mpre likely lower than the canonical valueof − s − (Cleeves et al. 2013a).In addition to varying the ionization rate, we conductedmodels in which the disc mass has been varied. We find thatthe surface density in the disc does not greatly affect thecomposition, as can be seen by comparing the models with M d = . M (cid:12) and . M (cid:12) . At M d = . M (cid:12) we see differ-ences, but these are due to the higher temperature in the disc, rather than higher density. Since these differences arelimited to within a few au from the star where the transporttime-scales are short, we conclude that disc mass does notmuch affect our conclusions.Finally, we tested the effects of our assumptions aboutthe grain size distribution. In the main text we compute theaverage grain area assuming a grain size distribution with n ( a ) d a ∝ a − . d a between 0.1 µ m and a max , where a max is setby the dust evolution model. To test this assumption wecompare our results with a model where the grain size usedin the gas-grain reactions is fixed at . µ m. Figure A1 showsthat the evolution is largely insensitive to the average grainarea. This paper has been typeset from a TEX/L A TEX file prepared bythe author. MNRAS , 1–14 (2019) hemistry & transport in discs -1 R [au]10 -8 -7 -6 -5 -4 -3 a b un d a n c e ³ = 10 ¡ COCO CH H O C H R [au] ³ = 10 ¡ R [au] ³ = 10 ¡ -1 R [au]10 -8 -7 -6 -5 -4 -3 a b un d a n c e a ® = = 0 : ¹ m R [au] : M ¯ R [au] : M ¯ Figure A1.
Dependence of the abundance of the main carbon and oxygen carriers on the disc properties in models without transport.Top: Models with varying cosmic-ray ionization rate, but fixed disc mass, M d = . M (cid:12) : Bottom Left: A model where the grain size hasbeen fixed at . µ m in the chemical model. Bottom Middle and Right: Models with varying disc mass. -9 -8 -7 -6 -5 -4 a b un d a n c e CO ³ = 10 ¡ CH ³ = 10 ¡ N ³ = 10 ¡ Figure A2.
Time evolution of the CO, CH and N abundanceat 5 au in models with varying cosmic-ray ionization rates.MNRAS000