© S. Savvidou et al. 2020

1 Introduction

Protoplanetary discs surround young stars for the first few million years after their formation and they are the birthplaces of planetary systems. The position of the iceline within the discs influences the formation and growth of planets. Planetesimal formation has been found to be enhanced or even initiated there because of water vapor that is diffused outwards from the hot, inner disc and recondenses after the iceline (Ros & Johansen 2013). This recondensation increases the abundances of icy pebbles, which have better sticking properties compared to dry aggregates (Supulver et al. 1997; Wada et al. 2009; Gundlach & Blum 2015), causing a pile-up near the iceline and triggering the streaming instability (Guillot et al. 2014; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017).

An increase in the dust surface density after the iceline can also aid in the growth of gas giant planet cores (Stevenson & Lunine 1988). The location and the evolution of the iceline location can be defining for the innermost boundary of gas giant formation and along with other parameters, such as the disc’s mass, it can also determine what kind of planets will be created (Kennedy & Kenyon 2008) and their masses (Morbidelli et al. 2015, 2016). In addition to that, the location of the iceline transition affects the composition of exoplanetary atmospheres (Madhusudhan et al. 2014, 2017; Cridland et al. 2016).

The location of the iceline is determined by the local temperature in the disc (Hayashi 1981; Sasselov & Lecar 2000; Podolak & Zucker 2004). The thermal structure of the discs is thus decisive for planetesimal and planet formation. It is, though, greatly affected by the dust content of the protoplanetary disc and the opacity that the dust grains provide. This complex interplay is caused by the influence between gas and dust. The relative velocity for each pair of grains is determined by the aerodynamic properties of the grains, namely the Stokes number, and the local properties of the gas, such as the temperature or the volume density (Ormel & Cuzzi 2007). The variety in the relative velocities results in different collisional outcomes between grain sizes, such as coagulation or fragmentation (Brauer et al. 2008; Zsom et al. 2011; Birnstiel et al. 2011, 2012). As a result, the dust content of the protoplanetary disc is described by a distribution of grain sizes, with number densities that are not necessarily equally distributed between all existing sizes.

Each grain size population has a different opacity, therefore having a distribution instead of a single grain size means that the disc’s total opacity will be affected and as a consequence it will affect the resulting structure of the disc. As stated in Birnstiel et al. (2016), the opacity, defines the observational characteristics of a protoplanetary disc by influencing the dust thermal continuum emission and the excitation conditions for the gas lines. Additionally, since opacity regulates the amount of light that can be absorbed by the disc, it determines its thermal structure. These reasons make opacity a very important factor of the structure and evolution of a protoplanetary disc.

The interplay between opacity and the thermal structure creates a feedback loop that we include in hydrodynamical simulations of equilibrium discs (Fig. 1). Even though the goal of theoretical models is to simulate protoplanetary discs as realistically as possible, typically they only include specific parts of the feedback loop, contrary to this work. In the following paragraphs we introduce what work has been done in parts of the feedback loop.

Fig. 1 Graphical illustration of the feedback loop. The thermal and density structure of a protoplanetary disc is determined by this loop: temperature and gas surface density affect the relative velocities for grains of different sizes. Through the relative velocities we find the outcomes of collisions between grains, therefore a grain size distribution is created. The grains are then vertically distributed according to their sizes and the turbulence strength. The spatial distribution of the grain sizes determines the opacity of the disc, which then affects its cooling rate and the stellar heating. This way the temperatureand density of the disc change and subsequently its whole structure.

Relative velocities and grain size distribution

Early on, Safronov (1969) worked on dust growth within the context of planet formation and on the time evolution equation for grain size distributions (often called Smoluchowski equation, Smoluchowski 1916). A lot of work was also done on dust dynamics and how they would affect collisional outcomes and, as a consequence, coagulation and fragmentation of dust particles (Weidenschilling 1980, 1984; Nakagawa et al. 1981; Brauer et al. 2008). A grain size distribution has been widely assumed to follow a power-law derived from the equilibrium between coagulation and fragmentation, inspired by the work of Dohnanyi (1969) on the number density distribution of objects in the asteroid belt. The number density, thus, can be approximated as n(s) ∝ s^ξ, where s is the grain size and ξ a constant. Several attempts were made in order to define this constant, mainly through analytical calculations combined with observational data for the interstellar medium grains (e.g., Mathis et al. 1977, MRN power-law), but also through experimental studies (e.g., Davis & Ryan 1990). It was shown by Tanaka et al. (1996) that the ξ constant is independent of the specific parameters of the collisional outcome model, as long as it is self-similar, which in this case means that the outcome of impacts between dust grains depends on the masses of two colliding particles only through their ratio.

However, such a description of a grain size distribution with only one power law is a simplification, since it only takes into account the coagulation/fragmentation equilibrium. More recently, the work on grain size distributions has been aided by laboratory experiments of dust collisions (review by Blum & Wurm 2008; Güttler et al. 2010). Such experiments determine what the collisional outcomes are between particles of equal or different size, for different relative velocities. They also help in creating models to simulate such collisions accurately and they can be used in the effort of understanding which processes are relevant within the context of planetesimal formation in protoplanetary discs (e.g., Zsom et al. 2010). If additional effects are also taken into account, such as cratering or different regimes due to size-dependent relative velocities, then the size distribution is described by broken power laws (Birnstiel et al. 2011). The studies discussed above focused on the local distribution of grains in a protoplanetary disc patch due to fragmentation and growth by coagulation, and typically assume that the gas disc does not evolve in time and the dust has no effect on the gas.

Opacity

As a first step, some work has been done on opacity alone within the context of protoplanetary discs (Miyake & Nakagawa 1993; Draine 2006; Cuzzi et al. 2014). The goal of those works is to create a simple opacity model that can describe as realistically as possible the dust opacity and can be then used in disc simulations (Bitsch et al. 2013a) or help in the interpretation of disc observations (Birnstiel et al. 2018). Alongside the theoretical models, several observations of the dust emission have been performed in order to connect opacity with the particle sizes present in the protoplanetary discs (Natta et al. 2007; Andrews & Williams 2005, 2007; Rodmann et al. 2006; Lommen et al. 2009; Ricci et al. 2010, 2011; Ubach et al. 2012).

Disc structure and grain size distribution

Several works in the recent years aimed to couple the dust and gas components of protoplanetary discs in simulations and in most of the cases such models include a grain size distribution. However, the models that we discuss here simulate the gas component of a protoplanetary disc and how the dust component is affected by the gas, but the solids do not influence the gas. Even without the back-reactions of dust on gas, modeling grain size distributions can be computationally challenging, given the long list of effects and parameters to be taken into account, especially using N-body like techniques to treat dust particles. As a consequence, some of the first attempts on this kind of models were made using the Monte-Carlo method (Ormel et al. 2007; Ormel & Spaans 2008) and the goal was to examine how the internal structure of dust affects the collisional evolution of the particles and the disc structure. The Monte-Carlo method has been also used in Zsom et al. (2010, 2011), while in their work the experimental collisional outcomes from Güttler et al. (2010) were implemented and the effect of the porosity and settling of the dust grains on the collisional outcomes was tested. Brauer et al. (2008) and Birnstiel et al. (2010) numerically solve the Smoluchowski equation for the coagulation/fragmentation equilibrium in vertically isothermal steady-state gas discs, while Okuzumi et al. (2012) studied the effects of the dust grain porosity on the dust evolution in a similar disc setup. In the works discussed above the feedback of the dust on the gas disc structure and especially its thermal part, is not taken into account.

Disc structure, opacity and cooling rate

The category of models that was described above neglected the effects of opacity, even though the dust opacity regulates the cooling rate of the disc, which affects the disc structure. In recent years, some studies tried to fill this gap by including the effect of dust opacity in disc simulations. Oka et al. (2011) performed 1+1D¹ simulation focusing on the effect that water-ice opacity has on the location of the iceline. In the aforementioned study the wavelength-dependent opacities of water-ice and silicates are directly used when calculating the radiative transfer. In Bitsch et al. (2013a, 2014, 2015a) the Bell & Lin (1994) opacity profile is followed (in Bitsch et al. 2013a constant opacity discs were also modeled) and 2D simulations (radial and vertical direction, assuming axisymmetry) are performed using the NIRVANA and the FARGOCA code adding radial heat diffusion and stellar irradiation. The effect of the water-ice to silicates ratio on the resulting thermal disc structures has also been studied recently (Bitsch & Johansen 2016) using the FARGOCA code and the opacity module from the RADMC-3D code to calculate the mean opacities (as in the present work), but the opacity differences for the water-to-silicate fractions considered are then translated into differences in the Bell & Lin (1994) opacity model. The Bell & Lin (1994) opacity model gives approximate values for the frequency averaged opacities within specific temperatureregimes (e.g., ice grains, evaporation of ice grains, metal grains, etc.) assuming micrometer-sized particles. The fixed opacity profile then gives the cooling rate and the stellar heating, therefore it defines the disc structure.

Even though including the opacity feedback in disc simulations is an important improvement, the aforementioned studies did not include the effect of grain growth and fragmentation, and thus, only employed opacities derived for single grain sizes. In addition to this, all of these studies assumed a uniform dust-to-gas ratio in the vertical direction of the disc which is in contrast to our approach in this work (see Sect. 2.4).

Disc structure, grain size distribution and opacity

Schmitt et al. (1997) coupled the dust and gas evolution in 1D simulations, while they also took into consideration the grain opacity. For the mean opacity calculations they followed the approach of Henning & Stognienko (1996), which is similar to the Bell & Lin (1994) opacity model approach. Moreover, the size distribution follows the Mathis et al. (1977) power-law. It was found that since grains determine the opacity, their evolution will subsequently change the opacity and therefore affect the structure and evolution of a protoplanetary disc. Prior to this study, Mizuno et al. (1988) and Mizuno (1989) included the dust component evolution in accretion discs and used the results to perform grain opacity calculations. In Suttner & Yorke (2001) the coagulation/fragmentation equilibrium is included in order to investigate how the dust emission is affected by the grain size distribution and its corresponding opacity. In this work the size distribution follows the Mathis et al. (1977) power-law and opacity was calculated using Mie theory. However, in the studies discussed above the back-reaction of the opacity onto the disc structure was not taken into account.

In the previous paragraphs some examples were given of the work that has been done in the context of grain growth within protoplanetary discs. Nevertheless, previous models were based on several simplifications, most important of which was that they neglected parts of the feedback loop (Fig. 1) that defines protoplanetary disc structures (e.g., Birnstiel et al. 2011) or used simplified assumptions for the opacity (e.g., Bitsch et al. 2013a). The few attempts that have been made to include the dust feedback on the gas of the disc, were 1D simulations or assumed an isothermal vertical structure for the gas, in contrast to the 2D hydrodynamical models that are presented here. Secondly, the opacity is either not included in the actual simulations or the opacities were included only for single fixed grain sizes.

The motivation for this project is to approach a more realistic model for disc structures and their evolution and more specifically to simulate the whole feedback loop including a detailed opacity module. We consider how grain dynamics and more specifically how grain size distributions affect the opacity and as a consequence the thermal structure of the disc in order to simulate the whole feedback loop. As far as the grain size distribution is concerned, two models were used for the simulations of this project. A simple power-law model following Mathis et al. (1977), hereafter MRN distribution and also a more complex model following Birnstiel et al. (2011), hereafter BOD distribution. Moreover, an opacity module was included in the 2D hydrodynamical disc simulations in order to more accurately calculate the opacity of the dust grain distribution and account for the back-reactions of dust to gas. In this opacity module, the Rosseland and Planck mean opacities as a function of temperature are used and they are calculated via Mie theory. The simulations were run until the disc reached thermal equilibrium. Such simulations offer us the opportunity to discuss the implications of the resulting disc structures to planet formation and could also serve as the basis to compare with observations (e.g., ALMA images) in future work.

The structure of this paper is as follows: in Sect. 2 we describe the energy equations used in the hydrodynamical simulations, the new opacity module and the two grain size distributions that we included in the code. Also, we mention the input parameters of the simulations and how the disc is set up. Then, in Sect. 3 we compare the discs utilizing the two different grain size distributions for one set of simulation with a nominal turbulence parameter and initial gas surface density. In Sect. 4, we discuss how the simulations change when we vary the turbulence parameter and how it changes when the surface density is different. The implications of the results are discussed in Sect. 5 and a summary follows in Sect. 6.

2 Methods

In this section, we discuss the different methods used in our work. We review the hydrodynamical equations in Sect. 2.1, the opacities, and how they are calculated in Sect. 2.2. In Sect. 2.3 we discuss the grain growth mechanism and in Sect. 2.4 we discuss the effects of vertical grain settling. We then finally describe our simulation setup in Sect. 2.5.

