© A. M. Arabhavi et al. 2022

1 Introduction

Protoplanetary disks are made up of gas and dust with different spatial distribution and temperature structures. The midplane is generally cold and dense for both dust and gas components (e.g. Bergin et al. 2007; Öberg & Bergin 2021). In those cold regions, which are well shielded from both stellar and interstellar ultra-violet (UV) radiation fields, gas phase species can accrete onto dust grain surfaces to form ices. The local thermo-chemical conditions in these regions determine the volume and chemical composition of the forming ice. This ice accretion affects the grain size and composition which are important for opacity calculations (Ossenkopf & Henning 1994; Allamandola et al. 1999; Boogert & Ehrenfreund 2004; Öberg 2009).

Observational data of dust in the diffuse interstellar medium, molecular clouds, and young stellar objects show the presence of ice covering the dust grains (Merrill et al. 1976; Whittet et al. 1983, 1998; Chiar et al. 2000, 2002; Boogert et al. 2008; Goto et al. 2018). For protoplanetary disks, the resolution, spectral range, and sensitivity requirements limit the observational capabilities. Yet some observations of H₂O ice in absorption at 3 µm (for example Terada et al. 2007; Honda et al. 2009, 2016; Terada & Tokunaga 2017) and emission at 47 µm and 63 µm (McClure et al. 2012, 2015; Min et al. 2016a) are available. Thus, disk modeling becomes important for understanding the properties of icy dust in disks as well as to design observing strategies for observing ices.

Observations, such as Terada et al. (2007) and Terada & Tokunaga (2012), indicate that the observed profile of an ice feature (for example the 3 µm H₂O ice feature) is a combination of diverse factors. In addition to the contributions of the ice species themselves, which include presence of different ice species such as CO, CH₄, NH₃, CH₃OH (Brooke et al. 1999; Dartois & d’Hendecourt 2001; Dartois et al. 2002), optical behavior of large grains (Smith et al. 1989; Boogert et al. 2015) and thermal state of the ices (Boogert et al. 2015; Sandford & Allamandola 1988; Hudgins et al. 1993; Palumbo & Strazzulla 1993) contribute to the shape of the ice feature. Grain growth is also expected in cold and dense environments (e.g. Terada & Tokunaga 2012; Terada et al. 2012; McClure et al. 2012).

Generally, disk models (e.g. Nomura & Millar 2005; D’Alessio et al. 2006; Semenov & Wiebe 2011; Bruderer et al. 2012; Woitke et al. 2016) assume a grain size distribution with a minimum and maximum grain size, and a power law with index a_pow that is either obtained by fitting the observations or set to 3.5 (MRN-distribution, Mathis et al. 1977). These grains are considered to be composed of refractory materials such as silicates and conducting materials. Some models consider the presence of ices by using ice opacities (e.g. Nomura & Millar 2005; Pontoppidan et al. 2005; McClure et al. 2012; Kamp et al. 2018). However, some of these models (Nomura & Millar 2005; D’Alessio et al. 2006; McClure et al. 2012, 2015) consider segregated grain distribution, that is, having separate size distribution of ice grains apart from refractory grains. Moreover, for calculating the opacities, these models (Nomura & Millar 2005; D’Alessio et al. 2006; McClure et al. 2015) either assume that ice exists everywhere below sublimation temperature, neglecting photodesorption, or do not consistently solve for photodesorption and ice abundance. Further, Nomura & Millar (2005), D’Alessio et al. (2006), McClure et al. (2012, 2015), and Kamp et al. (2018) consider only H₂O ice and Pontoppidan et al. (2005) consider H₂O, CO, and CO₂ ice species for opacity calculations. Segregated grain distribution do not account for increase in grain size due to ice accretion. Considering only H₂O ice underestimates the effect of ice accretion on size distribution of grains. A more recent study of ice opacities in disk models is presented in Ballering et al. (2021), where six ice species (H₂O, CO, CO₂, CH₃OH, NH₃, and CH₄) are considered along with two populations of refractory grains. The disk region is spatially classified into 25 zones which can be approximated with the same opacities. Thus avoiding the computationally intensive, yet more accurate, point-by-point opacity calculation. The effects of an interstellar UV field is also neglected. The model produces substantial ice absorption features including the 3 µm water ice band in disks that are not edge on. A more consistent treatment of the effect of ice on dust composition and size, hence, will help improve the understanding of ice observations in disks.

Ossenkopf & Henning (1994) show, for a fixed volume of water ice, that the opacity changes when ice mantles are present on dust grains, which is not taken into consideration by many disk models. Kamp et al. (2018) assume that H₂O ice is present in region defined by the snowline (Min et al. 2016a) which takes into account the temperature, pressure, and UV radiation field. The ice is distributed as ice mantles on top of underlying bare grain size distribution. Ices are assumed to be a constant mass fraction of the total dust mass, thus the ice abundance is not informed by the local chemical state. This leads to having too much ice in the upper midplane or too little ice in the inner midplane. Pontoppidan et al. (2005) use a single refractory grain size population, upon which ice mantle composed of H₂O, CO, and CO₂ ices is present. However, the model considers that ices evaporate at 90K (20K for CO) and also the local ice abundances are obtained by ice mixtures of fixed proportions (H₂O:CO₂:CO=1:0.2:0.03 and inclusions of CO₂:CO=1:0.7). In this paper, we present a model in which, for the first time, the ice chemistry, opacities, and radiative transfer are treated self-consistently.

We present here a computationally efficient method of including position dependent opacities in disk models, to account for ice abundances based on local thermo-chemical state, and discuss the implications for observations. In Sect. 2 the description of the disk model, PRODIMO, is provided in which the position dependent ice opacities are implemented. In Sect. 3, the method of calculating position dependent opacities is explained along with the verification of the calculation. The effect of the position dependent opacity on a typical T Tauri disk is discussed in Sect. 4. We also address two main factors that affect the ice features in our disk models, namely the relation between ice thickness and grain size, and the chemistry. We discuss the results of our models in the context of ice observations in Sect. 5. The method and results are summarized and concluded in Sect. 6.

2 Disk model and optical data

PRODIMO¹ or PROtoplanetary DIsk MOdel (Woitke et al. 2009; Kamp et al. 2010; Thi et al. 2011; Rab et al. 2018) is used in this paper to implement the position-dependent dust and ice opacities. The existing model, as explained in Woitke et al. (2009, 2016), uses a fixed grain size distribution with no ices to calculate the opacities. The model starts by setting up the physical state of the disk such as gas and dust densities including dust settling of the bare grains. The code then calculates the dust opacities and solves the continuum radiative transfer, thereby determining the dust temperature structure. Subsequently, the code computes the gas phase and ice chemistry and solves the heating & cooling balance to find the gas temperature structure. We use steady state chemistry in all the models presented in this paper. One model also includes grain surface reactions. More details about the chemistry and the disk parameters are given in Appendix B.

The bare grain opacities are explained in Woitke et al. (2016), based on Min et al. (2016b). The opacity calculations assume that the grains are spherical and the materials are well mixed. The refractory dust grains are assumed to be composed of 60% amorphous silicate (Dorschner et al. (1995), Mg_0.7Fe_0.3SiO₃) and 15% amorphous carbon (Zubko et al. 1996), with a porosity of 25%.

The computation method presented in this paper requires multiple iterations to converge. The first iteration uses bare grain opacities. The resulting position-dependent ice abundance and position-dependent settled bare grain size distribution are used to recompute the dust opacities including the ices in the next iteration. Further, the code continues to iterate between radiative transfer, chemistry with ice formation, and opacity calculation, until convergence is achieved. The models presented in this paper converged within three such iterations.

