CAI formation in the early Solar System
P. Woitke J. Drażkowska H. Lammer K. Kadam P. Marigo
 





Space Research Institute, Austrian Academy of Sciences,
           Schmiedlstr. 6, A-8042, Graz, Austria
           Max Planck Institute for Solar System Research,
           Justus-von-Liebig-Weg 3, 37077 Göttingen, Germany
           Department of Physics and Astronomy G. Galilei,
           University of Padova, Vicolo dell'Osservatorio 3,
           I-35122, Padova, Italy




Ca-Al-rich inclusions (CAIs) are the oldest dated solid materials in the solar system, found as light-coloured crystalline ingredients in carbonaceous chondrite meteorites. Their formation time is commonly associated with age zero of the Solar System. Yet, the physical and chemical processes that once led to the formation of these sub-millimetre to centimetre-sized mineral particles in the early solar nebula are still a matter of debate. This paper proposes a pathway to form such inclusions during the earliest phases of disc evolution. We combine 1D viscous disc evolutionary models with 2D radiative transfer, equilibrium condensation, and new dust opacity calculations. We show that the viscous heating associated with the high accretion rates in the earliest evolutionary phases causes the midplane inside of about 0.5 au to heat up to limiting temperatures of about 1500-1700 K, but no further. These high temperatures force all refractory material components of the inherited interstellar dust grains to sublimate – except for a few Al-Ca-Ti oxides such as Al_2O_3, Ca_2Al_2SiO_7, and CaTiO_3. This is a recurring and very stable result in all our simulations since these minerals form a natural thermostat. Once the Mg-Fe silicates are gone, the dust becomes more transparent and the heat is more efficiently transported to the disc surface, which prevents any further warm-up. This thermostat mechanism keeps these minerals above their annealing temperature for hundreds of thousands of years, allowing them to form large, pure and crystalline particles. These particles are dragged out by the viscously spreading disc, and once they reach a distance of about 0.5 au, the silicates re-condense on the surface of the Ca-Al-rich particles, adding an amorphous silicate matrix.  We estimate that this mechanism to produce CAIs works during the first 50 000 years of disc evolution. These particles then continue to move outward and populate the entire disc up to radii of about 50 au, before, eventually, the accretion rate subsides, the disc cools, and the particles start to drift inwards.



CAI formation in the early Solar System
    P. Woitke1
        J. Drążkowska2
        H. Lammer1
        K. Kadam1
        P. Marigo3

    Received 8 April 2024 / Accepted 23 April 2024
==============================================================================================








§ INTRODUCTION



New stars are born when the cold cores of molecular clouds collapse under the pull of their own gravity.

Due to the conservation of angular momentum during cloud collapse, a rapidly rotating disc is formed around the central protostar <cit.>.
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 <cit.>. 

Important information about the environmental conditions in the early Solar System, 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 <cit.>. 


The matrix is the background material into which the other constituents like the chondrules, the CAIs, and larger mineral grains are embedded <cit.>. 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 Myrs after the CAIs originated. Lesser amounts of sulphides, sulphates, carbonates, and several other minerals are also found, and some meteorites also contain significant amounts of phyllosilicates <cit.>. 

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 <cit.>. 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 <cit.>.  Four CAIs have been dated by using the Pb-Pb chronometer, which yields a Solar System age of 4567.30 ± 0.16 Myrs <cit.>.  Recently, <cit.> re-evaluated the age of the CAIs to 4568.7 Myrs through unified Al-Mg and Pb-Pb chronology.  


Chondrules are major constituents of most chondrites.  They are about millimetre-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 <cit.>, they crystallised in minutes to hours at temperatures of 1300-1800 K, prior to accretion into their parent bodies.  Chondrules can be separated in two types based on the Mg-number that is controlled by the oxygen fugacity of the chondrule-forming environment in the disc, where type I chondrules formed under more 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 <cit.>. 

