A&A 419, 203-213 (2004)
DOI: 10.1051/0004-6361:20034318

Constraints in the circumstellar density distribution of massive Young Stellar Objects

C. Alvarez1,2,3 - M. Hoare1 - P. Lucas4


1 - Physics and Astronomy Department, University of Leeds, Leeds LS2 9JT, UK
2 - Kapteyn Astronomical Institute, Postbus 800, 9700 AV Groningen, The Netherlands
3 - Max-Planck-Institut für Astronomie, Königstuhl 17, 69117 Heidelberg, Germany
4 - Dept. of Physical Sciences, University of Hertfordshire, College Lane, Hatfield AL10 9AB, UK

Received 12 September 2003 / Accepted 23 January 2004

Abstract
We use a Monte Carlo code to generate synthetic near-IR reflection nebulae that resemble those (normally associated with a bipolar outflow cavity) seen towards massive young stellar objects (YSOs). The 2D axi-symmetric calculations use an analytic expression for a flattened infalling rotating envelope with a bipolar cavity representing an outflow. We are interested in which aspects of the circumstellar density distribution can be constrained by observations of these reflection nebulae. We therefore keep the line of sight optical depth constant in the model grid, as this is often constrained independently by observations. It is found that envelopes with density distributions corresponding to mass infall rates of $\sim$10-4  ${M_{\odot}~{\rm yr}^{-1}}$ (for an envelope radius of 4700 AU) seen at an inclination angle of $\sim$ ${45}^{\circ }$ approximately reproduce the morphology and extension of the sub-arcsecond nebulae observed in massive YSOs. Based on the flux ratio between the approaching and receding lobe of the nebula, we can constrain the system inclination angle. The cavity opening angle is well constrained from the nebula opening angle. Our simulations indicate that to constrain the outflow cavity shape and the degree of flattening in the envelope, near-IR imaging with higher resolution and dynamic range than speckle imaging in 4 m-class telescopes is needed. The radiative transfer code is also used to simulate the near-IR sub-arcsecond nebula seen in Mon R2 IRS3. We find indications of a shallower opacity law in this massive YSO than in the interstellar medium, or possibly a sharp drop in the envelope density distribution at distances of $\sim$1000 AU from the illuminating source.

Key words: scattering - stars: formation - techniques: high angular resolution - radiative transfer - stars: winds, outflow

   
1 Introduction

Bipolar outflows appear to be a ubiquitous phenomenon during the formation of stars in all mass ranges (Beuther et al. 2002; Henning et al. 2000; Ridge & Moore 2001; Bally & Lada 1983). Low mass young stellar objects (YSOs) show highly collimated bipolar jets from a few 10 AU (Burrows et al. 1996) to several parsec (Eislöffel 2000; Reipurth et al. 1997) in length. These jets are thought to be magneto-hydrodynamically collimated in a wind formed at the inner star-disk system (e.g. X-wind, Shu et al. 1994). The jets are thought to drive the large scale molecular outflow (Masson & Chernin 1994).

The formation and collimation of outflows in massive YSOs is less well understood than in low mass YSOs. There appears to be a lack of highly collimated parsec-scale jets (Mundt & Ray 1994). In the near-IR, searches for shock-excited ${\rm H}_{2}$ show traces of jets in massive star forming regions, but probably driven by low mass young stars located in the same cluster (Wang et al. 2003; Davis et al. 1998). A recent search for optical shock-excited emission in the outer parts of the outflow, yielded no evidence of jet interaction (Alvarez & Hoare, in preparation). Very close to the driving source, there is also not clear evidence that jets are the rule in massive YSOs. In some cases, the free-free radio emission from the inner wind shows a jet morphology (e.g. HH80-81, Marti et al. 1993; Cep A, Torrelles et al. 1996). Such jets would have to be magneto-hydrodynamically driven, even though the OB stars themselves are not magnetically active. Magneto-hydrodynamics in the infalling rotating cloud could set up bipolar flows (Tomisaka 1998). In other cases, the ionised wind appears to be equatorial (e.g. Hoare & Muxlow 1996; Hoare et al. 1994; Hoare 2002). Theoretical models show that radiation pressure in massive young stars can drive gas off the surface of a disk, producing a predominantly equatorial wind (Drew et al. 1998; Drew & Proga 2000). Any initial flow maybe hydrodynamically collimated into a bipolar flow by the flattened surrounding cloud (e.g. Delamarter et al. 2000). These alternative theories will predict different morphologies for the base of the outflow cavities carved out. These variations in morphology occur at scales of a few 100 AU, which at the typical distances to massive YSOs of $\sim$1 kpc, correspond with angular sizes of $\sim$0 $.\!\!^{\prime\prime}$1. Therefore, high resolution techniques are fundamental to study the impact of the the outflow in the surrounding material.

In a related paper (Alvarez et al., in preparation, hereafter Paper I), we show high resolution near-IR speckle images which trace the circumstellar matter around massive YSOs at scales of a few 100 AU. The extended emission that is seen towards some of the sources can be interpreted as scattered light in an outflow cavity due to its monopolar morphology. This interpretation is supported in some cases by the blue colours of the nebula (e.g. Mon R2 IRS3, Paper I, Preibisch et al. 2002). Furthermore, polarimetric speckle imaging of the reflection nebula in the massive star forming region S140 IRS1 (Schertl et al. 2000) shows a centrosymmetric pattern which is typical of scattered light.

Intuitively, one can imagine that depending on the properties of the dust, the shape of the cavity, the density distribution and the orientation of the system with respect to the observer, the resulting reflection nebula will change. The morphology of the cavity is particularly important because it is shaped by the interplay between the infall and the outflow. For instance, it is expected that an equatorial or wide-angled wind will produce a cavity with a wide opening angle near the star. However, a jet is expected to open a rather narrow cavity.

