Late Delivery of Nitrogen to the Earth

Material in an accretion disk orbits the central star of mass, M_⋆, with Keplerian velocity at radius R with angular velocity ${\rm{\Omega }}=\sqrt{{{GM}}_{\star }/{R}^{3}}$ (Pringle 1981). The viscosity in a fully MRI turbulent disk is

$$\begin{eqnarray}&&\nu ={\alpha }_{{\rm{m}}}\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{\rm{\Omega }}},\end{eqnarray} \tag{ 2 }$$

where α_m is the viscosity parameter (Shakura & Sunyaev 1973) and the sound speed is ${c}_{{\rm{s}}}=\sqrt{{ \mathcal R }{T}_{{\rm{c}}}/\mu }$, where ${ \mathcal R }$ is the gas constant, μ is the gas mean molecular weight, and T_c is the midplane temperature. With a constant accretion rate, $\dot{M}$, through all radii, the surface density of a steady-state disk is

$$\begin{eqnarray}&&{\rm{\Sigma }}=\displaystyle \frac{\dot{M}}{3\pi \nu }\end{eqnarray} \tag{ 3 }$$

(Pringle 1981). The surface temperature is determined through the balance of heating from viscous effects and irradiation from the central star

$$\begin{eqnarray}&&\sigma {T}_{{\rm{e}}}^{4}=\displaystyle \frac{9}{8}\displaystyle \frac{\dot{M}}{3\pi }{{\rm{\Omega }}}^{2}+\sigma {T}_{\mathrm{irr}}^{4}\end{eqnarray} \tag{ 4 }$$

(e.g., Cannizzo 1993; Pringle et al. 1986), where σ is the Stefan–Boltzmann constant and the irradiation temperature is

$$\begin{eqnarray}&&{T}_{\mathrm{irr}}={\left(\displaystyle \frac{2}{3\pi }\right)}^{\tfrac{1}{4}}{\left(\displaystyle \frac{{R}_{\star }}{R}\right)}^{\tfrac{3}{4}}{T}_{\star }\end{eqnarray} \tag{ 5 }$$

(Chiang & Goldreich 1997), where T_⋆ is the temperature and R_⋆ is the radius of the star. The midplane temperature is related to the surface temperature through

$$\begin{eqnarray}&&{T}_{{\rm{c}}}^{4}=\tau {T}_{{\rm{e}}}^{4},\end{eqnarray} \tag{ 6 }$$

where the optical depth, τ, is given by

$$\begin{eqnarray}&&\tau =\displaystyle \frac{3}{8}\kappa \displaystyle \frac{{\rm{\Sigma }}}{2}\end{eqnarray} \tag{ 7 }$$

and the opacity is $\kappa ={{aT}}_{{\rm{c}}}^{{\rm{b}}}$. In our model, we take a = 0.0001 and b = 2.1 as dust dominates the absorption properties of matter in the disk (Armitage et al. 2001).

The snow line radius in the disk is found by solving ${T}_{{\rm{c}}}={T}_{\mathrm{snow}}$. In our model, we assume ${M}_{\star }=1\,{M}_{\odot }$, ${R}_{\star }=3\,{R}_{\odot }$, and T_⋆ =  4000 K. There remains some uncertainty in the value of the viscosity parameter (e.g., King et al. 2007). Observations of FU Ori show ${\alpha }_{{\rm{m}}}\approx 0.01$ (Zhu et al. 2007), while from observations of X-ray binaries and cataclysmic variables, it is estimated that α_m ≈ 0.1 ∼ 0.4 (King et al. 2007). Since we consider young stellar objects, we assume that α_m = 0.01 in a fully turbulent MRI disk.

We adapt the approximate analytical steady-state solution of the CO snow line in Martin & Livio (2014) and find analytic fits to the nitrogen snow line radius

