A&A 422, 171-176 (2004)
DOI: 10.1051/0004-6361:20040287

The vertical structure of T Tauri accretion discs

IV. Self-irradiation of the disc in the FU Orionis outburst phase

R. Lachaume

Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69, 53121 Bonn, Germany

Received 17 February 2004 / Accepted 8 April 2004

Abstract
I investigate the self-irradiation of intensively accreting circumstellar discs (backwarmed discs). This is modelled using the two-layer disc approach by Lachaume et al. (2003) that includes heating by viscous dissipation and by an external source of radiation. The disc is made of a surface layer directly heated by the viscous luminosity of the central parts of the disc, and of an interior heated by viscosity as well as by reprocessed radiation from the surface. This model convincingly accounts for the infrared excess of some FU Orionis objects in the range 1-200 $\mu$m and supports the backwarmed disc hypothesis sometimes invoked to explain the mid- and far-infrared excesses whose origins are still under debate. Detailed simulation of the vertical radiative transfer in the presence of backwarming is still needed to corroborate these results and spectroscopically constrain the properties of intensively accreting discs.

Key words: accretion, accretion disks - stars: individual: FU Ori - stars: individual: Z CMa - stars: individual: V1057 Cyg - stars: individual: V1515 Cyg - stars: circumstellar matter

1 Introduction

It is now widely accepted that circumstellar discs accompany the process of star formation, all the more as such discs have already been imaged in the millimetre (ex. DG Tau, Dutrey et al. 1994) and in the infrared (ex. HK Tau, Stapelfeldt et al. 1998). Matter accreting from the disc onto the star is supposed to build up solar-mass stars in the timescale of a Myr; yet recent studies of star forming regions indicate way too low accretion rates among T Tauri stars (TTS) - in the range $10^{-10}{-}10^{-7}~{M_\odot}/{\rm yr}$, as shown by Gullbring et al. (1998) in the Taurus region and Robberto et al. (2003) in the Trapezium. Part of the matter is supposed to be accreted in earlier phases of high accretion but young stellar objects (YSOs) are then so embedded that their disc cannot be optically observed (see André 1994, for the stages of evolution). The other part might also be accreted during brief, periodic phases of intense accretion of otherwise quiescent TTS (Hartmann & Kenyon 1996): FU Orionis objects (FUors), which feature accretion rates up to $10^{-4}~{M_\odot}/{\rm yr}$ and which are seen to undergo an increase in luminosity of more than 4 visual magnitudes in a few years (Herbig 1977,1966), may represent such a phase.

FUors are convincingly modelled by a self-heated viscous accretion disc (active disc) that overwhelms the stellar light and the properties of FUor outbursts are also well-studied (Bell et al. 1995; Hartmann & Kenyon 1985). Yet active disc models fail to explain the mid- and far-infrared excess (e.g. Simon & Joyce 1988) as well as weak silicate feature at $10~\mu$m sometimes present in emission (Hanner et al. 1998); these properties seem more typical of an irradiated disc. Alternative models have tackled this problem: Lodato & Bertin (2001) showed that the self-gravity of the disc can trigger an instability that produces additional warming and were able to reproduce the SED. However, this model cannot account for the silicate feature, which requires a temperature inversion at the surface of the disc. A more pleasing approach is the presence of a circumstellar envelope proposed by Adams et al. (1987): on one (observational) hand, Kenyon & Hartmann (1991) successfully fit the spectra of two FUors with such an envelope. On the other (theoretical) hand, the envelope could serve as a reservoir of infalling material, replenishing the disc between outbursts. On the third hand, the envelope can account for a part of the extinction observed in these objects.

Kenyon & Hartmann (1991) and Bell (1999) also proposed that the inner hot parts of the disc are bright enough to heat up the outer ones (backwarming) and produce an irradiated disc-like SED. Their simulation uses a black-body disc model, i.e. without vertical temperature profile, and convincingly explains the order of the excess at 30-100 $\mu$m, yet it fails to reproduce the SED at 10-30 $\mu$m and does not predict a silicate feature in emission. Malbet & Bertout (1991), Calvet et al. (1991) and Chiang & Goldreich (1997) showed that irradiation can produce a hot disc surface radiating at shorter wavelengths and account for emission features (CO bands, silicates). Using their results, Lachaume et al. (2003) developed a two-layer disc model in which irradiation by a central star and viscosity are taken into account. In this paper, I shall use this model with the disc itself as a source of radiation and investigate the backwarming in FUors.

