© The Authors 2024

1 Introduction

New stars are born when the cold cores of molecular clouds collapse under the pull of their own gravity. Because of the conservation of angular momentum during cloud collapse, a rapidly rotating disc is formed around the central protostar (Shu et al. 1987; Williams & Cieza 2011). The chemical composition of this proto-planetary disc, and the physical processes occurring therein, determine the elemental and molecular environment in which new planetary systems form (Öberg et al. 2023).

Important information about the environmental conditions in the early Solar System – that is, before the planets formed – can be deduced from the study of the material components found in primitive meteorites (chondrites): the matrix, the silicate-rich spherical chondrules, and the Ca-Al-rich inclusions (CAIs, Kita et al. 2012; Kruijer et al. 2017).

The matrix is the background material into which the other constituents – such as the chondrules, the CAIs, and larger mineral grains – are embedded (Brearley 1989). The matrix material consists mainly of anhydrous meteoritic phases, such as olivine, pyroxene, and Fe-Ni metals that have been produced by transient heating and rapid cooling about 2–4 Myr after the CAIs originated. Smaller amounts of sulphides, sulphates, carbonates, and several other minerals are also found, and some meteorites also contain significant amounts of phyllosilicates (Brearley 1989; Buseck & Hua 1993).

CAIs are thought to be condensation products of the gas near the star, which had a composition that is similar to that of the gas in the planet-forming regions (Bekaert et al. 2021). The analysis of short- and long-lived isotopes, together with petrologic studies of these chondritic components, suggests that the CAIs were the first solid matter; therefore, their age can be equated with the age of the cosmochemical origin of the Solar System (Lodders & Amari 2005; Chaussidon & Liu 2015; Bermingham et al. 2020; Connelly et al. 2012). Four CAIs have been dated using the Pb-Pb chronometer, which yields a Solar System age of 4567.30 ± 0.16 Myr (Amelin et al. 2010; Connelly et al. 2012). Recently, Piralla et al. (2023) re-evaluated the age of the CAIs to 4568.7 Myr through unified Al-Mg and Pb-Pb chronology.

Chondrules are major constituents of most chondrites; they are approximately millimetre(mm)-sized, round, and silicate-rich grains that are assumed to have been entirely or partly molten at some point in their history in the solar nebula. According to Scott & Krot (2005), the chondrules crystallised in minutes to hours at temperatures of 1300–1800 K prior to accretion onto their parent bodies. Chondrules can be separated into two types based on the Mg number, which is controlled by the oxygen fugacity of the chondrule-forming environment in the disc, where type I chondrules formed under more highly reducing conditions than type II chondrules. A comprehensive list of minerals that were discovered in chondrules along with all their chemical variations is given in the review of Brearley et al. (1998).

All internal rearrangement and phase-change processes within the chondrules require high temperatures – particularly the melting. However, refractory materials such as quartz and corundum require pressures of over 8 bar and 7000 bar, respectively, in order to have a liquid phase in a solar composition gas; see Appendix B. At lower pressures, there is only sublimation and resublimation, but no melting. Given that such large gas pressures are not available during disc evolution, we suggest the observations of ‘molten’ materials be interpreted as annealing. A true melting of refractory materials seems to only be possible inside larger, internally heated, gravitationally bound bodies such as planetesimals.

The chemical analysis of chondrules can be used for an estimation of the lifetime of the solar nebula. It is expected that most chondrules formed ≈1–3 Myr after the origin of the CAIs, in a period when the proto-Sun accreted more slowly (Pape et al. 2019). Bollard et al. (2017) analysed the 22 youngest chondrules and determined their ages from measurements of Pb isotopic ratios. Bollard et al. concluded that the production and melting of chondrules began contemporaneously with CAI condensation, and continued for about 4 Myr. This lifetime estimate is in good agreement with observations of protoplanetary discs, which show lifetimes of about 1–10 Myr (Montmerle et al. 2006; Mamajek 2009).

Isotopic data of chondritic components in the CAIs are an excellent tool for understanding the formation of the planetary building blocks as well as their mixing and transport in the disc. It is expected that CAIs formed in low-pressure environments where the temperatures are larger than about 1300 K in an ¹⁶O-rich reducing environment (Grossman & Larimer 1974; Krot et al. 2018; MacPherson et al. 2005). Recently, Nakashima et al. (2023) analysed CAIs that were discovered in returned samples of the Japanese mission Ryugu (Hayabusa 2). It was found that three chondrule-like objects with a size of < 30 μm are dominated by Mg-rich olivine, and are enriched in ¹⁶O compared to ¹⁷O by ~ (3–23)%_o, which is in agreement with what has been proposed as early generations of chondrules (Nakashima et al. 2023). According to Nakashima et al. (2023), the ¹⁶O-enriched objects are most likely melted amoeboid olivine aggregates that escaped incorporation into the ¹⁶O-poor precursor dust of the chondrules. From their analysis of the Ryugu samples, these authors further deduce that small objects experienced radial transport in the inner solar nebula to the location of the formation of the Ryugu parent body, which was farther from the Sun and scarce in chondrules. Moreover, Nakashima et al. (2023) suggest that the transported objects were most likely destroyed during aqueous alteration in the Ryugu parent body. Further results from the analysis of these chondritic components include nucle-osynthetic Ti isotope variations and their combination within different bulk chondrites, as well as oxygen isotope ratios.

Amoeboid olivine aggregates (AOAs) are another class of refractory objects that are related to CAIs; however, they escaped extensive melting and thus preserved the internal structure they acquired during their formation in the pristine solar nebula. Together with chondrules and matrix, CAIs and AOAs are the main refractory components of chondritic meteorites that apparently formed under different conditions and in different regions of the disc and therefore comprise isotopically different materials (Krot 2019; Amelin et al. 2010; Desch et al. 2022; Öberg et al. 2023).

Although there has been significant progress over recent years regarding the emergence and implications of CAIs, there are still many unanswered questions and long-standing problems related to the physical, chemical, and environmental conditions that lead to the formation of CAIs and their subsequent storage and transport in the disc (MacPherson et al. 2005; Krot et al. 2015; Che & Brearley 2021). To better understand these processes, we developed a model that combines disc evolution with 2D radiative transfer, equilibrium condensation, and new dust opacity calculations; see Sect. 2. In Sect. 3, we apply this model to the early solar nebula and present the results. In Sect. 4, we discuss the implications of our findings concerning the CAI formation. Section 5 concludes our work.

Fig. 1 Three modelling stages, passed quantities, and results.

2 The model

Our model for the early evolution of the Solar System is generated in three modelling stages, see Fig. 1. First, we run a one-dimensional disc evolution model developed by Drążkowska & Dullemond (2018, 2023) for an initial cloud mass of 1.15 M_⊙, which results in a stellar mass of about 1 M_⊙ after 3 Myr; see Sect. 2.1. Second, selecting a sequence of ages, we import the disc structures from stage 1 into the 2D radiation thermo-chemical disc model ProDiMo (Woitke et al. 2009, 2016, 2024) and run detailed radiative transfer calculations; see Sect. 2.2. And third, we import the 2D disc structures from stage 2 into the chemical and phase equilibrium code GGchem (Woitke et al. 2018) to investigate the thermal stability of the refractory dust and the ices. In this last modelling stage, we recompute the solid opacities and modify the temperature structure in the inner, viscously heated parts of the disc close to the midplane, until the opacities are consistent with diffusive radiative transfer. This reveals the thermostat regulation mechanism discovered by D’Alessio et al. (1998, 1999) and studied in detail by Min et al. (2011), where the sublimation of the main silicates leads to a loss of opacity, which determines how efficiently the viscous heat can be transported away, and thereby regulates the midplane temperatures. In the following sections, we explain the three modelling stages in detail. A discussion of the shortcomings of our modelling approach can be found in Sect. 4.4.

