THE ONSET OF DIFFERENTIATION AND INTERNAL EVOLUTION: THE CASE OF 21 LUTETIA

Asteroid 21 Lutetia plays an important role in our understanding of the origin and evolution of planetary objects. As noted by Bottke et al. (2005), the size of Lutetia allows it to survive disruption from impact. Lutetia accordingly preserves its original large-scale structure. Data provided by the Rosetta spacecraft suggest a high bulk density for Lutetia (3400 ± 300 kg m⁻³; Patzold et al. 2011). This fact, combined with the primitive nature of its crust (carbonaceous or enstatite chondrites; Coradini et al. 2011), is consistent with a scenario in which 21 Lutetia experienced partial differentiation with the formation of a metallic "core" surrounded by a primitive chondritic crust (Weiss et al. 2011).

The geophysical and thermophysical history of Lutetia depends strongly on its initial composition and global macroporosity; these values are greatly uncertain. Considerations of the global structure lead some authors to classify Lutetia as an asteroid with abundant fractures and joints (Asphaug 2009) for which the inferred macroporosity is in the range of ∼6%–40% (Consolmagno et al. 2008; Wilkinson et al. 2002). Furthermore, all other asteroids with similar sizes (except 20 Massalia) are thought to have macroporosities of >5%–10% and ranging up to ∼80% (Consolmagno et al. 2008). A strict upper limit on Lutetia's macroporosity (i.e., ∼52%) is provided by a model assuming that the entire asteroid is below a very thin chondritic surface layer and is made up of pure iron.

Weiss et al. (2011) proposed three different scenarios for producing the high bulk density of Lutetia via partial differentiation. In the first scenario, primordial Lutetia has the same size as the present asteroid and is initially undifferentiated. After radiogenic heating and subsequent internal melting, the metallic core and silicate mantle form with a decreased macroporosity and an increased bulk density. In the second scenario, primordial Lutetia has a larger radius and only a smaller fractional volume experiences melting. The undifferentiated outer layer is removed by subsequent impacts while the bulk density increases. In the third scenario, the differentiation of Lutetia occurs early, and at first the chondritic crust is not retained or it is removed by subsequent impacts. After this phase, chondritic debris is deposited on the asteroid.

In order to investigate which of these scenarios is the most plausible, we studied the thermal history of Lutetia. Before describing our model, however, we review the literature, as a wide range of thermal models of planetesimals assuming ²⁶Al as the main heat source has been developed in recent years.

Merk et al. (2002) analyzed the dependence of the accretion rate on the thermal evolution of the planetesimal, concluding that the accretion process had to be considered as long as the accretion time was not negligible with respect to the half-life of the radionuclides providing the energy. The authors also made use of the Stefan law formulation in order to incorporate the effect of latent heat into thermal evolution, but they neglected the role of sintering and convection. Furthermore, they did not use a radiation boundary condition (i.e., they assumed a fixed surface temperature).

Ghosh et al. (2003) also focused their attention on the importance of the accretion process in the thermal evolution of asteroids by studying the case of 6 Hebe. These authors investigated the differences between instantaneous and incremental accretion cases. They used the model developed by Ghosh & McSween (1998) with a moving boundary condition and a radiation boundary condition. The authors concluded that incremental models provide a link between theoretical models of measurable quantities in meteorites (i.e., peak temperatures, cooling rate, radioisotope closure times) that were determined by their thermal histories.

Hevey & Sanders (2006) incorporated convection in their thermal evolution models when the degree of partial melting exceeded 50% by volume. They analyzed the effects of sintering, whose onset was set at about 700 K, starting with a high porosity (50%) and a low thermal conductivity. They considered instantaneous accretion and, like Merk et al. (2002), used a fixed boundary temperature. The main result of their work was to constrain the accretion of the parent bodies of differentiated meteorites to within the first 1.5 Myr, or at most 2 Myr, from Ca–Al-rich inclusions (CAIs). These planetesimals therefore formed before most chondritic parent bodies, which accreted at a later time.