The opacity calculations are based on Mie theory as explained in Woitke et al. (2016). The dust grains (including those coated with ice) are considered to be spherical. Further, Woitke et al. (2016) use Mie theory and distributed hollow spheres (Min et al. 2005) for calculating opacities. However, the grain size distribution at each grid point in the disk changes and the distributed hollow spheres are computationally expensive, in this paper distributed hollow spheres are not used. In this paper, the dependence of ice thickness on grain size is generalized and made computationally easy to apply to disk models. The porosity of the refractory grains and ice layer is fixed to a certain constant value, that is, the possible enhancement or diminution of porosity upon ice accretion, grain collision, shattering, and other thermal events is neglected.

3 Opacity: Ice composition, thickness, and refractory materials

Extinction opacity $\left({\kappa }_{\lambda }^{\text{ext}}\right)$ can be defined as the total absorption $\left({\kappa }_{\lambda }^{\text{abs}}\right)$ and scattering $\left({\kappa }_{\lambda }^{\text{sca}}\right)$ cross sections of the dust distribution, as shown in Eq. (1) (Krugel 2002): $$\begin{array}{l}{\kappa }_{\lambda }^{\text{ext}}={\kappa }_{\lambda }^{\text{abs}}+{\kappa }_{\lambda }^{\text{sca}}\\ Q_{\lambda }^{\text{ext}}=Q_{\lambda }^{\text{abs}}+Q_{\lambda }^{\text{sca}}\\ {\kappa }_{\lambda }^{\text{ext}}={\displaystyle {\int }_{a_{\min }}^{a_{\max }}\pi }{a}^2{Q}^{\text{ext}}(a,\lambda )f(a)\text{d}a\end{array}$$(1)

where, Q^ext(a, λ), Q^abs(a, λ), Q^sca(a, λ) are the extinction, absorption, and scattering efficiencies defined as the ratio of respective cross sections to the geometric cross section of the grain, a is the grain size, a_min and a_max are the minimum and maximum grain sizes, f(a) is the considered size distribution ($=\frac{\text{d}N}{\text{d}a}$, i.e., number of grains, dN, per size bin da).

Due to ice accumulating on the bare (refractory) grains, both extinction efficiency Q(a, λ) and grain size distribution f(a) change. These directly affect the opacity κ_λ. The effect of ice on these two factors are discussed in the following subsections.

3.1 Effect of ice on size distribution

The abundance of ice varies throughout the disk (see appendix). This means that the thickness of ice layer differs at each point in the disk. The thickness of ice on the grains at any point depends on the local volume of ice available and the total grain surface area available for ice accretion at that point in the disk. A simple relation can be written as shown in Eq. (2). $$\begin{array}{l}\hfill V_{\text{IG}}=V_{\text{I}}+V_{\text{BG}}\hfill \\ n_{\text{d}}{\displaystyle {\int }_{a_{\min }}^{a_{\max }}\frac{4}{3}\pi }{\left(a+\Delta a\right)}^{3}f(a)\text{d}a=V_{\text{I}}+n_{\text{d}}{\displaystyle {\int }_{a_{\min }}^{a_{\max }}\frac{4}{3}\pi }{a}^{3}f(a)\text{d}a\end{array}$$(2)

where, V is volume of material per unit volume of space with subscripts IG for icy grains, I for ice, and BG for bare grains, a is the bare grain size, Δa is the ice layer thickness, f(a) [cm⁻¹] is the normalized size distribution function of the grains, a_min and a_max are the minimum and maximum sizes of the bare grains, and n_d [cm⁻³] is the number density of grains. We consider the following two cases. (1) If, in Eq. (2), Δa is constant for all grain sizes, then the ice/refractory volume ratio (V_I/V_BG) is significantly larger for small grains than for large grains. (2) Ice thickness can also depend on the grain size. This can be due to coagulation or shattering of icy grains in the midplane or heating and re-condensation events (the latter leads to growth of ice faster on larger grains, Kuroiwa & Sirono 2011). In this paper, we assume that the limiting case of such scenario is an uniform ice/refractory volume ratio for grains of different sizes.

In order to allow for a general approach, which includes both limiting cases, we assume that the resulting grains are still spherical. For an icy grain of certain refractory size a, if V_I, V_BG, and V_IG are volumes of ice, refractory material, and total volume of the icy grain, respectively, then the ratio V_I/V_IG and V_BG/V_IG are the volume fractions of ice and refractory material in an icy grain respectively. If a grain with refractory size a is present in an environment where there is no coagulation, shattering, or heating and recondensation events, then the ice thickness Δa follows some constant t and hence the volume fractions change across the size distribution. However, if the grain is present in an environment where the grain size distribution has been reached by coagulation, shattering, or heating and recondensation events, then the limiting case (according to the assumption) leads to an ice thickness Δa proportional to the size of refractory material, t · a, thus preserving the ice/refractory volume fraction. The following powerlaw approach is used to include both limiting cases $$\text{Δ}a=t\cdot a^q$$(3)

where, t and q are arbitrary coefficients that are indicative of local ice abundance and the local dominance of shattering, coagulation, heating and recondensation events respectively. The two limiting cases are (1) q = 0, where the thickness of the ice layer is constant and (2) q = 1, where the ice/refractory volume fraction is constant and hence the thickness linearly depends on the bare grain size. The values 0 < q < 1 relate to cases where the grain size distribution is affected by coagulation, shattering, or heating and recondensation events in varying levels. We assume that f(a) in Eq. (2) is the size distribution of the underlying bare grains.

The effect of the ice covering the dust on the size distribution (dust + ice) is shown in Fig. 1. For a constant ice thickness (q = 0), the relative increase in size of the smaller grains is several orders of magnitude and the new minimum grain size is determined by the ice layer thickness. The number density of grains of size close to the ice thickness increases significantly. However, the change in size for larger grain sizes is negligible. The overall effect on the grain size distribution is that it does no longer follows a powerlaw. For the other limiting case q = 1, the small grains have small ice thicknesses and large grains have large ice thicknesses. Both minimum and maximum grain sizes are changed, but this is not as significant as in the case of constant ice thickness. Further, the total grain surface area increases, still dominated by smaller grains, but the size distribution powerlaw is retained.

The total volume of ice at any point in the disk V_I is obtained from the chemical output of the disk code. If n_i [cm⁻³] is the number density at a grid point, µ_i [g] is the molecular mass, ρ_i [g cm⁻³] is the material mass density of ice species i, and V_i is the volume of ice species i, then the total volume of ice per unit volume of space, V_I, at a grid point is given by Eq. (4) for N ice species. $$\begin{array}{l}{V}_i=\frac{n_i{\mu }_i}{{\rho }_i}\\ V_{\text{I}}={\displaystyle \sum _{\text{i}=1}^N{V}_i}\end{array}$$(4)

Using Eq. (4), Eq. (2) is solved numerically for t at each grid point in the disk using Newton-Raphson method by approximating the root of the equation V_IG − V_I − V_BG = 0, for a given ice power law, q.

Fig. 1 Icy grain size distribution for limiting values of ice coating powerlaw, q. Bare size distribution: amin = 0.05 µm, amax = 3 mm, powerlaw index 3.5, VBG = 6.2894 × 10−13 cm3 cm−3, VI = 2.7577 × 10−12 cm3 cm−3.

Fig. 2 Extinction, absorption, and scattering opacity of bare grains (dashed) and ice coated grains (solid). A constant thickness of ice layers is assumed for all grains (q = 0). The composition of grains is Mg0.7Fe0.3SiO3, amorphous carbon, and porosity in the ratio 0.60:0.15:0.25, and that of ices: H2O, CH4, NH3 in the ratio 0.60:0.20:0.20. The volume of dust is VBG = 6.2894 × 10−13 cm3 cm−3 and thickness of ice is Δa = 0.246 µm. Refractory dust grain size distribution parameters are amin = 0.05 µm, amax = 3 mm, and powerlaw index 3.5. These opacities have been calculated using effective medium theory.

3.2 Effect of ice on extinction, absorption, and scattering efficiency