Radiative transfer simulations have been widely used to generate synthetic nebulae that resemble the observations (Wolf et al. 2002; Lucas & Roche 1997; Whitney & Hartmann 1992; Lazareff et al. 1990; Lucas & Roche 1998; Whitney et al. 1997; Fischer et al. 1996,1994; Kenyon et al. 1993). The work by Lazareff et al. (1990) was based on a ray-tracing code and it was focused mainly on the effect produced by different disc models on the synthetic nebulae. The authors compared the general features of the model images with previous seeing-limited images of the low mass systems HL Tau and L 1551 IRS5. Whitney & Hartmann (1992), Kenyon et al. (1993) and Whitney et al. (1997) developed a Monte Carlo code to investigate how the nebula morphology and the near-IR colours of the synthetic images vary with different model parameters. In particular, Whitney et al. (1997) used their code to constrain the colours of the central source, the dust model and the envelope density distribution in a sample of $\sim$20 low mass YSOs. Fischer et al. (1994) and Fischer et al. (1996) developed a new Monte Carlo scattering code and they focused on exploring the effect of different dust models in the synthetic images. Lucas & Roche (1997) and Lucas & Roche (1998) compared synthetic nebulae produced with their Monte Carlo code with high resolution multi-colour observations of reflection nebulae associated with low mass YSOs. From this comparison, they could constrain some parameters defining circumstellar density distribution as well as the dust model for several sources. Recently, radiative transfer Monte Carlo codes have been developed to simulate scattering by non-spherical dust particles Whitney & Wolff (2002), Wolf et al. (2002) and Lucas (2003).

These previous models have focused predominantly on low mass YSOs. Here, we apply the Monte Carlo code of Lucas & Roche (1998) to high mass YSOs, where the infall rates are much higher. We also adopt an observational approach, by presenting a grid of models in which as each parameter is varied, the overall density is scaled too to keep the optical depth along the line of sight constant. This is because the line of sight optical depth is often well constrained from other data such as the optical depth of the 9.7 $\mu$m silicate feature or the colour of the star. The models are decribed in Sect. 2. The grid of models is presented in Sect. 3. In Sect. 4, we use the models to constrain the density distribution in Mon R2 IRS3. Some concluding remarks are shown in Sect. 5.

   
2 Models

We used the Monte Carlo code of Lucas & Roche (1998,1997) with a set of parameters adapted to massive YSOs. The models consist of a central star surrounded by a dusty flattened envelope. A disc characterises the density distribution near the equator. The model also includes an empty bipolar cavity opened by the outflow in the circumstellar matter. The photons emerging from the central source are scattered off the dust grains in the envelope and disc. Each photon can suffer several scattering processes until it is either absorbed, or escapes. All photons traveling in a particular direction are binned, and projected onto the image plane.

The envelope is described by a density distribution resulting from the collapse of a slowly rotating cloud (Ulrich 1976; Terebey et al. 1984). The density ($\rho$) at any point (r, $\mu$) is given by Eq. (1),

 \begin{displaymath}% {\rho(r,\mu)} \sim \frac{\dot{M}}{8\pi r_{\rm c}({\rm GM})^... ...}} \cdot \frac{1}{(1 + \mu_{0})^{1/2}} \cdot \frac{1}{r^{1/2}} \end{displaymath} (1)


 \begin{displaymath}% \mu_{0}^{3} + \left(\frac{r}{r_{\rm c}} - 1\right)\mu_{0} - \frac{r}{r_{\rm c}}\mu = 0 \end{displaymath} (2)

where ${\dot{M}}$ represents the mass infall rate and M is the mass of the central source. ${r}_{{\rm c}}$ is the centrifugal radius, which determines the degree of flattening of the distribution. Flatter density distributions are characterised by larger centrifugal radii. $\mu = \cos \theta$, where $\theta$ is the position angle of each particle with respect to the polar axis. $\mu_0 = \cos \theta_0$ represents the initial position angle of each infalling particle. The equation of motion of the infalling particles (Eq. (2)) should be satisfied at every point (see Ulrich 1976, for details).

The disc plays a passive role in the models presented here. It absorbs and scatters the radiation from the central star, but it does not emit. We use Eq. (3) to describe the disc density structure,

 \begin{displaymath}% \rho_{\rm d}(R,z) = \rho_{0}(R/R_{*})^{-\alpha}{\rm e}^{-\frac{z^{2}}{2H(R)^{2}}} \end{displaymath} (3)

where R and z satisfy the relation r2 = R2 + z2, $\rho_{0}$ represents the density in the midplane at the surface of the star, and R* is the radius of the star. $H(R) = H_{0}(R/R_{*})^{\beta}$ is the disc scale height. $\beta$ parametrises the degree of flaring on the disc. We tried steady Keplerian discs ( $\alpha=15/8$ and $\beta=9/8$) as well as geometrically thin and optically thick discs ( $\alpha=3/4$ and $\beta=0$) (Lazareff et al. 1990; Whitney & Hartmann 1992).

An empty cavity represents the material evacuated by the outflow in the envelope. The shape of the cavity must be determined by the interplay between the infall and outflow processes. However, the lack of knowledge of these processes makes the choice of the cavity shape somewhat arbitrary. For some of the models a conical cavity was used,

 \begin{displaymath}% z = B \cdot (R - R_{{\rm cav}}) \end{displaymath} (4)

where B represents the tangent of the cavity half-opening angle. For other models, a parabolic cavity was used, which is represented by the expression,

 \begin{displaymath}% z = A \cdot (R^{2} - R_{{\rm cav}}^{2}) \end{displaymath} (5)

where the constant A determines the curvature of the parabola, and ${R_{{\rm cav}}}$ represents the radius of the cavity at the equator. The cavity opening angle ( ${\rm 2\theta _{lim}}$) is defined as twice the angle formed by the z axis and the line connecting the intersection between the parabola and the equator with the intersection between the parabola and the outer sphere.

For the dust, we used the mixture of Mathis et al. (1977). The values of the opacities are obtained using a dust to gas ratio by mass of $\sim$10-2. The optical constants and albedos for this mixture were chosen from Draine & Lee (1984) and Draine (1985). The values used for the opacity ${\rm\kappa^{ext}_\lambda}$ are 2.0, 3.8 and  $6.5~{\rm m^{2}~kg^{-1}}$, and the values for the albedo $\omega_\lambda$ are 0.22, 0.34 and 0.45, in the K, H and J bands respectively. The phase function that describes the scattering is within the Rayleigh approximation in the K band, and becomes gradually forward throwing towards the J band.

Table 1: Grid of models.

   
3 Grid of models

