Modelling and interpretation of SEDs

Circumstellar disks are mostly detected by larger continuum fluxes in the infrared to mm spectral regions as compared to naked stars (a flux excess). The analysis of the Spectral Energy Distribution (SED) of that flux excess was crucial for the development of the first theories about protoplanetary disks, and even nowadays, it is still one of the major tools to physically characterise the disks in terms of their mass, inner holes and gaps, vertical extension & shape, dust properties, and evolutionary state. In this chapter, we will review some of the early simple theories, show some examples, discuss the influence of typical disk shape and dust size parameters in modern SED analysis, and discuss how degenerate the results can be.


Introduction
First evidence for the existence of protoplanetary disks was derived from a "nebulosity" detected in optical images around young intermediate mass stars (Herbig 1960).However, at that time, the concept of protoplanetary disks wasn't yet established, and astrophysicists thought about star formation with a more simplistic spherical infall model, thus expecting pre-main sequence stars to be enshrouded in a spherical envelope of dust and gas.
With the launch of the IRAS satellite in 1983 (InfraRed Astronomical Satellite), it became evident that most pre-main sequence stars possess a strong far infrared (far-IR) excess.Spherical infall models have difficulties to explain both at the same time, the optical light received from the central objects as well as the strong IR excess, so the necessity to deviate from spherical geometry was becoming obvious.However, the IRAS telescope did not have sufficient spatial resolution to image the disks, and so the existence of protoplanetary disks was still a matter of debate.First images of mm dust emission in the late 1980s (e.g.Smith & Terrile 1984), and double-peaked line profiles at mm wavelength showed the extended size and velocity profiles expected from a disk in Keplerian rotation (see Dominik 2015).In the 1990s, NASA's Hubble Space Telescope (HST) clearly imaged the disks in reflected starlight, and in silhouette against the illuminated background of the Orion nebula (O'Dell & Wen 1994).These were the first times astronomers directly saw the protoplanetary disks where planets form.
To date, the observation of the Spectral Energy Distribution (SED) of the light emitted by premain sequence stars and their disks is by far the most frequently used astronomical technique to discover new, and to characterise the structure of protoplanetary disks, for example in terms of their dust properties, their disk masses, gaps and holes, flared or flat structure, etc..The results are often discussed in the context of possible evolutionary scenarios (e.g.Meeus et al. 2001).However, SED analysis is a highly degenerate problem, as we will see in this chapter, so the conclusions from fitting SEDs alone should be taken with a pinch of salt.To improve the reliability of SED analysis, other observations, such as visibilities, images, the spectral features of Polycyclic Aromatic Hydrocarbon (PAH) molecules and gas lines, must be folded in to break these degeneracies (e.g.Maaskant et al. 2014).Today's SED analysis involve advanced, at least 2D numerical radiative transfer techniques, which self-consistently determine the dust and PAH temperature structure in the disk.This wasn't always possible in the past, not only due to lacking computer power, but also because of only recent advances in the development of more powerful numerical algorithms, for example Monte-Carlo techniques.

Example SED
Figure 1 shows the SED of the Herbig Ae star HD 163296, fitted by a simple continuous disk model between 0.45 AU and 700 AU (Tilling et al. 2012).The SED can be subdivided into about 5 different spectral regions, where different physical emission mechanisms are at work, which we will study in this chapter: 1.The region λ ∼ < 1 μm (about ∼ < 2 μm for T Tauri stars) is dominated by direct star light, unless the disk is seen almost edge-on.At these wavelengths, the disk is simply too cold to make any significant contributions to the observed fluxes by thermal emission.Therefore, the analysis of the light in this spectral region does not contain much information about the disk, it is rather useful to characterise the central star.However, there is a certain contribution of the disk in form of starlight scattered on the surface of the disk.With high-resolution chronographs, taken at UV to near-IR wavelengths, it is possible to blend out the dominating star light, and so to image the residual, the so-called scattered light images, from which we can extract certain information about the structure of the disk, in particular about the flaring and radial extension.2. The near-infrared region (near-IR) 1 μm ∼ < λ ∼ < 5 μm is often featured by a strong excess with respect to the star, caused by thermal emission from the hottest grains in the disk, in particular from the inner rim.This disk region is still quite poorly understood.Adding some kind of blackbody emission with temperatures ∼ 1500 K is usually sufficient to fit that part of the SED.In the depicted case, Tilling et al. (2012) fitted the near-IR by assuming a tall inner rim, but others (e.g.Verhoeff et al. 2011) fit this region by assuming the existence of a spherically symmetric, optically thin "halo" of small grains close to the star.3. The silicate emission region 5 μm ∼ < λ ∼ < 30 μm often shows strong silicate dust emission features at 10 μm and 20 μm (the leading Si-O stretching and bending modes of the solid lattice), beside other spectroscopic fingerprints of crystalline silicates and PAHs.This is the spectral region most directly influenced by dust opacities, and beautifully covered by the ISO SWS and Spitzer IRS instruments.Accordingly, this region is most promising to retrieve information about the composition, size, and mineralogy of the dust grains (e.g.van Boekel et al. 2005).However, the photons detected at these wavelengths are almost entirely emitted from dust grains in the upper regions of the inner disk, and it is questionable whether the fitting results about the composition and size of the gains can be generalized to the whole disk.Protoplanetary disks are usually massively optically thick at these wavelengths.So why does the opacity play a role at all?We will discuss this question in Sect.3.3.
4. The multi-color region is characterised by a simple power-law SED between 30 μm ∼ < λ ∼ < 300 μm, and can be understood by optically thick black-body emission from a disk with a certain radial temperature gradient, see Sect.3.1.5. Eventually, in the sub-mm to cm region, the disks become optically thin, and the SED starts to decrease with a certain slope, mostly given by the slope of the dust absorption opacity, see Sect.3.2

