Rocklines as Cradles for Refractory Solids in the Protosolar Nebula

Meteorites and terrestrial planets show varying proportions of silicates and metallic iron, with Fe being distributed between Fe alloys and silicates (Nittler et al. 2000; Folco & Cordier 2015). Chondrites are the most common meteoritical bodies found on Earth, whose unaltered structure gives valuable information on the formation conditions of the building blocks of the solar system (Grossman 1972; Wasson & Kallemeyn 1988). Chondrules, which are round grains composed primarily of the silicate minerals olivine and pyroxene and found in the three families of chondrites (ordinary, carbonaceous, and enstatite chondrites), are then important to assess the thermodynamic evolution and the initial composition of the protosolar nebula (hereafter PSN) (Bouvier & Wadhwa 2010; Connelly et al. 2012; Bollard et al. 2017). Another source of meteoritical matter is cosmic spherules (CS), which are micrometeorites formed by the melting of interplanetary dust particles during atmospheric entry. CS are believed to sample a broader range of material than the collections of meteorites (Taylor et al. 2000). They are sorted by families and types, each one presenting its own structure and composition (Taylor et al. 2000; Alexander et al. 2002; Cordier et al. 2011; Sohl 2015). Neither chondrules, chondrites, or CS represent the earliest condensates of refractory phases in the PSN, most of which have been thermally processed before accretion on asteroids and planetary bodies. Their composition, however, can reasonably be used as a proxy of the composition of refractory solids in the PSN. The observed variations of composition of chondrites, chondrules, and CS suggest a compositional and redox gradient of refractory matter as a function of radial distance in the PSN and/or time. Recent studies explored the link between CS and chondrites, showing that it is possible to associate CS to their chondritic precursors, based on compositional analysis (Rudraswami et al. 2015; van Ginneken et al. 2017), potentially providing additional constraints to the thermodynamic conditions of the PSN.

The question of the origin of the composition differences among chondrites and chondrules is often answered by invoking the cooling of the inner disk where high temperature materials condense first (Grossman 1972; Taylor & Scott 2001, 2004). However, the influence of the temperature gradient within an evolving PSN, where more refractory materials form at closer distances from the Sun, has never been investigated to address this issue. Abundances of materials both in solid and gaseous phases are ruled by the chemical status of the disk and by the location of their condensation/sublimation fronts (Drażkowska & Alibert 2017; Öberg & Wordsworth 2019; Pekmezci et al. 2019). Significant increases of the abundances of solid materials can be generated at the location of these transition lines, due to the dynamics of vapors and grains (Cyr et al. 1998, 1999; Bond et al. 2010; Ali-Dib et al. 2014; Desch et al. 2017; Schaffer et al. 2018; Mousis et al. 2019). These processes potentially explain the metallicity of Jupiter through its formation near the water snowline (Mousis et al. 2019), as well as the high density of Mercury whose building blocks could have formed in regions where the abundances of Fe-bearing species are prominent (Taylor & Scott 2004; Vander Kaaden et al. 2019).

Here, we use a time-dependent coupled disk/transport model to investigate the role held by rocklines (the concept of snowlines extended to more refractory solids) of the most abundant solids into the shaping of Mg, Fe, and Si abundances profiles in the inner part of the PSN. The radial transport of solid grains through the different rocklines, coupled to the diffusion of vapors, leads to local enrichments or depletions in minerals and imply variations of the Mg-Fe-Si composition of dust grains in the inner regions of the PSN. We discuss our results in light of the relative abundances and compositions of minerals observed in meteoritic matter and planetary bulk compositions. Our approach can be used to derive the formation conditions of the primitive matter in the PSN, and to give insights on the origin of Mercury as well as Super-Mercuries (Santerne et al. 2018; Brugger et al. 2019).

Section 2 summarizes the disk/transport model used to compute the composition and the thermodynamic properties of the PSN. Our results are presented in Section 3. Section 4 is devoted to a discussion and our conclusions.

2.1. Disk Evolution

Our time-dependent PSN model is ruled by the following second-order differential equation (Lynden-Bell & Pringle 1974):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{g}}}}{\partial t}=\displaystyle \frac{3}{r}\displaystyle \frac{\partial }{\partial r}\left[{r}^{1/2}\displaystyle \frac{\partial }{\partial r}\left({r}^{1/2}{{\rm{\Sigma }}}_{{\rm{g}}}\nu \right)\right].\end{eqnarray} \tag{ 1 }$$