2.1 The evolutionary disc model

We use a hydrodynamic disc evolution model similar to the one described in Dullemond et al. (2006a,b); Drążkowska & Dullemond (2018, 2023). Our model considers the formation of a star and circumstellar disc from a single parent molecular cloud core with a mass of 1.15 M_⊙, constant temperature of 10 K, and rotation rate of 7 × 10⁻¹⁵ s⁻¹. The collapse proceeds from inside-out. We assume that angular momentum is conserved during the infall and thus the material from every shell of the parent cloud falls onto the disc inside of the centrifugal radius. The centrifugal radius moves outward with time, which strongly depends on the rotation rate of the parent cloud. The lower the initial cloud rotation rate, the smaller the centrifugal radius, and thus the material falls on the hot part of the disc for a longer time. The way the material is distributed inside of the centrifugal radius is described in Hueso & Guillot (2005). We assume that all material falling inside of the disc’s inner radius, set here to 0.1 au, is added directly to the central star, and thus the star is fed both directly from the molecular cloud and through the disc accretion. The inside-out infall leads to an initially compact and hot disc, which reaches its maximum mass of 0.3 M_⊙ after about 0.6 Myr.

The disc surface density evolution is calculated by taking into account the mass infall from the parent cloud and viscous evolution which we parameterise using the standard α-formalism and setting the global value to α = 10⁻³ (Shakura & Sunyaev 1973). During the infall, the compact disc vigorously expands and the outward-flowing gas drags mm-sized solid particles with it. Such disc models have been shown to produce the right conditions for outward transport of the CAIs necessary to explain their incorporation in meteorites formed in the outer parts of the Solar System (Cuzzi et al. 2003; Ciesla 2007; Jacquet et al. 2019; Jongejan et al. 2023). However, the formation of CAIs was not studied consistently in these models; instead, these were added afterwards.

2.2 The ProDiMo disc model

The Protoplanetary Disc Model (ProDiMo) developed by Woitke et al. (2009, 2016, 2024) is a state-of-the-art code that combines disc chemistry with radiative transfer, heating/cooling and high-energy physics. The considered chemical processes include 2-body and 3-body gas-phase chemistry (Kamp et al. 2017), X-ray ionisation (Aresu et al. 2011), FUV and EUV photo-processes (Rab et al. 2017), and ice chemistry. The ProDiMo models in this paper are based on the radial disc structures imported from the evolutionary disc model (Sect. 2.1) for a sequence of eight consecutive ages 0.001, 0.003, 0.01, 0.03, 0.1, 0.3,1, and 3 Myr, where time zero is the beginning of the cloud collapse. For each age we use the following data as input: the stellar mass M_*, the stellar mass accretion rate Ṁ, the disc surface density profile Σ(r) [g cm⁻²], and the vertically integrated viscous heating rate F_vis(r) [erg s⁻¹ cm⁻²].

2.2.1 Stellar setup

The stellar mass M_* and the mass accretion rate Ṁ_acc are taken from the evolutionary disc model (Sect. 2.1). The stellar luminosity L_* and the effective temperature of the star T_* are interpolated between the values given in the pre-main sequence evolutionary tracks of Siess et al. (2000) for given age and M_*; see Table 1. We note that there is considerable uncertainty in the stellar luminosity at the earliest stages since the protostar is deeply embedded within the natal cloud and thus, not directly observable (White et al. 2007). The derived values of L_* and T_* depend on the particular pre-main sequence evolutionary model selected to calculate the tracks. When comparing against the tracks of (Yorke & Bodenheimer 2008), we find a good agreement concerning T_*, but the stellar luminosities differ notably before 0.3 Myr, with luminosities being almost one order of magnitude lower at the earliest evolutionary stages. However, these disagreements diminish with age and differ only about 10% at later times.

The stellar radius R_* is calculated from the Stefan–Boltzmann law , where L_* represents the energy liberated by the contraction of the protostar and the beginning nuclear fusion. The accretion luminosity L_acc represents the energy liberated by the gas falling onto the star, from the inner disc radius R_in to the stellar radius R_*. It is calculated as

(1)

The inner radius of the disc R_in is set to where the total stellar irradiation causes the dust to sublimate, considering a dust temperature of about 1500K:

(2)

The value of 0.1 au for a total luminosity of 1.95 L_⊙ comes from previous experience with the ProDiMo disc models.

Equations (1) and (2) are solved iteratively until converged. Table 1 summarises the stellar setup parameters for the ProDiMo models. For comparison, the table also contains the viscous disc luminosity,

(3)

where F_vis is the viscous heating rate per surface area, and R_out is the outer disc radius. Our results indicate that L_vis < 0.1L_* at any evolutionary state considered here.

The stellar spectra required for the ProDiMo radiative transfer calculations are compiled from two components: (i) the photospheric spectrum, representing L_* and T_*, and (ii) a blackbody of T_acc = 10000 K, representing the accretion luminosity L_acc produced by the accretion shock. The total spectral intensities emitted from the star’s surface are calculated as

(4)

The photospheric intensities are interpolated from a grid of PHOENIX stellar atmosphere models by Brott & Hauschildt (2005), where the surface gravity is given by .

2.2.2 Disc structure, dust settling, and radiative transfer

Based on the input gas surface density structure Σ(r), ProDiMo computes an axisymmetric 2D density structure n_〈H〉(r, z), where n_〈h〉 [cm⁻³] is the hydrogen nuclei density, in vertical hydrostatic equilibrium as follows. For each column r, starting with a guess of the temperature structure Τ(r, z), ProDiMo numerically integrates the equation of hydrostatic equilibrium upwards (Woitke et al. 2009):

(5)

using the ‘simplified hydrostatic model’ (Woitke et al. 2016), where the isothermal sound speed is evaluated with Τ = T_d and a constant (H₂-rich) mean molecular weight. The resulting density structure ρ(r, z) is then integrated vertically and normalised to match the given Σ(r). The differences to a “full hydrostatic model”, where the gas temperature T_g and the proper mean molecular weight are used to calculate the gas pressure, are visualised in Fig. 11 of Woitke et al. (2016). Since we focus on the dust evolution in the midplane in this paper, where T_g = T_d is valid, this approximation has little effect on the results.

Next, ProDiMo sets up the settled dust structure and calculates the frequency-dependent dust opacities as explained in Woitke et al. (2016). We start with an unsettled size distribution function with a powerlaw between a minimum (a_min) and a maximum dust size (a_max). Using 100 dust size bins, the dust in each bin is then vertically redistributed (settled) in each column according to the given gas density structure, the turbulent mixing parameter a_settle, and the settling description of Riols & Lesur (2018); see Woitke et al. (2024) for details. We used the following parameter values in this paper: a_min = 0.05 μm, a_max = 3 mm, a_pow = 3.5, dust-to-gas mass ratio = 0.004, and a_settle = 10⁻³. The dust opacities are then calculated by Mie theory, assuming a constant mixture of 60% silicates, 15% amorphous carbon and 25% porosity, using the settled dust size distribution function f(a, r, z) at each point of the model.

