THE DUAL ORIGIN OF THE NITROGEN DEFICIENCY IN COMETS: SELECTIVE VOLATILE TRAPPING IN THE NEBULA AND POSTACCRETION RADIOGENIC HEATING

In our model, the volatile phase incorporated in cometesimals is composed of a mixture of pure ices, stoichiometric hydrates (such as NH₃–H₂O hydrate) and MG clathrates that crystallized in the form of microscopic grains at various temperatures in the outer part of the disk. Our model is based on the assumption that cometesimals have grown from the agglomeration of these icy grains due to collisional coagulation, implying no loss of volatile during their growth phase (Weidenschilling 1997). Here, the clathration process stops when no more crystalline water ice is available to trap the volatile species and then only pure condensates form at lower temperature in the disk. The process of volatiles trapping in icy grains is calculated using the equilibrium curves of hydrates and pure condensates, our model determining the equilibrium curves and compositions of MG clathrates, and the thermodynamic paths detailing the evolution of temperature and pressure between 5 and 30 AU in the protoplanetary disk. We refer to the work of Alibert et al. (2005) for a full description of the turbulent model of accretion disk used here.

The composition of the initial gas phase of the disk is defined as follows: we assume that the abundances of all elements (C, N, O, S, P, Ar, Kr, and Xe) are protosolar (Asplund et al. 2009) and that O, C, and N exist only under the form of H₂O, CO, CO₂, CH₃OH, CH₄, N₂, and NH₃. The abundances of CO, CO₂, CH₃OH, CH₄, N₂ and NH₃ are then determined from the adopted CO:CO₂:CH₃OH:CH₄ and N₂:NH₃ gas-phase molecular ratios. Once the abundances of these molecules are fixed, the remaining O gives the abundance of H₂O. Concerning the distribution of elements in the main volatile molecules, we set CO:CO₂:CH₃OH:CH₄ = 70:10:2:1 in the gas phase of the disk, values that are consistent with the interstellar medium (ISM) measurements considering the contributions of both gas and solid phases in the lines of sight (Mousis et al. 2009). In addition, S is assumed to exist in the form of H₂S, with H₂S:H₂ = 0.5 × (S/H₂)_☉, and other refractory sulfide components. We also consider N₂:NH₃ = 1:1 in the nebula gas phase, a value compatible with thermochemical calculations in the solar nebula that take into account catalytic effects of Fe grains on the kinetics of N₂ to NH₃ conversion (Fegley 2000). Note that CH₃OH is considered to be formed only as a pure condensate in our calculations since, to the best of our knowledge, no experimental data concerning the equilibrium curve of its associated clathrate have been reported in the literature.

The equilibrium pressure for the successive MG clathrates formed in the solar nebula can be expressed as (Hand et al. 2006)

$$\begin{equation} P_{{\rm eq,MG}} = \left[\sum _i \frac{y_i}{P_{{\rm eq},i}} \right]^{-1}, \end{equation} \tag{ 1 }$$

where y_i is the mole fraction of the component i in the fluid phase. The equilibrium pressure curves of each species are determined by fitting the available theoretical and laboratory data (Lunine & Stevenson 1985). Their equations are of the form log P_{eq, i} = A/T + B, where P_{eq, i} and T are the partial equilibrium pressure (bars) and temperature (K) of the considered species i, respectively.

The relative abundances of guest species incorporated in MG clathrate formed at a given temperature and pressure from the solar nebula gas phase are calculated following the method described by Lunine & Stevenson (1985) and Mousis et al. (2010), which uses classical statistical mechanics to relate the macroscopic thermodynamic properties of clathrates to the molecular structure and interaction energies. It is based on the original ideas of van der Waals & Platteeuw (1959) for clathrate formation, which assume that trapping of guest molecules into cages corresponds to the three-dimensional generalization of ideal localized adsorption. This approach is based on four key assumptions (Lunine & Stevenson 1985; Sloan & Koh 2008).

• 1.  
  The host molecule's contribution to the free energy is independent of the clathrate occupancy. This assumption implies in particular that the guest species do not distort the cages.
• 2.  
  (a) The cages are singly occupied. (b) Guest molecules rotate freely within the cage.
• 3.  
  Guest molecules do not interact with each other.
• 4.  
  Classical statistics are valid, i.e., quantum effects are negligible.

In this formalism, the fractional occupancy of a guest molecule K for a given type t (t = small or large) of cage can be written as

$$\begin{equation} y_{K,t}=\frac{C_{K,t}P_K}{1+\sum _{J}C_{J,t}P_J}, \end{equation} \tag{ 2 }$$

where the sum in the denominator includes all the species which are present in the initial gas phase. C_{K, t} is the Langmuir constant of species K in the cage of type t, and P_K is the partial pressure of species K. This partial pressure is given by P_K = x_K × P (we assume that the sample behaves as an ideal gas), with x_K the mole fraction of species K in the gas phase, and P the total gas pressure, which is dominated by H₂.