This equation describes the evolution of a viscous accretion disk of surface density Σ_g of dynamical viscosity ν, assuming hydrostatic equilibrium in the z direction. This equation can be rewritten as a set of two first-order differential equations coupling the gas surface density Σ_g field and mass accretion rate $\dot{M}$:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{{\rm{g}}}}{\partial t}=\displaystyle \frac{1}{2\pi r}\displaystyle \frac{\partial \dot{M}}{\partial r},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\dot{M}=3\pi {{\rm{\Sigma }}}_{{\rm{g}}}\nu \left(1+2Q\right),\end{eqnarray} \tag{ 3 }$$

where $Q=\mathrm{dln}({{\rm{\Sigma }}}_{{\rm{g}}}\nu )/\mathrm{dln}(r)$. The first equation is a mass conservation law, and the second one is a diffusion equation. The mass accretion rate can be expressed in terms of the gaz velocity field v_g as $\dot{M}=-2\pi r{{\rm{\Sigma }}}_{{\rm{g}}}{v}_{{\rm{g}}}$.

The viscosity ν is computed using the prescription of Shakura & Sunyaev (1973) for α-turbulent disks:

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

where ${{\rm{\Omega }}}_{{\rm{K}}}=\sqrt{{{GM}}_{\odot }/{r}^{3}}$ is the Keplerian frequency, with G as the gravitational constant, and c_s is the isothermal sound speed given by

$$\begin{eqnarray}&&{c}_{{\rm{s}}}=\sqrt{\displaystyle \frac{{RT}}{{\mu }_{{\rm{g}}}}}.\end{eqnarray} \tag{ 5 }$$

In this expression, R is the ideal gas constant and μ_g = 2.31 g.mol⁻¹ is the gas mean molar mass (Lodders et al. 2009). In Equation (4), α is a nondimensional parameter measuring the turbulence strength, which also determines the efficiency of viscous heating, hence the temperature of the disk. The value of α typically lies in the range 10⁻⁴–10⁻², from models calibrated on disk observations (Hartmann et al. 1998; Hueso & Guillot 2005; Desch et al. 2017).

The midplane temperature T of the disk is computed via the addition of all heating sources, giving the expression (Hueso & Guillot 2005):

$$\begin{eqnarray}&&\begin{array}{l}{T}^{4}=\displaystyle \frac{1}{2{\sigma }_{\mathrm{sb}}}\left(\displaystyle \frac{3}{8}{\tau }_{{\rm{R}}}+\displaystyle \frac{1}{2{\tau }_{{\rm{P}}}}\right){{\rm{\Sigma }}}_{{\rm{g}}}\nu {{\rm{\Omega }}}_{{\rm{K}}}^{2}\\ \,+\,{T}_{\odot }^{4}\left[\displaystyle \frac{2}{3\pi }{\left(\displaystyle \frac{{R}_{\odot }}{r}\right)}^{3}+\displaystyle \frac{1}{2}{\left(\displaystyle \frac{{R}_{\odot }}{r}\right)}^{2}\left(\displaystyle \frac{H}{r}\right)\left(\displaystyle \frac{\mathrm{dln}H}{\mathrm{dln}r}-1\right)\right]\\ \,+\,{T}_{\mathrm{amb}}^{4}.\end{array}\end{eqnarray} \tag{ 6 }$$

The first term corresponds to the viscous heating (Nakamoto & Nakagawa 1994), where σ_sb is the Stefan–Boltzmann constant, and τ_R and τ_P are the Rosseland and Planck mean optical depths, respectively. For dust grains, we assume τ_P = 2.4τ_R (Nakamoto & Nakagawa 1994). τ_R is derived from Hueso & Guillot (2005):

$$\begin{eqnarray}&&{\tau }_{{\rm{R}}}=\displaystyle \frac{{\kappa }_{{\rm{R}}}{{\rm{\Sigma }}}_{{\rm{g}}}}{2},\end{eqnarray} \tag{ 7 }$$

where κ_R is the Rosseland mean opacity, computed as a sequence of power laws of the form κ_R = κ₀ρ^aT^b, where parameters κ₀, a, and b are fitted to experimental data in different opacity regimes (Bell & Lin 1994) and ρ denotes the gas density at the midplane. The second term corresponds to the irradiation of the disk by the central star of radius R_⊙ and surface temperature T_⊙. It considers both direct irradiation at the midplane level and irradiation at the surface at a scale height H = c_s/Ω_K. The last term accounts for background radiation of temperature T_amb = 10 K.

