Shaping the CO snowline in protoplanetary disks

Context. Characterizing the dust thermal structure in protoplanetary disks is a fundamental task because the dust surface temperature can affect both the planetary formation and the chemical evolution. Because the temperature depends on many parameters, including the grain size, it can be challenging to properly model the grain temperature structure. Many chemistry disk models usually employ a sophisticated single dust structure designed to reproduce the effect of a realistic population presumably composed of a large diversity of sizes. This generally represents a good approximation in most cases. Nonetheless, it dilutes the effects of the complex radiative interactions between the different grain populations on the resulting dust temperature, and thus, the chemistry. Aims. We seek to show that the radiative interactions between dust grains of different sizes can induce a nontrivial dust temperature structure that cannot be reproduced by a single dust population and that can significantly affect the chemical outcome. Methods. The disk thermal structures were computed using the Monte Carlo radiative transfer code RADMC-3D. The thermal structures were postprocessed using the gas-grain code NAUTILUS to calculate the evolution of the chemical abundance. Results. We find that simultaneously using at least two independent dust grain populations in disk models produces a complex temperature structure due to the starlight that is intercepted by the upper layers of the disk. In particular, we find that micron-sized dust grains are warmer than larger grains and can even show a radial temperature bump in some conditions. This dust temperature spread between the grain populations results in the segregation of the CO snowline and in an unexpected CO gas hole in the midplane. We compare the results with observed close to edge-on class I/II disks. Conclusions. Our study shows that the size dependence of the dust temperature significantly impacts the chemistry, and that two dust populations at least are required to account for this property of the thermal structure in protoplanetary disk models over a wide range of disk masses and dust properties.


Introduction
Numerous studies have now confirmed the existence of grain growth in circumstellar disks, with grains of size up to millimeter detected (e.g., Testi et al. 2014;Ohashi et al. 2023).Some authors reported the presence of grain growth at a very early stage.For instance, Harsono et al. (2018) confirmed the presence of millimeter-sized grains in the very young disk (∼ 100,000 years) around the protostar TMC1A.These observations are of prior importance as they suggest that planet formation can initiate at a very early stage of protostellar evolution.
More specifically, it is known that the grain size distribution in protoplanetary disks is different from the one in the interstellar medium (ISM).Many authors have shown that the distribution can deviate from the typical Mathis, Rumpl, & Nordsieck (1977) model (hereafter MRN) found in the ISM (Ricci et al. 2010).In particular, the grain size range in disks is less narrow than in the ISM and is likely to be continuously distributed due to dust coagulation, fragmentation and radial drift (e.g., Dullemond & Dominik 2005;Brauer et al. 2008;Birnstiel et al. 2010Birnstiel et al. , 2018)), resulting in a typical range from nanometer sizes to centimeter sizes.The spatial distribution is also more complex.Grains with a small Stokes number are coupled to the disk gas (small Stokes numbers typically apply to small grains, and vice versa) and can be observed in near-infrared at very large vertical extent whereas grains with a larger Stokes number are more efficiently settled toward the disk midplane and occupy a thinner vertical extent (Fromang & Nelson 2009).Moreover, inward radial drift can occur for coagulated grains (e.g., Birnstiel et al. 2010;Lambrechts & Johansen 2014), strongly reshuffling the local dust density and size distributions.
Consequently, the extinction opacity in disks is also much more difficult to characterize than in the ISM.Most authors assume size-averaged opacities from a given local distribution (e.g., Pinte et al. 2016;Birnstiel et al. 2018).The opacities are commonly averaged over a typical dust size distribution from the MRN model, with the maximum grain size chosen to align with multi-wavelength observations.More precisely, most recent studies consider a dust model composed of two-grain populations, a small dust population and a large dust population (e.g., Du & Bergin 2014;Ballering et al. 2021;Schwarz et al. 2021;Zhang et al. 2021;Murillo et al. 2022) meant to fit observationally constrained evidence of grain growth and size-dependence to vertical and radial segregation (Gräfe et al. 2013;Pérez et al. 2015;Tazzari et al. 2016;Villenave et al. 2020).Typically, each population has its opacity averaged over a given size range, each with a different maximum grain size a max , the latter being a major parameter impacting the wavelength-dependent optical prop-erties (Birnstiel et al. 2018).A single resulting temperature is then extracted (density and temperature), either by averaging-out the two dust structures or by forcing dust thermal coupling during the radiative transfer simulations.The structure is later used for chemistry post-processes.However, since surface chemistry is strongly sensitive to dust surface temperatures, this treatment may overlook possible chemical effects by not considering the two dust structures independently.
Indeed, we saw that real disks are most likely composed of grains of multiple sizes with independent optical properties, resulting in polydispersed dust grains (Birnstiel et al. 2018).Note that the term 'polydispersion' (as opposed to monodispersion) in the present study is used in the strict sense as defined in Deirmendjian (1969) and Bohren & Huffman (1998), that is we speak of polydispersion when a population of scattering particles in a medium is uniform in shape and bulk material but not uniform in size.Grains are heated both by the external field and by the resulting scattered and re-emitted light from all the other grains.The absorption efficiency of grains is dependent on grain size and wavelengths which results in a) a complex interaction between the local radiation field and b) a size-dependent dust temperature structure (Heese et al. 2017;Gavino et al. 2021).It is therefore expected that utmost attention should be given to the selection of the dust population, dust absorption and scattering profiles in disk models.
In this work, we reproduce similar dust models as in recent studies (e.g., Zhang et al. 2021).However, we keep the dust structures independent so we can analyze (i) the precise dust temperature structures resulting from polydispersed dust grains computed with dust continuum radiative transfer simulations and (ii) the resulting chemical structure by treating the various dust populations simultaneously and independently in chemistry post-processes.We compare the results with test models composed of a single dust population meant to depict the classical treatment mentioned above.In particular, we will investigate the effect of polydispersion in the regions around the disk midplane and observe the impact on the CO snowline.
The paper is structured as follows: Section 2 describes the disk models used in the present study.Section 3 presents the resulting thermal structures and chemical post-process and, finally, in Sect. 4 we discuss in more details the origin of the temperature substructures and its potential implications on the CO gas distribution in observed disks.A summary is then provided in Sect. 5.

Model
In the present study, the models are assumed to be static protoplanetary disks, consisting of a typical 2D parametric smooth, axisymmetric and geometrically flared structure.Additionally, only passive heating is assumed.The disk is assumed to be in Keplerian rotation.We provide hereafter a concise description of the parametrization but the model is similar to Gavino et al. (2021).

