© ESO, 2012

1. Introduction

Protoplanetary disks are the main sites of planet formation. Within them, the small sub-micron sized dust grains grow to millimeter and centimeter sizes and settle to the midplane, where they eventually form kilometer-sized planetesimals and proto-planets that are the building blocks of planetary systems. How the dust grows and settles depends not only on the grain size and gas density, but also on the turbulent mixing strength of the gas that mixes small grains back up into the disk atmosphere and makes them collide in the midplane. Without turbulent mixing, all disks would be flat and dust would grow slower. Additionally, turbulent mixing is an important driver for the viscous evolution of the disk, for radial mixing of dust, and it may also play an important role in the later stages of planet formation.

Although turbulent mixing is fundamental for our understanding of disk evolution and planet formation, the cause of disk turbulence is not completely understood (e.g. Balbus & Hawley 1998; Armitage 2011). Turbulent mixing in disks is often implemented using the α prescription from Shakura & Sunyaev (1973) and Pringle (1981), but the value of the turbulent mixing strength α_turb is hard to predict from theory. Although it could be empirically derived from line-broadening with future observatories such as ALMA, this is challenging with current facilities (Hughes et al. 2011). The value typically used and inferred from disk models trying to match accretion rates is α_viscous = 0.01 (e.g. Hartmann et al. 1998). Comparison of steady-state accretion disk models to millimeter images yields α_viscous = 0.5...10^-4 (Isella et al. 2009).

The turbulent mixing strength manifests itself in the degree of dust settling. When dust grains grow to larger sizes, they decouple from the gas and start settling toward the midplane. At what size and height they decouple depends not only on the gas density and temperature but also on the turbulent mixing strength: a lower turbulent mixing strength leads to flatter disks (e.g. Dullemond & Dominik 2004b). As dust decouples, its relative velocity also increases, which leads to fragmentation, providing a natural stop to grain growth and replenishing the small dust grain population that is then mixed up into the disk surface (Dullemond & Dominik 2005).

There is by now a wealth of observational evidence supporting grain growth and dust settling in protoplanetary disks. Most notably, mid-infrared spectroscopic observations have shown that grains in the warm surface layers have grown beyond the size of a micron (e.g. Henning & Meeus 2009). In addition, the spectral index at (sub)millimeter and centimeter wavelengths shows that grains grow up to centimeter sizes in the cold disk midplane (e.g. Beckwith & Sargent 1991; Mannings & Emerson 1994). This picture has been confirmed with hydrostatic radiative transfer models by D’Alessio et al. (2006) and Furlan et al. (2005) for T Tauri stars, who require a depletion of small dust grains at the disk surface and a population of big grains in the disk midplane. Without grain growth and dust settling, hydrostatic disk models overpredict observed far-infrared fluxes.

Interestingly enough, this seems not to be the case for Herbig Ae and Be stars, the intermediate-mass counterparts of T Tauri stars. Hydrostatic disk models do well in explaining far-infrared fluxes without the need to deplete the upper layers of the disk (Dominik et al. 2003; Dullemond & Dominik 2004a; Acke et al. 2009), especially in Meeus group I sources¹. However, the (sub) millimeter slopes show that a population of millimeter-sized grains exist in these disks as well (e.g. Acke et al. 2004), indicating that the grain-size population is not that different from T Tauri stars. This opens up the possibility that these disks are less settled, and that disk properties – in particular the turbulent mixing strength – could vary strongly as a function of stellar mass. The recent discovery that disks around brown dwarfs appear to be more settled than T Tauri stars fits well into this picture (e.g. Lada et al. 2006; Szűcs et al. 2010), although this is not always found (Furlan et al. 2011), and indicates that the observed degree of settling could be inversely proportional to stellar mass.

However, the degree of settling does not directly reflect the turbulent mixing strength. Disks around T Tauri stars, Herbig stars and brown dwarfs vary in disk mass (e.g. Beckwith et al. 1990; Scholz et al. 2006; Vorobyov 2009) and stellar mass, and observations at the same wavelength probe different radial regions in the disk because of the changing stellar luminosity. Comparing the degree of settling using a fixed height or a lower scale height for bigger grains (e.g. D’Alessio et al. 2006; Sauter & Wolf 2011) is therefore less suitable for comparison along the stellar mass range. Radiative transfer models that describe dust settling using turbulent mixing exist (Dullemond & Dominik 2004b; Hasegawa & Pudritz 2010), but have not been used to compare to observations directly.

To constrain the turbulent mixing strength across the stellar mass range, we focused primarily on median spectral energy distributions (SEDs). We did this because measuring the degree of settling in individual disks is not straightforward because of a number of degeneracies in disk modeling (Sauter & Wolf 2011), and requires detailed multiwavelength modeling, including scattered-light images (Pinte et al. 2008), but conclusions can be drawn from larger samples (Furlan et al. 2005, 2006). They are available for both T Tauri stars (e.g. Furlan et al. 2005; D’Alessio et al. 2006), and brown dwarfs (Szűcs et al. 2010), and a substantial sample is available to construct them for Herbig stars as well (e.g. Acke et al. 2009; Juhász et al. 2010).

We will describe how dust settling is implemented in our radiative transfer code in Sect. 2, and compare the available approaches that use the turbulent mixing strength. In Sect. 3, we will use this model to constrain the turbulent mixing strength from the median SEDs of T Tauri stars, Herbig stars and brown dwarfs, which were constructed in Appendices A and B. We will discuss what the median disk model looks like in Sect. 4, and whether disks around lower mass stars are flatter or more settled than their high-mass counterparts. We summarize our findings in Sect. 5.

2. Model

2.1. Disk model

The disk model used in this paper is MCMax (Min et al. 2009), a 2D Monte Carlo radiative transfer code. It is based on the immediate temperature correction procedure from Bjorkman & Wood (2001), combined with the method of continuous absorption by Lucy (1999). The code has been benchmarked against other radiative transfer codes for modeling protoplanetary disks (Pinte et al. 2009). The radial grid around the inner radius and disk wall is refined to sample the optical depth logarithmically. The SED and images are calculated by integrating the formal solution to the equation of radiative transfer by ray-tracing.

In addition to radiative transfer, the code can explicitly solve for the vertical structure of the disk. With the implicit assumption that the gas temperature is set by the dust temperature, the vertical structure of the gas can be calculated by solving the equation of hydrostatic equilibrium. On top of that, the dust can be settled with respect to the gas, which we will describe in Sect. 2.2. After calculating the new structure for the dust, the radiative transfer code can be run to obtain a new dust temperature, and the whole procedure (radiative transfer, gas structure, dust settling) can be iterated until convergence is reached (Fig. 1). The temperature structure usually converges within five iterations. The exact number of iterations depends on disks parameters. In general, very settled or “flat” disks take longer to converge because they deviate more from the initial guess.

2.2. Dust settling