2 The model

I use the two-layer approach described in detail by Lachaume et al. (2003): the surface of the disc is directly irradiated by the incoming radiation with an optical thickness of unity in the visible along the slanted path of the incident beam. Its vertical optical thickness in the infrared, to which radiation is mostly reprocessed, is then much smaller than one. Therefore, it is at a higher temperature than its effective temperature. The interior of the disc is both backwarmed by the surface and heated by viscous dissipation, modelled with the $\alpha$ viscosity prescription by Shakura & Sunyaev (1973). A scheme of the model is presented in Fig. 1.

2.1 Irradiation

Since most of the flux originates from the first few solar radii of the disc, the irradiation source is seen as a central, infinitely flat ring from the the outer regions of the disc (1-100 AU). The radial flux distribution is, according to the standard disc model by Shakura & Sunyaev (1973),

\begin{displaymath}% F(r) = \frac{3}{8} \left(1-\sqrt{\frac{r_0}{r}}\right) \frac{G{M_\star}{\dot{M}}}{\pi r^3}, \end{displaymath} (1)

where ${M_\star}$ is the mass of the star, $r_0 < {R_{{\rm min}}}$ a constant, and  ${R_{{\rm min}}}$ the inner radius of the disc. The value of r0 depends on the torque exerted on the inner rim of the disc. If the latter extends down to the star, the standard model assumes that the former is linked to the deceleration of the matter falling from the supersonic disc to the subsonic equator of the star, and $r_0 \approx {R_\star}$, the stellar radius. However, the torque exerted on a disc with an inner gap remains unknown; it should depend on the magnetic coupling between the star and the disc. Without clear theoretical or observational hints, I decided to assume $r_0 = {R_\star}$, as in the standard model.


  \begin{figure} \par\includegraphics[width=6.5cm,clip]{0287f1.eps} \end{figure} Figure 1: Geometry of a backwarmed disc. The incidence angle $\phi $ of the incoming radiation onto the surface and the angle $\theta $ at which the inner regions are seen from it depend on the profile of the surface H(r).
Open with DEXTER

The total luminosity of the disc amounts to

                    $\displaystyle % {L_{{\rm disc}}}$ = $\displaystyle \int_{{R_{{\rm min}}}}^{+\infty} 2\pi r F(r) {\rm d} r$  
  = $\displaystyle \left(3-2 u \right) u^2 ~ \frac{G{M_\star}{\dot{M}}}{4{R_\star}},$ (2)

where $u = \sqrt{{R_\star}/{R_{{\rm min}}}} < 1$. At a given location on the disc's surface, the heating flux, ignoring the shadowing of some parts of the disc by the star, is:

 \begin{displaymath}% {F_{{\rm irr}}}= (2/3+\cos\theta) \sin\theta \sin\phi ~ \frac{L_{{\rm disc}}}{r^2 + H^2}, \end{displaymath} (3)

where $\theta $ is the angle at which the heating ring is seen, $\phi $the incidence angle onto the surface of the disc, and H the altitude of the surface above the mid-plane. In Eq. (3) the first term models the "limb-darkening'', while the second and third ones are projection effects. With an infinitely flat ring
    $\displaystyle \theta = {\rm Arctan} \frac Hr,$ (4)
    $\displaystyle \phi = {\rm Arctan} \left( {\frac{{\rm d}^{}H}{{\rm d}r^{}}} \right) - \theta.$ (5)

The characteristic temperature of the irradiation is needed to determine the optical thickness of the surface layer - it depends on the penetration of the incoming radiation. We use the flux-weighted average of the viscous effective temperature  ${T_{\rm V}}(r)$:
                             $\displaystyle % {T_{{\rm irr}}}^4$ = $\displaystyle \frac 1{L_{{\rm disc}}}\int_{{R_{{\rm min}}}}^{+\infty} 2\pi F(r) {T_{\rm V}}(r)^4 {\rm d} r$  
  = $\displaystyle \frac{45 - 80 u + 36 u^2}{3 - 2 u} u^6 ~ \frac{G{M_\star}{\dot{M}}}{160\sigma\pi{R_\star}^3}\cdot$ (6)