Moskovitz & Gaidos (2011) studied the thermal consequences of melt migration and, in particular, investigated how the redistribution of ²⁶Al from the interior into a crustal layer would affect the thermal evolution of planetesimals. They considered the case of instantaneous accretion excluding sintering, convection, and a radiation boundary condition, and they concluded that differentiation would be most likely for planetesimals larger than 20 km in diameter and those accreted within approximately 2.7 Myr from CAIs.

Sahijpal & Gupta (2011) performed numerical simulations of the processes involving both aqueous alteration and planetary differentiation. They used asteroids of 100 and 270 km as representative case studies, starting with high porosity and low thermal conductivity. They inserted the sintering and simulated the convection (at 50% of melting of silicates) by raising the thermal conductivity by three orders of magnitude. Their model assumes a fixed boundary temperature and does not take into account Darcy's law formulation for the segregation of silicates and melt. They concluded that a convective molten iron core is necessary to explain the remanent magnetism of the carbonaceous chondrite.

Finally, Elkins-Tanton et al. (2011) investigated the possibility that early radiogenic heating of planetesimals could create partially differentiated bodies with a primitive crust and magnetic core dynamos. The undifferentiated crust must be thick enough to prevent the majority of impacts from breaching it and persist through the internal magma ocean stage. They concluded that planetesimals accreting before about 1.5 Myr after CAIs are likely to fully differentiate through radiogenic heating, while bodies that accrete past about 1.5 Myr from CAIs would probably be characterized by a thick undifferentiated crust overlying a differentiated interior.

In the present work, we apply the thermal code we developed to study the thermal history of Vesta (Formisano et al. 2013) in order to constrain the formation time, size, and mass of the proto-core of Lutetia by varying the time delay in the injection of ²⁶Al in the solar nebula (a similar approach is also presented in Hevey & Sanders 2006, Moskovitz & Gaidos 2011, and Elkins-Tanton et al. 2011) and the post-sintering macroporosity. We define as proto-core a structure enriched in metals, formed as a consequence of the metal percolation, but still containing pristine silicates.

As in Hevey & Sanders (2006) and Elkins-Tanton et al. (2011), we assume instantaneous accretion and solve the heat equation with radiogenic heating provided by ²⁶Al and ⁶⁰Fe. As shown by Weidenschilling (2008), the formation of a planetesimal having the same size as Lutetia requires about 10⁵ yr or less. Moreover, as pointed out by Moskovitz & Gaidos (2011), recent dynamical studies treating the turbulent concentration of small particles in protoplanetary disks (Johansen et al. 2007, 2009; Cuzzi et al. 2010) show that planetesimals can grow "nearly instantaneously" in less than 100 yr to sizes of 100 km or larger. As a consequence, instantaneous accretion can be a reasonable approximation.

Similar to Moskovitz & Gaidos (2011), we do not analyze the sintering phase due to the large uncertainties associated with the assumed initial porosities and the corresponding thermal conductivities.

We include the treatment of the latent heat based on the Stefan law formulation, as in Merk et al. (2002), in order to incorporate its effect on the evolution, and we use Darcy's law model for metal percolation.

The developed scenarios will be useful in depicting a reliable geophysical and thermal history and constraining the internal structure of Lutetia. They also offer a complementary approach to the works present in the literature, like evidence for differentiation provided by gravitational features (Vincent et al. 2011) or by hydrocode modeling of the largest impact crater on Lutetia (Oklay et al. 2012).

2.1. Initial and Boundary Conditions

We consider 21 Lutetia as a spherical body with a radius fixed at 50 km with an initially homogeneous composition. We suppose that the object contains a metallic component (Y) (about 25% by volume) and a silicatic one (X) (about 75% by volume). This composition is similar to those of the H and L classes of ordinary chondrites. The post-sintering porosity ranges from 10%–30% by volume. The initial temperature (T₀), which is also the temperature of the solar nebula surrounding Lutetia, is fixed at 200 K (Lewis 1974). In Table 1, we report the physical parameters used in our model. We use a radiation boundary condition at the surface and a Neumann boundary condition at the center (heat flux equal to zero):