To settle the dust with respect to the gas, we used the formalism by Dullemond & Dominik (2004b). We will briefly describe the basic assumptions, and refer to the previously mentioned paper for details. We compare our approach to that of Hasegawa & Pudritz (2010) in Sect. 2.2.2, who used the midplane formalism by Dubrulle et al. (1995) in combination with radiative transfer models.

2.2.1. Theory

The vertical motion of dust particles within the disk is a balance of the downward motion caused by the gravity of the star versus the mainly upward motion caused by turbulent mixing. Comparing the timescales for both processes (the friction time t_fric versus the eddy turn-over time t_edd) results in an expression for the Stokes number St, which describes how well the dust is coupled to the gas at every location in the disk: $\mathrm{St}\mathrm{=}\frac{{\mathit{t}}_{\mathrm{fric}}}{{\mathit{t}}_{\mathrm{edd}}}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\frac{\mathit{m}}{\mathit{\sigma }}\frac{{\mathrm{\Omega }}_{\mathrm{k}}}{{\mathit{\rho }}_{\mathrm{gas}}{\mathrm{c}}_{\mathrm{s}}}\mathit{,}$(1)where m/σ is the mass-to-surface ratio of the dust grain, ${\mathrm{\Omega }}_{\mathrm{k}}\mathrm{=}\sqrt{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}\mathit{/}{\mathit{r}}^{\mathrm{3}}}$ the Keplerian frequency, ρ_gas is the local gas density and c_s = kT/μm_proton the local sound speed for the mean molecular weight μ. A Stokes number smaller than one means the dust is well coupled to the gas, whereas particles decouple from the gas at larger Stokes numbers. Using this expression for the dust-gas coupling and the standard α prescription from Shakura & Sunyaev (1973) in the form of Pringle (1981) (ν = αc_sH_p), the diffusion coefficient for a dust particle at every location in the disk can be written as²$\mathit{D}\mathrm{=}{\mathit{\alpha }}_{\mathrm{turb}}\frac{{\mathit{c}}_{\mathrm{s}}{\mathit{H}}_{\mathrm{p}}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}}\mathit{,}$(2)where H_p = c_s/Ω_k is the local pressure scale height and α_turb is the turbulent mixing strength.

