The effects of dust evolution on disks in the mid-IR

TL;DR

This study models how dust evolution in protoplanetary disks affects mid-infrared gas emission lines, revealing significant increases in line fluxes for most species as dust grains grow and settle over millions of years.

Contribution

It couples dust evolution modeling with thermochemical disk simulations to analyze the impact on observable mid-infrared spectral lines, a novel integrated approach.

Findings

Line fluxes of most species increase with dust evolution.

C2H2 line fluxes remain low despite detection, indicating complex chemistry.

CO2 flux increases faster than other species during disk evolution.

Abstract

In this paper, we couple together the dust evolution code two-pop-py with the thermochemical disk modelling code ProDiMo. We create a series of thermochemical disk models that simulate the evolution of dust over time from 0.018 Myr to 10 Myr, including the radial drift, growth, and settling of dust grains. We examine the effects of this dust evolution on the mid-infrared gas emission, focussing on the mid-infrared spectral lines of C2H2, CO2, HCN, NH3, OH, and H2O that are readily observable with Spitzer and the upcoming E-ELT and JWST. The addition of dust evolution acts to increase line fluxes by reducing the population of small dust grains. We find that the spectral lines of all species except C2H2 respond strongly to dust evolution, with line fluxes increasing by more than an order of magnitude across the model series as the density of small dust grains decreases over time. The C2H2 line fluxes are extremely low due to a lack of abundance in the infrared line-emitting regions, despite C2H2 being commonly detected with Spitzer, suggesting that warm chemistry in the inner disk may need further investigation. Finally, we find that the CO2 flux densities increase more rapidly than the other species as the dust disk evolves. This suggests that the flux ratios of CO2 to other species may be lower in disks with less-evolved dust populations.

Tables

Table 1. Table 1: A summary of important model parameters used in each step of the procedure. Some parameters, for example the disk mass, are important for two-pop-py , MCMax , and ProDiMo . These parameters are not listed for the subsequent codes because their effects are embedded in the data passed between the codes. The only parameter in this table that is passed through to FLiTs is the disk inclination. Of the two-pop-py parameters in this table, the disk mass and gas-to-dust ratio are only valid for the initial state of the disk. The other parameters in this table are constant in time. N grains subscript 𝑁 grains N_{\mathrm{grains}} is the number of grain sizes used for two-pop-py ’s reconstruction of the grain size distribution. The dust is a distribution of hollow spheres, where the maximum fractional volume filled by the central void is 0.8 (Min et al., 2005 , 2016 ) . The dust grain mixture in MCMax is 60 % percent 60 60\% amorphous \ce ​ M ​ g 0.7 ​ F ​ e 0.3 ​ S ​ i ​ O 3 \ce 𝑀 subscript 𝑔 0.7 𝐹 subscript 𝑒 0.3 𝑆 𝑖 subscript 𝑂 3 \ce{Mg_{0.7}Fe_{0.3}SiO_{3}} silicates (Dorschner et al., 1995 ) , 15 % percent 15 15\% amorphous carbon (Zubko et al., 1996 ) , and 25 % percent 25 25\% vacuum. The viscosity helps define the surface density Σ g ​ ( r ) subscript Σ g 𝑟 \Sigma_{\mathrm{g}}(r) , where Σ g ​ ( r ) ∝ ( r / r c ) − γ visc ​ exp ⁡ [ − ( r / r c ) 2 − γ visc ] proportional-to subscript Σ g 𝑟 superscript 𝑟 subscript 𝑟 c subscript 𝛾 visc superscript 𝑟 subscript 𝑟 c 2 subscript 𝛾 visc \Sigma_{\mathrm{g}}(r)\propto\left({r}/{r_{\mathrm{c}}}\right)^{-\gamma_{\mathrm{visc}}}\exp\left[-\left({r}/{r_{\mathrm{c}}}\right)^{2-\gamma_{\mathrm{visc}}}\right] with a characteristic radius r c = 60 ​ AU subscript 𝑟 c 60 AU r_{\mathrm{c}}=60\leavevmode\nobreak\ \mathrm{AU} (Birnstiel et al., 2012 ) .

Symbol	Quantity	two-pop-py	MCMax	ProDiMo
$R_{\ast }$	Stellar radius	$2.086R_{\odot }$	$2.086R_{\odot }$
$M_{\ast }$	Stellar mass	$0.70M_{\odot }$	$0.70M_{\odot }$	$0.70M_{\odot }$
$M_{\mathrm{disk}}/M_{\ast }$	Initial disk mass	$0.1$
$R_{\mathrm{taper}}$	Taper radius	$100\mathrm{AU}$	$100\mathrm{AU}$
${\gamma }_{\mathrm{visc}}$	Viscosity exponent	$1$
$N_{\mathrm{grains}}$	Number of dust grains	$150$	$150$
$a_{\min }$	Minimum dust grain size	${10}^{-5}\mathrm{cm}$	${10}^{-5}\mathrm{cm}$
$a_{\max }$	Maximum dust grain size		$1.9\times {10}^2\mathrm{cm}$
$f_{\mathrm{vac}}$	Dust grain porosity (vacuum fraction)	$0.25$
$f_{\max }$	Maximum hollow volume ratio		$0.8$
$t_{\max }$	Final timestep	${10}^7$ years
${\alpha }_{\mathrm{turb}}$	Turbulence	${10}^{-3}$	${10}^{-3}$
$R_{\mathrm{in}}$	Inner radius	$0.06835\mathrm{AU}$	$0.07\mathrm{AU}$
$R_{\mathrm{out}}$	Outer radius	$2000\mathrm{AU}$	$600\mathrm{AU}$
$g/d$	Initial gas-to-dust ratio	$100$
$T_{\ast }$	Stellar surface temperature	$4000K$	$4000K$	$4000K$
$N_{\mathrm{rad}}$	Number of radial grid points	$240$	$500$	$240$
$N_{\theta }$	Number of azimuthal grid points		$150$
$N_Z$	Number of vertical grid points		$160$
$N_{\mathrm{inner}}$	Number of radial points near inner rim		$100$
gasevol	Gas evolution	true
$D$	Distance of the disk		$140\mathrm{AU}$	$140\mathrm{AU}$
$i$	Inclination of the disk		${45}^{\circ }$	${45}^{\circ }$
${\chi }_{\mathrm{ISM}}$	Interstellar radiation field (Draine, 1978)	$1$	$1$
$H_0$	Gas scale height at $1\mathrm{AU}$		$0.05012\mathrm{AU}$
$\beta$	Flaring exponent		$1.15$
scset	Self-consistent dust settling		true
$L_{\ast }$	Stellar luminosity			$1.0L_{\odot }$
$f_{\mathrm{PAH}}$	PAH abundance (relative to ISM)		${10}^{-4}$
$f_{\mathrm{UV}}/L_{\ast }$	UV excess			$0.01$
$p_{\mathrm{UV}}$	UV powerlaw exponent			$1.3$
${\zeta }_{\mathrm{CR}}$	Cosmic ray $\text{ce}H2$ ionization rate		$1.7\times {10}^{-17}{s}^{-1}$
$L_{X-\mathrm{ray}}$	X-ray luminosity			${10}^{30}\mathrm{erg}{s}^{-1}$