In this section, a grid of models is presented (see Table 1) to illustrate how variations on the input parameters affect the morphology of the model images. All models shown have an outer radius ${R_{{\rm out}}}$ = 4700 AU, and are assumed to be located at a distance of 1 kpc. Since we aim to compare the model predictions with typical speckle observations in 4 m-class telescopes, synthetic images are generated with a size of $128\times128$ pixels at a pixel scale of 0 $.\!\!^{\prime\prime}$06. Our choice of  ${R_{{\rm out}}}$ keeps the outer boundary of the models outside the field of view of the synthetic images. The images were convolved with a Gaussian with a FWHM = 0 $.\!\!^{\prime\prime}$2, which is the typical resolution achieved with speckle imaging in 4 m-class telescopes (see Paper I).

The input parameters are varied with respect to a fiducial model (K01, in Table 1). The fiducial model consists of an envelope with a centrifugal radius of 50 AU and a mass infall rate $1.11\times 10^{-4}~{M_{\odot}~{\rm yr}^{-1}}$ (note that ${\dot{M}}$ $\sim$ 10-8-10-5  ${M_{\odot}~{\rm yr}^{-1}}$ are typically inferred for low mass YSOs; e.g., Whitney et al. 1997; Lucas & Roche 1997; Kenyon et al. 1993). The mass infall rates used in our grid are consistent with envelope models in massive stars (Maeder & Behrend 2002; Wolfire & Cassinelli 1986). This mass infall rates, though, would be smaller for models with a larger outer radius. The model also has a geometrically thin but optically thick disc ( $\alpha=3/4$ and $\beta=0$) of radius ${r_{\rm d}}$ = 250 AU and $\rho_{0}=2 \times 10^{-4}$ kg m-3. A conical cavity with a radius at the equator of 100 AU and an opening angle ${\rm 2\theta _{lim}}$ =  ${20}^{\circ }$ represents the effect of the outflow in the envelope. The central source is assumed to be a  $10~{M}_{\odot}$ star with a radius R* = 10  ${R}_{\odot}$. The number of input photons in model K01 is $5\times10^{6}$. Only 8.7% of these photons form part of the output (see Col. 9 in Table 1). The other 91.3% is absorbed either in the disc (9.3% of the total) or in the envelope (79.2%). There is a small fraction (2.8%) that is absorbed by the star itself after being scattered. The optical depth along the line of sight of the fiducial model in the K band ( ${\rm\tau _{\rm los}}$ = 8 in Table 1) corresponds to an extinction in the K band of AK = 8.6 mag. If the extinction law of He et al. (1995) is used, this corresponds to a visual extinction of $\sim$80 mag, i.e. an optical depth of the 9.7 $\mu$m silicate feature of ${\rm\tau_{Si}(9.7~\mu m)}$ $\sim$ 4, (where $\tau_{{\rm Si}}(9.7~\mu {\rm m}) = 0.053 A_{\rm v}$ from Draine & Lee 1984 was used), which is typically observed in massive YSOs.

3.1 Morphology


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi01.eps}\end{figure} Figure 1: K band images for models in which the infall rate has been varied with respect to the fiducial model (K01 in Table 1, see Fig. 2c). a) Model K02 with a  $0.87\times 10^{-4}~{M_{\odot}~{\rm yr}^{-1}}$. b) Model K03 ( $1.06\times 10^{-4}~{M_{\odot}~{\rm yr}^{-1}}$). c) Model K04 ( $1.30\times 10^{-4}~{M_{\odot}~{\rm yr}^{-1}}$). d) Model K05 ( $1.44\times 10^{-4}~{M_{\odot}~{\rm yr}^{-1}}$). The rest of the parameters remain unchanged with respect to model K01. In the four panels, the line of sight forms an angle of ${45}^{\circ }$ with respect to the cavity axis. The 1'' bar represents a length of 1000 AU, and the resolution used is 0 $.\!\!^{\prime\prime}$2. The greyscale varies from -5% (the lightest) to 20% of the maximum brightness (the darkest). The contours are at 5, 10, 15, 30 and 50% of the maximum brightness. These contour levels have been chosen for comparison with typical near-IR speckle imaging in 4 m-class telescopes (see Paper I). ${\rm\tau _{\rm los}}$ represents the optical depth along the line of sight. These and the other model images shown in this section are normalised to the brightness peak, unless otherwise stated.
Open with DEXTER

In this section, we investigate which parameters of the circumstellar density distribution can be constrained from the observed morphology of the nebula. Figure 1 shows the effect produced in the synthetic images by varying the overall optical depth in the envelope. This is done by changing the density scaling through the mass infall rate (${\dot{M}}$) in the Ulrich formula. At low optical depths (Fig. 1a), there is less dust available to scatter the light, and also to absorb the direct light from the star. Therefore, the nebula becomes relatively fainter than at higher optical depths, and the only contribution to the image at a 5% level is the central star. However, at higher optical depths, the nebula becomes much more extended, since the central star becomes highly obscured and there is more dust available for scattering. Hence, only a change of 50% in the mass infall rate has a dramatic change on the appearance of the system, due mainly to the exponential dependence on the line of sight optical depth of the central star brightness. Therefore, from an observational point of view, the line of sight optical depth is the most important parameter.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi02.eps}\end{figure} Figure 2: Synthetic K band images of reflection nebulae for the models at different inclination angles ($\iota $). The mass infall rate has been adjusted to yield the same optical depth along the line of sight ( ${\rm\tau _{\rm los}}$ = 8) at all four inclination. ${\dot{M}}$ is expressed in units of 10-4  ${M_{\odot}~{\rm yr}^{-1}}$ and $\iota $ in degrees. Panel  a) corresponds with a near-edge-on model (K06 in Table 1). The image in panel  b) corresponds with an intermediate inclination (model K07). Panel  c) shows the fiducial model (K01), which is seen at an inclination angle of  ${45}^{\circ }$. d) Model with a near-pole-on inclination (K08). The contours and the greyscale in all panels are defined as in Fig. 1.
Open with DEXTER