2.1 Hydrodynamical simulations

Calculations with mean opacities derived from single grain sizes were first introduced into the FARGOCA code by Lega et al. (2014) and Bitsch et al. (2014), who performed 2D and 3D radiation hydrodynamical simulations of discs and planet-disc interactions. The FARGOCA code solves the continuity and the Navier-Stokes equations, and uses the flux-limited diffusion approach to radiative transfer. More specifically, the time evolution of the energy profile of the protoplanetary disc is determined by $$\frac{\partial E_{\text{R}}}{\partial t}+\nabla \cdot F=\rho {\kappa }_{\text{P}}[B(T)-c{E}_{\text{R}}]$$(1) $$\frac{\partial ϵ}{\partial t}+\nabla \cdot (u\cdot \nabla )ϵ=-P\nabla \cdot u-\rho {\kappa }_{\text{P}}[B(T)-c{E}_{\text{R}}]+Q^{+}+S.$$(2)

The radiative energy density E_R is thus independent from the thermal energy density ϵ. In the expressions above the blackbody radiation energy is B(T) = 4σT⁴, where σ is the Stefan-Boltzmann constant, ρ is the gas density, κ_P the Planck mean opacity (further specified in Sect. 2.2), u the velocity, P is the thermal pressure, Q⁺ is the viscous dissipation or heating function and S is the stellar heating component (Levermore & Pomraning 1981; Dobbs-Dixon et al. 2010; Commerçon et al. 2011).

In our simulations we use the flux-limited diffusion (FLD) for the radiation flux F as described in Levermore & Pomraning (1981) $$F=-\frac{\lambda c}{\rho {\kappa }_{\text{R}}}\nabla E_{\text{R}}.$$(3)

In the flux-limited diffusion equation, c is the speed of light, α_R is the radiation constant, κ_R is the Rosseland mean opacity and λ the flux-limiter of Kley (1989). More details on the energy equations can be found in Bitsch et al. (2013a). The opacities that were introduced in the above equations are discussed in the following section.

The stellar heating density received by a grid cell of width Δr is defined as (Dobbs-Dixon et al. 2010): $$S=F_{\star }{\text{e}}^{-\tau }\frac{1-{\text{e}}^{-\rho {\kappa }_{\star }\Delta r}}{\Delta r},$$(4)

with $F_{\star }=R_{\star }^2\sigma T_{\star }^4/r^2$ being the stellar flux, R_⋆ the stellar radius, T_⋆ the stellar surface temperature, τ the radially integrated optical depth (up to each grid cell) and κ_⋆ the stellar opacity (further specified in Sect. 2.2).

Fig. 2 Rosseland, Planck, and stellar mean opacities (from left to right) as a function of temperature for grains of sizes 0.1, 1, 10, 100 μm, 1 mm, and 1 cm. They are independent of the gas density because they are dominated by the dust component. These values were calculated using RADMC-3D for a mixture of 50% silicates, 50% ice, and disc dust-to-gas ratio of 1%. The gray vertical line shows the location of the water iceline transition (170 K ± 10 K), causing a transition in opacity due to the evaporation/condensation of dust grains.

2.2 Opacity-temperature module

In the energy equations (Eqs. (1) and (2)) of the hydrodynamical simulation, we use the mean opacities that are averaged over all wavelengths. If we use the Planck black body radiation energy density distribution B_λ (λ, T) as a weighting function we can define the Planck mean opacity as $${\kappa }_{\text{P}}=\frac{{{\displaystyle \int }}_0^{\infty }{\kappa }_{\lambda ,ns}(T,\rho )B_{\lambda }(\lambda ,T)\text{d}\lambda }{{{\displaystyle \int }}_0^{\infty }{B}_{\lambda }(\lambda ,T)\text{d}\lambda }.$$(5)

Since the mean free path of thermal radiation in the disc is small compared to the disc’s scale height, the radiation field can be considered isotropic, blackbody emission.

The Rosseland mean opacity uses the temperature derivative of the Planck distribution as a weighting function and is defined as $${\kappa }_R^{-1}=\frac{{{\displaystyle \int }}_0^{\infty }{\kappa }_{\lambda ,s}^{-1}(T,\rho )(\partial B_{\lambda }(\lambda ,T)/\partial T)\text{d}\lambda }{{{\displaystyle \int }}_0^{\infty }(\partial B_{\lambda }(\lambda ,T)/\partial T)\text{d}\lambda }.$$(6)

It should be noted that scattering processes are neglected (subscript ns) when calculating the wavelength dependent opacities κ_λ for the Planck mean, but are included in the Rosseland mean opacity (subscript s).

We also consider the stellar radiation and define the stellar or optical opacity as $${\kappa }_{\star }=\frac{{{\displaystyle \int }}_0^{\infty }{\kappa }_{\lambda ,ns}(T,\rho )B_{\lambda }(\lambda ,T_{\star })\text{d}\lambda }{{{\displaystyle \int }}_0^{\infty }{B}_{\lambda }(\lambda ,T_{\star })\text{d}\lambda }.$$(7)

The stellar opacity is then the Planck mean opacity taking into consideration the stellar radiation temperature instead of the local disc temperature.

We calculate the mean Rosseland, Planck, and stellar opacities as a function of temperature using the RADMC-3D² code. We note that dust opacities are independent of the gas density, as opposed to gas opacities. The latter are not considered in this work as opacities in the disc are dominated by the dust component and the high temperature needed for dust evaporation is not reached within our simulations. The code utilizes Mie-scattering theory and the optical constants for water-ice (Warren & Brandt 2008) and silicates (Jaeger et al. 1994; Dorschner et al. 1995) in order to calculate the wavelength-dependent opacities, which are then averaged over all wavelengths. The main input parameters are the size of the grains and the dust-to-gas ratio of the disc. We can also choose the dust grain species, silicates, water ice, and carbon or the fraction between those in the dust mixture. In this work, we include a mixture of 50% water-ice and 50% silicates and the dust-to-gas ratio for the calculation of the opacities is 1%. Finally, the Rosseland, Planck and stellar mean opacities (see Bitsch et al. 2013a) are calculated.

In Fig. 2 we illustrate how each mean opacity scales with temperature for six different grain sizes, from 0.1 μm to 1 cm. The wavelength-dependent opacities and subsequently the mean opacities depend on the size parameter $x=\frac{2\pi s}{\lambda }$, but also on the refractive index of the given grain species, which is also itself dependent on wavelength (Movshovitz & Podolak 2008). By Wien’s law, the wavelength is inversely proportional to the temperature. Using the size parameter we find that the regime changes at approximately x = 1 and more specifically at x ≪ 1 we have the Rayleigh scattering, whereas at x ≫ 1 we have thegeometric optics regime (Bohren & Huffman 1998). Consequently, if the size of the particle is a lot smaller than the wavelength of the incident radiation, absorption dominates over scattering and the wavelength dependent opacities become independent of grain size. In the case of the larger grain sizes, or when x ≫ 1, the opacities become independent of wavelength (and consequently temperature) but depend on the grain size. Most of the regions though lie somewhere in between, which means that calculating the opacity depends on both the grain size with its individual refractive index and the given wavelength or temperature.

The Rosseland mean opacities (Fig. 2) for the largest particles of the set (100 μm, 1 mm and 1 cm) are almost flat, except for the transition around the iceline at 170 K ± 10 K. At this temperature, ice sublimates, and the opacity is then only determined by silicates. For those large particle sizes, the size parameter is greater than 1, therefore we are in the geometric optics regime and the Rosseland mean opacity is independent of temperature. However, we note that for a grain size of 100 μm and temperaturebelow 20 K the opacity depends on the temperature. In this region, the regime has changed and the opacity is determined by Rayleigh scattering. Equally, the size parameter is well below 1. The same trend can be seen for the smaller particles, namely 0.1, 1 and 10 μm before the iceline. The opacity of the 10 μm grain sizes goes into the geometric optics regime after the iceline and tends to become independent of temperature. In the region after the iceline for the smallest particles (0.1 and 1 μm) the size parameter is closer to 1, so the opacities are also influenced by the refractive indices.

The Planck opacities have a weaker dependency on temperature compared to the Rosseland mean opacities. The stellar opacities depend only on the stellar temperature and grain sizes, but not on the disc temperature, except for the transition at the water iceline, when water-rich particles evaporate. Both the Planck and stellar opacities are calculated using only the absorption coefficient, which does not have a strong dependency on wavelength and consequently temperature, as opposed to the extinction coefficient. The Planck mean opacities are calculated taking into account the temperature of the disc, while the stellar opacities, use the temperature of the star, which is constant.

Using RADMC-3D several files are created with the mean opacity values as a function of temperature. These files are then used in the hydrodynamical code (FARGOCA). The opacity calculations from these files are interpolated and in this way, we get in the code the appropriate opacity values given the temperature of the grid cell. The reason why the interpolation is done instead of directly calculating opacity using Mie theory is that the computational time would be very long. We include the direct opacity-temperature calculations for at least ten grain sizes, from 0.1 μm to 1 mm or 2 cm. We then create size bins and as a simplification, each grain size within a bin shares the same opacity-temperature calculations (corresponding to the logarithmic mean size of that bin). We note here that these bins are different and for computational reasons wider than the bins used for the calculations of the vertically integrated dust surface densities (see Sect. 2.3).

As a comparison, we also use the frequency-averaged Bell & Lin (1994) opacity law. The opacity, in this case, depends on the local temperature and density. There are several transitions in this opacity regime caused by the processes which dominate each temperature region, such as the evaporation of ices interior to the iceline, which is also present in the prescription we are using for the discs with the grain size distributions. However, the greatest difference between the two opacity regimes is that the Bell & Lin (1994) opacity law is based on micrometer-sized dust and does not take the opacity provided by all of the dust sizes present in the disc into account. The Bell & Lin (1994) law also considers the gas opacities, but these are relevant for high temperatures that are not reached in the simulations presented here.

2.3 Grain size distributions

The collision between two dust grains can result in various outcomes. The possible outcomes are coagulation, fragmentation, cratering, and bouncing (review by Blum & Wurm 2008). The outcome of a collision is determined by the relative velocities of the collidingbodies and their mass ratio (Weidenschilling 1977; Brauer et al. 2008).

The relative velocities between grains are determined by the mass of the particles, but they are also greatly affected by the local temperature and the gas scale height. Dust dynamics involve not only collisions between grains, but also with the molecules of the protoplanetary disc’s gas. These collisions with the gas cause a lag to the dust particles that leads to relative velocities between themselves.

We compare in this work two different grain size distribution models. These models provide the vertically integrated surface density of dust as a function of the grain size. The first and simple model (hereafter MRN) is inspired by the groundwork on dust distributions (Dohnanyi 1969; Mathis et al. 1977; Tanaka et al. 1996). At a given distance to the star, the equilibrium between fragmentation and coagulation results in a steady-state size distribution, where the number density of the particles can be written as $$n(m)\text{d}m\propto m^{-\xi }\text{d}m$$(8)

or $$n(s)\text{d}s\propto s^{2-3\xi }\text{d}s,$$(9)

with m the particle mass, s the particle size, and ξ a constant.

The mass of a specific size, within a size bin [s_i −ds′, s_i + ds″] is $$M_{s_i}={\displaystyle {\int }_{s_i-\text{d}{s}^{\prime }}^{s_i+\text{d}{s}^{″}}m}\cdot n(s)\text{d}s\propto {\left[\frac{s^{6-3\xi }}{6-3\xi }\right]}_{s_i-\text{d}{s}^{\prime }}^{s_i+\text{d}{s}^{″}},$$(10)

assuming 5-3ξ≠ − 1. The grain sizes for this project are distributed over a logarithmic grid, so s_i − ds′ is $\sqrt{s_i\cdot s_{i-1}}$ and s_i + ds″ is $\sqrt{s_i\cdot s_{i+1}}$. The vertically integrated surface density of each grain size bin is then $${\Sigma }_{\text{d},s_i}\propto f_{\text{DG}}{\Sigma }_{\text{g}}\left[{\left(\sqrt{s_i\cdot s_{i+1}}\right)}^{6-3\xi }-{\left(\sqrt{s_i\cdot s_{i-1}}\right)}^{6-3\xi }\right],$$(11)

where f_DG is the dust-to-gas ratio and Σ_g is the gas surface density.

We use a grain size grid, such as s_i+1 = c ⋅ s_i and the assumption that ξ = 11/6 (Dohnanyi 1969; Williams & Wetherill 1994), so then the expression for the unnormalized vertically integrated surface density for each grain size bin can be simplified to $${\Sigma }_{\text{d},s_i}\propto s_i^{1/2}{f}_{\text{DG}}{\Sigma }_{\text{g}}.$$(12)

