Abundances of Ordinary Chondrites in Thermally Evolving Planetesimals

Planetesimals are some of the central objects that play crucial roles in planet formation (e.g., Pollack et al. 1996; Ida & Lin 2004; Mordasini et al. 2009). It has been suggested for a long time that planetary cores were built from planetesimal collisions and resulting mergers in protoplanetary disks (e.g., Wetherill & Stewart 1989; Kokubo & Ida 1998; Morishima et al. 2010; Raymond et al. 2014). The presence of planetesimals in planet-forming regions can also affect the orbital evolution of protoplanets there, which is referred to as planetesimal-driven migration (e.g., Levison et al. 2010; Kominami et al. 2016). More recently, planetesimals have been regarded as the main contributors to metal enrichment of Jovian planets both in the solar system and extrasolar planetary systems (e.g., Saumon & Guillot 2004; Miller & Fortney 2011; Mordasini et al. 2016). These are all regulated by the dynamics of planetesimals, and hence are determined by the properties of planetesimals such as their abundances and masses. Thus, it is of fundamental importance to identify and characterize when and how planetesimal formation takes place in protoplanetary disks for fully exploring the formation and migration histories of planets.

Understanding of planetesimal formation is still far from complete despite recent progress (e.g., Johansen et al. 2014). There are two reasons for this. One is that, theoretically, no force (i.e., surface adhesion force or gravity) is effective for kilometer-sized bodies. Many theoretical studies thereby attempt to establish high density regions of solids in protoplanetary disks that eventually become gravitationally unstable and lead to direct and instant formation of planetesimals due to self-gravity (e.g., Cuzzi et al. 2001; Johansen et al. 2007; Kretke & Lin 2007). The other is that, observationally, planetesimals are invisible at almost all wavelengths in the nearby star-forming regions. This arises from the fact that emission from gas and dust and their scattered light become generally dominant in the observations. Investigating planetesimals is therefore not straightforward even for nearby, young stars.

The solar system stands out in exploring planetesimal formation because we can obtain key information from asteroids and meteorites. It is reasonable to consider that asteroids are the remnants of planetesimals that survived a number of constructive and destructive events occurring over the age of the solar system (e.g., Bottke et al. 2005a, 2005b; Morbidelli et al. 2009; Zheng et al. 2017; Tsirvoulis et al. 2018). In other words, asteroids are unique objects that contain a wealth of information about planetesimals. There are completed and ongoing sample return missions for investigating the composition of asteroids: Hayabusa, OSIRIS-REx, and Hayabusa 2. For instance, the Hayabusa revealed that the chemical and isotopic compositions of particles on the asteroid Itokawa are similar to those of meteorites falling onto the Earth (e.g., Nakamura et al. 2011; Yoshikawa et al. 2015). Asteroids can also be used as a probe to specify the size distribution of planetesimals. Indeed, theoretical studies suggest that primordial planetesimals would be larger than 100 km in diameter in order to explain the size distribution of current main belt asteroids (e.g., Bottke et al. 2005a; Morbidelli et al. 2009; Tsirvoulis et al. 2018). Thus, asteroids give a number of invaluable insights about planetesimals. It is however important to realize that these insights are derived from the current status of asteroids, and hence one has to take into account the full history of a number of events that occurred for asteroids.

Meteorites act as important and complimentary objects for understanding planetesimal formation (e.g., McSween 1999). This is because they literally contain the fossil record of how the solar system formed. For example, chondrites are known as the most primitive meteorites and contain the first condensed materials in the solar system, called Ca-Al-rich inclusions (CAIs; e.g., MacPherson 2005). It is well known that the age of the solar system is measured from CAIs, which is 4567.30 ± 0.16 million years (e.g., MacPherson et al. 1995; Amelin et al. 2010; Connelly et al. 2012; Davis & McKeegan 2014). Also, the detailed chemical and isotopic analyses of meteorites suggest that planetesimals should have formed soon after CAI formation and continued for about a few million years (Davis et al. 2014; Gail et al. 2014; Blackburn et al. 2017). Moreover, meteorites contain another primitive component called chondrules that are considered to have formed within 5 Myr after CAI formation (e.g., Connelly et al. 2012; Bollard et al. 2015). The formation mechanisms of chondrules are currently under active investigation, one of which is highly relevant to planetesimal collisions (e.g., Johnson et al. 2015; Hasegawa et al. 2016a, 2016b; Matsumoto et al. 2017; Wakita et al. 2017a).