Extinction $\left(Q_{\lambda }^{\text{ext}}\right)$, scattering $\left(Q_{\lambda }^{\text{sca}}\right)$, and absorption $\left(Q_{\lambda }^{\text{abs}}\right)$ efficiencies are related as given in Eq. (1). These efficiencies are complex functions of size parameter (x = 2πα/λ) and complex refractive index m(λ).

The default Mie-efficiency calculations in PRODIMO, as explained in Sect. 2, assume that the same grain material is present everywhere in the disk (i.e. position-independent). In this work, however, at all locations in the disk where the local physical state allows for the presence of ice on grains, the opacity calculations include the increase of size and change of material composition of the grains due to ice formation. Two factors that should be considered for efficiency calculations are the different ice species and their abundances. The former affects the complex refractive index and the latter affects both the size and the refractive index of the icy grains. Figure 2 shows the effect of these two factors on extinction, absorption, and scattering opacities by assuming a grain size distribution and fixed ice thickness (q = 0) and composition. The figure shows the opacity cross sections per gram of refractory material [cm² g⁻¹]. Scattering opacity increases both at shorter and longer wavelengths and absorption increases in the mid- to far-infrared range. This leads to an extreme albedo of >99% at the optical to near-infrared wavelength region. Estimating the effects of these two factors, the different ice species and their abundances, on the Mie efficiencies at each grid point are individually explained in the following subsections.

Fig. 3 Extinction efficiency for a grain size 0.0559 µm (with respect to the bare grain geometrical cross section) with different ice thicknesses composed of H2O, NH3, CO, CO2, CH4, CH3OH in ratio 0.45:0.20:0.10:0.05:0.05:0.15 with 20% of porosity. Also shown is the validity of the interpolation for a constant ice thickness layer of 0.515 µm. The embedded panel shows the relative error between the interpolation and explicit calculation of efficiency for the icy grain.

3.2.1 Effect of ice thickness on Mie efficiencies

To estimate the position dependent opacities, a Mie efficiency matrix (Q-matrix) is first generated. Each point in this multidimensional matrix corresponds to Mie efficiency defined by individual ice species, ice thickness, ice power law, bare grain size, and wavelength. For a given bare grain size distribution and the wavelength bins, a finite number of ice thickness and ice power law values (Eq. (3)) for each ice species are used to populate this matrix. Figure 3 shows how extinction efficiency changes as the thickness of ice layer increases. Since it is computationally expensive to calculate the efficiencies based on thickness at each grid point using Mie routines, log interpolation (Eq. (5)) between points in the pre-calculated g-matrix is done.

In Eq. (5), θ and ϕ are the interpolation weights which define the closeness of the local thickness to the nearest two thickness grid points. However, for interpolation between the grid points of the smallest ice thickness and bare grain opacities, a simple linear interpolation is used. Figure 3 presents the values obtained by such interpolation. It should be noted that in Fig. 3, the efficiency, Q_ext, is the ratio of extinction cross section of ice covered grains to the geometrical cross section of bare grains. Here, bare grain geometrical cross section is considered because it makes the log interpolation of efficiency simple. $$\begin{array}{l}{Q}_{\text{a},\lambda }(r,z)=Q_{\text{a},\lambda }{\left(t_1\right)}^{\theta }\cdot Q_{\text{a},\lambda }{\left(t_2\right)}^{\varphi }\\ \theta =\frac{\log (t_2/t_{r,z})}{\log (t_1/t_2)}\varphi =\frac{\log (t_1/t_{r,z})}{\log (t_1/t_2)}\end{array}$$(5)

where, r, z are the radial distance from the rotation axis and vertical distance from the midplane of a point in disk, t_r,z is the thickness coefficient (see Eq. (3)) of ice layer at the location (r,z) in the disk for which the opacity is being estimated, t₁ and t₂ are the two ice thickness values considered in the pre-calculated Q-matrix closest to t_r,z, using which the efficiencies Q_a,λ(t₁) and Q_a,λ(t₂) are obtained from the matrix, and θ, ϕ are the weights for log interpolation. Finally, the values calculated from Eq. (5) are scaled by $\left({\left(a\text{/}\overline{a}\right)}^2\right)$ to obtain the actual efficiencies with respect to ice covered grains ($\overline{a}$ is size of the grain with ice. i.e. $\overline{a}=a+\Delta a=a+t_{r,z}\cdot a^{q_{r,z}}$). However, the local ice thickness power law (q_r,z) requires extensive solutions of shattering, coagulation physics, and heating-recondensation events, which cannot be obtained by disk chemistry. Hence, trivial values (q = 0,1) are considered in Sect. 4 independent of position (i.e. same throughout the disk).

3.2.2 Effect of ice composition on Mie efficiencies

The protoplanetary disks are generally dominated by H₂O ice (Öberg et al. 2011). Though the direct observation of ices in disks are scarce, observations of molecular clouds and young stellar objects (YSOs) show the presence of H₂O, CO, CO₂, NH₃, CH₃OH, CH₄ amongst other ice species (Gibb et al. 2000, 2004; Boogert et al. 2008; Pontoppidan et al. 2008; Noble et al. 2013). These six ice species are considered for opacity calculations in this paper. Our sources for the refractive indices are listed in Table A.1. In this paper, the composition of ice at any point in the disk (the number density of ice species, see Eq. (2)) is obtained from the chemistry results of the disk code, which is explained in Kamp et al. (2017). The ice composition is determined by accretion and desorption processes (including thermal, UV, and cosmic ray processes). The effect of grain surface chemistry is discussed in Sect. 4.

It is found that the efficiencies based on local ice composition, for a given ice thickness and a given bare grain size distribution, can be approximated by a simple volume fraction weighted summation of efficiencies from the Q-matrix for grains coated with different pure ices. This is shown in Eq. (6). This method is computationally efficient because the number of Mie calculations required becomes independent of the spatial grid size of the disk, instead depends on the number of ice species, N, considered. $$Q_{\text{a},\lambda }\left({\text{Bare, Ice}}_1\text{,\hspace{0.17em}}{\text{Ice}}_2,\dots \text{,\hspace{0.17em}}{\text{Ice}}_n\right)={\displaystyle \sum _{i=1}^N\left({\eta }_i{Q}_{\text{a},\lambda }\left({\text{Bare+Ice}}_i\right)\right)}$$(6)

where Q_a,λ is the efficiency, the left hand side of the equation is the efficiency of the grains coated with local mixed ice, composed of species 1 to N, η_i is their respective ice abundance volume fraction (η_i = V_i/V_I, Eq. (4)). It should be noted that in Eq. (6), the Q_a,λ on left and right side of the equation are based on the same refractory grain composition and size. Further, Q_a,λ on the right is calculated using the same thickness of ice as in the left except that the ice layer is assumed to be entirely composed of Ice_i. The validity of this approach is shown in Fig. 4, which compares the efficiencies calculated for grain coated with a mixed ice layer (1) by using effective medium theory for ice mixture and (2) by weighted summation as given in Eq. (6). The relative error plotted in the bottom panel of the figure shows that the interpolation is in good agreement with calculation using effective medium Mie routine. These approximations (Eqs. (5) and (6)) introduces uncertainties (of similar order as the relative errors) in the ice features in the disk spectra obtained from the solution. However, these uncertainties are ~1 order of magnitude smaller than the uncertainties introduced due to difference in the opacities of a mixed and a differentiated core-mantle grain. Hence, these approximations are sufficiently accurate for our models but the uncertainties have to be taken into account while fitting observations. The following section discusses the effective medium and core-mantle Mie theory.

Fig. 4 Extinction efficiency for a grain size 0.0559 µm with 0.515 µm ice thickness composed of H2O, NH3, CO, CO2, CH4, CH3OH in ratio 0.16:0.16:0.16:0.2:0.16:0.16 with 20% of porosity. The bottom panel shows the relative error between the interpolation and explicit calculation of efficiency for the icy grain.