To compute T_d(r, z), ProDiMo uses its in-built frequency-dependent ray-based 2D radiative transfer module, which has been thoroughly tested against the leading Monte-Carlo radiative transfer programs by Pinte et al. (2009); see also Woitke et al. (2016) and Appendix A in Oberg et al. (2022). In the optically thick and viscously heated midplane, most relevant for this paper, the conservation of the frequency-integrated radiative flux can be formulated as

(6)

where Γ is the local net heating rate converting non-radiative to radiative energy. The input model (Sect. 2.1) provides the viscous heating rate per unit area F_vis [erg cm⁻² s⁻¹]. To calculate the heating rate per volume Γ, we assume that the heating occurs proportional to the gas density ρ as

(7)

where the factor 1/2 is because of the mirror symmetry at z = 0. In the diffusion approximation, the bolometric radiation flux,

(8)

is given by the gradient of the bolometric mean intensity J(r, z) = σ Τ⁴/π, where σ is the Stefan-Boltzmann constant. The Rosseland-mean opacity is defined as

(9)

where is the dust extinction coefficient at position (r, z) in the disc, Τ = T_dust(r, z) is the temperature, and B_v(T) is the Planck function. ProDiMo solves these diffusion equations in 2D for given opacity structure in the optically thick core of the midplane, using the temperature results from the proper 2D frequency-dependent ray-based radiative transfer as upper boundary condition. We assume that the vertical component of the flux is zero in the midplane, because of symmetry, and use this condition as lower boundary condition; see Appendix A in Oberg et al. (2022).

Once the 2D disc radiative transfer problem is solved, and the dust temperature structure T_d(r, z) is determined, ProDiMo calculates the chemical composition in time-independent mode and the gas temperature structure T_g(r, z) in heating/cooling balance. These results are not relevant to the problem of dust stability in the midplane, and will hence not be discussed any further in this paper. However, T_d(r, z) is required for the computations of the vertical hydrostatic equilibrium (Eq. (5)), and therefore, ProDiMo performs global iterations between radiative transfer and hydrostatic equilibrium, which converges after about ten iterations.

2.3 The GGchem phase equilibrium model

The basic idea for this paper is to consider the process of refractory dust sublimation as post-processing of the ProDiMo disc models, using GGchem (Woitke et al. 2018). GGchem is a publicly available thermo-chemical equilibrium code to calculate the concentrations of all neutral and single ionised atoms, electrons, molecules, molecular ions, and condensates. GGchem is based on the principle of minimisation of the total Gibbs free energy, including condensation, and applicable in a wide temperature range between 100 and 6000 K. For this paper, we use the code down to 50 K; for temperatures below 50 K we use Τ = 50 K.

The sole purpose of GGchem here is to predict the amount and volume composition of all solids, for given gas density n_〈h〉, temperature Τ, and total element abundances prior to condensation, where we use the solar abundances of Asplund et al. (2009). The GGchem setup in this paper is the same as in Woitke et al. (2018). We select 26 elements (H, He, Li, C, N, O, F, Na, Mg, Al, Si, P, S, Cl, K, Ca, Ti, V, Cr, Mn, Fe, Ni, Cu, Zn, Zr, and W) and include condensates, ions and charged molecules, which leads to a selection of 595 gaseous chemical species and 248 condensed species.

For this paper, we extended GGchem to calculate the opacity of the resulting mixture of solid materials as explained in Sect. 2.3.3. However, this opacity structure changes the diffusive radiative transfer in the optically thick core of the disc, and therefore, dust stability and radiative transfer become a coupled problem. The following section describes how we solve this issue using a simplified vertical diffusive radiative transfer, which results in a thermostat regulation mechanism, which changes Τ(r, z) in the optically thick midplane, and only there.

2.3.1 The thermostat regulation mechanism

When the midplane is very optically thick and heated by fric-tional processes, the radiative flux is almost perfectly directed upwards in the vertical direction. In this case, we can simplify Eqs. (6) and (8) as

(10)

and hence the vertical flux F(z) is known because we know Γ(z). Considering a vertical grid {z_k | k = 0,…, K}, taken over from the ProDiMo disc model, Eq. (10) is written numerically as

(11)

where z_k−1/2 = (z_k−1 + z_k)/2. Equation (11) can be solved downwards (k → k − 1) to obtain the consistent temperature and opacity structure. For each disc column at radius r, we start the integration at a height z_K over the midplane, where the vertical and horizontal Rosseland optical depth are both >10, to make sure we are safely in the diffusive regime. We then read off Τ(z_K) at that point as upper boundary condition, and the density structure n_〈h〉(z) below (z ≤ z_K) from the ProDiMo disc model. Then, we

1. guess the temperature Τ(z_k−1),

2. compute the amount of dust and its material composition for given total element abundances, temperature Τ(z_k−1) and hydrogen nuclei particle density n_〈H〉(z_k−1) with GGchem,

3. calculate the dust size distribution function f(a) at z_k−1 (Sect. 2.3.2),

4. calculate the frequency-dependent dust extinction coefficient according to Eq. (15) by effective medium and Mie theory (Sect. 2.3.3),

5. calculate the Rosseland mean opacity κ_R(z_k−1) according to Eqs. (9) and (18),

6. calculate Τ(z_k−1) according to Eq. (11),

7. go back to step 2 if not yet converged.

Using this procedure, we implemented a stable thermostat regulation mechanism as discovered by Min et al. (2011). If the temperature Τ(z_k−1) becomes too high, the dust starts to sublimate, which lowers the Rosseland opacity κ_R, and hence limits the downward temperature increase. In all cases where both the viscous heating rate F_vis and the vertical optical depths are large enough, the solution will always lead to silicate sublimation in the deeper layers, but rarely to complete dust evaporation, as some of the most stable condensates will remain: the Al-Ca-Ti oxides.

2.3.2 Dust size distribution

For the dust size distribution function f(a), we consider a reference dust size distribution function f_ref(a) ∝ a^−3.5 (in units cm⁻¹/H-nucleus) between 0.05 and 3000 μm, which is normalised to a dust/gas mass ratio of 0.004, assuming an interstellar dust material density of 2 g cm⁻³. The value 0.004 results from complete condensation of a solar composition gas (Asplund et al. 2009) before ice formation; see Woitke et al. (2018) for details. The reference dust size distribution function is the same as the unsettled dust size distribution function in the ProDiMo disc model. The total reference volume of condensed species is

(12)

In the phase equilibrium simulations, GGchem provides a different total volume of condensed species V per H-nucleus, depending on which condensed materials are found to be stable, and how much of them. If V < V_ref, which happens, for example, when the silicates sublimate, we assume that all particles have shrunk by a common factor in radius

(13)

The numerical grid for the size distribution function {(a_i, f_i) | i = 0,…, I} is then modified as a_i = γ × a_ref,i, while the number of dust particles between size gridpoints a_i is assumed to remain the same. However, since f(a) is the number of dust particles per size interval, we have f_i = f_ref,i/γ. If ice condenses on the existing surfaces, leading to V > V_ref, we assume that all grains increase in size accordingly, which leads to increased dust opacities per hydrogen nucleus according to Eq. (15), first and foremost because of the increased condensed volume per H-nucleus.

2.3.3 Dust and gas opacities

The phase equilibrium model GGchem provides the volume mixing ratios {V_s | s = 1,…,S}, where s is an index for the stable solids; see examples in Table 2. We add 25% porosity. Using the effective medium theory of Bruggeman (1935), we calculate the effective optical constants of the mixed material as