$$\begin{eqnarray} T(r,t = 0) &=& T_0, \nonumber\\ \left[ {\frac{{\partial T}}{{\partial r}}} \right]_{r = 0} &=& 0, \\ \left[ \frac{\partial T}{\partial r}\right]_{{\rm surf}} &=& - \frac{\varepsilon \sigma }{K} \left(T_{{\rm surf}}^4 - T_{0}^4\right), \nonumber \end{eqnarray} \tag{ 1 }$$

where T_surf is the temperature of the surface, ε is the emissivity, K is the thermal conductivity, and σ is the Stefan–Boltzmann constant (see Table 1).

Table 1. Physical Parameter Values Used in this Work

Quantity	Value	Unit	Reference
Final primordial radius (R)	50 × 103	m	Weiss et al. (2011)
Density of metal (ρmet)	7800	kg m−3	Sramek et al. (2012)
Density of silicate (ρsil)	3200	kg m−3	Sramek et al. (2012)
Specific heat of metal (solid) (cmet, sol)	600	J kg−1 K−1	Sahijpal et al. (2007)
Specific heat of metal (liquid) (cmet, liq)	2000	J kg−1 K−1	Sahijpal et al. (2007)
Specific heat of silicate (solid) (csil, sol)	720	J kg−1 K−1	Sahijpal et al. (2007)
Specific heat of silicate (liquid) (csil, liq)	720	J kg−1 K−1	Sahijpal et al. (2007)
Latent heat of metal (Lmet)	270	KJ kg−1	Ghosh & McSween (1998)
Latent heat of silicate (Lsil)	400	KJ kg−1	Ghosh & McSween (1998)
Metal solidus ($T_{{\rm sol}}^{{\rm met}})$	1213	K	Ghosh & McSween (1998)
Metal liquidus ($T_{{\rm liq}}^{{\rm met}})$	1233	K	Ghosh & McSween (1998)
Silicate solidus ($T_{{\rm sol}}^{{\rm sil}}$)	1425	K	Taylor (1992)
Silicate liquidus ($T_{{\rm liq}}^{{\rm sil}}$)	1850	K	Taylor (1992)
Thermal conductivity of metal (Kmet)	50	W m−1 K−1	Sramek et al. (2012)
Thermal conductivity of silicate (Ksil)	3	W m−1 K−1	Sramek et al. (2012)
Initial metal volume fraction (Y)	25%
Initial silicate volume fraction (X)	75%
Post-sintering porosity (ϕ)	10%–30%
Temperature of solar nebula (T0)	200	K	Lewis (1974)
Stefan–Boltzmann constant (σ)	5.67 × 10−8	W m−2 K−4
Emissivity (ε)	1
Half-life of 26Al (τAl)	0.717	Myr	Castillo-Rogez et al. (2009)
Specific heat production of 26Al	0.355	W kg−1	Castillo-Rogez et al. (2009)
Initial isotopic abundance of 26Al in ordinary chondrites ([26Al]0)	6.20 × 10−7	ppb	Castillo-Rogez et al. (2009)
Half-life of 60Fe (τFe)	2.62	Myr	Rugel et al. (2009)
Specific heat production of 60Fe	0.068–0.074	W kg−1	Castillo-Rogez et al. (2007)
Initial isotopic abundance of 60Fe in ordinary chondrites ([60Fe]0)	(22.5–225) × 10−9	ppb	Castillo-Rogez et al. (2007)

Download table as:  ASCIITypeset image

2.2. Physical and Numerical Description

We numerically study the thermal evolution of 21 Lutetia, solving the following coupled equations:

$$\begin{equation} {\left({\rho c} \right)_m}\frac{{\partial T}}{{\partial t}} = \vec{\nabla }\cdot ({{K_m}\vec{\nabla }T} ) + H, \end{equation} \tag{ 2 }$$

$$\begin{equation} \frac{{\partial Y}}{{\partial t}} + \vec{v} \cdot \vec{\nabla }Y = 0. \end{equation} \tag{ 3 }$$