Taking the gas density and temperature structure from the radiative transfer and hydrostatic structure calculation, we can now solve for the vertical structure of a dust grain of certain size ρ_a at every location in the disk by solving a vertical diffusion equation: $\frac{\mathit{\partial }{\mathit{\rho }}_{\mathrm{a}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{\mathit{\partial t}}\mathrm{=}\frac{\mathit{\partial }}{\mathit{\partial z}}\left({\mathit{D}}_{\mathrm{a}}\mathrm{\left(}\mathit{z}\mathrm{\right)}{\mathit{\rho }}_{\mathrm{gas}}\frac{\mathit{\partial }}{\mathit{\partial z}}\left(\frac{{\mathit{\rho }}_{\mathrm{a}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}{{\mathit{\rho }}_{\mathrm{gas}}}\right)\right)\mathrm{-}\frac{\mathit{\partial }}{\mathit{\partial z}}\left({\mathit{\rho }}_{\mathrm{a}}{\mathit{v}}_{\mathrm{sett}\mathit{,}\mathrm{a}}\mathrm{\left(}\mathit{z}\mathrm{\right)}\right)\mathit{,}$(3)where v_sett,a(z) is the local settling speed for that dust grain in the absence of turbulent mixing. This equation is then solved in a time-dependent way using implicit integration. Because the timescales on which the disk settles toward equilibrium are shorter than one million years (Dullemond & Dominik 2004b), we jump toward equilibrium in ten relatively large time steps.

We assumed that the turbulent mixing strength is constant throughout the disk. There is no a priori reason to do so, nor a computational one, and there are indications for variations in both the vertical (Fromang & Nelson 2009; Simon et al. 2011) and radial (Isella et al. 2009) direction. However, because we cannot constrain these variations from the data used in this paper, we kept α_turb constant.

2.2.2. Impact on disk structure

The main difference between this approach and that by Dubrulle et al. (1995) is that we used the local sound speed and gas density to calculate the dust-gas coupling. Their approach is an analytical solution for modeling the dust distribution in the isothermal midplane, and uses only the midplane density and temperature. Hence, we will refer to the latter as midplane settling, though it is sometimes applied outside this regime (Hasegawa & Pudritz 2010). To compare the different approaches, we have displayed the vertical structure at the first iteration in a vertical slab centered at 10 AU in Fig. 2. Because this model has not been iterated, the vertical structure of the gas is calculated from the optically thin dust temperature, and corresponds to a Gaussian.

The well-mixed model has a uniform dust-to-gas ratio of 0.01, whereas the settled models show a clear trend from an increased dust-to-gas ratio near the midplane to a decreased ratio near and above the surface, which is located around z/r ~ 0.18 (radial τ = 1 surface in the optical). This is caused by stratification, the settling of different dust species to different heights in the disk. In the midplane, the midplane settling approach agrees very well with the self-consistent settling, which is also the region that the former has been constructed for (Dubrulle et al. 1995).

Closer to the disk surface, the two settled models start to deviate. With increasing height above the midplane, the gas density decreases. For the self-consistent settling, this means that smaller particles will decouple from the gas. For each particle size, there is a characteristic height above the midplane where the coupling is so weak that they are almost fully depleted (Dullemond & Dominik 2004b). For particles down to a size of roughly 1 micron, this decoupling takes place below the disk surface, but smaller grains can still be found above the disk surface.

The midplane settling approach shows a different behavior, because the dust-gas coupling is calculated in the midplane only. Particles do not decouple from the gas at a specific height, but keep following a Gaussian distribution all the way through the disk surface. This gives an almost constant dust-to-gas ratio for the smallest particles (<1 μm), and a much weaker depletion for the intermediate sizes.

Settling the dust also impacts the temperature structure of the disk (Fig. 3), and therefore its density structure. Non-settled disk models have a more or less constant temperature distribution below the disk surface (Fig. 3, dashed line). Big grains, however, are efficient coolers, and settling them leads to a colder midplane and warmer surface. The stratification in dust sizes also leads to a gradual increase in temperature from the midplane to the disk surface (see also Hasegawa & Pudritz 2010). Because the temperature in and above the midplane is no longer isothermal, it becomes important to solve for the vertical structure of the disk using the local scale height, which increases with height. Using only the midplane temperature to calculate the vertical structure of the disk underestimates the scale height within the midplane by up to a factor of two (Fig. 3). This leads to flatter disks with – in turn – even lower midplane temperatures (Fig. 3, dotted line).

2.3. Dust grain properties

Dust grains in our model are characterized by three different parameters: composition, size and shape. Different grain sizes settle to different heights in the disk (see Fig. 2), and their opacities are displayed in Fig. 4.

2.3.1. Composition and shape

Our dust grains consist mainly of amorphous silicates, measured for the composition of interstellar dust toward the Galactic center by Min et al. (2007), co-added with with a continuum opacity source, for which we used amorphous carbon. The ratio of amorphous carbon/amorphous silicates influences the peak strength of the 10-micron silicate feature. We used a mass fraction of 15% carbon, a value found from modeling this ratio in Herbig Ae/Be stars (Meijer et al. 2005, 2009). The final composition, with reference to the optical constants, is 11.73% MgFeSiO₄ (Dorschner et al. 1995), 32.6% Mg₂SiO₄ (Henning & Stognienko 1996), 36.5% MgSiO₃ (Dorschner et al. 1995), 1.5% NaAlSi₂O₆ (Mutschke et al. 1998), 15% C (Preibisch et al. 1993). The shape of our particles is irregular, and approximated using a distribution of hollow spheres (DHS, Min et al. 2005). We used a vacuum fraction of 0.7, though this choice does not influence our results because we do not focus on detailed feature shapes.

2.3.2. Size

Observations of protoplanetary disks show evidence for (sub)micron and millimeter-sized or larger grains. Intermediate sized grains are not observed, but are predicted by dust coagulation models (e.g. Weidenschilling 1984; Birnstiel et al. 2010a). We therefore used a grain size distribution ranging from sub-micron to millimeter sizes that follows a power law, f(a) = a ^− q, and used a default index of q = 3.5 similar to the Mathis, Rumpl and Nordsieck (MRN) distribution (Mathis et al. 1977). To allow settling of different grain size to a different height in the disk, we divided our size distribution into two bins per order of magnitude. We define the grain radius for the settling calculation with the logarithmic mean of the minimum and maximum grain sizes to (i.e. a = 1.78 micron for a grain size bin between 1.00 and 3.16 micron). The opacities of these grains were calculated using 10 sizes per bin to avoid strong resonant features that could result from using a single grain size.

2.4. Disk parameters

Because our disks are in vertical hydrostatic equilibrium, the only free fitting parameter to control the vertical structure for a given mass is the turbulent mixing strength α. The radial distribution of dust and gas is parametrized by a power-law as Σ(r) = r ^− p with the commonly used p = 1 for a steady-state accretion disk. The inner radius is placed at a dust evaporation temperature of 1500 K, whereas the outer radius is fixed to 300 AU in accordance with D’Alessio et al. (2006). The total disk mass is a fixed fraction of 1 percent of the stellar mass (e.g. Scholz et al. 2006; Williams & Cieza 2011), assuming a dust-to-gas ratio of 1:100.

3. Fits to median SEDs

To compare our disk model to the three different data sets, we used median SEDs. Instead of fitting all individual sources or comparing them to a model grid, we only fitted our model to one median SED to extract “typical” fit parameters. Furlan et al. (2005, 2006) have shown that this approach works well, taking into account the observed spread in mid-infrared colors.

When fitting the data, we tried to keep as many parameters fixed as possible, see Table 1. Parameters that could not be fixed because of the different stellar properties are given in Table 2. These are the inner radius, which we kept at 1500 Kelvin and the dust mass, which is a fixed fraction of the stellar mass of 10^-4 M_∗.

3.1. T Tauri stars

The Taurus median SED has been fitted using a hydrostatic disk model by D’Alessio et al. (2006), which demonstrated the need for grain growth and dust settling. Our goal is to reproduce these results, but now in the context of a more self-consistent settling approach to constrain the turbulent mixing strength, and not an arbitrary height for the midplane layer of big grains. We have extended the existing data set from Furlan et al. (2006) with sub-millimeter fluxes at 850 μm using the catalog from Andrews & Williams (2005). The data from D’Alessio et al. (1999) include some non-detections and upper limits (see Appendix A). The updated millimeter median is well fitted by a dust mass of 5 × 10^-5 M_⊙, consistent with the result from Andrews & Williams (2005).

To fit the median SED in the near-infrared, we needed to include viscous heating (Fig. 5). We did not self-consistently solve the radial structure of the disk, but calculated the viscous heating with a mass accretion rate of 3 × 10^-8 M_⊙/yr, the same as Furlan et al. (2006). To fit the mid-infrared flux, we flattened the disk by settling the dust. A turbulent mixing strength of α = 10^-2 – a value commonly used for viscous accretion disks – overpredicts the mid-infrared fluxes longward of 20 micron, which arises from the region around  ~10 AU. We need to lower the turbulent mixing strength to achieve a good fit, which results in α = 10^-4.

Because the 20 micron flux mainly probes the height of the disk surface, there are also other ways to lower it, see also Sect. 4.5. We can reduce the number of small grains, which will lower the optical depth in the upper layers and decrease the visible disk surface. Decreasing the total disk mass would be inconsistent with the millimeter flux, but we can modify the grain size distribution because it is dominated by big grains. Decreasing the power-law index of the grain size distribution removes small grains from the disk surface, but does not significantly affect the mass in big grains near the disk midplane. Changing the power-law index to 3.25 is sufficient to fit the median SED with a turbulent mixing strength of α = 10^-2.

3.2. Herbig stars

Spectral energy distributions for a large number of Herbig stars are available, but a median SED has not been constructed before. We constructed one in Appendix B, which we present in Fig. 6. To test the hypothesis that more massive stars have less settled disks, we retained as many disks parameters from the T Tauri stars as possible. We changed (see also Table 2) the disk mass – which scales with the stellar mass and is increased with a factor 5 with respect to the T Tauri star, consistent with the submillimeter data – and the inner radius – which also increases to be consistent with dust evaporation.

Although accretion does not play a role in the SEDs of Herbig stars because of the increased stellar luminosity and larger inner radius, the near-infrared fluxes (2–8 micron) of hydrostatic disk models underpredict the observations (Vinković et al. 2006, see also Fig. 6). Different mechanisms have been proposed to explain the difference, from an increased rim scale height (Verhoeff et al. 2010; Acke et al. 2009) to halos (Vinković et al. 2006; Verhoeff et al. 2011) to material within the dust evaporation radius (e.g. Kraus et al. 2008; Tannirkulam et al. 2008a). Although an artificially increased scale height works well for the mid-infrared SED (Acke et al. 2009), it does not work well when fitting simultaneously the near-infrared region of the Herbig median SED (see Appendix C). Because an inner gaseous disk is beyond the current capabilities of our code, we modeled the near-infrared flux deficit of the hydrostatic model with an optically thin spherical halo. We emphasize that this halo is a parametrization of optically thin emission from the inner regions, including dust or gas within the inner rim, and refer the interested reader to the appendix.

The optically thin halo plays the role of accretion in T Tauri stars: it provides near-infrared flux without affecting the flux longward of 20 micron, in contrast to the puffed-up inner rim. We used a small halo of 0.3...1 micron grains between 0.14 and 0.3 AU³ that has an optical depth of 0.14 in the visual and a dust mass of 4 × 10^-12 M_⊙ (Fig. 6). As a result of the fitting procedure, we found the same turbulent mixing strength for the disk of α = 10^-4 as for the T Tauri stars. A higher mixing strength of α = 10^-2 overpredicts the Spitzer IRS fluxes from 20 to 40 micron and the IRAS flux at 60 micron.

3.3. Brown dwarfs

Disks around brown dwarfs are more difficult to detect owing to their low luminosity. Complete samples down to photospheric values only exist up to 8 micron, and are incomplete at 24 micron (e.g. Lada et al. 2006; Scholz et al. 2007; Luhman et al. 2008). Furthermore, very few millimeter detections exist (Klein et al. 2003; Scholz et al. 2006). We used the median SED constructed by Szűcs et al. (2010) for Chameleon I, where the median 24 micron flux should be treated as an upper limit. For convenience, we also plot the millimeter upper limit from Klein et al. (2003), scaled to the distance of the Chameleon I star-forming region.

Unlike the Herbig and T Tauri stars, hydrostatic disk models of brown dwarfs do not show a flux deficit in the near-infrared compared to the observed median. We scaled the disk mass and inner radius to the brown dwarf regime (see Table 2). The stellar photosphere was fitted by a NextGen stellar atmosphere (Hauschildt et al. 1999), and we adjusted the luminosity down to 0.025 L_⊙. We fitted the median shortward of 8 micron without any additional adjustments (see Fig. 7). The assumed dust mass of 8 × 10^-6 M_⊙ is consistent with the upper limit from Klein et al. (2003), and amounts to only 2–3 earth masses of solid material.

We found that a turbulent mixing strength of α = 10^-4 is consistent with the 24 micron upper limit, whereas a mixing strength of α = 10^-2 clearly overpredicts this upper limit. We found the same turbulent mixing strength as for the earlier spectral types.

4. Discussion

Although we found the same turbulent mixing strength for the three different samples, this does not necessarily mean that these disks are equally flat. The vertical structure is a product of more than just the turbulent mixing strength, most notably the local density, temperature and Keplerian frequency. In the following section we will describe the disk properties, and show that disks around lower mass stars appear more settled at a specific radius, consistent with Szűcs et al. (2010), but that the structure in the mid-infrared emitting region is remarkably self-similar.

We will also discuss why our choice of grain size distribution does not affect our conclusion that the turbulent mixing strength does not vary along the stellar mass range in Sect. 4.3 and its implication for the planet-forming potential of the disk in Sect. 4.4. The dependence of the fitted turbulent mixing strength on other disk parameters is discussed in Sect. 4.5.

4.1. Disk properties at fixed radius

To compare the dust and gas distribution in our best-fit models, we display their pressure scale height and the radial τ = 1 surfaces in Fig. 8. The pressure scale height is highest in the brown dwarf disk and lowest in the Herbig star, reflecting the differences in stellar mass. Although temperatures are also lower for less massive stars, this is not enough to compensate the lower gravitational potential in which the disk resides, and the gas remains more flaring.

The dust, on the other hand, follows a different path: the τ = 1 surface in the outer disk reaches a relative height (z/r) of 0.17, 0.15 and 0.14 for the Herbig star⁴, T Tauri star and brown dwarf, respectively. This indicates that disks around less massive stars are a little flatter, though the difference remains within 20%. Compared to the much higher pressure scale height at the same radius, the disk of lower mass stars therefore appears more settled, because the disk surface reaches down to one pressure scale height for the Herbig star, but to less than half a pressure scale height in the brown dwarf. This can be explained by lower (surface) densities and temperatures (see Figs. 10a, 11a and Table 3) that make the dust decouple closer to the midplane, as well as a lower optical depth.

The shape of the disk surface can also be viewed in scattered-light images. We display the radial profiles of synthetic scattered-light images in the visual in Fig. 12a. The radial profiles are scaled to the photospheric flux at the same wavelength, such that it traces only the shape of the disk surface, and not the properties of the central star. A clear trend is visible at all radii, with higher fluxes and shallower slopes for Herbig stars and lower fluxes and steeper slopes for brown dwarfs (see also Table 3). This indicates that the disk surface at a specific location is less flaring.

4.2. Self-similar solutions in the MIR emitting region

Because the mid-infrared emission probes a different region in all three samples, it is also interesting to compare regions with similar temperatures, which are located at different radii for different stars. These regions can be identified by using that the distance r, at which a specific flux from the star is received, scales with the square root of the stellar luminosity L _∗ . By scaling the radius to the value of the stellar flux, ${\mathit{r}}_{\mathrm{T}}\mathrm{=}\mathit{r}\mathit{/}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}$, we can construct dimensionless equivalents of our disk properties (denoted by a ^), and check if they are self-similar across the stellar mass range.

For example, the optically thin dust temperature for the T Tauri star can be described as (see Fig. 9a) ${\mathit{T}}_{\mathrm{dust}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{307} \mathit{K}\mathrm{\ast }\mathrm{(}\mathit{r}\mathit{/}\mathrm{AU}{\mathrm{)}}^{-0.5}\mathit{.}$(4)By scaling the radius to the intensity of the radiation field, we can retrieve a self-similar solution for the optically thin dust temperature profile $\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}{\mathit{T}}_{\mathrm{dust}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathrm{307} \mathrm{K}\mathrm{\ast }{\left({\mathit{r}}_{\mathrm{T}}\mathit{/}\mathrm{AU}\mathrm{\ast }\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}\right)}^{-0.5}\begin{array}{}\end{array}\mathrm{=}\mathrm{315} \mathrm{K}\mathrm{\ast }\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathit{/}\mathrm{AU}{\mathrm{)}}^{-0.5}\mathrm{=}\mathit{T̂}\mathrm{dust}\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathrm{)}\mathit{,}\end{array}\begin{array}{}\end{array}\end{array}$(5)where the temperature at 1 AU for the T Tauri star is close to the selfconsistent solution because its luminosity is close to 1 L_⊙. As can be seen in Fig. 9b, these scaled solutions agree to within 20% for the T Tauri star, Herbig star and brown dwarf.