Table 2. Table 2: The emission line of each species chosen for analysis, including upper level energies E up subscript 𝐸 up E_{\mathrm{up}} and the Einstein A 𝐴 A coefficient. The description of the ro-vibrational lines of \ce ​ C ​ O ​ 2 \ce 𝐶 𝑂 2 \ce{CO2} , \ce ​ C ​ 2 ​ H ​ 2 \ce 𝐶 2 𝐻 2 \ce{C2H2} , \ce ​ H ​ C ​ N \ce 𝐻 𝐶 𝑁 \ce{HCN} , and \ce ​ N ​ H ​ 3 \ce 𝑁 𝐻 3 \ce{NH3} is an abbreviated form of that described in Jacquemart et al. ( 2003 ); Rothman et al. ( 2005 ) , where v j subscript 𝑣 𝑗 v_{j} are the normal mode vibrational quantum numbers, l j subscript 𝑙 𝑗 l_{j} are the vibrational angular momentum quantum numbers, and l 𝑙 l is the absolute value of the sum of l j subscript 𝑙 𝑗 l_{j} . The final entry, for example R ​ 11 ​ e 𝑅 11 𝑒 R11e , denotes that it is an R 𝑅 R -branch transition, the lower-state rotational energy level is 11 11 11 , and e 𝑒 e or f 𝑓 f denotes the symmetry for l 𝑙 l -type doubling.

Species	$\lambda \left(\mu m\right)$	Transition	$E_{\mathrm{up}}\left(K\right)$	$A\left(s^{-1}\right)$	Reference
$\text{ce}CO2$	$14.98299$	$v_1{v}_2{l}_2{v}_{3}r=01101\to 00001$, $Q6e$	$983.85$	$1.527$	Bosman et al. (2017)
$\text{ce}C2H2$	$13.20393$	$v_1{v}_2{v}_3{v}_4{v}_{5}l\pm =000011\to 000000$, $R11e$	$1313.1$	$3.509$	Woitke et al. (2018)
$\text{ce}HCN$	$14.03930$	$v_1{v}_2{l}_2{v}_3=0110\to 0000$, $Q6e$	$1114.1$	$2.028$	Bruderer et al. (2015)
$\text{ce}o-H2O$	$17.75408$	$J^{\prime }=6\to J^{\prime \prime }=5$	$1278.5$	$0.002869$	Notsu et al. (2017)
$\text{ce}NH3$	$10.33756$	$v_1{v}_2{v}_3{v}_4=0100\to 0000$, $J^{\prime }=3\to J^{\prime \prime }=3$	$1515.3$	$11.57$
$\text{ce}OH$	$20.11506$	$J^{\prime }=13.5\to J^{\prime \prime }=12.5$	$5527.2$	$50.47$	Woitke et al. (2018)

Full text

11institutetext: 1 Kapteyn Astronomical Institute, University of Groningen, Postbus 800, 9700 AV Groningen, The Netherlands

2 SRON Netherlands Institute for Space Research, P.O. Box 800, 9700 AV Groningen, The Netherlands

3 Anton Pannekoek Institute for Astronomy, University of Amsterdam, PO Box 94249, 1090 GE Amsterdam, The Netherlands

4 SUPA, School of Physics & Astronomy, University of St. Andrews, North Haugh, St. Andrews KY16 9SS, UK

5 Max Planck Institute for Extraterrestrial Physics, Gie\textbetaenbachstra\textbetae 1, 85741 Garching, Germany

11email: greenwood@astro.rug.nl

The effects of dust evolution on disks in the mid-IR

A.J. Greenwood1

  

I. Kamp1

  

L.B.F.M. Waters2,3

  

P. Woitke4

  

W.-F. Thi5

In this paper, we couple together the dust evolution code two-pop-py with the thermochemical disk modelling code ProDiMo. We create a series of thermochemical disk models that simulate the evolution of dust over time from $0.018\leavevmode\nobreak\ \mathrm{Myr}$ to $10\leavevmode\nobreak\ \mathrm{Myr}$, including the radial drift, growth, and settling of dust grains. We examine the effects of this dust evolution on the mid-infrared gas emission, focussing on the mid-infrared spectral lines of $\ce{C2H2}$, $\ce{CO2}$, $\ce{HCN}$, $\ce{NH3}$, $\ce{OH}$, and $\ce{H2O}$ that are readily observable with Spitzer and the upcoming E-ELT and JWST.

The addition of dust evolution acts to increase line fluxes by reducing the population of small dust grains. We find that the spectral lines of all species except $\ce{C2H2}$ respond strongly to dust evolution, with line fluxes increasing by more than an order of magnitude across the model series as the density of small dust grains decreases over time. The $\ce{C2H2}$ line fluxes are extremely low due to a lack of abundance in the infrared line-emitting regions, despite $\ce{C2H2}$ being commonly detected with Spitzer, suggesting that warm chemistry in the inner disk may need further investigation. Finally, we find that the $\ce{CO2}$ flux densities increase more rapidly than the other species as the dust disk evolves. This suggests that the flux ratios of $\ce{CO2}$ to other species may be lower in disks with less-evolved dust populations.