Spatial distribution of gas
The radial gas surface density follows a simple truncated powerlaw, where p = 3/2 (Shakura & Sunyaev 1973;Hersant et al. 2009;Guilloteau et al. 2011;Le Gal et al. 2019) and R 0 is the reference radius chosen to be 100 au in the study.
The surface density at the given reference radius R 0 is derived using where M disk is the total disk mass and its inner radius and outer radius are R in and R out , respectively.Following Dartois et al. (2003), the gas temperature is imposed by analytical laws.The gas kinetic temperature in the midplane T mid is given by a power law where T mid,0 is the gas temperature in the midplane at R 0 , and we allow for a warmer disk atmosphere using the formulation of Williams & Best (2014) where σ = 2 is a stiffness parameter of the vertical profile.Above z atm , the temperature is vertically constant.T atm (r) follows a power law similar to T mid (r) and z atm is taken as being 4 times the mid-plane gas scale height, H g , which at the the reference radius R 0 is given by, where µ is the mean molecular weight equals to 2.4, m H is the proton mass, k B the Boltzmann constant, and G the gravitational constant.The mid-plane gas scale height then follows a radial power-law varying as r 3/2−q/2 , a typical flared geometry.Our assumed gas temperature and density profiles were selected for comparison with previous works and differ from those of a steady-state viscous disk with radially constant viscosity, which would have a r −1.1 dependence.Note, however, that a selfconsistent model would have a more complex radial density profile even under the simplified viscosity assumption because of the regions of increased temperatures that result from the treatment of different dust populations.
We solve the gas density by considering the (non-isothermal) vertical hydrostatic pressure equilibrium, that we solve iteratively following Hersant et al. (2009), where ρ g(z i ) is the gas mass density at altitude z i , Ω the Keplerian rotation at the given radius.The vertical gas density structure thus slightly deviates from a Gaussian profile at altitudes z > 2 − 3 H g .
Although this imposed density and gas temperature structure are not exactly self-consistent, we verified that the gas temperature T g (r, z) is not too different from the area-weighted dust temperature (calculated from Monte-Carlo simulations as described in Sect.2.2), as would be expected when the gas is mostly heated by thermal contact with the dust grains, rather than by the radiation field.
Here, a j represents the grain size of bin j and n d is the dust number density (details are given in the next section).Note that the chemical models can be computed using the gas temperature as given by Eq. 4 or derived from the area-weighted dust temperature (Eq.7).We systematically verified that both cases lead to very similar abundances for the most common molecules, as expected since the rates of the main gas phase chemical reactions are not strong functions of the temperature1 .We thus only present chemistry results from the model using the gas temperature as given by Eq. 4, the other choice resulting in practically identical figures.

Spatial distribution of dust
The vertical density distribution of any dust component is defined by the following equation: where σ d(r,a) and H d(r,a) are the surface density and scale height, respectively, of the grain population of size a at radius r.The dust scale height is defined following (see Dubrulle et al. 1995;Youdin & Lithwick 2007;Fromang & Nelson 2009;Dong et al. 2015), so that the dust scale height represents a fraction of the gas scale height, depending on the dimensionless stopping time, or Stokes number, St. S c is the Schmidt number and α the turbulent viscosity coefficient, set to 1 and 0.01, respectively, in the whole study.The coefficient α sets the degree of turbulent viscosity and therefore the degree of settling (Shakura & Sunyaev 1973).The settling factor depends on the dust grain size a via the Stokes number, which can be defined near the midplane as where ρ m is the dust material density and under the assumption of Epstein regime and spherical particles (see Cuzzi et al. 2001;Birnstiel et al. 2012).
The dust-to-gas mass ratio varies from an altitude to another due to the settling, but the dust-to-gas surface density ratio is ζ at all radii.The model does not include radial drift.
We point out here that dust settling is not expected to result in a simple Gaussian distribution.For instance, Fromang & Nelson (2009)'s simulations and analytical approach (their Eq. 19) show that in a vertically isothermal disk, the dust settling results in a dust mass distribution that is approximately Gaussian near the mid-plane (z < 2H), but falls more steeply than a Gaussian at high heights.Nevertheless, Eq. 19 of Fromang & Nelson (2009) is only valid for a Gaussian gas density vertical distribution.In our case, the increased temperature at high heights above z/H > 2, where the dust profile starts to be steeper, leads to a profile that has broader wings than a Gaussian (our Eq. 6).The two effects partially cancel and the settled dust distribution is thus quite likely better approximated by a simple Gaussian than in the vertically isothermal case.In order to make sure that the extra term in Eq. 19 of Fromang & Nelson (2009) does not cause major differences, we performed additional simulations using their approach.These tests have showed very minor changes, we therefore only present the results using our Eq.8.

Dust mass and size distribution
In all the paper, the size distributions follow the MRN distribution with a power-law dn(a) ∝ a −d and with d = 3.5.It is possible to divide the distribution into N d intervals (logarithmically distributed in all the study).From this, a mass and a size can be defined for each grain bin j.Following Gavino et al. (2021), the effective discretized size value of the jth interval is given by the mass-weighted average, where a +, j and a −, j are the maximum and minimum cutoff values of the jth interval, respectively.The total dust mass in the disk of the population j of size a j is a fraction x(a j ) of the total disk dust mass M dust : with ζ the dust-to-gas mass ratio and M disk the total disk mass (gas+dust).The term x(a j ) derives from the MRN distribution and is a fraction of the full size distribution where a min and a max are the minimum and maximum cutoff values, respectively, of the full size distribution and such that The fiducial physical parameters used in all disk models are summarized in Table 1.We note that these parameters are chosen to provide an interpretative framework for the resulting thermal structure and are not intended to directly compare with observations.However, the models are chosen to be similar to a typical low-mass T Tauri disk.1.675 g.cm −3 α (turbulence coefficient) 10 −2 T mid,0 (chosen gas temperature at R 0 ) 25 K q (radial temperature profile exponent) 0.4 p (surface density exponent) 1.5