2.2 Location of the surface

Chiang et al. (2001) provide a method to consistently find the location of the irradiated surface in an irradiated two-layer disc. However, the method is more difficult to implement in models that also include viscous heating, in which viscosity-dominated regions may present self-shadowing. Instead, we used a simpler approach proposed by Chiang & Goldreich (1997) and state that irradiation-dominated regions meet

 \begin{displaymath}% H \propto h, \end{displaymath} (7)

where h is the vertical pressure scale-height of the disc. It appears that such a prescription results in a surface profile H(r) close to that predicted by the self-consistent method in the irradiation-dominated regions (Lachaume et al. 2003). Although one expects it to fail in viscosity-dominated regions, the impact on the disc structure is low, since these regions are by definition not affected by the irradiation.

We now look for a self-similar solution $H/r = r^\gamma$. Within the approximation $H(r)/r \ll 1$, Eq. (3) reads

 
                            $\displaystyle % {F_{{\rm irr}}}$ $\textstyle \approx$ $\displaystyle 2/3 \left(\frac{H}{r}\right) \left((\gamma-1) \frac{H}{r}\right) \frac1{r^2}$  
  $\textstyle \propto$ h2 r-4. (8)

The pressure scale height is obtained from the vertical gravity field of the star:

\begin{displaymath}% h \approx r^{3/2} \sqrt{\frac{{k_{\rm B}}{T_{\rm c}}}{G{M_\star}{m_{\rm H}}}} \end{displaymath} (9)

where

\begin{displaymath}% {T_{\rm c}}= ({F_{{\rm irr}}}/2)^{1/4} \end{displaymath} (10)

is the temperature of the interior of the disc in irradiation-dominated regions. The last two equations, together with $H \propto h$, yield

 \begin{displaymath}% H \propto r^{3/2} {F_{{\rm irr}}}^{1/8}. \end{displaymath} (11)

By subsituting Eqs. (8) into (11) one finally derives

\begin{displaymath}% H/r \propto r^{1/3}. \end{displaymath} (12)

This determination does not provide a estimate for the proportionality constant, which I leave as a free parameter in the model.


  \begin{figure} \par\mbox{\includegraphics[width=7.3cm,clip]{0287f2a.eps} \hspace... ...\hspace*{1.3cm} \includegraphics[width=7.3cm,clip]{0287f2d.eps} } \end{figure} Figure 2: The SED of FUors. Solid line: two-layer disc model including viscous heating and backwarming. Dashed line: black-body disc model including viscous heating and backwarming. Dotted line: disc model with viscous heating alone.
Open with DEXTER

3 Results and discussion

3.1 Strategy

I apply this model to unembedded FUors, since more embedded objects cannot be properly described without considering reprocessing of the disc emission by the surrounding dusty material, which is beyond the scope of this paper. As a criterion for "unembeddedness'' I choose a decreasing SED in the infrared, which leads me to consider three confirmed FUors, FU Ori, V1057 Cyg, and V1515 Cyg, and one FUor candidate, Z CMa (Kenyon & Hartmann 1991). Their unembeddedness is confirmed by their limited extinction $A_{\rm V} \la 4$ mag (Bell et al. 1995; Cohen & Kuhi 1979). Photometry data are taken from Kenyon & Hartmann (1991) for V1057 Cyg and V1515 Cyg, from Thiébaut et al. (1995) for Z CMa, and from the Gezari et al. (1999) catalogue for FU Ori.

The models have been computed using wavelength-dependent grain opacities for dust aggregates by Henning & Stognienko (1996), needed for the determination of the optical thickness of the surface layer and for the computation of the SED. At higher temperatures, I used grey gas opacities by Bell & Lin (1994). Model parameters  ${R_{{\rm min}}}$, ${M_\star}{\dot{M}}$, and  $A_{\rm V}$ can be determined by first fitting an active disc model to the SED in the range 0.3-3 $\mu$m; they do not depend on the strength of the mid-IR excess to be modelled. The constant $H/r\ (10~\rm AU)$ is then adjusted to fit the backwarmed disc model to the mid- and far-IR excess. ${M_\star}{\dot{M}}$ and $H/r\ (10~\rm AU)$ are unambiguouisly constrained (by the 1-3 $\mu$m flux and the 10-50 $\mu$m SED respectively), while  ${R_{{\rm min}}}$ and  $A_{\rm V}$ are derived from the SED at $\lambda \la 0.5~\mu$m.