The Langmuir constant depends on the strength of the interaction between each guest species and each type of cage, and can be determined by integrating the molecular potential within the cavity as

$$\begin{equation} C_{K,t}=\frac{4\pi }{k_B T}\int _{0}^{R_c}\exp \Big (-\frac{w_{K,t}(r)}{k_B T}\Big)r^2dr, \end{equation} \tag{ 3 }$$

where R_c represents the radius of the cavity assumed to be spherical, k_B is the Boltzmann constant, and w_{K, t}(r) is the spherically averaged Kihara potential representing the interactions between the guest molecules K and the H₂O molecules forming the surrounding cage t. For a spherical guest molecule, this potential w(r) can be written as (McKoy & Sinanoğlu 1963)

$$\begin{eqnarray} w(r) &=& 2z\epsilon \Big [\frac{\sigma ^{12}}{R_c^{11}r}\Big (\delta ^{10}(r)+\frac{a}{R_c}\delta ^{11}(r)\Big)\nonumber\\&& - \frac{\sigma ^6}{R_c^5r}\Big (\delta ^4(r)+\frac{a}{R_c}\delta ^5(r)\Big)\Big ], \end{eqnarray} \tag{ 4 }$$

with

$$\begin{equation} \delta ^N(r)=\frac{1}{N}\Big [\Big (1-\frac{r}{R_c}-\frac{a}{R_c}\Big)^{-N}-\Big (1+\frac{r}{R_c}-\frac{a}{R_c}\Big)^{-N}\Big ]. \end{equation} \tag{ 5 }$$

In Equation (4), z is the coordination number of the cell. This parameter depends on the structure of the clathrate (I or II) and on the type of the cage (small or large). The Kihara parameters a, σ, and for the molecule–water interactions employed in this work are listed in Table 1. Note that our parameters are different from those used by Iro et al. (2003), who instead of using potential parameters describing the guest-clathrate interaction in their statistical model, employed potential parameters corresponding to guest–guest interactions (case of pure solutions). When comparing different data sets determined for the same X–H₂O interaction (with X = CH₄, N₂, and Xe), we opted for the most recent ones because they are fitted to a larger range of experimental data. In the case of Ar, Kr, CO, and PH₃, the listed sets are the unique ones that we found in the literature. For the H₂S–H₂O and CO₂–H₂O interaction parameters, we used those of Parrish & Prausnitz (1972) and Thomas et al. (2009; the set of CO₂–H₂O interaction parameters provided by Thomas et al. 2009 is an update of the one determined by Parrish & Prausnitz 1972 and very close to the one recently derived by Herri & Chassefière 2012). We also investigated the interaction parameters of Sloan & Koh (2008) but these parameters yield unphysical results when used in our models. We suspect this is because the code of Sloan & Koh (2008) is optimized for industrial purposes and the parameters are not appropriate for our temperature–pressure regime. We have tested our code against the one of Herri & Chassefière (2012) so we do not believe the problem resides in our code. Finally, the mole fraction f_K of a guest molecule K in a clathrate can be calculated with respect to the whole set of species considered in the system as

$$\begin{equation} f_K=\frac{b_s y_{K,s}+b_\ell y_{K,\ell }}{b_s \sum _J{y_{J,s}}+b_\ell \sum _J{y_{J,\ell }}}, \end{equation} \tag{ 6 }$$

where b_s and b_l are the number of small and large cages per unit cell, respectively, for the clathrate structure under consideration, and with ∑_K f_K = 1. Values of z, R_c, b_s, and b_l are taken from Parrish & Prausnitz (1972). Each time a species is fully trapped in an MG clathrate, its gaseous abundance is set to zero in the disk. This change of the disk's gas-phase composition induces the formation of a new MG clathrate at different temperature and pressure conditions and with a distinct composition. As a consequence, the compositions and formation conditions of successive MG clathrates are calculated n times in the nebula, with n corresponding to the total number of fully enclathrated species. Once all the water budget has been used for clathration, the volatiles remaining in the gas phase form pure condensates at lower temperatures. The equilibrium curves of these pure condensates derive from the compilation of laboratory data given in the CRC Handbook of Chemistry and Physics (Lide 2002).

Table 1. Adopted Set of Parameters for the Kihara Potential

Mol.	σ	/kB	a	Reference
	(Å)	(K)	(Å)