1 Introduction

Throughout the few-million-year lifetime of a protoplanetary disk, its evolution is dominated by the gas and dust which accrete onto the star, coalesce into larger bodies, or are dissipated by photoevaporation and stellar winds. The dust in these disks is a significant factor in their evolution. Although dust particles (assuming typical ISM abundances) make up only a small portion of the total disk mass, small grains composed of dust and ice dominate the radiative transfer of the disk and provide “seeds” for planets to form.

It is crucial that we can determine the distribution of this dust. In the eras of ALMA, JWST, and E-ELT we have both the sensitivity and the angular resolution to spatially resolve nearby protoplanetary disks, to spectrally resolve the emission lines from molecules, and to detect the broad emission features from silicate grains. The distributions of gas and dust in the disk are not necessarily co-spatial, and

nor is the ratio constant with height, because dust grains are expected to grow and fall towards the midplane (Blum & Wurm, 2008; Dominik & Dullemond, 2008) and the scale heights of the dust and gas disks may differ (de Gregorio-Monsalvo et al., 2013a; Avenhaus et al., 2018). T Tauri stars and brown dwarfs have flatter dust disks and relatively thicker gas disks than Herbig stars (Mulders & Dominik, 2012). A crucial step that we must take is to understand the physics behind these effects, and how we can incorporate them into a thermochemical disk model.

In the infrared, we can observe the spectral lines of many species that trace the warm upper layers of the inner disk. For example, typical mid-infrared water line fluxes in T Tauri disks have been observed to range between about $2\times 10^{-13}\leavevmode\nobreak\ \mathrm{erg}\leavevmode\nobreak\ \mathrm{cm}^{-2}\leavevmode\nobreak\ \mathrm{s}^{-1}$ and $2\times 10^{-15}\leavevmode\nobreak\ \mathrm{erg}\leavevmode\nobreak\ \mathrm{cm}^{-2}\leavevmode\nobreak\ \mathrm{s}^{-1}$ (Pontoppidan et al., 2010). The authors also note that there appears to be a sub-class of disks which have strong $14.98\leavevmode\nobreak\ \mu\mathrm{m}$ $\ce{CO2}$ emission, but with no detectable contributions in the $14$ to $16\leavevmode\nobreak\ \mu\mathrm{m}$ range from water, $\ce{HCN}$, and $\ce{C2H2}$. Previous results by Meijerink et al. (2009) and Bosman et al. (2017) suggest that a high gas-to-dust ratio of $1000{:}1$ in the line-emitting regions might be necessary to reproduce the mid-infrared line fluxes that we have observed with Spitzer. This ratio has often been adopted by the modelling community (Bruderer et al., 2015; Bosman et al., 2017).

To justify this ratio, we need to explain what has happened to the dust: if the disk formed from a molecular cloud with a gas-to-dust ratio of $100{:}1$, this dust cannot disappear without reason. We know that dust settles in disks (Dubrulle et al., 1995), and that especially in older disks, the midplanes are typically more dust-rich than the upper layers (Dullemond & Dominik, 2005). We need to account for this dust settling in order to quantify the gas-to-dust ratio of the line-emitting regions. In this paper, we couple a dust evolution code with a thermochemical modelling code and for the first time analyze the effects of dust evolution on the mid-infrared spectral lines.

Accurately determining the gas and dust masses of protoplanetary disks is difficult. Primarily, because much of the line and continuum emission is optically thick: there is no known way to derive gas and dust masses from line and continuum observations that is free of significant degeneracies. Although observing disks with gaps can help to measure scale heights and thus estimate midplane densities, there can be a lot of optically-thick gas and dust in the midplane that cannot be observed directly – for example, Pinte et al. (2016) find that the dust density in the gaps of HL Tau is at least a factor of 10 lower than in the rings surrounding the gaps. However, it is difficult to place an upper limit on the dust density in the rings.

Likewise, it is difficult to accurately observe the dust grain size distribution as a function of radius and height. For example, although there is progress to spatially resolve the dust surface density of disks (Pinilla et al., 2014), the measured spectral index is dependant upon assumptions about the dust temperature. Scattered light imaging using SPHERE has been used to observe the surface of dust disks (Ginski et al., 2016; Stolker et al., 2016; Avenhaus et al., 2018; Muro-Arena et al., 2018). Particularly, the survey of eight T Tauri disks introduced by Avenhaus et al. (2018) is a valuable set of observations. Although a detailed analysis of each source is left for forthcoming papers, Avenhaus et al. (2018) provide insights such as showing that all of the observed T Tauri disks appear to have $\tau=1$ surface flaring indices between $\sim 1.1$ and $\sim 1.6$, and shows striking geometric structures such as gaps, rings, and even the lower disk surface of three disks. If an inclined disk has rings or gaps in the dust structure, the height of the dust scattering surface can be estimated by geometrically fitting these structures (Ginski et al., 2016; Avenhaus et al., 2018).

In comparison to micron-sized dust grains, decimetre-sized bodies have very little effect on the radiative transfer of the disk in the mid-infrared regions. However, their effect on the disk’s evolution is much more significant. Because they are less affected by processes such as viscosity and turbulence, decimetre-sized bodies are more likely to be found close to the disk’s midplane, and are a crucial step towards the eventual formation of planets.

To analyze the effect of dust evolution and migration on the mid-infrared spectral lines of a T Tauri protoplanetary disk, we model the evolution of dust over time, and create a thermochemical model of the disk at discrete time steps in the evolution. In this way, we incorporate a self-consistent description of dust evolution into our thermochemical disk models, where the gas-to-dust ratio and dust size distribution vary not only at each timestep, but also radially and vertically through the disk. By analyzing the mid-infrared spectral lines of each individual model, we can build an understanding of how dust evolution affects these lines.