4.2.1. Self-similar temperature and surface density

With the same procedure, we can also compare different quantities to see if they are also self-similar in all three samples. We find that also the midplane temperatures for all three samples (see Fig. 10b) can be described with a self-similar function: $\mathit{T̂}\mathrm{mid}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}\mathrm{=}\mathrm{90} \mathrm{K}\mathrm{\ast }\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathit{/}\mathrm{AU}{\mathrm{)}}^{-0.44}\mathit{,}$(6)which is a strong indication that this must be a consequence of a self-similar vertical structure across the stellar mass range, which will we explore below in more detail. Interestingly, this scaling law also works for the surface density (Fig. 11), such that regions with similar temperatures also have similar surface densities: $\mathrm{\Sigma ̂}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}\mathrm{=}\mathrm{0.25} {\mathrm{g}\mathit{/}\mathrm{cm}}^{\mathrm{2}}\mathrm{\ast }\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathit{/}\mathrm{AU}{\mathrm{)}}^{-1}\mathit{.}$(7)This is because across the three samples, the luminosity scales roughly with stellar mass squared⁵. Since the disk mass is a fixed fraction of the stellar mass, the surface density scales as $\mathrm{\Sigma ̂}\mathrm{1}\mathrm{AU}\mathrm{=}\left(\frac{{\mathit{M}}_{\mathrm{\odot }}}{{\mathit{M}}_{\mathrm{\ast }}}\right) {\mathrm{\Sigma }}_{\mathrm{1} \mathrm{AU}}\mathrm{=}\sqrt{\frac{{\mathit{L}}_{\mathrm{\odot }}}{{\mathit{L}}_{\mathrm{\ast }}}} {\mathrm{\Sigma }}_{\mathrm{1} \mathrm{AU}}\mathit{.}$(8)However, the region corresponding to the same temperature lies farther away from the star, compensating for the change in disk mass: $\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}\mathrm{\Sigma ̂}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}\mathrm{=}\mathrm{\Sigma ̂}\mathrm{1} \mathrm{AU} {{\mathit{r}}_{\mathrm{T}}}^{-1}\mathrm{=}\left(\frac{{\mathrm{\Sigma }}_{\mathrm{1} \mathrm{AU}}}{\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}}\right){\left(\frac{\mathit{r}}{\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}}\right)}^{-1}\begin{array}{}\end{array}\mathrm{=}{\mathrm{\Sigma }}_{\mathrm{1} \mathrm{AU}} {\mathit{r}}^{-1}\mathrm{=}\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}\end{array}\begin{array}{}\end{array}\end{array}$(9)thus the self-similar solutions for the surface density agree to within 20%.