Another remarkable feature of meteorites is that more than 40,000 of them have been found. Interestingly, 90% of the analyzed meteorites are ordinary chondrites (Meteoritical Bulletin website http://www.lpi.usra.edu/meteor/metbull.php; Grady et al. 2014). Thanks to such abundant samples, statistical studies have become possible, including classification. In general, meteorites are grouped by their chemical composition and texture (e.g., Van Schmus & Wood 1967; Krot et al. 2005; Scott & Krot 2005; Weisberg et al. 2006). We focus here on the petrologic type that links the degree of thermal metamorphism. The degree of metamorphism is considered to reflect the maximum temperatures that the parent bodies of meteorites experienced (e.g., Dodd 1981; Scott & Krot 2005; Huss et al. 2006). Type 3 chondrites are the least metamorphosed ones and would not pass through temperatures higher than ∼600 °C. Type 6 chondrites are the most highly metamorphosed, possibly in the environment above ∼800 °C, and type 4 and 5 chondrites would be between them (Keil 2000; Scott & Krot 2005; Wlotzka 2005; Huss et al. 2006). Since the experienced temperatures depend on the position inside planetesimals, previous studies make use of thermal evolution models of planetesimals (e.g., McSween et al. 2002; Henke et al. 2013; Gail et al. 2014). Attempts to reproduce the fall statistics of ordinary chondrites (Ghosh & McSween 2000; McSween et al. 2002) have not worked well under the assumption that parent bodies of meteorites are constituted of single-sized planetesimals.

It is, however, highly possible that ordinary chondrites originate from a variety of planetesimals (e.g., Monnereau et al. 2013; Vernazza et al. 2015; Blackburn et al. 2017). Furthermore, the internal temperature of planetesimals is sensitive not only to the location within the planetesimals, but also to their formation times and initial sizes (e.g., Gail et al. 2014; Wakita et al. 2014; Lichtenberg et al. 2016; Ricard et al. 2017). In this paper, therefore, we extensively investigate the thermal evolution of planetesimals with a wide range of model parameters such as formation times and sizes of planetesimals. By comparing our numerical results with the fall statistics of ordinary chondrites in each petrologic type, we will specify the most plausible values of the formation times and size distributions of planetesimals.

The plan of this paper is as follows. In Section 2, we describe how we calculate the thermal evolution of planetesimals with model parameters. Our numerical results are shown in Section 3, where we discuss the dependencies of the model parameters on the results. In Section 4, we compare our results with the fall statistics of ordinary chondrites and discuss implications for planetesimal formation. Our conclusions are given in Section 5.

2.1. Thermal Evolution Models of Planetesimals

We adopt thermal evolution models of planetesimals that were constructed in our previous work (Wakita et al. 2014, 2017b). We assume that planetesimals are spherical and have a spatially uniform chemical composition. It is also assumed that, once planetesimals form with a radius of R_pl at a given time (t = t_pl), they do not experience any further growth and destructive processes afterward. Then, thermal evolution of a planetesimal is numerically calculated by solving the equation of heat conduction,

$$\begin{eqnarray}&&\rho c\displaystyle \frac{\partial T(r,t)}{\partial t}=\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}K\displaystyle \frac{\partial T(r,t)}{\partial r}\right)+A\exp (-\lambda t),\end{eqnarray} \tag{ 1 }$$

where t is the time measured since the formation of the planetesimal, r is the distance from the center of the planetesimal, T(r, t) is the temperature of materials located at r in the planetesimal at a time t, A is the radiogenic heat generation rate per unit volume, and λ is the decay constant of the radionuclides. In our simulations, decay heat arising from the short-lived radionuclide, ²⁶Al is taken as the heating source of planetesimals (e.g., Miyamoto et al. 1982). The initial temperature of planetesimals is set at T(r, t = t_pl) = −123 °C (150 K), which is low enough not to affect the results of our calculations. We assume that physical quantities, such as ρ, c, and K, are constant for simplicity; ρ = 3300 kg m⁻³, c = 910 J kg⁻¹ K⁻¹, and K = 2 J s⁻¹ m⁻¹ K⁻¹ are the adopted values of the bulk density, specific heat, and thermal conductivity, respectively (Yomogida & Matsui 1983; Opeil et al. 2010).

If planetesimals contain volatiles such as water, CO₂, and organics, then different values should be adopted for the above three quantities. Furthermore, if this were the case, the temperature distribution of planetesimals would be regulated not only by thermal conductivity, but also by transport of such volatiles. However, these complexities should be taken into account in exploring thermal evolution of volatile-rich bodies such as parent bodies of carbonaceous chondrites (e.g., Grimm & McSween 1989; Travis & Schubert 2005; Bland & Travis 2017). In this paper, we focus on parent bodies of ordinary chondrites that accreted in very volatile-poor environments, and most of them are (highly) thermally metamorphosed. Therefore, the impact of volatiles can be negligible, except for the ones that contain a small amount of volatiles and/or experience the least metamorphism. In this sense, the physical quantities used here are reasonable, and considering heat conduction is valid for thermal evolution of planetesimals in the context of this paper.

2.2. Key Model Parameters

The main purpose of this study is to carry out a parameter study to determine the fundamental quantities of planetesimals. In this paper, we introduce five parameters (see Table 1): the radius of planetesimals (R_pl), the power-law index (α) for the size distribution of planetesimals, the initial time (t_int) at which planetesimal formation begins, the final time (t_fin) when planetesimal formation ends, and the formation timescale (τ_f) of planetesimals. In the following, we describe what values are considered for these parameters.

Table 1.  Summary of Parameters in This Study

Parameters of Planetesimals	Symbol	Value
Radius (km)	Rpl	1–500
Power-law index for size distribution	α	0.0, 2.8, 3.5, 4.5
Initial time of formation (Myr)	tint	0.1, 1.0, 1.5, 2.0, 2.5, 3.0
Final time of formation (Myr)	tfin	tint, 5.0 or 7.0
Timescale for the formation rate (Myr)	τf	1.0, 4.0a

Note.

^aThe formation rate of planetesimals is assumed to be constant (${\tau }_{f}=\infty$) for all of the simulations, except in Section 3.4.4.

Download table as:  ASCIITypeset image

The typical size of planetesimals and their size distribution are hardly constrained. Accordingly, we take the broad range from R_pl = 1 km to R_pl = 500 km. We assume that the size distribution can be described as a power law in terms of their radius. Then, the number (n_pl) of planetesimals between R_pl and R_pl + dR_pl can be written as

$$\begin{eqnarray}&&{n}_{{pl}}={n}_{0}{({R}_{\mathrm{pl}}/{R}_{\mathrm{pl},0})}^{-\alpha },\end{eqnarray} \tag{ 2 }$$