In all models, the viscosity parameter is $\alpha = 0.1$ and the albedo of the surface ${\mathcal{A}}= 0.2$. $\alpha$ has almost no influence on the SED at $\lambda \la 100~\mu$m because these discs are optically thick, and the presence of the albedo increases the visible and near-IR SED by a few percent while decreasing the mid- and far-IR excess by up to 20%. Z CMa and FU Ori, which exhibit observational hints of not beeing seen pole-on, are assumed to be seen with $i = 40{^{\rm o}}$, other targets are modelled as pole-on ( $i = 0{^{\rm o}}$). The inclination has little influence for angles $i \la 45{^{\rm o}}$ and mostly scales down the visible and near-IR SED by a factor $\cos i$.

Table 1: Fundamental parameters of FUor models: accretion rate, inner radius of the disc, flaring H/r, type of dust opacity  $\kappa _\lambda $, and extinction $A_{\rm V}$.

3.2 Results

The fits are displayed in Fig. 2 and their fundamental parameters in Table 1. For comparison a standard viscous disc model and a black-body backwarmed disc, that is with no superheated surface, as in Kenyon & Hartmann (1991), are shown along with the two-layer model fit. It appears that the black-body model gives a slightly higher far-infrared excess and a smaller mid-infrared excess than the present model and does not produce the 10 $\mu$m silicate feature.

3.2.1 FU Ori

Up to $10~\mu$m the SED of FU Ori can be modelled by a standard viscous disc (Fig. 2, dotted line of the top-left panel). The accretion rate ${\dot{M}}{M_\star}= {3.2\times10^{-5}}~{M_\odot}^2/{\rm yr}$ corresponds to typical values in the literature ($\approx$ $3/\cos i~{M_\odot}^2/{\rm yr}$ where i is the inclination of the disc according to Hartmann & Kenyon 1996; Malbet et al. 1998). Nevertheless, the 5-12 $\mu$m SED is about 30% weaker as expected from the fit of 0.5-3 $\mu$m fluxes.

A black-body backwarmed disc (see Kenyon & Hartmann 1991) can explain the far-IR excess but it either overestimates the 30-100 $\mu$m fluxes (same figure, dashed line) or underestimates the 20-30 and 100-200 $\mu$m SED (not shown). If the backwarmed disc model comprises a hot surface resulting from direct irradiation, then the SED can be convingly fitted in the range 20-200 $\mu$m (solid line) with $H/r \approx 0.2$ at 10 AU from the central star. This fit, however, overestimates the M and N fluxes by a factor 2.5. Interestingly, the model predicts a $10~\mu$m silicate feature with the same intensity as in the observations, but with an overestimated continuum.

That this centro-symmetric disc model does not accurately reproduce the observations should not be taken too seriously as an argument against backwarming, since the target appears to be a rather complicated object: Malbet et al. (1998) found by means of optical long-baseline interferometry that the disc should present a hot spot located 15 AU from the star with a flux ratio of 4-5 mag in K; Wang et al. (2004) also detected a companion at 200 AU (0.5'') with a similar flux ratio.

3.2.2 V1057 Cyg

V1057 Cyg is well described with an accretion disc featuring an accretion rate ${M_\star}{\dot{M}}= {2.0\times10^{-5}}~{M_\odot}^2/{\rm yr}$ (in the range predicted by Kenyon et al. 1988), as shown in Fig. 2 (top-right panel, dotted line). Backwarming is needed to account for the SED in the range 10-100 $\mu$m with a flaring $H/r \approx 0.35$ at 10 AU from the star. The black-body approach either overestimates the 20-50 $\mu$m fluxes by a factor 3 (dashed line) or underestimates the 50-100 $\mu$m SED by a similar factor (not shown). The two-layer disc model fits the data well in the range 1-100 $\mu$m (solid line). It is worth noticing that the model does produce a strong silicate feature while the flux in the N-band seems to hint to the contrary - it is at the level of the continuum predicted by the model.

