MID-INFRARED VARIABILITY OF THE BINARY SYSTEM CS Cha

We assume that the edge of the circumbinary disk is a vertical cylindrical wall illuminated by the stars. To the best of our knowledge, there are no calculations on the shape of the wall in circumbinary disks. For accretion disks around single stars, it has been proposed that the inner wall is curved, because the sublimation temperature of the dust grains depends on density (Isella & Natta 2005) and on grain size (Tannirkulam et al. 2007). However, in the case of a circumbinary disk, the inner disk truncation is produced by dynamical effects and not dust sublimation. For the case of disks with embedded planets, there are detailed calculations by Varnière et al. (2004) and Crida et al. (2006) of the surface density profile of the gap, that in principle could be used to define a shape of the inner wall. However, the cases discussed in these papers correspond to low-mass companions (planets) and it is not clear how one can scale this to binary systems. In addition, the density profile depends on viscosity, the disk aspect ratio, etc., which are unknown properties for CS Cha.

We assume that the dust in each surface element of the inner edge of the circumbinary disk is in radiative equilibrium with the impinging radiation field from the stars. The flux received at each surface element of the wall is the addition of the contributions of both stars, taking into account the geometrical dilution of the stellar radiative flux, for the different distances of the stars at different orbital positions. Also, we assume that each star is accreting from a small circumstellar accretion disk that receives mass from the circumbinary disk (Günther & Kley 2002). It is usually assumed that for T Tauri stars surrounded by accretion disks, accretion shocks exist at the stellar surface that emit UV and optical radiation. We include the emission of these shocks as a source of irradiation (Muzerolle et al. 2003). The luminosity liberated in these shocks is the accretion luminosity, ∼L_acc, and we assume an effective temperature for the shocks of ∼8000 K (Calvet & Gullbring 1998). Thus, the SED of each star is a combination of their photospheric and shock emission.