Equation (2) is the heat equation with the radiogenic heat source (H) (Castillo-Rogez et al. 2007), where (ρc)_m and K_m are the overall heat capacity and the overall thermal conductivity, respectively. Equation (3) is the advection equation which controls the metal percolation in the silicatic matrix, with the migration velocity being a function of the permeability of the silicatic medium, the density contrast between molten metals and solid silicates, and the gravitational acceleration and of the viscosity of the molten iron (Yoshino et al. 2004). As in Ghosh & McSween (1998), when the melting temperature of Fe–FeS is reached (1213 K), the percolation of the metals takes place through the silicatic matrix and the formation of the proto-core occurs. To take into account in a simple way the latent heat during the phase transition, the specific heat is modified through the Stefan coefficient:

$$\begin{equation} {\rm Ste} = \frac{L}{c}\frac{{d\chi }}{{dT}} = \frac{L}{c}\frac{1}{{T_{{\rm liq}} - T_{{\rm sol}}}}, \end{equation} \tag{ 4 }$$

where χ is the degree of melting:

$$\begin{equation} \chi = \frac{{T - T_{{\rm sol}}}}{{T_{{\rm liq}} - T_{{\rm sol}}}}, \end{equation} \tag{ 5 }$$

and so

$$\begin{equation} \bar{c} = c(1 + {\rm Ste}), \end{equation} \tag{ 6 }$$

assuming a linear growth of χ with increasing temperature.

The numerical method is based on one dimension, forward time-central space (FTCS), and a finite difference scheme in the radial direction with a Lax correction. A spatial grid (Δr) of 500 m is used. To ensure the stability of our numerical approach, we use an adaptive time increment according to the Courant–Friedrichs–Lewy stability conditions for each of the physical processes we consider in our work (heat diffusion, metal percolation, and the radiation boundary condition). Following Toksoz & Solomon (1973), thermal conduction imposes the following critical time step:

$$\begin{equation} \Delta t_{{\rm cond}} = \frac{\left(\rho c\right)_m \Delta r^2}{2K_m}. \end{equation} \tag{ 7 }$$

Analogous to Toksoz & Solomon (1973), we define the following critical time step associated with the radiation boundary condition:

$$\begin{equation} \Delta t_{{\rm rad}} = \frac{\left(\rho c\right)_m \Delta r T_{{\rm surf}}}{\sigma \left(T^4_{{\rm surf}} - T^4_0\right)}. \end{equation} \tag{ 8 }$$

Finally, while the percolation of metals is taking place, we need to solve also Equation (3) and introduce a third critical time step:

$$\begin{equation} \Delta t_{{\rm perc}} = \frac{\Delta r}{v}, \end{equation} \tag{ 9 }$$

where v is the velocity of the metal percolation. The Courant–Friedrichs–Lewy stability condition requires that the time step used in our model satisfies the following criterion:

$$\begin{equation} \Delta t < {\rm min}(\Delta t_{{\rm cond}}, \Delta t_{{\rm rad}}, \Delta t_{{\rm perc}}). \end{equation} \tag{ 10 }$$

Therefore, at each temporal iteration of the program, we select the actual time step based on the minimum critical time step among those we computed. As the stability condition requires the actual time step to be lower than the critical one, we choose to use a value equal to 90% of the smallest critical time step as a compromise between the competing needs for stability and performance. So, our time step is defined as

$$\begin{equation} \Delta t = 0.9 \times {\rm min}(\Delta t_{{\rm cond}}, \Delta t_{{\rm rad}}, \Delta t_{{\rm perc}}). \end{equation} \tag{ 11 }$$