2 Modelling strategy

A total of four separate codes were used in a set order to produce the models in this paper. The underlying disk model is based upon the “DIANA standard” T Tauri disk described in Woitke et al. (2016), except for the dust structure. Table 1 describes the parameters of our T Tauri disk model. This section serves as a technical description of what each code does and how they are interfaced.

First, we use the two-population dust evolution code two-pop-py111two-pop-py is available on Github at http://birnstiel.github.io/two-pop-py/. (Birnstiel et al., 2012, 2015, 2017). The two-pop-py code produces a one-dimensional description of the dust at each age increment. It is a simplified version of the dust evolution code described in Birnstiel et al. (2010), calibrated with the use of these more complex models and based upon the use of two dust grain populations: a small monomer population, and a large grain population. The small grain population is coupled to the gas, not affected by drift velocities, and is constant in time and space. The population of large grains can grow in size, and is affected by the limiting mechanisms of radial drift and grain fragmentation. The reconstruction of the grain size distribution from the simple two-population model is described by Birnstiel et al. (2015).

From two-pop-py, we obtain the gas and dust surface densities as functions of radius, and the dust mass fraction as a function of both radius and grain size (which records the fraction of the total mass at each radius that is accounted for by each grain size). For each age increment, we calculate a separate disk model. For this, we use the code MCMax (Min et al., 2009) to create a two-dimensional model of the disk structure and calculate the dust continuum radiative transfer for the dust temperature. We ensure that the grain sizes are conserved between two-pop-py and MCMax.

There are two standard methods in MCMax to calculating the scale heights of the gas and dust: to have parameterized scale heights, or to solve for hydrostatic equilibrium. We take a hybrid approach, first forcing the scale heights of the gas to remain parameterized, then solving the dust for hydrostatic equilibrium and allowing the dust to settle in a self-consistent manner (described by Mulders & Dominik 2012). We do not argue for the physical correctness of this procedure (which is unknown), rather for the simplicity it brings. The first argument is that we want to remain as close to the DIANA standard of models as possible, where only the dust distribution deviates from the standard. The second argument is that it simplifies the results and analysis: although the total gas mass in the disk decreases with age, the gas scale heights remain constant. This allows us better to analyse the effects that different dust distributions have on the mid-infrared molecular lines, under the assumption that everything else is equal.

The third code in the pipeline is ProDiMo (Woitke et al., 2009; Kamp et al., 2010; Aresu et al., 2011), which takes the MCMax output and calculates a full thermochemical disk model. We take the disk structure directly from MCMax, then run the disk chemistry, taking the reaction rates for 235 species from the UMIST2012 chemical network (McElroy et al., 2013) to calculate the disk’s thermal and chemical equilibrium state. Because the overall opacity and density structure of the disk is fixed, there is no need for global disk structure iterations.

The final element is FLiTs, which calculates a high-resolution spectrum from a ProDiMo model (Woitke et al., 2018) across the infrared wavelength range. ProDiMo provides the gas and dust temperatures, the stellar spectrum and local radiation field, opacities, molecular number densities, and level populations. FLiTs reads this file and calculates the disk’s spectrum for the desired molecules and wavelength range.

3 Dust migration and disk surface densities

The effects of dust migration over time are visible in both the surface densities and dust grain size distributions, as computed by two-pop-py. These are one-dimensional results, while we later use MCMax to compute the two-dimensional disk structure. The precise timescales are not strictly relevant for this discussion because the ages are very uncertain, both in the models and any observed disks. It is better to compare disks with different surface densities than to compare disks of different ages, because these parameters are more readily observable. For example, disks in Ophiuchus are thought to be less than $1\leavevmode\nobreak\ \mathrm{Myr}$ old (Furlan et al., 2009), but basic properties such as the disk mass may vary by about two orders of magnitude (Testi et al., 2016). We need many different models in order to accurately represent such a diverse disk population. Nevertheless, for brevity we refer to each model by its two-pop-py age.

In each ProDiMo model we use steady-state chemistry: the ProDiMo chemistry may run for longer than the “age” of the two-pop-py model. However, the timescales for disk surface chemistry (on the order of years) are much shorter than the timescale of dust evolution ($\gtrsim 10^{4}\leavevmode\nobreak\ \mathrm{years}$): it is only when considering longer wavelengths such as the sub-mm that we might need to calculate time-dependent chemistry.

Figure 2 show the gas and dust surface densities from $t_{0}=0.018\leavevmode\nobreak\ \mathrm{Myr}$ to $t_{\mathrm{max}}=10\leavevmode\nobreak\ \mathrm{Myr}$. Although the gas mass does slowly decrease over time due to viscous evolution, the loss of dust due to radial drift is much greater222To simplify the analysis of comparing different dust structures, we do not include the accretion luminosity in our MCMax and ProDiMo models., and this causes the total gas-to-dust ratio of the disk to increase over time, from $100$ at $t=0$ to over $6000$ at $10\leavevmode\nobreak\ \mathrm{Myr}$. For the column of gas at a radius of $1\leavevmode\nobreak\ \mathrm{AU}$, the youngest model has a gas-to-dust ratio of about $40$, whereas the oldest model reaches about $4\times 10^{4}$.

The youngest two-pop-py model is close to the initial state of the system, and the low gas-to-dust ratio in the initial state of the inner disk is the result of a high level of inwards dust mass flux. As the disk ages, radial drift continues to deplete the total dust mass in the disk. One other notable trend in Fig. 2 is the “bump” in gas-to-dust ratio, which becomes significant at $0.32\leavevmode\nobreak\ \mathrm{Myr}$ and a radius of $40\leavevmode\nobreak\ \mathrm{AU}$, and moves inwards with time. This bump is due to a localized depletion in the dust, reported by Birnstiel et al. (2012) as a pile-up effect caused by increases in the disk temperature due to viscous heating, similar to the grain pile-up in Youdin & Chiang (2004).