(14)

where (n_s, k_s) are the real and imaginary parts of the optical constants of the pure materials s. The dust extinction, absorption and scattering opacities [1/cm] are then calculated as

(15)

(16)

(17)

where a is the radius of a dust particle and x = 2π a/λ is the size parameter. f (a) [1/cm/H-nucleus] is the size distribution function per hydrogen nucleus; see Sect. 2.3.2. The extinction, absorption and scattering efficiencies are calculated according to Mie theory. We use the routine MIEX (Voshchinnikov 2002) in the implementation of Sebastian Wolf. The opacity computations are based on a collection of optical data by Marigo et al. (2024), with additional data collected by Carsten Dominik and Michiel Min¹. The complete selection of relevant materials and opacity sources is summarised in Table A.1.

Finally, the total Rosseland opacity used in Eq. (11) is assumed to be given by the dust Rosseland opacity according to Eqs. (9) and (15), and a gas Rosseland opacity

(18)

The gas Rosseland mean opacity is interpolated from 2D tables provided by Paola Marigo, using the same opacity data sources and techniques explained in Marigo et al. (2022, 2024), but without the dust opacities. In those works we updated and expanded molecular absorption to include 80 species, predominantly utilising the recommended line lists currently accessible from the ~~ExoMol~~~ (Tennyson & Yurchenko 2012) and ~~HITRAN~~~ (Gordon et al. 2022) databases. Additionally, in response to a recent investigation, we revised the H⁻ photodetachment cross section, incorporated the free-free absorption of negative ions of atoms and molecules, and refreshed the collision-induced absorption due to H₂/H₂, H₂/H, H₂/He, and H/He pairs. To construct the Rosseland mean opacity tables, we integrated the ӔSOPUS 2.0 code (Marigo et al. 2024) with the GGchem code (Woitke et al. 2018). We adopted a solar mixture following Asplund et al. (2009), with a hydrogen abundance (in mass fraction) X = 0.7374, and metallicity Z = 0.0134. We consider the number density of hydrogen (range: 6 ≤ log(n_〈H〉) ≤ 20 in steps of 0.1 dex) and the temperature (range: 1.7 ≤ log(T/K) ≤ 4.0 in steps of 0.01 dex) as independent variables.

The gas opacity provides a minimum base opacity, mostly provided by water vapour line opacity, which is reached when all dust has sublimated. The inclusion of the gas opacity means that there is an upper limit for the thermostat regulation mechanism; see Sect. 4.3.

3 Results

We start the presentation of our modelling results with Fig. 2, which shows some selected disc properties along the midplane, before we turn to the 2D vertical structure in Fig. 3. Figure 2 shows the gas and dust column densities in our models as evolutionary sequence. The earlier snapshots are similar to the first one depicted after 0.03 Myr. The disc mass increases from about 0.005 M_⊙ after 1000 yr to about 0.26 M_⊙ after 1 Myr, before it starts to decrease again; see Table 1. The radial disc extension, as measured by the radius that produces 90% of the continuum flux at 1.3 mm, according to the ProDiMo dust radiative transfer model, increases from about 5 to 130 au.

The dust column densities are equal to 0.004 × Σ(r) by assumption in the ProDiMo models, but Fig. 2 shows them after dust sublimation and ice formation. In the early phases, there is an inner region (≲0.5 au) where the silicates have sublimated in the midplane regions (only there), leading to a significant decrease of the dust column densities. One can also see a smaller, sudden increase in the dust-to-gas ratio, which is due to water ice formation. The snowline is situated at about 3 au after 0.03 Myr, 4 au after 0.1 Myr, 6 au after 0.3 Myr, 3 au after 1 Myr, and about 1.5 au after 3 Myr.

The middle column of plots in Fig. 2 shows the calculated midplane temperatures. While the ProDiMo models (dashed cyan lines), which assume fully condensed refractory dust, result in maximum midplane temperatures as high as 4000 K, the GGchem models (after dust sublimation) shows a temperature limit of about 1500–1700 Κ that is never exceeded, indicated by the black thin dashed line: the thermostat mechanism. Fitting the GGchem results in the inner regions, where the silicates have sublimated and this limiting temperature is reached, we find

(19)

where is the scale height derived from the midplane sound speed c_T. In the following, we denote T_sil as the silicate stability temperature. The physical idea behind the fit in Eq. (19) is that the sublimation temperature increases with pressure, and Σ/Η_ρ is roughly proportional to the midplane gas pressure. The silicate stability temperature is reached inside of about 0.4 au after 0.03 Myr, 0.5 au after 0.1 Myr, and 0.5 au after 0.3 Myr. At later evolutionary stages, the viscous heating of the disc subsides, the disc midplane becomes much cooler, and the silicate stability temperature is not reached anymore. Beyond the snowline, the GGchem midplane temperatures (after ice formation) are slightly higher than in the ProDiMo models, because the ProDiMo models do not include ice opacities. Consequently, the vertical optical depths are larger beyond the snowline in the GGchem models, and so the viscous heating has a stronger impact on the midplane temperature.

The right column in Fig. 2 shows the computed material composition of the solid particles in terms of volume fractions in the disc midplane. The materials are only plotted when they exceed 8% of the total solid volume anywhere in the midplane. We see a repetitive pattern of Ca-Al-Ti compounds in the innermost hot regions, followed by silicates with a complex composition. The most important Ca-Al-Ti oxides are corundum (A1₂O₃), gehlenite (Ca₂Al₂Si0₇), spinel (MgAl₂O₄), and perovskite (CaTiO₃). There is also baddeleyite (ZrO₂). The most abundant silicate species are fosterite (Mg₂SiO₄), enstatite (MgSiO₃), the feldspar end-members anthorite (CaAl₂Si₂O₈) and albite (NaAlSi₃O₈), mixed with either some solid iron at higher temperatures, or fayalite (Fe₂SiO₄) at lower temperatures. As the midplane temperature decreases to 650 K, sulphur starts to condense as well, mostly in form of troilite (FeS), before a number of phyllosilicates replace the silicates below about 350 K, step by step, most prominently lizardite (Mg₃Si₂O₉H₄). Below about 150 Κ water ice freezes out, and below about 85 Κ ammonia ice freezes out in addition.

Figure 3 shows some 2D results of our disc model after modelling stage 3 for an early and a late evolutionary phase. After 0.1 Myr, the strong viscous heating leads to a highly puffed-up inner region which puts the outer disc into a shadow. At later evolutionary phases (1 Myr) the inner disc flattens and we have a more textbook-like flaring disc. The T_d(r, z) structure can be subdivided into three distinct vertical layers: (1) the optically thin zone which is directly illuminated by the star, where we find vertically constant temperatures , (2) a transition zone that is radially optically thick but vertically optically thin, where the temperature drops quickly as we enter the disc shadow in the downward direction, and (3) the optically thick midplane region where viscous heating leads to another temperature increase.

Figure 3 also shows the local dust-to-gas ratio after modelling stage 3, which does not include dust settling. The generic value of 0.01 is only reached in the cold icy regions. The input value 0.004 is present in most other disc regions. However, the dust-to-gas ratio drops to much lower values of order 10^−3.5–10^−4.5 in the core of the inner disc, where the silicates have sublimated; the viscous heating burns a “hole” into the midplane from the left, very similar to Fig. 1 in Min et al. (2011). However, Min et al. only considered one dust species. In our models, the dust never sublimates completely, and the Al-Ca-Ti oxides remain.