Figure 2 shows the K band images for different inclination angles of the line of sight with respect to the system axis (i.e. cavity axis). The overall density scaling has been adjusted in each of the four models to yield the same optical depth along the line of sight as the fiducial model ( ${\rm\tau _{\rm los}}$ = 8). For views near edge-on (panel a in Fig. 2), the nebula is clearly bipolar at the 5% level. The receding lobe appears less bright than the approaching lobe. In this case (model K06 in Table 1), a mass infall rate 2.3 times lower than for the fiducial model was used. Even at this near-edge-on inclination, it is possible to see the central star due to the low overall density scaling. At intermediate inclinations (panels b and c), the receding lobe is not detected any more at the 5% level. At low inclinations, (panel d) a faint monopolar nebula can still be seen.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi03.eps}\end{figure} Figure 3: a) Variation of the approaching to receding lobe flux ratio with system inclination angle, where the line of sight optical depth is kept constant. The system inclination angle varies from near-edge-on ($\mu = 0.1$) to near-face-on ($\mu = 0.9$). b) Variation of the approaching to receding lobe flux ratio with overall opacity.
Open with DEXTER

The contrast between the approaching and receding nebula ( ${F_{{\rm app}}/F_{{\rm rec}}}$) is a useful quantity to estimate the system inclination angle. Figure 3a shows the variation of the ${F_{{\rm app}}/F_{{\rm rec}}}$ with the inclination angle. Each point in the plot represents a model whose overall density scaling has been adjusted to yield an optical depth along the line of sight of 8 at its corresponding inclination. The flux ratio has been calculated using aperture photometry with an aperture radius of 1'' in the synthetic images. The aperture on the approaching lobe was centered at 0 $.\!\!^{\prime\prime}$8 from the star along the cavity axis and it includes the star itself. The aperture on the receding lobe was centered at 1 $.\!\!^{\prime\prime}$2 from the star also along the cavity axis and it does not include the star. This avoids sensitivity to the actual location of the apertures. For near-edge-on views (small $\mu$'s), the approaching and receding lobe have roughly the same brightness. Therefore, the ${F_{{\rm app}}/F_{{\rm rec}}}$ is nearly 1. As $\mu$ increases, the approaching lobe becomes relatively brighter, and hence the value of  ${F_{{\rm app}}/F_{{\rm rec}}}$ increases. The figure also shows that the ratio ${F_{{\rm app}}/F_{{\rm rec}}}$ becomes 100 for values of $\mu$ in the range between 0.70 and 0.80 (i.e. for inclination angles in the range ${35}^{\circ}$ and  ${45}^{\circ }$). Hence, a 1% upper limit in the detection of the receding lobe indicates that the system is seen under an inclination angle $\le$ ${45}^{\circ }$.

Figure 3b shows the change of ${F_{{\rm app}}/F_{{\rm rec}}}$ due to variations in the optical depth (i.e. variations in the ${\dot{M}}$) for an inclination angle of  ${45}^{\circ }$. At very low optical depths ( ${\rm\tau _{\rm los}}$ $\la 5$) the star dominates the flux. As the opacity increases ($6 \la$  ${\rm\tau _{\rm los}}$ $\la 9$), the star is increasingly obscured, and the counter-lobe starts to show up. At larger optical depths ( ${\rm\tau _{\rm los}}$ $\ga 10$), the approaching lobe starts to dominate the emission, while the receding lobe hardly changes. Hence, the ratio ${F_{{\rm app}}/F_{{\rm rec}}}$ increases again. The contrast between the nebula lobes is less sensitive to changes in the mass infall rate than to changes in the inclination angle. Hence, in principle, the approaching to receding lobe flux ratio can be used to constrain the inclination angle.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi04.eps}\end{figure} Figure 4: K band images resulting from models with different centrifugal radius seen at an inclination angle of  ${45}^{\circ }$. The mass accretion rate for each model was set to keep the ${\rm\tau _{\rm los}}$ = 8 along the line of sight at  ${45}^{\circ }$ in the K band. The rest of the parameters remain unchanged with respect to the fiducial model. In all panels, the value of  ${r}_{{\rm c}}$ is expressed in AU and ${\dot{M}}$ is expressed in units of 10-4  ${M_{\odot}~{\rm yr}^{-1}}$. a) Model K09 in Table 1, which represents a spherical envelope (i.e. the limiting case of a negligible centrifugal radius). b) Model K10, which has a ${r}_{{\rm c}}$ = 100 AU. c) Model K11, which has a  ${r}_{{\rm c}}$ = 150 AU. d) Model K12 with a  ${r}_{{\rm c}}$ = 200 AU. The contours and the greyscale in all four panels are defined as in Fig. 1.
Open with DEXTER

In Fig. 4, we investigate whether the shape of the envelope can be derived. The images show models with an increasing centrifugal radius (i.e. flattness of the envelope), from panel a) to panel d). The mass infall rate for each model has again been set such that the optical depth along a line of sight at  ${45}^{\circ }$ remains the same as for the fiducial model (i.e. ${\rm\tau _{\rm los}}$ = 8 in the K band). The fiducial model, with a centrifugal radius of 50 AU, is shown in Fig. 2c. Figure 4a represents the case of a spherical density distribution. The nebula becomes less bright as the centrifugal radius increases (panels b-d in Fig. 4). Since the material is predominantly concentrated on the equatorial plane and there is less dust available in the polar regions of the envelope, where a large fraction of the scattered light is generated. However, this is a very subtle change compared to that for the line of sight optical depth o inclination angle. Hence, the reflection nebula tells us little about the degree of flattening of the envelope.

Figure 5 shows the effect of varying the cavity opening angle ( $10^{\circ} \le 2\theta_{{\rm lim}} \le 50^{\circ}$) on the synthetic nebulae. The mass infall rate was increased with the cavity angle to keep the same optical depth along the line of sight in all images. Unsurprisingly, the nebula opening angle appears to be larger for models with a wide-angled cavity. The significant changes on the nebula shape occur for the external contours, while the inner contours remain nearly unchanged. This is consistent with the fact that variations in the cavity opening angle will affect more the external regions of the envelope than the regions close to the cavity base. At wider cavity opening angles, more stellar photons can escape the system without being scattered. This is shown in Col. 9 of Table 1, where the fraction of photons (scattered and stellar) that leaves the system ( ${\rm\mathcal{F}_{out}}$) is listed.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi05.eps}\end{figure} Figure 5: K band images corresponding to models with different cavity opening angles seen at an inclination angle of ${45}^{\circ }$. The ${\dot{M}}$ for each model (which is shown in units of 10-4  ${M_{\odot}~{\rm yr}^{-1}}$) has been set to yield a K band optical depth of 8 along the line of sight. The conical cavity opening angle increases from panel  a) ( ${10}^{\circ }$) to panel  d) ( ${50}^{\circ }$) (i.e. models from K13 to K16 in Table 1). The greyscale and contours are as in Fig. 1.
Open with DEXTER