At each time-step, Σ_g is evolved with respect to Equation (2). Then Equation (6) is solved iteratively with Equations (4), (5), and (7) to produce the new thermodynamic properties of the disk. The new velocity field is then computed following Equation (3). The time-step is computed using the diffusion timescale in each bin ${\rm{\Delta }}t=0.5\,min({\rm{\Delta }}{r}^{2}/\nu )\,\simeq 0.1\,{\rm{yr}}$, where Δr is the spatial grid size, and the factor 0.5 is taken for safety. The spatial grid is formed from N = 500 bins scaled as ${r}_{i}=({r}_{\min }-{r}_{\mathrm{off}})\times {i}^{\beta }+{r}_{\mathrm{off}}$, where $\beta \,=\mathrm{log}\left(\tfrac{{r}_{\max }-{r}_{\mathrm{off}}}{{r}_{\min }-{r}_{\mathrm{off}}}\right)/\mathrm{log}(N)$, giving a nonuniform grid, with r_min and r_max as the limits of the grid, and r_off is a parameter that allows us to control the spatial distribution of points in the simulation box. We compute the sum of the mass lost from both limits of the simulation box ($\int \dot{M}({r}_{\min }){dt}$ and $\int \dot{M}({r}_{\max }){dt}$) along with the mass of the disk itself at each time. This quantity remains constant within the precision of the machine.

The initial condition is the self-similar solution ${{\rm{\Sigma }}}_{{\rm{g}}}\nu \propto \exp \left(-{\left(r/{r}_{{\rm{c}}}\right)}^{2-p}\right)$ derived by Lynden-Bell & Pringle (1974). Choosing $p=\tfrac{3}{2}$ for an early disk, we derive the initial surface density and mass accretion rate as follows:

$$\begin{eqnarray}&&\left\{\begin{array}{l}{{\rm{\Sigma }}}_{{\rm{g}},0}=\displaystyle \frac{{\dot{M}}_{\mathrm{acc},0}}{3\pi \nu }\exp \left[-{\left(\displaystyle \frac{r}{{r}_{{\rm{c}}}}\right)}^{0.5}\right],\\ {\dot{M}}_{0}={\dot{M}}_{\mathrm{acc},0}\left(1-{\left(\displaystyle \frac{r}{{r}_{{\rm{c}}}}\right)}^{0.5}\right)\exp \left[-{\left(\displaystyle \frac{r}{{r}_{{\rm{c}}}}\right)}^{0.5}\right].\end{array}\right.\end{eqnarray} \tag{ 8 }$$

To compute the numerical value of ${{\rm{\Sigma }}}_{{\rm{g}},0}$, we solve iteratively Equation (6) with the imposed density profile in the first member of Equation (8). r_c regulates the size of the disk, and from the second member of Equation (8) we see that it matches, at t = 0 to the centrifugal radius.

The value of r_c is computed by dichotomy, assuming a disk mass M_d of 0.1 M_⊙. For α = 10⁻³, r_c is equal to 1.83 au, and 99% of the disk's mass is encapsulated within ∼100 au. We assume the initial mass accretion rate onto the central star ${\dot{M}}_{\mathrm{acc},0}$ to be ${10}^{-7.6}$ M_⊙· yr⁻¹ (Hartmann et al. 1998). The resulting density and temperature profiles are shown in Figure 1 at different epochs for α = 10⁻³.

Zoom In Zoom Out Reset image size

2.2. Size of Dust Particles

The size of dust particles used in our model is determined by a two-population algorithm derived from Birnstiel et al. (2012). This algorithm computes the representative size of particles through the estimate of the limiting Stokes number in various dynamical regimes. We assume that dust is initially present in the form of particles of sizes a₀ = 10⁻⁷ m, and grow through mutual collisions on a timescale (Lambrechts & Johansen 2014):

$$\begin{eqnarray}&&{\tau }_{\mathrm{growth}}=\displaystyle \frac{a}{\dot{a}}=\displaystyle \frac{4{{\rm{\Sigma }}}_{{\rm{g}}}}{\sqrt{3}{\epsilon }_{{\rm{g}}}{{\rm{\Sigma }}}_{{\rm{d}}}{{\rm{\Omega }}}_{{\rm{K}}}},\end{eqnarray} \tag{ 9 }$$

where a is the size of dust grains and Σ_d is the dust surface density. We set the growth efficiency parameter _g equal to 0.5 (Lambrechts & Johansen 2014). The size of particles that grow through sticking is thus given by ${a}_{\mathrm{stick}}={a}_{0}\exp (t/{\tau }_{\mathrm{growth}})$. However, this growth is limited by several mechanisms preventing particles from reaching sizes greater than ∼1 cm. The first limit arises from fragmentation, when the relative speed between two grains due to their relative turbulent motion exceeds the velocity threshold u_f. This sets a first upper limit for the Stokes number of dust grains, which is (Birnstiel et al. 2012):

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{frag}}=0.37\displaystyle \frac{1}{3\alpha }\displaystyle \frac{{u}_{{\rm{f}}}^{2}}{{c}_{{\rm{s}}}^{2}},\end{eqnarray} \tag{ 10 }$$

where we set u_f = 10 m.s⁻¹ (Birnstiel et al. 2012; Mousis et al. 2019).