The contributions from each grain size are then summed up. In order to get the normalized surface density values we divide each contribution by the aforementioned sum $${\tilde{\Sigma }}_{\text{d},s_i}=\frac{s_i^{1/2}{f}_{\text{DG}}{\Sigma }_{\text{g}}}{{\sum }_i{s}_i^{1/2}}.$$(13)

The second and more complex model (Birnstiel et al. 2011, hereafter named BOD) takes into account fragmentation, coagulation, and also cratering, where only part of the mass of the target body is excavated after the collision with a small impactor. The input parameters for this model are the dust and gas surface densities (Σ_d,0 and Σ_g,0), the local disc temperature (T), the alpha turbulence parameter (α), the volume density of the particles (ρ_s) and finally the fragmentation velocity (u_f), which is the critical velocity above which all collisions lead to either fragmentation or cratering. The logarithmic grid for the sizes of both distributions is defined as s_i+1 = 1.12s_i, while the smallest grain size is 0.025 μm. As mentioned in Sect. 2.2 this grain size grid is finer than the size grid we use to determine the opacities.

Considering that different particle sizes lead to different collision outcomes, this recipe takes into account the relative velocities that particles of different sizes will develop in order to create different regimes for each size. These regimes are created according tosize boundaries, within which different power-laws apply for the fit to the size distribution. It should be mentioned that these size boundaries are defined by the corresponding relative velocities of the dust grains. The smallest particles of the distribution follow Brownian motion, which means their motion is affected by collisions with the gas molecules, there is no preferred direction and they do not have angular momentum. The next regime regards larger particles that start to get affected by turbulent mixing. It was also found (Ormel & Cuzzi 2007) that when particles have stopping times approximately equal or larger compared to the turn-over time of the smallest eddy of the gas, they start to decouple from the gas, so they follow a different regime. Finally, the distribution has an upper end or a fragmentation barrier above which particles can no longer grow and only fragmentation occurs.

Between two size boundaries, the distribution is described by a power-law $n(m)\cdot m\cdot s=s_i^{{\delta }_i}$ of different powers δ_i, depending on the regime (Brownian motion or turbulent mixing). Within each one of the regimes, the power-law indices are different if the grains are affected by settling, given their sizes and the disc parameters (see Sect. 2.4 for a discussion on the vertical distribution of grains). The powers for each regime are found in Table 1 and using these we can create a first fit f(s_i). It is necessary then to include a bump caused by cratering and the cut-off effects of the distribution that cause an increase in the fit for large enough particles. This boundary effect is caused by the fact that large particles near the upper boundary of the grain sizes grid do not have larger particles to collide with, but the mass transfer from one size bin to the other needs to be constant to keep a steady-state grain size distribution. Therefore the number density is increased. Similarly, erosion by small impactors slows down the growth of large particles, and an increase in the number density is needed to keep the flux constant. More details for this recipe can be found in Birnstiel et al. (2011).

Finally, the fit is normalized according to the total dust surface density at the given location (Fig. 3, also see Sect. 5.2 in Birnstiel et al. 2011) as in the first model. This fit represents the vertically integrated dust surface densities per logarithmic bin of grain size, N(s) ⋅ m ⋅ Δ logs, where

Fig. 3 Vertically integrated dust surface density distribution per logarithmic bin of grain size as a function of grain size for the two distributions used here, after Birnstiel et al. (2011; BOD) and Mathis et al. (1977; MRN), at 10 AU, for the simulation with α = 5 × 10−3 (see Fig. A.1 for the distributions over all orbital distances). The dust to gas ratio is 1% and the gas surface density is 1000 g cm−2 at 1 AU in both cases. For both of the distributions we additionally used uf = 1 m s−1 and ρs = 1.6 g cm−3. In the BOD, we can distinguish three regions. Small particles follow Brownian motion (BM), then as they grow they followthe turbulent gas motions (T1), and finally they are affected by the turbulence, but get decoupled from the gas as their stopping times are much larger than the turn-over time of the eddies (T2). The barrier sBT is the grain size limit for Brownian motion, while the barrier s12 separates the two turbulent regimes. The bump near the end of the distribution is caused by cratering, since large particles only lose part of their mass this way, while small particles can only coagulate and form larger grains. This causes the distribution to be top-heavy. The MRN distribution is a single power-law (Eq. (13)). Both grain size distributions have a maximum value determined by the same fragmentation limit (Eq. (15)), but they are not the same for the BOD and the MRN distribution because of the self-consistently calculated temperature.

$$N(s)={\displaystyle {\int }_0^{z_{\text{max}}}n}(s) \text{d}z$$(14)

is the vertically integrated number density.

The maximum grain size or in other words the fragmentation barrier is defined in the simulations with both of the grain size distributions as $$s_{\text{max}}\simeq \frac{2{\Sigma }_{\text{g}}}{\pi \alpha {\rho }_s}\frac{u_{\text{f}}^2}{c_s^2},$$(15)

with ρ_s = 1.6 g cm⁻³* the density of each particle, u_f = 1 m s⁻¹ the fragmentation threshold velocity, k_B the Boltzmann constant, m_p the proton mass and μ = 2.3 the mean molecular weight in proton masses. The sound speed is given by $$c_s=\sqrt{\frac{k_{\text{B}}T}{\mu m_{\text{p}}}}.$$(16)

The threshold velocity u_f ~ 1 m s⁻¹ corresponds to the threshold after which collisions between silicates always lead to fragmentation (Poppe et al. 2000). However, it has also been experimentally found that water-ice shows a higher threshold velocity, u_f ~ 10 m s⁻¹ (Gundlach & Blum 2015). We choose to use only the lower fragmentation threshold in the here presented work, but the composition dependency will be studied in future work.

Because of Eq. (15), which applies to both distributions, and the different regime boundaries in BOD, which depend on the local disc parameters (see Fig. 3 and Table 1), there is not a global size distribution, but rather a self-consistent spatial distribution of grain sizes both radially and vertically (see also Sect. 2.4).

It is noteworthy that even though we consider the coagulation/fragmentation equilibrium and the effects of cratering and settling, we neglect in the following work the drift of grains and the effect of bouncing. However, in the simulations presented here we find that the fragmentation barrier is always smaller than the drift barrier. This means that the particles have already fragmented and replenished the smaller pieces before they would have the chance to experience drift. The small particles are less affected by radial drift (Weidenschilling 1977) and since they coagulate, an equilibrium forms that drives the grain size distribution. The fragmentation barrier decreases with increasing α-viscosity parameter, which is expected since an increased α leads to increased turbulent relative velocities. The maximum possible grain size also decreases when the fragmentation threshold velocity decreases. Drift is an important effect acting on dust grains in protoplanetary discs, but it is a reasonable simplification to neglect it for the chosen parameters of our simulations. In future work where, for example, the fragmentation threshold velocities are increased or the composition dependency is included, drift is an effect that needs to be taken into consideration.

2.4 Vertical distribution of grains

The grains of a given size are vertically distributed according to the following $${\rho }_{\text{d}}={\rho }_{\text{d},0}\exp \left(-\frac{z^2}{2H_{\text{d}}^2}\right),$$(17)

with $${\rho }_{\text{d},0}=\frac{{\Sigma }_{\text{d}}}{\sqrt{2\pi }{H}_{\text{d}}}$$(18)

the dust density at midplane and the dust scale height derived by Dubrulle et al. (1995) $$H_{\text{d}}=H_{\text{g}}\sqrt{\frac{\alpha }{\alpha +\text{St}}},$$(19)

where H_d and H_g is the dust and gas scale height respectively. The Stokes number of the particles in the Epstein regime is $$\text{St}={\tau }_{\text{f}}{\Omega }_{\text{K}}=\sqrt{\frac{\pi }{8}}\frac{{\rho }_{s}s}{{\rho }_{\text{g}}{H}_{\text{g}}},$$(20)

with τ_f the stopping time of the particles, Ω_K the Keplerian velocity, andρ_g the gas volume density. At midplane, where ρ_g = ρ_g,0 (resembling Eq. (18)), the Stokes number is $$\text{St}=\frac{\pi }{2}\frac{{\rho }_{s}s}{{\Sigma }_{\text{g}}}.$$(21)

The vertically integrated dust surface densities Σ_d as a function of orbital distance are determined by the grain growth and fragmentation equilibrium prescriptions that were introduced in Sect. 2.3. The BOD grain size distribution has already taken into account the effect of settling (to calculate the distribution itself), depending on the grain sizes and the disc parameters (see Sect. 2.3 and Table 1). We then distribute in our model the grains vertically according to their sizes and how much they are expected to be affected, in a fashion consistent with the assumptions made in BOD (Eqs. (17)–(21)). However, it has been shown that small particles can get trapped in lower altitudes by the concentration of larger grains due to settling (Krijt & Ciesla 2016). This effect is not taken into account here as it is beyond the scope of this work but could be an improvement in future work.

We use the volume density of dust within a grid cell to find the opacity through $$\overline{\kappa }={\sum }_i\left(\frac{{\rho }_{\text{d},i}}{{\rho }_{\text{g}}}100\right) {\kappa }_i,$$(22)

where κ_i is the opacity of each grain size i, as shown in Fig. 2, and ρ_d,i ∕ρ_g the dust-to-gas volume density ratio for a given grain size i. In the case of single grain sizes summing is not needed. The dust-to-gas term for the volume densities includes the settling effect (Eq. (17)). In this expression, we multiply the volume density dust-to-gas ratios by 100 to account for the fact that in the calculations of κ_i (with the module from RADMC-3D) a dust-to-gas ratio of 1% was assumed. This way we multiply this κ_i with the appropriate factor depending on the volume density dust-to-gas ratio.

As an example of the effect of settling, we show in Fig. 4 the dust density as a function of height at 3 AU for 5 different grain sizes. In this simulation, the α value is 10⁻⁴, the initial gas surface density at 1 AU is Σ₀ = 1000 g cm⁻² and the grain size distribution used is the BOD. As a reference, we also plot the gas density to indicate the different volume density dust-to-gas ratios depending on the grain size and the volume density dust-to-gas ratio.

Considering the α prescription for viscosity of Shakura & Sunyaev (1973) the turbulence must be α ≤ St for settling to become important (Armitage 2009). We see from Eq. (19) that the larger the Stokes number of a particle, the more it will be affected by settling for a given α. Additionally, the lower the α value is, the more effective settling will be for even smaller dust particles. For this reason, we choose to show an example in Fig. 4 of a simulation with α = 10⁻⁴ which is the lowest α value we used in our simulations.

Not only the size but the location within the disc matters, because the Stokes number for a given grain size depends on the gas density which decreases with the increasing orbital distance. As a consequence, the same particles experience less settling the closer they are to the inner boundary of the disc.

The smallest particles which are shown in Fig. 4, namely 1 and 10 μm are not affected by settling, despite the low turbulence strength. Their dust-to-gas ratio remains the same at all heights, so they are well coupled to the gas. Then, the larger the particle, the more effective settling is. The 100 μm sized dust particles are already affected by settling, but beyond this grain size the difference is even larger. The cm-sized dust, which is nearly the maximum grain size in this simulation, is almost constrained at the midplane. The dust-to-gas ratio (dashed gray line), ρ_d ∕ρ_g, is well below 1% above z = 0.05 AU, but reaches 4% at midplane.

The main difference between various grain sizes is their different opacities as a function of temperature. However, as illustrated in Fig. 4, settling is another important effect, with a distinct efficiency depending on the grain size. Test simulations (presented in Appendix C) show that a significant settling changes significantly the thermal structure of the disk. Indeed, without it a constant dust-to-gas ratio leads to overestimated opacities above midplane, hence more “puffed-up” inner discs with higher temperatures, which cause a shadowing of the outermost region and prevent it from reaching an equilibrium state. Thus, settling is an important effect that needs to be taken into account in models in order to accurately study the thermal structures of protoplanetary discs.

Fig. 4 Dust density as a function of height for grains of five representative grain sizes within a disc with the BOD, α = 10−4 and Σg,0 = 1000 g cm−2, at 3 AU. The black line shows the gas density of the disc. The dashed gray line shows the dust-to-gas ratio as a function of height for the whole grain size distribution. The dust-to-gas ratio of the smallest particles of the sample remains constant, but as grains get larger they get more affected by settling. Consequently, the largest grain (1 cm) is mostly concentrated at the midplane.

2.5 Simulations setup

The stellar mass used in the simulations is M_⋆ = 1 M_⊙, the temperature is 4370 K, and the radius is R_⋆ = 1.5 R_⊙. The total dust-to-gas ratio is f_DG = 1%. Viscosity in the simulations follows an α prescription (Shakura & Sunyaev 1973), where $$\nu =\alpha \frac{c_s^2}{{\Omega }_{\text{K}}}.$$(23)