We explored several scenarios characterized by different strength of the energy sources (the radiogenic heat due to the decay of ²⁶Al) and values of post-sintering macroporosity (10%, 20%, and 30% by volume). The scenarios are labeled N0 (instantaneous accretion, Δt_d = 0), N1 (Δt_d ≃ 0.3 Myr), and N2 (Δt_d ≃ 0.7 Myr). The main results are shown in Table 2, in which we report the size, the time of formation, the mass of the proto-core (i.e., a structure enriched in metals and containing pristine silicates), and the maximum temperature reached after 5 Myr. In all the cases we analyzed, the maximum degree of silicate melting (∼10% by volume, corresponding to about 1450 K; Taylor 1992) is reached only in a limited region of Lutetia. Since the degree of silicate melting is very low and, following Hevey & Sanders (2006), the onset of convection should require that the silicate melt fraction exceed 50%, we conclude that heat transport via convection is negligible. In fact, even in those regions where a significant temperature difference is present, the Rayleigh number does not overcome the critical values lying between 1000 and 2000.

In the case of post-sintering porosity of 10% by volume, we observe that in all scenarios Lutetia does not completely differentiate and only a proto-core forms (see the maximum temperature versus time profile of Figure 1(f)). In Figure 1(a), after 0.1 Myr, in the three scenarios the temperature is lower than the solidus temperature of silicates and the asteroid is homogeneously heated. After 0.5 Myr (see Figure 1(b)), N0 enters in the melting temperature of silicates, while in N1 and N2 the temperature are less than 1425 K. After 1 Myr (see Figure 1(c)), in N1 the temperature reaches the solidus temperature of silicates, while in N2 the values are still low. In Figure 1(d), after 3 Myr, we observe a slight general increase of the temperature in all scenarios and then, after 5 Myr (see Figure 1(e)), in N0 the temperature has almost the same value while in N1 and N2 we observe a general decrease. The proto-core size ranges from 25 to 31 km and the time of formation ranges from 1.2 to 2.2 Myr.

Zoom In Zoom Out Reset image size

The general trend, for a porosity of 20% by volume, is similar to the previous case, but the temperature values reached are lower because of the lower amount of material, and therefore energy sources, per unit volume (see Figure 2(a)). As we can observe in Figure 2(b), after 0.5 Myr, the temperatures overcome the liquidus temperature of the metallic component in N0 and N1, while in N2 the temperature is in the window of melting of metals. After 1.5 Myr (see Figure 2(c)), Lutetia is in the heating phase for all the scenarios and, after 3 Myr (see Figure 2(d)), the general trend is the same. In Figure 2(e), after 5 Myr, we observe a general decrease of the temperatures for N1 and N2. In Figure 2(f), the maximum temperature versus time profile is shown. We observe that the maximum temperature is reached in the hottest scenario, characterized by instantaneous accretion (N0), and that the time of formation of the proto-core ranges from 1.7 to 3.2 Myr, while the proto-core size ranges from 14 to 21 km.

Zoom In Zoom Out Reset image size

If we choose a value of porosity of 30% by volume, we can observe that, after an isothermal phase (see Figure 3(a)) for all scenarios, at 0.5 Myr the temperature overcomes the liquidus melting temperature of metals in N0, while in N1 it is in the windows of melting of metals and in N2 it is less than 1213 K (see Figure 3(b)). After 1.5 Myr (see Figure 3(c)), N1 overcomes the liquidus temperature of metals and N2 enters in the windows of melting of metals. The general trend is the same after 3 and 5 Myr (see Figures 3(d) and (e), respectively). In no scenario does the temperature reach the solidus temperature of silicates (see Figure 3(f)). The formation of the proto-core occurs from 2.3 to 4 Myr. Figure 4 summarizes the results obtained for the three cases of porosity explored.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 4(a)–(c), the maximum temperatures are reached in the middle region of the asteroid (from 30 to 40 km from the center) as a consequence of partial differentiation; in this region, in fact, there is less mass to be heated and so the temperatures are higher. In Figure 4(d)–(f), we can see that the general trend is the same as in Figure 4(a)–(c) but the temperatures reach lower values than those in the previous case because the porosity is increased (i.e., 20% by volume). In N0, the melting of silicates is possible, while in N1 and N2 only the melting of metals occurs. Figure 4(g)–(i) shows that high values of porosity (i.e., 30% by volume) prevent the temperature from reaching the melting temperature of silicates. In particular, when the delay in the injection of ²⁶Al is greater (i.e., 0.72 Myr), the melting of metals is possible only in a narrow region of Lutetia, ranging from 5 to 15 km.