for ${R}_{\mathrm{pl},\min }\leqslant {R}_{\mathrm{pl}}\leqslant {R}_{\mathrm{pl},\max }$, where n₀ is the normalization constant, and we set R_pl,0 = R_pl,min. For the power-law index (α), we adopt four values: α = 0, 2.8, 3.5, and 4.5. In our fiducial model, we choose α = 2.8, following Johansen et al. (2015) and Simon et al. (2016), who propose that it would represent the primordial distribution of planetesimals.

The formation time of planetesimals is also uncertain since their formation mechanisms are still unclear. It nonetheless might be reasonable to consider that planetesimal formation would take place when the (gaseous) solar nebula was present. Then, we assume that formation time (t_pl) is in the range of 0.1 Myr ≤ t_pl ≤ 7.0 Myr. In this parameter study, we consider that the onset time (t_int) of planetesimal formation is between 0.1 Myr and 3.0 Myr. For the termination time (t_fin) of planetesimals, we consider two cases. For one case, we assume t_pl = t_int = t_fin to focus on the single generation of planetesimals (see Sections 3.1 and 3.2). For the other case, in which multiple generations of planetesimals are examined, t_fin = 5.0 Myr or 7.0 Myr (see Sections 3.3 and 3.4). Note that in the following sections, the time is all measured since the formation of CAIs: CAIs are regarded as the first condensates in the solar system that formed 4567 Myr ago with the initial ratio of ²⁶Al/²⁷Al = 5.25 × 10⁻⁵ (e.g., Connelly et al. 2012). Given that the abundance of ²⁶Al decreases exponentially with a half-life of 0.72 Myr, the earlier formed planetesimals have a higher abundance of the heating source.

The formation rate of planetesimals is another important quantity that can affect global thermal histories of forming planetesimals. In this paper, we simply assume that the formation rate (ϕ) of planetesimals is constant with time for most cases. We will also consider the case that the formation rate of planetesimals decreases exponentially with time. Namely, it is given as ϕ = ϕ₀ exp (−t/τ_f), where ϕ₀ is the hypothetical formation rate of planetesimals in number at t = 0, and τ_f is the formation timescale (see Section 3.4.4).

2.3. Normalization Factors for Mass Fractions

3.1. Thermal Evolution of a Single-sized Planetesimal Formed at a Single Epoch

In this section, we demonstrate the fundamental properties of thermal evolution of a single-sized planetesimal. The results are obtained, using the following set of parameters: R_pl = 100 km, and t_pl = t_int = t_fin = 2.0 Myr or 2.5 Myr, that is, a single generation of planetesimals is considered.

Figure 1 shows the results for the cases of t_pl = 2.0 Myr and t_pl = 2.5 Myr on the left and the right panels, respectively. We first discuss the case of t_pl = 2.0 Myr. The top-left panel depicts temperature evolutions at the positions of r = 0, 80, and 90 km. Initially, the temperature increases with time, as a result of the continuous energy input from radionuclides ²⁶Al. The central region in the planetesimal (r = 0 km) achieves the highest temperature, which exceeds 1000 °C at t = 5.5 Myr. On the other hand, the temperature at r = 80 km (r = 90 km) reaches 920 °C (680 °C) at t = 5.5 Myr (4.2 Myr). And then it starts to decrease because the heating rate due to decay of ²⁶Al becomes lower than the cooling rate at the position. This means that each position of the planetesimal has a peak temperature (T_peak(r)) at t = t_peak(r). The bottom-left panel of Figure 1 shows the temperature evolution at all the positions of the planetesimal. We can see that the temperatures at r ≤ 65 km follow the same evolution and can increase to more than 1000 °C.

Zoom In Zoom Out Reset image size

For the case of t_pl = 2.5 Myr (see right panels of Figure 1), T_peak(r) at the center is 600 °C, demonstrating that the planetesimal formed later has the lower T_peak at the same position than the earlier formed one. This indicates that the thermal evolution is significantly affected by the formation time of planetesimals.

3.2. The Maximum Temperature and the Corresponding Mass Fraction

As discussed above, the peak temperatures inside planetesimals that are heated by radioactive nuclei are different at different locations. For the analysis that will be developed for understanding the origin of the petrologic type of ordinary chondrites, we manipulate the computed values of the temperature (T(r, t)) as follows.

In this work, we are interested in the maximum temperatures that planetesimals experience through thermal evolution. Given that the peak temperatures (T_peak(r)) at position r are achieved at different time (t = t_peak(r)), the maximum temperature (T_max (r, t)) is described as