Recent observations of protoplanetary discs find α values from 10⁻⁴ to 10⁻² or even 10⁻¹ (Hueso & Guillot 2005; Andrews et al. 2009; Rafikov 2016; Ansdell et al. 2018), but such large α would cause discs to rapidly expand to great extent (Hartmann et al. 1998) in contrast to observations. We use in this work five sets of simulations with α = 10⁻², 5 × 10⁻³, 10⁻³, 5 × 10⁻⁴, and 10⁻⁴ in order to testthe effect of turbulence on the thermal structure of the disc. We choose these values in order to include a simulation with α = 5 × 10⁻³ so that one can directly compare with the work in Bitsch et al. (2015a) and a simulation with α =10⁻³ to allow comparison with the discs in Bitsch et al. (2013a). In the simulations with α =10⁻² and 5 × 10⁻³ the grid cells are 480 × 70 (radial-vertical direction) and the disc extends from 2 to 50 AU, while in the simulations with 10⁻³, 5 × 10⁻⁴ and 10⁻⁴ the grid cells are 150 × 35 and the discextends from 0.1 to 3.1 AU.

The gas surface density follows a profile $${\Sigma }_{\text{g}}={\Sigma }_{\text{g},\text{0}}\cdot {(r/\text{AU})}^{-p},$$(24)

with p = 1∕2 and we test two different initial surface densities, Σ_g,0 = 100 and 1000 g cm⁻² for every α value that was mentioned above. We run more combinations of different initial gas surface densities and total dust-to-gas ratios, however these are mainly used in order to produce a fitting for the iceline position as a function of the three parameters, α, Σ_g,0 and f_DG (see Sect. 5.1, Appendices A and B).

Since we simulate equilibrium discs, the surface density profile does not evolve significantly during the simulation, because the thermal equilibrium is reached faster than the viscous evolution equilibrium. At the top of the disc we manually set T = 3 K, the temperature of the interstellar medium, so that the disc can be cooled by the upper boundary (as described in Bitsch et al. 2013a). The simulations run for some hundreds of orbits (typically 200–1000 orbits) until they reach thermal equilibrium. Nevertheless, some of the simulations might show signs of convection (Bitsch et al. 2013b), which means that they will remain unstable regardless of integration time.

At first, we perform simulations with single grain sizes in order to see the difference in the disc structures between them. Dust grains affect the hydrodynamical simulation through the opacity in each grid cell. Every simulation has a different grain size and the opacity values for this specific size are used (see Fig. 2). The simulations of single sizes offer the chance to examine the extent to which different grain sizes affect the disc’s evolution and equilibrium structure and predict how much grain growth or a grain size distribution will change the outcome.

In the next step, we also consider settling and how it affects large grains. For these simulations we only use single grain sizes and the dust surface density is assumed to be Σ_d = f_DGΣ_g or specifically Σ_d = 0.01Σ_g as before. The difference between discs without and with settling is further discussed in Appendix C.

We, then, include the two grain size distribution models that were discussed in Sect. 2.3. The distributions are self-consistently calculated in the code using as input parameters in each time-step and grid cell, the gas surface density (for both distributions) and the temperature (only for the BOD). The upper size boundary for the MRN power law model can be either fixed or follow the same fragmentation barrier formula as the second, more complex model (Birnstiel et al. 2011), but in this work, we use the latter.

3 Grain size distributions

3.1 Comparison between the two grain size distributions

In this section we compare the simulations utilizing the two different grain size distributions (MRN and BOD). At first, we focus only on the case of α = 5 × 10⁻³, Σ_g,0 = 1000 g cm⁻² and f_DG = 1%.

In Fig. 5 we can see the aspect ratio as a function of orbital distance from the star for simulations of different grain sizes along with the two grain size distributions. The gray areas correspond to the iceline transition (T = 170 ± 10 K). In this area, the change in opacity is responsible for the bumps in the aspect ratios. We first focus on the simulation with 0.1 μm, which roughly corresponds to an unevolved dust population as found in the interstellar medium. The simulation with particles of 0.1 μm results in an increasing aspect ratio as a function of orbital distance up to 6 AU, where it reaches a maximum and then decreases up to approximately 15 AU. Using 1 μm-sized particles wesee a similar disc structure. The aspect ratio is almost constant for the first few AU and has a small increase around 6 AU. Then it converges with the simulation of 0.1 μm up to the minimum around 15 AU. If the grain size increases to 10 μm, then the aspect ratio also increases with distance, features a bump closer to 7.5 AU and decreases with the same slope as the previous two simulations. The larger particles have distinct profiles. Specifically, the aspect ratio of the simulation with particles of 0.1 mm is a monotonically increasing function of orbital distance with a small bump at 3.5 AU. The same can be seen for the simulation with the largest particles, namely 1 mm, but in this case, a bump is not visible at any part of the aspect ratio profile since the iceline transition does not exist within the boundaries of the simulation.

The gradients in the aspect ratios for the inner region of the discs are determined by the opacity of the disc. We can compare the opacity gradients in Fig. 2 with the aspect ratio gradients keeping in mind that the temperature decreases as we move further away from the star. Depending on the temperature at each orbital distance of the disc, we see that the opacity gradients are responsible for the dips and bumps in the aspect ratio profiles of the discs.

The outer region of the discs is dominated by stellar irradiation, which causes the flaring of the discs. The simulation with the BOD shows influence from the smaller particles at the inner parts of the disc, but moving outwards the aspect ratio gets affected by larger particles, around 100 μm. A more detailed analysis of the grain sizes that contribute the most to the opacity of the disc is presented in Sect. 3.4. The simulation with the MRN distribution shows similar aspect ratio gradients. For this case, with high α, there is onlya minimal difference between the dust surface densities, which leads to similar opacities in total. Both discs are affected by grains of similar sizes and for this reason, the aspect ratios are almost the same there.

In Fig. 5b we show the midplane temperature as a function of orbital distance for the same set of simulations. We can see that the changes in the temperature gradients correspond to the changes in the aspect ratio gradients. The gray horizontal line is again the iceline transition. The simulation with the millimeter-sized particles does not reach the iceline temperature within the extent of the simulated disc and thus does not feature the aspect ratio bump.

We can compare in Fig. 6 the resulting aspect ratios and midplane temperatures of the simulations with the BOD and MRN distributions with a simulation utilizing the opacities from Bell & Lin (1994). This opacity profile is based on micrometer-sized particles, hence it is expected that it resembles the simulation with only micrometer-sized particles. It should be pointed out that the Bell & Lin (1994) prescription includes the gas opacities as well, but these are relevant for high temperatures that are not reached within the extent of the discs here. We notice in Fig. 6 that the gradient after the iceline of the simulation with the Bell & Lin (1994) prescription or only micrometer-sized particles is much steeper than the corresponding gradient in the simulations of the full distributions. This is an important difference as the gradients in the aspect ratio affect the migration speed of planets that could be embedded in such a protoplanetary disc (see Sect. 5.2).

In conclusion, we find that including either the BOD or the MRN distribution leads to comparable results. The differences between the two grain size distributions tested in this work for different values of the turbulence parameter α and differentsurface densities are discussed in Sect. 4. Prior to that, a more extended discussion on the dust surface densities, dominant grain sizes, opacities, and temperatures follows in the next paragraphs.

Fig. 5 Aspect ratio (left plot) and midplane temperature (right plot) as a function of orbital distance in AU, for discs with five different single grain sizes from 0.1 μm to 1 mm (see Fig. 2 for the opacities of those five grain sizes). All of the simulations include viscous heating and stellar irradiation, have α = 5 × 10−3 as the turbulence parameter in viscosity, Σg,0 = 1000 g cm−2 as the initial gas surface density, and the dust-to-gas ratio is fDG = 1%. In this set of simulations, we also consider settling so that we can compare with the simulations that include the grain size distributions. The gray areas in the plot indicate the water iceline transition. Overplotted with dashed lines are the discs with the MRN distribution in reddish pink and the BOD distribution in dark blue. The simulations with the distributions show influence from small particles at the inner part of the disc, while the outer parts are mostly affected by larger particles, around 0.1 mm (see Fig. 8). The small differences in the aspect ratio and temperature of the discs with the distributions stem from the different slope of the vertically integrated dust surface densities of the two distributions (see Figs. 3 and A.1).

Fig. 6 Aspect ratio (left) and midplane temperature (right) as a function of orbital distance for the discs with the BOD, the MRN distribution and a disc that utilizes the Bell & Lin (1994) opacities. All of the simulations have α = 5 × 10−3, Σg,0 = 1000 g cm−2, and fDG = 1%. We also show the aspect ratio of the simulations with 1 and 100 μm for reference. The Bell & Lin (1994) opacities are based on micrometer-sized particles resulting in comparable aspect ratios. The main differences with the discs including the full grain size distributions are the steeper radial temperature (and thus aspect ratio) gradient for the disc with the Bell & Lin (1994) opacities and the reversed slope within 4 AU. These influence the migration speed and direction of planets embedded in those protoplanetary discs (see Figs. 13–15). The gray areas correspond to the water iceline transition.

3.2 Opacities and temperature

The opacity and temperature within the disc for the BOD and MRN distribution are plotted in Fig. 7 as a function of orbital distance on the x-axis and height on the y-axis. The total opacity of the disc is determined by accounting for the contribution of each grain size according to Eq. (22).

The highest opacity values in the figures correspond to the iceline as it can be also seen in the temperature plot (gray band). Almost every particle size has its highest opacity at the iceline, consequently, the total opacity of the disc is the highest at the iceline, as it can be seen already in Fig. 2. It was already briefly mentioned in Sect. 3.1 (detailed discussion in Sect. 3.4) that the dominant grain sizes at the inner disc are small, therefore we see the same pattern in the opacity of the disc around the iceline as the opacities of the small particles, with the bump around the transition. For this reason, we are tracing the iceline at the opacity plot.

The total opacity is scaled down in the simulation with the MRN distribution compared to the one with the BOD, which is what we also find in the aspect ratio. Since the opacities have the same pattern and since the bumps in the aspect ratio are caused by opacity transitions it is expected to find there the same gradient.

Viscous dissipation is the dominant source of heating for the inner parts of the disc, while stellar irradiation becomes important at larger orbital distances and more importantly for the upper layers of the disc (Dullemond et al. 2001; Bitsch et al. 2013a). Since the upper layers are heated up, the opacities are also higher there. If we move vertically up, the iceline moves inwards as viscous heating becomes weaker. The specific radius and height at which stellar heating begins to affect the structure of the disc are, among other parameters, influenced by the strength of turbulence. The dependence of the disc’s thermal structure on the turbulence strength is discussed in Sect. 4 and more opacity plots as a function of orbital distance and height can be found in Appendix B.

In Fig. 7 the τ = 1 line is overplotted (light blue line). The optical depth τ is defined as $$\tau ={\displaystyle {\int }_{z_{\text{max}}}^0{\kappa }_{\text{R}}}{\rho }_{\text{g}}\text{d}z,$$(25)

therefore it increases as the height z is decreasing. The τ = 1 line marks the difference between the optically thick and the optically thin medium. When τ ≥ 1, then the disc is optically thick, which means that the mean free path of the photons is much smaller than the length scale over which temperature changes. In optically thin parts of the disc, photons can “freely” travel out of the disc. The τ = 1 line thus marks the region of the disc where cooling becomes efficient. A τ = 1 line close to midplane corresponds to lower opacities, which results in a cooler disc. Even though the regions above this line are optically thin and cool down very efficiently, the uppermost layers are directly heated by stellar irradiation and we also see an increase in the opacity.

The transition from the optically thin part of the disc to the optically thick (moving from the top layers towards the midplane) is also where the boundary for observations would be if these were integrated over all wavelengths. Mid-infrared wavelength observations of the optically thick disc probe the temperature of the dust “photosphere”, the effective surface layer of the disc. Observations are in general carried out at various wavelengths, thus probing different grain sizes and different information for the disc (Andrews 2015). The optical depth relevant for such observations might differ for individual grain sizes.

3.3 Dust surface densities

The vertically integrated dust surface densities per orbital distance and grain size are presented in Appendix A. The maximum grain sizes in both of the distributions depend on the local temperature and gas surface density, which change with time. Additionally, all of the boundary sizes of the BOD distribution depend on the local gas surface density.

We stress here the loop that is created; the dust surface densities play a major role in determining the opacity of the disc (Eq. (22)), which then influences the cooling rate and the stellar heating and thus changes the temperature. The shift in the temperature leads to a new fragmentation barrier (and regime boundaries for the BOD), hence the dust surface densities change and so forth. Given the fact that this loop exists between the dust and the gas, it is important to consider the self-consistent calculations of the dust surface densities in the simulations.