Figure 1 describes how the gas mass, dust mass, and gas-to-dust ratio vary over time, integrated over the entire disk. Figure 3 shows vertical cuts of the grain size distribution at radii of $0.1$, $1$, $10$, and $100\leavevmode\nobreak\ \mathrm{AU}$. At radii of $10\leavevmode\nobreak\ \mathrm{AU}$ or less, the trends are similar: the grain size distribution shows both an overall decrease in surface densities, and a disproportionate depletion of larger grains (which are transported inwards due to radial drift). As a consequence, in the inner disk at $0.1\leavevmode\nobreak\ \mathrm{AU}$ there are many large grains up to $10\leavevmode\nobreak\ \mathrm{cm}$ in size at an age of $10\leavevmode\nobreak\ \mathrm{Myr}$. However, at $1\leavevmode\nobreak\ \mathrm{AU}$ this cut-off is at a grain size of $1\leavevmode\nobreak\ \mathrm{cm}$, and about $0.05\leavevmode\nobreak\ \mathrm{cm}$ at $10\leavevmode\nobreak\ \mathrm{AU}$. This happens because a population of large grains can only be maintained so long as there is an inwards flux of grains coming from larger radii. The situation is slightly different at $100\leavevmode\nobreak\ \mathrm{AU}$: although there is no significant population of cm-sized, grains, the population of mm-sized grains grows until an age of $0.56\leavevmode\nobreak\ \mathrm{Myr}$.

4 Size of the dust disk

One consequence of our dust settling method is that the dust disk is less vertically-extended than the gas disk. This is not an artefact of parameterizing the gas structure, but a result of the self-consistent description of dust settling. Disks with processed dust and dust disks that are thinner than the gas disk have been reported in previous literature (de Gregorio-Monsalvo et al., 2013a; Ginski et al., 2016; Muro-Arena et al., 2018; Avenhaus et al., 2018), but the number of observations is still relatively few.

The dust disk calculated using our self-consistent settling is thinner than a disk model that uses Dubrulle settling (Dubrulle et al., 1995). The main differences between the self-consistent settling method (Mulders & Dominik, 2012) and Dubrulle setting is that Dubrulle settling uses only the midplane density and temperature. Away from the midplane, their analytical solution becomes less accurate. Self-consistent settling agrees well with Dubrulle settling in the midplane, but sub-micron particles can also decouple from the gas in the upper disk (Mulders & Dominik, 2012), contrary to Dubrulle settling. The distribution of $<0.1\leavevmode\nobreak\ \mu\mathrm{m}$ particles is nearly uniform throughout a Dubrulle-settled disk, and larger grain sizes are affected less by settling than the self-consistent model (Mulders & Dominik, 2012).

All of our models have a dust disk that is thinner than the gas disk. Figure 3 shows vertical cuts of the gas-to-dust ratio through the disks at each age: in every case, in the upper layers of the disk where some tenuous gas still exists, the gas-to-dust ratio can easily be $10^{5}$ or more. In the midplane of the inner disk, particularly at $0.1\leavevmode\nobreak\ \mathrm{AU}$, the gas-to-dust ratio of models younger than $0.56\leavevmode\nobreak\ \mathrm{Myr}$ can be less than $10$. This is in contrast to at $100\leavevmode\nobreak\ \mathrm{AU}$, where most of the disk at $0.018\leavevmode\nobreak\ \mathrm{Myr}$ has a gas-to-dust ratio of $100$. At young ages, there are few large grains in the outer disk: the dust only decouples significantly from the gas in older models, once mm-sized dust grains have formed. However, at small radii, there is already a significant population of large grains at $t=0.18\leavevmode\nobreak\ \mathrm{Myr}$. This gives rise to the low gas-to-dust ratio in the midplane at small radii.

When comparing the gas and dust disk thickness with observations, our youngest models are comparable to disks such as IM Lup that appear to be less settled (Cleeves et al., 2016), while the older models are more comparable to disks such as HD 163296, where the dust disk may be significantly settled in comparison to the gas disk (de Gregorio-Monsalvo et al., 2013b).

The settled dust structure causes every disk to have some regions where the gas-to-dust ratio is $1000$ or more, which as we find in Sect. 5 is the minimum gas-to-dust ratio from which all species except $\ce{C2H2}$ tend to emit.

Both dust evolution and the new self-consistent method of settling contribute to changes in the mid-infrared line fluxes, because both of these processes affect small micron-sized dust grains. Dust evolution affects the surface density of these grains in the disk, while settling affects the vertical distribution of those grains. Figures 4 and 5 shows the effects of dust evolution and the settling method on dust densities, absorption coefficient, and gas temperatures, while Fig. 6 shows those same effects on the mid-infrared line fluxes.

Counter-intuitively, the mid-infrared line fluxes can in fact increase due to Dubrulle settling. Even though Dubrulle settling leads to lower gas temperatures in the mid-infrared line-emitting region for each species, in Fig. 6 we do see that models with Dubrulle settling tend to have somewhat stronger lines than those same models with self-consistent settling. Although the total vertical optical depth remains the same in models with both Dubrulle and self-consistent settling, the dust in the latter models is more concentrated towards the midplane. Thus, reduced optical depths in the upper layers of models with Dubrulle settling can allow the mid-infrared line-emitting regions to grow in size. However, it is important to note that although the method of settling can affect the line fluxes, the effects of dust evolution remain more significant.

5 Mid-infrared spectra

The migration and evolution of dust has a profound effect on the mid-infrared spectra and line-emitting regions of our disk model. To discuss the line-emitting regions of each species, it is necessary to define the exact molecular lines we are using, and how the line-emitting regions themselves are derived.

We define the line-emitting region as the area from which $70\%$ of the flux originates, in both the radial and vertical directions. The lower limit, $\mathrm{x15}$, is defined so that $85\%$ of the total line flux of that spectral line is emitted at radii greater than $\mathrm{x15}$. The upper limit, $\mathrm{x85}$, is the opposite – only $15\%$ of the total line flux of that spectral line is emitted at a radius greater than $\mathrm{x85}$. At each radius, $\mathrm{z15}$ and $\mathrm{z85}$ define the heights in the disk above which $15\%$ and $85\%$ of the flux is emitted respectively. The properites of the line-emitting region for a given species (such as $T_{\mathrm{gas}}$) are averaged and weighted with the volume density of that species across the line-emitting region.