$$\begin{eqnarray}&&{T}_{\max }(r,t)=\left\{\begin{array}{l}T(r,t)\mathrm{at}\,t\,\lt \,{t}_{\mathrm{peak}}(r)\\ {T}_{\mathrm{peak}}(r)\mathrm{at}\ t\,\geqslant \,{t}_{\mathrm{peak}}(r),\end{array}\right.\end{eqnarray} \tag{ 6 }$$

from the temperature evolution as given in Figure 1. The time evolutions of T_max(r, t) are shown in the top and middle panels of Figure 2. As is obvious from Equation (6), T_max becomes constant after t = t_peak = 5.5 Myr (4.2 Myr) at r = 80 km (r = 90 km) for t_pl = 2.0 Myr.

Zoom In Zoom Out Reset image size

We then introduce the mass fraction (f(t, T_th)) of planetesimals that is determined by a threshold temperature (T_th). When the single generation of planetesimals with the single size is considered, a fraction of mass that experiences the temperatures equal to and higher than T_th can be computed as

$$\begin{eqnarray}&&f(t,{T}_{\mathrm{th}})=\displaystyle \frac{{N}_{1}}{{M}_{1}}\displaystyle \frac{4\pi \rho }{3}{r}^{3}({T}_{\max }(t)\geqslant {T}_{\mathrm{th}}),\end{eqnarray} \tag{ 7 }$$

where t is explicitly labeled both in f and T_max for clearly representing that both quantities are functions of time.

The bottom panels of Figure 2 show the time evolutions of f(t, T_th) for different T_th. Our results show that f(t, T_th) with higher T_th rises up as time goes on. This is a direct reflection of the fact that the temperature at each position increases with time (see the middle panels of Figure 2). Then, any values of f(t, T_th) become almost constant after t ≃ 6 Myr for t_pl = 2.0 Myr. This is because, at this time, almost all positions already pass through T_peak, and T_max no longer changes. Also, our results clearly show that a mass fraction of the planetesimals that undergo a high temperature is reduced for ones formed at later times (see the bottom-right panel of Figure 2): f(t, T_th = 600°C) = 0.3 and f(t, T_th = 800°C) = 0 at t ≳ 7 Myr for the planetesimal formed at t_pl = 2.5 Myr. On the contrary, the planetesimal formed at t_pl = 2.0 Myr has f(t, T_th = 600°C) = 0.76 and f(t, T_th = 800°C) = 0.65 at t ≳ 7 Myr.

This clearly demonstrates that a mass fraction of planetesimals that undergo the temperature above 600 °C strongly depends on the formation time of planetesimals.

3.3. Multiple Generations of Planetesimals with a Single Size

We have so far considered a single generation of planetesimals, that is, t_pl = t_int = t_fin. In general, however, planetesimals are expected to form over some period. Here, we present the results from multiple generations of planetesimals (t_int ≤ t_pl ≤ t_fin) with single sizes. As an example, we adopt the following set of model parameters: R_pl = 100 km, t_int = 0.1 Myr, and t_fin = 7.0 Myr with a constant formation rate.

We first see the effect of planetesimal formation time on f(t, T_th). Figure 3 plots f(t, T_th) at t = 7.0 Myr ≡ t_ms as a function of t_pl(t_int ≤ t ≤ t_fin). We confirm that the mass fraction becomes constant at t_ms ≥ t_fin. We find that as the formation time shifts later, f(t_ms, T_th) becomes lower. For instance, the planetesimals formed after t_pl = 2.3 Myr cannot reach 800 °C. Also, f(t = t_ms, T_th = 0 °C) = 0 at t_pl > 4.3 Myr, indicating that the planetesimals formed later than t_pl = 4.3 Myr cannot be heated above 0 °C even at the center.

Zoom In Zoom Out Reset image size

We then examine the behavior of time-integrated mass fractions of planetesimals formed at multiple epochs. The fraction of mass F(t, T_th) that experienced T_max ≥ T_th for all of the formed planetesimals is given as

$$\begin{eqnarray}&&F(t,{T}_{\mathrm{th}})=\displaystyle \frac{1}{{M}_{2}(t)}\displaystyle \sum _{{t}^{{\prime} }={t}_{\mathrm{int}}}^{t}\phi ({t}^{{\prime} })\displaystyle \frac{4\pi \rho }{3}{r}^{3}({T}_{\max }({t}^{{\prime} })\geqslant {T}_{\mathrm{th}}).\end{eqnarray} \tag{ 8 }$$

Figure 4 shows the time evolution of F(t, T_th). It can be seen that F(t, T_th) for T_th ≥ 800 °C increases until t ≃ 1.5 Myr. This stems from the fact that it takes about 1 Myr for earlier formed planetesimals to reach T_peak(r). On the other hand, F(t, T_th) for T_th ≥ 800 °C starts to decreases after t = 1.5 Myr because a small, or no, part of later formed planetesimals can reach temperatures higher than 800 °C (also see Figure 3). Our results also show some spiky structures at the early times. These are numerical noises, but do not affect our conclusions that are derived from the later stage of the evolutions.

Zoom In Zoom Out Reset image size

In this study, we examine the abundance of ordinary chondrites moderately metamorphosed in thermally evolving planetesimals. In order to perform direct comparison among our results, we define the ratio of F(t, T_th) as follows,

$$\begin{eqnarray}&&{F}_{\mathrm{ratio}}(t,{T}_{\mathrm{th}})\equiv \displaystyle \frac{F(t,{T}_{\mathrm{th}}-100\,^\circ {\rm{C}})-F(t,{T}_{\mathrm{th}})}{F(t,{T}_{\mathrm{th}}=0\,^\circ {\rm{C}})-F(t,{T}_{\mathrm{th}}=1000\,^\circ {\rm{C}})},\end{eqnarray} \tag{ 9 }$$

where T_th varies from 100 °C to 1000 °C.

Figure 5 shows the time evolutions of F_ratio(t, T_th) for different T_th. We find that F_ratio(t, T_th) becomes almost constant after t = 5.0 Myr. Planetesimals formed later than 5.0 Myr cannot contribute to the change of F_ratio(t, T_th), because they do not experience the temperature higher than 100 °C. Accordingly, we focus on planetesimals that form earlier than 5.0 Myr in the following sections, and denote F_ratio(t = t_ms = 5.0 Myr, T_th) as F_ratio(T_th) for brevity.

Zoom In Zoom Out Reset image size

We are now in a position to examine the effect of R_pl on F_ratio(T_th), which is shown in Figure 6. In this figure, we adopt the following parameters: R_pl = 5–500 km, t_int = 0.1 Myr, t_fin = 5.0 Myr, and the constant formation rate is considered. We find that F_ratio(T_th) for larger planetesimals with R_pl ≥ 100 km are very similar, and F_ratio(T_th ≥ 600 °C) is 0.29. This means that even if we consider the arbitrary size distribution for planetesimals larger than 100 km, the results do not change very much. On the other hand, F_ratio(T_th) of planetesimals with R_pl ≤ 50 km are different. Planetesimals of R_pl = 50 km have F_ratio(T_th ≥ 600 °C) of 0.27. For R_pl = 10 and 20 km, F_ratio(T_th ≥ 600 °C) are 0.21 and 0.25, respectively. This indicates that smaller planetesimals produce lower F_ratio(T_th ≥ 600 °C) than larger ones, reflecting the fact that they cannot reach higher temperatures than the larger ones due to more effective cooling from their surfaces.

Zoom In Zoom Out Reset image size

It should be stressed that the radius of the planetesimal (R_pl) plays an important role in thermal evolution of planetesimals and hence the resulting ratio (F_ratio(T_th)). Also, the results are not affected by planetesimals that form later than 5.0 Myr.

3.4. Multiple Generations of Planetesimals with the Size Distribution

In this section, we present the results that are obtained from the most realistic case where planetesimals form within the time interval between t = t_int and t = t_fin and they have the size distribution from R_pl,min to R_pl,max.

Under this situation, the mass fraction of planetesimals that experience a threshold temperature (T_th) or higher is written as

$$\begin{eqnarray}&&{\mathscr{F}}(t,{T}_{\mathrm{th}})=\displaystyle \frac{1}{{M}_{3}(t)}\displaystyle \sum _{{t}^{{\prime} }={t}_{\mathrm{int}}}^{t}\phi ({t}^{{\prime} })\\ &&\times \,{\displaystyle \int }_{{R}_{\mathrm{pl},\min }}^{{R}_{\mathrm{pl},\max }}{{dR}}_{\mathrm{pl}}{n}_{\mathrm{pl}}({R}_{\mathrm{pl}})\displaystyle \frac{4\pi \rho }{3}{r}^{3}({T}_{\max }({t}^{{\prime} })\gt {T}_{\mathrm{th}},{R}_{\mathrm{pl}}).\end{eqnarray} \tag{ 10 }$$

Note that r is now a function of R_pl. Then, the ratio of ${\mathscr{F}}(t,{T}_{\mathrm{th}})$ can be given as

$$\begin{eqnarray}&&{{\mathscr{F}}}_{\mathrm{ratio}}(t,{T}_{\mathrm{th}})\equiv \displaystyle \frac{{\mathscr{F}}(t,{T}_{\mathrm{th}}-100\,^\circ {\rm{C}})-{\mathscr{F}}(t,{T}_{\mathrm{th}})}{{\mathscr{F}}(t,{T}_{\mathrm{th}}=0\,^\circ {\rm{C}})-{\mathscr{F}}(t,{T}_{\mathrm{th}}=1000\,^\circ {\rm{C}})},\end{eqnarray} \tag{ 11 }$$

for 100 °C ≤ T_th ≤ 1000 °C. Armed with this formulation, we carry out a parameter study and examine how ${{\mathscr{F}}}_{\mathrm{ratio}}(t,{T}_{\mathrm{th}})$ at t = t_ms is influenced by model parameters of planetesimals. For brevity, ${{\mathscr{F}}}_{\mathrm{ratio}}(t={t}_{\mathrm{ms}}$ = 5.0 Myr, T_th) is denoted as ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}})$. Unless otherwise mentioned, we take the power-law index α = 2.8 and consider that planetesimals form from t_int = 0.1 Myr up to t_fin = 5.0 Myr with a constant formation rate.