Another interesting result of our models is that the stability of the silicates is not guaranteed in the inner disc high above the midplane. Since solid stability depends not only on temperature, but also on pressure, and the pressure drops by orders of magnitude high above the midplane, first the silicates and then even the Al-Ca-Ti oxides sublimate, and the inner disc rim takes on a rounded shape. Again, these results are very similar to Fig. 1 in Min et al. (2011), but more smooth as we include more than just one solid species. We mark the vertical Rosseland optical depth τ_ross with contour lines for values 1 and 10 in our plots, where the latter marks the height used as boundary condition for the stage 3 modelling procedure with GGchem as explained in Sect. 2.3.

The bottom row of plots in Fig. 3 give an impression of the spatial distribution of the different grain materials in the disc. We show six simplified categories of grain materials ordered by temperature: (1) the Al-Ca-Ti oxides, (2) ordinary silicates, (3) sulphur-containing silicates, (4) phyllosilicates, (5) water ice, and (6) ammonia ice. Each category includes materials from the previous category. Examples of the complex material composition from different points of a selected disc model are listed in Table 2.

The absorption opacities of these categories of condensates are shown in Fig. 4. We see that (1) the ordinary silicates have about one order of magnitude larger opacities than the Al-Ca-Ti oxides, (2) the inclusion of the conductive FeS into the grains leads to a substantial increase in the absorption at near-IR and mm wavelengths, and (3) ice formation leads to another substantial increase in the dust absorption opacities per gas mass, in particular in the UV and the far-IR, with the occurrence of additional near-IR, mid-IR, and far-IR ice features.

Figure 5 shows the evolution of the spectral energy distribution (SED) over time of our disc model up to modelling stage 2. The maximum luminosity (L_* = 3.7 L_⊙, L_acc = 0.89 L⊙) is reached after about 300 000 yr (black line), when the disc is featured by a tall ‘puffed-up’ inner disc. The strong internal viscous heating in this evolutionary phase causes a strong temperature inversion in the inner midplane, leading to high vertical extensions of the disc inside of about 1 au. This puffed-up inner disc (not an puffed-up inner rim) intersects a lot of starlight, which is then converted into near-IR and mid-IR dust emission. In contrast, the outer disc is cold, being situated in the shadow of the tall inner disc, which leads to relatively little far-IR emission and an overall steep slope of the SED. Once the mass accretion rate diminishes, and the viscous heating of the inner disc subsides, the inner midplane regions get a lot cooler, and the scale heights decrease by a factor of about two. The increase in the stellar mass with time also contributes to this evolution of the scale height, yet with less significance. Therefore, in the later evolutionary stages, we find a flaring disc, and a slightly warmer outer disc. The flat inner disc intersects only little starlight, leading to much less near-IR emission up to 10 μm. The shape of the SED in later phases resembles transitional discs with a reincrease in the flux beyond 20 μm.

Fig. 2 Evolution of gas and dust column densities (left column) and midplane temperatures (middle column) with ages indicated on the left. The green dots on the left mark the popular MMSN-value of 1700 g cm−2 at 1 au. The dashed cyan lines are the resulting midplane temperatures from the respective ProDiMo models, and the thick magenta lines are the midplane temperatures after dust sublimation. The dotted black line indicates the silicate stability temperature; see main text. The right column shows the material composition of the solid particles in the midplane.

Fig. 3 Selected 2D results of disc models after stage 3 as function of radius r and relative height over the midplane z/r. Upper row: hydrogen nuclei density n〈H〉(r, z). Second row: dust temperature structure Td(r, z). Third row: local dust-to-gas ratio. Bottom row: classes of solid mixtures(r, z); see Table 2. Additional dashed contour lines mark the vertical Rosseland optical depth τRoss = 1 and τRoss = 100. The black areas in plots in the two lower rows marks the end of the GGchem modelling domain n〈H〉 > 107 cm−3.

Fig. 4 Calculated dust absorption opacities per gas mass for the solid mixtures typically found in our stage 3 GGchem models; see examples in Table 2. The two red stars represent dust absorption opacities commonly used to derive disc masses from mm fluxes: 10 cm2/g(dust) at 1000 GHz (Beckwith et al. 1990) scaled to 3.5 cm2/g (dust) at 850 μm.

Fig. 5 Spectral energy distributions as function of age calculated by ProDiMo (after modelling stage 2).

4 Discussion

After having established how the disc temperature structure evolves with time, and which minerals are thermodynamically stable as function of radius and time, we performed a timescale analysis in order to explore how the grains are expected to move, grow, and change internal structure and material composition. We discuss this analysis below. The objective is to find a natural pathway to explain the occurrence of mm-sized pure Al-Ca-rich grains embedded in a matrix composed of several amorphous silicate-rich materials.

4.1 Annealing timescales

Annealing is the process of internal rearrangement of a lattice structure by solid diffusion, where atoms (or groups of atoms) hop between neighbouring lattice places, which allows to purify the material and eliminate imperfections in the lattice. The timescale for annealing according to Gail & Sedlmayr (1999) is given by

(20)

where ℓ is the total distance travelled by the atoms via a random walk, v is a typical oscillation frequency of the lattice, and λ is the spatial step-size, i.e. the distance between two neighbouring lattice places. E_b is the energy barrier involved in moving the group of atoms from on lattice place to the next. For a solid layer of thickness ℓ to develop a well-ordered crystalline structure, as observed for the CAIs, the annealing timescale must be short compared to the residence time t, i.e. τ_anneal ≪ t. Since we are interested here in the question whether or not a particle of size a can become crystalline throughout, we use a = ℓ.

For silicates, Gail & Sedlmayr (1999) estimated ν ≈ 2 × 10¹³ Hz from an average vibrational frequency of SiO₄ tetrahedrons, and assumed E_b/k ≈ 41000 Κ based on annealing experiment of Nuth & Donn (1982), later confirmed by experiments on silicate annealing by Hallenbeck et al. (1998).

We consider Al₂O₃ as a representative material for the Al-Ca-Ti oxides. We estimate ν ≈ 2.4 × 10¹³ Hz and E_b/k ≈ 48 000 Κ from the general observation that annealing temperatures roughly scale with sublimation temperatures. In a solar composition gas at 1 bar, the sublimation temperatures of A1₂O₃ and Mg₂SiO₄ have been found to be 1960 and 1660 K, respectively (Woitke et al. 2018), suggesting a correction factor of 1.18. A1₂O₃ is known to have a hexagonal rhomboidal lattice structure with lattice constants 4.76 and 12.99Å, so the annealing becomes dependent on direction, also because E_b is likely to depend on direction. However, here we simply assume λ ≈ 1 nm. The resulting annealing timescales are shown in Fig. 6 for different particle sizes a, and will be discussed together with the coagulation timescales at the end of the next section.

4.2 Coagulation timescales

A particle of size α will undergo π (a + a′)²Δν n_〈H〉 f(a′) da′ collisions per second with particles of sizes a′…a′+da′, where f(a) is the dust size distribution function [cm⁻¹/H-nucleus]. Δv(a, a′) is the turbulence-induced relative velocity between particles of sizes a and a′ (Ormel & Cuzzi 2007) and π(a + a′)² is the geometrical cross section. We define the timescale of coagulation as the time required to increase a particle’s volume (or mass) by a factor of exp(l)