CO	3.1515	133.61	0.3976	Mohammadi et al. (2005)
CO2	2.9681	169.09	0.6805	Thomas et al. (2009)
CH4	3.14393	155.593	0.3834	Sloan & Koh (2008)
N2	3.13512	127.426	0.3526	Sloan & Koh (2008)
H2S	3.1558	205.85	0.36	Parrish & Prausnitz (1972)
PH3	3.771	275.0	0	Vorotyntsev & Malyshev (1998)
Ar	2.9434	170.50	0.184	Parrish & Prausnitz (1972)
Kr	2.9739	198.34	0.230	Parrish & Prausnitz (1972)
Xe	3.32968	193.708	0.2357	Sloan & Koh (2008)

Note. σ is the Lennard–Jones diameter, is the depth of the potential well, and a is the radius of the impenetrable core.

Download table as:  ASCIITypeset image

The composition of the volatile phase incorporated into cometesimals formed at a given distance from the Sun is finally given by the intersection of the disk's thermodynamic path at this location with the equilibrium curves of the different formed ices. The amount of volatile, i (relative to water), used to form the ice j in the nebula (either in the form of stoichiometric hydrate, MG clathrate, or of pure condensate), can be determined by the relation

$$\begin{equation} {m_{i,j} = \frac{X_{i,j}}{X_{{\rm H}_2{\rm O}}} \frac{\Sigma (r; T_j, P_j)}{\Sigma (r; T_{{\rm H}_2{\rm O}}, P_{{\rm H}_2{\rm O}})}}, \end{equation} \tag{ 7 }$$

where X_{i, j} is the mass mixing ratio with respect to H₂ of the volatile i used to form ice j in the nebula. $X_{{\rm H}_2{\rm O}}$ is the mass mixing ratio of H₂O with respect to H₂ in the nebula. Σ(R; T_j, P_j) and $\Sigma (R; T_{{\rm H}_2{\rm O}}, P_{{\rm H}_2{\rm O}})$ are the surface density of the disk at distance R from the Sun at the epoch of formation of ice j, and at the epoch of condensation of water, respectively. The global volatile, i, to water mass ratio incorporated in cometesimals is then given by M_j = ∑_{j = 1, k}m_j, with k corresponding to the number of solid phases in which volatile i is incorporated. Note that a thermodynamic path has been arbitrarily selected at the heliocentric distance of 5 AU to determine the composition of the formed ices. The adoption of any other distance range between ∼5 and 30 AU in the nebula would not affect the composition of the ices because it remains almost identical irrespective of (1) their formation distance and (2) the input parameters of the disk, provided that the initial gas-phase composition is homogeneous (Marboeuf et al. 2008).

Figure 1 represents the composition of the volatile phase incorporated in cometesimals agglomerated from ices condensed in the outer nebula and Figure 2 displays the molar ratios between incorporated volatiles of interest. Both figures show calculations expressed as a function of the formation temperature of cometesimals in the primordial nebula. Observations of active comets when their activity is driven by the sublimation of water ice—i.e., for heliocentric distances smaller than about 2 AU—suggest that the mean CO production rate relative to water is typically a few percent, with extreme values comprised between 0.4% and 30% (Bockelée-Morvan et al. 2004). Assuming that the production rates relative to water measured in comets are representative of the bulk composition of cometesimals, Figure 2 shows that their formation temperature must be lower than ∼47 K in the nebula to get CO/H₂O > 1% in their interior. The figure also shows that a formation temperature lower than ∼25 K would increase the CO/H₂O ratio up to ∼13% in cometesimals, as a result of the incorporation of pure CO condensate at this lower nebular temperature.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Formation temperatures of ices lower than ∼22 K in the nebula allow trapping of N₂ with N₂/H₂O of the order of 1% in the interior of planetesimals, a value far above the highest estimates (∼0.02%) inferred from observations (Wyckoff et al. 1991; Cochran et al. 2000). On the other hand, assuming a formation temperature of cometesimals in the 22–47 K range in the nebula would allow N₂/H₂O ∼0.01% and CO/H₂O ∼6.6% in their interior, values that are consistent with observations. A simple explanation for the formation of cometesimals would then argue that they need to be agglomerated from grains formed in the 22–47 K range in order to match the observations of the cometary production rates. However, this mechanism is not consistent with the fact that Pluto and Triton possess N₂-rich ice covers while they are expected to be formed in the same region as comets. An alternative possibility would be that comets accreted from planetesimals formed at very low temperature in the disk, but that subsequent radiogenic heating enabled their partial devolatilization through their pores. In the next section, we explore how radiogenic heating can induce important losses of N₂ in comets after their formation.

We use a fully three-dimensional model of heat transport described by Guilbert-Lepoutre et al. (2011) to compute the time evolution of the temperature distribution within cometary nuclei, accounting for the heating due to the radioactive decay of short-lived nuclides ²⁶Al and ⁶⁰Fe. We assume that these objects formed by cold accretion, i.e., no accretional heating is accounted for. The model computes the temperature distribution and its evolution with time, by solving the heat diffusion equation