4.2.2. Self-similar scale height

If we now look at a vertical slab in the disk located at ${\mathit{r}}_{\mathrm{T}}\mathrm{=}\mathit{r}\mathit{/}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}$, we find that two of the main parameters for the settling calculation in Eqs. (1)and (2)– the gas surface density (8)and temperature (6)– are self-similar. The third important parameter, the Keplerian frequency, is not the same, but scales with $\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}$$\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}\mathrm{\Omega ̂}\mathrm{k}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}\mathrm{=}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{\ast }}}{{\mathit{r}}_{\mathrm{T}}^{\mathrm{3}}}}\mathrm{=}\sqrt{\mathit{G}\frac{{\mathit{M}}_{\mathrm{\odot }}}{\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}}{\left(\frac{\mathit{r}}{\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}}\right)}^{-3}}\begin{array}{}\end{array}\mathrm{=}\sqrt{\frac{\mathit{G}{\mathit{M}}_{\mathrm{\odot }}}{{\mathit{r}}^{\mathrm{3}}}}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}\mathrm{=}{\mathrm{\Omega }}_{\mathrm{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)} \sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}\mathit{,}\end{array}\begin{array}{}\end{array}\end{array}$(10)which is equivalent to r/r_T. The gas pressure scale height at this location is therefore not self-similar, but the relative scale height H_p/r is: $\frac{\mathit{Ĥ}\mathrm{p}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}}{{\mathit{r}}_{\mathrm{T}}}\mathrm{=}\frac{\mathit{ĉ}\mathrm{s}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}}{{\mathit{r}}_{\mathrm{T}}\mathrm{\Omega ̂}\mathrm{k}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{{\mathit{r}}_{\mathrm{T}} \mathrm{\left(}{\mathrm{\Omega }}_{\mathrm{k}}\mathrm{\right(}\mathit{r}\mathrm{)} \mathit{r}\mathit{/}{\mathit{r}}_{\mathrm{T}}\mathrm{)}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{s}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{r}{\mathrm{\Omega }}_{\mathrm{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{=}\frac{{\mathit{H}}_{\mathrm{p}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathit{r}}\mathrm{·}$(11)This means that if we express the height above the midplane (z) in terms of the relative height (z/r), the local scale height, and therefore the vertical structure of the gas, is self-similar: $\mathit{Ĥ}\mathrm{p}\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathit{,}{\mathit{z}}_{\mathrm{T}}\mathit{/}{\mathit{r}}_{\mathrm{T}}\mathrm{)}\mathrm{=}{\mathit{H}}_{\mathrm{p}}\mathrm{(}\mathit{r,z}\mathit{/}\mathit{r}\mathrm{)}\mathit{.}$(12)

4.2.3. Self-similar dust settling

The vertical structure of the dust in the mid-infrared emitting region at r_T = 10   AU is shown in Fig. 13. Like that of the gas, the vertical structure of the dust also appears self-similar when expressed in terms of its relative height z/r. This is because also the Stokes number (Eq. (1)) $\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}\mathrm{St̂}\mathrm{(}{\mathrm{r}}_{\mathrm{T}}\mathit{,}{\mathrm{z}}_{\mathrm{T}}\mathit{/}{\mathrm{r}}_{\mathrm{T}}\mathrm{)}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\frac{\mathit{m}}{\mathit{\sigma }}\frac{\mathrm{\Omega ̂}\mathrm{k}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}}{\mathit{\rho ̂}\mathrm{gas}\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathit{,}{\mathit{z}}_{\mathrm{T}}\mathit{/}{\mathit{r}}_{\mathrm{T}}\mathrm{)} {\mathit{c}}_{\mathrm{s}}\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathit{,}{\mathit{z}}_{\mathit{T}}\mathit{/}{\mathit{r}}_{\mathrm{T}}\mathrm{)}}\begin{array}{}\end{array}\mathrm{=}\frac{\mathrm{3}}{\mathrm{4}}\frac{\mathit{m}}{\mathit{\sigma }}\frac{{\mathrm{\Omega }}_{\mathrm{k}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{/}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}}{{\mathit{\rho }}_{\mathrm{gas}}\mathrm{(}\mathit{r,z}\mathit{/}\mathit{r}\mathrm{)}\mathit{/}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}} {\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathit{r,z}\mathit{/}\mathit{r}\mathrm{)}}\begin{array}{}\end{array}\mathrm{=}\mathrm{St}\mathrm{(}\mathit{r,z}\mathit{/}\mathit{r}\mathrm{)}\end{array}\begin{array}{}\end{array}\end{array}$(13)is self-similar because the gas density scales with 1/H_p. The same is true for the diffusion coefficient (2)$\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}\mathit{D̂}\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathit{,}{\mathit{z}}_{\mathrm{T}}\mathit{/}{\mathit{r}}_{\mathrm{T}}\mathrm{)}\mathrm{=}{\mathit{\alpha }}_{\mathrm{turb}}\frac{\mathit{ĉ}\mathrm{s}\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathit{,}{\mathit{z}}_{\mathrm{T}}\mathit{/}{\mathit{r}}_{\mathrm{T}}\mathrm{)}\mathit{Ĥ}\mathrm{p}\mathrm{(}{\mathit{r}}_{\mathrm{T}}\mathit{,}{\mathit{z}}_{\mathrm{T}}\mathit{/}{\mathit{r}}_{\mathrm{T}}\mathrm{)}}{\mathrm{1}\mathrm{+}\mathrm{St̂}\mathrm{2}\mathrm{(}{\mathrm{r}}_{\mathrm{T}}\mathit{,}{\mathrm{z}}_{\mathrm{T}}\mathit{/}{\mathrm{r}}_{\mathrm{T}}\mathrm{)}}\begin{array}{}\end{array}\mathrm{=}{\mathit{\alpha }}_{\mathrm{turb}}\frac{{\mathit{c}}_{\mathrm{s}}\mathrm{(}\mathit{r,z}\mathit{/}\mathit{r}\mathrm{)}{\mathit{H}}_{\mathrm{p}}\mathrm{(}\mathit{r,z}\mathit{/}\mathit{r}\mathrm{)}}{\mathrm{1}\mathrm{+}{\mathrm{St}}^{\mathrm{2}}\mathrm{(}\mathit{r,z}\mathit{/}\mathit{r}\mathrm{)}}\begin{array}{}\end{array}\mathrm{=}\mathit{D}\mathrm{(}\mathit{r,z}\mathit{/}\mathit{r}\mathrm{)}\end{array}\begin{array}{}\end{array}\end{array}$(14)because it is only a product of self-similar quantities when expressed in units of r_T and z_T/r_T. Consequently, the vertical structure of both gas and dust in the mid-infrared emitting regions of protoplanetary disks are self-similar.