3.4.1. Dependence on Size Ranges of Planetesimals

Planetesimals are likely to form with some size range (size distribution) but not with a single size. Therefore, we consider four size ranges; small (1–50 km in radius), middle (60–100 km in radius), large (100–500 km in radius), and all (1–500 km in radius).

Figure 7 shows the resulting values of ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}})$ for the four size ranges. We find that the results of the middle and large size ranges ((c) and (d) in Figure 7) look similar. Both of them have ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\geqslant 600\,^\circ {\rm{C}})$ = 0.29. On the other hand, the small size range (Figure 7(b)) takes a smaller value of ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\geqslant 600\,^\circ {\rm{C}})$ = 0.22. As we see in the previous section, only a small portion of the planetesimal can reach higher temperatures in small sized planetesimals. Figure 7(a) shows the ratio of all planetesimals combined from the other three ranges.

Zoom In Zoom Out Reset image size

From this, we can conclude that the cumulative mass fraction ratio of planetesimals that experience certain temperatures depends on their size range. In particular, there is a significant difference between the smaller size range (R_pl ≤ 50 km) and the larger one (R_pl ≥ 60 km).

3.4.2. Dependence on Onset of Planetesimal Formation

One of our aims in this paper is to constrain the formation time of planetesimals. In previous sections, we have considered that planetesimals start to form at t_int = 0.1 Myr. In this subsection, we see the dependence of the results on the formation timing (t_int) of planetesimals.