Fig. 7 Mean Rosseland opacity values for the disc that includes the BOD distribution (top) and the MRN (bottom). These opacities correspond to the discs with the grain size distributions which were presented in Figs. 5 and 6. The highest opacity values are found at the iceline transition (gray band in the right plots). Moving outwards, the disc gets colder and the opacity of the disc decreases. The light blue line is the location where the vertically integrated optical depth reaches unity (τ = 1), so it divides the disc in the optically thick lower region and the optically thin upper region. Above this line, the opacities decrease due to the efficient cooling of the disc. The uppermost layers show increased opacities caused by the direct stellar heating, which increases the temperature. The same gradients can be seen in the temperature plots (right).

3.4 Dominant grain sizes

Depending on the local gas disc parameters, the grain size which plays the role of the dominant opacity source will change. We find the individual contribution of each grain size to the total opacity of the disc through Eq.(22). For each set of particles of the same size we calculate its contribution, $\frac{{\rho }_{\text{d},i}}{{\rho }_{\text{g}}}100{\kappa }_i$ and we present it as a function of orbital distance in Fig. 8 for the nominal case of α =5 × 10⁻³,  Σ_g,0 = 1000 g cm⁻². In order for a grain size to dominate the opacity, it needs a combination of high dust-to-gas volume density ratio (for this specific particle size) and high opacity at the given part of the disc.

We can see in the same figure, plotted as a black solid line, that the inner disc with either the BOD or the MRN distribution is influenced by very small particles, around 3 μm. Specifically, those small particles have the maximum contribution to the total opacity; their opacity dominates the opacity of the disc. The dominant grain sizes in the disc with the BOD grain size distribution feature a jump at around 20 AU. This jump is caused by the dip in the distribution in the transition between the first turbulence regime and the second (see Fig. 3), after which particles arelarge enough to get completely decoupled from the gas. After this jump, the dominant grain size is near 200 μm. The dominant grain size in the disc with the MRN distribution smoothly increases in the inner regions and then remains constant at around 90 μm exterior to 20 AU.

At the inner, hot parts of the disc, the grains around the micrometer size have surface densities that are around an order of magnitude lower than the larger grains (see Appendix A). However, they have the highest opacity by several orders of magnitude (see Fig. 2) at these high temperatures. This results in them being the dominant particles in that region. Farther out, the temperature of the disc decreases and as a consequence, larger particles carry the highest opacity. The decreased temperature means that the opacities of the larger particles get comparable with those of the smaller particles. At the same time, the difference between the very small and the largest grain sizes slightly increases, aiding the dominance of the largest grains.

In Fig. 8 the grain sizes below the dark red line give 25% of the contribution to the total opacity, thus the grain sizes above this line give 75% contribution. The line above this divides the grain sizes into two groups, each contributing50% to the opacity of the disc. In the same way, the uppermost line has grain sizes with 75% of the contributionbeneath it and sizes with 25% contribution above it. These lines show the general trend of the contribution to the opacity, which mainly comes from the small grains in the inner disc and by the large grains in the outer disc.

For a given grain size, the contribution to the opacity shows similar patterns as the dust surface densities. The maximum opacities per grain size are seen at approximately 6.5 AU. This location corresponds to the iceline and it is expected to have the highest opacity contribution by almost all of the grain sizes (Fig. 2). In conclusion, what defines the grain size with the maximum contribution to the total opacity is the combination of the dust surface densities and the opacity of each grain size at a specific orbital distance (which is determined by the local temperature). Once more it is evident that the self-consistent calculations within the feedback loop (Fig. 1) are crucial to the disc structure evolution.

Fig. 8 Contribution to the midplane mean Rosseland opacity per grain size for the simulation with the BODon the top and the MRN distribution on the bottom, for the nominal disc parameters α = 5 × 10−3, Σg,0 = 1000 g cm−2 and fDG = 1%. The black line indicates the grain sizes that contribute the most to the opacity of the disc as a function of orbital distance. The BOD distribution causes a jump in these dominant grain sizes because of the dip in the dust surface density (see Fig. 3) in the transition between the two turbulence regimes. The total opacity at each radius is the sum of the contributionfrom each grain size, or in other words, it is the sum of the corresponding column. For each grain size, the maximum opacities can be found at around 6–6.5 AU, where the iceline is located. The red lines show the percentage of the contribution from the grains below the corresponding line. This clearly illustrates that the dominant grain size does not necessarily determine the whole opacity. At the same time, we can see that the smallest grain sizes (up to roughly 50 μm beyond 20 AU) contribute the least to the total opacity.

4 Dependence on viscosity, gas surface density and dust-to-gas ratio

Decreasing the α-viscosity values is expected to affect the discs in two ways. First of all, the viscous heating decreases, therefore the discs will cool down and their aspect ratio will be lower (see Fig. 9). Secondly, the lower turbulence means that particles will face less destructive collisions and thus they will be able to grow to larger sizes. The larger grains have in general lower opacities, which means that the discs experience an additional cooling because of the change in the opacities. The general trend that we show in Fig. 9 is that the aspect ratio indeed decreases as the turbulence parameter decreases. The location of the iceline also moves further in.

The aspect ratios and corresponding midplane temperatures in the low α models (bottom plots in Fig. 9) show some wiggles due to convection. Convection is caused by the vertical temperature gradient which depends on the opacity gradient in the vertical direction and is present in the optically thick regions of the disc (Bitsch et al. 2013b). This also implies that as the grains are more vertically diffused in the higher α case and the vertical temperature gradients are less steep, the effect should be less strong. Indeed, convection is also present in the high α simulations, but only at the inner, hotter regions of the disc. All regions that are affected by convection experience some sort of instability, so that reaching a steady state is very hard, if not impossible.

The vertically integrated dust surface densities as a function of orbital distance and grain size for the simulations with the two grain size distributions and the rest of the α values (10⁻², 10⁻³, 5 × 10⁻⁴, 10⁻⁴) are presented in Fig. A.1 in the appendix. In Sect. 3 we presented the vertically integrated dust surface densities for the nominal value of α = 5 × 10⁻³.

We also plotted the 50% contribution line for the same α value in Fig. 8. This line divides the grain sizes into two groups that contribute equally to the total opacity at the given orbital distance. In Fig. 10 we present the comparison of the 50% contribution lines as a function of orbital distance and α-viscosity. We plot also, in thicker lines of the same color, the maximum grain size in each disc so that the two groups of grain sizes with equal contribution to the opacity of the disc can be seen. The lower boundary of this plot corresponds to the minimum grain size (constant in all of the simulations).

We find that the position of this 50% contribution line is similar almost independently of the grain size distribution utilized in the model. Similarly the position of the 50% contribution line within the grain size range (vertical axis in Figs. 8 and 10) does not change significantly as α decreases (see Fig. 10), but at the same time the maximum grain size increases. The very large particles (≥ 100 μm) have significantly lower opacities and thus each order of magnitude added in grain size only adds a small contribution to the total opacity of the disc.

If we decrease the gas surface density then the total dust surface density also scales down. The reduced surface density results in a colder disc because of two effects. On one hand, the viscous heating decreases and on the other hand the radiative cooling increases, as it is inversely proportional to the disc’s density. This is what we find in the simulations with the lowest initial gas surface density we used, namely Σ_g,0 = 100 g cm⁻² (Fig. 11). The discs with lower surface density are much colder, therefore have a significantly lower aspect ratio. In this case, the difference in the aspect ratio of the discs with the two distributions almost vanishes completely, therefore the position of the iceline is also practically the same for the two discs. This is explained by the fact that the dust surface densities of the two distributions are comparable, leading to contributions to the opacity from similar grain sizes (see Fig. A.1).

On the contrary, if we increase the dust-to-gas ratio the total dust surface density is by definition enhanced. This means that the viscous heating is higher and the cooling rate decreased, resulting in hotter discs with higher aspect ratios. In addition to that, the opacity increases because of the enhanced total dust surface densities and as a consequence the optically thick region of the discs extends to higher heights compared to the discs with lower dust-to-gas ratio (see an example of a disc with f_DG = 3% in Appendix A).

The opacities as a function of orbital distance and height for the discs with the two distributions and the rest of the α values (10⁻², 10⁻³, 5 × 10⁻⁴, 10⁻⁴) can be found in Fig. A.2 in the appendix. The turbulence strength affects the viscous heating and the orbital distance and height at which stellar heating takes over. This direct influence on the thermal structure of the disc leads to different opacities at each position in the disc, depending on the α-viscosity value. A decrease in the gas surface density, as mentioned, decreases the viscous heating and increases the cooling rate, but the lower temperaturedecreases the opacity of the disc, and cooling is enhanced. A comparison between the opacities in discs with all of the α values can be found in Fig. A.2. We also show in Appendix A the opacity of the disc with the lowest initial gas surface density Σ_g,0 = 100 g cm⁻² (with α = 10⁻², f_DG = 1%) and the opacity of the disc with the highest dust-to-gas ratio, f_DG = 3% (with α = 10⁻² again and Σ_g,0 = 1000 g cm⁻²).

Fig. 9 Aspect ratio as a function of orbital distance for the discs with high α-viscosity values (top left), namely α = 10−2, 5 × 10−3 and low α-viscosity values (bottom left), namely α = 10−3,  5 × 10−4 and 5 ×10−4, for the two grain size distributions (BOD in dark blue, MRN in reddish pink). The iceline moves inwards as α decreases, due to reduced viscous heating. The wiggles that can be seen in parts of the low α discs are caused by convection. Right plots: temperature as a function of orbital distance for the discs with high α-viscosity values (top), namely α = 10−2, 5 × 10−3 and for low α-viscosity values (bottom), namely α = 10−3,  5 × 10−4 and 10−4.

Fig. 10 50% contribution lines as a function of orbital distance for the discs with the high α values at the top and the low α values at the bottom. The group of thicker lines at the top corresponds to the maximum grain sizes of each disc so that the 50% contribution lines below divide the grain sizes into two groups, which contribute equally to the total opacity of the disc (see the red lines in Fig. 8). The difference between the two grain size distributions is small. As α-viscosity decreases, the maximum grain size increases, and more influence to the total opacity comes from larger grains.

5 Implications

5.1 Iceline

The location of the iceline can be theoretically calculated by considering the viscous and stellar radiation heating and the partial pressure of water vapor (e.g., Podolak & Zucker 2004; Davis 2005; Ciesla & Cuzzi 2006). But it is also well predicted by setting a single sublimation temperature (e.g. Hayashi 1981; Min et al. 2011) as we do here as well. Nevertheless, the location where this sublimation temperature is reached depends on several parameters, as described by the feedback loop (Fig. 1).

Planet formation studies indicate that the iceline in protoplanetary disc models should be outside of 1 AU, the Earth’s orbit. Otherwise, mechanisms like blocking the inward flow of pebbles from the outer disc by growing planets need to be invoked to keep the planets in the inner system dry (Morbidelli et al. 2016). If we take into consideration the observations of the composition of the bodies in the asteroid belt, the iceline at the time of formation of the asteroids would be at ~2.7 AU. But if icy planetesimals were formed at such small orbital distances and contributed to the formation of the terrestrial planets, we would observe larger amounts of water on Earth than what we observe today. The composition of asteroids from the inner region of the asteroid belt suggests that at their time of formation they should have been interior to the iceline.

Evolving disc models indicate that at the time of the formation of terrestrial planets the iceline has already moved towards 1 AU (Bitsch et al. 2015a). In Oka et al. (2011) it is suggested that in order to reach a better conclusion about the location of the iceline a grain size distribution is required, rather than uniform dust grain sizes. In Garaud & Lin (2007) the decoupling of dust particles from gas is discussed as a potential influence on the thermal structure of the disc. Lecar et al. (2006) argue that in order to move the iceline outside 3 AU one possible solution would be to increase the opacity used in their model. Here we do not include the opacity of a single grain size, but use an evolving grain size distribution that regulates the opacity of the disc more realistically because we take into account their individual contributions.

All of these suggested effects have been therefore taken into account in the here presented work where we study the influence of α-viscosity, initial gas surface density, and total dust-to-gas ratio on the position of the iceline. In Appendix B we present the simulations which were used and the procedure that was followed to do the fitting of the iceline position as a function of those three disc parameters. We find that the location of the iceline is independent of the grain size distribution which was utilized in the disc and it follows