The result of varying the cavity radius at the equator from 50 to 250 AU is shown in Fig. 6. A radius of 50 AU corresponds approximately to the sublimation radius for dust in OB stars. Radiation pressure or wind interactions could increase the size of the hole. The nebula appears fainter and slightly wider for larger values of the ${R_{{\rm cav}}}$. In this case, the differences between nebulae can also be observed in the innermost contours. This is caused by the fact that changes in ${R_{{\rm cav}}}$ will have a stronger impact in the regions of the envelope closer to the equator than in the outer regions. The total number of output photons increases by a factor of 1.4 from the model with ${R_{{\rm cav}}}$ = 50 AU to the model with ${R_{{\rm cav}}}$ = 250 AU (see Col. 9 in Table 1). For larger values of  ${R_{{\rm cav}}}$, similar number of photons are scattered of a larger surface. Hence, the surface brightness of the nebula decreases from Figs. 6a to d.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi06.eps}\end{figure} Figure 6: Variation of the nebula shape with the cavity radius at the equator (models K17 to K20 in Table 1). The mass infall rate was adjusted to yield the same optical depth along the line of sight in all four images. The values for  ${R_{{\rm cav}}}$ are expressed in AU. All panels represent a line of sight at  ${45}^{\circ }$. The greyscale and contours are as in Fig. 1.
Open with DEXTER

We now investigate the effect of the detailed shape of the base of the cavity. A parabolic cavity might be expected if the central wind is initially equatorial, whilst a conical cavity would arise from a jet-driven flow in low mass stars (Bachiller et al. 1995). Figure 7a shows the K band image for a model with a parabolic cavity compared with a conical nebula in Fig. 7b. The other parameters have the same values as for the fiducial model, except the overall opacity scaling, which has been enhanced by a factor of 1.4 to yield the same optical depth along the line of sight as in the fiducial model. The parabolic shape is clearly seen in the resulting nebula, as expected (Fig. 2). A larger fraction of photons escape from the parabolic model than from the conical model because the parabolic cavity is broader near the star than the conical cavity. The concave shape of the cavity walls, favours the scattering in the outer regions of the envelope, which also contributes to make the parabolic nebula more extended.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi07.eps}\end{figure} Figure 7: a) Model with a parabolic cavity (K26 in Table 1) compared with the fiducial model at an inclination of  ${45}^{\circ }$. The cavity opening angle is the same as for the conical cavity in the fiducial model ( ${\rm 2\theta _{lim}}$ =  ${20}^{\circ }$). b) Fiducial model at an inclination of  ${45}^{\circ }$. This is the same image as shown in Fig. 2c. The greyscale and contours are as in Fig. 1. The density scaling is a factor of 1.4 larger in the model with the parabolic cavity than in the conical cavity model.
Open with DEXTER

All models shown up to this point included an optically thick, and geometrically thin flat disc of radius ${r_{\rm d}}$ = 250 AU. To investigate any dependence of the model images on the disc radius, we have generated models with the same parameters as the fiducial model (K01, Fig. 2a) but with varying disc radius ( $100~{\rm AU} \le r_{\rm d} \le 500~{\rm AU}$). No significant differences were found between the images resulting from these models and the fiducial model. We also investigated whether the introduction of a flaring disc may change the morphology of the reflection nebulae. Model K25 in Table 1 has the same parameters as the fiducial model except for a flaring in the inner disc. The values used for $\alpha$ and $\beta$ in the disc equation (Eq. (3)) were 15/8 and 9/8 respectively (i.e. a Keplerian disc). No relevant differences were appreciated between these images and the fiducial model. The reason is that the disc flaring angle for the value of $\beta=9/8$ is small compared with the line of sight inclination angle. The envelope density dominates the disc density for all inclinations $\le$ $77^{\circ}$ ( $\mu \ge 0.22$).

   
3.2 Dependence with wavelength


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi08.eps}\end{figure} Figure 8: Multi-band models. The two panels represent the fiducial model in the J and H bands, seen at an inclination angle of  ${45}^{\circ }$. In Table 1 they are labeled as J01 (panel  a)), H01 (panel  b)). $\tau $ is the optical depth at an inclination of  ${45}^{\circ }$. The greyscale and the contours are the same as in Fig. 1, and the grey cross indicates the position of the star.
Open with DEXTER