3.3 Effective medium vs. core-mantle Mie theory

The composition of the icy grains can be thought of as (1) a mixture of ices and refractory material or (2) a multi-layered structure with refractory material as the core and an ice mantle on top. The opacity calculations for these two cases can be carried out using either Mie theory with effective medium theory or core-mantle Mie theory, respectively. However, the computational effort required for core-mantle solutions is a factor 4 larger than the other.

It can be argued that in the limiting case q = 1 , the coagulated or shattered grains do not retain a core-mantle structure. Further, it is debatable whether icy grains are layered or mixed in dense environments (Noble et al. 2017). For example, ices could be organized as layers starting with a layer of polar ices such as H₂O, NH₃, and CH₃OH and on top a layer of non-polar ices such as CO, N₂, and O₂, creating an onion-like structure (Öberg 2009). Ice covered grains can undergo processing due to thermal processes, ultraviolet photons, and cosmic ray particles, called energetically processed ice (Boogert et al. 2015), which can mix the ices in the mantle and as well cause segregation. Impinging energetic particles and irradiation with X-rays and UV photons may cause ice desorption (Palumbo et al. 2006; Andersson & van Dishoeck 2008), mantle explosion (Ivlev et al. 2015), or the formation of new ice species (Allamandola et al. 1988; Muñoz Caro et al. 2002; Palumbo et al. 2006), possibly followed by molecular ejections (Dulieu et al. 2013) which can lead to mixed ices. The mixed ice components can undergo processing leading to ice restructuring and segregation (Fayolle et al. 2011; Öberg et al. 2009a). Further, thermal and radiative processes can also form larger and more complex molecules (Garrod et al. 2008; Gerakines et al. 1996; Öberg et al. 2009b; Muñoz Caro et al. 2002). In the present model, neither ice segregation nor restructuring is considered. Moreover, using Eq. (6) does not take into account the ice features that result from interactions between different ice species, for example CO-H₂O at 4.68 µm, 4.647 µm (Tielens et al. 1991; Sandford et al. 1988), NH₃-H₂O at 3.4 µm (Bertie & Morrison 1980). Figure 5 shows the ice opacities calculated for the limiting ice powerlaws with effective medium theory and core-mantle Mie theory (here, the mantle is a composed of a mixture of ices). The differences between the opacity values calculated by these two approaches are small. Due to these reasons our opacity calculations are based on effective medium theory.

Fig. 5 Extinction opacity of bare grains (green) and ice covered dust grains calculated using the effective medium theory (red) and Core-Mantle Mie theory (blue) for ice powerlaws q = 0 (solid lines) and q = 1 (small dashed lines). The ice composition is H2O:NH3:CH4:Vacuum = 0.51:0.40:0.04:0.05. The bare grain size distribution, composition, and volume of ice is the same as in Fig. 2.

4 Impact of ice opacities on the disk models

In the following, we discuss four protoplanetary disk models with identical stellar, disk shape and mass, dust settling, and bare grain parameters, see Appendix B. The models only differ with respect to the treatment of the ice opacities and grain surface chemistry. All four models use steady state chemistry and the fourth model includes grain surface chemistry. The four models are:

• bare grain: bare grain opacity.

• const ice thickns: ice opacities with a constant ice thickness, q = 0 (Eq. (3)).

• const ice vol frac: ice opacities with a constant ice volume fraction, q = 1 (grain size dependent ice thickness).

• const ice thickness + surf chem: same as the constant ice thickness model, that is, q = 0, but includes grain surface chemistry.

The position-dependent icy grain opacity calculations are based on the abundances of ice species as provided by the chemical network at each point in the disk. The iteration between chemistry, opacity calculation, and continuum radiative transfer, as described in Sect. 2, is continued until the opacity computations converge to a tolerance of less than 5%.

4.1 Extinction and opacity

The inclusion of the ice opacities significantly changes the extinction properties of the disk in the model, particularly in the midplane, where ices are abundant. Figure 6 shows the vertical visual extinction, radial visual extinction, dust temperature, gas temperature, spectral energy distribution, and disk mass averaged opacity of the entire T Tauri disk (in cm² g⁻¹ (refractory material)) for the four disk models. The assumptions of constant ice thickness (q = 0) and constant ice volume fraction (q = 1) modifies the grain size distribution differently (Fig. 1), which in turn affects the mean disk opacity. Panel F of Fig. 6 shows that including ice opacities significantly increases the disk mean opacity compared to the bare grain opacity (black) model. The constant ice thickness model (red) shows an increase in the opacity at shorter wavelengths (<3 µm) by more than an order of magnitude. At longer wavelengths (>100 µm), the opacity change with respect to bare grain model is negligible. The mean opacity of the constant ice volume fraction model (blue) is larger than the bare grain model at all wavelengths. Particularly, at shorter wavelengths the increase in opacity is by a factor ~2 and at longer wavelengths by a factor of ~1.3. At intermediate wavelengths, the opacities of both constant ice thickness and constant ice volume fraction models are of the same order of magnitude. However, the constant ice volume fraction model produces weaker ice features. The model with constant ice thickness and grain surface chemistry (green) shows that an increase in the ice abundance (due to grain surface reactions) leads to a larger opacity compared to all the other models in the mid-IR wavelengths. For a given volume of ice at one location in a disk, the effect of the ice thickness to grain size relation on the opacity can also be seen in Fig. 5. The difference in the opacities between the models is an interplay of changing grain size distribution due to the ice thickness to grain size relation and the absolute abundance of ice due to chemistry which is explained in the following paragraph.

The model with constant ice thickness produces ice features that are stronger than for the model with constant ice volume fraction by an order of magnitude. However, the ice abundances are similar in these models. This implies that ice distribution across the grain sizes become important for the ice features. Figure 7 shows the contribution of different grain sizes to the disk mean opacity for the four models. Grains of particular sizes (a) have opacity maxima at wavelengths regions corresponding to 2πa/λ ~ 1. At wavelengths corresponding to 2πa/λ ≫ 1, the Mie extinction efficiencies reach a constant value (Q_ext → 2). At wavelengths corresponding to 2πa/λ ≪ 1, the extinction is dominated by absorption (Q_ext ≃ Q_abs) and falls off inversely proportional to the wavelength (Q_abs ∝ λ⁻¹, Krugel 2002). This behavior is seen in the bare grain model (left). For the assumed refractory grain size distribution, the opacity in the wavelength range 2 µm–100 µm is dominated by grains of sizes 1.0 µm ≤ a < 10 µm and 10 µm ≤ a < 100 µm. Further, the small grains (<1 µm) have strong silicate (and ice for other models) features. In the constant ice thickness model, small grains (<1 µm) have large volume fractions of ice (~90%) due to larger surface area. This leads to change in grain size distribution as shown in Fig. 1, where there is a pile up of grains around the size corresponding to ice thickness and removal of grains smaller than the ice thickness. Due to these changes, the opacity contribution of grains smaller than local ice thickness drastically reduces (~1 order) and that of grains of size corresponding to the ice thickness drastically increases (~1 order, Fig. 7). That is, the main opacity carriers in the infrared wavelengths changes from 1.0 µm ≤ a < 100 µm to 0.1 µm ≤ a < 10 µm. This effect is more pronounced in the constant ice thickness model with grain surface chemistry (right) in which the relevant grain sizes are 0.1 µm ≤ a < 1 µm. In contrary, in the constant ice volume fraction model most of the ice resides on larger grains (>100 µm) which do not contribute to opacity in the near- and mid-infrared wavelengths and the change in size distribution is uniform across grain sizes. Thus, models with constant ice thickness show stronger ice features. The opacity curve of the bare grain model shows two silicate features at 10 µm and 20 µm. The constant ice volume fraction model shows small features around 2.95 µm (NH₃), 3 µm (H₂O), 4 µm (CO₂), 26 µm (NH₃), and the silicate features. Due to the changes in the size distribution, the models with constant ice thickness (with and without grain surface chemistry) show significantly enhanced ice features, prominently showing H₂O, NH₃, CO₂, CO, and CH₄. Hence, the average opacity of the disk strongly depends on 1) the ice power law, q, in Eq. (3) and 2) total abundance of ices.