Analytical Approximations
In the following, we will study some simple analytical theories about SEDs, valid in certain spectral regions only, if at all, but important for basic understanding.The sections are ordered according to complexity and time at which these ideas were developed.

Multi-color black-body SEDs
The most basic scenario is a flat circumstellar disk seen under inclination angle i, see Fig. 2, where i = 0 o means face-on, and i = 90 o means edge-on.Let r and dr denote the radius and the width of an annulus of the disk which, as seen from the Earth, occupies solid angle where d is the distance.We assume that the vertical structure in the disk is isothermal1 with T (r) being the temperature of the dust grains at radius r.We furthermore assume that the disk is vertically optically thick, in which case the emitted spectral intensity is given by I ν = B ν T (r) , where B ν (T ) is the Planck function.The observed flux from that annulus is F ν = I ν dΩ, i.e. the total flux is given by By means of Eq. ( 2), the SED can be easily calculated if the radial temperature profile T (r) in the disk is known.Let's assume that T (r) is given by a power-law To roughly determine the power-law index q, let's put T (R in = 0.1 AU) ≈ 1500 K, and estimate the temperature at R out ≈ 200 AU.From observations we know that the disk should be as cold as T min ≈ (5 ... 15) K at the outer radius, which provides an estimate2 of q ≈ (0.6 ... 0.75).For such q-values, we find the following situation in the multi-color spectral region as depicted in Fig. 3.At some wavelength λ, the flux is dominated by a certain disk annulus of radius r, where the temperature T (r) is such that the maximum of its local Planck function coincides with λ.The disk regions smaller than r are unimportant, because their solid angles ∝ r dr ∝ r 2 are comparably small, and the more distant disk regions are unimportant because they are too cool to emit at λ.In other words, considering a spectral sweep of increasing λ, in the multi-color region, there is aways a much bigger outer disk region to come which is yet too cold to contribute to the flux at λ, see Fig. 3.One can use these considerations to derive the SED-slope in the multi-color region, see (e.g.Beckwith et al. 1990).The location of the maximum of the Planck function is at ν ∝ T (Wien's displacement law), the maximum of the Planck function is max{νB The multi-color region ends at the "turn-over" wavelength λ turnover , where the coolest temperature in the disk T min is reached, i.e.where the maximum of the coolest Planck function is located at λ turnover .We can use Wien's displacement law to estimate λ turnover , assuming again T min ≈ (5...15) K λ turnover = 2900 μm T min /K ≈ (200 ... 600) μm . (5) Beyond the turn-over wavelength λ > λ turnover , all disk regions contribute to the spectral flux νF ν with the descending r.h.s. of their Planck functions.Eventually, for λ λ turnover , all contributing Planck functions scale as in the Rayleigh-Jeans limit, so we find Equations ( 4), ( 5) and ( 6) give an astonishingly realistic first picture of disk SEDs.In fact, observed SEDs tend to kink down at a few 100 μm (see, for example, Fig 1).And, assuming q ≈ (0.6 ... 0.75), the predicted SED slope in the multi-color region is 2/q − 4 ≈ − (0.66 ... 1.33) which resembles real far-IR observations quite well.Note, however, that the SEDs would be flat for q = 0.5, and even increasing for smaller q, i.e. for values we are nowadays used to find in modern disk simulations.

Low optical depths in the mm-cm region
One modification of the above derived picture of SEDs is needed at mm − cm wavelengths, where the disk becomes more and more transparent, so our assumption of optically thick black-body emission will break down at some wavelength λ thin .Unfortunately, this just happens at about the same wavelength as predicted by λ turnover .Thus, from observations, it is not clear whether the SED kinks down Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰ because of lacking opacities, or because even the coolest disk parts approach the Rayleigh-Jeans limit.
An improved SED-model can be formulated as where τ ν (r) = κ abs ν (r)Σ(r) is the vertical optical depth across the disk, Σ(r) [g/cm 2 ] is the dust column density and κ abs ν (r) [cm 2 /g] is the dust absorption opacity per dust mass at radius r.Let's consider the case λ max{λ turnover , λ thin } for a moment, where the Rayleigh-Jeans limit B ν (T ) ≈ 2ν 2 kT/c 2 is valid.The total mass of dust in the disk can be written as ν is furthermore independent of position in the disk, we can write where T = T (r) Σ(r) r dr Σ(r) r dr is the average dust temperature in the disk.Modern disk simulations result in values T ≈ (10 ... 50) K, depending on stellar luminosity, outer radius, opacities, and shape of the disk.The mm − cm dust absorption opacities are usually featureless, and can be conveniently fitted with a power-law like κ abs ν ∝ λ −β , with β having typical values between 2.0 (small non-conducting grains) to about 1.0 or even lower (very large particles, or elongated particles of conducting materials, see Min 2015), but the value of κ abs ν at a particular wavelength is highly uncertain.
Equation ( 9) is nevertheless widely used in the literature (e.g.Andrews & Williams 2007) to estimate the total dust masses of disks (and total gas masses, assuming M disk = 100 × M dust ) from the observed mm-fluxes, assuming to know the mm dust opacity.It is also not clear whether the condition λ max{λ turnover , λ thin } is already satisfied at λ = 1.3 mm in all cases.Beckwith et al. (1990) improved Eq. ( 9) by introducing a "Δ-correction" for optical depths effects, i.e. for wavelengths λ ∼ < λ thin .The derivation of Beckwith et al. (1990), however, still requires λ λ turnover , and assumes a vertically isothermal disk with radial power-laws for T (r) and Σ(r).From an observational point of view, the method is difficult to apply properly, because it requires the knowledge of the radius at which the disk becomes optically thin at the selected wavelength.

Flared disks, and the surface layer
The motivation to go beyond the SED models discussed above is at least twofold.(i) Physically, only the upper "surface" disk layers can directly absorb star light, see Fig. 5.A simple consideration of the energy budget of surface layer and disk interior shows that these regions must have vastly different temperatures.(ii) Observationally, the (5 − 30) μm region of SEDs often shows strong 10 μm and 20 μm silicate emission features, although the disks are certainly massively optically thick in that spectral region.If the disk was truly vertically isotherm, the opacity wouldn't have any effect on the vertical emission.We need a warm, optically thin layer of dust over a cold thick midplane to explain the observed spectral emission features!Chiang & Goldreich (1997) and Chiang et al. (2001) developed a two-layer disk model which accounts for the physically different conditions of surface layer and disk interior, see also Dominik et al. 2003;Dullemond & Dominik 2004;Dullemond et al. 2001).The optical depth for photons emitted by the star is given by the integral over the extinction opacity κ ext (r , z ) [1/cm] along a radial EPJ Web of Conferences Figure 5.Only the blue shaded disk region, the "surface layer", is directly illuminated by the star.The stellar photons are absorbed in the surface layer, which re-emits half of that energy away from the disk, and half of the energy downward, heating the disk interior.That energy is thermalised by the disk interior and eventually escapes it in form of long wavelength photons toward the observer.Studying the energy budget of midplane and surface layer reveals that the surface layer must be much warmer that the disk interior.ray where, in cylinder coordinates, z/r = z /r = tan(α) = const.Thus, at any position in the disk (r, z), the radial optical depth is given by The spectral index is deliberately omitted here, as we consider implicitly a "typical wavelength" for stellar photons around the stellar maximum, ∼ 400 nm for a Herbig Ae star, and ∼ 1 μm for a T Tauri star.We can now implicitly define the "surface height" of the disk h s (r) by At any given radial distance r, there will be one disk height z = h s (r) where the radial optical depth approaches unity, i.e.where the majority of stellar photons are absorbed or scattered fir the first time, which we identify as the location of the "surface layer".According to the definition of the radial optical depth (Eq.10), h s (r)/r is a strictly monotonic increasing function.Anticipating a power-law behaviour, we write where ξ > 0 is the flaring parameter.The basic idea of Chiang & Goldreich (1997) is to study the energy balance of a vertical column in the disk, between radii r 1 and r 2 , as irradiated indirectly by the surface layer, see Fig. 6.The angle α, as sketched in Fig. 6, is given by tan α = h s (r)/r, and the solid angle of the flared surface of the disk between r 1 and r 2 , as seen from the star, is where Δr = r 2 − r 1 .The photon energy received by that piece of surface per time is L ΔΩ/(4π).
Assuming that half of that energy is re-emitted or scattered toward the midplane, the heating of the column is Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰ 00007-p.7 Since the interior of the disk is assumed to be optically thick, radiative transfer will be diffusive, the column will be vertically roughly isothermal, and the considered column will cool via black-body emission from its surface where σ is the Stefan-Boltzmann constant.Balancing L heat and L cool , we get the temperature of the disk interior where we have used the Stefan Boltzmann law L = 4π R 2 T 4 , with R and T being the stellar radius and effective temperature.Equation ( 16) is a remarkable result.It states that the temperature of the disk interior is triggered by the degree of surface flaring.In a self-shadowed disk, we have virtually ξ = 0, which would imply T = 0.This must be wrong, of course, showing the shortcomings of the assumption made above.In reality, there is also multiple scattering and diffusive radial radiation transport, both of which will prevent T = 0. Using h s /r ∝ r ξ (Eq.12), Eq. ( 16) predicts a temperature power-law T ∝ r −q with index q = 0.5−ξ/4, Let us now turn to the surface layer itself.The condition of dust radiative equilibrium states an equality between the amount of absorbed and emitted radiation energy per time where J ν is the mean intensity.The surface layer is approximately optically thin.Ignoring the flux emitted upward by the disk below, the mean intensity in the surface J surf ν is approximately given by a narrow pencil of direct stellar radiation J surf ν = 1 4π F ν , where the direct stellar flux is , and L ν is the specific luminosity of the star, approximated by a Planckian of effective temperature.Introducing the Planck mean dust absorption opacity EPJ Web of Conferences 00007-p.8 and utilising B ν (T ) dν = σ π T 4 , the temperature structure in the surface layer is found to be Equation ( 19) is an implicit equation for the determination of T surf (r).For very large grains, for example, which have flat opacities (κ abs ν ≈ const), the Planck-mean opacities would cancel, and Eq. ( 19) would result in T = T √ R /(2r), i.e. we get a q = 0.5 radial power-law.However, small grains are warmer, because their opacities generally decrease with increasing wavelength, such that κ Pl abs (T ) κ Pl abs (T surf ), and this effect tends to decrease the value of q.For example, for a power-law opacity κ abs ν ∝ λ −p , the result is q = 2/(4 + p) ≈ 0.5− p/8.By combining Eqs. ( 16) and ( 19), the temperature contrast between surface and interior is By order of magnitude, at r ∼ 100 AU, we have ξ ∼ 0.3, h s /r ∼ 0.3, κ Pl abs (T )/κ Pl abs (T surf ) ∼ 100, predicting a temperature contrast between surface and interior of about a factor of 5.This is a remarkable result, however, modern disk models result in even larger temperature contrasts, of up to a factor of 10, see Fig. 7, because additional physical effects are taken into account, in particular dust settling.Chiang & Goldreich (1997) noticed that Eq. ( 16) offers a way to self-consistently determine the flaring index ξ when combined to the equation of vertical hydrostatic equilibrium in the nearly isother-Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰ 00007-p.9mal disk interior.To do that, we assume3 that the unknown surface height h s is proportional to the hydrostatic scale height where χ is a free fit parameter ∼ 1 ... 6.Assuming the disk to be in Keplerian rotation with orbital frequency Ω= GM /r 3 and using c T = kT/μ = H Ω, the scale height is given by where c T is the isothermal sound speed, μ ≈ 2.3 amu is the mean molecular weight and M is the stellar mass.Taking Eq. ( 22) to the 8 th power, and then substituting T 4 by the central term in Eq. ( 16), and using Eq. ( 21), we find which shows that the scale height should vary with a power-law as H ∝ r 9/7 , and hence the flaring parameter (introduced as h s /r ∝ r ξ , Eq. 12) is ξ = 2/7 ≈ 0.29.Using χ ≈ 2.5, and values for a T Tauri star as M = 0.7 M and L = 1 L , Eq. ( 23) predicts H ≈ 10 AU and h s ≈ 25 AU at radius r = 100 AU.In view of modern disk models, these predictions for the scale heights are excellent, but ξ = 0.29 seems a bit too high.The resulting SEDs would be featured by an enormous far-IR excess, increasing in the multi-color region (see Fig. 7 of Dullemond et al. 2001), which may be appropriate for some extremely flared Herbig Ae/Be disks like AB Aur or HD 100546, but typical T Tauri stars have decreasing SEDs in the multi-color region, and their disks seem to be less flared than predicted by Eq. ( 23), with ξ ≈ 0.05 ... 0.2 being more typical fit results.Dust settling might be the key to understand these discrepancies.Dullemond et al. (2001) noted that Eq. ( 16) suggests a "flaring instability", see Fig. 8.If the disk is flared, it intersects more star light, thereby warming up and extending higher up, which again amplifies the flaring, etc.. On the other hand side, if the flaring would stop for a reason, the irradiation would be reduced, the disk would get cooler, and would start to shrink vertically, reducing the flaring, and so on.This could produce a vertically collapsed, very cold, and self-shadowed disk.Dust settling could trigger this transformation.
The above discussed "passive" disk models do not take into account the effect of viscous heating due to frictional processes involved in the mass accretion process.These effects are important, in particular, for classical T Tauri stars with high mass accretion rate.Viscous heating can change the temperature structure in the inner disk regions substantially and can provide an additional near-IR excess.D' Alessio et al. (1999Alessio et al. ( , 1998) ) have developed disk models which take these effects into account in the frame and basic approach of the Chiang & Goldreich (1997) models.

The inner rim
The inner rim of a protoplanetary disk adds an important feature to the SED: the near-IR excess around (1 − 5) μm.In case of Herbig Ae/Be stars, this excess often simply looks like a Planckian of a single temperature of order 1500 K (see, for example, Figure 2 in Dullemond & Monnier 2010).This suggests that, in the near-IR, we see dust grains at a certain temperature, probably just at the edge of their thermal stability -glowing particles which define the position of, and are shaping, the inner rim.Unfortunately, SED analysis cannot provide much more information, in particular about where these particles are.Near-IR and mid-IR interferometers (e.g.PIONIER) have just started to achieve a performance (coverage in spatial frequency domain and signal/noise) that will allow us to deduce more information about the shape of the inner rim in the future.
Therefore, the physical structure of the inner rims of protoplanetary disks is still poorly understood, and is not much more than a theoretical concept so far, see for example Fig. 9. Consequently, SED-modellers are reluctant to assume anything complicated here.The most popular model, by far, is a simple sharp inner edge, a discontinuity in surface density at r = R in , beyond which all dust grains are assumed to be stable, ignoring the possible existence of dust-free gas inside of R in , because it is not required to explain the SED.
Understanding the near-IR excess is principally easy.Just read off the surface height h s from (Eq. 11) at a suitable small radius, e.g.r 0 = 1 AU, compute the solid angle which the inner rim occupies as seen from the star via Ω inner = 2 × 2π sin arctan(h s /r 0 ) , and compute the near-IR excess as where a is the dust albedo at a typical wavelength for the stellar irradiation, ∼ 400 nm for a Herbig Ae star, and ∼ 1 μm for a T Tauri star.Equation ( 24) assumes that the direct radiation absorbed by the inner rim is thermalised and re-emitted in the near-IR.We can estimate the inner disk radius R in , where the dust reaches its sublimation temperature, by balancing the incoming absorbed stellar flux inner rim position is where the r.h.s.formula is valid for albedo a ≈ 0.2 and T sub ≈ 1500 K.In practice, Eq. ( 25) leads to slightly too small inner radii, because the inner rim is not solely irradiated by the star, but also from the far side of the inner rim itself, see review by Dullemond & Monnier (2010).
There is one interesting problem that remains with the near-IR excess.SED models which are based on the assumption of vertical hydrostatic equilibrium tend to under-predict the near-IR excess, in particular for T Tauri disks.The spectral shape of the near-IR excess produced by the model is often OK, but it is too faint: hydrostatic inner rims are too flat.Hydrostatic models, despite their puffed-up inner rims, tend to produce too small h s to achieve the observed level of near-IR excess.
This mismatch has inspired modellers to search for additional/alternative ways to explain the near-IR excess, other than "tweaking" the model, e.g. by lowering the gravity at the inner rim.One idea was presented by Thi et al. (2011), who suggested to use the gas temperature instead of the dust temperature (and the properly calculated mean molecular weight -the gas at the inner rim is not necessarily H 2 -rich) to calculate the vertical structure of the inner rim, which then results to be 2-3 times taller.An even more radical solution of the problem is simply to prescribe the shape of the disk entirely, such that h s can be tuned directly (Woitke et al. 2011).Also, using dust materials with very low albedo at the inner rim, like amorphous carbon, helps to absorb more stellar radiation, which then produces a stronger near-IR excess (Carmona et al. 2014).Another idea assumes the presence of hot grains close to the star in a spherical fashion, entirely disconnected from the disk, to form a "dust halo" (e.g.Maaskant et al. 2014;Verhoeff et al. 2011).One way to interpret this "halo" is to assume that we see small grains entrained in disk winds close to the star, or could we have larger pebbles on off-plane, possibly eccentric orbits?
The uncertain physical structure of the inner disk is hence connected to many other unknowns in this scientific field, such as hydrodynamics and magnetic fields.Unfortunately, improving the disk models for a better treatment of the inner rim will mean a lot of cumbersome numerical efforts to obtain converged results (see, e.g., Kama et al. 2009).But some physical effects connected to the gassolid phase-transition seem so well-understood and robust that future models should aim at including them, because the structure of the inner rim does affect all other disk regions (compare Fig. 9): • The temperatures at the directly irradiated inner rim are much higher that in the shadowed regions behind it: the inner rim should be puffed-up (Dullemond et al. 2001).• The sublimation temperature of any condensate is pressure-dependent.With increasing height over the midplane (decreasing pressure), there must come a point where any condensate becomes thermodynamically unstable: the (dust-) shape of the inner rim should be rounded (Kama et al. 2009).• Dust grains in disks are not made of a single solid material.Different materials evaporate at different temperatures: The inner rims of protoplanetary disks should have a layered dust structure, similar as in brown dwarf atmospheres or AGB star winds (Woitke 2006).• There must be a connection between the dusty inner rim of Herbig Ae/Be stars (located at R in ∼ 0.5 AU) and the region where magnetospheric accretion can work (of the order of a few stellar radii ∼ 0.05 AU).Mass must flow from the dusty disk through this interface, which cannot contain dust: There must be gas inside the dust inner rim of Herbig Ae/Be disks.

00007-p.12
Table 1.Estimated stellar and disk properties as function of spectral type for stellar age = 2 Myr. SpTyp

Disk properties of function of spectral type
In Table 1 we list some expected disk properties as function of spectral type at a fixed age, here 2 Myrs, as resultant from the analytical approximations presented so far.We have used the pre-main sequence stellar evolutionary models of Siess et al. (2000) to determine effective temperatures, stellar masses and luminosities as function of spectral type for 2 Myrs old stars.We then roughly set the disk mass to be M disk = 0.01 M , use Eq. ( 25) to determine the inner radius, and Eq. ( 23), with ξ = 2/7 and χ = 2.5, to determine the scale height.Results will have large uncertainties because of the assumptions involved.Looking at the distribution of stellar masses in Table 1 shows that, for an age of 2 Myrs, most stars should be either M-type or K-type T Tauri stars, with very little G and F stars, and only a few high-mass young Herbig Be stars.The corresponding SEDs are computed from the disk models outlined in Sect. 5 and shown in Fig. 10.

Modern disk models
In modern (passive) disk models, it is no more necessary to make any of the approximations described in Sect.3, we just solve the equation of local dust radiative equilibrium (Eq.17) together with the 3D radiative transfer equation for an axisymmetric disk configuration.It is hence sufficient to simply prescribe (i) the stellar and interstellar irradiation and (ii) the spatial structure of gas, dust and opacity in the disk.Anything else is a result of the model, in particular the dust temperature structure T (r, z), and all predicted continuum observations like the SED, visibilities, and images.Using highly efficient new computational algorithms, in particular Monte Carlo techniques (see Pinte 2015), combined with nowadays computer power, it is possible to compute one complete 2D disk model within a few tens of CPU minutes.
The new Monte-Carlo codes include MCFOST (Pinte et al. 2006) and MCMax (Min et al. 2009), whereas ProDiMo (Woitke et al. 2009) uses a deterministic ray-based method, which is less powerful.The codes have been benchmarked against each other and against other codes (Pinte et al. 2009), showing a reassuring degree of agreement despite large optical depths, in particular with regard to the SED computations.What modern disk models distinguish from each other is rather the way how the disk structure is treated, which dust sizes and opacities are used, and whether or not, and how, dust settling is taken into account.In the DIANA consortium, see http://www.diana-project.com,we have agreed on a "minimum set" of common assumptions and approximations, aiming at setting new standards for disk modelling.These assumptions shall be quickly outlined here, without discussion, to introduce a number of parameters that we need for the discussion later on.
Disk mass and column density structure: The gas column density structure Σ(r) [g/cm 2 ] is given by a radial power-law with index , modified by an exponential tapering off where R tap is the tapering-off radius.The default choice for the exponent is γ = 2 − (self-similar solution).Radial integration over Eq. ( 26), from R in to R out , results in the total disk gas mass M disk , which is used to fix the proportionality constant in Eq. ( 26).The inner rim is sharp and positioned at R in .The outer radius R out R tap is chosen large enough such that Σ(R out ) becomes vanishingly small.

Vertical gas stratification:
We assume a Gaussian vertical gas distribution with a simple parametric prescription of the scale height as function of radius where H is the gas scale height, H 0 a reference value thereof at radius r 0 , and β is the flaring exponent.Vertical integration of ρ(r, z) results in Σ(r) which is used to fix the proportionality constant in Eq. ( 27).
Dust size distribution: We assume a power-law dust size distribution function f 0 (a) = dn 0 (a)/da [cm −4 ] as function of particle radius a [cm] as Equation ( 28) prescribes the dust size distribution function "before settling".Since settling only moves the grains vertically, but doesn't change the total number of dust grains in a complete vertical column, we have f 0 (a, r) = f (a, r, z) dz / dz.Therefore, the local dust size distribution f (a, r, z) deviates from f 0 (a, r).The dust mass density [g/cm 3 ], before settling, is given by 4π 3 ρ d a max a min f 0 (a) a 3 da, amorph.carbon 15% porosity 25% where ρ d is the assumed dust material density, and must equal the gas density ρ times the assumed dust/gas mass ratio, which is used to fix the proportionality constant in Eq. ( 28).

EPJ Web of Conferences
Dust settling: Dust settling is included according to Dubrulle et al. (1995), assuming an equilibrium between upward turbulent mixing and downward gravitational settling.The result is a size and density-dependent reduction of the dust scale heights H(r, a) with respect to the gas scale height H(r), where Ω is the Keplerian orbital frequency, γ 0 ≈ 2 for compressible turbulence, and τ f is the frictional timescale, ρ mid the midplane gas density, and c T is the midplane sound speed.To avoid iterations involving the midplane temperature, we use c T = H Ω where H is gas scale height from Eq. ( 27).α settle is the dimensionless viscosity parameter describing the strength of the turbulent mixing.The l.h.s. of Eq. ( 29) is smoothly limited to a maximum value of one.Technically, in every disk column, Eq. ( 29) is computed for a number of (about 100) dust size bins.Starting from the unsettled dust size distribution, the dust particles in each bin are re-distributed in z-direction accordingly, building up a numerical representation of the local dust size distribution function in every point in the disk f (a, r, z).

Dust opacities:
We assume an effective mix of laboratory silicates and amorphous carbon, with 25% porosity, and use a distribution of hollow spheres (DHS) to simulate shape deviations of the dust particles from spherical symmetry.The details are explained in Min (2015).For every size bin, the dust opacities can be pre-computed, for example the absorption efficiency Q abs (λ, a).The local dust opacities are finally computed, at every point in the disk, by integrated over the local dust size distribution function, for example κ abs ν (r, z) = A Fortran-90 package to compute the DIANA standard dust opacities is available at http://www.diana-project.com/data-results-downloads.
Altogether, a single-zone disk model is fully characterised by the parameters listed in Table 2, not mentioning all stellar and other irradiation parameters here, and disregarding further gas parameters not needed for the computation of SEDs.
The values listed in Table 2 define our "reference disk model".Concerning the star, we consider a "typical" class II T Tauri star with effective temperature T = 4000 K and stellar luminosity L = 1 L .

Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰
These values correspond to spectral type K7, a stellar mass of M = 0.7 M and an age of about 1.6 Myrs according to the pre-main sequence stellar evolutionary models of Siess et al. (2000).The star is assumed to have a power-law UV excess with total luminosity L UV = 0.01 L , which is added to the stellar spectrum in form of a power-law L λ ∝ λ p UV with p UV = 1.3.The reference model has been designed to predict continuum and line observations that roughly resemble the observations of real class II T Tauri stars, in particular • a near-IR excess (2 to 7 μm) of about 0.11 L , • clearly visible silicate dust emission features around 10 μm and 20 μm, • a descending SED-slope in the multi-color region, • a 1.3 mm flux of about 70 mJy at 140 pc, and

Effects of Parameters
We will now study the influence of some of the introduced disk and dust parameters on the calculated SED, to investigate how far disk properties can be determined by SED analysis.In Fig. 11, we plot the computed SED of the reference model (black line, identical in all plots).Each of the 9 part figures show the influence of a single model parameter, where in each case, two additional models have been calculated, usually one with a larger value and one with a smaller value of that parameter, as compared to the reference model.The resulting range of SED results are highlighted by the green shaded areas.Further details of the model results are shown in Table 3.
The dust mass M dust basically shifts the SED up and down at long wavelength, where the disk is predominantly optically thin, according to Eq. ( 9).Its influence diminishes at λ ∼ < 200 μm where we enter the multi-color region in the reference model.However, a change of M dust still produces noticeable changes at wavelengths as short as ∼ 20 μm.The latter can be understood in terms of the relation between H and the surface height h s .For more massive disks, h s lies higher in the disk albeit H remains constant.This refers to the parameter χ as introduced in Eq. ( 21), which should be mass-dependent.
The reference scale height H 0 affects the SED at all shorter wavelength λ ∼ < 200 μm in a similar fashion, but not the Rayleigh-Jeans tail of the SED.Its influence can be understood by Eq. ( 16) which states that T 4 ∝ h s .Larger scale heights mean to intersect more star light, producing a warmer disk interior which re-emits more thermal radiation from all optically thick disk regions.The "+"/"-" models refer to an increase/decrease of the indicated parameter to the value annotated in Fig. 11.The IR excess is measured between 2 μm and 7 μm after subtraction of the naked stellar flux.The 10 μm emission amplitude is computed as F ν (9.6 μm)/ 1 2 F ν (6.8 μm) + F ν (13.1 μm) .The mm and cm slopes, α SED = −Δ log F ν /Δ log λ, are measured between (0.85 − 1.3) mm, and (5 − 10) mm, respectively.
The flaring index β tends to rotate the SED around a point at λ ≈ 20 μm in this case.Large β values mean that we have a flared disk with a low inner rim but with tall outer regions, which produce less near-IR but more far-IR excess.Small β lead to a "self-shadowed" disk structure with a very cold disk in the outer parts.
The inner radius R in regulates the maximum temperature of the dust grains at the inner rim, see Eq. ( 25).Larger R in therefore result in a deficit of near-IR emission ("transitional disks").However, according to Eq. ( 24), the amount of excess luminosity produced by the disk stays the same since the same amount of stellar radiation is still being absorbed.For large R in , the luminosity excess merely shifts from the near-IR to the mid-IR region, and beyond.Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰ 00007-p.17Dust settling, see parameter α settle , almost exclusively affects the long wavelength parts of the SED (λ ∼ > 20 μm).According to the Dubrulle prescription (see Eq. 29), dust settling is much more effective at large radii where the densities are low, in which case the dust grains cannot be easily dragged along turbulent gas motions.Consequently, the outer disk parts become flat "as seen in dust", although the gas still extends high up.When looking at the effects of β and α settle in Fig. 11, there seems to be no way to disentangle flaring from settling by SED analysis of observations at long wavelengths.In order to discriminate these two physically entirely different mechanisms, we need to observe both, gas and dust.However, the amplitude of the silicate emission features at 10 μm and 20 μm is also affected by settling, see Table 3.For a well-mixed dust/gas mixture, the mm-grains cover all spectral features produced by the small grains with their flat, grayish opacity, which tends to wash out the 10 μm and 20 μm silicate emission features in the SED.Dust settling removes the large grains from the disk surface which seems necessary in many cases to produce the observed amplitude of the silicate emission features.
The maximum dust particle size a max has a similar influence as M dust at long wavelengths.Increasing a max effectively means to put more and more dust mass into very large particles which have almost no opacity at shorter wavelengths.However, beyond 1 mm in this model, where the largest particles do provide the dominating opacities, the SED starts to change slope depending on the value of a max .This happens at about dust size parameter x = 2πa max /λ ∼ < 1, i.e. for λ ∼ > 2π a max .
The dust size distribution power-law index a pow regulates the mixture of small and large dust particles in the disk.It thereby changes, by order of magnitude, the transparency of the dusty gas at small, in particular UV wavelengths, which changes the UV penetration depths important for the gas modelling.At long wavelengths, a pow controls the opacity slope which translates into a change of the "mm-slope", although the effect only fully develops at cm-wavelengths, see Table 3.
The volume fraction of amorphous carbon has a surprisingly large impact on the SED at all wavelengths.Pure laboratory silicates (both amorphous and crystalline) are almost completely transparent to optical and near-IR radiation.Such particles are efficient scatterers, keeping the stellar radiation out of the disk, but they will hardly absorb it.Therefore, disks made of pure silicates emit much less in the near and far IR, as well as much less mm flux.At mm and cm wavelengths, the conducting amorphous carbon inclusions (if properly treated with deviations from spherical symmetry as in these models) show the "antenna-effect" (see Min 2015) which significantly increases the dust absorption opacities at long wavelength.Large volume fractions of amorphous carbon also tend to suppress the 10 μm and 20 μm silicate emission features, see Table 3.
The surface density power-law index has practically no influence on the SED, the same with the tapering-off radius R tap , the outer radius R out (not depicted, see e.g.Bouy et al. 2008) and with the minimum dust size a min (as long as a min ∼ < 0.5 μm, not depicted).The dependencies of the SED on the latter, non-depicted parameters is less than for .Other observations must be used to determine these parameters, like images, visibilities, or line observations.Figure 12 visualises which spatial parts of the disk are mainly responsible for the observable fluxes at different λ, according to the reference model.At near-IR wavelengths as short as 1 μm, we have extended "emission" in form of scattered star light from the uppermost disk regions.For slightly longer near-IR wavelengths, we mostly see thermal emission from the inner rim and innermost disk regions < 0.3 AU at high altitudes.For wavelengths λ ∼ > 10 μm, the emission region detaches from the inner rim, moves to more distant regions and deeper layers, before it finally falls down at mm-wavelengths to a narrow region centred around the midplane, where most of the mm-grains are concentrated via dust settling.The red box is located almost symmetrically around the midplane, suggesting that the mm-emission is optically thin, but its curved shape shows that, even at λ = 1 mm, there are still some optical depths effects at work.Modellers should be aware of this complication.SED analysis is a good example for the "inverse problem" in astrophysics (Lucy 1994).Tuning the free model parameters to fit the SED is relatively easy, but what have we actually learnt?What we would like to know is, in how far this process can be inverted, starting with the observations, and ending in values for the physically meaningful model parameters with the so-called confidence intervals.

Limits and degeneracies in SED analysis
Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰ In fact, very satisfying SED-fits can often be obtained either by varying the dust opacity parameters or by varying the disk shape parameters and zone setup.This should make us suspicious.Just because it fits doesn't mean it's correct.
Another important conclusion from Figure 11 is that some of the introduced physical parameters have almost no direct impact on any spectral region in the SED, for example the disk outer and tapering-off radii, R out and R tap , the column density power-law index , and the disk inclination i (unless the disk is seen close to edge-on).In order to determine these parameters, we need spatially resolved data, or line data.For some objects, (sub-)mm and cm flux measurements might be entirely missing.In those cases, it seems impossible to determine certain disk properties from the observed SED, such as the disk mass and the maximum dust particle size.Knowing all this, one should be reluctant to treat those parameters as "free" when fitting just SEDs, because the outcome will be more or less random.It can improve the fits for sure, but in spurious ways.
In order to break these degeneracies, we need to include other types of observational data, such as resolved images and visibility data.Furthermore, using "holistic disk models", which predict not only continuum but also line observations, allows us to utilise the large body of information contained in line fluxes, line velocity profiles, and line maps/visibilities (e.g.Carmona et al. 2014).The drawback of this approach is that these models are much more complicated, and that the gas physics involved is much less well-understood as compared to dust radiative transfer physics.In addition, it makes the job of the modeller obviously much harder.A not unlikely result of such efforts could be that we simply cannot understand the complete set of observational data based on the current model.
In order to deal with these degeneracies, Bayesian analysis is the method of choice (see e.  errors, the likelihood of the model θ is defined via the joint probability to fit the data as where p(y) is an unimportant normalisation constant, and the "raw" χ 2 (not reduced by factor 1/N) is defined as Pure χ 2 -minimisation methods simply optimise the model parameters θ until the deviations between model predictions and measurements are minimum, but these methods do not determine the confidence intervals, i.e. the statistical uncertainties of the derived model parameters.Bayesian analysis, in contrast, is designed exactly for that purpose.The starting point of Bayesian analysis is the prior distribution p(θ), which can be any given probability density distribution of where the notation beautifully expresses the need to invert the modelling process.Here, p(θ | y) is called the posterior distribution, i.e. the conditional distribution of the parameters θ after having seen the data y.The essence of Bayesian analysis is that the posterior distribution p(θ | y) is a proper probability distribution in model parameter space, which can be used, for instance, to compute expectation values and standard deviations as where θ k are the most probable parameter values (the posterior means) and σ θ k are the confidence intervals (the posterior standard deviations).These definitions are so-called frequentists definitions (see VanderPlas 2014).There are some subtle differences between the frequentist and the Bayesian way to understand probabilities, in particular if the prior distribution p(θ) is not flat.However, in the most simple case of a uniformly sampled parameter space, we have p(θ) = const and the differences disappear.Figure 14 shows an example for a computed Bayesian posterior distribution p(θ | y), from which one can easily read off the confidence intervals by eye, or according to whatever recipe is favoured.One could, for example, determine the parameter intervals that contain 68% of the probabilities, or the parameter ranges with probabilities above exp(−1/2) of the maximum probability.The second definition is more safe, as it will contain the complete interval if the posterior probability is flat, for instance.
Concerning SED fitting, a painful conclusion from Bayesian analysis is that, in general, the more parameters we allow to vary freely in the model, the wider the confidence intervals become, which simply shows how degenerate SED fitting can be.For the case depicted in Figure 14, nine parameters have been varied to fit the SED of BP Tau, and the derived parameter ranges are almost "generic", i.e. similar to the ranges of parameters considered as "typical" for a whole class of T Tauri stars.Picking just the most probable parameter values θ k does not even guarantee to yield a well-fitting SED model.Just think about Fig. 13 for a moment.
Another practical challenge involved in determining the confidence intervals, e.g. according to Eqs. ( 35) and ( 36), is that we need to sample all relevant parts of the parameter space, i.e. we need to compute a sufficient number of models which fill in all volumina in the parameter space that are relevant for the computation of Eqs. ( 35) and (36).We need not only good fitting models, but we need a large number of badly, but not too badly fitting models.There are well-established numerical algorithms which fill in those parts in parameter space automatically, in particular the MCMC (Monte-Carlo Markov Chain) algorithm, but these techniques have been developed for "models" that take a split-second to run.In our case, we need at least a couple of CPU minutes to run one single disk model (which predicts the SED) or even several CPU hours, if we want to also predict images, visibilities, and gas lines.Therefore, although desirable, running tens of thousands of such disk models, just to get the errorbars, does not seem to be appropriate in all cases.

Figure 1 .
Figure 1.The SED of the Herbig Ae star HD 163296 (Tilling et al. 2012, reproduced with permission c ESO).The blue, red, black and green circles with errorbars are photometric measurement points from different instruments, the green line shows a de-reddened observed UV spectrum, and the orange line shows an ISO/SWS spectrum.The full black line is the disk model, and the thin red line indicates the photospheric model.The SED can be roughly sub-divided into 5 different spectral regions, as discussed in the text.

Figure 2 .Figure 3 .
Figure 2. The spectral flux received from an annulus of a flat disk, seen under inclination i.

Figure 4 .
Figure 4. Steeper SEDs in the mm-region because the dust becomes increasingly optically thin.

Figure 6 .
Figure 6.Stellar irradiation of a flared disk surface

Figure 7 .
Figure 7. Dust temperature structure in a T Tauri star disk from the (1+1) D model of Chiang & Goldreich (1997, c AAS, reproduced with permission) (l.h.s.), and from modern 2D disk models (r.h.s.).Note the effect of the emission from the disk interior on the surface temperature structure at large radii on the right side.There is also interstellar irradiation on the r.h.s., which prevents the midplane dust temperature to decrease further beyond r ∼ > 100 AU.
in 2 σT 4 with the emitted cooling flux F cool = σT 4 sub , where T sub is the dust sublimation temperature.The resulting Summer School "Protoplanetary Disks: Theory and Modeling Meet Observations‰ 00007-p.11

Figure 11 .
Figure11.Effect of dust and disk parameters on model SED at distance 140 pc and inclination 45 • .The thick full black line is the reference model (identical in every part figure), whereas the green shaded area indicates the effect of a single parameter on the SED, where the dashed on dotted lines correspond to the changed parameter values as annotated.Top row: dust mass M dust , scale height H 0 , and flaring exponent β. 2 nd row: inner radius R in , column density power-law index , and dust settling α settle .3 rd row: maximum grain size a max , dust size power-law index a pow , and volume ratio of amorphous carbon amC.The dependency of the SED on the outer and tapering-off radii, R out and R tap (not shown) are less than the one shown for .

Figure 12 .
Figure 12.Spatial regions in the disk, as function of radius r and height z, which emit most of the spectral fluxes seen at different wavelengths (for the reference model, see Table 2, in face-on orientation).The underlying greyshaded contours represent the assumed disk gas density structure, where n H (r, z) is the hydrogen nuclei density [1/cm 3 ].The colored boxes identify the disk regions emitting (and scattering) about 50% of the total flux.The dashed lines indicate vertical optical extinctions A V = 1 and A V = 10.

Figure 11
Figure 11 indicates that certain observed SED characteristics can be explained in different ways, for example by a lack of disk flaring or by dust settling.These "degeneracies" often render it difficult to confidently determine a disk properties by means of SED analysis.Another example of such degeneracies is shown in Fig 13.Known degeneracies include • disk dust mass -mm dust opacity • distance -luminosity • flaring -settling • dust opacities -nearly all disk shape parameters • one zone?two zones?three zones?...

Figure 13 .
Figure 13.Example for degeneracies in SED analysis.The contours show the likelihood of fitting the SED, as function of inner rim radius R in and maximum dust size a max .The SED of this transitional disk can be fitted about equally well if the inner rim is located at about 50 AU and the disk has only ≤ 1 μm dust particles, or, if the inner rim is located at about 15 AU and the disk has up to 1 mm dust particles.However, values around R in ≈ 25 AU seem most likely.Figure kindly provided by C. Pinte.

Figure 14 .
Figure 14.Bayesian analysis, based on 30 000 disk models, to fit the SED of the classical T Tauri star BP Tau.The model parameters have meanings as explained in Sect.5, however "aexp" is a pow , "surf" is − and "strat" is a parameter involved in a different way to prescribe dust settling.Figure kindly provided by C. Pinte.

Table 2 .
Disk model parameters, and values of the reference model.

Table 3 .
Changes of some calculated disk and SED properties.ref.model M disk + H 0 + β− R in + α settle − a max − a pow + amC+