Figures 2c, 8a and 8b show the variation of the synthetic images with the wavelength (J, H and K band respectively). The geometry and density distribution for the J and H band models are the same as for the fiducial model. However, the opacities, the albedos and the scattering matrix change with wavelength. The opacities correspond with an extinction law ( ${\rm\kappa^{ext}_{\lambda} \propto \lambda^{-\gamma}}$) with an exponent $\gamma = 2.3$. The number of photons in the input spectrum at H ( $\mathcal{N}_{{H}}^{{\rm in}}$) is a factor of 2.00 times the number of photons at K ( $\mathcal{N}_{{K}}^{{\rm in}}$), and the number of photons at J ( $\mathcal{N}_{{J}}^{{\rm in}}$) is 2.45 times the number of photons at K. The number of input photons at each wavelength was calculated using the model stellar atmospheres from Kurucz (1979). The spectral energy distribution (SED) of model OB main sequence stars were integrated in the J, H and K bands using the transmission profiles of the filters J98, H98 and K98 at UKIRT. The ratio between the number of photons emitted at two given bands approaches asymptotically to a constant value for earlier spectral types. This limit corresponds with a slope in the stellar SED of -2.4 ( $\mathcal{N}_\lambda^{{\rm in}} \propto \lambda^{-2.4}$), which is shallower than the theoretical Rayleigh-Jeans limit ( $\mathcal{N}_\lambda^{{\rm in}} \propto \lambda^{-3}$). This shallower SED is a better representation of the colours of OB main sequence stars than the Rayleigh-Jeans approximation (see the UKIRT web-page: http://www.jach.hawaii.edu/JACpublic/UKIRT/). Figures 2c, 8a and 8b show that the nebula appears more extended at short wavelengths, while the star becomes totally obscured, due to the increase in the opacity. It can also be seen that the separation between the star and the nebula apex decreases towards longer wavelengths, i.e. it is possible to probe the circumstellar density distribution closer to the star at longer wavelengths.

We now address the question of how the colours of the synthetic images vary with the different model parameters. The H-K colour for the models can be estimated using the following expression,

 \begin{displaymath}% (H - K)_{{\rm mod}} = -2.5\log\left(\frac{\mathcal{N}^{{\rm... ...ac{\lambda_{K}}{\lambda_{H}}\frac{f_{K}^{0}}{f_{H}^{0}}\right) \end{displaymath} (6)

where the $\mathcal{N}$'s represent the number of photons within a certain aperture on each band. The $\lambda$'s are the central wavelengths of the filters, and the f0's are the zero magnitude fluxes. The values used for $\lambda$ are 1.65 and  $2.20~{\rm\mu m}$ at H and K respectively, and $1.12\times10^{-9}$ W m$^{-2}~\mu$m-1 and $4.07\times10^{-10}$ W m$^{-2}~\mu$m-1 are the values for fH0 and fK0 respectively (UKIRT web-page: http://www.jach.hawaii.edu/JACpublic/UKIRT/). For simplicity, the colours shown here assume that the foreground extinction is negligible. Note that the presence of unknown foreground extinction is an added complication when models are compared with observations.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi09.eps}\end{figure} Figure 9: a) Variation of the model H-K colour (Eq. (6)) with the line of sight inclination angle for the fiducial model. b)- d) The inclination angle is fixed at  ${45}^{\circ }$ and the overall opacity, centrifugal radius and cavity opening angle are varied. In panels  c) and  d) the mass infall rate was also changed to keep the opacity along the line of sight constant.
Open with DEXTER

Figure 9a shows the change in the H-K colour with the line of sight inclination angle for the fiducial model. An aperture radius of 1'' centered at the image centre was used. The H-K colour at all inclinations, except for the very close to pole-on, is clearly redder than the H-K colour of the input spectrum ((H-K) $_{{\rm inp}}$ = 0.04). At edge-on inclinations ( $\mu \sim 0$), no star is seen in either band, and the nebula at K is less extincted than the nebula at H. At intermediate inclinations, the nebula at H becomes brighter because the stellar light passes through a less dense part of the envelope. The star begins to appear at K but is not yet seen at H. The overall effect is that the H-K colour becomes bluer. For $\mu > 0.65$ the star brightens quickly at K, dominating the flux in this band, while the star just begins to appear at H. Hence, the H-K colour becomes redder. At face-on views (i.e. $\mu \ge 0.94$) the H-K colour tends to the (H-K) $_{{\rm inp}}$ because the star is now seen directly through the cavity. The behaviour of the H-K colour with the inclination angle in our simulations is different to that found by Kenyon et al. (1993) and Whitney et al. (1997). Their models become bluer at high inclination angles, while ours become redder due to the higher envelope optical depth (i.e. higher mass infall rate) and due to the fact that, at high inclinations, no star is seen in the H band.

Figure 9b illustrates how the H-K colour varies with the overall opacity at a fixed inclination angle of  ${45}^{\circ }$. At large opacities, the H-K colour is bluer because the reflection nebula becomes relatively brighter than the star. In Fig. 9c, we represent models with different centrifugal radii seen at a fixed inclination of  ${45}^{\circ }$. The mass accretion rate was adjusted so that all of them have an optical depth along the line of sight of 8. In the models with a larger centrifugal radius (i.e. a flatter density distribution) there is less dust available in the outer parts of the envelope, since the dust is mainly concentrated at the equator. This density enhancement in the equatorial region allows more direct light from the central star to escape (predominantly in the K band due to the low opacity) than in models with a small  ${r}_{{\rm c}}$. The consequence is that the H-K colour becomes slightly redder for models with large ${r}_{{\rm c}}$.

Figure 9d shows variations of the H-K colour with the cavity opening angle. All models represented in this plot are seen at an inclination of  ${45}^{\circ }$, and the overall density scaling was adjusted to yield an ${\rm\tau _{\rm los}}$ = 8. The general trend is that at wide cavity opening angles the H-K colour becomes bluer because more scattered photons in the outer parts of the envelope (at H, not at K) can escape from the system. This general trend is favoured by the fact that the direct light from the star is equally extincted at all cavity opening angles shown in Fig. 9d. The figure also illustrates that the H-K colour for models with an ${\rm 2\theta _{lim}}$ >  ${50}^{\circ }$ becomes even bluer than the H-K colour of the input spectrum. In summary, the four panels in Fig. 9 show that the change in the H-K colour is not higher than 1 mag within the ranges of the parameter space investigated. Hence, the integrated colour is not particularly sensitive to the density distribution, and is unlikely to yield unique solutions for the model parameters.

   
3.3 Dust in the cavity

Up to now, we have assumed that the cavity evacuated by the outflow is empty. However, in a more realistic situation, one would expect that some dust may remain inside the cavity. We estimated the expected  ${A}_{{\rm v}}$ due to dust within the cavity from the typical ${\rm H}_{2}$ column densities in massive outflows as follows. Ridge & Moore (2001) show that the integrated CO (J=1-2) intensity in molecular outflows from 11 massive YSOs is, on average, $\sim$200 K  ${\rm km~s}^{-1}$. If we use the expression $N({\rm H}_2)~({\rm cm}^{-2}) = 3$ $\times$ $10^{20}\int{T_{{\rm mb}}~{\rm d}v}$, where $\int{T_{{\rm mb}}~{\rm d}v}$ represents the integrated CO (J=1-2) intensity (Osterloh et al. 1997; Henning et al. 2000), typical ${\rm H}_{2}$ column densities in the outflow of $6\times10^{22}$ cm-2 are obtained. This yields typical visual extinctions of $\sim$3.5 mag, if the expression $A_{\rm v}/N({\rm H}_2)~({\rm cm}^{-2}) = 5.88$ $\times$ 10-23from Bohlin et al. (1978) is used. This implies an $A_{K} \sim 0.4$ mag which is about 200 times smaller than the AK = 8, used in our fiducial model.

To see the effect that dust in the cavity has in our simulations, we have run three models with the same parameters as our fiducial model but with an AK along the cavity axis of 0.1, 1 and 10 mag. A dusty density distribution within the cavity given by the expression $\rho(z) = (1 + \vert z\vert/R_{{\rm cav}})^{-2}$, which is expected in a disc wind with a constant flow density distribution in a conical cavity, was chosen. The output from the models with AK = 0.1 and 1 are hardly distinguishable from the fiducial model. For AK = 10, only a few photons from the star can scape the system, and the reflection nebula is not seen at all. However, the latter value of AK is far larger than our previous estimate of extinction due to dust in the cavity. Therefore, for the typical extinction expected due to dust within the cavity in massive YSOs, no noticeable effect on the observable properties of our scattering models is detected. Note that this result depends on the selected shape for the density distribution within the cavity. If an optically thin uniform density distribution is chosen, the resulting reflection nebula is expected to be more extended than for an empty cavity (see Lucas & Roche 1996).

   
4 Model for Mon R2 IRS3 S

In Paper I, multi-colour (H and K band) speckle images of a pair of outflow cavities in Mon R2 IRS3 were presented. In this section, the data for IRS3 S are used to find observational constraints to the density distribution using the models presented in the previous section. In Figs. 10a and b, we show the H and K band images of Mon R2 IRS3 S shown from Paper I, with the presumed cavity axis oriented along the vertical direction. These images have a pixel scale of 0 $.\!\!^{\prime\prime}$057 and a resolution of 0 $.\!\!^{\prime\prime}$19, which is comparable to the values used the simulations. The procedure used to find a good fit model for these data is the following. A sub-set of models that match the morphology of IRS3 S at K are chosen amongst the grid shown in Sect. 3. A model is searched amongst these that also matches the morphology of the source in the H band. Finally, the H-K colour of the model is compared with the observations.


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi10.eps}\end{figure} Figure 10: a), b) H and K speckle images of Mon R2 IRS3 S. These images were obtained by rotating the images shown in Paper I by an angle of  ${198}^{\circ }$. The orientation is such that the presumed cavity axis is along the vertical axis. The greyscale varies from -5% (white) to 40% (black) of the brightness peak. The contour levels are at 5, 10, 15, 30 and 50% of the peak. The images are normalised to the peak brightness. c), d) Images for models K01 and HH01 seen at an inclination angle of  ${45}^{\circ }$. These images are also normalised to their brightness peak. The resolution, contours and greyscale are the same as on panels  a) and  b). The cross in the H band images indicate the location of the embedded star.
Open with DEXTER