This result is consistent with that of Furlan et al. (2011), who also founnd no variation in the degree of settling between brown dwarfs and T Tauri stars, but partly contradicts Szűcs et al. (2010). These authors also found the same relative scale height (H/r) for the dust in the mid-infrared emitting region. However, the required reduction in dust scale height compared to their fully flared models implies a difference in gas scale height between brown dwarfs and T Tauri stars, which we do not find.

The self-similarity of vertical structure is further illustrated by synthetic scattered-light images that trace the shape of the dust disk surface (Fig. 12). If we plot the radial intensity profile as a function of the scaled radius r_T, corrected for the relative height⁶ (F_scat/F_∗ ∗ L_∗/L_⊙), we see that the radial profiles for the three different stars line up nicely (Fig. 12b), and that the shape of the disk surface is indeed self-similar.

Since the amount of radiation reprocessed at every radius is set by the covering fraction (z/r) – which is the same for all stars – even the SED in the mid-infrared becomes self-similar (Fig. 14). The scaled model SEDs line up to within 20% in the mid-infrared at 20–30 micron where we fitted the SEDs. The similarity extends into the far-infrared (within 40%) and even the near-infrared up to 3 micron (within a factor of 2). The latter is surprising and probably coincidental, given the different nature of the emission mechanism at that wavelength: viscous heating for the T Tauri star, optically thin emission for the Herbig star and mostly photospheric flux for the brown dwarf.

4.3. Grain size distributions

Throughout this article, we have assumed that there are no radial gradients in the grain size distribution. Models of grain growth predict that grain sizes decrease with radius (e.g. Birnstiel et al. 2010b), which has also been observed by Guilloteau et al. (2011) and Banzatti et al. (2011). Additionally, disks around later type stars are colder and less massive, leading to smaller grain sizes. Modeling the grain size distribution is beyond the scope of this paper, though a recipe for doing so is available (Birnstiel et al. 2011). For now, we will only discuss if we expect to see differences in the grain size distribution in the mid-infrared emitting region between the three different samples.

Because the vertical structure in the mid-infrared region is self-similar, we investigated if this is also true for the grain-size distribution and its planet-forming potential. Grains in protoplanetary disks are in a coagulation-fragmentation equilibrium (Weidenschilling 1984), driven by turbulence. According to Birnstiel et al. (2011), the most important parameters that define the size distribution are the turbulent mixing strength, the gas surface density, the temperature, and three parameters that refer to the microphysics of dust: the fragmentation velocity u_f, the solid density of dust ρ_s and the size distribution of fragments ξ. We assumed the last three to be the same across the three samples.

We find no variations in the turbulent mixing strength (Sect. 3). Regions with similar temperature also have a similar surface density (Sect. 4.2). The grain size distribution should then also be self-similar. We can see why this is the case by looking at the expression for the maximum grain size for turbulence-driven growth⁷ – as given by Eq. (12) of Birnstiel et al. (2009): $\mathit{â}\max \mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}\mathrm{=}\frac{\mathrm{\Sigma ̂}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}{\mathit{u}}_{\mathrm{f}}^{\mathrm{2}}}{\mathit{\pi }{\mathit{\alpha }}_{\mathrm{turb}}{\mathit{\rho }}_{\mathrm{s}}\mathit{ĉ}\begin{array}{}\mathrm{2}\\ \mathrm{s}\end{array}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}}\mathrm{=}\frac{\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}{\mathit{u}}_{\mathrm{f}}^{\mathrm{2}}}{\mathit{\pi }{\mathit{\alpha }}_{\mathrm{turb}}{\mathit{\rho }}_{\mathrm{s}}{\mathit{c}}_{\mathrm{s}}^{\mathrm{2}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}\mathrm{=}{\mathit{a}}_{\max }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{,}$(15)which is also self-similar when looking at regions with similar temperature (${\mathit{r}}_{\mathrm{T}}\mathrm{=}\mathit{r}\mathit{/}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}$). This means that also the grain size distribution in the mid-infrared emitting region is not expected to vary across the stellar mass range.

4.4. Planet-forming potential

If the same disk properties are found at a different location in the disk around heavier stars, does this mean that the resulting planetary systems are simply scaled versions of each other?