Panels A and B show the effect of the ice opacities on the radial and vertical visual extinctions. When ice opacities are included, the vertical extinction increases significantly toward the midplane, but only once the disk has become optically thick. This increase is more than an order of magnitude in the models that use a constant ice thickness (in which the opacity increases significantly as well) closer to the midplane, like an opacity wall or barrier. For example, the vertical visual extinction at 30 au in the midplane changes from about 10 to 100 when q = 0 ice opacities are considered. The observable disk extent can be defined as the vertical optical depth be equal to unity at optical wavelengths. This observable disk extent increases from about 120 au for the bare grain model to about 300 au for the models with constant ice thickness. This behavior is mostly a size effect, because ice formation increases the sizes of the small grains in particular, which are the carriers of the opacity at optical wavelengths. However, the impact of the icy opacities in the q = 1 model is less pronounced. Hence, the visual extinction profile in a disk strongly depends on the relation between the ice thickness and grain size, that is, the ice power law index q in Eq. (3).

Fig. 6 Comparison of disk properties of the models. Panels A–F show the vertical visual extinction, radial visual extinction, dust temperature, gas temperature, spectral energy distribution, and disk mass averaged opacity of the entire T Tauri disk (in cm2 g−1 (refractory material)) respectively for four disk models. The four models are (1) bare grains, (2) icy grains with constant thickness of ice layer (q = 0), (3) icy grains with constant ice/refractory volume ratio (q = 1), (4) icy grains with constant ice thickness (q = 0) with grain surface reactions. All the models use steady state chemistry.

Fig. 7 Contribution of different grain sizes (refractory + ice) to the disk mean opacity for bare grain model, constant ice thickness model, constant ice volume fraction, and constant ice thickness with grain surface chemistry model. The legend at the top indicates the curves corresponding to different grain size bins.

4.2 Spectral energy distribution and temperature profile

Although the mean opacities in the disk increase drastically when ice coating is included, the main disk properties such as temperature, pressure, and chemical abundances, etc., do not change significantly in the models. This is because the ice only builds up in the (radially and vertically) optically thick and shielded regions in the disk. In particular, there is almost no ice in the model at locations where stellar or interstellar UV photons can penetrate, even if the local dust temperatures are low. This is because of the UV photodesorption process, which removes ice species as they interact with UV photons. The changes of opacity in the optically thick regions have no significant influence on the local mean intensities (J_v), as J_v ≈ B_v(T_d) is valid anyway (B_v is the Planck or black-body radiation), and hence there is no significant influence on the radiation that escapes from these regions. Panels C and D of Fig. 6 show the computed dust and gas temperatures in the four disk models. The surface layer of the disk re-emits or scatters radiation received from the star toward the observer and toward the disk midplane (Natta 1993). The latter largely determines the temperatures in the disk (Chiang & Goldreich 1997; Chiang et al. 2001). Since the ice is formed in regions that are already optically thick, there is almost no ice in the upper regions of the disk that can scatter incoming radiation. Thus, the temperature structure does not change much when position dependent ice opacities are included. For example, the largest change in temperature at any grid point compared to the bare grain model is 18%; this is for the model with grain surface chemistry. The overall temperature changes are less than 1.5% for the models without grain surface chemistry and less than 5% for the model with grain surface chemistry.

The spectral energy distributions (SEDs) of the four models are presented in panel E of Fig. 6 for disk inclinations 0°, 45°, 75°, and 85°. Since the visual extinction and the albedo in the disk surface, and the temperature in the disk largely remain unaffected, we do not see substantial change in the shape of the SED. However, at higher inclinations due to large optical depths in the midplane, the flux density at shorter wavelengths (<3 µm) for the models with constant ice thickness is about an order of magnitude smaller compared to the bare grain model. Similarly, the increase in opacity at longer wavelengths for the model with constant ice volume fraction (for example, the opacity at 10³ µm is 1.17 times that of bare grain model) results in slightly lower flux density in the longer wavelength region of the SED. This is encapsulated in Fig. 8 which shows the spectral indices (n_a–b) for 1–4 µm, 35–300 µm, and 0.85–1.3 mm for the four models at different disk inclinations. The spectral index is defined by the following relation:

$$n_{a-b}=\frac{{\log }_{10}\left({\lambda }_a{F}_a\right)-{\log }_{10}\left({\lambda }_b{F}_b\right)}{{\log }_{10}\left({\lambda }_a\right)-{\log }_{10}\left({\lambda }_b\right)}$$(7)

where λF is the flux density and λ is the wavelength. The three spectral indices probe different regions of the disk. n_{1–4 µm} corresponds to the hot inner disk surface, n_{0.85–1.3 mm} corresponds to the cooler outer disk midplane, and n_{35–300 µm} corresponds to warm disk regions. The spectral range for the latter is large because multiple ice features are present between these wavelengths and determination of the slope becomes difficult for smaller wavelength ranges, see Fig. 9.

The spectral indices of the models at near- to mid-infrared wavelengths (n_{1–4 µm}) remain constant across the models for a given inclination. This is because the radiation at these wavelengths is emitted from the surface layers close to the inner rim of the disk. These regions are devoid of ices due to thermal desorption. At millimeter wavelengths, the radiation arises from large part of the disk and hence can be used to estimate the dust mass, size of the grains, and their evolution in the disk (Natta et al. 2007; Testi et al. 2014). Thus, spectral indices at millimeter wavelengths (n_{0.85–1.3 mm}) are used to quantify the dust properties (e.g. Woitke et al. 2019; Villenave et al. 2020). The spectral indices for any given model depends on the inclination (Fig. 8). In other words, the model with the same dust properties exhibit different mm-spectral indices at different inclinations. Secondly, for a given disk inclination, the spectral index at millimeter wavelengths can depend on the ice opacity and chemistry model. Since the radiation at these wavelengths arise from the bulk of the outer and cooler disk, the disk mass averaged opacity can be used to understand the behavior of the spectral indices. Panel F of Fig. 6 shows that the model with a constant ice volume fraction has a steeper opacity curve at these wavelengths and the model with a constant ice thickness has essentially the same opacity curve as the bare grain model. This is also reflected in the spectral indices at lower inclinations (0°, 45°). However, at higher inclinations, the orientation effects play a role leading to certain parts of the disk being less visible and some being more visible. At mid- to far-infrared wavelengths, the spectral index (n_35–300µm) becomes more negative (steeper slope) for models with ice opacities without grain surface chemistry. The model with constant ice volume fraction has a lower spectral index than the model with constant ice thickness at all disk inclinations. This shows that the ice power law has an effect on the slope of mid- to far-infrared SED. The model with constant ice thickness with grain surface chemistry shows the spectral indices similar to that of the model without grain surface chemistry. The emission region corresponding to these wavelengths spans across large range of temperatures (due to large wavelength range), hence it is difficult to qualitatively analyze the complex interaction of ice opacities and the spectral index. The indices for the model with grain surface chemistry for disk inclination of 85° are unreliable due to substantial ice features, hence are not shown. The ice power law determines the change in grain size distribution (Fig. 1) and the distribution of ice across the grain sizes. Thus, the ice power law influences the slope of the SED at longer wavelengths.

Fig. 8 Spectral indices for 1–4 µm (left), 35–300 µm (middle) and 0.85–1.30 mm (right) for the four models presented in Sect. 4. Different markers represent the spectral indices for different disk inclinations as indicated in the legend on the right.

4.3 Ice emission and absorption features