Radiation source and radiative transfer modeling
For the stellar radiation we consider a pre-main sequence star emitting as a blackbody with an effective temperature of T star = 4100 K, a stellar luminosity of L star = 1 L ⊙ , and a stellar mass of M star = 1 M ⊙ .As for the interstellar radiation field (ISRF), we consider a typical distribution from Draine (1978) with an extension of van Dishoeck & Black (1982) for wavelengths > 200 nm.We use these parameters for all models described in this study.The stellar irradiation will be the dominant source of thermal heating for the dust and gas.On the other hand, the photochemistry of the disk surface will be mainly controlled by the FUV spectra of both the ISRF and stellar radiation field.
To compute the dust temperature we perform radiative transfer calculations using the RADMC-3D code (Dullemond et al. 2012).Each dust population has its own density distribution and are allowed to be thermally decoupled from each other.The radiative transfer calculations include the absorption opacity, scattering opacity, and the Henyey-Greenstein g parameter for the treatment of anisotropic scattering.
The structure is a spherical grid centered at the star's location, with 300 logarithmically distributed radii, and 180 angles from 0 to π.We use 100 logarithmically distributed wavelengths from 0.1 µm to 2 mm for the radiative transfer calculation.We use 10 7 photon packages in each simulation.

Chemical network and gas-grain simulation
We use the three-phase Nautilus Multi-Grain Code (NMGC) (Ruaud et al. 2016;Iqbal & Wakelam 2018;Gavino et al. 2021), derived from the NAUTILUS gas-grain code (Hersant et al. 2009).NMGC is specifically designed for chemistry simulations of models embedding multiple, independent grain populations.All dust populations are chemically active and all chemical compounds can interact with each population independently.The chemistry model is a 1D+1 structure spanning from 4 au to 360 au from the star.All 1D structures are composed of 64 vertical spatial points from the disk midplane up to 4 H g , defined as the surface of the disk.
NMGC uses the rate equation approach (Hasegawa et al. 1992;Hasegawa & Herbst 1993) in which the gas phase, the grain surfaces, and the mantles are chemically active (Ruaud et al. 2016).The gas phase can exchange species with the surfaces but not directly with the mantle.For surface chemistry, adsorption, thermal desorption, cosmic-ray induced desorption, Following the prescription of Gavino et al. (2021), the H 2 formation rates are calculated by consistently interpolating rate curves which are initially solved by using a stochastic treatment that considers fluctuations of hydrogen atoms on the grain surfaces by Bron et al. (2014).
We use the KInetic Database for Astrochemistry network (kida.uva.2014) 2 (Wakelam et al. 2015), to which we include the additional reaction rates for the H 2 formation on grain surfaces.The chemical evolution is simulated over a period of t chem = 2 Myrs, well above the CO freeze-out timescale in the whole disk.Considering the relatively large simulation run-time, we use atomic initial conditions (except for H 2 ) as presented in Table 2. See Wakelam et al. (2019) for more information on the impact of initial conditions on the chemical evolution.
For the species other than H 2 and CO, the photoprocess rates are calculated using the following equation: Where G 0 is the unattenuated local UV flux at each radius relative to the standard ISM (Draine) value.Note that the γ value are taken as the default values appropriate for the ISM and that for many species γ = 1 in the Kida network.The A v is computed using the local UV flux calculated from the Monte-Carlo radiative transfer simulation.
For a given wavelength λ and given coordinates (r, z), the local flux is the unattenuated flux F 0 reduced by the dust extinction: where the opacity is related to the extinction by the usual relation

Dust models and optical properties
To illustrate the importance of polydispersion, we consider four dust models.Model S (for Single) is composed of a simple single dust distribution.Model C and D are models using two dust distributions.In Model D (for Dual), each dust population has its own temperature, while for Model C (for Common), an averaged temperature is derived from the two temperatures of Model D following Eq. 7. The difference between Model C and D therefore lies entirely in the chemical computation, as we will see.
Lastly, Model M16 (for Multi) is composed of 16 dust grain size bins.The four models are described below, while the thermal and chemical results are shown in Sects.3.2 and 3.3, respectively.A summary of the four models is given in Table 3.More specifically, the number of grain size bins used during radiative transfer and chemistry simulations as well as the effective grain size values are indicated in columns 2, 3, and 4, respectively.

Model S, D, and C
We adopt the opacities from the DSHARP collaboration by using the standard size-averaged smoothed optical properties of the dsharp_opac package3 as described in Birnstiel et al. (2018) with the mixed mie coefficients used for the DSHARP project.
For Model D, we reproduce the same distribution as in Zhang et al. ( 2021) (who uses the DSHARP opacities), for each of the two populations we use a size-averaged opacity derived from a dust grain size distribution following the MRN distribution with a power-law index q = 3.5.The dust grains are composed of a mixture of water ice (Warren & brandt 2008), astronomical silicates (Draine 2003), troilite and refractory organics (Henning & Stognienko 1996).We consider two populations, a small population D small with an upper grain radius a max,s = 1 µm, and a large population D large with an upper grain radius a max,l = 1 mm.Both populations have a lower grain radius of a min = 0.005 µm (e.g., Schwarz et al. 2018;Ballering et al. 2021;Zhang et al. 2021).The wavelength range is a logarithmic distribution of 210 values between λ min = 0.1 µm and λ max = 10 4 µm.Because the small population is globally less protected from photodesorption, we set its ice mass fraction to zero (following the prescription in Zhang et al. (2021)).The resulting mass-weighted dust absorption opacities κ abs , dust scattering κ sca and the g parameter of anisotropy are shown in Fig. 1.These profiles are to reproduce as closely as possible the ones in Zhang et al. (2021).The thermal calculations as well as the chemistry post-process are performed using the two populations simultaneously.
For Model S, we only use the small population opacity D small (dashed lines in Fig. 1) for the computation of the dust temperature.A single dust structure is thus considered at all steps, including the chemistry post-process.
Model C is the same as Model D as it uses the two dust size distributions and optical properties simultaneously for the dust temperature computation, but differs during the chemistry postprocess step.Indeed, the two resulting dust temperatures are averaged over the dust surface area using Eq.7 so a single dust structure is used as input to solve the chemical evolution.Model C is therefore close to the prescription that has been commonly used in most disk models so far.