To compare the later stages of planet formation, we looked at the Hill sphere and feeding zone of a forming planet, to see how much mass is available for its formation, in a way similar to Raymond et al. (2007). The Hill sphere is given by ${\mathit{R}}_{\mathrm{H}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathrm{=}\mathit{r}\mathrm{3}\sqrt{\frac{{\mathit{m}}_{\mathrm{P}}}{\mathrm{3}{\mathit{M}}_{\mathrm{\ast }}}}\mathit{,}$(16)where m_P is the mass of the planet. Expressed in the same units that result in self-similar solutions for the vertical structure (${\mathit{r}}_{\mathrm{T}}\mathrm{=}\mathit{r}\mathit{/}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}$) it becomes $\mathit{R̂}\mathrm{H}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}\mathrm{=}\mathit{r}\mathit{/}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}\mathrm{3}\sqrt{\frac{{\mathit{m}}_{\mathrm{P}}}{\mathrm{3}{\mathit{M}}_{\mathrm{\ast }}\mathit{/}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}}}\mathrm{=}\frac{{\mathit{R}}_{\mathrm{H}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{3}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}}\mathit{,}$(17)so the size of the Hill sphere is not self-similar, though its dependence on luminosity is weak. It decreases in size for less luminous stars, despite the lower gravity of the central star. The total mass available in the feeding zone deviates more because it scales with an extra factor 2πr_T: $\begin{array}{}\begin{array}{}\end{array}\begin{array}{}\begin{array}{}\end{array}\mathit{M̂}\mathrm{feed}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}\mathrm{=}\mathrm{2}\mathit{\pi }{\mathit{r}}_{\mathrm{T}}\mathrm{\Sigma ̂}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}\mathrm{\ast }\mathrm{2}\mathit{R̂}\mathrm{H}\mathrm{\left(}{\mathit{r}}_{\mathrm{T}}\mathrm{\right)}\begin{array}{}\end{array}\mathrm{=}\mathrm{2}\mathit{\pi }\mathrm{\Sigma }\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{r}\mathit{/}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}\mathrm{\ast }\mathrm{2}\frac{{\mathit{R}}_{\mathrm{H}}\mathrm{\left(}\mathit{r}\mathrm{\right)}}{\mathrm{3}\sqrt{{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}}}\begin{array}{}\end{array}\mathrm{=}{\mathit{L}}_{\mathrm{\ast }}\mathit{/}{\mathit{L}}_{\mathrm{\odot }}^{\mathrm{-}\mathrm{5}\mathit{/}\mathrm{6}}{\mathit{M}}_{\mathrm{feed}}\mathrm{\left(}\mathit{r}\mathrm{\right)}\mathit{.}\end{array}\begin{array}{}\end{array}\end{array}$(18)The size of the feeding zone therefore scales almost linearly with luminosity. This means that although the initial conditions for planet formation are the same across the steller mass range, this is not true for the later stages of planet formation, in agreement with Raymond et al. (2007).

4.5. Dependence on disk parameters

The method we present here is useful for constraining variations in the turbulent mixing strength, but it is more difficult to constrain its absolute value. The mid-infrared flux with which we fitted our models traces the height of the disk surface, which can also be influenced by other factors than the turbulent mixing strength. However, because our solutions for the mid-infrared emitting region are self-similar, these uncertainties will affect our estimate of the turbulent mixing strength in the same way for all stars, and do not influence our conclusion that it does not vary across the stellar mass range.

In this section we will explore how other parameters can affect our estimate of the turbulent mixing strenght. In particular, we will explore how to keep the turbulent mixing strength fixed at α_turb = 0.01 and fit the three median SEDs by varying these other parameters. As already mentioned in Sect. 3.1, small changes in the grain size distribution also affect the height of the disk surface. Using a flatter slope for this distribution removes small grains from the disk surface, and we can fit the median SEDs with a power-law index of 3.25. Another parameter that directly influences the SED in the mid infrared is the gas-to-dust ratio. A lower ratio leads to a weaker coupling of the dust and hence to flatter disks. We fitted a gas-to-dust ratio of 1 for α_turb = 0.01.

One parameter that we cannot change is the dust mass, because this would be inconsistent with the millimeter photometry. As a consequence, all parameters that influence the millimeter opacity require a refitting of the dust mass. For example, we can reduce the opacity in the upper layers by reducing the carbon content of our dust grains, but this also reduces the opacity at millimeter wavelengths. If we then increase the dust mass to fit the flux at 850 μm, we again increase the opacity in the upper layers. The two effects cancel out, and there is no net decrease in the mid-infrared flux. In a similar fashion, the maximum grain size a_max does not influence the estimate of α_turb because it decreases the opacity in the upper layers and at millimeter wavelengths simultaneously.

The same is true for parameters that influence the distribution of mass in the outer disk: the outer radius and surface density powerlaw. Although a smaller disk has a smaller covering fraction for the same surface density, a fit to the millimeter data requires a higher surface density, compensating the decrease in covering fraction. The same is true for steeper surface density powerlaws.

4.6. Spread in disk parameters

Furlan et al. (2005) have shown that the spread in observed disk colors in the mid infrared is well reproduced by a spread in inclination, accretion and surface layer depletion. However, millimeter observations show that there is also a considerable spread in disk mass (Andrews & Williams 2005), which will affect the amount of mid-infrared emission through changing dust surface densities and stronger or weaker dust-gas coupling. We investigated if this is sufficient to explain the observed spread in disk colors, without the need for a spread in surface layer depletion.

We ran a series of disk models based on the best fit to the Taurus median SED with α_turb = 10^-4, varying only the disk mass over two orders of magnitude. The result can be found in Fig. 15. Although the spread of this grid reproduces the region between the lower and upper quartile at 20–40 micron well, it is much too large for the quartiles in the millimeter. This means that the real spread in disk mass must be smaller, and therefore cannot be solely responsible for the spread in observed colors. Albeit slightly smaller than previously estimated, a spread in surface layer depletion must still be present⁸.

This result also serves to validate one of the key assumptions of this paper: that a model with median fit parameters is a good representation of the median SED. Because the model SEDs do not overlap (which is also the case when varying most other disk parameters) the median SED is exactly equal to the SED with the median fit parameters.

5. Conclusions

We have investigated the turbulent mixing strength in protoplanetary disks across the stellar mass range by comparing radiative transfer models with dust settling to median SEDs of three samples of T Tauri stars, Herbig stars and brown dwarfs. Our conclusions are:

• We find no significant variations of the turbulent mixing strengthacross the stellar mass range from Herbig stars down to browndwarfs.

• We find a relatively low turbulent mixing strength of α_turb = 10^-4 for an MRN-like grain size distribution up to a millimeter and a gas-to-dust ratio of 100, whereas a size distribution with fewer small grains or a gas-to-dust ratio of 1 is required for a higher mixing strength of α_turb = 10^-2.

• The mid-infrared emitting regions of protoplanetary disks are remarkably self-similar, with a similar surface density, vertical structure of gas and dust, and even grain size distribution. This provides comparable environments for the first stages of planet formation, though the later phases dominated by the planet’s Hill sphere are likely different.

• When viewed at the same distance from the central star, disks around later type stars appear more settled. They have a flatter disk surface as traced in scattered light. However, their gas scale height is much higher.

• We have constructed a median SED for Herbig stars, analogous to its Taurus and Chameleon I counterparts. This median SED points to high inner-rim temperatures of  ~1700 K, and is only poorly fitted with a puffed-up inner rim, but is more consistent with emission from optically thin material from above the rim or from within it.

---

Acknowledgments