As shown in Figure 3, the mass fraction with higher temperatures decreases with increasing the formation time. If the formation time is later than 2.0 Myr after CAI formation, there is no chance to reach 1000 °C. Figure 8 shows ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}})$ for the four size ranges with onset times of planetesimal formation, t_int = 0.1, 1.0, 1.5, 2.0, 2.5, and 3.0 Myr. When the onset time of planetesimal formation becomes later, ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\,\geqslant 600\,^\circ {\rm{C}})$ decreases. In particular, ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\geqslant 600\,^\circ {\rm{C}})$ in the small size range (1–50 km) rapidly decreases as t_int goes later (Figure 8(b)). This is because smaller planetesimals cannot reach temperatures higher than 600 °C when they form later (t_pl = 1.5 Myr for R_pl = 10 km and t_pl = 2.5 Myr for R_pl = 50 km).

Zoom In Zoom Out Reset image size

We should emphasize that the onset time of planetesimal formation should be earlier than 2.2 Myr (1.2 Myr) to contain the volume that experiences T_max = 800 °C for R_pl = 50 km (R_pl = 10 km). This is important to compare our results with fall statistics of ordinary chondrites (see Section 4).

3.4.3. Dependence on Size Distributions of Planetesimals

Next we see the effect of the size distribution of planetesimals. So far, we take α = 2.8 as the power-law index of the size distribution. This can be regarded as the initial size distribution of planetesimals (Johansen et al. 2015; Simon et al. 2016). Meanwhile, α = 3.5 is known to be the size distribution to explain the current asteroid belt (Morbidelli et al. 2009). We also consider α = 0 and 4.5 as extreme cases.

The results for different size distributions with α = 0, 3.5, and 4.5 are shown in Figure 9. When we compare α = 2.8 (Figure 8(a)) and 3.5 (Figure 9(b)), we find a moderate difference between them: ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\geqslant 600\,^\circ {\rm{C}})$ at t_int = 2.0 Myr are 0.21 and 0.13 for α = 2.8 and 3.5, respectively. For α = 0.0, which means that the number of any size of planetesimals is the same, ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}})$ is mainly determined by larger planetesimals: that is ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\geqslant 600\,^\circ {\rm{C}})$ = 0.28 at t_int = 2.0 Myr (Figure 9(a)). This is comparable with ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\geqslant 600\,^\circ {\rm{C}})$ = 0.27 at t_int = 2.0 Myr for α = 2.8 and R_pl = 100–500 km (Figure 8(d)). When we take α = 4.5, smaller planetesimals regulate the result and ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\,\geqslant 600\,^\circ {\rm{C}})$ is lower: ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\geqslant 600\,^\circ {\rm{C}})$ at t_int = 2.0 Myr is 0.05 for α = 4.5 at t_int = 2.0 Myr (Figure 9(d)).

Zoom In Zoom Out Reset image size

The value of ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}})$ depends largely on the size distribution in the extreme case of α = 0.0 or 4.5. However, the effect of the size distribution is moderate, as long as α ∼3.0.

3.4.4. Dependence on Formation Rate of Planetesimals

In the previous sections, we have assumed that the formation rate of planetesimals is constant. Here, we consider that the formation rate decreases on the timescale of τ_f (see Section 2.2 and Table 1).

The results for τ_f = 1.0 and 4.0 Myr are shown in Figure 10, where the radius of planetesimals ranges from 1 to 500 km with the power-law index of α = 2.8. When τ_f = 1 Myr, ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\geqslant 600\,^\circ {\rm{C}})$ gets larger than that for the constant rate (see Figures 8(a) and 10(a)). This means that, when the formation rate of planetesimals rapidly decreases with time, the contribution from earlier formed planetesimals, which experience higher temperatures, becomes stronger. For τ_f = 4 Myr, the results are close to those with the constant formation rate (see Figures 8(a) and 10(b)).

Zoom In Zoom Out Reset image size

From these results, we conclude that there is the distinguishable difference between the results for the constant formation rate and the decreasing formation rates with short τ_f.

In this section, we compare our results with the abundance of each petrologic type of ordinary chondrites and discuss the formation and evolution histories of planetesimals.

The fall statistics of ordinary chondrites indicates that the number abundance of each petrologic type is 5.2% for type 3, 16.5% for type 4, 38.5% for type 5, 37.9% for type 6, and 1.9 % for others (Grady et al. 2014). The petrologic types of chondrites are divided by chemical compositions and texture, as well as their peak metamorphic temperature. Thus, the peak metamorphic temperature would not necessarily be adequate to separate each type rigorously; especially because the estimated peak metamorphic temperatures are overlapped between types 4 and 5 (Scott & Krot 2005). Nevertheless, in comparison to our results, we assign the following temperatures to each type: 100–600 °C to type 3, 600–800 °C to types 4 and 5, and 800–1000 °C to type 6 (Scott & Krot 2005; Huss et al. 2006). In our simulations, we have computed the mass fraction of planetesimals (${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}})$) for a variety of initial radii and formation times of planetesimals (see Equation (11)). Although the fall statistics is given by the number of meteorites (by percentages), the statistics may represent the original (mass) abundance in their parent bodies, as a results of (thanks to) the huge number of meteorite samples. Hence, we compare ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}})$ with the fall statistics.

First, we find that, for all models considered in this study, ${{\mathscr{F}}}_{\mathrm{ratio}}(100^\circ {\rm{C}}\leqslant {T}_{\mathrm{th}}\lt 600\,^\circ {\rm{C}})$ is higher than ${{\mathscr{F}}}_{\mathrm{ratio}}(600\,^\circ {\rm{C}}\,\leqslant {T}_{\mathrm{th}}\lt 800\,^\circ {\rm{C}})$ and ${{\mathscr{F}}}_{\mathrm{ratio}}(800\,^\circ {\rm{C}}\leqslant {T}_{\mathrm{th}}\leqslant 1000\,^\circ {\rm{C}})$, meaning that type 3 ordinary chondrites are the most abundant in our models. This result is consistent with previous work (McSween et al. 2002), but it is opposite to the fall statistics of type 3 chondrites, whose mass fraction is much lower than those of types 4 and 5 and type 6 (Grady et al. 2014). We will discuss this point later.