(21)

where V = (4π/3) a³ is the particle volume and (4π/3) a′³ is the increase in the particle’s volume by one sticking collision. α_S is the sticking efficiency. We use min{1, a′³/a³} in the integration kernel instead of just a′³ /a³ to avoid interpreting the case dV > V as super fast growth. Instead, in the case of small grains colliding with big grains (a < a′), we use just their collisional frequency.

To calculate the turbulence-induced relative velocities Δv(a, a′), we use the semiempirical formulae of Ormel & Cuzzi (2007), which depend on the local sound speed, the turbulence parameter α, the two Stokes numbers of the colliding particles, and the gas Reynolds number. At gas densities of about 10¹⁵ cm⁻³, typical for these midplane regions, the Stokes numbers of mm-sized particles are of order St ≈ 10⁻⁴, so we mostly have the limiting case of small Stokes numbers (small particles) as described by Eq. (26) in Ormel & Cuzzi (2007), where Δv ≈ C_T (3 α St)^1/2; see also Birnstiel (2023). In that case, close inspection of Eq. (21) shows that for all sizes a, the integral on the right side is dominated by the largest collisional partners as these collisions have the largest Δv. For small particles a, the coagulation timescale turns out to be size-independent, whereas for the largest particles, it is slightly longer due to the lack of even larger particles.

In the close midplane regions, where the silicates sublimate, we find a dust-to-gas mass ratio of approximately 2 × 10⁻⁴, corresponding to a size-reduction factor of γ ≈ 0.3 (see Eq. (13)). Therefore, according to our assumption a_max = 3 mm, the largest Al-Ca-Ti oxide particles are initially about 1 mm in size after silicate sublimation, and with α = 10⁻³, we find relative velocities Δv between any particles and mm-sized particles between a few 10 cm s⁻¹ to a few m s⁻¹.

These relative velocities are already close to the fragmentation threshold velocity (Blum 2018; Birnstiel 2023). Furthermore, for large particles colliding with large particles, the sticking efficiency is likely to be much smaller than α_S = 1 because of the bouncing barrier (e.g. Güttler et al. 2010). When mm-sized particles collide with mm-sized or even larger particles, the contact (van der Waal’s) forces during the collisions become negligible in comparison to their inertia, and therefore, such particles simply do not stick together when they collide, they just bounce. For simplicity, we calculate the integral in Eq. (21) anyway up to a_max (few mms), with α_S = 1, to get an order of magnitude estimation of the coagulation timescale, keeping in mind that the actual timescales can be significantly longer because of the fragmentation threshold and the bouncing barrier.

Our timescale analysis is shown in Fig. 6 for the disc model after 0.1 Myr. We compare our coagulation timescale (Eq. (21)) with the simplified growth timescale proposed by Birnstiel (2023)

(22)

where Z is the dust-to-gas mass ratio and Ω is the Keplerian rotation frequency, showing good agreement. Thus, creating mm-sized particles by coagulation is a fast process taking only about 100yr. The growth by coagulation will stop when the particles have reached the fragmentation barrier. Once the particles have reached this barrier, fragmentation replenishes the small particles (Blum & Wurm 2008; Birnstiel 2023), and an equilibrium between coagulation and fragmentation establishes a power-law size distribution up to the fragmentation barrier.

In the close midplane, where temperatures are regulated to remain close to the silicate sublimation temperatures ≈1500–1700 K (Eq. (19)), annealing a 1 mm Al-Ca-Ti oxide particle takes about 10⁴ yr, which is still short compared to the lifetime ~10⁵ yr of the T-regulated midplane zone, where the silicates are predominantly in the gas phase. However, already beyond 0.3 au, where the temperatures are just slightly lower, 1400–1500 K, the annealing timescales increase significantly. However, it is still possible to form 100 μm annealed Al-Ca-Ti oxide particles there, which can then coagulate quickly.

The situation changes substantially outside of about 0.45 au in this model after 0.1 Myr, where the silicates are thermally stable and the thermostat mechanism breaks down. The temperature drops quickly, the dust-to-gas ratio increases by almost 2 orders of magnitude, the coagulation timescales (Eq. (22)) become even shorter, and the annealing timescales become huge. In this case, we expect the silicates to resublimate on the existing Al-Ca-Ti oxide surfaces to form an amorphous mantle mainly composed of Mg₂SiO₄ and MgSiO₃, without having the time to form an ordered crystal structure. Other materials like feldspar-components, iron and nickel, and later FeS, will be incorporated in this mantel; see Table 2. These core-mantle grains will continue to collide, fragment and coagulate. The expected physical properties of this mantel correspond well to the ‘matrix’ as observed in chondritic materials. Laboratory analysis has shown that the matrix contains olivines, pyroxenes, Fe-Ni metals, lesser amounts of sulphides, sulphates, carbonates, and in some cases significant amounts of phyllosilicates (Brearley 1989; Buseck & Hua 1993; Scott & Krot 2005).

Fig. 6 Annealing and coagulation timescales for Al-Ca-Ti oxide particles of different sizes in the disc midplane as a function of r for an age of the solar nebula of 0.1 Myr. The silicates have sublimated inside of about 0.45 au. Inside of about 0.2 au, only Al2O3 and ZrO2 remain, whereas outside of 0.2 au, Ca2Al2SiO7 and CaTiO3 are present as well, which creates another small step in the dust-to-gas ratio and hence in the coagulation timescales.

4.3 Range of validity of the thermostat mechanism

As discussed in Sect. 2.3.3, the thermostat regulation mechanism enforces the midplane temperature to stay close to the silicate sublimation temperature (Eq. (19)), provided that the viscous heating is strong enough and the disc is optically thick enough. However, there is a limit of this regulation mechanism, namely when the dust opacity becomes smaller than the gas opacity. In that case, a further increase in the local temperature does not result in a significant lowering of the total Rosseland opacity (Eq. (18)), and hence the thermostat breaks down.

In order to study the range of validity of the thermostat mechanism, we calculated a series of stage-3 (GGchem) models where we artificially modified the viscous heating rate F_vis (Eq. (7)) by a factor ranging from 0.005 to 200; see Fig. 7. We selected a representative column at r = 0.2 au from our disc model after 0.1 Myr for this investigation. Since F_vis is proportional to the mass accretion rate, we plot the results over Ṁ_acc. The nominal values of the unmodified model are Ṁ_acc = 1.04 × 10⁻⁶ M_⊙ yr⁻¹ and F_vis = 4.14 × 10⁶ erg cm⁻², which corresponds to the centre of the plot, where we reach T ≈ 1600 K in the midplane.

As we increase Ṁ_acc from very low values, we reach the silicate sublimation temperature at about 7 × 10⁻⁸ M_⊙ yr⁻¹, beyond which the dust opacity decreases rapidly, until only the Al-Ca-Ti oxides remain in the midplane. This dust opacity, however, is still larger than the gas opacity by about one order of magnitude. Increasing Ṁ_acc further, to a critical value of about 4 × 10⁻⁶ M_⊙ yr⁻¹, has little effect on the midplane temperature; this is the range where the thermostat mechanism is at work. However, for even larger Ṁ_acc, the Al-Ca-Ti oxides eventually start to sublimate, too, and once their opacity becomes negligible, the gas opacity takes over the regulation of the further increase in the midplane temperature with increasing Ṁ_acc.