We thank Kees Dullemond for helping with the implementation of self-consistent settling and for discussions during the development of the paper, and Til Birnstiel for useful comments on the manuscript. The authors would also like to thank Bram Acke, Koen Maaskant and Michiel Min for help and discussion in constructing the Herbig median SED. We would like to thank the not-so-anonymous referee, Phil Armitage, for his concise and constructive referee report. This research project is financially supported by a joint grant from the Netherlands Research School for Astronomy (NOVA) and The Netherlands Institute for Space Research (SRON).

References

Appendix A: Taurus median

As mentioned in Sect. 3.1, we have added the 850 micron flux point to the existing Taurus median from Furlan et al. (2006). Although the original Taurus median from D’Alessio et al. (1999) is based on the same sources and includes millimeter photometry, its coverage is incomplete. It is strongly biased toward brighter – and thus heavier – disks and was “meant to be indicative rather than definitive”. Andrews & Williams (2005) have created a nearly complete census of the Taurus star-forming region at millimeter wavelengths, allowing us to complete the median fluxes at 850 micron. The median flux we derive is log (νF_ν [erg/s/cm²] ) =  −12.70, with the upper and lower quartiles at  − 13.05 and  − 12.13, respectively. These fluxes are roughly an order of magnitude lower than those included in the original median, but are consistent with the average disk mass in Taurus found by Andrews & Williams (2005).

Appendix B: Herbig median

The Herbig median SED is based on the samples of Acke et al. (2009) and Juhász et al. (2010). Because knowing the dust mass is crucial in measuring the degree of settling in the outer disk, we selected only those sources with a millimeter detection. Following Acke et al. (2009), we also excluded transitional disks, since their mid-infrared SEDs should be explained in the context of gaps and inner holes, rather than dust settling.

The sample consists of 32 sources, and has a complete wavelength coverage in bands, B, V, J, H, K, L, Spitzer IRS longward of 10 micron, IRAS 12, 25 and 60 and millimeter wavelengths. U- and M-band photometry were also included because they lacked only one and two measurements, respectively. The Spitzer spectra were reduced to narrow-band photometry following the same method as in Furlan et al. (2006). At the shortest wavelengths (5.7, 7.1, 8.0, 9.2 μm) the coverage is sometimes incomplete (lacking 5, 5, 3 and 3 measurements resectively), and was supplemented with ISO data where available. We only included the IRAS 60 micron fluxes, because the 12 and 25 micron regions is covered by Spitzer, and results in the same median fluxes. Although all sources have millimeter detections, not all are measured at the same wavelength. We therefore constructed one photometric point at 850 μm, and interpolated fluxes according to νF_ν(mm) = νF_ν(obs)·(λ_obs/850   μm)^-4, averaging if more than one measurement was available.

The median star in our sample is a Herbig A6 star, with a median effective temperature of 8500 K, luminosity of 21 L_⊙ and mass of 2 M_⊙. Because these stars are not in a single star-forming region, we scaled their SEDs to a distance of 140 parsec and their luminosity to the median luminosity before constructing the median. The median including quartiles are given in Table B.1, and are displayed in Figs. 6 and C.1. This approach results in an SED that has a small spread in photometry at stellar wavelengths, and a similar spread as the Taurus median SED at disk wavelengths. Typical features like the 2 micron bump and silicate features are clearly visible.

Appendix C: The inner regions of Herbig stars

As mentioned in Sect. 3.2, the inner regions of Herbig stars deserve special attention in disk modeling. The strong 2 micron bump is a prominent feature in most Herbig Ae and Be stars (Meeus et al. 2001), and has been explained as the fully illuminated surface of a disk truncated at the dust evaporation radius by Dullemond et al. (2001). Because of the large fraction of reprocessed light in the near-infrared, this would require the rim to be “puffed up”, a natural cause of its high temperature down to the midplane (Natta et al. 2001).

However, Vinković et al. (2006) showed that a rim in hydrostatic equilibrium does not puff up far enough to explain the near-infrared fluxes. To fit the SED, the inner rim has to be puffed up beyond hydrostatic equilibrium by a factor 2–3 (Verhoeff et al. 2010; Acke et al. 2009). Other explanations have also been proposed, ranging from dust halos (Vinković et al. 2006; Verhoeff et al. 2011) to optically thin gas or dust within the inner rim (Kraus et al. 2008; Tannirkulam et al. 2008a), as well as theoretical motivations for an increased scale height in the inner rim (Thi et al. 2011) or the presence of such a halo (Krijt & Dominik 2011). For the context of this paper, it is important whether or not the added emission in the inner disk influences the geometry and emission of the outer disk, and could affect our estimate of the turbulent mixing strength. A puffed-up inner rim casts a shadow on the outer disk, decreasing the flux at longer wavelengths (Acke et al. 2009), whereas material within the inner rim or a halo does not cast a shadow (Mulders et al. 2010).

Figure C.1 shows the SEDs of these alternative inner disk geometries. When we increased the scale height in the inner rim by a factor 2.5 to fit the near infrared part of the SED (Fig. C.1, solid line), it cast a shadow over almost the entire outer disk, and we were unable to fit the far-infrared flux even with a turbulent mixing strength of α = 10^-2 or higher. We were able fit the far infrared flux by increasing the amount of small grains, which we did by increasing the slope of the grain size distribution to q = 4.0 (dotted line). However, the strong effect of the inner rim shadow remains, supressing the mid-infrared flux. We tried several prescriptions to modify the inner rim structure, such as rounding it off, but they did not provide enough mid infrared flux. The reason why the puffing up the inner rim does not work well in our model – in contrast to e.g. Acke et al. (2009) and Thi et al. (2011) – is that the median SED requires additional flux at 2 micron – whereas previous work focused mainly on the 8 micron region. This requires a hot inner rim, 1500–1700 K, which is consistent with observations and dust evaporation, but which does not produce enough MIR emission.

Alternatively, dust or gas within the inner rim can also contribute to the near-infrared flux. Additional components originating in that region have been observed for a number of Herbig stars, such as MWC 147 (Kraus et al. 2008; Bagnoli et al. 2010), MWC 275 and AB Aurigae (Tannirkulam et al. 2008a,b), HD 163296 (Benisty et al. 2010) and HR599 (Benisty et al. 2011). Material within the inner rim can increase the near-infrared flux if it generates extra luminosity from accretion (Kraus et al. 2008) or if it captures additional radiation that would otherwise not reach the disk.

For the second option, the material needs to extend sufficiently close to the star, such that starlight is captured that crosses the midplane and would normally not be reprocessed by the disk. We modeled this by extending the largest grain component in our model down to two stellar radii, or 0.02 AU (Fig. C.1, dashed line). Because even the largest grains in our model at this location become extremely hot, we were unable to fit the shape of the near-infrared SED this way, but it does capture roughly sufficient energy at shorter wavelengths. Because this solution to the near-infrared flux problem is similar to a halo in the sense that it captures additional energy, we did not explore it in more detail. We limited ourselves to modeling a halo, where we note that this could also represent a solution with material within the inner disk.

All Tables

All Figures