It is however worth noticing that this disc, though comparable to FU Ori in terms of accretion properties, presents a much larger flaring. If, for some unknown reason, it features a lower viscosity, this could account for a larger disc mass and a more flared structure. On the other hand, it cannot be ruled out that part of the far-IR excess is not to be imputed to backwarming and could result from other phenomena, like an envelope or self-gravity-related energy dissipation, in which case the flaring requirement is only an upper limit.

3.2.3 V1515 Cyg

V1515 Cyg presents an SED very similar to that of V1057 Cyg at $\lambda \la 20~\mu$m but shifted down by a factor 5, which the larger distance roughly explains. It is therefore not surprising to find a similar accretion rate ${M_\star}{\dot{M}}= {1.3\times10^{-5}}~{M_\odot}^2/{\rm yr}$ with a standard disc model (Fig. 2, dotted line in the bottom-left panel). Despite of the lack of data at $\lambda \ga 20~\mu$m, the 10-20 $\mu$m fluxes require backwarming with a flaring $H/r \approx 0.35$ at 10 AU from the star. However, it is not possible to convincingly discriminate between the black-body and the two-layer disc models.

3.2.4 Z CMa

Z CMa is a close (0.1'') binary (Koresko et al. 1989) YSO with FUor characteristics. The total flux can be described as emerging from an intensively accreting FUor ( ${\dot{M}}\sim 10^{-3}~{M_\odot}/{\rm yr}$) for $\lambda \la 10~\mu$m, yet it fails at reproducing the almost flat far-IR SED. A closer look at the SED of the individual components (see Thiébaut et al. 1995, and references therein) hint that the primary is a cool YSO in its primary phase of evolution (it peaks at 5 $\mu$m) while the secondary features an FU Ori-like SED. Z CMa B is fairly well described with an accretion disc ${M_\star}{\dot{M}}= {1.4\times10^{-5}}~{M_\odot}/{\rm yr}$, but the lack of data at $\lambda \ga 5~\mu$m rules out checking the backwarming hypothesis (Fig. 2, bottom-right panel).

3.3 Reliability of the model

This two-layer model proves that backwarming is a credible explanation for the mid- and far-IR excess among a few FUors, without the need to invoke the contribution of a circumstellar envelope. It does not self-consistently tackle the determination of the vertical structure, so I discuss here the possible caveats linked to these simplifications.

3.3.1 Simplification of the vertical structure