A second limitation for dust growth is due to the drift velocity of grains, i.e., when grains drift faster than they grow, setting another limit for the Stokes number (Birnstiel et al. 2012):

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{drift}}=0.55\displaystyle \frac{{{\rm{\Sigma }}}_{{\rm{d}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\displaystyle \frac{{v}_{{\rm{K}}}^{2}}{{c}_{{\rm{s}}}^{2}}{\left|\displaystyle \frac{\mathrm{dln}P}{\mathrm{dln}r}\right|}^{-1},\end{eqnarray} \tag{ 11 }$$

where v_K is the Keplerian velocity, and P is the disk midplane pressure.

Equation (10) only considers the relative turbulent motion between grains that are at the same location, but the fragmentation threshold u_f can also be reached when dust grains drift at great velocities, and in the process collide with dust grains that are on their path. In that case, we obtain a third limitation for grains' size (Birnstiel et al. 2012):

$$\begin{eqnarray}&&{\mathrm{St}}_{\mathrm{df}}=\displaystyle \frac{1}{(1-N)}\displaystyle \frac{{u}_{{\rm{f}}}{v}_{{\rm{K}}}}{{c}_{{\rm{s}}}^{2}}{\left(\displaystyle \frac{{\rm{d}}\mathrm{ln}P}{{\rm{d}}\mathrm{ln}r}\right)}^{-1},\end{eqnarray} \tag{ 12 }$$

where the factor N = 0.5 accounts for the fact that only particles of bigger size fragment during collisions.

The relation between Stokes number and dust grain size depends on the flow regime in the disk (Johansen et al. 2014):

$$\begin{eqnarray}\mathrm{St}=\left\{\begin{array}{ll}\sqrt{2\pi }\displaystyle \frac{a{\rho }_{{\rm{b}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}} & \mathrm{if}\,{\rm{a}}\leqslant \displaystyle \frac{9}{4}\lambda \\ \displaystyle \frac{8}{9}\displaystyle \frac{{a}^{2}{\rho }_{{\rm{b}}}{c}_{{\rm{s}}}}{{{\rm{\Sigma }}}_{{\rm{g}}}\nu } & \mathrm{if}\,{\rm{a}}\geqslant \displaystyle \frac{9}{4}\lambda ,\end{array}\right.\end{eqnarray} \tag{ 13 }$$

where ρ_b is the bulk density of grains. The first case corresponds to the Epstein regime, occurring in the outermost region of the disk, and the second case corresponds to the Stokes regime. The limit between the two regimes is set by the mean free path $\lambda =\sqrt{\pi /2}\,\nu /{c}_{{\rm{s}}}$ (computed by equating both terms of Equation (13)) in the midplane of the disk. Combining Equations (9)–(13), the limiting Stokes number and representative particle size are computed for both regimes independently, giving a_E and St_E in the Epstein regime, and a_S and St_S in the Stokes regime. Then, as shown by Figure 2, the representative size and Stokes number is given by

$$\begin{eqnarray}&&a=\min \left({a}_{{\rm{E}}},{a}_{{\rm{S}}}\right),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\mathrm{St}=\max \left({\mathrm{St}}_{{\rm{E}}},{\mathrm{St}}_{{\rm{S}}}\right).\end{eqnarray} \tag{ 15 }$$

Zoom In Zoom Out Reset image size

In the following, two end cases are considered. In case A, we assume that all trace species are entirely independent, i.e., a run with several species is equivalent to several runs with a single trace species at a time. In case B, we assume that at each orbital distance, dust grains are a mixture of all available solid matter at that distance. For this case, the considered dust surface density is the sum over all surface densities ${{\rm{\Sigma }}}_{{\rm{d}}}={\sum }_{i}{{\rm{\Sigma }}}_{{\rm{d}},i}$. The bulk density ρ_b of resulting grains is also the mass-average of the bulk densities of its constituents (given in Table 1):

$$\begin{eqnarray}&&{\rho }_{{\rm{b}}}=\displaystyle \frac{{\displaystyle \sum }_{i}{{\rm{\Sigma }}}_{{\rm{d}},i}{\rho }_{{\rm{b}},i}}{{\displaystyle \sum }_{i}{{\rm{\Sigma }}}_{{\rm{d}},i}}.\end{eqnarray} \tag{ 16 }$$

Case B is favored from a dynamical point of view, as dust grains of different composition mix over long timescales, whereas case A is favored from a thermodynamic point of view, since sublimation and condensation tend to separate species into their pure forms. Since we focus on rocklines, whose positions are all at distances ≤1 au, we consider our disk volatile-free. In our model, the closest iceline would be that of H₂O, which is located at ∼4 au. This iceline is far enough to ignore its impact on processes at play around rocklines.

Table 1.  Main Refractory Phases Present in the Disk with Corresponding Initial Abundances and References for Saturation Pressures P_sat

Type	Name	Formula	ni = (Xi/H2)	ρb (g cm−3)	Psat
Pyroxene	Ferrosilite	FeSiO3	9.868 × 10−6	3.95	Nagahara et al. (1994)
	Enstatite	MgSiO3	4.931 × 10−5	3.20	Tachibana et al. (2002)
Olivine	Fayalite	Fe2SiO4	2.897 × 10−6	4.39	Nagahara et al. (1994)
	Forsterite	Mg2SiO4	1.452 × 10−5	3.22	Nagahara et al. (1994)
Alloys	Iron sulfide	FeS	3.266 × 10−5	4.75	Ferro et al. (1989)
	Kamacite	Fe0.9Ni0.1	1.875 × 10−5	7.90	Alcock et al. (1984)
	Nickel	Ni	1.875 × 10−6	8.90	Alcock et al. (1984)

Download table as:  ASCIITypeset image

2.3. Evolution of Vapors and Dust

We follow the approaches of Desch et al. (2017) and Drażkowska & Alibert (2017) for the dynamics of trace species in terms of motion and thermodynamics, respectively. The disk is uniformly filled with seven refractory species considered dominant (see Table 1), assuming protosolar abundances for Fe, Mg, Ni, Si, and S (Lodders et al. 2009), similar Fe/Mg ratios in olivine and pyroxene, and that half of Ni is in pure metallic form while the remaining half is in kamacite.

Sublimation of grains occurs during their inward drift when partial pressures of trace species become lower than the corresponding vapor pressures. Once released, vapors diffuse both inward and outward. Because of the outward diffusion, vapors can recondense back into solid form, and condensation occurs either until thermodynamic equilibrium is reached or until no more gas is available to condense. Over one integration time-step Δt, the amount of sublimated or condensed matter is (Drażkowska & Alibert 2017)

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{subl},i} & = & {\dot{Q}}_{\mathrm{subl},i}{\rm{\Delta }}t\\ & = & \min \left(\displaystyle \frac{6\sqrt{2\pi }}{\pi {\rho }_{{\rm{b}}}a}\sqrt{\displaystyle \frac{{\mu }_{i}}{{RT}}}{P}_{\mathrm{sat},i}{{\rm{\Sigma }}}_{{\rm{d}},i}{\rm{\Delta }}t;{{\rm{\Sigma }}}_{{\rm{d}},i}\right),\end{array}\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{{\rm{\Sigma }}}_{\mathrm{cond},} & = & {\dot{Q}}_{\mathrm{cond},i}{\rm{\Delta }}t\\ & = & \min \left(\displaystyle \frac{2H{\mu }_{g}}{{RT}}\cdot \left({P}_{{\rm{v}},i}-{P}_{\mathrm{sat},i}\right)\,;\,{{\rm{\Sigma }}}_{{\rm{v}},i}\right),\end{array}\end{eqnarray} \tag{ 18 }$$

where μ_i is the molar mass of a given trace species, P_sat,i is its saturation pressure, and ${P}_{{\rm{v}},i}$ is its partial pressure at a given time and place in the PSN. The second term of the min function ensures that the amount of sublimated (resp. condensed) matter is at most the amount of available matter to sublimate (resp. condense), i.e., the dust (resp. vapor) surface density Σ_d,i (resp. Σ_v,i). We define the rockline as the location at which the surface density Σ_d,i of solid and Σ_v,i and vapor of a given species i are equal. Species exist mostly in solid forms at greater heliocentric distances than their rocklines while they essentially form vapors at distances closer to the central star. Figure 3 shows the locations of the considered rocklines as a function of time in the PSN for both cases. No gas phase chemistry is assumed in the disk.

Zoom In Zoom Out Reset image size

The motion of dust and vapor, which coexist as separate surface densities Σ_d,i and Σ_v,i, is computed by integrating the 1D radial advection-diffusion equation (Birnstiel et al. 2012; Desch et al. 2017):

$$\begin{eqnarray}&&\displaystyle \frac{\partial {{\rm{\Sigma }}}_{i}}{\partial t}+\displaystyle \frac{1}{r}\displaystyle \frac{\partial }{\partial r}\left[r\left({{\rm{\Sigma }}}_{i}{v}_{i}-{D}_{i}{{\rm{\Sigma }}}_{{\rm{g}}}\displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{{{\rm{\Sigma }}}_{i}}{{{\rm{\Sigma }}}_{{\rm{g}}}}\right)\right)\right]+{\dot{Q}}_{i}=0.\end{eqnarray} \tag{ 19 }$$

This equation holds for both vapor and solid phases since the motion is determined by the radial velocity v_i and the radial diffusion coefficient D_i of species i, $\dot{Q}$ being the source/sink term. When a species i is in vapor form, we assume v_i ≃ v_g and D_i ≃ ν. When this species is in solid form, the dust radial velocity is the sum of the gas drag induced velocity and the drift velocity (Birnstiel et al. 2012):

$$\begin{eqnarray}&&{v}_{{\rm{d}}}=\displaystyle \frac{1}{1+{\mathrm{St}}^{2}}{v}_{{\rm{g}}}+\displaystyle \frac{2\mathrm{St}}{1+{\mathrm{St}}^{2}}{v}_{\mathrm{drift}},\end{eqnarray} \tag{ 20 }$$

where the drift velocity is given by (Weidenschilling 1977; Nakagawa et al. 1986):

$$\begin{eqnarray}&&{v}_{\mathrm{drift}}=\displaystyle \frac{{c}_{{\rm{s}}}^{2}}{{v}_{{\rm{K}}}}\displaystyle \frac{\mathrm{dln}P}{\mathrm{dln}r}.\end{eqnarray} \tag{ 21 }$$

This expression usually holds for a population of particles sharing the same size. However, we work here with the two-population algorithm from Birnstiel et al. (2012), i.e., the dust is composed of a mass fraction f_m of particles of size a, and a mass fraction $1-{f}_{{\rm{m}}}$ of particles of size a₀. The dust radial velocity can be then approximated by a mass-weighted velocity:

$$\begin{eqnarray}&&{v}_{i}={f}_{{\rm{m}}}{v}_{d,\mathrm{size}{\rm{a}}}+(1-{f}_{{\rm{m}}}){v}_{d,\mathrm{size}{{\rm{a}}}_{0}},\end{eqnarray} \tag{ 22 }$$

where ${v}_{d,\mathrm{size}{\rm{a}}/{{\rm{a}}}_{0}}$ is calculated with respect to Equation (20), St is computed for both populations at each heliocentric distance, and f_m depends on the size limiting mechanism (Birnstiel et al. 2012):

$$\begin{eqnarray*}{f}_{{\rm{m}}}=\left\{\begin{array}{ll}0.97 & \ \mathrm{if}\ \,{\mathrm{St}}_{\mathrm{drift}}=\min \left({\mathrm{St}}_{\mathrm{frag}},{\mathrm{St}}_{\mathrm{drift}},{\mathrm{St}}_{\mathrm{df}}\right)\\ 0.75 & \mathrm{otherwise}.\end{array}\right.\end{eqnarray*}$$

Due to the low Stokes number of dust grains (St < 1), we make the approximation ${D}_{i}=\tfrac{\nu }{1+{\mathrm{St}}_{i}^{2}}\simeq \nu$.

Figures 4 and 5 represent the time evolution of the composition of refractory particles evolving throughout the PSN in a Mg-Fe-Si ternary diagram and the time evolution of the Fe abundance profile (in wt%) as a function of heliocentric distance, and in case A and case B, respectively. Ternary diagrams display the composition of cosmic spherules and chondrules. Solid particles start with a protosolar composition, which changes during their drift throughout the innermost regions of the PSN, due to the successive sublimation of minerals. Because alloys, which contain only Fe, S, and Ni, are the first to sublimate (see Figure 3), solid particles loose a substantial amount of iron (in the ∼20%–50% range), but not Mg and Si. This is illustrated by the composition profiles that shift toward the Mg-Si axis with an unchanged Mg/Si ratio on the ternary diagrams. Closer to the Sun, ferrosilite and fayalite begin to sublimate until no Fe remains in solid form. In this case, the composition profiles in ternary diagrams are located on the Mg-Si axis, and because some silicon is vaporized too, the Mg/Si ratio increases. Finally, if the temperature is high enough to sublimate forsterite, only enstatite remains in solid particles, with an atomic ratio (Mg/Si)_at = 1 corresponding to 46 wt% of Mg. In case A, lighter grains are more subject to the pressure gradient and drift inward faster than dense Fe-rich grains. As a result, the Fe wt% increases by ∼10% in the whole disk. This leads to a slightly wider range of possible compositions in the PSN than for case B, but the differences between the two runs are minor from a compositional point of view.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 represents the composition of PSN refractory particles in case A (as in Figure 4) superimposed with the mean bulk compositions of chondrite groups. The figure also shows the composition of a random selection of chondrules and matrix (Hezel & Palme 2010). The mean bulk composition of chondrites is close to the protosolar value, and its spread seems to follow the profiles derived from our model. The same behavior can be observed for the matrix. On the other hand, chondrules exhibit a very low amount of bulk Fe, and the average composition seems to be close to the one computed by our model.

Zoom In Zoom Out Reset image size

Top panels of Figure 7 show the radial profiles of the disk's metallicity (defined as Z = Σ_d/Σ_g) at different epochs of its evolution. As expected, solid matter is concentrated at the position of rocklines. The composition of the PSN around the first cluster of rocklines (iron sulfide, kamacite, and nickel) corresponds to the 30–60 wt% Fe part of curves in ternary diagrams and Fe wt% profiles, which matches the S-BO type (barred olivine) spherules compositions. The composition of the PSN around the second cluster of rocklines (fayalite and ferrosilite) corresponds to the 10–30 wt% Fe part of curves in the ternary diagram and Fe wt% profiles, matching the S–V type (glass) spherules compositions. In the same manner, we would expect chondrules to be formed in the innermost regions if a sufficient amount of material were present due to continuous drift from the outer disk. At t = 10⁵ yr and r = 0.67 au (rockline of iron sulfide), the PSN has 56 wt% and 58 wt% of Fe in case A and case B, respectively. This increase of the Fe wt% leads to compositions of the PSN richer in Fe than the protosolar value.

Zoom In Zoom Out Reset image size

Finally, bottom panels of Figure 7 show the Stokes number of dust grains in the disk. Because case A has many independent species evolving, we only follow forsterite, fayalite, and iron sulfide, namely the most, least, and intermediate refractory materials considered in our particles. In case B, all grains are mixed together, giving a single Stokes number at each heliocentric distance. Using Equations (10), (11), and (12), we expect (i) St_frag ∝ 1/T, (ii) St_drift ∝ Z/(rT), and (iii) ${\mathrm{St}}_{\mathrm{df}}\propto 1/(\sqrt{r}T)$ (assuming ${\left|\tfrac{\mathrm{dln}P}{\mathrm{dln}r}\right|}^{-1}\propto 1$). In the innermost region, dust size is limited by fragmentation up to ∼5 au. In the 5–10 au range, a competition between drift and drift-limited fragmentation sets the dust grains' size. Beyond 10 au dust is in the growth phase.

In case A, variations in limiting sizes only come from the difference in bulk densities of grains. As a result, different species display the same Stokes number during most of their drift throughout the PSN (see the bottom left panel of Figure 7). However, because the amount of solid matter decreases below 10 au, as a result of sublimation and/or radial drift, the Stokes number diminishes as well. In case B, minor species embedded in large grains are transported more efficiently toward their rocklines. Hence, higher metallicities are found around the rocklines in case B compared to case A, at early epochs. In turn, the PSN becomes depleted in solid matter in shorter timescales in case B compared to case A.

As indicated in Section 2.1, α is a free parameter whose value can change with time and heliocentric distance (Kalyaan et al. 2015). For this disk model, an increase of the α value leads to a centrifugal radius r_c located at higher heliocentric distance, which in turn leads to a larger disk. This also leads to a larger diffusion coefficient D_i = ν for the trace species. As a consequence, vapors diffuse outward faster and enrich the solid phase more evenly. This results in peaks of abundances that are wider and less prominent than those observed in the right panels of Figures 4 and 5. For the extreme case α = 10⁻², the peaks of abundance are not observable anymore. This shows that the choice of α is critical for both the PSN and trace species evolutions. However, results of simulations with nonuniform α show that this quantity is increasing with heliocentric distance, and takes values of ∼10⁻³ at 1 au (Kalyaan et al. 2015). Since we are mostly interested in the dynamics of the inner PSN, we adopted this value for our α parameter.

Our model shows that the diversity of the bulk composition of cosmic spherules and chondrules can be explained by their formation close to rocklines, suggesting that solid matter is concentrated in the vicinity of these sublimation/condensation fronts. The slightly lower Fe content of S-type (ordinary) chondrites than for C-type (carbonaceous) observed in Figure 6 could be due to a partial sublimation mechanism.

Several transport mechanisms can explain the presence of these processed minerals in the Main Belt, as well as in the outer regions of the PSN. For example, small particles with sizes within the 10⁻⁶–1 m range and formed in the inner nebula can diffuse radially toward its outer regions (Bockelée-Morvan et al. 2002; Desch et al. 2017). Figure 8 illustrates the efficiency of diffusion in our model by representing the time evolution of the radial distribution of particles initially formed in the 0–1 au region of the PSN. In this simulation, a protosolar dust-to-gas ratio is assumed within 1 au at t = 0 and the effects of rocklines are not considered. The figure shows that particles can diffuse well over 10 au after 10⁵ yr of PSN evolution and thus fill the Main Belt region. Photophoresis is another mechanism that can be at play when the inner regions of the PSN become optically thin. In this case, the disk still has a reasonable gas content and allows particles with sizes in the 10⁻⁵–10⁻¹ m range to receive light from the proto-Sun and be pushed outward by the photophoretic force beyond 20 au in the PSN (Krauss et al. 2007; Mousis et al. 2007; Moudens et al. 2011; Loesche et al. 2016). These two mechanisms would have taken place before the formation of larger planetesimals in the PSN.

Zoom In Zoom Out Reset image size

Interestingly, these high concentrations do not allow the triggering of streaming instability. The smallest metallicity Z_c required to trigger a streaming instability for low Stokes number regimes $\left(\mathrm{St}\lt 0.1\right)$ is (Yang et al. 2017):

$$\begin{eqnarray}&&\mathrm{log}{Z}_{{\rm{c}}}=0.10{\left(\mathrm{logSt}\right)}^{2}+0.20\mathrm{logSt}-1.76.\end{eqnarray} \tag{ 23 }$$

Figure 7 shows that the highest disk metallicity computed with our model, i.e., at t = 10⁵ yr and r = 0.1 au in case B, is Z = 4.4 × 10⁻³, way below Z_c = 2. 7 × 10⁻², suggesting our model does not allow the triggering of a streaming instability. However, the decoupling of dust grains from gas at greater Stokes numbers or the back-reaction of dust onto the gas could slow down the drift of particles in the pile-up regions, thus increasing the local dust-to-gas ratio.

Although the spread in bulk compositions of CS can be well explained by alteration during atmospheric entry (Rudraswami et al. 2015), the two scenarios are not mutually exclusive as they happen at very different times. The grains composition computed by our model can be seen as the average composition at each time and heliocentric distance, and the deviation from the mean composition could be the result of full dynamics of grain growth and interaction (the same ideas can be found in Horner et al. 2007). The effect of rocklines could be among the first processing mechanisms altering the uniform composition of refractory matter in the PSN. Moreover, condensation/sublimation fronts for refractory matter may have other implications. Close to the host star, the pressure at the midplane can be high enough (up to 1 bar) to melt partially or entirely solid grains. This would highly affect the collisional dynamics of grains (or droplets), and allow grains to overcome the meter barrier (Boley et al. 2014) and form Fe-rich planetesimals that later gave birth to Mercury.

Although our model relies a lot on the number of considered species and the availability of thermodynamic data governing state change, it suggests that rocklines played a major role in the formation of small and large bodies in the innermost regions of the PSN. For example, even if the large amount of iron in Mercury (83 wt% in the ternary diagram) cannot be explained with this model alone, the increased proportion of Fe in the PSN (62 wt% at most in the vicinity of rocklines; see Figures 4 and 5) from the protosolar value (47 wt%) might have contributed to the accretion of the planet's large core by forming Fe-rich regions. As our model only tracks the evolution of dust grains in the early PSN, it is compatible with current formation mechanisms of terrestrial planets (Raymond et al. 2014; Izidoro & Raymond 2018). The relevant PSN composition in terms of age and heliocentric distance must then be chosen accordingly to the considered formation scenario.

In its current state, the model fails to reproduce the extreme enrichments in Fe needed to account for the formation of Mercury. However, giant impact simulations performed by Chau et al. (2018) show that the resulting Mercury-like planets display core mass fractions (CMF) in the 0.5–0.7 range (see their Figure 3), when starting with a protosolar CMF of 0.3. If the initial CMF of Mercury was 10% higher due to the formation of its building blocks close to rocklines, as suggested by our findings, the post-collision CMF would lie in the 0.6–0.8 range, which is in better agreement with the estimated value of 0.7 (Stacey 2005; Righter 2007; Hauck et al. 2013). The combination of multiple scenarios to explain the large CMF of Mercury seems more likely. At greater heliocentric distances the PSN composition becomes again protosolar (mainly for case B), which is in agreement with the bulk composition of Earth and Venus (close to protosolar) derived from interior structure models (Stacey 2005; Rubie et al. 2015; Dumoulin et al. 2017).

Our study suggests that Mercury-like planets should exist in other planetary systems. More than ∼2000 small planets in the 1–3.9 ${R}_{\oplus }$ range have been discovered so far at close distances to their host stars.⁴ Most of the measured densities are poorly determined and the detection of sub-Earth planets remains challenging, implying it remains difficult to quantify the size of the population of Mercury-like planets. Finally, the presence of Mercury-like planets should be ruled by the amount of available matter to form Fe-rich building blocks. In the case of very massive and hot circumstellar disks, rocklines and Mercury-like planets would be located at much higher distances to the host star. In contrast, less massive and colder disks could impede the formation of Fe- rich planets because rocklines would be located too close to their host star.

O.M. acknowledges support from CNES. We thank the anonymous referee for useful comments that helped improving the clarity of our paper.