$$\begin{eqnarray}\begin{array}{rcl}{R}_{{{\rm{N}}}_{2},\mathrm{snow}} & \approx & 3.6{\left(\displaystyle \frac{{\alpha }_{{\rm{m}}}}{0.01}\right)}^{-\tfrac{2}{9}}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{\tfrac{1}{9}}{\left(\displaystyle \frac{\dot{M}}{{10}^{-8}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{\tfrac{2}{9}}\\ & & \times \ {\left(\displaystyle \frac{{T}_{{{\rm{N}}}_{2},\mathrm{snow}}}{58{\rm{K}}}\right)}^{-0.64}{\left(\displaystyle \frac{{R}_{\star }}{3{R}_{\odot }}\right)}^{\tfrac{2}{3}}{\left(\displaystyle \frac{{T}_{\star }}{4000{\rm{K}}}\right)}^{\tfrac{8}{9}}\,\mathrm{au}\end{array}\end{eqnarray} \tag{ 8 }$$

and the ammonia snow line radius

$$\begin{eqnarray}\begin{array}{rcl}{R}_{{\mathrm{NH}}_{3},\mathrm{snow}} & \approx & 2.1{\left(\displaystyle \frac{{\alpha }_{{\rm{m}}}}{0.01}\right)}^{-\tfrac{2}{9}}{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{\tfrac{1}{9}}{\left(\displaystyle \frac{\dot{M}}{{10}^{-8}{M}_{\odot }{\mathrm{yr}}^{-1}}\right)}^{\tfrac{2}{9}}\\ & & \times \ {\left(\displaystyle \frac{{T}_{{\mathrm{NH}}_{3},\mathrm{snow}}}{131{\rm{K}}}\right)}^{-0.64}{\left(\displaystyle \frac{{R}_{\star }}{3{R}_{\odot }}\right)}^{\tfrac{2}{3}}{\left(\displaystyle \frac{{T}_{\star }}{4000{\rm{K}}}\right)}^{\tfrac{8}{9}}\,\mathrm{au}.\end{array}\end{eqnarray} \tag{ 9 }$$

These are valid for relatively low accretion rates, $\dot{M}\,\lesssim {10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$, where irradiation dominates the viscous heating term.

The numerical values of the snow line radii in Equations (8) and (9), for the standard parameters chosen, are uncomfortably small. The ammonia snow line lies inside of the main asteroid belt, while the nitrogen snow line lies within it. Figure 1 shows the full numerical solution for the snow line radii for nitrogen and ammonia in a fully turbulent disk as a function of the steady-state accretion rate (dashed lines) found by solving T_c = T_snow along with Equations (4)–(6). This exposes the problem with fully turbulent disks, in that the snow line moves in too close to the host star during the low accretion rate phase when the disk is near the end of its lifetime (Garaud & Lin 2007; Martin & Lubow 2011; Oka et al. 2011; Martin & Livio 2014). This contradicts observations of the composition of asteroids in the asteroid belt (e.g., DeMeo & Carry 2014). To explore possible solutions to this problem, we consider a more realistic disk—one with a dead zone—in the following section.

Zoom In Zoom Out Reset image size

The dashed lines in Figure 2 show the surface density at the snow line radius as a function of the accretion rate. For low accretion rates, the surface density at the snow line radius becomes very small, ${\rm{\Sigma }}\lesssim 10\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. Assuming a gas to solid ratio of 100 (e.g., Bohlin et al. 1978; Hayashi 1981), the surface density of solid material is $\lesssim 0.1\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. In the fully turbulent disk model there is little solid material available for planetesimal formation in the outer parts of the disk.

Zoom In Zoom Out Reset image size

The MRI requires a sufficiently high level of ionization in the disk for it to operate. The inner regions of a protoplanetary disk are hot enough to be thermally ionized where ${T}_{{\rm{c}}}\geqslant {T}_{\mathrm{crit}}$. The critical temperature, T_crit, is thought to be around 800 K for activating the MRI (Umebayashi & Nakano 1988). However, farther from the host star, external sources such as cosmic rays and X-rays from the central star ionize only the surface layers up to a critical surface density, Σ_crit. The value of the critical surface density is not well determined (Martin et al. 2012a, 2012b). For T Tauri stars, cosmic rays ionize a critical surface density of ${{\rm{\Sigma }}}_{\mathrm{crit}}\approx 100\,{\rm{g}}\,{\mathrm{cm}}^{2}$ (Gammie 1996). However, this value could be much lower if we consider the chemical reaction of charged particles (Sano et al. 2000). Bai & Goodman (2009) and Turner & Sano (2008) found that the MRI could operate in the surface layer of a disk ionized by stellar X-rays with ${{\rm{\Sigma }}}_{\mathrm{crit}}\approx 1\mbox{--}30\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. Where the surface density is larger than this, ${\rm{\Sigma }}\geqslant {{\rm{\Sigma }}}_{\mathrm{crit}}$, a dead zone forms at the disk midplane where the MRI cannot operate (e.g., Gammie 1996; Gammie & Menou 1998).

In the dead zone layer, material builds up due to the lower viscosity and the small mass transport rate. The disk becomes self-gravitating when the Toomre (1964) parameter, $Q={c}_{{\rm{s}}}\,{\rm{\Omega }}/\,\pi G{\rm{\Sigma }}$, is less than the critical, Q_crit = 2. Then, a second type of turbulence, gravitational turbulence, is driven, with a viscosity given by

$$\begin{eqnarray}&&{\nu }_{{\rm{g}}}={\alpha }_{{\rm{g}}}\displaystyle \frac{{c}_{{\rm{g}}}^{2}}{{\rm{\Omega }}},\end{eqnarray} \tag{ 10 }$$

and the alpha parameter is

$$\begin{eqnarray}&&{\alpha }_{{\rm{g}}}={\alpha }_{{\rm{m}}}\exp (-{Q}^{4})\end{eqnarray} \tag{ 11 }$$

(e.g., Zhu et al. 2010). The disk can reach a steady state with a self-gravitating dead zone (e.g., Martin & Lubow 2013, 2014; Rafikov 2015). Martin & Livio (2012) found that the evolution of the snow line is significantly altered in a time-dependent disk with a dead zone compared with a fully turbulent disk model. The build up of material can lead to the disk being gravo-magneto unstable (see the disk structure sketches in Martin & Livio 2013). The temperature in the dead zone increases sufficiently for the MRI to be triggered there, leading to an accretion outburst. However, the outer parts of the disk that we are interested in here remain in a steady state. Thus, we consider only steady-state self-gravitating disk solutions in this work.

We follow Martin & Livio (2012) to obtain an analytic solution for the snow line radii assuming that both the ammonia and nitrogen snow lines are in the self-gravitating part of the dead zone. The disk has MRI active layers over a self-gravitating dead zone, but we assume that ${\rm{\Sigma }}\gg {{\rm{\Sigma }}}_{\mathrm{crit}}\approx 0$. The self-gravitating disk has a surface density

$$\begin{eqnarray}&&{\rm{\Sigma }}=\displaystyle \frac{{c}_{{\rm{g}}}{\rm{\Omega }}}{\pi {GQ}}.\end{eqnarray} \tag{ 12 }$$

The accretion rate of the steady disk is

$$\begin{eqnarray}&&\dot{M}=3\pi {\nu }_{{\rm{g}}}{\rm{\Sigma }}=\left(\displaystyle \frac{3{c}_{{\rm{g}}}^{3}{\alpha }_{{\rm{m}}}}{G}\right)\displaystyle \frac{\exp (-{Q}^{4})}{Q}.\end{eqnarray} \tag{ 13 }$$

The term in brackets is constant for a fixed snow line temperature. This expression depends sensitively on Q, but Q is approximately constant in the range of accretion rates we use. We scale the variables to ${T}_{{{\rm{N}}}_{2},\mathrm{snow}}^{{\prime} }={T}_{{{\rm{N}}}_{2},\mathrm{snow}}/58\,{\rm{K}}$, ${T}_{{\mathrm{NH}}_{3},\mathrm{snow}}^{{\prime} }={T}_{{\mathrm{NH}}_{3},\mathrm{snow}}/131\,{\rm{K}}$, ${\alpha }_{{\rm{m}}}^{{\prime} }={\alpha }_{{\rm{m}}}/0.01$, ${M}_{\star }^{{\prime} }={M}_{\star }/{M}_{\odot }$, ${R}_{\star }^{{\prime} }={R}_{\star }/3\,{R}_{\odot }$, and $R^{\prime} =R/\mathrm{au}$. Then, we can solve Equation (13) to calculate the scaled Toomre parameter

$$\begin{eqnarray}&&Q^{\prime} =\displaystyle \frac{Q}{{Q}_{\mathrm{crit}}}=0.69{\left[\displaystyle \frac{W(x)}{W({x}_{0})}\right]}^{\tfrac{1}{4}},\end{eqnarray} \tag{ 14 }$$

where we define

$$\begin{eqnarray}&&{x}_{{{\rm{N}}}_{2}}={x}_{0,{{\rm{N}}}_{2}}\displaystyle \frac{{\alpha }_{{\rm{m}}}^{{\prime} 4}{T}_{{{\rm{N}}}_{2},\mathrm{snow}}^{{\prime} 6}}{{M}_{\star }^{{\prime} 4}}\end{eqnarray} \tag{ 15 }$$

with ${x}_{0,{{\rm{N}}}_{2}}=3.86\times {10}^{10}$ and

$$\begin{eqnarray}&&{x}_{{\mathrm{NH}}_{3}}={x}_{0,{\mathrm{NH}}_{3}}\displaystyle \frac{{\alpha }_{{\rm{m}}}^{{\prime} 4}{T}_{{\mathrm{NH}}_{3},\mathrm{snow}}^{{\prime} 6}}{{M}_{\star }^{{\prime} 4}}\end{eqnarray} \tag{ 16 }$$

with ${x}_{0,{\mathrm{NH}}_{3}}=5.13\times {10}^{12}$. The equation

$$\begin{eqnarray}&&x=W(x)\exp [W(x)]\end{eqnarray} \tag{ 17 }$$

defines the Lambert function, W.

The steady-state energy equation of a disk with a self-gravitating dead zone is

$$\begin{eqnarray}&&\sigma {T}_{{\rm{e}}}^{4}=\displaystyle \frac{9}{8}{\nu }_{{\rm{g}}}{\rm{\Sigma }}{{\rm{\Omega }}}^{2}+\sigma {T}_{\mathrm{irr}}^{4}.\end{eqnarray} \tag{ 18 }$$

The midplane temperature of the disk is related to the disk surface temperature with Equations (5)–(7). In Figure 1, the two solid lines show the snow line radii that are the solutions to T_c = T_snow with Equation (18). The snow line radius in the disk model with a dead zone is farther from the central star than in the fully turbulent disk model because the temperature of the disk is higher. With a dead zone, for low accretion rates the snow line radius is insensitive to the accretion rate. The self-gravitating dead zone prevents the snow lines from moving too close to the central star in the late stages of disk evolution.

For lower accretion rates ($\dot{M}\leqslant {10}^{-8}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$), irradiation dominates over the viscous heating term, and we approximate Equation (18) with T_e = T_irr and find the analytic snow line radius for nitrogen to be

$$\begin{eqnarray}&&{R}_{{{\rm{N}}}_{2},\mathrm{snow}}=12.6\,{M}_{\star }^{{\prime} \tfrac{1}{9}}{R}_{\star }^{{\prime} \tfrac{2}{3}}{T}_{\star }^{{\prime} \tfrac{8}{9}}{T}_{{{\rm{N}}}_{2},\mathrm{snow}}^{{\prime} -0.31}{\left[\displaystyle \frac{W(x)}{W({x}_{0})}\right]}^{-\tfrac{1}{18}}\,{\rm{a}}{\rm{u}}\end{eqnarray} \tag{ 19 }$$

and for ammonia to be

$$\begin{eqnarray}&&{R}_{{\mathrm{NH}}_{3},\mathrm{snow}}=9.8\,{M}_{\star }^{{\prime} \tfrac{1}{9}}{R}_{\star }^{{\prime} \tfrac{2}{3}}{T}_{\star }^{{\prime} \tfrac{8}{9}}{T}_{{\mathrm{NH}}_{3},\mathrm{snow}}^{{\prime} -0.31}{\left[\displaystyle \frac{W(x)}{W({x}_{0})}\right]}^{-\tfrac{1}{18}}\mathrm{au}.\end{eqnarray} \tag{ 20 }$$

These values are given approximately by

$$\begin{eqnarray}\begin{array}{rcl}{R}_{{{\rm{N}}}_{2},\mathrm{snow}} & = & \,12.6{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{\tfrac{1}{9}}{\left(\displaystyle \frac{{R}_{\star }}{3{R}_{\odot }}\right)}^{\tfrac{2}{3}}{\left(\displaystyle \frac{{T}_{\star }}{4000{\rm{K}}}\right)}^{\tfrac{8}{9}}\\ & & \times \ {\left(\displaystyle \frac{{T}_{{{\rm{N}}}_{2},\mathrm{snow}}}{58{\rm{K}}}\right)}^{-0.31}\,\mathrm{au}\end{array}\end{eqnarray} \tag{ 21 }$$

and

$$\begin{eqnarray}\begin{array}{rcl}{R}_{{\mathrm{NH}}_{3},\mathrm{snow}} & = & \,9.8{\left(\displaystyle \frac{{M}_{\star }}{{M}_{\odot }}\right)}^{\tfrac{1}{9}}{\left(\displaystyle \frac{{R}_{\star }}{3{R}_{\odot }}\right)}^{\tfrac{2}{3}}{\left(\displaystyle \frac{{T}_{\star }}{4000{\rm{K}}}\right)}^{\tfrac{8}{9}}\\ & & \times \ {\left(\displaystyle \frac{{T}_{{\mathrm{NH}}_{3},\mathrm{snow}}}{131{\rm{K}}}\right)}^{-0.31}\,\mathrm{au}.\end{array}\end{eqnarray} \tag{ 22 }$$

The numerical values of Equations (21) and (22) for our standard parameters are consistent with the solid lines in Figure 1 at low accretion rates.

The solid lines in Figure 2 show the surface density of the disk with a self-gravitating dead zone at the snow line radius, for nitrogen (blue line) and for ammonia (green line). The disk with the dead zone has a roughly constant surface density with accretion rate, particularly for low accretion rates. The surface density is much higher than even the maximum possible value for the critical surface density that may be MRI active, ${{\rm{\Sigma }}}_{\mathrm{crit}}\approx 200\,{\rm{g}}\,{\mathrm{cm}}^{-2}$. Thus, our assumption of ${\rm{\Sigma }}\gg {{\rm{\Sigma }}}_{\mathrm{crit}}$ is justified. Furthermore, the surface density at the snow line radii is significantly higher in a disk with a dead zone, thus allowing for more solid bodies to form.

For the low accretion rates expected at the end of the disk lifetime, the snow line radii for nitrogen and ammonia in the solar system are both farther from the Sun than the main asteroid belt. The ammonia snow line lies at about 9 au, while the nitrogen snow line is at 12 au. Thus, comets in the Kuiper Belt are expected to form with significantly more nitrogen than asteroids in the main asteroid belt. In the next section, we consider observational constraints from the isotopic ratio of nitrogen and estimate the fraction of the present Earth's atmosphere that could have been delivered by comets.

In the context of the Kuiper Belt, interior mean-motion resonances with Neptune generally produce scattering events from eccentricity excitations compared with exterior mean-motion resonances that provide resonant capture of debris (Beauge & Ferraz-Mello 1994). Thus, we narrow our investigation of a nitrogen transport mechanism to secular resonances rather than including mean-motion resonances, as there are no interior mean-motion resonances within the Kuiper Belt. There are two types of secular resonances, apsidal and nodal. An apsidal resonance excites eccentricity, while a nodal resonance excites inclination (Froeschle et al. 1991). Apsidal resonances are the most important secular resonances for transporting material to the Earth, as the comets remain in the orbital plane of the system's solid large bodies with an increased eccentricity and are more likely to collide with the Earth. Here, we apply a first-order secular resonance model to the solar system to find all of the apsidal resonance locations. We are particularly interested in secular resonances that are outside of the ammonia snow line (at around 9 au at the end of the gas disk lifetime; see Section 2) as objects outside of this location may contain large amounts of nitrogen.

A secular resonance arises when the free precession rate of a test particle is equal to the eigenfrequency of one of the planets (e.g., Froeschle & Scholl 1986; Yoshikawa 1987; Minton & Malhotra 2011; Smallwood et al. 2018a, 2018b). The planetary eigenfrequencies are obtained by calculating the eigenvalues from the Laplace–Lagrange equations in the context of the generalized secular perturbation theory (see, for example, Murray & Dermott 1999). The eigenfrequencies for each of the solar system planets are shown in Table 1. The free precession rate of a test particle with orbital frequency, n, and semimajor axis, a, in the potential of the solar system is

$$\begin{eqnarray}&&\dot{\omega }=\displaystyle \frac{n}{4}\displaystyle \sum _{j=1}^{N}\displaystyle \frac{{m}_{j}}{{m}_{* }}{\alpha }_{j}{\bar{\alpha }}_{j}{b}_{3/2}^{(1)}({\alpha }_{j})\end{eqnarray} \tag{ 27 }$$

(Murray & Dermott 1999), where m_j is the mass of the solar system planets, m_* is the mass of the Sun, b^(j)_s is the Laplace coefficient, and the coefficients α_j and ${\bar{\alpha }}_{j}$ are defined as

$$\begin{eqnarray}{\alpha }_{j}=\left\{\begin{array}{ll}{a}_{j}/a, & \mathrm{if}\ {a}_{j}\lt a,\\ a/{a}_{j}, & \mathrm{if}\ {a}_{j}\gt a,\end{array}\right.\end{eqnarray} \tag{ 28 }$$

and

$$\begin{eqnarray}{\bar{\alpha }}_{j}=\left\{\begin{array}{ll}1, & \mathrm{if}\ {a}_{j}\lt a,\\ a/{a}_{j}, & \mathrm{if}\ {a}_{j}\gt a.\end{array}\right.\end{eqnarray} \tag{ 29 }$$

The top panel of Figure 4 shows the eigenfrequencies of the planets as horizontal dashed lines. The test particle free precession rate is represented by the solid lines. Where the solid lines cross a horizontal line there is an apsidal secular resonance.

Zoom In Zoom Out Reset image size

We now consider the strength of the apsidal resonances. In the bottom panel of Figure 4, we show the forced eccentricity of a test particle in the solar system (for details on these calculations see Murray & Dermott 1999 and Minton & Malhotra 2011). Planetary debris that falls within the high eccentricity parts of the gray shaded regions undergoes eccentricity growth that leads to ejection or collision with a larger object. In the context of nitrogen transport to the Earth, the source of nitrogen may originate from the Kuiper Belt, which extends roughly from 30 to 50 au (Duncan et al. 1995). The ν₈ secular resonance is the most prominent secular resonance within the Kuiper Belt (Froeschle & Scholl 1986). It is produced when the free precession rate of the comets is close to the eigenfrequency of Neptune. According to our linear secular resonance model, the ν₈ resonance is located at 40.9 au (the resonance location is based on first-order approximations), which is consistent with the actual location of the resonance (Nagasawa & Ida 2000).

Populations of Kuiper Belt objects that are in resonance suggest that Uranus and Neptune underwent outward migration (Fernandez & Ip 1984; Malhotra 1993, 1995; Hahn & Malhotra 1999, 2005; Chiang & Jordan 2002; Chiang et al. 2003; Murray-Clay & Chiang 2005; Tsiganis et al. 2005; Morbidelli et al. 2007; Nesvorný 2015; Nesvorný & Vokrouhlický 2016). During this dynamical phase, resonant bodies in the Kuiper Belt were destabilized and there was a sudden massive delivery of planetesimals to the inner solar system (e.g., Levison & Morbidelli 2003; Gomes et al. 2005; Nesvorný et al. 2017). In particular, the ν₈ secular resonance swept through the inner parts of the Kuiper Belt (Nagasawa & Ida 2000) and caused a large fraction of the comets to acquire free precession frequencies similar to that of Neptune. Comets that originated or got "bumped" into this sweeping resonance would have undergone increased eccentricity excitations leading to collisions with Earth, allowing for the transport of nitrogen.

To study the delivery of comets to the Earth, we use the N-body simulation package, REBOUND, to model the Kuiper Belt including the Earth, the four giant planets, and the Sun in the solar system (Rein & Liu 2012). We use a hybrid integrator called mercurius. This integrator uses WHFast, which is an implementation of a Wisdom–Holman symplectic integrator to simulate the long-term orbit integrations of the planetary system. For close encounters of test particles and planets, it switches smoothly to IAS15, which is a 15th-order integrator, to simulate gravitational dynamics (Rein & Spiegel 2015; Rein & Tamayo 2015). Thus, we can calculate the evolution of the Kuiper Belt for a duration of one hundred million years (Myr) quickly and precisely. In our simulation, we take a snapshot every 0.1 Myr.

Comets interact with the Earth, the gas giant planets, and the Sun due to gravitational forces, while there is no interaction between comets. The interaction timescale of some of the largest asteroids is of the order of magnitude of the age of the solar system (Dohnanyi 1969).

The obit of each comet is defined initially by six orbital elements, $a,i,e,\omega ,{\rm{\Omega }},\nu$. In our model, we assume the semimajor axis, a, is distributed uniformly in the range of a_min = 38 to a_max = 45 au, including the region of the ν₈ secular resonance (see Figure 4). The uniform distribution is given by

$$\begin{eqnarray}&&{a}_{\mathrm{comet}}=({a}_{\max }-{a}_{\min })\times \,{\chi }_{r}+{a}_{\min }\end{eqnarray} \tag{ 30 }$$

(Lecar & Franklin 1997; Smallwood et al. 2018a), and χ_r is a randomly generated number between 0 and 1. The inclination, i, is randomly allocated in the range from 0° to 10° and the eccentricity, e, is randomly allocated in the range from 0 to 0.1. The argument of periapsis, ω, the longitude of the ascending node, Ω, and the true anomaly, ν, are all uniformly distributed in the range from 0° to 360°.

The present-day orbital elements for each of the planets are set as the initial parameters for the planets because the solar system is stable over long timescales (Duncan & Lissauer 1998; Ito & Malhotra 2006). REBOUND queries those data automatically from the NASA HORIZONS database. A test particle that has a semimajor axis beyond 200 au is considered to be ejected.

In our simulation, we inflate the size of the Earth to its Hill radius. The Hill radius is given by

$$\begin{eqnarray}&&{R}_{{\rm{H}}}={a}_{{\rm{p}}}{\left(\displaystyle \frac{{M}_{\mathrm{planet}}}{3{M}_{\odot }}\right)}^{1/3},\end{eqnarray} \tag{ 31 }$$

where M_planet is the mass of the planet and a_p is the semimajor axis of the planet (Delsanti & Jewitt 2006). The Hill radius for the Earth is 2 × 10⁶ km. In principle, 3 Hill radii is the range of the planet's gravitational force, and comets cannot remain on stable orbits within this region (Gladman 1993; Chatterjee et al. 2008; Morrison & Malhotra 2015). The increased size of the Earth artificially increases the rate of collisions with the Earth, allowing us to simulate fewer particles than would have been present in the early Kuiper Belt.

Figure 5 shows the initial eccentricity and semimajor axis distribution of 50,000 test particles. There are a total of 104 test particles that collide with the Earth within 100 Myr, and the initial conditions for these are shown by the blue dots. Seventy of these test particles originate from the region of the ν₈ resonance around 40–42 au. The 3:2 and 5:3 mean-motion resonances of Neptune that are at 39.5 and 42 au also contribute a few of the collisions with the Earth.

Zoom In Zoom Out Reset image size

On the other hand, a test particle that has a semimajor axis a > 43 au is quite stable, the region is regarded as having classic Kuiper Belt objects, and the final distribution of eccentricity and semimajor axis of the comets in Figure 6 is roughly consistent with that observed for Kuiper Belt objects. However, there are still many particles in the inner region where a < 40 au. Current observations show that this region is quite clear (see Figure 7 in Jewitt et al. 2009). As we do not consider the outward migration of Uranus and Neptune in our simulation, we do not reproduce this clearing. The region of the ν₈ resonance is not completely cleared out of comets yet. About 50% of the test particles that were initially distributed in the region of the ν₈ resonance are still in this region.

Zoom In Zoom Out Reset image size

Resonances in the Kuiper Belt excite the eccentricities of the test particles and this increases the chance of a close encounter with a planet. These close encounters result in outward or inward migration of test particles. Once a test particle comes into the inner solar system, the Earth has an opportunity to capture it. Figure 7 shows that particles begin colliding with the Earth after about 25 Myr, and there are more than four impact particles every 5 Myr after 30 Myr. In the next section, we use the results of our N-body simulation to estimate the amount of nitrogen that was delivered to the Earth from the Kuiper Belt.

Zoom In Zoom Out Reset image size

We now consider a simple calculation to scale our simulation result to calculate how much mass may have been delivered to the Earth by comets. Most of the impacts come from the ν₈ resonance region and we focus on these. The particles that begin in the region we are interested in, 40–42 au, assume a uniform distribution of the number of test particles in semimajor axis. However, we inflated the size of the Earth to increase the possibility of a collision. Smallwood et al. (2018a) calculated that the number of collisions is reduced by a factor of 800/5 if we shrink the size of the Earth from its Hill's radius to the radius of the Earth. Thus, the collision fraction should be 3 × 10⁻⁵ collisions per particle per 100 Myr.

The initial mass of the Kuiper Belt was at least 10 times the mass of Earth, ${M}_{\oplus }=6\times {10}^{27}\,{\rm{g}}$ (Delsanti & Jewitt 2006). The amount of mass that originated in the ν₈ resonance and collided with the Earth is approximately 5 × 10²³ g in 100 Myr. Because at the end of our simulation there are still 50% of the original particles in the region of the ν₈ resonance, we multiply the number of collisions by a factor of 2 to find the mass to be ≈1 × 10²⁴ g. For comparison, the total mass that was delivered to the Earth during LHB is estimated to be around $2\times {10}^{23}\mbox{--}2\times {10}^{24}\,{\rm{g}}$ by modeling the thermal effects of the LHB (Abramov et al. 2013). However, according to the constraint of atmospheric ³⁶Ar, Marty et al. (2016) estimated the cometary contribution to be about 2.0 × 10²² g. This value is 2 orders of the magnitude lower than our estimate. However, we have not taken into account the fact that a comet starts to sublime volatile when it enters the inner solar system. For example, the dynamical lifetime of the comet Halley is about 3000 perihelion passages or 0.2 Myr (Olsson-Steel 1988). It has lost about 0.4% of the total mass of the nucleus over the last 2200 yr (Schmude 2010). From observations, we have already found some small comets that have a radius <0.4 km that have lost of most their mass (see Table 9.4 in Fernández 2006). If we consider a comet that has a nucleus radius of 1 km, the bulk density is $0.5\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ and the perihelion distance is 1 au, and the dynamical lifetime is less than 1000 orbits (see Figure 9.2 in Fernández 2006). Due to the mass loss of comets, comets only possess 1% to 10% of their original mass when they collide with the Earth. Thus, the mass delivered by comets in our model is in the approximate range of 10²²–10²³ g.

To calculate how much nitrogen could have been delivered in these comets, we use data of 16 comets from Le Roy et al. (2015) and Paganini et al. (2012) to estimate the mean values of the highest relative abundances of NH₃ and HCN to be 0.7% and 0.28% with respect to water. Besides, in a recent study, ROSINA measured the N₂/CO ratio by mass to be 0.570% on the comet 67P/Churyumov–Gerasimenko, and CO in this comet is about 10% with respect to water (Le Roy et al. 2015). Thus, we use 0.057% for the relative abundance of N₂ with respect to water. H₂O takes up 50% of the composition of the comet by mass, and thus the total mass of nitrogen from comets is estimated to be in the range of 3.9 × 10¹⁹–3.9 × 10²⁰g.

The total mass of the Earth's atmosphere is 5.136 × 10²¹g (Weast et al. 1989). Nitrogen makes up 78% of the composition of the current Earth's atmosphere so that the mass of nitrogen in the atmosphere is 3.8 × 10²¹ g (Verniani 1966). If 10% of nitrogen was delivered by comets (see Section 3), the mass of cometary nitrogen is 4 × 10²⁰g, and our simulation is consistent with this value.

Note also that here we only consider the contribution from the ν₈ region. There would have been additional collisions as a result of mean-motion resonances of Neptune (see Figure 5) and the migration of Uranus and Neptune, which we did not include in our simulation. In the Archean era, the partial pressure of nitrogen was at least 1.4  ∼  1.6 times higher than today's atmosphere. Mallik et al. (2018) estimated this value by subtracting the net degassing flux of N from arcs and back-arc basins (Hilton et al. 2002), mid-ocean ridges and intraplate settings (Sano et al. 2001) from the influx of N at subduction zones today. Due to N recycling efficiency, some nitrogen in the Earth's atmosphere was brought to the subduction zones and the deep mantle (Mallik et al. 2018). However, there are some contradictory results that the partial pressure of the nitrogen in Archean atmosphere is similar or less than the present atmosphere's based on different analyses of Archean samples (Som et al. 2016; Avice et al. 2018). Consequently, it leaves some uncertainties, and our simulations cannot solve this problem.

Because comets also contain a large amount of water, they deliver both nitrogen and water, and the D/H ratio in the Earth may also be changed. Since the D/H ratios of comets are much higher than terrestrial water and enstatite chondrites do not contain large amounts of water, it could be used to constrain the contribution of water from the comet to the Earth and it should be less than about 10% (Dauphas et al. 2000; Bockelée-Morvan et al. 2004). Since water takes up 50% of a comet's composition, we estimate the amount of water that was delivered from comets, and this value is $5\times {10}^{21}\sim \ 5\times {10}^{22}$ g. The total mass of the water on the Earth is 1.35 × 10²⁴g. Therefore, the water contribution from comets is just about 0.4% ∼ 4%. This value is more in agreement with Marty et al. (2016), who state the cometary contribution of water is ≤1%, whatever the reservoir considered.