The zeroth and first moments of the radiative transfer equation are solved as in Calvet et al. (1991; see also D'Alessio et al. 2005). The irradiation flux interacts with the dust in the wall, a fraction is absorbed and another fraction is scattered. The scattered fraction produces a diffuse radiation field, which is finally absorbed at deeper layers. The transfer equation is integrated along the radial direction (parallel to the disk midplane) for each surface element of the disk inner wall, and a temperature is calculated assuming radiative equilibrium. The temperature is a function of radial optical depth. In this calculation, we assume that the density of the wall at the edge of the circumbinary disk is high enough to allow thermal equilibrium between grains of different sizes and gas. Thus, a unique temperature describes the wall's thermal emission for each wall location.

For the calculation of the temperature on the surface of the wall [T(τ = 0)], we use the optical properties of the dust at this temperature. The implicit equation for this temperature is then solved with an iterative procedure. The opacity evaluated at the wall surface temperature is used to calculate the temperature as a function of τ. Finally, with the temperature radial distribution of each wall pixel, we calculate the emergent intensity, and taking into account the geometry of the wall and occultation effects, we calculate the contribution of the wall to the SED.

The models of gap formation in circumbinary disks show a decrease in mass surface density between an inner hole formed due to gravitational interactions of the binary and the disk (Artymowicz & Lubow 1994) and the outer disk, which suggests the gap is optically thin. Note that planets immersed in the disk are able to increase the size of the hole (Zhu et al. 2011; Dodson-Robinson & Salyk 2011). Observations of circumbinary disks at millimeter wavelengths, which resolve the inner hole (Guilloteau et al. 1999 for GG Tau; Guilloteau et al. 2008 for HH30; Andrews et al. 2010), confirm that it has a low density and low optical depth. This is the case for transitional disks (Hughes et al. 2007 for TW Hya; Hughes et al. 2009 for GM Aur; Andrews et al. 2011). The inner hole might have optically thin dust that contributes to the near- and mid-IR SED (e.g., TW Hya, Akeson et al. 2011; Calvet et al. 2005; Espaillat et al. 2007a, 2007b, 2010).

We consider that the distribution of material inside the hole is not axisymmetric (see Section 3.2.1), following the results from simulations of the interaction between the binary system and the disk. Also, since the dust in the hole is optically thin, with very low density, we calculate a different temperature for each grain size, assuming radiative equilibrium between the grain absorption of stellar and shock radiation and emission of its own thermal radiation. Finally, we calculate the contribution to the SED of this optically thin material.

3.2.1. Dust Spatial Distribution

The gravitational interaction between the disk and the binary system sets constraints on the possible configurations of the disk material. A binary system is able to create a hole in the circumbinary disk (Lin & Papaloizou 1979; Artymowicz & Lubow 1994), because there is an inner region without stationary orbital configurations. The stable orbits for the particles are located very close to each star (circumstellar disks) and further out, in trajectories around both stars (circumbinary disk; Pichardo et al. 2005; Nagel & Pichardo 2008). As one would expect, the closer the material is to one of the stars, the less important its interaction with the other one. However, this configuration is perturbed when hydrodynamical interactions are taken into account. In this case, the material can lose angular momentum, and the resulting orbits end up connecting their point of origin in the circumbinary disk with each one of the circumstellar disks (Artymowicz & Lubow 1996; Günther & Kley 2002). In order to define a structure consistent with the results of these hydrodynamical simulations, we model the system as a circumbinary disk, two circumstellar disks and two streams of material, crossing the circumbinary disk hole, and connecting the outer circumbinary disk with the inner circumstellar disks. Note that the larger hole present in CS Cha (Espaillat et al. 2007a, 2011) cannot be explained only with the binary; thus, one requires another mechanism to produce it, perhaps another stellar mass companion or a multiple planetary system (Zhu et al. 2011; Dodson-Robinson & Salyk 2011). This will change the details of the structure in the hole, which we will not consider here.

Artymowicz & Lubow (1996) present smoothed particle hydrodynamics simulations (Monaghan 1992) which follow the material that is traveling from the circumbinary disk to the stars. They note that the accretion is modulated in time with a periodicity given by the orbital period of the binary. The simple picture is that the material in the edge of the disk is perturbed near the apocenter of the binary, but requires a time delay to fall to one of the stars, arriving around pericenter. Of course, the precise evolution depends on the mass ratio and eccentricity of the binary system (Artymowicz & Lubow 1996). A result of this is that the time of evolution of the material through the gap is of the order of the orbital period (P = 2482 days). This timescale is smaller than a typical dust coagulation timescale (∼10⁵ years for R ∼ 10 AU; Weidenschilling 1977; Birnstiel et al. 2010), so we assume that the grain composition and size distribution do not change during travel.

The circumbinary material in the outskirts of the disk initially evolves slowly, characterized by a viscous timescale evaluated at the radius of the circumbinary disk. Then, the material arrives to a radius defined by the gravitational interaction of the binary system with the disk (Artymowicz & Lubow 1994; Pichardo et al. 2005). Starting from this radius, the matter loses its axisymmetry and flows in two streams. The two streams are launched from two points called the saddle points. As noted in the last paragraph, the characteristic time for the evolution of the material in the gap is of the order of the binary orbital period. Thus, when the material arrives at the saddle points, it moves toward the inner circumstellar disks within a dynamical timescale.

The saddle points are analogous to the Lagrangian points in the three-body circular restricted problem in classical mechanics (Murray & Dermott 1999). For a circular binary, one can define a rotating coordinate system where the stars are at rest. Thus, there is a constant of motion, the Jacobi constant, and specifically two linear Lagrangian points (Murray & Dermott 1999), that define the locus where material with enough value for this constant is able to move from orbits around the binary system to orbits connecting the circumbinary disk with circumstellar disks. In order to extend the analogy to an eccentric system, we note that the gravitational potential can be expanded as an infinite sum of terms. In many cases, one can describe the potential with only a few of them. For an eccentric binary system, the potential can be simplified using the terms (m, l) = (0, 0) and (m, l) = (2, 1) of its expansion (Artymowicz & Lubow 1994). Moreover, this potential rotates at a rate Ω_b/2, where Ω_b is the mean angular velocity of the binary. For this simplified case, two points can be defined (the corresponding Lagrangian points), which corotate with the potential. These are the launching points for the streams mentioned in the last paragraph, i.e., the saddle points. Analogous to the circular case, the material that accretes from the external regions of the disk is able to cross the radius where the saddle points are located, primarily at these locations. The saddle points are located at a radius R_s = 1.587a, where a is the semimajor axis. Thus, we can conclude that at larger radii, the orbital backbone of the material orbits around the binary system (Pichardo et al. 2005). At smaller radii the material is located in two streams starting at the saddle points (Artymowicz & Lubow 1996; Günther & Kley 2002) and ending in the inner disks.

The configuration used in this work, which is consistent with the orbital restrictions given by the binary system presence, is a ring of material between R_s and R_wall, where R_wall is the position of the inner edge of the outer disk. At R < R_s, the material is located in the streams that connect the circumbinary and circumstellar disks. This is shown schematically in Figure 2. The ring corresponds to a region where the material evolves viscously, because this is consistent with the gravitational disk–stars interaction as described previously. The difference between this and the actual outer disk is that the ring is optically thin. As we mentioned above, in the case of CS Cha, for which we have Spitzer mid-IR fluxes, we do not include the contribution of the circumbinary disk to the SED, because it is too cold to emit in this wavelength range (Espaillat et al. 2007a, 2011). Also, we assume that the circumstellar disks, streams, and outer ring are all optically thin.

Zoom In Zoom Out Reset image size

3.2.2. Dust Temperature

We assume that the optically thin dust in the inner hole is in thermal equilibrium with the stellar radiation. However, the density in the gap is so low, that thermal equilibrium between gas and dust grains is not expected. Therefore, each grain with a different size would have a different temperature.

The critical gas density for thermal equilibrium between dust grains and gas can be estimated by assuming that the collisional timescale is equal to the thermal timescale (e.g., Chiang & Goldreich 1997; Glassgold et al. 2004). The collisional timescale is given by

$$\begin{equation} \tau _{{\rm col}}={1\over n_{{\rm H}}\bar{v}\pi a^{2}}, \end{equation} \tag{ 1 }$$

where n_H is the number density of the molecular hydrogen, $\bar{v}$ is the mean velocity of a molecule, and a is the radius of the grain. We consider that the velocity is thermal, however, it is possible that there is a turbulent component, which amounts to a fraction of this velocity. Using the fact that the typical velocity associated with turbulence in the α-parameterization of the viscosity (Shakura & Sunyaev 1973) in protoplanetary disks (Hawley et al. 1995) is around 0.01 of the sound speed, we can conclude that including the turbulent velocity does not change our following conclusions. The thermal timescale can be estimated as the time that a particle of dust takes to radiate all its thermal energy at temperature T_d(a),

$$\begin{equation} \tau _{{\rm ther}}={kT_{d}(a)\over 4\pi a^{2}\sigma _{{\rm SB}}T_{d}(a)^{4}}, \end{equation} \tag{ 2 }$$

where k is the Boltzmann constant and σ_SB is the Stefan–Boltzmann constant.

Thus, the critical density is (Glassgold et al. 2004)

$$\begin{equation} n_{{\rm cr}}={4\sigma _{{\rm SB}}T_{d}(a)^{3}\over k\bar{v}}, \end{equation} \tag{ 3 }$$

and substituting typical numerical values, it can be written as

$$\begin{equation} n_{{\rm cr}}=(1.133\times 10^{8}\ {\rm cm}^{-3}){T_{d}^{3}\over T_{g}^{1/2}}, \end{equation} \tag{ 4 }$$

where T_g is the gas temperature.

Next, we would estimate the gas density associated with the optically thin material to compare it with this critical density. If we assume T_g ≈ T_d, and take T_d ≈ 100 K as characteristic of the streams, then n_cr = 3.34 × 10⁻¹¹ g cm⁻³. Note that a change in T_g of an order of magnitude corresponds to a factor of three change in n_cr (see Equation (4)), which does not modify our conclusions.

In order to get an order of magnitude for the density in the streams, we note that the area covered by the streams is around 0.1πR²_s (Artymowicz & Lubow 1996). We assume that the material falls from the inner edge of the ring toward the star in a few times the free-fall time, which is 1.5 years. This is consistent with the fact that in Artymowicz & Lubow (1996), the streams change on the order of an orbital period of the binary, T_orb ∼ 7 years. Thus, according to a mass accretion rate of $\dot{M}=10^{-8}\,M_{\odot }\ {\rm yr}^{-1}$ (Espaillat et al. 2007a), the mass in the streams should be around 7 × 10⁻⁸ M_☉. Note that in the hydrodynamical simulations of a disk around a binary system, Günther & Kley (2002) assume a typical height in terms of the radius R in the disk, without dependence on the azimuthal angle, which is given by

$$\begin{equation} H(R)=\left(\Sigma _{i=1,2}{GM_{i}\over c_{s}^{2}|R-R_{i}|^{3}}\right)^{-1/2}, \end{equation} \tag{ 5 }$$

where G is the gravitational constant, M_i is the mass of each of the stars, c_s is the sound speed, and finally R_i is the distance of each star with respect to the origin. For Case 1 and using a typical temperature of 100 K for molecular hydrogen and |R − R_i| = R_s, H(R) = 0.049 R_s. Thus, a typical value for H is 0.1 R_s, which we consider here. Differences in H do not change the following conclusions.

Given the mass in the streams, the area covered and the height calculated previously, the characteristic density in the streams is around n_str = 10⁻¹⁴ g cm⁻³, which means that n_str ≪ n_cr. This implies that the gas and the dust grains are not in thermal equilibrium, and since the absorption coefficient of the grains depends on grain size, a, grains with different sizes have different temperatures, T(a).

We divide the grain sizes into 100 bins equally spaced in log (a) between a_min = 0.005 μm and a_max = 4 μm, distributed with a power-law profile with exponent p = −3.5. We take this value of a_max according to the detailed modeling of Espaillat et al. (2011). For each size, we calculate the temperature from radiative equilibrium between the energy of the stellar and shock radiation absorbed by the grain and the energy lost by radiation. As expected, the temperature of the smallest grains is higher than the temperature of the larger grains (see discussion in Section 5.2 and Figure 3). We also estimate a mean temperature, as a unique temperature for the whole grain distribution, assuming the grains have a Mathis–Rumpl–Nordsieck (MRN) power-law size distribution, and use this temperature as a reference. This temperature is also calculated assuming radiative equilibrium with the radiation fields of the stars and their accretion shocks, but using mean absorption coefficients. As expected, at the same position in space, the temperature of the smallest grains is higher and the temperature of the larger grains is lower than this mean temperature, due to a higher heating efficiency. Figure 3 illustrates this by showing temperature versus distance in the primary's stream (Case 1, GO observations) for the smallest grain (0.005 μm), the largest grain (4 μm), and for a power-law distribution of grain sizes between 0.005 μm and 4 μm.

Zoom In Zoom Out Reset image size

Remember that the estimation of the temperature is done for all the locations where there is material, and in each case the temperature of each grain size is calculated. The grain composition adopted here is amorphous olivine, consistent with Espaillat et al. (2011), where 95% is olivine and 5% is crystalline silicates. In addition, we include organics and troilite. The abundances are ζ_sil = 0.0034, ζ_org = 0.001, and ζ_troi = 0.000768, for silicates, organics, and troilite, respectively.

For the best model for CS Cha in Espaillat et al. (2007a) and Espaillat et al. (2011), the inner edge of the circumbinary disk, which we refer to as the "wall," is located at R_wall = 43 AU and R_wall = 38 AU, respectively. The difference found in R_wall for CS Cha was due to using different dust opacities. The models of Espaillat et al. (2011) use a different distribution of silicates with respect to the models for CS Cha in Espaillat et al. (2007a). In the first modeling, 60% of the dust is silicates compared with 88% in the other case. For both cases, crystalline silicates correspond to a small amount. Here we use the most recent value, R_wall = 38 AU, because it corresponds to a better fit of the 10 μm band. At such a large distance, there is almost no difference between a model with one star and a model with two stars, because a ≪ R_wall (a = 3.703 AU for Case 1, a = 3.587 AU for Case 2). Thus, it is safe to assume the wall parameters to be the same inferred by Espaillat et al. (2011), including the wall height (h = 7 AU) and the inclination angle between the disk axis and the line of sight (i = 60°).

As noted in Section 2, we fix the value estimated by Espaillat et al. (2007a) for the mass accretion rate at 1.2 × 10⁻⁸ M_☉ yr⁻¹. We distributed this value, 0.4 × 10⁻⁸ M_☉ yr⁻¹ and 0.8 × 10⁻⁸ M_☉ yr⁻¹, for the primary and secondary, respectively. This is in agreement with Artymowicz & Lubow (1996), where the simulations with e = 0.1 show that the mass accretion rate to the secondary is a factor of two larger than the mass accretion rate to the primary star. Notice that the accretion luminosities are an order of magnitude lower than the star luminosity. As a result, we expect that the precise value of the former does not make a substantial difference in the resulting SED. The stellar parameters for Cases 1 and 2 are summarized in Table 1.

The streams start at the saddle points of the simplified potential (Artymowicz & Lubow 1996; see Section 3.2.1) and end at each circumstellar disk. For simplicity, their shape is a line that connects both positions. Although the thickness of the streams is a free parameter, one condition should be held: the flows fill around 10% of the surface of the hole (Artymowicz & Lubow 1996). Following the simulations of Artymowicz & Lubow (1996), we can safely assume that the density of streams and the outer ring are one tenth of the values associated with the compact circumstellar disks. We assume that both disks have a minimum radius of R_min = 0.1 AU and a maximum radius of R_max = 0.5 AU. As we mentioned above, the mass accretion rate to the secondary is about two times the value associated with the primary following Artymowicz & Lubow (1996, see also Bate 2000). For the modeling, we translate this as a condition on the filling factor inside the hole associated with each stream. Thus, we assume that the area occupied by the primary star stream is half the area of the stream reaching the secondary. The last requirement consistent with Artymowicz & Lubow (1996) is that the mass from the circumbinary disk that falls into any of the streams at apocenter is around twice the value at pericenter. Finally, in order to characterize the mass of the streams for every orbital location, we modulated the area of the streams (equivalently the mass) with a sinusoidal function in the azimuthal angle. The advantage of using this information to define the stream structure is that their geometrical parameters are fixed. The density on each stream is constant, and it is fixed using the fit for the observed spectrum (see Section 5).

The eccentricity is a parameter that is not restricted by the observations of Guenther et al. (2007). Artymowicz & Lubow (1996) point out that the variability of the stream mass depends on the eccentricity. The amount of material associated with the streams depends on the distance of the stars to the wall at apocenter, and this strongly depends on the eccentricity. Note that for a circular binary, the minimum distance of each star to the disk does not change for the axisymmetric wall; thus, the structure is stationary and there is no variability in the amount of material in the hole. Thus, in order to explain the variability, the system must be eccentric. Also, the observed SED variations suggest significant changes in the amount of material in the streams (17% according to Espaillat et al. 2011). For the modeling, in order to look for a best fit we take two values for the eccentricity, e = 0.1 and e = 0.2.

The remaining parameter needed to fully characterize the binary system is ϕ. It is the angular position of the stars along the orbit at the time we are observing the system: ϕ = 0.5π corresponds to pericenter and ϕ = 1.5π to apocenter. The values we try in order to find the best fit for the earlier observation are ϕ₁ = (0.0, 0.5, 1.0, 1.5)π. The set of models tried is shown in Table 2. The second epoch corresponds to ϕ₂, constrained by the time elapsed between the two observations, i.e., 1065 days. We solve the equations of motion of the binary system, looking for configurations separated by this time.

Finally, the dust mass surface density for the region assigned to the streams is a free parameter. The density is increased from zero until the flux level in the SED is reached.

The fit is done independently for both epochs. In order to find the best-fitting model, we search for the model that minimizes χ². Because this is a small set of the parameter space, we cannot claim that we find a unique model. Due to the large number of assumptions and unknowns (for example, the exact configuration of the streams), degeneracy of the SED fitting is expected. Thus, one can only claim that we find models consistent with the observations.

The fit for the GTO observation in 2005 is shown in Figure 4 for Cases 1 and 2. The fit for the GO observation in 2008 is shown in Figure 5, also for both cases. The parameters for the best-fit model for Cases 1 and 2 are the same, this is e = 0.1 and ϕ = (0.5, 1.38)π. The shape of the streams' structure for the fit is shown in Figure 6 for Case 1, compared to the GTO and GO observation of CS Cha. The streams for Case 2 are almost identical to the ones for Case 1, thus, we do not show them here.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The fit in the last section is done for both epochs, but consistently taking into account the time difference between the observations. Thus, the variability can be explained with the previous fits. In Figure 7, we show the Case 1 SED fits for both observations, which were previously presented in Figures 4 and 5. Figure 8 is the same but for Case 2, showing the fits also presented in Figures 4 and 5. We note that the best fit for either model with e = 0.1 or e = 0.2 is ϕ₁ = 0.5. In order to check this tendency for every e, we ran models for e = 0.5 and e = 0.9. The tendency is corroborated, but the best models according to χ² are worst for larger e. We ran models for e = 0.1 for other values of ϕ₁, such that the whole set is ϕ₁ = (0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75). Even for this larger set, the best fit is found for ϕ = 0.5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Given the observational constraints and the properties we have fixed for the inner wall based on our previous modeling (Espaillat et al. 2011), we find that any eccentricity that satisfies e > 0 and ϕ₁ around 0.5 produces a configuration consistent with the observed SEDs in the two epochs. This means that the eccentricity e is not constrained by our models and it appears that only observations of the detailed orbit of the binary system would allow us to quantify this parameter.

The SEDs of the present models (for both, Case 1 and Case 2) are similar to the SED of a single star model presented by Espaillat et al. (2011). In the single star model the SED variation is related to a change in the mass of the optically thin dust in the inner hole. In the present model with two central stars, we also consider a change in dust mass, but besides this there is a source of variability given by the change in illumination produced by the stars, which is absent in a single star model.

As we show, the variability can be explained by the differences for each orbital configuration of the binary system, in the illumination of the material in the circumstellar disks and streams. Also the amount and the geometrical configuration of the dust is important in terms of the orbital configuration. A visualization of this effect is given in Figure 9, where for the Case 1 fit, the intensity for a wavelength of 10 μm is presented. One can clearly see that the contribution to the emission is larger for the configuration associated with the GO observation (Figure 6). Notice that in Figure 9, the variation in area for the streams means a larger amount of mass, and the gray scale corresponds to changes in the illumination.

Zoom In Zoom Out Reset image size

The associated dust mass in circumstellar disks and streams is around 10⁻¹² M_☉, of the same order of the estimate by Espaillat et al. (2011), where the material was located only in one circumstellar disk. The dust surface density on each stream is uniform and can be calculated with the total mass of dust and the area associated with each stream. In Table 3, we present the dust surface density for the primary and the secondary streams, for Cases 1 and 2, and for the configurations associated with both epochs. In order to estimate from the model a mass accretion rate we do the following. We assume that the gas attached to the dust in the disks and streams is responsible for the mass accretion rate estimated in Espaillat et al. (2007a). Using the dust abundances, the mass of gas is of the order of 10⁻¹⁰ M_☉, if we now use a free-fall time from the saddle point as the fastest time for the accretion, the mass accretion rate is around 10⁻¹⁰ M_☉ yr⁻¹. This value is two orders of magnitude smaller than the estimate in Espaillat et al. (2007a). A possible explanation for this is that the dust-to-gas ratio is lower than usually expected. A way to do this is through a filtering effect in the inner edge of the circumbinary disk. Rice et al. (2006) note that the pressure gradient in this location acts as a filter, only allowing small dust grains to trespass it, holding the larger ones. However, the typical critical size is around 10 μm, which is larger than the maximum value taken in our work, a_max = 4 μm. Another idea is that the mass accretion rate inferred from the U-band excess has uncertainties due to the estimation process in itself and due to the fact that the system is a binary. An improvement of this estimate requires a detailed description of the evolution of the particles (solid and gaseous), since they leave the inner edge of the circumbinary disk to finally accrete at one of the stars (Günther & Kley 2002). Note that a U-band excess is not always indicative of accretion, because one has to check the shape of the lines profiles to safely assume that the system is accreting (Ingleby et al. 2011). However, a complete study able to elucidate the dynamics of the hole material is beyond the scope of this paper.

In previous works (D'Alessio et al. 2005; Nagel et al. 2010; Espaillat et al. 2011), the dust temperature in the hole and the emissivity are grain-size-mean values calculated using an opacity, given by the average opacity over the grain size distribution. In here, since the dust in the hole has such a low density and cannot be considered in thermal equilibrium with the gas (see Section 3.2.2; Chiang & Goldreich 1997; Glassgold et al. 2004), we calculate a dust temperature and emissivity as a function of grain radii. For the sake of evaluating the differences between both approaches, we calculate a model with the grain-size-mean opacity to fit the SED. We find that the resulting SEDs are very similar, but for the model with the grain-size-mean opacity assumption, the mass of the optically thin material is 19% of the mass of the models where no thermal equilibrium is assumed. This difference between both approaches means that it is worthwhile to do the proper estimate.