The models produce ice absorption and emission features in the SED. At inclinations lower than the edge-on case, the SEDs do not show any ice absorption features in any of the four models. However, weak emission features at far-infrared wavelengths can be observed corresponding to H₂O, NH₃ ice, and CO₂ ice. Close to edge-on (85°), models with constant ice thickness show ice absorption features. Particularly, the model without grain surface chemistry shows only a strong CO ice absorption feature and a weak NH₃ ice absorption feature at 4.6 µm and 26 µm respectively. The model with grain surface chemistry shows significantly more ice absorption features (and stronger) in addition to those seen in previous model, such as H₂O, CO₂, and CH₃OH. This model also shows stronger ice emission features at far-IR wavelengths. This shows that the chemical network strongly affects the strength of the ice absorption and emission features. Therefore, observations of such ice features can lay strong constraints on the type of chemistry in disks.

Figure 9 shows the zoomed-in SED (top panels) and optical depths (bottom panels) of the ice emission features (at 75° inclination) from 40–300 µm for the models discussed before (Sect. 4). The black curves in the SED plots in all the top panels show the pseudo continuum used to obtain the optical depths of the ice features. This continuum is obtained by spline interpolation using python library scipy (Virtanen et al. 2020) and then correcting spline of other models for systematic error produced by the spline for bare grain model. The top and bottom plots in Panel A show the SED and optical depth of the four models. Due to higher ice abundances, the model with grain surface chemistry produces stronger features (by about a factor 2 in optical depth) than the model without grain surface chemistry. The peaks of the optical depths of emission features for the constant ice thickness and the constant ice volume fraction models are quite similar, however, the constant ice volume fraction model shows larger optical depths (by a factor less than 2) between the peaks of these features. The figure also shows the ices that correspond to the features, the prominent ices being H₂O, NH₃, CO₂, and CH₃OH. Unlike other ices, CH₃OH ice has a broad feature that spans from about 50 µm to 160 µm and hence does not show up as a single peak instead contributes to the optical depth across the wavelength range. Panel B of Fig. 9 shows that the optical depths are quite similar across the inclinations except for the close to edge-on case, where some of the features (such as 47 µm, 70 µm H₂O and NH₃ ice features) can be seen in absorption. Thus, the strength of the ice features depend on the chemistry producing the ices, the relation governing ice distribution over different grain sizes and the disk inclination.

The SED shows absorption features only at high disk inclinations. Hence, it is interesting to explore the observability of absorption features when the disk is spatially resolved. Figure 10 shows images of the four models at 26 µm in the top panels. The bottom of the figure shows the spectral flux density of these models obtained from the region marked by the cyan rectangles. The cyan rectangles are 0.2″ × 0.3″ in size². This is a simple study to estimate the fluxes for spatially resolved observations and is not an attempt to simulate any specific instrument. Hence, point spread function, noise, errors, sensitivity, and saturation limits are not considered. Simulating observations and identifying best observation strategies are left to a future work.

The images in the top panels show that inclusion of ice opacities produces a distinct band of absorption around the midplane, subsequently called the dark lane. This dark lane is more prominent for models with constant ice thickness and also depends on the chemical network. These differences reflect the increase in the extinction in the disks as discussed in Sect. 4. Panel A of Fig. 6 shows that the inclusion of ice opacities increases the extinction radially, more strongly for models with constant ice thickness compared to model with constant ice volume fraction. The ice power law has a strong effect on the appearance of the disks, particularly on the size of the dark lane. The spectral flux density (in the bottom panel) extracted in the dark lane shows that ice absorption features can be seen by spatially resolving the disks at inclinations less than the edge-on case. Further, the constant ice volume fraction model shows substantially weaker features compared to the models that consider a constant ice thickness. The absorption coefficients of ices and the refractory silicate used are also shown at the bottom of the figure to help identify the features observed in the spectra. More information regarding the optical constants used in this paper is described in Appendix A. The models show features corresponding to five of the six ices included in the opacity calculation. Though the local abundance of CH₄ ice can be as large as 1.4 × 10⁻⁴ in the disk, no feature corresponding to it is observed. This is because, though the local CH4 abundance can be large, the column density along the line of sight of the spectra extracted in Fig. 10 for other ices such as H₂O, NH₃ and CH₃OH are larger due to the spatial extent of these ices (see appendix). Moreover, the CH₄ features lie in the wavelength regions with features from more dominating ices. This can also be seen in the disk mean opacity of the model with grain surface chemistry as shown in panel F of Fig. 6. The model with grain surface chemistry shows stronger absorption features than the model without grain surface chemistry by about 2 orders of magnitude. Ice emission features are also present at longer wavelengths but are substantially weaker (by several orders of magnitude) compared to the absorption features. Hence, spatially resolving the disk produces significant ice absorption and emission features that are not seen in the disk integrated spectra, particularly strong in the absorption by few orders of magnitude.

Figure 11 shows the optical depth of the 3 µm ice feature obtained from spatially resolved spectral density plot of the model with constant ice thickness from Fig. 10, along with the absorption coefficients of different ice species that have a feature in that wavelength range. The 3 µm profile is a combination of H₂O, NH₃, and CH₃OH ices. The sharp peak at 2.95 µm in the profile is due to NH₃ ice. H₂O ice dominates the profile of the 3 µm feature, however, the longer wavelength side (>3.3 µm) shows a contribution of CH₃OH ice. NIRSpec IFU should be able to spatially resolve the disks and hence we can expect observations of such features from JWST.

Fig. 9 Spectral energy distributions and optical depths of ices in different models. The top panels show the spectral energy distributions of different models between 40 µm and 300 µm. The black curves indicate a pseudo-continuum used to obtain the optical depths of the ice features. In each panel, the SED of the different models are arbitrarily shifted to provide better visualization. The bottom panels show optical depths of the ice features in the SEDs in the top panels. Panel A shows the ice features in the four models at a disk inclination of 75°. Panel B shows the ice features for the model with constant ice thickness without grain surface chemistry at different disk inclinations. Note that the orange curve in each panel is the same model. The ice features are marked using colored rectangles in the optical depth plots.

Fig. 10 Simulated images and spatially resolved spectral flux densities of the disk models presented in this work. The top panels show simulated images (intensity maps) at 26 µm and inclination of 75°. The cyan rectangles (0.2″× 0.3″) indicates the region whose spectral flux density is plotted in bottom panel. The spectral flux density plot also shows the absorption coefficients of the ices and silicates.

5 Discussion

Ice accretion on the grains modifies the grain size distribution. In Sect. 3.1, we showed that a constant ice thickness leads to a deviation of the size distribution from the powerlaw and a constant ice volume fraction increases the minimum and maximum grain sizes. Further, these changes can be enhanced by models that produce more ice, such as the model with grain surface chemistry. The emission at (sub)millimeter wavelengths probe the bulk of the dust in the disk. This corresponds to the large dust grains which are settled in the disk midplane. In these environments grain evolution processes such as grain growth are important as they affect the emission at these wavelengths. Natta et al. (2007) show that the millimeter slope is strongly affected by the maximum grain size and the dust distribution powerlaw. Particularly, when the maximum grain size is of the order of millimeters, the slope is very sensitive to the changes in maximum grain size. Since, the emission in these wavelengths can be used to quantify the dust properties and their evolution, it is important to consider the effect of ice on the grain size distribution and the opacity while interpreting these spectral indices.

Ices have been detected in protoplanetary disks, as discussed in Sect. 1. These detections are from the total disk integrated fluxes (SED) and are not spatially resolved. Our models show weak emission features at far-IR wavelengths at all inclinations, which can be, however, very weak and barely observable for some model setups. Absorption features can only be seen in the edge-on orientation in the total disk spectra. Even when seen edge-on, most of the absorption features are seen only when grain surface reactions are included, for example, the 3 µm H₂O feature. However, spatially resolving the disk produces a multitude of absorption features across the spectral range.