Firstly, a reasonable match to the K band image for IRS3 S was found, which is given by the fiducial model used in the grid seen at an inclination angle of  ${45}^{\circ }$ (Fig. 2c). A zoom into the central 2'' of this model is shown in Fig. 10d. The H band image for the fiducial model using the opacity ${\kappa^{{\rm ext}}_H}$ = $3.8~{\rm m^{2}~kg^{-1}}$ (cf. Sect. 2) appears to be too extended (see Fig. 8b). Other combinations of inclination angle and mass infall rate were also too extended since the optical depth at H was always too high. Therefore, a shallower opacity law was used in an attempt to achieve a better fit to the IRS3 S H band image. The new value used for ${\rm\kappa^{ext}_\lambda}$ at H was $2.9~{\rm m^{2}~kg^{-1}}$ which corresponds with an opacity law ( $\kappa_\lambda \propto \lambda^{-\gamma}$) with $\gamma = 1.3$ between $1.65~{\rm\mu m}$ and $2.2~{\rm\mu m}$ (instead of the previous value of $\gamma = 2.3$). The exponent 1.3 was chosen to yield an opacity law slightly flatter than than the $\gamma = 1.7$inferred from studies on the interstellar extinction in the near-IR (e.g. He et al. 1995). The model using the shallower opacity law (HH01 in Table 1) is shown in Fig. 10. It can be seen that this law still yields a rather extended nebula compared with the H band image of IRS3 S. A shallower opacity law than the one given by Draine & Lee (1984) has been used by other authors (e.g. Lucas & Roche 1998) to explain the small variation with wavelength of reflection nebulae in low mass YSOs. This is possibly due to a different dust composition or grain size distribution in the circumstellar matter of YSOs than in the interstellar medium. Another possible explanation relies on the fact that speckle imaging acts as a spatial filter for diffuse extended features of the same order and larger than the seeing. Therefore, even if the nebula in the H band was extended for a few arcseconds, the reconstructed speckle image would only pick up the structure closer to the nebula peak, which is dominated by the high spatial frequencies. It may also be the case that the envelope density distribution is truncated at a distance of $\sim$1000 from the central source, causing both H and K images to have the same extent.