Although our model cannot reproduce the fall statistics of type 3 ordinary chondrites, it is interesting to compare our results with the observed relative abundance of type 6 to types 4 and 5 chondrites. This is because these types of chondrites are dominant (90% of all chondrites). Additionally, when thermal evolution of planetesimals is the main cause to generate types 4, 5, and 6 chondrites in their parent bodies, these chondrites should have been present deep inside of the planetesimals (see Figure. 1). Consequently, the relative abundance of type 6 to types 4 and 5 chondrites could be a key to specifying formation times of planetesimals and their primordial sizes. Figure 11 represents the mass ratios of type 6 to types 4 and 5 chondrites obtained from our results as a function of the onset time of planetesimal formation (see legends). In this plot, we consider some representative cases. Note that ${{\mathscr{F}}}_{\mathrm{ratio}}({T}_{\mathrm{th}}\geqslant 600\,^\circ {\rm{C}})$ becomes zero at t_int > 2.5 Myr, thus we cannot plot the result of t_int ≥ 3.0 Myr. Table 2 summarizes our results for the case of t_int = 1.5 Myr, together with a set of parameters.

Zoom In Zoom Out Reset image size

One of our important findings is that the parent bodies of ordinary chondrites must begin to form before t ∼ 2.0 Myr, to produce type 6 chondrites that undergo thermal metamorphism with the temperature of 800 °C or higher. In other words, if formation of planetesimals would be initiated after t ∼ 2.0 Myr, no planetesimal could be heated high enough to contain type 6 chondrites, which contradicts the observational fact.

It should be emphasized that the calculated mass ratios of type 6 to types 4 and 5 are in agreement with the fall statistics within the uncertainty of 30% for most of our models when t_int ≤ 2 Myr. The planetesimal formation time of 2.0 Myr is consistent with the estimation from our previous work for the parent bodies of the asteroid Itokawa (Wakita et al. 2014). However, adopting t_int = 2.0 Myr, the mass ratios of type 6 to types 4 and 5 from the small size range (R_pl = 1–50 km) of planetesimals become lower than the fall statistics (see the legend Fig. 8(b) in Figure 11 and Table 2). The result with α = 4.5, wherein small planetesimals dominate the total mass, also underproduces the fall statistics (see the legend Fig. 9(c) in Figure 11 and Table 2). On the other hand, when planetesimal formation begins at t ≤ 1.5 Myr, these cases can successfully explain the fall statistics of type 6 ordinary chondrites. Therefore, our results suggest two scenarios for the formation condition of planetesimals; (1) planetesimals should have started to form at t ≤ 2.0 Myr after CAIs when their primordial sizes were larger than 60 km and/or (2) planetesimals whose size distribution was weighted toward small radii of 1–50 km should have started to form at t  ≤ 1.5 Myr after CAI formation. These conclusions are schematically depicted in Figure 12.

Zoom In Zoom Out Reset image size

Figure 11 also shows that the results for size distributions with α = 0 and 3.5 lead to a slightly worse agreement with the fall statistics than that for the fiducial case with α = 2.8 (see the legends Fig. 8(a), 9(a), and 9(b) in Figure 11 and Table 2). However, the difference is not so large that we cannot impose any constraint on the primordial size distribution of planetesimals. Meanwhile, our results can provide some implications to the formation rate of planetesimals (see the legends Fig. 8(a), 10(a), and 10(b) in Figure 11 and Table 2). When the formation rate of planetesimals decreases slowly (τ_f = 4.0 Myr), the mass ratios of type 6 to types 4 and 5 can reproduce the fall statistics relatively well (see the legend Fig. 10(b) in Figure 11 and Table 2). Even for the rapid (τ_f = 1.0 Myr) case whose result is slightly higher than the fall statistics, it is in a reasonable range (see the legend Fig. 10(a) in Figure 11 and Table 2). Although the preferable formation rate of planetesimals can be a constant one or a relatively long formation timescale (τ_f = 4.0 Myr), we cannot rule out the rapid formation case (τ_f = 1.0 Myr).

Our results show that the mass ratios of type 6 to types 4 and 5 are produced by the thermal evolution model of planetesimals. However, the mass fraction of type 3 ordinary chondrites cannot be explained. As shown in Figures 8–10, even if planetesimals start to form as early as 0.1 Myr, about 70% of the total mass of the planetesimals cannot experience temperatures higher than 600 ^◦C. Therefore, to lower the abundance of type 3 chondrites than those of types 4–6 ones, planetesimals have to be heated more than the present calculations. However, we do not know other internal sources and effective mechanisms of heating the whole planetesimals (except for the decay of short-lived radionuclides).