In contrary, the 3 µm water ice feature has been reported in total disk flux of some objects (mostly edge-on) such as HK Tau B (Terada et al. 2007), d132-1832, d216-0939 (Terada et al. 2012). Further, emission features are observed in a few protoplanetary disks (one such example is VW Cha, McClure et al. 2015). The presence of emission features and absence of absorption features in our models can be attributed to the limited parameters explored. Objects differ, for example in terms of the parent star, disk mass, the inclination, the flaring angle, inner and outer disk radii, gaps, silhouette disks, etc. All these parameters can strongly affect the ice features. Consequently, each disk has to be modeled separately and our models are not representative of every disk. Further, dynamic physical processes such as vertical and radial transport are not included in our models and can significantly alter the ice features. For example, vertical mixing can significantly change the ice composition in the disk (Furuya & Aikawa 2014) and readily produces ice features that are stronger compared those from our models by factors of a few in terms of optical depth (Woitke et al. in prep.). Similarly, models with gaps can also show strong ice features in the SED. The goal of this paper is to present the consistent position dependent opacities, and exploring such parameters and simulating observations is left to future work. For example, the disks with ice features in the far-IR reported by McClure et al. (2015) can be modeled with our ice opacity approach.

Ballering et al. (2021) also modeled ices in disks and they find absorption features of multiple ices in the total disk flux even at inclinations lower than the edge-on case. However, these models differ from our models significantly and hence a direct comparison is not possible. For example, Ballering et al. (2021) do not consider the interstellar UV field and use a cosmic ray ionization rate that is an order of magnitude lower than our models. These assumptions can significantly affect the ice features. The vertical extent of ice in the disk is largely determined by the interstellar UV field, because this radiation can photodesorb the ices back into the gas phase. This in turn inhibits the presence of ice in the emitting regions of the disk. The ice destruction mechanism by photodesorption is consistently treated in our models, with chemistry-dependent ice opacities and fully coupled UV radiative transfer. Further, Ballering et al. (2021) use two dust populations with significantly different settling conditions and the treatment of ice opacities is not fully consistent with local conditions in each grid point of the model.

The strengths of the ice features strongly depend on the ice power law and the chemistry network. Spatially resolving the disk presents multiple ice features and can probe the vertical and radial distribution of ice in the disk. In the end, future observations will allow us to test these various models and to refine the assumptions.

Fig. 11 Optical depth of the 3 µm H2O ice feature derived from the spectral flux density obtained from model 2 in Fig. 10. The optical depth is shown in dark gray thick line and the absorption coefficients of ice species that have a feature in the plotted wavelength range are shown as colored thin lines.

6 Conclusions

A computationally efficient way of computing position-dependent ice and bare grain opacities in protoplanetary disks has been implemented in PRODIMO, which consistently couples chemistry, ice formation, opacity computations, and radiative transfer in protoplanetary disks. The conclusions drawn are the following:

1. A weighted sum of the opacities of grains coated with the different pure ices and log interpolation of opacities between pre-calculated points in ice thickness space provide a computationally efficient and sufficiently accurate method to calculate the icy opacities in protoplanetary disks.

2. Locally, the continuum opacity can change significantly by ice formation. The opacity changes are particularly large at shorter wavelengths, due to the increased size of the small grains, when assuming a constant ice thickness on all grains. Opacity changes are less pronounced, and more uniform in wavelength space, when the grains are assumed to have a constant ice/refractory volume ratio, which is expected if coagulation and shattering are important. In both cases, the primary reason for those opacity changes is the altered size distribution of the icy grains.

3. The position-dependent ice opacities do not alter the physical disk structure significantly, because the ice opacities are significant only in the optically thick disk regions (A_V ≳ 10) in our models.

4. The spectral energy distributions (SEDs) of the disks show faint ice emission features in the far-infrared regions beyond 35 µm for all our models with ice opacities. However, absorption features are seen only at high inclinations, close to edge-on conditions. In particular, the 3 µm water ice feature is seen only in case of the model with grain surface chemistry in close to edge-on orientation.

5. The assumption made on how the ice is distributed across the grain size distribution (ice power law) influences the far-infrared and millimeter slope of the SED.

6. The strengths of the ice features are influenced by the ice power law (about a factor of 2 in optical depths). The number of ice features and their strength increase (a factor 2 in optical depth) in models having grain surface chemistry, which is directly related to an increase in ice abundance.

7. Our models show a prominent dark lane in resolved images. This dark lane is caused by icy grains situated just above the disk midplane, and the size of the dark lane depends on the assumptions used to calculate the modified dust size distribution and the chemical network.

8. To detect the ice absorption features, extracting spatially resolved spectra from the dark lane is more promising.

The models presented here are useful to model existing (e.g. Herschel) and upcoming observations (e.g. JWST) to better understand the ice composition in disks. And further improvements, such as transport processes, will help us understand the role of these processes in disk environments.

Acknowledgements

CHR acknowledge the support of the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 325594231. CHR is grateful for support from the Max Planck Society.

Appendix A Ice Optical Constants

The source of optical data, the measurement temperature, and structure of six types of ice species and refractory material used are tabulated in Table A.1. Ice structures change with temperature, so does the refractive index and absorption coefficient. Though the abundances of the different ices follow the local thermo-chemical state of the disk, the optical constants used are not consistent with temperature. Table A.1 summarizes the state and temperature of ices at which the optical data is obtained. It should be noted that the ice optical data do not cover the entire spectrum (0.0912 µm – 10⁴ µm). Hence, the data has to be expanded extensively for all species to enable studying all the ice features. For absence of optical data corresponding to a wavelength λ which is bounded by optical data at λ₁ and λ₂ (such that λ₁ < λ < λ₂), we perform a log interpolation to estimate the optical data at λ. In an unbounded case (extrapolation), for longer wavelength side (λ > λ₁) we assume: n(λ) = n(λ₁) and k(λ) = k(λ₁) · λ₁/λ; for shorter wavelength side (λ < λ₂) we assume: n(λ) = n(λ₂) and k(λ) = k(λ₂). The optical constants used are shown in Fig. A.1 and Fig. A.2.

Fig. A.1 Refractive indices (n) of ices species used in the models.

Fig. A.2 Absorption coefficients (k) of ices species used in the models.

Appendix B Disk Model

We present a fiducial T Tauri disk ProDiMo model. Table B.1 presents the disk properties considered for all the models presented in this paper. For more details regarding the parameters, meanings and their values, refer Woitke et al. (2009) and Woitke etal. (2016).

We use the large DIANA chemical network composed of 235 species with 13 elements. The network is described in detail in Kamp et al. (2017). The adsorption energies of different ices considered in this paper are given in Table B.2. The grain surface chemistry used in our models is explained in (Thi et al. 2020), we add two species on top of the large DIANA network: H# and H₂#, allowing for a minimal surface chemistry, e.g, hydrogenation on surfaces. These models only use physisorption and no chemisorption for the surface chemistry. The purpose of including the surface chemistry is to assess the impact of different chemical networks, and is not to explain the ice observations and the compositions in these objects. We also ran time-dependent chemistry model with grain surface chemistry where the initial composition was determined by a molecular cloud model. We evolved the model up to 2Myr. We noticed that the ice abundances, hence the opacities and the ice features, are very similar to that of the steady state chemistry with grain surface reactions, hence the model is not included in Sect. 4.

Fig. B.1 Distribution of the six ice species (H2O, CO, CO2, NH3, CH3OH, and CH4) in the presented four models with different ice opacity (q) and chemistry settings (See Sect. 4). Each row corresponds to the different ice species and each column corresponds to the one of the four models. The dust temperature and vertical visual extinction contours are also shown with solid and dashed lines respectively.

Fig. B.2 Mean disk abundance of ice species in the four models.

Fig. B.3 Ice thickness distribution in the constant ice thickness model without grain surface chemistry. The green diamond indicates the location of largest ice thickness (0.61 µm)