Therefore, we expect the Al-Ca-Ti oxides to sublimate when the mass accretion rate exceeds the order of ≈10⁻⁶ to 10⁻⁵ M_⊙ yr⁻¹. The exact value of that critical mass accretion rate will depend on the surface density, which is about 10⁴g cm⁻² in our model at 0.2 au after 0.1 Myr. If the surface density is smaller, the optical depths are smaller, leading to less pronounced temperature inversions, and the critical mass accretion rate will be higher.

Fig. 7 Range of validity of the thermostat regulation mechanism. The upper plot shows the calculated midplane temperatures in our disc model at r = 0.2 au after 0.1 Myr for a sequence of simulations where we modified the mass-accretion rate Ṁacc from its nominal value of 1.04 × 10−6 M⊙ yr−1. The dashed line shows the resulting midplane temperatures when we use the fixed dust opacities from the ProDiMo model. The lower plot shows the corresponding dust and gas Rosseland opacities in the disc midplane, taking into account dust sublimation.

4.4 Shortcomings of the model, and future improvements

Although our model gives a clear picture of how dust grains move in the early solar nebula, what their thermal history is, which we determined using detailed radiative transfer calculations, and how the different minerals sublimate and resublimate as a function of time, we are far from a fully consistent, time-dependent 3D treatment of the CAI formation problem. In each modelling stage (Fig. 1), we made slightly different assumptions about dust opacities and temperatures that could have an impact on the previous modelling stages. For example, in the second modelling stage (ProDiMo) we assume a constant dust material mixture with a settled size distribution between constant minimum and maximum grain sizes, whereas in the third modelling stage (GGchem) we assume a varying material composition and an unsettled powerlaw size distribution with adjusted maximum and minimum size according to the resulting dust/gas mass ratio. However, these differences only play a role in the optically thick hot core of the disc where the silicates can sublimate. We checked that at the boundary τ_Ross ~ 10, below which the GGchem models are applied, the Rosseland opacities in both models agree with each other to within a factor of two. An implementation of the dust opacity including dust sublimation in the phase-1 hydro-simulatiomns might be feasible in future applications, but goes beyond the scope of this paper.

The midplane temperatures (see Fig. 2) can indeed change substantially when silicate sublimation is taken into account, but this is relevant only in the optically thick midplane (τ_Ross > 10). The temperature deviations in those deep layers, however, do not affect the upper layers, as the upper layers only ‘see’ the boundary temperature, which remains unchanged in modelling phase 3. There is one exception though, and that is the vertical hydrostatic structure, which is more extended when the mid-plane is warmer. Here, when we compute the hydrostatic vertical structure in ProDiMo, we implemented a cap of the local dust temperature by the silicate sublimation temperature (Eq. (19)) to make sure we do not overpredict the vertical scale heights.

In addition we note that all dust stability computations with GGchem are time-independent by design. They are based on the principle of minimisation of the total Gibbs free energy, so we assume that all internal rearrangement processes necessary to reach the thermodynamically most favourable solid state happen in due time, which becomes increasingly questionable when the temperatures become low; see Herbort et al. (2020, their Fig. 11). However, close to the stability temperatures of all minerals, sublimation and resublimation are actually quite fast processes, and we have the annealing timescales to discuss the rearrangements. In fact, we do not follow the grains in a time-dependent way and calculate how their position, size, material composition and internal structure change, but we base all our conclusions on the trajectories found in the first modelling phase, and the condensation, growth, and annealing timescales found in third modelling phase, assuming that an equilibrium between dust coagulation and fragmentation holds and always installs a powerlaw size distribution quickly. A more detailed treatment of these processes in the frame of the phase-1 hydrodynamical disc evolution model would be desirable.

5 Summary and conclusions

This paper proposes a pathway to form mm-sized, pure, and annealed Al-Ca-Ti oxide grains embedded in a matrix of amorphous silicates in the earliest evolutionary phases of the solar nebula. We associate these particles with the Ca-Al-rich inclusions (CAIs) found in many chondritic meteorites.

Our model combines (1) the 1D viscous disc evolutionary models of Drążkowska & Dullemond (2018, 2023) with (2) the 2D radiation thermo-chemical disc models of Woitke et al. (2024), (3) the equilibrium condensation models of Woitke et al. (2018), and (4) new dust opacity calculations.

The models reveal a thermostat mechanism in the disc midplane regions close to the star (Min et al. 2011), which keeps the local temperature stable at about 1500–1700 K, just above the silicate sublimation temperature. All dust grains inherited from the molecular cloud, which temporarily reside in this region, will loose their silicates and all other more volatile components quickly, making the Al-Ca-Ti oxides the only condensates that remain in this region. The thermostat mechanism works by lowering the dust opacity by silicate sublimation, until the viscous heat accumulating in the midplane can escape vertically.

These particular, regulated temperature conditions remain relatively stable for hundreds of thousands of years, which allows the Al-Ca-Ti particles to grow by coagulation and create a semiordered crystal structure by annealing. The growth of these particles by coagulation only stops once the size of the Al-Ca-Ti particles reaches the fragmentation barrier of a few mm, creating the size distribution that we can still observe today by measuring the sizes of the CAIs in chondritic materials; see Nakamura et al. (2007) for example.

To form the CAIs, that is, to embed these particles in a silicate matrix, it is essential to have an outward motion of the dust grains during the earliest evolutionary phases - as is the case in the Drążkowska & Dullemond (2018, 2023) models -, which is sketched as scenario 2 on the right side of Fig. 8. The Al-Ca-Ti particles are dragged out with the viscously expanding gas to radii beyond 0.5 au, where the thermostat regulation mechanism breaks down, and the silicates resublimate on the existing surfaces of the Al-Ca-Ti oxide particles. When this happens, temperatures are already too low for annealing, and so this creates an amorphous silicate mantel. In contrast, according to the more common scenario 1 in Fig. 8, where all matter is accreted outside-in through the disc, at no time or space in the disc would one expect to find large, pure, and crystalline Al-Ca-Ti rich particles embedded in an amorphous silicate matrix.

Our scenario for CAI formation in the early Solar System therefore requires two ingredients: (1) the existence of an inner disc zone where the temperature is regulated to stay just above the silicate sublimation temperature, and (2) an outward flow of solid particles through this zone. According to our disc model, both of these requirements are met, as shown in Fig. 9 by the cyan overlap region, up to a disc age of about 50 000 yr. Eventually, after the mass-accretion rate has fallen to values of $ 10″⁷ M_e yr″¹ and the viscous heating of the disc has subsided, the thermostat mechanism breaks down, and the dust particles from the molecular cloud that continue to join the disc do so mostly unaltered, that is, without annealing and without losing their silicates of other refractory components on the way.

Fig. 8 Two scenarios for the early evolution of the solar nebula. In the common scenario 1, the disc is fed by the molecular cloud at relatively large radial distances, and the gas and the dust grains move inwards. In our favoured scenario 2, the disc is fed at much smaller radii, and the dust grains are dragged outward with the viscously spreading gas, which provides a natural explanation for the formation of CAIs embedded in a silicate matrix.

Fig. 9 CAI formation in time and space. The critical radius (black line) divides this diagram into a region where millimetre grains move inward (due to inward gas motion and radial drift) and a region where these grains are dragged outward with the expanding gas. The light red line encircles the hot T-regulated region with thermally unstable silicates, leading to pure Al-Ca-Ti oxide particles. The cyan shaded area is the region of lasting CAI-formation. The Al-Ca-Ti oxide particles formed in this region will eventually reside in a silicate matrix and populate the disc. Four trajectories for different particle sizes are plotted with green lines.