Type 3 chondrites are divided into petrologic subtypes from 3.0 to 3.9 based on their characteristics of thermoluminescence sensitivity and chemical compositions of minerals, such as iron content of olivines (Sears et al. 1980; Scott & Jones 1990; Scott et al. 1994; Grossman & Brearley 2005; Bonal et al. 2006). Although our model cannot reproduce the abundance of type 3 chondrites relative to types 4 and 5 and type 6, it is still interesting to see what implications our results have for the fall statistics of type 3 subtypes. It is suggested that types 3.0–3.2 chondrites experienced much lower peak temperatures than the other type 3s and that type 3.9 would have experienced almost the same temperature as type 4 (e.g., Huss et al. 2006; Cody et al. 2008; Vernazza et al. 2014). It is very hard to distinguish the subtypes only by their peak metamorphic temperatures, but we here roughly divide type 3 subtypes into three subgroups, assigning peak metamorphic temperatures as follows: 100–400 °C for types 3.0–3.2, 400–500 °C for types 3.3–3.6, and 500–600 °C for types 3.7–3.9. The relative fall statistics of type 3 subtypes are derived based on the Meteoritical Bulletin Database (http://www.lpi.usra.edu/meteor/metbull.php): 0.132, 0.398, and, 0.470 for types 3.0–3.2, 3.3–3.6, and 3.7–3.9, respectively. Note that we use only the samples of type 3 ordinary chondrites that are clearly classified as subtypes. We summarize the normalized fall statistics of type 3 subtypes and our results for the case of t_int = 1.5 Myr in Table 3. Our results show that the lower numbers of type 3 subtypes are more abundant than the higher numbers of them (i.e., types 3.0–3.2 > 3.3–3.6 > 3.7–3.9), which is opposite to the trend that appeared for the fall statistics (types 3.0–3.2 < 3.3–3.6 < 3.7–3.9). In particular, our models predict much higher abundances of types 3.0–3.2 (>0.8) than its fall statistics (0.13, see Table 3). Even if we take the ratio of types 3.3–3.6 to types 3.7–3.9, almost all of our results overestimate the observed values out of the uncertainty of 30% (Figure 13). This disagreement might partially stem from the facts that the fall statistics of type 3 subtypes used here are far from complete (only 1/3 of type 3 ordinary chondrites is sorted out into subtypes) and/or that we do not consider the heat transport through volatile materials, which would affect the abundances of the least metamorphosed ones (types 3.0–3.2).

Zoom In Zoom Out Reset image size

One strong possibility to account for the small fall statistics of type 3 ordinary chondrites (as well as the tendency that less metamorphosed type 3 subtypes are less abundant) is impact events for planetesimals. In this study, we assume that planetesimals do not undergo any destructive processes once they formed. However, planetesimals that were not incorporated into (proto)planets would have suffered from destructive collisions with their surrounding planetesimals over the age of the solar system. In fact, destructive collisions are one of the important ingredients for reproducing the size distribution of the current main asteroid belt (Morbidelli et al. 2009). When planetesimals experience destructive collisions, materials located in their surface regions are ejected. Such materials can potentially serve as the parent bodies of type 3 chondrites due to a low level of metamorphism. It nonetheless can be expected that they would not contribute to the (current) fall statistics very much. This is because they would probably be accreted onto nearby (proto)planets and/or be gone beyond the solar system on a relatively short timescale. Furthermore, the removal of surface materials via destructive collisions would increase the chance for the core region of planetesimals to be exposed to subsequent collisions. Given that materials in the core region can become the parent bodies of types 4–6 chondrites, destructive collisions can eventually enhance the possibility that types 4–6 chondrites fall onto the Earth as meteorites. It is interesting that the asteroid Itokawa would be mainly composed of types 5–6 chondrites, and its parent body is considered to have been originally larger than 20 km in radius before being catastrophically destroyed by a big impact (e.g., Nakamura et al. 2011; Wakita et al. 2014). Thus, the destructive collisions that planetesimals have underwent during their long life may effectively decrease the amount of type 3 chondrites. In other words, the fall statistics of type 3 chondrites may hold the key information on the collision histories of planetesimals.

In this paper, we systematically examine the thermal evolution of planetesimals with a wide range of initial sizes and formation times. Our numerical results show that planetesimals can reach higher temperatures when they form at earlier times and/or have larger initial radii. There is a considerable difference in maximum temperatures between a small size range (from 1 to 50 km) and a larger size range (from 60 to 90 km and/or from 100 to 500 km) (see Figure 8).

We also examine the effect of the formation rate of planetesimals. We found that, in the case of the continuous formation of planetesimals, their mass fraction is significantly affected by the onset time of planetesimal formation (Figure 8). This trend is prominent when we focus on the planetesimal mass that experiences higher temperatures (≥600 °C). Furthermore, our results show that planetesimals should form within 2.0 Myr after formation of Ca-Al-rich inclusions (CAIs), in order to produce type 6 ordinary chondrites, which would have experienced the peak metamorphic temperature of 800 ∼ 1000 °C. We also found that the ratio of type 6 to types 4 and 5 chondrites in fall statistics can be well explained by thermal evolution models of planetesimals. This strongly suggests that these types of chondrites are produced mainly by thermal metamorphism inside the parent bodies. Preferable scenarios to account for the fall statistics are that planetesimals with radii of ≥60 km start to form around 2.0 Myr and/or that planetesimals of ≤50 km start to form at ≤1.5 Myr after the CAIs formation.

On the other hand, our thermal evolution models of planetesimals cannot explain the fall statistics of type 3 ordinary chondrites; we predict that type 3 chondrites are unavoidably more abundant than any of the types 4 to 6 chondrites. This is against the fact that type 3 chondrites are the least abundant in the fall statistics. The overabundance of type 3 chondrites could be resolved by taking into account subsequent destructive processes such as collisions between planetesimals. In our future work, we will explore the effect of impact and subsequent destruction for a better understanding of the formation and evolution histories of planetesimals.

We thank the referee, Andrew M. Davis, for helpful comments and suggestions. Numerical computations were carried out on the PC cluster at the Center for Computational Astrophysics, National Astronomical Observatory of Japan. T.N. has been supported in part by a JSPS Grant-in-Aid for Scientific Research (26400223). Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. Y.H. is supported by JPL/Caltech.