Even though a reasonably good fit was found to the K band speckle image of IRS3 S, the synthetic nebula shows some differences with respect to the observed nebula. In particular, the observed nebula departs clearly from the axi-symmetry, probably due to foreground extinction or the presence of a clumpy envelope. Besides, the brightness of the synthetic nebula (Fig. 10d) drops faster near the star than in the case of the observed nebula (Fig. 10b). Further out, at $\sim$1'' from the star, the brightness of the synthetic nebula appears to drop more smoothly than observed in IRS3 S. The sub-arcsecond nebula is the inner region of a large reflection nebula (with a diameter of $\sim$15'') that is seen in the near-IR seeing limited images of Mon R2 IRS3 (see Yao et al. 1997; Aspin & Walther 1990). The large scale nebula is a halo of reflected light around the whole cloud but its polarisation pattern indicates enhanced scattering along an axis that coincides with the sub-arcsecond nebula axis (PA $\sim$  ${198}^{\circ }$) .


  \begin{figure} \par\includegraphics[width=8.8cm,clip]{0318fi11.eps}\end{figure} Figure 11: K band profile along the cavity axis for IRS3 S. The speckle data (filled triangles) are obtained from Fig. 10b. The solid circles represent the data from the contour map shown by Aspin & Walther (1990). The solid line corresponds with an extension of the fiducial model`up to a radius of ${9.4}\times 10^{4}$ AU. The speckle data and the model are flux calibrated following the procedure explained in the text. The offsets are in arcseconds from the position of the star. The errors for the speckle data are smaller than the size of the plotting symbol.
Open with DEXTER

To investigate if the model that matches the sub-arcsecond morphology of IRS3 S fits the outer part of the reflection nebula, the outer radius of the fiducial model was doubled while the mass infall rate was slightly reduced to  $1.03\times 10^{-4}~{M_{\odot}~{\rm yr}^{-1}}$ to keep a K band optical depth of 8 at an inclination of  ${45}^{\circ }$. The profile along the cavity axis for this model was compared with the profile for the K band reconstructed image of IRS3 S and with the profile along the same direction for the reflection nebula studied by Aspin & Walther (1990). This comparison is shown in Fig. 11. The data for the large scale nebula were taken from the flux-calibrated contour map shown in Fig. 3 of Aspin & Walther (1990). The error bars correspond to half of the contour level separation in their map ( ${\rm0.25\ mag/arcsec^2}$). The data points for the observed sub-arcsecond nebula were calculated from a 0 $.\!\!^{\prime\prime}$2 wide line along the cavity axis passing through the star in IRS3 S on the normalised speckle image. The zero point offset for the magnitude scale was calculated by assigning a K magnitude of 6.6 (obtained from the 2MASS survey) to the total number of counts within an aperture of radius 2 $.\!\!^{\prime\prime}$5 that included both IRS3 N and S. The statistical error in the surface brightness is $\la$ $0.2~{\rm mag/arcsec}^2$.

The profile for the model was calculated from 0 $.\!\!^{\prime\prime}$2 wide line along the cavity axis in the synthetic image normalised to the brightness peak. The zero point offset for the magnitude scale was chosen such that it yields the same brightness at the peak as the flux-calibrated speckle data.

Figure 11 shows that there may be a discrepancy of at least ${\rm 1\ mag/arcsec^2}$ between the data from Aspin et al. and the speckle data. A possible explanation is that speckle imaging filters out any background smooth emission, and hence less counts in the speckle image than in the seeing limited image of Aspin et al. appear to account for the same flux. The profiles plotted in Fig. 11 show that the model predicts emission from the outer parts of the nebula (at a distance from the star >1'') that is about 2 mag brighter than it is observed. This may be caused by a drop in the density distribution of IRS3 at a radius $\sim$1000 AU that is not included in the models (cf. Figs. 10a and b).

Willner et al. (1982) estimated an optical depth for the silicate absorption feature at  $9.7~{\rm\mu m}$ of ${\rm\tau_{Si}(9.7~\mu m)}$ = 4.30, using a 27'' aperture centered at Mon R2 IRS3. The same value is also inferred from the ISO spectrum of the region (Jackie Keane private communication). This value can be used as an estimate of the total extinction (foreground and through the envelope) towards IRS3 S. This ${\rm\tau_{Si}(9.7~\mu m)}$ leads to an optical extinction of ${A}_{{\rm v}} \sim 82~{\rm mag}$ (Draine & Lee 1984) . This corresponds with an extinction in the K band AK = 8.6 mag (i.e. $\tau_{{K}} \sim 7.9$) using the reddening law of He et al. (1995), which is in good agreement with the K band optical depth for the fiducial model ( ${\tau_{K} = 8}$).

We found that the H-K colour of IRS3 S is redder than the H-K colour of any of the model images. This is probably caused by the fact that we used the SED of an OB main sequence star to define the input spectrum. However, it is very likely that the illuminating source has an excess in the K band due to the presence of a possible accretion disc, and/or hot dust near the star. The inclusion of an excess of K band photons in the input spectrum may solve the discrepancy between the observed and modeled colours. High resolution photometry at longer wavelengths is required to constrain the input spectrum in the simulations.

   
5 Conclusions

Radiative transfer Monte Carlo simulations have been used to investigate the density distribution in massive YSOs at scales where the outflow is generated. The assumed density distribution consists of a central massive star within a flattened dusty envelope, with a cavity and an inner optically thick disc. It is found that envelopes with density distributions corresponding to typical mass infall rates of $\sim$10-4  ${M_{\odot}~{\rm yr}^{-1}}$ seen at an inclination angle of $\sim$ ${45}^{\circ }$ approximately reproduce the morphology and extension of the sub-arcsecond nebulae observed in massive YSOs. The inclination angle can be constrained by the measurement of the contrast between the approaching and the receding nebular lobe, although observations with a high dynamic range are required (e.g. adaptive optics). The cavity opening angle is well constrained by the nebula opening angle. The simulations indicate possibly some constraints on cavity shape and radius at the equator, which could have implications for the initial angle of the outflow (e.g. jet, wide-angled, equatorial). However, higher resolution than provided by speckle imaging in 4 m-class telescopes is needed to achieve better constraints of these two quantities. The models do not provide significant constraints on the flattening in the envelope or the size of the equatorial disc, which require direct observations with millimetre interferometry.

The Monte Carlo code was also applied to the near-IR sub-arcsecond reflection nebula seen in Mon R2 IRS3 S. An envelope with a mass infall rate of 10-4  ${M_{\odot}~{\rm yr}^{-1}}$ that includes a conical cavity with an opening angle of  ${20}^{\circ }$ seen at an inclination angle of  ${45}^{\circ }$ provides a reasonable match for the K band image. However, no set of input parameters was found that reproduces both the H and K band images of IRS3 S. An opacity law with an exponent $\gamma = 2.3$ (Draine & Lee 1984; Draine 1985) yields H band nebulae that are too extended with respect to the observations. This would also be the case for the observed interstellar extinction law ( $\gamma = 1.7$; He et al. 1995). A shallower opacity law ( $\gamma = 1.3$) yields a better match to observed H band nebula, although still too extended. This indicates that the dust in the circumstellar envelope of massive YSOs may have a rather different optical properties to the dust that forms part of the interstellar medium. However, a truncated density distribution could also explain the data.

Overall, this work shows that future high resolution ($\sim$0 $.\!\!^{\prime\prime}$05) high dynamic range (>100) near-IR imaging has the potential to constrain the inclination angle and shape of the base of the outflow cavity. In turn, this could test hydrodynamic models of the interplay between the infall and outflow in massive YSOs.

Acknowledgements
C.A. would like to thank to the Physics and Astronomy Department at Leeds University for their support. C.A. is also deeply grateful to Kapteyn Astronomical Institute for allowing him to use their facilities during the realization of this work. We would like to thank B. Whitney, the referee of this work, for her comments and suggestions.

References


Copyright ESO 2004