Fig. 11 Aspect ratio as a function of orbital distance for two different surface densities (Σg,0 = 100 g cm−2 and Σg,0 = 1000 g cm−2). The aspect ratio decreases when the gas surface density decreases (while keeping the α-viscosity constant at the nominal value of 5 × 10−3) due to reduced viscous heating and increased cooling, which is inversely proportional to the density of the disc. As a result, the position of the water iceline moves inwards, close to 1 AU. This inward movement of the water iceline due to a reduced gas surface density follows the trend of previous disc evolution simulations that show an inward movement of the water icelineas the disc evolves and the gas surface density reduces (e.g., Oka et al. 2011; Bitsch et al. 2015a; Baillié et al. 2015). The low gas surface density results in almost the same aspect ratio profile for the discs with the two distributions.

$$r_{\text{ice}}=9.2\cdot {\left(\frac{\alpha }{0.01}\right)}^{0.61}\cdot {\left(\frac{{\Sigma }_{\text{g},\text{0}}}{1000 \text{g}{\text{cm}}^{-2}}\right)}^{0.8}\cdot {\left(\frac{f_{\text{DG}}}{0.01}\right)}^{0.37} \text{AU}\text{.}$$(26)

In order to investigate the theoretical background of the above power law fit we can start, as in Bitsch et al. (2013a), by considering the heating and cooling balance, Q⁺ = Q⁻, which means $$2\sigma T_{\text{eff}}^4=\frac{9}{4}{\Sigma }_{\text{g}}\nu {\Omega }_{\text{K}}^2.$$(27)

Given that the midplane temperature can be expressed as $T_{\text{mid}}={\left(\frac{3{\tau }_{\text{d}}}{4}\right)}^{1/4}{T}_{\text{eff}}$, the above equation becomes $$\frac{8\sigma }{3{\tau }_{\text{d}}}{T}_{\text{mid}}^4=\frac{9}{4}{\Sigma }_{\text{g}}\nu {\Omega }_{\text{K}}^2.$$(28)

We also substitute viscosity with the α prescription (Eq. (23) and Eq. (16) for the sound speed) and express the vertical optical depth as ${\tau }_d=\frac{1}{2}{\Sigma }_{\text{g}}{f}_{\text{DG}}\kappa$, so we get $$T_{\text{mid}}^3=\left(\frac{27}{64}\frac{k_{\text{B}}}{\sigma \mu m_{\text{H}}}\right){\Sigma }_{\text{g}}^2 f_{\text{DG}} \kappa \alpha {\Omega }_{\text{K}}.$$(29)

The surface density profile follows ${\Sigma }_{\text{g}}={\Sigma }_{\text{g},\text{0}}{\left(\frac{r}{\text{AU}}\right)}^{-1/2}$ and ${\Omega }_{\text{K}}=\sqrt{\frac{G{M}_{*}}{r^3}}$, so we obtain $$T_{\text{mid}}^3\propto {\Sigma }_{\text{g},\text{0}}^2 f_{\text{DG}} \kappa \alpha r^{-5/2}.$$(30)

We can then solve for the position of the iceline r = r_ice where T_mid = T_ice, $$r_{\text{ice}}\propto {\Sigma }_{\text{g},\text{0}}^{4/5} f_{\text{DG}}^{2/5} {\kappa }^{2/5} {\alpha }^{2/5}.$$(31)

We thus find that the power-law indices for the dependencies on Σ_g,0 and f_DG are very similar to what we find in our fitting (Eq. (26)). Comparing Eqs. (26) and (31) suggests that at the iceline κ ∝ α^1∕2, but is almost independent of Σ_g,0 and f_DG. The reason for this dependency has no easy analytical explanation, but it appears to be the outcome of the feedback between the disc structure and the dust evolution. This further illustrates, as we also discuss in Sect. 3.4, that we cannot rely on single grain size opacities to accurately describe the disc structure.

The position of the iceline as a function of the α-turbulence parameter and the initial gas surface density Σ_g,0 from our fitting formula is presented in Fig. 12, for discs with constant f_DG = 1%. The iceline transition is defined as the location where T = (170 ± 10) K. Increasing values of either one of the three parameters, α, Σ_g,0, f_DG, leads to hotter discs so that the iceline moves closer to the star (see Sect. 4). In the models with any of the grain size distributions, the iceline is located outside 1 AU for α ≥2.6 × 10⁻⁴ and exterior to 2.7 AU only when α ≥ 1.4 × 10⁻³ for a disc with Σ_g,0 = 1000 g cm⁻² and f_DG = 1%. However, this also depends on the disc’s surface density and total dust-to-gas ratio. If the surface density reduces in time, the disc becomes cooler and the iceline moves inwards, even for the high viscosity cases. For example, for Σ_g,0 = 100 g cm⁻² the iceline is located outside 1 AU for α ≥ 5.4 × 10⁻³ and exterior to 2.7 AU only when α ≥ 2.8 × 10⁻² (again with a constant f_DG = 1%). In contrast, if the total dust-to-gas ratio is increased to 3%, then we find the iceline outside 1 AU even for α ≥ 1.4 × 10⁻⁴ and outside 2.7 AU for α ≥ 6.9 × 10⁻⁴.

We can conclude that utilizing Mie theory for the opacities of the grains and taking into account a distribution of grain sizes helps in keeping the iceline sufficiently far out from the star, especially for high α and Σ_g,0 values. In general, the location of the iceline might also depend on the composition of the grains, which will be examined in future work. In contrast to Oka et al. (2011), who suggest that grain size distributions might help to keep the iceline at larger distances compared to single grain size discs, we actually find the opposite. Including a distribution in a disc simulation results in a similar position of the iceline to the discs with the smallest grain sizes (0.1 and 1 μm). At low viscosities, the opacity is dominated by larger grains and thus the disc becomes colder. Unrealistic single grain opacities (typically of micrometer size particles) result in discs that are too hot. This implies that potentially other heating sources are needed to keep the iceline at large orbital distances, especially if viscous heating is low (Mori et al. 2019).

Fig. 12 Iceline position as a function of α turbulence and initial gas surface density at 1 AU with a constant fDG = 1% (Eq. (26)). The iceline transition is defined as T = 170 ± 10 K. The black lines mark rice = 0.5, 1, 2.7 and 5 AU. Higher viscosity or gas surface density leads to hotter discs, with the iceline located at greater distances from the star. The same applies to higher total dust-to-gas ratio. The gray dashed lines mark rice = 0.5, 1, 2.7 and 5 AU for a disc with fDG = 3%.

5.2 Planet migration

The protoplanetary disc’s structure also affects planet migration. Very roughly said, if the aspect ratio increases with orbital distance then planets migrate inwards (Type-I migration, Paardekooper et al. 2011). On the contrary, if the aspect ratio is a decreasing function of orbital distance, planets will migrate outwards (Bitsch et al. 2013a, 2014, 2015a), if the viscosity is large enough (Baruteau & Masset 2008).

We focus here on the results with α = 5 × 10⁻³, a viscosity large enough to trigger outward migration by the entropy driven corotation torque. For the discs with the grain size distributions, we see an aspect ratio which is a decreasing function of orbital distance beyond the iceline, therefore in those disc regions planets could migrate outwards. Interior to the iceline the aspect ratio is an increasing function of orbital distance which means that planets embedded in this region of the disc would only migrate inwards. The minima in the aspect ratio are locations where planets could get trapped and if (as it is more likely) more than one planet existed, they could get into resonance and remain at those fixed orbital distances until the local parameters of the discs changed sufficiently to force them to migrate again (Horn et al. 2012; Cossou et al. 2013, 2014; Izidoro et al. 2017, 2019).

Simulations with single particle sizes larger than 100 μm show a monotonically increasing aspect ratio with orbital distance (Fig. 5), implying that planets in these colder discs would migrate inwards. The speed of the migration scales with the inverse of the square of the aspect ratio. If the aspect ratio decreases, then the migration is faster (for a review see Baruteau et al. 2014). Therefore, generally speaking, migration would be faster in the single grain size models with large particles and the exact migration speed of planets depends on their exact position in the protoplanetary disc.

To predict more precisely the migration of planets in the discs presented here, we include the migration maps for two discs with the distributions (Figs. 13 and 14). Migration rates are derived from the torque formula of Paardekooper et al. (2011). For comparison, we also show the migration map of the disc with the Bell & Lin (1994) opacities (Fig. 15), as these were the opacities used in Bitsch et al. (2015a) in Fig. 18. The black solid line in Figs. 13–15 encircles the regions of outward migration. Planets with masses less than approximately 10 M_⊕ always migrate inwards in the simulations with the distributions, whereas for the simulation using the Bell & Lin (1994), pure micrometer opacities, inward migration is the only possibility for planets less than 6 M_⊕.

The innermost region of outward migration in Fig. 15 corresponds to the area where the aspect ratio decreases as the orbital distance increases (Fig. 6). We can also see that this area has a steeper temperature gradient (Fig. 6) both interior and exterior to the water iceline transition. The increased torques and consequently migration speeds at the outer region of outward migration between 5 and 20 AU in the Bell & Lin (1994) disc simulations are caused by the steep increase of opacity for temperatures larger than 170 K, which is not the case for the BOD and MRN distribution, as these distributions are not dominated by micrometer grains, in contrast to the Bell & Lin (1994) opacity model. This illustrates that the grain sizes which dominate or contribute the most to the opacity of the disc have significant implications, not only for the disc structure itself but also indirectly for the planets embedded in that disc. In addition, this has important effects on the formation of planetary systems, because the migration rates determine how close planets can migrate towards each other, which sets the stability of the planetary system (Matsumoto et al. 2012; Chambers 2006).

Fig. 13 Torque acting on planets with different masses for the disc utilizing the BOD distribution for the nominal viscosity of α = 5 × 10−3. The black line encircles the regions of outward migration and corresponds to the region of the disc where the aspect ratio decreases as a function of the orbital distance. The temperature in the same region shows the steepest gradient.

Fig. 14 Same as Fig. 13 for the disc with the MRN distribution. The difference to the BOD distribution is small regarding the size of the region of outward migration, however, the torque is weaker for the MRN distribution. This could lead to different migration and growth behavior of planets forming in the outer disc.

Fig. 15 Same as Fig. 13 for the disc with the Bell & Lin (1994) opacity profile. In contrast to the discs with the BOD or the MRN grain size distribution, the disc simulated with the Bell & Lin (1994) opacity law shows an inner region of outward migration, which is caused by another bump in the H/r profile in the inner disc (Fig. 5).

5.3 Implications for planet formation and protoplanetary disc simulations

Our simulations presented here are the first step towards more self-consistent protoplanetary disc structure and evolution simulationsas well as planet formation simulations. Planet formation in the pebble assisted core accretion scenario relies crucially on the pebble sizes and distributions (e.g., Ormel & Klahr 2010; Lambrechts & Johansen 2012; Johansen & Lambrechts 2017, for review), as well as on the disc structure to calculate the planet migration rates as the planets grow (Bitsch et al. 2015b). The here presented model opens the avenue to simulations with self-consistent disc structures and pebble sizes, which can then be accreted onto planets. This can increase the accuracy of future planet formation simulations by pebble accretion.

The here used FARGOCA code also allows for 3D hydrodynamical simulations with embedded planets. A combination with the presented model of thermal structures calculated from full grain size distributions allows a very detailed comparison with observations, which are more advanced that the most frequently used 2D isothermal simulations followed by 3D radiative transfer (e.g., Zhang et al. 2018). This could potentially change our interpretation of observed protoplanetary discs featuring substructures potentially caused by planets.

6 Summary

We perform 2D hydrodynamical simulations including the whole feedback loop shown in Fig. 1. Specifically, we include and test two full grain size distributions and mean opacities (calculated via Mie theory) and study their influence on the disc structure. The particles have a minimal size of 0.025 μm and the upper boundary is regulated by the fragmentation barrier (Eq. (15)). We test five different α-viscosity values (10⁻², 5 × 10⁻³, 10⁻³, 5 × 10⁻⁴, 10⁻⁴), five values of initial gas surface density Σ_g,0 (100, 250, 500, 750 and 1000 g cm⁻²) and five values of dust-to-gas ratio f_DG (1, 1.5, 2, 2.5, 3%). We also perform simulations with only single grain sizes and with the Bell & Lin (1994) opacity law for comparison and in order to understand to greater extent the influence of grains of different sizes on the thermal structures of protoplanetary discs.

The dust component in protoplanetary discs is believed to follow a size distribution, regulated by a coagulation-fragmen- tation equilibrium (Brauer et al. 2008; Birnstiel et al. 2011). We utilize and compare two different grain size distributions. The first and simple model (MRN; Dohnanyi 1969; Mathis et al. 1977; Tanaka et al. 1996) results from the equilibrium between fragmentation and coagulation, whereas the second and more complex model (Birnstiel et al. 2011, BOD) takes into account fragmentation, coagulation and also cratering and adjusts the dust surface densities according to the grain sizes and how they compare to the size of the gas molecules and the gas turbulent eddies.