Model M16
It is expected, based on models of dust coagulation and fragmentation, that a dust population is not composed of only two distinct populations but rather of a continuous size distribution (Birnstiel et al. 2010(Birnstiel et al. , 2018)).To address this, we consider Model M16, which mimics a continuous grain size distribution by splitting the MRN distribution into 16 logarithmically distributed grain size intervals, each with its own size (Eq.11) and associated single-grain-size opacities.Model M16 is therefore meant to be more realistic than the other models.
Similarly to the other models, we adopt a classical powerlaw grain size distribution ranging from a min = 5 nm to a max = 1 mm with a size distribution index q = 3.5.However, instead of having a total absorption opacity averaged over the total size distribution (as it is done in the other models), the resulting opacities are discretized into the 16 bins (see Table 4).
To do so, we perform calculations for each grain bins using the Fortran subroutine included in the dsharp_opac package (Birnstiel et al. 2018), itself derived from the original Mie code from Bohren & Huffman (1998).To avoid unrealistic interference effects due to the spherical treatment of the Mie theory, we calculate the opacity for 30 linearly distributed grain sizes around each of the 16 bin sizes.
To simplify the discussion in the following, we differentiate three sub-populations from the 16 sizes, each of which characterized by similar optical properties.We call small grains all grains of size < 0.1 µm (4 bins), intermediate grains all grains of sizes between 0.1 and 10 µm (6 bins) and large grains all grains of size > 10 µm (6 bins).
The resulting mass-weighted absorption opacities κ abs (a) = 2πa 2 Q abs (a)/m(a) [g/cm 2 ] and scattering albedos ω are shown in Fig. 2. The absorption opacities (Fig. 2a) of the large grains (solid lines) are roughly flat in short wavelengths (< 100 µm) whereas they all asymptotically follow the same decreasing slope beyond.The large grains also prefer to emit at wavelengths comparable to their size.Intermediate (dashed lines) and small grain (dotted lines) opacities are strongly wavelength-dependent and show structured patterns due to their material composition.These various profiles result in various equilibrium grain temperatures and thermal emissions.
The small grains and intermediate grains show little albedo or are completely opaque (Fig. 2b) at wavelengths λ > 10 µm.The albedo ω is at its maximum when 2πa ∼ λ (Birnstiel et al. 2018) and can be as high as 0.9 between ∼ 0.1 to 10 µm in the case of the intermediate grains.The albedo of the largest grains (bin 12 to bin 16) remains globally flat for λ < 100 µm at a fairly large albedo (0.5 < ω < 0.6).However, since their relative surface area is very small (< 1 % in the midplane), their contribution to scattering events leading to significant effects on the thermal structure is assumed to be negligible.

Comparing models
The total disk mass is kept constant in all models.Thus, when comparing models, their respective optical depths can differ since the dust mass is not distributed among the same dust populations.Typically, Model D and C are expected to exhibit a slightly smaller vertical optical depth than Model S due to the presence of the large dust population in the former.However, these opacity differences are rather small and the differences between models can be mainly attributed to their different dust temperature structures.Note that we tested various dust models and 10 1 10 0 10 1 10 2 10 3 [ m]   4 for description of bin sizes.

Results
In this section, we evaluate the effects of different assumptions on the dust models on the thermal and chemical outcomes.First, we describe the 2D dust structures.

Dust structures
Figure .3 shows the 2D dust area densities of Model D. The choice of showing the area density is motivated by the fact that the dust surface area is one of the main parameters that drives the gas-grain interactions during chemistry.As detailed in Sec. 2, the dust populations are computed independently to one another so the two dust populations are to be shown separately.
Using Eq. 11 and Eq. 12, we can define two discretized subranges representing two populations with their own size and mass fraction values.The small dust population (Fig. 3a) is composed of grains of size a small = 0.02 µm and, despite representing less than 5 % of the total dust mass, is the larger contributor to the extinction as the grains represent most of the dust surface area.The large dust population is composed of grains of size a large = 10.4 µm and represents 95.5 % of the total dust mass (Fig. 3b).
The settling factor is defined using Eq. 9.The large population is more confined toward the midplane than the small population due to its significantly larger settling factor.Indeed, the Stokes number ratio of the two populations is S t,large /S t,small = a large /a small = 5.2×10 2 .While the settling factor of the small population is close to unity at all radii (the small dust scale height follows that of the gas), the ratio between the large dust and small dust scale height equals to 0.92 at 10 au, 0.35 at 100 au and 0.20 at 500 au.The two populations exhibit notable differences.In particular, the 25 K isothermal is located at ∼ 30 au from the star (in the midplane) for the large grains (Fig. 4, right), which is a typical expected distance where the CO snowline is located in a T-Tauri disk.On the other hand, the small grain temperature drops below 25 K near 180 au only (Fig. 4, left).There is therefore a significant temperature spread with radial distance of 150 au between the two 25 K isothermal.

Comparison between Model D, S, and C
To have a better understanding of the temperature structures we further investigate the disk midplane by comparing the radial dust temperatures of Model D with that of Model S and C. Figure 5 shows the radial profiles in the midplane of the three models.In Model D, the two populations are thermally decoupled from r ∼ 30 au and the small grains are warmer than the large grains of approximately 13 K at all radii.Moreover, the small dust population temperature shows a surprising bump between ∼ 30 au and ∼ 80 au.This structure therefore deviates from an expected monotonously decreasing temperature along the radial distance.
In Model S, the dust optical properties are the same as the small population of Model D. We can therefore expect that their respective dust temperature exhibit similar profiles but Fig. 5 shows that it is not the case.The temperature profiles of the two models overlap in the optically thick (τ >> 1) inner region (r < 30 au) due to the very effective temperature coupling.Further away, the two populations in Model D become spread out as they become thermally decoupled.Such temperature spread does not exist in Model S since the dust structure is unique.This also implies that there can be no bump in Model S.
The optical depth at both visible and infrared (IR) domains between the star and any radial distance in the midplane is too large (τ > 10 4 at 1 µm), thus the re-heating of the small dust population cannot be generated by direct starlight emission.Since this bump does not exist in Model S, we can assume that it originates from the interaction of the re-emission and scattered light between the two dust populations.This assumption is also supported by the vertical optical depth.The vertical solid, dashed, and dotted green lines represent the radial distances from which the vertical optical depth at 1 µm becomes smaller than 10, 5, and 1, respectively (in Model D).We see that the dust temperature starts to increase as the vertical optical depth goes down to values close to 10, where dust radiative interactions between upper layers and lower layers start to be possible.
As for Model C, the temperature profile follows closely the profile of the small dust population of Model D because most of the dust surface-area is from the small grain populations (> 95 %).