Table 2 details the molecular lines chosen for analysis, one for each species. For consistency, where possible, we analyze spectroscopic lines that have previously been analyzed in other literature (see references in Tab. 2). Whenever an individual molecular line is referenced in this paper, it refers to the line in this table.

The flux of a particular molecular line is driven primarily by three factors: the gas-to-dust ratio, the gas temperature in the line-emitting region, and the optical depth in the line-emitting region $\langle\tau_{20\leavevmode\nobreak\ \mu\mathrm{m}}\rangle$. Figure 7 shows how the average gas-to-dust ratio $\langle\mathrm{g}/\mathrm{d}\rangle$, average gas temperature $\langle T_{\mathrm{gas}}\rangle$, and the line fluxes change with age in the disk. Although the gas-to-dust ratio of the entire disk begins at $100{:}1$, the self-consistent dust settling is not informed by the “age” of the model. It solves for the vertical dust structure without considering the settling timescales, which may be a few $\times 10^{5}$ years (Dullemond & Dominik, 2004).

We measure fluxes in two different ways. Where we mention the flux of an individual spectral line, these fluxes are calculated using a vertical escape probability method (Woitke et al., 2009). These line fluxes are accurate only for a face-on disk, because no detailed radiative transfer is taken into account. In contrast, where discuss FLiTs spectra, these have been calculated for an inclined disk and account for both the radial and vertical optical depth. They are also convolved to an instrumental resolution, so that we measure the peak flux of a complex of lines and not the integrated flux of a single line. The purpose of the first method is to measure how line emission responds to changes in the dust distribution, by picking these representative lines to measure the properties of the line-emitting region. The purpose of the second method is to put these measurements into an observational context. By accounting for the many spectral lines that are in the bandheads of some species and convolving the resulting spectrum, we can relate the peak fluxes of the spectra to observations and sensitivities of Spitzer or JWST.

In comparison to observed Spitzer spectra (e.g. Pontoppidan et al. 2010; Salyk et al. 2011; Pascucci et al. 2013), the FLiTs spectra of all of our species except $\ce{C2H2}$ can produce spectral lines that are at least as bright as we have observed. The spectra are highly responsive to changes in the dust distribution, where the loss of small dust particles over time acts to reduce the opacity of the disk and dramatically increases line fluxes. The only exception to this trend is $\ce{C2H2}$, where our models do not produce lines that are as bright as Spitzer observations (despite $\ce{C2H2}$ being readily detected), suggesting that the chemical networks of $\ce{C2H2}$ may be inaccurate or incomplete.

5.1 Line fluxes: C2H2 in comparison to other species

Throughout the dust evolution simulation, the line flux of each species except $\ce{C2H2}$ increases by at least an order of magnitude, whereas $\ce{C2H2}$ remains more or less constant in line flux. For all species except $\ce{C2H2}$, the gas-to-dust ratio in the line-emitting region increases significantly as the dust evolves.

On the other hand $\ce{C2H2}$ simply does not exist in these models in the upper layers of the disk. It is easily destroyed when the dust evolves and the disk becomes more optically thin. As the disk ages the remaining $\ce{C2H2}$ is concentrated further towards the midplane, and the line-emitting region is pushed both towards the midplane and out to larger radii (up to $\sim 0.6\leavevmode\nobreak\ \mathrm{AU}$). In these models, we see the same bifurcated structure of $\ce{C2H2}$ discussed in Greenwood et al. (submitted) and seen also in Walsh et al. (2015). A significant proportion of the line-emitting population of $\ce{C2H2}$ is below the $\mathrm{A}_{V}=10$ line, meaning that the dust continuum in this region is very optically thick in the mid-infrared. As the dust evolves and becomes more optically thin, the $\ce{C2H2}$ responds and shifts closer towards the midplane. Hence, there is no chance for the line flux to increase. Greenwood et al. (submitted) show that a more dust-rich disk model has stronger $\ce{C2H2}$ lines, where for all other species, more dust-rich models have weaker lines. Our models from this work do not predict $\ce{C2H2}$ line fluxes that are observable with any current spacecraft. However, there are clear Spitzer detections of $\ce{C2H2}$ in T Tauri disks with similar detection rates to the other species (Pontoppidan et al., 2010). More work is needed to understand this particular species.

5.2 Properties of the line-emitting regions

Figure 7 shows changes in the gas-to-dust ratio, gas temperature, vertical dust optical depth, escape probability line flux, and line-emitting area of each line over time. Although for each species the gas-to-dust ratio in the line-emitting region increases steadily over time, $\ce{C2H2}$ has a consistently lower ratio than the other species because its line-emitting area is closer towards the mid-plane.

The trends in gas temperature are more complex. For $\ce{C2H2}$, at $1.8\leavevmode\nobreak\ \mathrm{Myr}$ the gas temperature of the line-emitting region drops dramatically, from about $900\leavevmode\nobreak\ \mathrm{K}$ to $350\leavevmode\nobreak\ \mathrm{K}$. The temperature drops even more as the disk ages further. This trend is also reflected in the fact that $\ce{C2H2}$ is in absorption for the oldest few disk models (see Fig. 9). The reason for this is changes in the line-emitting area: the $\ce{C2H2}$ abundances in the oldest disk models drop significantly, and the line-emitting area moves outwards and towards the mid-plane. The gas temperatures of the $\ce{CO2}$ and $\ce{NH3}$ line-emitting regions stay relatively stable over time. For $\ce{HCN}$ and $\ce{H2O}$, we see an increase in gas temperature, particularly for the oldest models. This is because the line-emitting regions become larger and move towards slightly warmer, higher layers in the disk.