The dust surface densities are calculated as dictated by the aforementioned grain size distributions and the dust grains are also vertically distributed according to their sizes taking into account the effect of settling. We also have a spatial distribution radially, since the size distribution depends on the local disc parameters and changes self-consistently. In conclusion, a whole loop of growth, fragmentation, and settling of the resulting grains for each vertical slice of the disk is modeled in our simulations and updated at every timestep according to the local disc parameters.

We show disc structures calculated with the full grain size distributions and single grain sizes in Figs. 5a and b. Additionally, we show that the grain sizes which dominate or contribute the most to the opacity of the disc are not the same at all orbital distances of the disc (Figs. 8 and 10). As a consequence, the opacity prescriptions which assume a single dust size lead to inaccurate calculations of the thermal structures of the discs.

It is also important to stress that the dust surface densities, or in other words the distribution of mass among the grain sizes, play a major role in determining the disc opacity (Eq. (22)), which in turn influences the cooling rate and the stellar heating and changes the temperature and surface density of the gas. This shift in the local disc parameters leads to a new fragmentation barrier (and regime boundaries for the BOD), therefore the dust surface densities change and so on. For this reason, it is important to include the self-consistent calculations of the dust surface densities in the simulations.

The two grain size distributions show minimal differences in the dust surface densities (Fig. A.1). The reason for this is that both of the grain size distributions we have used in the discs feature the same fragmentation barrier. Therefore the grain size range in the discs with either one of the distributions is similar. Any difference between the discs with the BOD distribution and the discs with the MRN distribution comes mainly from the difference in the surface densities as a function of grain size (see Figs. 3 and A.1), which is usually smaller than an order of magnitude. The dominant grain sizes (Sect. 3.4) might not be the same, because of the small differences in the dust surfacedensities per grain size, but the total opacity of the disc is similar independently of the grain size distribution.

With this accurate prescription, we investigate the dependency of the iceline position on the α-viscosity, the initial gas surface density and the dust-to-gas ratio, where we see the effect of the feedback loop and find $r_{\text{ice}}\propto {\alpha }^{0.61}{\Sigma }_{\text{g},\text{0}}^{0.8}{f}_{\text{DG}}^{0.37}$ (Eq. (26), Fig. 12) independently of the grain size distribution utilized in the disc model. Specifically, for high gas surface density (Σ_g,0 = 1000 g cm⁻²) the position of the iceline is exterior to 1 AU for α ≥ 2.64 × 10⁻⁴ and exterior to 2.7 AU only when α ≥ 1.35 × 10⁻³. For higher values than the nominal f_DG = 1% we find that the iceline moves closer to the star as it is expected by the enhanced dust surface densities and the consequent hotter discs. However, for the nominal f_DG = 1%, lowering thegas surface density results in colder discs and the iceline is below 2.7 AU, even for the high viscosity models.

The changes in the aspect ratio gradient as a function of orbital distance affect the regions where outward migration is possible for planets that could be embedded in the disc (Figs. 13–15). Utilizing an α-viscosity of 5 × 10⁻³, Σ_g,0 = 1000 g cm⁻² and f_DG = 1% we find thatthe regions where outward migration could be possible in the discs with the two distributions are similar to the one present in a disc with the Bell & Lin (1994) opacities, that feature only micrometer-sized grains, at around 5–15 AU for planets with masses greater than 10 M_⊕. However, the region is more extended for the disc with the BOD distribution (up to 20 AU) and the disc with the Bell & Lin (1994) opacities has one more outward migration regions, near the inner boundary (2–3 AU), which is not present in the discs with the grain size distributions.

We can hence conclude that given the complexity and computational expense of the BOD distribution and the fact that it does not take into account radial drift or bouncing of the dust particles it is not necessary to prefer it over a simple MRN-like power-law distribution.

As the iceline can be the starting point for planetesimal formation (Guillot et al. 2014; Schoonenberg & Ormel 2017; Drążkowska & Alibert 2017) it is important to have as realistic models as possible, therefore include the feedback loop of grain growth and thermodynamics in hydrodynamical models (Fig. 1). Given also the fact that dust in protoplanetary discs follows a size distribution regulated by a coagulation-fragmentation equilibrium, the opacity prescription of a single grain size is not able to accurately calculate the thermal structures of discs.

The here presented model has some limitations that we wish to further investigate in future work. Both of the distributions tested here neglect radial drift (Brauer et al. 2008; Birnstiel et al. 2012) and bouncing (Zsom et al. 2010; Lorek et al. 2018) which can be detrimental for grain growth and in general affect the dust dynamics and subsequently the disc’s thermal structure. Also the onset of convection in some regions and for a subset of α and Σ_g,0 values might change the vertical distribution of the grains beyond settling and turbulent stirring by viscosity, as taken into account here. A very important future step is to model accretion discs instead of equilibrium discs and in this way, we will be able to also study different evolutionary steps of the protoplanetary disc. The particle composition and abundances are also determinants for dust dynamics and opacities, so it is important to consider a population that is as realistic as possible and use more accurate fragmentation velocities depending on our dust composition. Similarly, we are assuming a dust population consisting of 50% silicates and 50% water-ice, but we can relax this assumption and test different fractions (as done for example, in Bitsch & Johansen 2016).

It is, therefore, evident that the prescription that we used and present for this work opens up new avenues for protoplanetary disc simulations and planet formation. The inclusion of the feedback loop of grain growth and disc thermodynamics leads to more self consistent simulations of protoplanetary accretion discs and planet formation simulations in the pebble accretion scenario. Eventually, such models target a more precise comparison of protoplanetary disc observations to simulations that allow us to move away from simple 2D isothermal models with post-processing of radiation transfer.

Acknowledgements

B.B. and S.S. thank the European Research Council (ERC Starting Grant 757448-PAMDORA) for their financial support. M.L. thanks the Knut and Alice Wallenberg Foundation (grant 2017.0287, PI: A. Johansen). S.S is a Fellow of the International Max Planck Research School for Astronomy and Cosmic Physics at the University of Heidelberg (IMPRS-HD). We would like to thank the anonymous referee for the valuable comments and suggestions that helped us improve the manuscript.

Appendix A Dust surface densities for different α-viscosity values

We present here the vertically integrated dust surface densities and the opacities for the simulations using the rest of the α values and the two grain size distributions. The maximum grain sizes of the BOD and MRN discs are approximately the same because they follow the same fragmentation barrier formula (Eq. (15)). The vertically integrated dust surface densities arearound one order of magnitude lower in the discs with the BOD grain size distribution compared to the discs with the MRN distribution for small particles (Fig. A.1). On the other hand, they are always comparable for the largest particles in the discs. This is already evident by the shape of the two grain size distributions (Fig. 3).

As α decreases, the surface densities of the smallest particles in the discs with the BOD distribution are several orders of magnitude lower than these of the largest particles. The gradients are smoother in the discs with the MRN distribution as α decreases. The fragmentation barrier depends on the initial gas surface density so when the latter decreases, the maximum grain size gets smaller (Eq. (15)).

In Fig. A.2 we show the opacities as a function of orbital distance and height for a selection of the simulated discs for this work. The τ = 1 line is located at similar heights in all of the discs with high α-viscosity (around 2 AU at the outer edge). The same applies to the low α-viscosity discs where the τ = 1 line is always around 0.15 AU at the outer edge. As expected the optically thick region is extended towards higher altitudes with a higher total dust-to-gas ratio and contained near the midplane for low gas surface densities. Above the τ = 1 line opacity always decreases as cooling is more efficient. However, the uppermost layers show increased opacities because of the stellar irradiation that directly heats them up.

Fig. A.1 Dust surface densities as a function of orbital distance and grain size for the different α values used here, and for additional simulations with the lowest gas surface densities and with the highest dust-to-gas ratio that we have tried. When the turbulence is reduced, the maximum grain size increases (since less destructive collisions are expected). In addition to that, the reduced α-viscosity allows the discs to become cooler, so the opacity of the larger grains (~mm) becomes comparable or larger than that of the smaller (~μm) grains. For α = 10−4 we find thatthe grains grow up to a few centimeters in the discs with either one of the distributions. The dashed lines divide the grain sizes into two groups which contribute equally to the total opacity of the disc (see Figs. 8and 10).

Fig. A.1 continued.

Fig. A.2 Mean Rosseland opacities as a function of orbital distance and height for the different α values, the lowest gas surface densities, and the highest dust-to-gas ratio. Due to the larger grain sizes for the MRN distribution (Fig. A.1), the opacities for the MRN distribution are also generally lower compared to the BOD distribution. The light blue line corresponds to optical depth τ = 1 integrated vertically starting from infinity towards midplane, so it divides the optically thin (above) and thick region (below).

Fig. A.2 continued.

Appendix B Iceline position as a function of α-turbulence, initial gas surface density and dust-to-gas ratio

In Sect. 5.1 we present the position of the iceline as a function of the α-viscosity parameter, the initial gas surface density, Σ_g,0, and the total dust-to-gas ratio f_DG. In order to do this fitting, we used five simulations for each parameter (see Table B.1). In Fig. B.1 we show the individual fitting over each one of the parameters. The fit to the three parameters writes $$f=C\cdot {\left(\frac{\alpha }{0.01}\right)}^{p_1}\cdot {\left(\frac{{\Sigma }_{\text{g},\text{0}}}{1000 \text{g}{\text{cm}}^{-2}}\right)}^{p_2}{\left(\frac{f_{\text{DG}}}{0.01}\right)}^{p_3},$$(B.1)

with C = 9.20 ± 0.05 AU, p₁ = 0.61 ± 0.03, p₂ = 0.77 ± 0.03, p₃ = 0.37 ± 0.01. The resulting fit is the same regardless of the grain size distribution used in the disc (solid line in Fig. B.1).

Fig. B.1 Iceline position as a function of α-viscosity (top), initial gas surface density (middle), and dust-to-gas ratio (bottom) for the discs with the BOD distribution (left column) and the discs with the MRN distribution (right column). The iceline transition is defined as T = (170 ± 10) K. The specific parameters used for the simulations presented in this plot are shown in Table B.1. The solid lines are the fits to each parameter and are the same for the discs with either one of the distributions.

Appendix C The effect of settling

We implement in our work the effect of settling for the grains in the disc as described in Sect. 2.4, in order to vertically distribute the grains according to their sizes and the local disc parameters. This implies that the disc structure can be affected both by a change in the grain size due to the different opacities that each size provides (Fig. 1) and by a change in the settling efficiency of the given grain size. In order to test if both of these effects are significant factors that define the disc structure, we run one simulation where the disc only contains millimeter grains and compare with a disc which also contains only millimeter grains but does not take settling into account. Thus in this latter case, the millimeter grains are vertically distributed according to a constant dust-to-gas ratio throughout the whole disc. Additionally, we run a simulation where the opacities of the grains correspond to millimeter grains, but we assume that the grains are vertically distributed as micrometer grains (so we assume s = 1 μm in the equations describing settling, Eqs. (17)–(21)). We choose α = 10⁻⁴ for which the settling of millimeter grains will be very effective. However, micrometer dust grains are not expected to be affected by settling even at this low α-viscosity. The models also have an initial gas surface density of Σ_g,0 = 1000 g cm⁻² and total dust-to-gas ratio f_DG = 1%.

In Fig. C.1 we present the aspect ratio of these three disc models. The aspect ratio as a function of orbital distance for the disc where grains are vertically distributed as micrometer-sized dust resembles the one of a disc where the dust-to-gas ratio is constant all over the disc. This is expected because micrometer-sized dust is not significantly affected by settling even at the low α-viscosity of 10⁻⁴ (see Fig. 4). However, we find that the aspect ratio is lower in the inner regions of discs when the millimeter grains are allowed to settle with their corresponding properties. When settling is included, the millimeter grains are mainly concentrated near the midplane (see also Fig. 4), while at higher altitudes the opacity diminishes. Without settling or with reduced efficiency of settling the opacity is similar at all altitudes, which leads to a reduced cooling rate and higher aspect ratio in the inner regions.

Fig. C.1 Comparison of the aspect ratio as a function of orbital distance for discs with mm grains with settling, without settling and with the settling of μm grains. The disc structure is almost the same when no settling is implemented and when the opacities correspond to mm grains, but the grains are distributed as micrometer-sized dust (Eqs. (17)–(21)). When settling is included the millimeter grains are concentrated near the midplane, leaving the upper layers with diminished opacities and enhanced cooling rate. Without efficient settling, the inner discs get hotter. However, this creates a shadow that prevents theefficient heating of the outermost regions from stellar irradiation.