Model M16
Let us go further and test the effect of polydispersion in a more realistic structure by considering Model M16. Figure .6a shows the resulting radial temperature profiles of the 16 grain species in the midplane.There is a temperature spread just like in Model D, with the intermediate grains being warmer than the small grains and the large grains, and with a difference of about 17 K between the warmest grain and the coldest one.The model also shows a temperature bump inside 100 au.It is clear that the grain populations can be divided into two groups, one group that undergoes a bump effect (small and intermediate grains) and one that does not (large grains).The division in two groups is in fact analogous to the one in Model D. This suggests that a two-population Model Such as Model D could be sufficient to properly approximate the effect of polydispersion taking place in a realistic dust grain population.In order to help us understand which of the absorption, emission, or scattering processes as well as which of the grain populations contribute to the bump effect, we successively set the scattering opacity to zero inside wavelength intervals for each grain population and run radiative transfer simulations.We find that the radiation of wavelengths between 0.7 and 2 µm scat-  We can conclude that when scattering is included (Fig. 6a), the temperature degeneracy of small and intermediate grains is raised, the spread is larger, and a bump is generated.Note that a similar test with Model D shows that setting the albedo to zero also suppresses the temperature bump.This effect will be discussed in more details in Sect.4.1.2.
Another look at the temperatures shows that this temperature spread is located right in the vicinity of the CO freeze-out temperature (typically around 20-30 K in the disk midplane).We test the effect of the temperature spread on the chemical evolution in the following sections.

Limitations
It is easy to anticipate the effects of the temperature bump on the chemical abundance of light species, particularly that of CO.However, computation of the chemical evolution using as many grain sizes as in Model M16 can be highly cpu-consuming (a single 1D structure needs ∼ 60 h in cpu-time to be solved on a typical laptop compared to around ∼ 1 h for a model using a single grain) and can generate convergence issues due to the sizes and very small abundances of the largest grains.To duck this obstacle, we choose to use Model D instead of M16 to compute chemistry, based on the assumption previously made on the equivalent thermal structure in models D and M16.Said other- wise, we assume that the thermal spread from these models will have similar chemical outcomes.In addition, to fully appreciate the effect of the temperature spread, we also compute chemistry in Model S and C.