Figure 7 shows that where we compute the spectra of our models using FLiTs, we expect the fluxes of every species except $\ce{C2H2}$ to increase significantly as the dust in the disk evolves. This is illustrated by Figs. 8 and 9, which show that between the youngest and oldest disks models, the peak flux densities for every species except $\ce{C2H2}$ increase by factors ranging from 30 for $\ce{H2O}$ to 200 for $\ce{OH}$. The escape probability line fluxes in Fig. 7 increase by similar amounts, showing that measuring the escape probability line flux of these single lines forms a reasonable proxy for measuring trends in the overall strength of the line-emitting region.

For comparison to infrared observatories, the median $1\leavevmode\nobreak\ \sigma$ mid-infrared line sensitivity of Spitzer is about $8.5\times 10^{-19}\leavevmode\nobreak\ \mathrm{W\leavevmode\nobreak\ m}^{-2}$, thus many individual spectral lines are far too faint for Spitzer to detect.333irsa.ipac.caltech.edu/data/SPITZER/docs/irs/. However, the ro-vibrational bandheads of species such as $\ce{CO2}$ contain many lines that at a low spectral resolution are blended together, which is why we can see these species with Spitzer. In some of the more evolved disks, even individual $\ce{OH}$, $\ce{o-H2O}$, $\ce{NH3}$, and $\ce{HCN}$ lines are significantly above this sensitivity limit (see Fig. 7). Following the flux densities of the convolved spectra presented in Figs. 8 and 9, and assuming a $5\leavevmode\nobreak\ \sigma$ Spitzer sensitivity of $5\leavevmode\nobreak\ \mathrm{mJy}$ in $512\leavevmode\nobreak\ \mathrm{s}$, we expect that only $\ce{C2H2}$ would be undetectable at every age. $\ce{OH}$, $\ce{H2O}$, and $\ce{NH3}$ are only detectable in more evolved disks, while $\ce{HCN}$ and $\ce{CO2}$ should be detectable in every model. At $R=600$, the complexes of water lines can easily become blurred at Spitzer’s low spectral resolution, making it difficult to determine the continuum flux (in comparison to $R=2800$ spectra). The $10\leavevmode\nobreak\ \sigma$ continuum sensitivity of JWST’s MIRI instrument is around $0.1\leavevmode\nobreak\ \mathrm{mJy}$ at $15\leavevmode\nobreak\ \mu\mathrm{m}$, with an exposure time of $10^{4}\leavevmode\nobreak\ \mathrm{s}$ (Glasse et al., 2015). Thus with JWST, we can expect that in every model, all species except $\ce{C2H2}$ would be detectable with exposure times much less than $10^{4}\leavevmode\nobreak\ \mathrm{s}$.

$\ce{C2H2}$ even goes into absorption at $5.6\leavevmode\nobreak\ \mathrm{Myr}$ and $10\leavevmode\nobreak\ \mathrm{Myr}$. The line emission features in disk models up to $1.8\leavevmode\nobreak\ \mathrm{Myr}$ old appear to be dominated by emission from near the inner wall, while in the older disk models the line-emitting region extends out to $0.6\leavevmode\nobreak\ \mathrm{AU}$ and is close to the mid-plane. The absorption lines in the $5.6\leavevmode\nobreak\ \mathrm{Myr}$ and $10\leavevmode\nobreak\ \mathrm{Myr}$ old models are likely attributable to the fact that the $\ce{C2H2}$ emission comes from close to the mid-plane. As described in Greenwood et al. (submitted), there exists an inversion in the gas temperature around the $\mathrm{A}_{V}=1$ line, such that there may exist clouds of colder gas above the warmer gas in the disk. Absorption lines may result if these clouds of colder gas exist between our line of sight and the line-emitting region of the species.

$\ce{CO2}$ is the brightest species by a significant margin, suggesting that our younger models may be representative of the sub-class discovered by Pontoppidan et al. (2010), where in six disks only $\ce{CO2}$ was unambiguously detected. This sub-class may consist of disks where the dust is significantly settled in comparison to the gas, similar to the structure produced in our models by the self-similar dust settling. In this case, we would expect $\ce{CO2}$ to be the brightest species in the mid-infrared, and thus other species may remain undetected. The reason that $\ce{CO2}$ fluxes can be so high is that it is a robust species which is able to survive in relatively optically thin regions because it has a small cross-section for photodissociation (van Dishoeck et al., 2006). As the dust structure evolves and becomes more optically thin, the line-emitting area of $\ce{CO2}$ also grows: Fig. 7 shows that $\ce{CO2}$ can have a line-emitting area close to $1000\leavevmode\nobreak\ \mathrm{AU}^{2}$ (the $\ce{CO2}$ emission extends to about $15\leavevmode\nobreak\ \mathrm{AU}$ at $1\leavevmode\nobreak\ \mathrm{Myr}$).

6 Discussion

As the dust evolves, grains are depleted from the disk and the mid-IR molecular line emission grows substantially for every species except $\ce{C2H2}$. This happens because dust particles (particularly small micron-sized grains) are the main carriers of opacity in the disk, so reduced continuum optical depths lead directly to a larger line-emitting area and increased molecular line fluxes.

The lack of $\ce{C2H2}$ emission is because this species exists closer towards the mid-plane than the other molecules. Thus, it is more affected by both the optical depth of the dust in the disk and the gas temperature inversions around $\mathrm{A}_{V}=1$. The T Tauri spectra observed by Carr & Najita (2011) have peak $\ce{C2H2}$ flux densities of up to about $20\leavevmode\nobreak\ \mathrm{mJy}$, and $\ce{C2H2}$/$\ce{HCN}$ line flux ratios between about $0.08$ and $0.8$. As Figs. 7, 8 and 9 show, when we compare the peak flux densities of the FLiTs spectra of our strongly-evolved models, we see much lower ratios: although the ratio between the spectral flux densities of $\ce{C2H2}$ and $\ce{HCN}$ is $\sim 0.02$ at $0.18\leavevmode\nobreak\ \mathrm{Myr}$, it decreases with age to as low as $5\times 10^{-5}$ at $10\leavevmode\nobreak\ \mathrm{Myr}$ (with the caveat that for the oldest models, $\ce{C2H2}$ is in weak absorption). The fact that $\ce{C2H2}$ behaves so differently to other species suggests either that our chemical network misses some key $\ce{C2H2}$ formation pathways, or that the dust structure of our models is not representative of disks with $\ce{C2H2}$ detections. Although it is possible that uncertainties in the chemical network and reaction rates are leading to significant inaccuracies in the distribution of $\ce{C2H2}$ calculated by ProDiMo, our $\ce{C2H2}$-deficient models may also be a natural result of the dust structure. Walsh et al. (2015) find that $\ce{C2H2}$ is significantly less abundant in Herbig Ae disks due to increased levels of photodissociation, which also increases when the population of small dust grains is reduced. This suggests that the low $\ce{C2H2}$ line fluxes are a result of an evolved dust structure, but it does not adequately explain the fact that $\ce{C2H2}$ is all but absent in even our youngest $0.018\leavevmode\nobreak\ \mathrm{Myr}$ model.