Without settling of the millimeter grains according to their size properties, the outer regions are not sufficiently heated. Due to the increased aspect ratio in the inner disc, stellar irradiation to the outer disc is diminished, creating a shadow that cools down the outer region. At the same time, the millimeter grains have very low opacity so they cannot absorb the stellar heating efficiently inthe outer disc. The disc at the outermost radii might keep cooling down until it reaches the temperature of the surroundings (Dullemond & Dominik 2004). Hence including settling is very important to avoid such complications and inconsistencies in the disc structures.

Different grain sizes lead to different disc structures even without any settling implemented. The distinctive structures of discs with different grain sizes come mainly by their individual opacities (Fig. 2). However, without settling the disc opacity above midplane is overestimated so the discs are hotter in the inner regions and thus do not allow the stellar irradiation to heat the outer regions. In order to consistently take into account the influence of a grain size distribution to the resulting disc structures, it is important to include settling.

References

All Tables

All Figures

Fig. 1 Graphical illustration of the feedback loop. The thermal and density structure of a protoplanetary disc is determined by this loop: temperature and gas surface density affect the relative velocities for grains of different sizes. Through the relative velocities we find the outcomes of collisions between grains, therefore a grain size distribution is created. The grains are then vertically distributed according to their sizes and the turbulence strength. The spatial distribution of the grain sizes determines the opacity of the disc, which then affects its cooling rate and the stellar heating. This way the temperatureand density of the disc change and subsequently its whole structure.

In the text

Fig. 2 Rosseland, Planck, and stellar mean opacities (from left to right) as a function of temperature for grains of sizes 0.1, 1, 10, 100 μm, 1 mm, and 1 cm. They are independent of the gas density because they are dominated by the dust component. These values were calculated using RADMC-3D for a mixture of 50% silicates, 50% ice, and disc dust-to-gas ratio of 1%. The gray vertical line shows the location of the water iceline transition (170 K ± 10 K), causing a transition in opacity due to the evaporation/condensation of dust grains.

In the text

Fig. 3 Vertically integrated dust surface density distribution per logarithmic bin of grain size as a function of grain size for the two distributions used here, after Birnstiel et al. (2011; BOD) and Mathis et al. (1977; MRN), at 10 AU, for the simulation with α = 5 × 10−3 (see Fig. A.1 for the distributions over all orbital distances). The dust to gas ratio is 1% and the gas surface density is 1000 g cm−2 at 1 AU in both cases. For both of the distributions we additionally used uf = 1 m s−1 and ρs = 1.6 g cm−3. In the BOD, we can distinguish three regions. Small particles follow Brownian motion (BM), then as they grow they followthe turbulent gas motions (T1), and finally they are affected by the turbulence, but get decoupled from the gas as their stopping times are much larger than the turn-over time of the eddies (T2). The barrier sBT is the grain size limit for Brownian motion, while the barrier s12 separates the two turbulent regimes. The bump near the end of the distribution is caused by cratering, since large particles only lose part of their mass this way, while small particles can only coagulate and form larger grains. This causes the distribution to be top-heavy. The MRN distribution is a single power-law (Eq. (13)). Both grain size distributions have a maximum value determined by the same fragmentation limit (Eq. (15)), but they are not the same for the BOD and the MRN distribution because of the self-consistently calculated temperature.

In the text

Fig. 4 Dust density as a function of height for grains of five representative grain sizes within a disc with the BOD, α = 10−4 and Σg,0 = 1000 g cm−2, at 3 AU. The black line shows the gas density of the disc. The dashed gray line shows the dust-to-gas ratio as a function of height for the whole grain size distribution. The dust-to-gas ratio of the smallest particles of the sample remains constant, but as grains get larger they get more affected by settling. Consequently, the largest grain (1 cm) is mostly concentrated at the midplane.

In the text

Fig. 5 Aspect ratio (left plot) and midplane temperature (right plot) as a function of orbital distance in AU, for discs with five different single grain sizes from 0.1 μm to 1 mm (see Fig. 2 for the opacities of those five grain sizes). All of the simulations include viscous heating and stellar irradiation, have α = 5 × 10−3 as the turbulence parameter in viscosity, Σg,0 = 1000 g cm−2 as the initial gas surface density, and the dust-to-gas ratio is fDG = 1%. In this set of simulations, we also consider settling so that we can compare with the simulations that include the grain size distributions. The gray areas in the plot indicate the water iceline transition. Overplotted with dashed lines are the discs with the MRN distribution in reddish pink and the BOD distribution in dark blue. The simulations with the distributions show influence from small particles at the inner part of the disc, while the outer parts are mostly affected by larger particles, around 0.1 mm (see Fig. 8). The small differences in the aspect ratio and temperature of the discs with the distributions stem from the different slope of the vertically integrated dust surface densities of the two distributions (see Figs. 3 and A.1).

In the text

Fig. 6 Aspect ratio (left) and midplane temperature (right) as a function of orbital distance for the discs with the BOD, the MRN distribution and a disc that utilizes the Bell & Lin (1994) opacities. All of the simulations have α = 5 × 10−3, Σg,0 = 1000 g cm−2, and fDG = 1%. We also show the aspect ratio of the simulations with 1 and 100 μm for reference. The Bell & Lin (1994) opacities are based on micrometer-sized particles resulting in comparable aspect ratios. The main differences with the discs including the full grain size distributions are the steeper radial temperature (and thus aspect ratio) gradient for the disc with the Bell & Lin (1994) opacities and the reversed slope within 4 AU. These influence the migration speed and direction of planets embedded in those protoplanetary discs (see Figs. 13–15). The gray areas correspond to the water iceline transition.

In the text

Fig. 7 Mean Rosseland opacity values for the disc that includes the BOD distribution (top) and the MRN (bottom). These opacities correspond to the discs with the grain size distributions which were presented in Figs. 5 and 6. The highest opacity values are found at the iceline transition (gray band in the right plots). Moving outwards, the disc gets colder and the opacity of the disc decreases. The light blue line is the location where the vertically integrated optical depth reaches unity (τ = 1), so it divides the disc in the optically thick lower region and the optically thin upper region. Above this line, the opacities decrease due to the efficient cooling of the disc. The uppermost layers show increased opacities caused by the direct stellar heating, which increases the temperature. The same gradients can be seen in the temperature plots (right).

In the text

Fig. 8 Contribution to the midplane mean Rosseland opacity per grain size for the simulation with the BODon the top and the MRN distribution on the bottom, for the nominal disc parameters α = 5 × 10−3, Σg,0 = 1000 g cm−2 and fDG = 1%. The black line indicates the grain sizes that contribute the most to the opacity of the disc as a function of orbital distance. The BOD distribution causes a jump in these dominant grain sizes because of the dip in the dust surface density (see Fig. 3) in the transition between the two turbulence regimes. The total opacity at each radius is the sum of the contributionfrom each grain size, or in other words, it is the sum of the corresponding column. For each grain size, the maximum opacities can be found at around 6–6.5 AU, where the iceline is located. The red lines show the percentage of the contribution from the grains below the corresponding line. This clearly illustrates that the dominant grain size does not necessarily determine the whole opacity. At the same time, we can see that the smallest grain sizes (up to roughly 50 μm beyond 20 AU) contribute the least to the total opacity.

In the text

Fig. 9 Aspect ratio as a function of orbital distance for the discs with high α-viscosity values (top left), namely α = 10−2, 5 × 10−3 and low α-viscosity values (bottom left), namely α = 10−3,  5 × 10−4 and 5 ×10−4, for the two grain size distributions (BOD in dark blue, MRN in reddish pink). The iceline moves inwards as α decreases, due to reduced viscous heating. The wiggles that can be seen in parts of the low α discs are caused by convection. Right plots: temperature as a function of orbital distance for the discs with high α-viscosity values (top), namely α = 10−2, 5 × 10−3 and for low α-viscosity values (bottom), namely α = 10−3,  5 × 10−4 and 10−4.

In the text

Fig. 10 50% contribution lines as a function of orbital distance for the discs with the high α values at the top and the low α values at the bottom. The group of thicker lines at the top corresponds to the maximum grain sizes of each disc so that the 50% contribution lines below divide the grain sizes into two groups, which contribute equally to the total opacity of the disc (see the red lines in Fig. 8). The difference between the two grain size distributions is small. As α-viscosity decreases, the maximum grain size increases, and more influence to the total opacity comes from larger grains.

In the text

Fig. 11 Aspect ratio as a function of orbital distance for two different surface densities (Σg,0 = 100 g cm−2 and Σg,0 = 1000 g cm−2). The aspect ratio decreases when the gas surface density decreases (while keeping the α-viscosity constant at the nominal value of 5 × 10−3) due to reduced viscous heating and increased cooling, which is inversely proportional to the density of the disc. As a result, the position of the water iceline moves inwards, close to 1 AU. This inward movement of the water iceline due to a reduced gas surface density follows the trend of previous disc evolution simulations that show an inward movement of the water icelineas the disc evolves and the gas surface density reduces (e.g., Oka et al. 2011; Bitsch et al. 2015a; Baillié et al. 2015). The low gas surface density results in almost the same aspect ratio profile for the discs with the two distributions.

In the text

Fig. 12 Iceline position as a function of α turbulence and initial gas surface density at 1 AU with a constant fDG = 1% (Eq. (26)). The iceline transition is defined as T = 170 ± 10 K. The black lines mark rice = 0.5, 1, 2.7 and 5 AU. Higher viscosity or gas surface density leads to hotter discs, with the iceline located at greater distances from the star. The same applies to higher total dust-to-gas ratio. The gray dashed lines mark rice = 0.5, 1, 2.7 and 5 AU for a disc with fDG = 3%.

In the text

	Fig. 13 Torque acting on planets with different masses for the disc utilizing the BOD distribution for the nominal viscosity of α = 5 × 10−3. The black line encircles the regions of outward migration and corresponds to the region of the disc where the aspect ratio decreases as a function of the orbital distance. The temperature in the same region shows the steepest gradient.
In the text

	Fig. 14 Same as Fig. 13 for the disc with the MRN distribution. The difference to the BOD distribution is small regarding the size of the region of outward migration, however, the torque is weaker for the MRN distribution. This could lead to different migration and growth behavior of planets forming in the outer disc.
In the text

	Fig. 15 Same as Fig. 13 for the disc with the Bell & Lin (1994) opacity profile. In contrast to the discs with the BOD or the MRN grain size distribution, the disc simulated with the Bell & Lin (1994) opacity law shows an inner region of outward migration, which is caused by another bump in the H/r profile in the inner disc (Fig. 5).
In the text

Fig. A.1 Dust surface densities as a function of orbital distance and grain size for the different α values used here, and for additional simulations with the lowest gas surface densities and with the highest dust-to-gas ratio that we have tried. When the turbulence is reduced, the maximum grain size increases (since less destructive collisions are expected). In addition to that, the reduced α-viscosity allows the discs to become cooler, so the opacity of the larger grains (~mm) becomes comparable or larger than that of the smaller (~μm) grains. For α = 10−4 we find thatthe grains grow up to a few centimeters in the discs with either one of the distributions. The dashed lines divide the grain sizes into two groups which contribute equally to the total opacity of the disc (see Figs. 8and 10).

In the text

	Fig. A.1 continued.
In the text

Fig. A.2 Mean Rosseland opacities as a function of orbital distance and height for the different α values, the lowest gas surface densities, and the highest dust-to-gas ratio. Due to the larger grain sizes for the MRN distribution (Fig. A.1), the opacities for the MRN distribution are also generally lower compared to the BOD distribution. The light blue line corresponds to optical depth τ = 1 integrated vertically starting from infinity towards midplane, so it divides the optically thin (above) and thick region (below).

In the text

	Fig. A.2 continued.
In the text

Fig. B.1 Iceline position as a function of α-viscosity (top), initial gas surface density (middle), and dust-to-gas ratio (bottom) for the discs with the BOD distribution (left column) and the discs with the MRN distribution (right column). The iceline transition is defined as T = (170 ± 10) K. The specific parameters used for the simulations presented in this plot are shown in Table B.1. The solid lines are the fits to each parameter and are the same for the discs with either one of the distributions.

In the text

Fig. C.1 Comparison of the aspect ratio as a function of orbital distance for discs with mm grains with settling, without settling and with the settling of μm grains. The disc structure is almost the same when no settling is implemented and when the opacities correspond to mm grains, but the grains are distributed as micrometer-sized dust (Eqs. (17)–(21)). When settling is included the millimeter grains are concentrated near the midplane, leaving the upper layers with diminished opacities and enhanced cooling rate. Without efficient settling, the inner discs get hotter. However, this creates a shadow that prevents theefficient heating of the outermost regions from stellar irradiation.

In the text