Although we would expect an increase in temperature with increasing porosity due to the lowering of the thermal conductivity (Opeil et al. 2010), the dominant effect of the increase in porosity is the decrease of the concentration of the radioactive source (²⁶Al) linked to the silicates. Therefore, the hottest scenarios are those characterized by the lowest porosity, as shown in Table 2.

We also observe that in all the scenarios analyzed, a surface layer ranging from 2 to 5 km is below the Curie temperature of carbonaceous chondrites (corresponding to about 873 K). This observation means that a remanent magnetization of this body is possible (Richter et al. 2012). The remanent magnetization of Lutetia is not for certain—it could be externally generated by sources in the primordial solar system (Weiss et al. 2010) or it could be due to the presence of an internal core dynamo that could confirm the scenario of partial differentiation.

We can make a direct comparison between our temperature profiles and those of other models. Comparing our results with Hevey & Sanders (2006), in the case of a 50 km body and instantaneous or late accretion (0.75 Myr), the maximum temperature reached (about 8000 and 4000 K, respectively) are very different from our values (in the hottest scenario of the current paper, i.e., N0 and 10% by volume porosity, we obtain 1455 K). Our thermal profiles also show a maximum in the "mantle," due to the differentiation and the affinity of ²⁶Al for silicates. These global differences probably depend on the different methodology and initial conditions adopted in the two works. The authors concluded that the thin crust that formed is destined to be processed by the magma flowing beneath due to small impacts or convective drag. On the contrary, Elkins-Tanton et al. (2011) concluded that an undifferentiated chondritic crust survives through the internal magma ocean phase. If we define the formation of the magma ocean at 50% by volume of the silicate melting, the small degree of silicate melting reached in our scenarios prevent this situation and preserve the primitive unmelted crust. Sahijpal & Gupta (2011) also showed the existence of a chondritic crust, but they assumed a fixed boundary temperature instead of a more realistic radiation boundary condition (so the extent of their chondritic crust was likely overestimated). Since the methodology assumed in Moskovitz & Gaidos (2011) is similar to that of our work, their maximum temperature values are more compatible with our own than those of Hevey & Sanders (2006), even if their global profiles are different.

Observational data do not provide stringent constraints about the internal structure. Currently we know that Lutetia possesses a chondritic crust (carbonaceous or enstatitic) and that its high bulk density has been interpreted as an indication of the presence of a metallic core. The results of our model suggest that partial differentiation is possible. In fact, the maximum degree of silicate melting is about 10% by volume in a limited region of the "mantle." This is consistent with the scenarios proposed by Weiss et al. (2011), if the current macroporosity (10%–30% by volume) is the same as the post-sintering one. In all scenarios, only the formation of a proto-core, a structure enriched in metals, occurs. The proto-core formation takes from 1 to about 4 Myr and proto-core size ranges from 6 to 30 km. The relative proto-core mass ranges from about 1 to about 36% of the total mass. Our results suggest that the accretion time of Lutetia should not exceed 0.7 Myr from CAIs and the post-sintering macroporosity does not exceed 30% by volume.

In all scenarios we considered, a primitive, undifferentiated crust survives thermal evolution (Elkins-Tanton et al. 2011; Sahijpal & Gupta 2011) and it possibly could be reduced by subsequent impacts.

A possible remanent magnetization is retained, if we assume a carbonaceous surface composition, since in all the scenarios we explored a surface layer of about 2–5 km is below the Curie temperature.

The main source of energy is represented by ²⁶Al, while the contribution due to ⁶⁰Fe and other possible sources (e.g., accretion and differentiation processes) is negligible.

We thank Guy J. Consolmagno and an anonymous referee for their helpful comments and Romolo Politi for his numerical analysis assistance. M.F. thanks his friend Demetra De Cicco for the revision of the paper. The computational resources used in this research have been supplied by INAF-IAPS through the project HPP-High Performance Planetology.