Another notable result is that $\ce{CO2}$ is much brighter than all other species in our models, suggesting that our older models may be representative of the disks found by Pontoppidan et al. (2010) where only $\ce{CO2}$ was detected: if these sources have an evolved and settled dust structure, our results show that we would expect $\ce{CO2}$ to be much brighter than other species. The sources where $\ce{CO2}$ was not the only detected species may have a less evolved and settled dust structure, giving lower $\ce{CO2}$ line fluxes and more equal ratios between the line fluxes of different species. In all of our models, the $\ce{CO2}$ spectrum has a brighter flux density than all other species by at least a factor of ten. The explanation for the relative brightness of $\ce{CO2}$ likely lies in the dust and temperature structure of the models.

While Pontoppidan et al. (2010) note a strong correlation between spectral type and line detections, there also appears to be a dichotomy in their results between B, A, and F type objects and G, K, and M spectral types. For spectral types earlier than G, there is a dramatic and as-yet unexplained drop in the detection rates of all species except $\ce{CO}$. This may appear to be in conflict with Mulders & Dominik (2012), who suggest that late-type (brown dwarf and T Tauri) disks have flatter dust disks and relatively thicker gas disks than earlier (Herbig) disks do.444One caveat to their analysis is that the best-fit brown dwarf, T Tauri, and Herbig disk models have identical outer radii. Their brown dwarf model has a very low surface density and is strongly flared, while the outer radius of $400\leavevmode\nobreak\ \mathrm{AU}$ is not representative of a typical brown dwarf disk (Bate et al., 2003). In this paper we suggest that a highly-evolved and strongly-settled dust structure may lead to $\ce{CO2}$ emission being much brighter than other species. However, if younger disks have less-evolved dust (thus more small dust grains) and a less-settled dust structure, and because our models show that $\ce{CO2}$ line fluxes increase with age more than other species, we suggest that the flux ratios between $\ce{CO2}$ and the other species that we have modelled will be lower in early-type disks.

The work by Mulders & Dominik (2012) is based upon SEDs, and they state that “regions with the same temperature have a self-similar vertical structure independent of stellar mass”. It is difficult to distinguish changes in the SEDs of our dust evolution models at wavelengths below $30\leavevmode\nobreak\ \mu\mathrm{m}$; thus we are probing different parameter spaces. The lack of detections around early-type stars could be because the stronger UV radiation of early-type stars is enough to prevent species that are vulnerable to photodissociation from forming in sufficient quantities. Such species include $\ce{C2H2}$, $\ce{H2O}$, $\ce{NH3}$, and $\ce{HCN}$, which all have photodissociation cross sections $>10^{-17}\leavevmode\nobreak\ \mathrm{cm}^{-2}$ (van Dishoeck et al., 2006). Walsh et al. (2015) find that weaker UV fields do indeed allow for more molecule-rich disk atmospheres. Thermochemical models by Antonellini et al. (2015) have found that strong UV radiation does not weaken water line fluxes and that likely explanation for the lack of water detections around Herbig stars is noisy spectra combined with high continuum flux levels. Future modelling and observational work is required in order thoroughly to understand the interplay between the disk dust structure, the radiation field, and the mid-infrared molecular lines. What appears most clear is that as the main opacity carrier in the disk, the evolution and settling of dust can have a substantial impact.

7 Conclusion

We have coupled the dust evolution code two-pop-py with the radiative transfer code MCMax, the thermochemical disk modelling code ProDiMo and the line-tracing code FLiTs. We have produced a series of models that simulate the evolution of dust over time in a $2$D thermochemical disk model, and have created infrared spectra of these models including $\ce{C2H2}$, $\ce{CO2}$, $\ce{HCN}$, $\ce{NH3}$, $\ce{OH}$, and $\ce{H2O}$ mid-infrared line emission.

This paper shows that dust evolution has a very clear and straightforward effect on the mid-infrared line fluxes, where the evolving and diminishing dust population in the surface acts to decrease continuum optical depths and thus increase line fluxes by up to a factor of $100$. Even in the least-evolved disk models, the mid-infrared lines prefer to emit from a gas-to-dust ratio of at least $1000{:}1$. This provides a physical explanation for the need in previous literature to use such ratios in order to produce enough line flux to match models with observations.

This is a double-edged sword: although the dust structure is a significant degeneracy when fitting models to observations, if we can observe the dust directly (for example, by observing the $\tau=1$ surface with scattered-light imaging), the dust structure also becomes a powerful diagnostic in understanding the mid-infrared emission of T Tauri disks. Upcoming observatories such as JWST and E-ELT will significantly improve our ability to observe these disks, and combined with more sophisticated models, we can build a two- or three-dimensional understanding of dust structure in the inner disk and the mid-infrared spectral lines.

8 Acknowledgements

We would like to thank Lucia Klarmann, Gabriela Muro Arena, Paola Pinilla, Til Birnstiel, and Michiel Min for their discussions and comments on the manuscript. We also thank the Center for Information Technology of the University of Groningen for their support and for providing access to the Peregrine high performance computing cluster.