All internal rearrangement and phase-change processes within the chondrules require high temperatures, in particular the melting.  However, refractory materials such as quartz and corundum require pressures over 8 bar and 7000 bar, respectively, to have a liquid phase in a solar composition gas, see Appendix <ref>.  At lower pressures, there is only sublimation and resublimation, but no melting.  Since such large gas pressures are not available during disc evolution, we suggest to interpret the observations of “molten” materials by annealing.  A true melting of refractory materials seems to be possible only inside of internally heated, larger 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 between ≈ 1-3 Myr after the origin of the CAIs, in a period when the proto-Sun accreted more slowly <cit.>.  <cit.> analysed the 22 youngest chondrules and determined their ages from the 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 Myrs. This lifetime estimate is in good agreement with observations of protoplanetary discs, which show lifetimes of about 1-10 Myrs <cit.>.



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 ^16O-rich reducing environment <cit.>. 
Recently, <cit.> 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 ^16O compared to ^17O by ∼(3-23), which is in agreement with what has been proposed as early generations of chondrules <cit.>. According to <cit.> the ^16O-enriched objects are most likely melted amoeboid olivine aggregates that escaped from the incorporation into the ^16O-poor precursor dust of the chondrules. These authors further suggest from their analysis of the Ryugu samples 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, <cit.> suggest that the transported objects have been most likely destroyed during aqueous alteration in the Ryugu parent body. 
Further results from the analysis of these chondritic components include nucleosynthetic 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 their internal structure since 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 <cit.>. 









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 <cit.>. 
To better understand these processes, 

we have developed a model that combines disc evolution with 2D radiative transfer, equilibrium condensation and new dust opacity calculations, see Sect. 2.


In Sects. 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.












§ THE MODEL


Our model for the early evolution of the solar system is generated in three modelling stages, see Fig. <ref>.
First, we run a one-dimensional disc evolution model developed by <cit.> for an initial cloud mass of 1.15 M_⊙, which results in a stellar mass of about 1 M_⊙ after 3 Myrs, see Sect. <ref>.
Second, selecting a sequence of ages, we import the disc structures from stage 1 into the 2D radiation thermo-chemical disc model ProDiMo <cit.> and run detailed radiative transfer calculations, see Sect. <ref>.
And third, we import the 2D disc structures from stage 2 into the chemical and phase equilibrium code GGchem <cit.> 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 <cit.>, 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. <ref>. 






 §.§ The evolutionary disc model







We use a hydrodynamic disc evolution model similar to the one described in <cit.>. 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^-15 s^-1. 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 <cit.>. 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, 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 Myrs.  


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^-3 <cit.>. During the infall, the compact disc vigorously expands and the outward-flowing gas drags millimetre-sized solid particles with it.  Such disc models have been proven 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 <cit.>, however, the formation of CAIs was not studied consistently in these models, instead these were added afterwards.






 §.§ The ProDiMo disc model




The Protoplanetary Disc Model (ProDiMo) developed by <cit.> 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 <cit.>, X-ray ionisation <cit.>, FUV and EUV photo-processes <cit.>, and ice chemistry. The ProDiMo models in this paper are based on the radial disc structures imported from the evolutionary disc model (Sect. <ref>) for a sequence of eight consecutive ages 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, and 3 Myrs, 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^2], and the vertically integrated viscous heating rate F_ vis(r) [erg/s/cm^2].   



  §.§.§ Stellar setup


The stellar mass M_⋆ and the mass accretion rate Ṁ_ acc are taken from the evolutionary disc model (Sect. <ref>). 
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 <cit.> for given age and M_⋆, see Table <ref>.
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 <cit.>.
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 <cit.>, 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–Boltz­mann law L_⋆=4π R_⋆^2 σ T_⋆^4, 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

    L_ acc = G M_⋆Ṁ_ acc/2 R_⋆(1-R_⋆/R_ in)  .

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 1500 K 

    R_ in = (0.1  au)
                 √(L_⋆+L_ acc/1.95 L_⊙) .

The value of 0.1 au for a total luminosity of 1.95 L_⊙ comes from previous experience with the ProDiMo disc models.  Equations (<ref>) and (<ref>) are solved iteratively until converged. Table <ref> summarises the stellar setup parameters for the ProDiMo models.  For comparison, the table also contains the viscous disc luminosity 

    L_ vis = ∫_R_ in^R_ out 
                      F_ vis(r) 2π r dr  ,

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.1 L_⋆ 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=10000K, representing the accretion luminosity L_ acc produced by the accretion shock.  The total spectral intensities emitted from the star's surface are calculated as

    I_ν^⋆ =  
      I_ν^ phot(T_⋆,log g) 
       + L_ acc/4π R_⋆^2 σ T_ acc^4 B_ν(T_ acc)  .

The photospheric intensities I_ν^ phot are interpolated from a grid of PHOENIX stellar atmosphere models by <cit.>, where the surface gravity is given by g=GM_⋆/R_⋆^2. 





  §.§.§ Disc structure, dust settling, and radiative transfer