Figure B.1 shows the abundances of H₂O, NH₃, CO, CO₂, CH₃OH, and CH₄ ices in the four models discussed in Sect. 4. The figure also shows the dust temperatures and vertical extinction in the disk models. Each column of panels in the figure shows the abundance of the ice species in each model. Models with ice opacities show an increase in water ice in the mid-plane which extends radially. The outer radii of water ice region extends from 7 au in bare grain model to 25 au, 10 au, and 400 au in models with constant ice thickness, constant ice volume fraction, and constant ice thickness with grain surface chemistry respectively. Fig. B.2 shows the mean abundances of the ices in the models presented in this paper. The models with ice opacities are generally more abundant in ices, with exception of CO and CH₄ ices. The ice abundances are slightly higher in the model with constant ice thickness than the model with constant ice volume fraction by about a factor 1.5. The model with grain surface chemistry has several orders of magnitude higher mean ice abundances, except for CO ice which is about two orders lower. In general, including ice opacities results in higher ice abundances except for CO ice. It should be noted that these are mean abundances of the entire disk, hence, essentially are model properties and should not be mistaken for the ice abundances that can be obtained from observations. Figures B.3 and B.4 show the distribution of ice thickness in the disk for the models with constant ice thickness without and with grain surface chemistry respectively. The mean ice thickness for entire the disk (which includes regions without any ice) in these models are 0.045 µm and 0.094 µm respectively. Figure B.5 shows the distribution of ice volume fraction across the disk for the constant ice volume fraction model. The mean ice volume fraction for the entire disk in this model is 24.2%. Unlike constant ice thickness models, the constant ice volume fraction model shows lower ice volume fraction in the midplane and larger ice volume fraction in the upper midplane. Dust settling leads to an increase in the large grain population in the midplane. Since ice largely resides on large grains in this model, the regions with little large grains yet having significant ice abundances have larger ice volume fractions.

Fig. B.4 Ice thickness distribution in the constant ice thickness model with grain surface chemistry. The green diamond indicates the location of largest ice thickness (0.59 µm)

Fig. B.5 Ice volume fraction distribution in the constant ice volume fraction model. The green diamond indicates the location of largest ice volume fraction (0.99)

References

All Tables

All Figures

	Fig. 1 Icy grain size distribution for limiting values of ice coating powerlaw, q. Bare size distribution: amin = 0.05 µm, amax = 3 mm, powerlaw index 3.5, VBG = 6.2894 × 10−13 cm3 cm−3, VI = 2.7577 × 10−12 cm3 cm−3.
In the text

Fig. 2 Extinction, absorption, and scattering opacity of bare grains (dashed) and ice coated grains (solid). A constant thickness of ice layers is assumed for all grains (q = 0). The composition of grains is Mg0.7Fe0.3SiO3, amorphous carbon, and porosity in the ratio 0.60:0.15:0.25, and that of ices: H2O, CH4, NH3 in the ratio 0.60:0.20:0.20. The volume of dust is VBG = 6.2894 × 10−13 cm3 cm−3 and thickness of ice is Δa = 0.246 µm. Refractory dust grain size distribution parameters are amin = 0.05 µm, amax = 3 mm, and powerlaw index 3.5. These opacities have been calculated using effective medium theory.

In the text

Fig. 3 Extinction efficiency for a grain size 0.0559 µm (with respect to the bare grain geometrical cross section) with different ice thicknesses composed of H2O, NH3, CO, CO2, CH4, CH3OH in ratio 0.45:0.20:0.10:0.05:0.05:0.15 with 20% of porosity. Also shown is the validity of the interpolation for a constant ice thickness layer of 0.515 µm. The embedded panel shows the relative error between the interpolation and explicit calculation of efficiency for the icy grain.

In the text

	Fig. 4 Extinction efficiency for a grain size 0.0559 µm with 0.515 µm ice thickness composed of H2O, NH3, CO, CO2, CH4, CH3OH in ratio 0.16:0.16:0.16:0.2:0.16:0.16 with 20% of porosity. The bottom panel shows the relative error between the interpolation and explicit calculation of efficiency for the icy grain.
In the text

	Fig. 5 Extinction opacity of bare grains (green) and ice covered dust grains calculated using the effective medium theory (red) and Core-Mantle Mie theory (blue) for ice powerlaws q = 0 (solid lines) and q = 1 (small dashed lines). The ice composition is H2O:NH3:CH4:Vacuum = 0.51:0.40:0.04:0.05. The bare grain size distribution, composition, and volume of ice is the same as in Fig. 2.
In the text

Fig. 6 Comparison of disk properties of the models. Panels A–F show the vertical visual extinction, radial visual extinction, dust temperature, gas temperature, spectral energy distribution, and disk mass averaged opacity of the entire T Tauri disk (in cm2 g−1 (refractory material)) respectively for four disk models. The four models are (1) bare grains, (2) icy grains with constant thickness of ice layer (q = 0), (3) icy grains with constant ice/refractory volume ratio (q = 1), (4) icy grains with constant ice thickness (q = 0) with grain surface reactions. All the models use steady state chemistry.

In the text

	Fig. 7 Contribution of different grain sizes (refractory + ice) to the disk mean opacity for bare grain model, constant ice thickness model, constant ice volume fraction, and constant ice thickness with grain surface chemistry model. The legend at the top indicates the curves corresponding to different grain size bins.
In the text

	Fig. 8 Spectral indices for 1–4 µm (left), 35–300 µm (middle) and 0.85–1.30 mm (right) for the four models presented in Sect. 4. Different markers represent the spectral indices for different disk inclinations as indicated in the legend on the right.
In the text

Fig. 9 Spectral energy distributions and optical depths of ices in different models. The top panels show the spectral energy distributions of different models between 40 µm and 300 µm. The black curves indicate a pseudo-continuum used to obtain the optical depths of the ice features. In each panel, the SED of the different models are arbitrarily shifted to provide better visualization. The bottom panels show optical depths of the ice features in the SEDs in the top panels. Panel A shows the ice features in the four models at a disk inclination of 75°. Panel B shows the ice features for the model with constant ice thickness without grain surface chemistry at different disk inclinations. Note that the orange curve in each panel is the same model. The ice features are marked using colored rectangles in the optical depth plots.

In the text

Fig. 10 Simulated images and spatially resolved spectral flux densities of the disk models presented in this work. The top panels show simulated images (intensity maps) at 26 µm and inclination of 75°. The cyan rectangles (0.2″× 0.3″) indicates the region whose spectral flux density is plotted in bottom panel. The spectral flux density plot also shows the absorption coefficients of the ices and silicates.

In the text

	Fig. 11 Optical depth of the 3 µm H2O ice feature derived from the spectral flux density obtained from model 2 in Fig. 10. The optical depth is shown in dark gray thick line and the absorption coefficients of ice species that have a feature in the plotted wavelength range are shown as colored thin lines.
In the text

	Fig. A.1 Refractive indices (n) of ices species used in the models.
In the text

	Fig. A.2 Absorption coefficients (k) of ices species used in the models.
In the text

Fig. B.1 Distribution of the six ice species (H2O, CO, CO2, NH3, CH3OH, and CH4) in the presented four models with different ice opacity (q) and chemistry settings (See Sect. 4). Each row corresponds to the different ice species and each column corresponds to the one of the four models. The dust temperature and vertical visual extinction contours are also shown with solid and dashed lines respectively.

In the text

	Fig. B.2 Mean disk abundance of ice species in the four models.
In the text

	Fig. B.3 Ice thickness distribution in the constant ice thickness model without grain surface chemistry. The green diamond indicates the location of largest ice thickness (0.61 µm)
In the text

	Fig. B.4 Ice thickness distribution in the constant ice thickness model with grain surface chemistry. The green diamond indicates the location of largest ice thickness (0.59 µm)
In the text

	Fig. B.5 Ice volume fraction distribution in the constant ice volume fraction model. The green diamond indicates the location of largest ice volume fraction (0.99)
In the text