Acknowledgements

P.W. acknowledges funding from the European Union H2020-MSCA-ITN-2019 under Grant Agreement no. 860470 (CHAMELEON). J.D. acknowledges funding from the European Union under the European Union’s Horizon Europe Research & Innovation Programme 101040037 (PLAN-ETOIDS). P.M. acknowledges funding from the Italian Mnisterial Grant PRIN 2022, “Radiative opacities for astrophysical applications”, no. 2022NEXMP8, CUP C53D23001220006.

Appendix A Optical constants

Table A.1 lists the names, chemical formula and references for optical constants for the solid materials that occur frequently and abundantly in our disc model.

Electronic files containing the optical constants (real and imaginary parts of the refractory index (n, k)) can be collected from https://zenodo.org/records/8221362 and https://github.com/cdominik/optool.

Appendix B Phasediagrams for Al₂O₃ and SiO₂

In order to discuss the availability of liquid phases in a protoplanetary disc, we performed simple GGchem simulations with only 4 elements (H, O, Si, Al) for hydrogen-rich and hydrogen-poor environments. For the hydrogen-rich model, we used solar abundances (Asplund et al. 2009): ϵ_H = 12, ϵ_O = 8.69, ϵ_Al = 6.45, ϵ_Si = 7.51. For the hydrogen-poor model we used ϵ_H = 1 instead while leaving the other element abundances unchanged. The resulting phasediagrams are shown in Fig. B.1. The coloured areas show the regions in the (p, T)-plane where a certain molecule, solid [s] or liquid [l] contains most of an element. We read off the triplepoints for quartz in hydrogen-rich and hydrogen-poor environments as (8 bar, 1696 K) and (1.5 × 10⁻⁷ bar, 1696 K), respectively. For corundum, we find (7000 bar, 2327 K) and (2.5 × 10⁻⁵ bar, 2327 K), respectively.

In the H-rich case, H₂ is the dominating molecule, providing the gas pressure, and most of the oxygen is locked in H₂O. Both effects reduce the stability and amounts of Si-oxide and Al-oxide molecules in the gas phase required for the condensation of quartz and corundum, and hence large gas pressures are required for condensation. In the H-poor case, the most abundant molecule is O₂, providing the gas pressure, and Si-oxide and Al-oxide molecules form easier and in larger quantities. In that case, condensates can form already at much lower gas pressures – the difference is about 8 orders of magnitude in pressure. This finding has been described as “chemical sputtering” by Gail & Sedlmayr (1999): At the surface of a hot oxide grain, the collisions with H₂-molecules and H₂O tend to destabilise and decompose the refractory materials; oxides are less stable in a reducing atmosphere. At millibar pressures, the solids decompose in the presence of H₂ and sublimate before the melting point is reached.

The H-poor case resembles most situations studied in laboratory experiments or observations on Earth, where it is possible to create a melt of a refractory material by heating. However, gas pressures as high as 8 bar and 7000 bar (for quartz and corundum, respectively) are required for the existence of liquid phases in a solar composition gas. Such pressures are not available in protoplanetary discs, where the gas pressures barely reach a few millibars in the closest midplane during disc evolution. Therefore, our conclusion is that liquid phases are not available during disc evolution. A true melting of oxide materials only seems possible in the interior of gravitationally bound bodies such as planetesimals, where such pressures can be reached, and where the melt is not in direct contact with H₂.

Fig. B.1 Phase diagrams for quartz (SiO2) and corundum (Al2O3) in H-rich (left) and H-poor (right) environments. The black circles mark the triple points of quartz and corundum in these cases.

References

All Tables

All Figures

	Fig. 1 Three modelling stages, passed quantities, and results.
In the text

Fig. 2 Evolution of gas and dust column densities (left column) and midplane temperatures (middle column) with ages indicated on the left. The green dots on the left mark the popular MMSN-value of 1700 g cm−2 at 1 au. The dashed cyan lines are the resulting midplane temperatures from the respective ProDiMo models, and the thick magenta lines are the midplane temperatures after dust sublimation. The dotted black line indicates the silicate stability temperature; see main text. The right column shows the material composition of the solid particles in the midplane.

In the text

Fig. 3 Selected 2D results of disc models after stage 3 as function of radius r and relative height over the midplane z/r. Upper row: hydrogen nuclei density n〈H〉(r, z). Second row: dust temperature structure Td(r, z). Third row: local dust-to-gas ratio. Bottom row: classes of solid mixtures(r, z); see Table 2. Additional dashed contour lines mark the vertical Rosseland optical depth τRoss = 1 and τRoss = 100. The black areas in plots in the two lower rows marks the end of the GGchem modelling domain n〈H〉 > 107 cm−3.

In the text

	Fig. 4 Calculated dust absorption opacities per gas mass for the solid mixtures typically found in our stage 3 GGchem models; see examples in Table 2. The two red stars represent dust absorption opacities commonly used to derive disc masses from mm fluxes: 10 cm2/g(dust) at 1000 GHz (Beckwith et al. 1990) scaled to 3.5 cm2/g (dust) at 850 μm.
In the text

	Fig. 5 Spectral energy distributions as function of age calculated by ProDiMo (after modelling stage 2).
In the text

Fig. 6 Annealing and coagulation timescales for Al-Ca-Ti oxide particles of different sizes in the disc midplane as a function of r for an age of the solar nebula of 0.1 Myr. The silicates have sublimated inside of about 0.45 au. Inside of about 0.2 au, only Al2O3 and ZrO2 remain, whereas outside of 0.2 au, Ca2Al2SiO7 and CaTiO3 are present as well, which creates another small step in the dust-to-gas ratio and hence in the coagulation timescales.

In the text

Fig. 7 Range of validity of the thermostat regulation mechanism. The upper plot shows the calculated midplane temperatures in our disc model at r = 0.2 au after 0.1 Myr for a sequence of simulations where we modified the mass-accretion rate Ṁacc from its nominal value of 1.04 × 10−6 M⊙ yr−1. The dashed line shows the resulting midplane temperatures when we use the fixed dust opacities from the ProDiMo model. The lower plot shows the corresponding dust and gas Rosseland opacities in the disc midplane, taking into account dust sublimation.

In the text

Fig. 8 Two scenarios for the early evolution of the solar nebula. In the common scenario 1, the disc is fed by the molecular cloud at relatively large radial distances, and the gas and the dust grains move inwards. In our favoured scenario 2, the disc is fed at much smaller radii, and the dust grains are dragged outward with the viscously spreading gas, which provides a natural explanation for the formation of CAIs embedded in a silicate matrix.

In the text

Fig. 9 CAI formation in time and space. The critical radius (black line) divides this diagram into a region where millimetre grains move inward (due to inward gas motion and radial drift) and a region where these grains are dragged outward with the expanding gas. The light red line encircles the hot T-regulated region with thermally unstable silicates, leading to pure Al-Ca-Ti oxide particles. The cyan shaded area is the region of lasting CAI-formation. The Al-Ca-Ti oxide particles formed in this region will eventually reside in a silicate matrix and populate the disc. Four trajectories for different particle sizes are plotted with green lines.

In the text

	Fig. B.1 Phase diagrams for quartz (SiO2) and corundum (Al2O3) in H-rich (left) and H-poor (right) environments. The black circles mark the triple points of quartz and corundum in these cases.
In the text