$$\begin{equation} \rho _{{\rm bulk}}c \frac{\partial T}{\partial t} + \nabla (- \kappa \overrightarrow{\nabla }T) =\mathcal {Q}_{{\rm rad}} \end{equation} \tag{ 8 }$$

with ρ_bulk [kg m⁻³] the bulk density, c [J kg⁻¹ K⁻¹] the heat capacity, κ [W m⁻¹ K⁻¹] the thermal conductivity, and $\mathcal {Q}_{{\rm rad}}$ [W m⁻³] the internal power production per unit volume due to the decay of short-lived radioactive nuclides. As an example, we show in this paper the results for an idealized "Hale–Bopp" with a size of 35 km (Szabo et al. 2011; Lamy et al. 2004; Fernandez 2000; Weaver & Lamy 1997). We consider this object as a sphere, made of a porous mixture of crystalline water ice and dust, homogeneously distributed within the ice matrix. We use generic parameters to describe the bulk material: we consider that it can be described with a dust to water ice mass ratio X_d/X_w = 1 (with X_d and X_w the mass fractions of dust and water, respectively), a porosity ψ = 50%, and a bulk density ρ_bulk = 700 kg m⁻³ (typical values of the comet literature; see Huebner et al. 2006). The heat capacity and thermal conductivity of such a mixture are c = 800 J kg⁻¹ K⁻¹ and κ = 0.09 W m⁻¹ K⁻¹, respectively. The thermal evolution is computed with the formation time t_F—time for an object to grow to its final size—as a reference for time zero, from an initial temperature of 20 K.

Heating due to radioactive decay is described by

$$\begin{equation} \mathcal {Q}_{{\rm rad}}= \sum _{{\rm rad}}\rho _{d} X_{{\rm rad}}(t_F) H_{{\rm rad}}\frac{1}{ \tau _{{\rm rad}}} \exp \left(\frac{-t}{\tau _{{\rm rad}}}\right), \end{equation} \tag{ 9 }$$

with ρ_d the dust bulk density, X_rad(t_F) the mass fraction of the given radioactive isotope after a formation time t_F, H_rad the heat released per unit mass upon decay (4.84 × 10¹² J kg⁻¹ and 5.04 × 10¹² J kg⁻¹ for ²⁶Al and ⁶⁰Fe, respectively; see Guilbert-Lepoutre et al. 2011), τ_rad its mean lifetime (using 1.05 Myr for ²⁶Al and 3.78 Myr for ⁶⁰Fe; Rugel et al. 2009), and t the time. The object formation time t_F is accounted for in the model through the decay of radioactive isotopes that occurs during accretion. No additional accretional heating is accounted for, though. Decay during the formation time t_F results in a decrease of each nuclide initial abundance in the simulations:

$$\begin{equation} X_{{\rm rad}}(t_F)=X_{{\rm rad}}(0) e^{-t_F/\tau _{{\rm rad}}}, \end{equation} \tag{ 10 }$$

with X_rad(0) the initial mass fraction of each nuclide. We consider the initial nuclide ratios ²⁶Al/²⁷Al = 5 × 10⁻⁵ (MacPherson et al. 1995) and ⁶⁰Fe/⁵⁶Fe = 1 × 10⁻⁸ (Spivak-Birndorf et al. 2011), and mass fractions from Wasson & Kallemeyn (1988) for CV chondrites.

Although the model has many free parameters, only a few can actually strongly affect the resulting thermal history, such as the formation time of the object, its composition or the thermal conductivity. We also point out that the object's size plays an important role in its thermal evolution, owing to the surface to volume ratio. Figure 3 represents the influence of the formation time, through the abundance of radiogenic nuclides. It shows the evolution of the object's central temperature as a function of time after formation, for different formation times. The central temperature corresponds to the maximum temperature achieved inside the object: since the surface is in thermal equilibrium with the surrounding nebula at these stages of evolution, the heating due to radioactive decay is balanced by a cold load coming from the surface. Figure 4 displays the time evolution of the temperature profile within the body, calculated from a formation delay of 3 Myr. This figure shows that, at a given epoch of its evolution and for the adopted set of input parameters, the condition that the temperature profile must be in the 22–47 K range can be satisfied in almost the whole body, thus enabling the release of N₂ from this zone. Only the very close subsurface (down to a few hundreds meters deep) remains at lower temperatures, but this region should be rapidly depleted in volatiles once the comet follows subsequent orbits closer to the Sun.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The combined effects of composition and porosity are presented in Figure 5. For those simulations, we also considered a plausible formation time of 3 Myr and a fixed density of 700 kg m⁻³. Figure 6 represents the influence of thermal conductivity for the standard parameters previously mentioned. All together, these simulations show that it is possible to heat planetesimals at temperatures high enough to simultaneously imply the loss of N₂ and the preservation of CO, between 22 and 47 K.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size