Freeing the flaring parameter H/h may seem hazardous, since its value strongly affects the amount of irradiation caught by the disc and therefore the strength of the IR excess (ex. Kenyon & Hartmann 1991, for backwarmed discs), so a self-consistent determination of the irradiated surface by more elaborate simulations (Dullemond et al. 2002; D'Alessio et al. 1998) may seem more able to corroborate or invalidate the backwarming scenario. Yet the latter simulations hide a flaring hypothesis, since they depend on the up-to-now poorly constrained mass column in the disc; in other words, they do not help to rule out or support the backwarming scenario better than an empirical approach, as far as the diagnosis is based on the SED. In the future, they should rely on better constrained masses and provide a reliable diagnostic, when the Atacama Large Millimetre Array (ALMA) probes the disc mass at the 10 AU-scale or when IR long-baseline interferometry measures flaring-related asymmetries at the AU-scale - for instance with AMBER on the Very Large Telescope Interferometer (VLTI).

The two-layer approximation may also seem problematic, but Dullemond & Natta (2003) showed that its SED predictions stand close to those of full transfer simulations. On the contrary, it is not fit for spectroscopic studies (Aikawa et al. 2002) since spectral lines are very sensitive to the temperature and density profiles; in particular, it cannot be used to issue a diagnostic using the strength and shape of the silicate feature.

3.3.2 Ignoring the disc's own gravity

FUor disc models predict large masses, as shown in Table 2, so that the gravity of the disc can overwhelm the stellar gravity in the outer, massive parts of the disc (self-gravitating discs, see Bertin & Lodato 1999; Huré 2000, and references therein). These works single out three effects:

1.
 an additional vertical gravitational field that flattens the outer parts of the disc and should prevent them from being illuminated by the central source of radiation (Huré 2000), which would include backwarming from inner disc regions;
2.
  an additional horizontal gravitational field so that the disc is no longer in Keplerian rotation. This affects the viscous dissipation and the distribution of the matter. (Bertin & Lodato 1999);
3.
  Possible instabilities in self-gravitating regions. The disc is then more likely to depart from a symmetrical geometry and presents an additional dissipation of energy, hence an enhanced IR excess. Lodato & Bertin (2001) proposed this phenomenon to explain the far-IR excess among FUors.
Phenomenon 1 is critical when determining the SED of a backwarmed disc. For each model fit, I computed an estimate of the radius ${R_{\rm g}}$ at which the disc undergoes a transition toward self-gravity, namely where the vertical gravity field produced by the disc at its surface equals that of the central star. As shown in Table 2, the three successfully fit FUors should have ${R_{\rm g}}\approx 50~\rm AU$ with the hypothesis of a high viscosity coefficient $\alpha = 0.1$. In this case, self-gravity should mostly affect the coolest regions and therefore the SED at $\lambda \ga 50~\mu$m. Phenomena 2, 3 mostly affect viscous dissipation and can account for the far-IR SED for $\alpha \sim 0.01$ (Lodato & Bertin 2001). It is not impossible that the nature of the IR excess should be imputed to contributions of both backwarming and self-gravity.

It should be emphasised that no study of self-gravity included the influence of irradiation. The latter possesses a conterweight effect: by heating the outer parts of the disc, it increases the flaring and - at least within the $\alpha$ viscosity prescription - decreases the amount of matter needed to sustain a given accretion rate. Furthermore, it is not clear whether FU Ori discs are as massive as expected from standard steady-disc models: it might well happen during the outbursts that only the innerer parts undergo a sudden increase of the accretion rate, while the outer parts of the disc conserve TTS-like properties (Bell et al. 1995).

Table 2: Disc mass in FUor models: column density $\Sigma $ at 1 and 10  AU from the star, cumulated disc mass  ${M_{\rm d}}$ within 10 and 100 AU, and radius ${R_{\rm g}}$ at which the disc's vertical gravity field equals that of the star. Masses are lower values and radii ${R_{\rm g}}$ upper values, obtained with a high viscosity parameter $\alpha = 0.1$.

4 Conclusion

Using the viscous and irradiated disc model by Lachaume et al. (2003), it is possible to account for the strength of the mid- and far-IR excess among a few FUors using the viscous luminosity of the inner parts of the disc as an irradiation source. Yet the model makes several approximations in the determination of the structure: the temperature profile is assimilated to two isothermal layers, a surface and an interior, which should not prove critical in the SED diagnostic; the flaring of the disc has been left as a free parameter, but it proves no less relevant than full numerical simulations that assume the unknown amount of material in accretion discs at a scale $\la$ $100~\rm AU$; more critical is my leaving aside the self-gravity of the disc, which could "unflare'' the disc and prevent irradiation, but the simulation of self-gravity in the presence of irradiation has not yet been performed, so its influence still remains speculation.

This work can be seen as a feasibility study for a future vertical structure simulation of irradiated discs that I am developing using the radiative transfer formalism presented in Malbet et al. (2001). Such a simulation is needed to issue a spectral diagnostic of irradiation (e.g. silicate feature in emission), that forthcoming IR long-baseline insterferometers will ease with their AU-scale resolution for the closest FUors. In particular, MIDI on the VLTI will be able to measure the silicate feature while disentangling the contribution of an envelope from that from an irradiated disc, which an SED diagnostic cannot do (Vinkovic et al. 2003). One also expects to obtain constraints on the flaring with IR closure phases (AMBER on the VLTI) and on the disc mass with ALMA; their determination is not model-independent and also requires a reliable model.

Acknowledgements
This work has made use of NASA's Astrophysics Data System Bibliographic Services and of CDS's Vizier Catalogue Database. Computations and graphics have been done with free software, in particular Yorick by D. Munro. I also wish to thank C. P. Dullemond for helpful comments that improved the quality of the paper. Language corrections have been suggested by K. Smith.

References

 
Copyright ESO 2004