CO in the gas phase
Figure 7 shows the integrated surface density of CO in all 3 models.While Model S and C result in similar profiles, Model D shows a clear depression in CO near 150 au, as a result of the combined effects of the small and large grains.
In order to understand where those differences in column densities come from, we show in Fig. 8 the resulting 2D CO gas number density [cm −3 ] structures in (r,z) coordinates, from r = 2 au to 360 au using Model S (Fig. 8a), Model C (Fig. 8b), and Model D (Fig. 8c).The solid white lines marks the transition where the ratio n g (CO)/(n g (CO)+n s (CO)) becomes smaller than 0.5 or said otherwise the CO snowline.
The interpretation of the CO gas structure in Model S is straightforward.The disk surface is vulnerable to the external UV radiation and the photodissociation rate of CO is large enough to prevent efficient formation.In the molecular layers the H 2 shielding and CO self-shielding become sufficiently protective and CO gas becomes the most abundant molecule (after H 2 ) with an abundance close to the interstellar value 10 −4 .The grain surfaces are still too warm for effective CO condensation.Closer to the midplane, the grain surface temperatures drop and CO becomes massively depleted onto the grain surfaces.The CO snowline delimits the molecular layer and the cold midplane.Model C shows a similar pattern with a larger molecular layers.This is not surprising since the dust temperature is globally larger than in Model S. This 'V-shape' structure is typically observed in protoplanetary disks (see van't Hoff et al. 2020, for a detailed description of the CO gas phase distribution in disks).
The CO gas map in Model D shows a more intriguing structure.Indeed, although the structure in the upper layers is similar to Model S and C, we note the presence of a CO 'hole' from ∼ 50 au to ∼ 200 au and extending up to ∼ 10 au vertically.We then observe a CO gas rise between 200 au and 300 au before diminishing beyond 300 au.The CO snowline exhibits a total different   shape as well, with a gas-ice transition displaced toward the midplane between 150 au and 250 au.
By looking at Fig. 5, it appears that the CO gas rise is located well outside the temperature bump (around ∼ 70 au) and the CO depression is located right where the dust temperature is at its maximum.This striking discrepancy between CO depression and dust temperature is due to the dust distribution producing multiple condensation fronts, as shown in the following section.

CO ice and snowline
Figure 9 shows the 2D number density maps of the total CO ice in Model S, C, and D. All models present the typical dual structure split by the snowline, below which CO is massively depleted onto the grains.With a global higher dust temperature, the snowline of Model C is shifted toward lower layers of the disks as compared to Model S.
In Model D, the total CO ice content is the sum of CO ice adsorbed onto the small and large grain populations.It is clear that CO ice lies under the snowline profile, and although the snowline is clearly different than in the other models, the CO ice distribution seems to follow a similar trend as the other models.Additionally, we also present the snowlines of each population taken separately (red dashed lines).The snowlines appear to be divided into two individual surfaces, overlapping near 250 au.To better understand, let us provide a formal definition of the snowline in the case of multiple discretized dust grain sizes.The freeze-out temperature is tantamount to solving the balance between the desorption and adsorption rate, and depends on multiple parameters among which the dust temperature and dust and gas densities (see for instance van 't Hoff et al. 2017).Each grain population having its own density and temperature, we can define a snowline for each of them.More precisely, if n small (CO) and n large (CO) represent the CO ice number density on the small and large population, respectively, and n g (CO) the number density of CO gas, the two dashed white lines thus correspond to the location where the ratios n g (CO)/(n g (CO) + n small (CO)) and n g (CO)/(n g (CO) + n large (CO)) equal 0.5.On the other hand, the total snowline (white line) is the location where In the case of Model D, the number of grain bins N = 2 and Eq.19 writes n ice tot (CO) = n small (CO) + n large (CO).While Model S and Model C exhibit an expected structure of the snowline, the shape of the snowline in Model D is more complex.This is discussed in more details in Sect.4.2.

Thermal decoupling in the disk midplane
Figure 10 shows the 2D map (grey scale) of the optical depth at 1 µm in Model D computed from the star to all points in the disk τ ⋆ (r, z).We see that the optical depth values located below one scale height z < H g (white dashed curve) are very large (τ ⋆ » 10 2 ) such that the disk midplane is always optically thick in regards to the star.The two dust populations are therefore expected to be strongly thermally coupled: the warmer population is cooled by the other one, and vice versa.However, outside 30 au, the grains become thermally decoupled (the small grains become warmer), even though the optical depth remains very large.
To take a closer look, the blue line marks the altitude where the vertical optical depth calculated from the midplane equals 1.This way we can investigate from which altitude radiation contributes to the heating of the midplane.We see that at small radii (r < 30 au), the radiation reaching the midplane is exclusively due to thermal emission from nearby dust with very similar temperature, hence the temperature coupling.At larger radii (r > 30 au), we see that the altitude of the blue line suddenly increases to reach the uppermost layers of the disk at around r = 100 au.On the other hand, the red lines corresponds to the surface where the optical depth from the star τ ⋆ equals to 1.The layers above that surface are therefore directly heated by the star.At the coordinates (r = 70 au, z = 11 au), the blue line crosses the red line, meaning that the optical depth from the star to this point and from there to the midplane equals 2, a value much smaller than the optical depth calculated from the star through the midplane.The midplane can thus be effectively heated by stellar light coming from the disk surface.
However, since each grain has a different absorption/emission opacity, they are not heated with the same efficiency.Figure 11 shows the ratios of the absorption cross section for wavelengths at which the star typically radiates to the absorption cross section for wavelengths at which the dust typically emits (∼ 20 µm) as a function of the grain size.In Model D, the small population has a ratio close to 4 whereas the large population has a ratio close to 1.2.The starlight at 1 µm is thus more effective in heating the small grains than the large grains.Likewise, in Model M16, the intermediate grains have the largest ratio values, so the heating is maximum for these grains.

The 'Bump effect'
On top of that, we described a surprising effect where the temperature re-increases, reaches a maximum and then monotonously decreases again.This is what we called the bump effect.
In Sec.3.2 we demonstrated that scattering is at the origin of the bump effect in the dust temperatures.
To illustrate the importance of the contribution of scattering of stellar radiation, we indicated the wavelength λ max for the temperature at which the star emits according to Wien's law in Fig. 2a.
In Model M16, we saw that the bump effect was caused by the scattering of stellar radiation by the intermediate grains where a fraction of the starlight intercepted in the upper layers is scattered vertically toward the midplane.In the case of Model D, this is the small grains that scatter the starlight vertically (we note that the large grains are strongly settled so they cannot effectively intercept the starlight in the upper layers.They can, however, interact with the starlight already scattered by the small dust population in the intermediate layers).
More generally, the crucial combination for a bump to occur is the simultaneous presence of at least one scattering dust grain population with a large scattering albedo accompanied with the presence of at least one other grain population with a large absorption to emission opacity ratio.
In Model S, no bump effect is visible, even in the optically thin regions.This comes obvious given the mechanisms needed to generate a bump, as we saw that at least two dust species with different opacities are needed, so that the various dust populations can have different temperatures in regions where temperature coupling is not relevant.Therefore, a bump effect cannot occur in a dust model using a single opacity profile, regardless of how sophisticated the model is.

Sensitivity to other parameters
Thermal decoupling and the bump effect originate from different mechanisms.While thermal decoupling always occurs in optically thin environments and does not require radiative interaction between different grains, a bump effect is generated on certain grains by the scattered light coming from other grains in their vicinity.It is therefore sensitive to the chosen optical properties and in particular to parameters that affect the albedo: the material composition, the size range, etc.To illustrate this, we performed again dust continuum radiative transfer simulations for Model D but this time by using the DIANA standard opacities in Appendix A.1, and by using another size distribution in Appendix A.2.We also tested the sensitivity to different disk masses (Appendix A.3) and different turbulence parameters (Appendix A.4).As expected, while the magnitude of the bump is sensitive to the opacity properties, a temperature dichotomy exists in all cases, though at different radii.
Finally, the reader is referred to Appendix B where we address the question of the possible numerical origin of the bump effect arising in optically thick regions.Fig. 11: ratios of the absorption cross section for wavelengths at which the star radiates (∼ 1 µm) and the absorption cross section for wavelengths at which the dust emits (∼ 20 µm) for the grains in Model D (red points) and M16 (black points).

Snowline shape in Model D
To understand the snowline shape requires to first consider the separate contributions of the two grain populations of Model D.
Figure 12 shows the CO ice on each population in Model D. Figure 12a represents the distribution of CO adsorbed onto the small grains only whereas Fig. 12b  (Eq.18) whereas the red dashed lines are the snowlines of the individual populations.
In the midplane, the CO snowline for the small dust population is located close to 250 au where the dust temperature has sufficiently decreased.Inside 250 au, CO ice on the small grains is virtually non-existent because the temperature is above the freeze-out threshold for this population.For the large grains, the CO snowline in the midplane is located at a much inner radius (∼ 25 au) where the dust temperature drops below the freeze-out temperature for this population.We can straightforwardly conclude then that only the large dust population is at the origin of the CO gas hole visible in Fig. 8c.
On the other hand, we see that the large grain population also exhibits another midplane snowline (near ∼ 250 au, also where the midplane snowline of the small grains is located), so that the overall snowline of the large grains forms a closed Beyond that bubble, we can notice that CO is virtually non-existent on the large grains (Fig. 12b), even though their surface temperature is low enough for CO to adsorb.This behavior can be explained by multiple factors, but the most crucial one is the competition of interactions with CO molecules between the large and small grains.Indeed, the dust surface-areas (and therefore the number of accessible sites) is never the same between the two populations.For instance, at 250 au, the ratio between the small grain surface-area to the large grain surfacearea is n small (CO)a 2 small /n large (CO)a 2 large ∼ 6.5, meaning that a CO molecule has a higher probability (about 6.5 times higher) to interact with a small grain than a large grain.Consequently, Once the small grains are below the freeze-out temperature, they retain most of CO at the expense of the large grains.Hence, the sharp delimitation between the two snowlines.
The rise in CO gas observed at the delimitation (Fig. 8c) is also a consequence of this competition.Between r = 150 au to r = 250 au, the desorption to adsorption rate ratio progressively increases, and reaches the value 0.5 close to 250 au (for the small grains).At that location, CO interacts more with the small grains, but the latter still being slightly too warm, is returned to the gas phase ultimately.In other words, the freeze-out timescale around 250 au is longer than elsewhere, and the snowline assumes this peculiar S-shape with a minimum height at 250 au.
Note the presence of a small isolated 'island' of CO ice on small grains (Fig. 12a) at ∼ 30 au around the midplane, which corresponds to the location of the temperature dip, where the large grains have cooled enough the smaller grains.This happens because of the specific dust mass and structure of the disk.Should the disk be slightly less massive, this small island would likely not exist.
Lastly, with two grain populations and therefore two temperatures, Model D is an efficient simplification of the dust distribution to understand the impact of a proper dust temperature onto the chemistry and in particular the behavior of snowlines.Nevertheless, a more realistic distribution with more grain populations and more related temperature would necessarily produce more "individual" (one per grain population and temperature) snowlines, making the resulting CO snowline even more complex.Although the resulting CO ice segregation should be hard to probe in observed protoplanetary disks, it requires an interpretation on its own that may have significant implications for the chemical evolution and planet formation.

Radial drift and planet formation
Our model does not include any radial drift so the CO snowline shape is solely due to the different dust temperatures, though its shape and precise location is also affected by dust settling.
Cleeves (2016) has shown, using a simple parametric model, that radial drift can reshape the CO distribution and create multiple snowlines.In particular, the removal of a large fraction of dust grains from the outer disk can allow the region to become warmer by reprocessed stellar radiation from the upper layers.We showed that such reprocessing can occur even without a significant change in dust-to-gas mass ratio and only requires the presence of different dust populations.We can therefore expect that the combination of the two mechanisms (i.e., temperature spread and radial drift) should build up the radial segregation of the snowlines.
Radial drift occurs when gas drag acts on the dust grains (the gas disk 'feels' its own pressure and is therefore sub-Keplerian).As a consequence, the dust grains drift inward.The drift efficiency increases with the Stokes number St < 1 under Epstein drag law (Birnstiel et al. 2012).If, for instance, we consider the grain distribution in Model M16, the large grains will drift more efficiently than the intermediate grains, whereas the small grains remain mostly coupled to the gas and will barely drift.Although the drift time is expected to be relatively large for the whole grain size range considered (St < 10 −2 for all grains of size < 1 mm at 10 au in our model), radial size-dependent segregation can occur with enough orbital periods.We can make an estimate of the radial segregation between all grain populations by considering the following drift timescale equation (Birnstiel et al. 2012), where V k is the Keplerian velocity at radius r, c s is the isothermal sound speed, and γ = | − p − q 2 − 3 2 | is the absolute value of the power-law index of the gas pressure.The drift timescales at 100 au of all 16 grains are given in the fourth column of Table 4.These timescales can be compared with our chemical evolution time t chem = 2 Myr.We can see that all grains larger than about 10 µm have a drift timescale < t chem , whereas the small grain populations virtually do not drift (τ drift (100 AU) >> t chem ).We can therefore expect a strong radial size-dependent segregation within the chemical evolution time.The colder and larger grains should be found within a smaller radial extent than the smaller grains, such that the surface-area ratios become even more in favor of the smaller grains at large radii.
Accordingly, we expect that the effect found in Model D, which depends on a competition between small and large dust grain populations, will be moved inward by the radial drift of large grains, and one might expect an even more complex shape for the 2D snowline curve than in our Model D. Likewise, this effect should be enhanced by inward drift in a more realistic distribution like in Model M16.This could have observational implications (see Sect. 4.4).
Our study ignores grain-grain collisions, which may lead to more homogeneous ice compositions for grains of different size.However, grain-grain collisions are a complex process.While it is conceivable that slow velocity collisions between a small grain and a large one effectively results in the ice coverage of the small grain being incorporated in the ice coverage of the larger one, faster collisions, those which replenish the small size distribution from the large size one, are unlikely to have such a simple outcome.Grain growth has been shown by Van Clepper et al. (2022) to affect the chemistry, with an impact on the C/O gas phase ratio as a function of disk age, but their model did not consider disruptive collisions.Our chemical model also ignores diffusion, in particular vertical mixing, which can homogeneize the gas phase abundances Hersant et al. (2009).
In all cases, however, accounting for the size dependence of dust temperatures leads to higher temperatures than using a common temperature, so that this process should be considered when evaluating the chemical evolution.
On another hand, our mechanism also has implications for the chemical evolution during planet formation.Indeed, Model D shows that the small grains are virtually 'naked' from CO ice inside the large grain snowline whereas the large grains are naked from CO inside the small grain snowline.This suggests that the process of CO ice toward complex organic molecules (COMs) takes place only on the large grains in the inner regions where planet formation takes place, whereas it occurs only on the small grains in the outermost regions, which later act as a reservoir of solid particles after they slowly grew toward drifting larger grains (Youdin & Shu 2002;Lambrechts & Johansen 2014).
Additional chemistry models which take into account sizedependent radial drift are going to be investigated in future work.

CO holes in observed disks
We found that Model D shows a CO depression of up to 3 orders of magnitude in density and a factor of 10 in surface density, forming a 'hole' between 50 au and 200 au, followed by a rise beyond.Such behavior is reported in the case of the observations of edge-on disks.Dutrey et al. (2017) analyzed the temperature structure of the edge-on Class II disk 2MASS J16281370-2431391 (the 'Flying Saucer').They found a CO and CS temperature drops inside 200 au followed by a rise between 200 and 300 au.They argued that an effective rise of UV penetration can explain the heat up of the disk in the outer region (200 au being the outer radius of the 'mm' dust disk).Similarly, Flores et al. (2021); Villenave et al. (2022) reported in the edge-on Class II disk around SSTC2D J163131.2-242627(Oph 16313) a rise in the CO temperature in the outermost region.They also argued that an external UV field can heat the outer disk where UV photons can penetrate more effectively.
More recently, Lin et al. ( 2023) detected a similar behavior in both the 12 CO and 13 CO emissions in the close-to edge-on Class I disk IRAS 04302+2247 (the 'Butterfly ' Wolf et al. 2003) which closes off the freeze-out region, forming a hole.However, contrary to the two other disks, the rise does not coincide with the outer radius of millimeter-continuum emission.Furthermore, Gräfe et al. (2013) showed that in this disk the population of large grains has drifted inward whereas the modeled scattered light indicates that the small grain distribution has a larger radial extent invisible to mm-wavelength.
Hence, in these objects, the rise of CO appears beyond a more or less important drop of the dust-to-gas ratio (although it is less obvious in the case of the younger Butterfly star where the dust is less evolved).Such a behavior can be qualitatively explained by Cleeves (2016), where a CO rise occurs naturally because the dust-to-gas ratio drops quickly beyond the 'mm' dust disk edge allowing for a raise of temperature.
The Cleeves (2016) model requires an effective change in the (small grain) dust-to-gas ratio, due to radial drift.On the contrary, our model shows that even for a radially relatively constant value of this ratio, the temperature difference between small and large grains can result in a similar feature.This will be investigated in a future work.
Finally, Fig. 7 shows that the CO surface density, although of similar amplitude, clearly differs for Model D compared to the others.Model D, with two-grain temperatures, shows a dip between 70 and 230 au which reflects the CO depression observed in Fig. 12. Observing a disk with a moderate inclination, this feature can apparently mimic an unresolved gap which may be then erroneously attributed to planet formation.
In any case, we note that the observation of the CO snowline shape is only feasible for disks seen at high inclination, ideally close to edge-on.In less inclined disks, the substantial CO opacity from the top layer will hide the empty region behind the upper layers.The surface density effect shown in Fig. 7 should however be visible using isotopologues, principally 13 CO.
Other manifestations of the temperature gradient may be found with other molecular probes.The differences in temperature and timescales will also affect the surface chemistry because of the temperature dependence of the atoms and molecule mobility on the dust surfaces, and thus, in return, the chemical composition of the gas phase.Grids of models may be required to quantitatively evaluate the impact and possibly find out more sensitive indicators of this behavior than the rotational lines of CO.

Summary
We tested the impact of polydispersion in a medium composed of multiple independent grain populations on the resulting dust temperatures and chemical evolution in protoplanetary disk models.More specifically, we compared disk models using 1, 2, and a more realistic 16 individual grain populations with sizedependent vertical settling, densities, and optical properties and showed the effect on the CO distribution in the disks.The main results can be summarized as follows: -The presence of multiple dust populations with different optical properties generates a systematic dust temperature spread, where the smaller grains are always warmer than the larger grains, in regions of the disk where dust can be thermally decoupled.This is true even in the midplane of massive disks.Interestingly, the temperature of micron-sized grains can show a radial bump: it slightly re-increases with the radial distance due to stellar scattered light coming from the upper layers when the vertical optical depth decreases.-Our work shows that opacities derived from a size distribution are not well-suited to reproduce the temperature spread that exists in a dust population composed of multiple sizes because these opacities dilute the effects of polydispersion, regardless of how sophisticated the opacity models are.In particular, they underestimate the temperature of the small grains.This makes impossible the choice of a 'fit-all' dust opacity model to compute a realistic dust temperature.However, given the computational limitations of a dust opacity model using single-grain size opacities like Model M16, we argue that a satisfactory compromise is to use two (maybe three) independent size-averaged dust opacities to approximate the dust thermal spread of a realistic dust grain population, providing a good choice of composition and size range.-The grains of different sizes, having different temperatures and different coupling timescales to the gas phase, do not interact in the same way with the gas phase species, in particular when CO condensation/evaporation balance is considered.-The multiple temperatures significantly affects the CO gas phase distribution.Our model using two dust populations shows the presence of a CO depletion from 50 au up to 250 au, followed by a CO gas rise at larger radii.This rise does not require additional external UV irradiation nor a change of dust-to-gas mass ratio.The effect cannot be reproduced in chemical models that use a single dust temperature structure.-The CO ice distribution is also affected and the snowline is reshaped.The dust temperature spread creates multiple CO condensation fronts, inducing a radial snowline splitting.This does not require radial drift, although inward drift will affect the snowline even more.A more physical model would need a larger number of grain sizes and related temperature, implying a more complex 2D CO snowline.
-These results have potential implications for observations and planet formation.On one hand, the presence of CO holes followed by a rise in observed Class I/II disks has been recently reported.We speculate that the origin of these CO holes could be linked to the combination of dust temperature spreads and size-dependent dust radial segregation.On another hand, the CO snowline splitting has implications on the chemical evolution during planet formation.Indeed, the formation of COMs via CO surface processes is in turn also dust-size segregated such that it affects the way COMs are brought toward planetary cores.
Overall, this work illustrates the necessity for modelers to consider simultaneously multiple dust grain populations with independent optical properties, both for the thermal and chemical simulations.

Fig. 2 :
Fig.2: Absorption opacity and scattering albedo profiles for the 16 bin sizes in Model M16.The dotted, dashed, and solid lines represent the small, intermediate, and large grains, respectively.The dashed black vertical line indicates the wavelength λ max at which the star radiates according to Wien's law.The grain sizes are indicated in µm.Refer to Table4for description of bin sizes.
Figure. 4 shows the dust temperature structures of Model D computed with the thermal Monte Carlo simulation in RADMC3D.The two populations exhibit notable differences.In particular, the 25 K isothermal is located at ∼ 30 au from the star (in the midplane) for the large grains (Fig.4, right), which is a typical expected distance where the CO snowline is located in a T-Tauri disk.On the other hand, the small grain temperature drops below 25 K near 180 au only(Fig.4, left).There is therefore a significant temperature spread with radial distance of 150 au between the two 25 K isothermal.

Fig. 3 :
Fig. 3: Panels a and b: 2D dust area density [cm −1 ] of Model D. Panel c: ratio of the large grains area density to the small grains area density (note the different z scale).Panel d: mass ratio between the large and small grains.The green lines mark 1 × H g and 2 × H g .

Fig. 4 :Fig. 5 :
Fig. 4: 2D dust temperature structures of Model D. Left: temperature of the small grains.Right: temperature of the large grains.

Fig. 6 :Fig. 7 :
Fig. 6: Midplane radial dust temperature profiles of Model M16.The dotted curves represent the dust temperatures of the small dust grains, the dashed curves represent the temperature of the intermediate grains and the solid curves are the dust temperatures of the large dust grains.The grain sizes are given in µm.Left: temperatures with scattering.Right: temperatures with scattering opacity of the intermediate grains set to zero.

Fig. 8 :Fig. 9 :
Fig. 8: 2D maps of CO gas number density [cm −3 ].Left: CO in Model S. Middle: same in Model C. Right: same in Model D. The solid white line corresponds to the CO snowline where n g (CO)/(n gas (CO)+n ice (CO)) = 0.5.The black line is the maximum scale height below which chemistry is computed.

Fig
Fig. 10: 2D map of the disk.The grey scale represent the values of the optical depth computed from the star to all points in the disk τ ⋆ (r, z).The red lines represent the surface where τ ⋆ (r, z) = 1.The green lines mark the location where the total vertical optical depth at the midplane becomes smaller than 10, 5, and 1.The solid blue line marks the altitude at which the optical depth in the z-direction toward the upper layers becomes 1.The line diverges to infinity beyond the line τ v = 1.All opacities are given at 1 µm.
Fig. 12: 2D maps of CO ice number density [cm −3 ] in Model D. Left: CO ice adsorbed onto the small dust population only.Right: CO ice adsorbed onto the large dust population only.The solid grey line corresponds to the total CO snowline.The dashed red lines correspond to the snowline of each population independently.The black line is the maximum scale height below which chemistry is computed.

Table 2 :
Adopted elemental initial abundances relative to H.

Table 3 :
Summary of the dust disk models.The grain sizes in column 4 and the mass fraction in column 5 are calculated using Eq.11 and Eq. 12, respectively.