Based on the input gas surface density structure Σ(r), ProDiMo computes an axisymmetric 2D density structure (r,z), where [cm^-3] is the hydrogen nuclei density, in vertical hydrostatic equilibrium as follows.  For each column r, starting with a guess of the temperature structure T(r,z), ProDiMo numerically integrates the equation of hydrostatic equilibrium upwards <cit.>

    1/ρdp/dz = - z GM_⋆/(r^2+z^2)^3/2 ,

using the “simplified hydrostatic model” <cit.>, where the isothermal sound speed c_T = √(p/ρ) = √(kT/(2.3  amu)) is evaluated with T= and a constant (H2-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  and the proper mean molecular weight are used to calculate the gas pressure, are visualised in Fig. 11 of <cit.>. Since we focus on the dust evolution in the midplane in this paper, where = 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 <cit.>.  We start with an unsettled size distribution function with a powerlaw f(a)∝a^- between a minimum () and a maximum dust size (). 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 , and the settling description of <cit.>, see <cit.> for details.  We used the following parameter values in this paper: =0.05 μm, =3mm, =3.5, dust-to-gas mass ratio =0.004, and =10^-3.  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 (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 <cit.>, see also <cit.> and Appendix A in <cit.>. 

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

    ÷F⃗ = Γ ,

where Γ is the local net heating rate converting non-radiative to radiative energy.








The input model (Sect. <ref>) provides the viscous heating rate per unit area F_ vis [erg/cm^2/s]. To calculate the heating rate per volume Γ, we assume that the heating occurs proportional to the gas density ρ as

    Γ(r,z)  = 1/2 F_ vis(r) ρ(r,z)/∫_ 0^∞ρ(r,z) dz ,

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

    F⃗(r,z) = -4 π/3(r,z)  J(r,z)

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

    (r,z) =
        ∫dB_ν(T)/dT dν  / ∫1/κ^ ext_ν(r,z)dB_ν(T)/dT dν

where κ^ ext_ν(r,z) is the dust extinction coefficient at position (r,z) in the disc, T=T_ dust(r,z) is the temperature, and B_ν(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 <cit.>.

Once the 2D disc radiative transfer problem is solved, and the dust temperature structure (r,z) is determined, ProDiMo calculates the chemical composition in time-independent mode and the gas temperature structure (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, (r,z) is required for the computations of the vertical hydrostatic equilibrium (Eq. <ref>), and therefore, ProDiMo performs global iterations between radiative transfer and hydrostatic equilibrium, which converges after about 10 iterations.





 §.§ 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 <cit.>. 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 K and 6000 K. For this paper, we use the code down to 50 K; for temperatures below 50 K we use T=50K.

The sole purpose of GGchem here is to predict the amount and volume composition of all solids, for given gas density , temperature T, and total element abundances prior to condensation, where we use the solar abundances of <cit.>. The GGchem setup in this paper is the same as in <cit.>. 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 have extended GGchem to calculate the opacity of the resulting mixture of solid materials as explained in Sect. <ref>.  However, that 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 T(r,z) in the optically thick midplane, and only there.



  §.§.§ The thermostat regulation mechanism



When the midplane is very optically thick and heated by frictional processes, the radiative flux is almost perfectly directed upwards in the vertical direction.  In this case, we can simplify Eqs. (<ref>) and (<ref>) as

    F(z)  = ∫_0^z Γ(z') dz'
            =  -4 σ/3(z) dT^4/dz ,

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. (<ref>) is written numerically as

    F(z_k-)    
       =  -4 σ/3√((z_k-1) (z_k)) T(z_k)^4-T(z_k-1)^4/z_k-z_k-1 ,

where z_k- = (z_k-1+z_k)/2. Equation (<ref>) 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 T(z_K) at that point as upper boundary condition, and the density structure (z) below (z≤z_K) from the ProDiMo disc model. Then, we



  1. guess the temperature T(z_k-1),

  2. compute the amount of dust and its material composition for given total element abundances, temperature T(z_k-1) and hydrogen nuclei particle density (z_k-1) with GGchem,

  3. calculate the dust size distribution function f(a) at z_k-1 (Sect. <ref>),

  4. calculate the frequency-dependent dust extinction coefficient κ^ ext_ν(z_k-1) according to Eq. (<ref>) by effective medium and Mie theory (Sect. <ref>), 

  5. calculate the Rosseland mean opacity (z_k-1) according to Eqs. (<ref>) and (<ref>),

  6. calculate T(z_k-1) according to Eq. (<ref>),

  7. go back to step 2 if not yet converged.
  
This way, we have implemented a stable thermostat regulation mechanism as discovered by <cit.>.  If the temperature T(z_k-1) becomes too high, the dust starts to sublimate, which lowers the Rosseland opacity , 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.



  §.§.§ 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^-1/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^3. The value 0.004 results from complete condensation of a solar composition gas <cit.> before ice formation, see <cit.> 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

    V_ ref =    ∫_0.05 μ m^3 mm4π/3a^3 f_ ref(a) da 
     =      4.53×10^-27  .

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

    γ = (V/V_ ref)^1/3 .

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. (<ref>), first and foremost because of the increased condensed volume per H-nucleus.



  §.§.§ 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 <ref>. We add 25% porosity. 

Using the effective medium theory of <cit.>, we calculate the effective optical constants of the mixed material as 

    (,) = F(V_s,n_s,k_s)

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

    κ^ ext_ν =     ∫_^ f(a)
                           π a^2 Q^ ext_ν(x,,) da
    
    κ^ abs_ν =     ∫_^ f(a)
                           π a^2 Q^ abs_ν(x,,) da 
    κ^ sca_ν =     ∫_^ f(a)
                           π a^2 Q^ sca_ν(x,,) da

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. <ref>. The extinction, absorption and scattering efficiencies are calculated according to Mie theory. We use the routine MIEX <cit.> in the implementation of Sebastian Wolff. The opacity computations are based on a collection of optical data by <cit.>, with additional data collected by Carsten Dominik and Michiel Min (<https://github.com/cdominik/optoo>). The complete selection of relevant materials and opacity sources is summarised in Table <ref>.



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

    = ^ dust(,V,V_s,T)  + ^ gas(,T)  .

The gas Rosseland mean opacity ^ gas is interpolated from 2D-tables provided by Paola Marigo, using the same opacity data sources and techniques explained in <cit.>, 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  <cit.> and  <cit.> 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_2/H_2, H_2/H, H_2/He, and H/He pairs.
To construct the Rosseland mean opacity tables, we integrated the  code <cit.> with the  code <cit.>. We adopted a solar mixture following <cit.>, 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() ≤ 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. <ref>.





§ RESULTS



We start the presentation of our modelling results with Fig. <ref>, which shows some selected disc properties along the midplane, before we turn to the 2D vertical structure in Fig. <ref>.  Figure <ref> 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 Myrs.

The disc mass increases from about 0.005 M_⊙ after 1000 yrs to about 0.26 M_⊙ after 1 Mys, before it starts to decrease again, see Table <ref>. 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 au to 130 au.  

The dust column densities are equal to 0.004×Σ(r) by assumption in the ProDiMo models, but Fig. <ref> shows them after dust sublimation and ice formation. In the early phases, there is an inner region (0.5au) 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 of the dust-to-gas ratio, which is due to water ice formation.  The snowline is situated at about 3 au after 0.03 Myrs, 4 au after 0.1 Myrs, 6 au after 0.3 Myrs, 3 au after 1 Myrs, and about 1.5 au after 3 Myrs.



The middle column of plots in Fig. <ref> 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 K 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 

    ln(T_ sil [K])  ≈ ln(1700) 
       +  0.12 ln(Σ [10^4 g/cm^2]/H_p [0.01 au])  ,

where H_p=(r^3 c_T^2/GM_⋆)^1/2 is the scale height derived from the midplane sound speed c_T. In the following, we will denote T_ sil as the silicate stability temperature.

The physical idea behind the fit in Eq. (<ref>) is that the sublimation temperature increases with pressure, and Σ/H_p is roughly proportional to the midplane gas pressure. 
The silicate stability temperature is reached inside of about 0.4 au after 0.03 Myrs, 0.5 au after 0.1 Myrs, and 0.5 au after 0.3 Myrs. 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.<ref> 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 (Al2O3), gehlenite (Ca2Al2SiO7), spinel (MgAl2O4), and perovskite (CaTiO3). There is also baddeleyite (ZrO2).  The most abundant silicate species are fosterite (Mg2SiO4), enstatite (MgSiO3), the feldspar end-members anthorite (CaAl2Si2O8) and albite (NaAlSi3O8), mixed with either some solid iron at higher temperatures, or fayalite (Fe2SiO4) 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 (Mg3Si2O9H4).  Below about 150 K water ice freezes out, and below about 85 K ammonia ice freezes out in addition.  

Figure <ref> shows some 2D-results of our disc model after modelling stage 3 for an early and a late evolutionary phase.  After 0.1 Myrs, 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 Myrs) the inner disc flattens and we have a more textbook-like flaring disc.  The (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 (r,z)∝√(r), (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 a temperature re-increase.    

Figure <ref> 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 to 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 <cit.>.  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 <cit.>, but more smooth as we include more than just one solid species.  We have marked 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. <ref>.

The bottom row of plots in Fig. <ref> 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 <ref>.



The absorption opacities of these categories of condensates are shown in Fig. <ref>. 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 of the absorption at near-IR and millimetre wavelengths, and (3) ice formation leads to another substantial increase of 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 <ref> 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 yrs (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 of 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 re-increase of the flux beyond 20 μm.






§ 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 will now perform a timescale analysis to discuss how the grains are expected to move, grow, change internal structure, and material composition. The objective is to find a natural pathway to explain the occurrence of millimetre-sized pure Al-Ca rich grains embedded in a matrix composed of several amorphous silicate-rich materials.



 §.§ 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 <cit.> is given by

    τ_ anneal^-1( ℓ )  = λ^2 ν exp(-E_b/kT)/3 ℓ^ 2 ,

where ℓ is the total distance travelled by the atoms via a random walk, ν 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, <cit.> estimated ν≈2×10^13Hz from an average vibrational frequency of SiO4 tetrahedrons, and assumed E_b/k≈ 41000K based on annealing experiment of <cit.>, later confirmed by experiments on silicate annealing by <cit.>.  

We consider Al2O3 as a representative material for the Al-Ca-Ti oxides.  We estimate ν≈2.4×10^13Hz and E_b/k≈ 48000K from the general observation that annealing temperatures roughly scale with sublimation temperatures. In a solar composition gas at 1 bar, the sublimation temperatures of Al2O3 and Mg2SiO4 have been found to be 1960 K and 1660 K, respectively <cit.>, suggesting a correction factor of 1.18. Al2O3 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 λ≈1nm. 
The resulting annealing timescales are shown in Fig. <ref> for different particle sizes a, and will be discussed together with the coagulation timescales at the end of the next section. 










 §.§ Coagulation timescales



A particle of size a will undergo π (a+a')^2Δ v  f(a') da' collisions per second with particles of sizes a'... a'+da', where f(a) is the dust size distribution function [cm^-1/H-nucleus]. Δ v(a,a') is the turbulence-induced relative velocity between particles of sizes a and a', see <cit.>, and π(a+a')^2 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(1)  

    τ_ coag^-1(a)  = dV/dt/V = ∫min{1,a'^3/a^3} π (a+a')^2 Δ v α_S f(a') da'  ,

where V=(4π/3) a^3 is the particle volume and (4π/3) a'^3 is the increase of the particle's volume by one sticking collision. α_S is the sticking efficiency. We use min{1,a'^3/a^3} in the integration kernel instead of just a'^3/a^3 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 semi-empirical formulae of <cit.>, 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^15 cm^-3, typical for these midplane regions, the Stokes numbers of mm-sized particles are of order ≈10^-4, so we mostly have the limiting case of small Stokes numbers (small particles) as described by equation (26) in <cit.>, where Δ v≈ c_T (3 α )^1/2, see also <cit.>. 
In that case, close inspection of Eq. (<ref>) 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^-4, corresponding to a size-reduction factor of γ≈0.3 (see Eq. <ref>).  Therefore, according to our assumption =3mm, the largest Al-Ca-Ti oxide particles are initially about 1 mm in size after silicate sublimation, and with α=10^-3, 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 <cit.>. 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 <cit.>. 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. (<ref>) anyway up to  (few millimetres), 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. <ref> for the disc model after 0.1 Myrs. We compare our coagulation timescale (Eq. <ref>) with the simplified growth timescale proposed by <cit.>

    τ_ growth^-1 = Z Ω ,

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 100 yrs. 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 <cit.>, 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-1700K (Eq. <ref>), annealing a 1 mm Al-Ca-Ti oxide particle takes about 10^4yrs, which is still short compared to the lifetime ∼10^5yrs 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-1500K, 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 Myrs, 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. <ref>) become even shorter, and the annealing timescales become huge.  In this case, we expect the silicates to re-condense on the existing Al-Ca-Ti oxide surfaces to form an amorphous mantle mainly composed of Mg2SiO4 and MgSiO3, 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 <ref>.  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 <cit.>.




 §.§ Range of validity of the thermostat mechanism



As discussed in Sect.<ref>, the thermostat regulation mechanism enforces the midplane temperature to stay close to the silicate sublimation temperature (Eq. <ref>), 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 of the local temperature does not result in a significant lowering of the total Rosseland opacity (Eq. <ref>), 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. <ref>) by a factor ranging from 0.005 to 200, see Fig. <ref>. We selected a representative column at r=0.2au from our disc model after 0.1 Myrs 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^-6 M_⊙/yr and F_ vis=4.14×10^6 erg/cm^2, which corresponds to the centre of the plot, where we reach T≈1600K in the midplane.




As we increase Ṁ_ acc from very low values, we reach the silicate sublimation temperature at about 7×10^-8 M_⊙/yr, beyond which the dust opacity κ_ R^ dust 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^-6 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 of 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^-6 to 10^-5 M_⊙/yr.   The exact value of that critical mass accretion rate will depend on the surface density, which is about 10^ 4 g/cm^2 in our model at 0.2 au after 0.1 Myrs.  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.



 §.§ 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 have determined using detailed radiative transfer calculations, and how the different minerals sublimate and re-sublimate 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. <ref>), 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 disk 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. <ref>) 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, no 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 midplane is warmer.  Here, when we compute the hydrostatic vertical structure in ProDiMo, we have implemented a cap of the local dust temperature by the silicate sublimation temperature (Eq. <ref>) 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 <cit.>. However, close to the stability temperatures of all minerals, sublimation and re-sublimation 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.




§ SUMMARY AND CONCLUSIONS


This paper proposes a pathway to form millimetre-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 Cal-Al-rich inclusions (CAIs) found in many chondritic meteorites.

Our model combines (1) the 1D viscous disc evolutionary models of <cit.> with (2) the 2D radiation thermo-chemical disc models of <cit.>, (3) the equilibrium condensation models of <cit.>, and (4) new dust opacity calculations.  


The models reveal a thermostat mechanism in the disc midplane regions close to the star <cit.>, 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 allow the Al-Ca-Ti particles to grow by coagulation and create a semi-ordered 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 millimetres; creating a size distribution that we can still observe today by measuring the sizes of the CAIs in chondritic materials, see e.g. <cit.>.

To form the CAIs, i.e. 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 <cit.> models, which is sketched as scenario 2 on the right side of Fig. <ref>.  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 re-condense on the existing surfaces of the Al-Ca-Ti oxide particles.  When this happens, 
temperatures are already too low for annealing, so this created an amorphous silicate mantel.  In contrast, according to the more common scenario 1 in Fig. <ref>, where all matter is accreted outside-in through the disc, at no time and space in the disc one would expect to find large, pure and crystalline Al-Ca-Ti rich particles embedded in an amorphous silicate matrix.

Our scenario for the CAI-formation in the early solar system hence 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 requirements are met, as shown in Fig. <ref> by the cyan overlap region, up to a disc age of about 50 000 yrs.  
Eventually, after the mass accretion rate has fallen to values 10^-7 M_⊙/yr and the viscous heating of the disc subsides, the thermostat mechanism breaks down, 
and the dust particles from the molecular cloud that continue to join the disc will do so mostly unaltered, i.e. without annealing and without losing their silicates of other refractory components on the way.



   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 (PLANETOIDS). 
   P. M. acknowledges funding from the Italian Ministerial Grant PRIN 2022, “Radiative opacities for astrophysical applications”, no. 2022NEXMP8, CUP C53D23001220006.






aa




§ OPTICAL CONSTANTS


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





§ PHASEDIAGRAMS FOR AL_2O_3 AND SIO_2





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 <cit.>: ϵ_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. <ref>.  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^-7 bar,1696 K), respectively.  For corundum, we find (7000 bar,2327 K) and (2.5×10^-5 bar,2327 K), respectively. 

In the H-rich case, H2 is the dominating molecule, providing the gas pressure, and most of the oxygen is locked in H2O.  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 O2, 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 <cit.>:  At the surface of a hot oxide grain, the collisions with H2-molecules and H2O 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 H2 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 H2.

