Radiation Transfer of Models of Massive Star Formation. II. Effects of the Outflow
We present radiation transfer simulations of a massive () protostar forming from a massive () protostellar core, extending the model developed by Zhang & Tan (2011, Paper I). The two principal improvements are (1) developing a model for the density and velocity structure of a disk wind that fills the bipolar outflow cavities, based in part on the disk-wind model of Blandford & Payne (1982); and (2) solving for the radially varying accretion rate in the disk due to a supply of mass and angular momentum from the infall envelope and their loss to the disk wind. One consequence of the launching of the disk wind is a reduction in the amount of accretion power that is radiated by the disk. We also include a non-Keplerian potential appropriate for a growing, massive disk. For the transition from dusty to dust-free conditions where gas opacities dominate, we now implement a gradual change as a more realistic approximation of dust destruction. We study how the above effects, especially the outflow, influence the spectral energy distributions (SEDs) and the synthetic images of the protostar. Dust in the outflow cavity significantly affects the SEDs at most viewing angles. It further attenuates the short-wavelength flux from the protostar, controlling how the accretion disk may be viewed, and contributes a significant part of the near- and mid-IR fluxes. These fluxes warm the disk, boosting the mid- and far-IR emission. We find that for near face-on views, i.e. looking down the outflow cavity (although not too close to the axis), the SED from the near-IR to about 60 m is very flat, which may be used to identify such systems. We show that the near-facing outflow cavity and its walls are still the most significant features in images up to 70 m, dominating the mid-IR emission and determining its morphology. The thermal emission from the dusty outflow itself dominates the flux at m. The detailed distribution of the dust in the outflow affects the morphology, so the model can be constrained by considering detailed intensity profiles along and perpendicular to the outflow axis. For example, even though the outflow cavity is wide, at 10 to 20 m, the dust in the disk wind can make the outflow appear narrower than in the near-IR bands.
Subject headings:ISM: clouds, dust, extinction, jets and outflows — stars: formation
There is still no consensus on how massive stars form. Possible mechanisms include collapse from massive cores supported by turbulence and magnetic fields, as in the Turbulent Core model (McKee & Tan 2003, hereafter MT03), competitive accretion at the centers of forming star clusters (Bonnell et al. 2001), and stellar collisions (Bonnell et al. 1998). Theoretically, if a massive star forms similarly to its low-mass counterpart, i.e. from a core with a typical ratio of rotational to gravitational energy of a few percent (Goodman et al. 1993), a massive accretion disk with a size of several hundreds to 1000 AU and bipolar outflows are expected to appear around the massive protostar (Tan 2005). Such disks and outflows have been seen in hydrodynamic simulations with varying implementations of magnetic field, turbulence and radiation (e.g., Krumholz et al. 2007, Kuiper et al. 2010, Seifried et al. 2011, 2012, Peters et al. 2011, Vaidya et al. 2011). However, resolving this issue observationally is difficult as massive stars tend to form in distant, crowded, and highly obscured environments. Hot cores with collimated molecular outflows and disk-like structures have been found around massive young stars, but so far these observed disk-like structures (or toroids) are much larger and more massive than expected for the rotationally-supported disk of an individual protostar, and there is little direct evidence for Keplerian or near-Keplerian disks (e.g., Cesaroni et al. 2007). The outflows have also been observed in near- and mid-IR continuum emission, and it has been argued that the outflow cavity may dominate the mid-IR morphology (De Buizer 2006).
In order to better understand the properties of these observed massive young stellar objects (MYSOs), a number of radiative transfer models have been developed to compare with observations. One of the most popular is the model grid by Robitaille et al. (2006), but it is focussed on low-mass star formation. Massive star-forming cores usually reside in clumps within giant molecular clouds with large mass surface densities, (MT03), so they are themselves relatively dense and collapse with high accretion rates. The expected properties of the fiducial cores of MT03 lie outside the parameter space of the Robitaille et al. (2006) grid. Also the components in this radiative transfer model are relatively simple — e.g., the density in the outflow is usually assumed to be either a constant value or a power-law. Offner et al. (2012) compared this model grid with results of numerical simulations of low-mass star formation, and found possible discrepancies between the parameters inferred from fitting with the radiative-transfer model grid and those of the protostars in the simulation, especially for deeply embedded objects. Other radiative-transfer model grids, such as those of Molinari et al. (2008), are for MYSOs, but the components are also relatively simple. So it is worth including a more realistic density distribution in the outflow cavity and constructing a fully self-consistent model for all the components of the star-forming core.
This is the second of a series of papers on modeling the spectral energy distributions (SEDs) and images of massive star-forming cores. In the previous paper (Zhang & Tan 2011, hereafter Paper I), we studied a core with an initial mass of 60 that contained an 8 protostar forming at the center, and showed the effects of an active accretion disk and an empty bipolar outflow cavity on the SEDs and the images. The stellar mass chosen is high enough that the star has a high luminosity () but low enough that the ionizing luminosity is small; we ignore the latter here. We included gas opacities inside the dust destruction front in the disk since hot gas becomes optically thick. However, the outflow cavities in Paper I were left empty (an optically thin limit), which is not realistic since dust is expected to survive in these outflows; furthermore, under certain conditions of high mass outflow rates, the gas in the outflow may become optically thick near the star. In this paper, we shall study the effects on the SED of the gas and dust in the wind filling the outflow cavities. We shall also consider how the outflow removes angular momentum and accretion energy from the disk. In the next section, we discuss our model setup. We present the SED and image results of our models in Section 3, and summarize our main conclusions in Section 4. In future papers, based on the fiducial model developed here, we shall present the SEDs and images of an evolutionary sequence of MYSOs, and also examine their dependence on the initial conditions of the core.
We now briefly describe the model setup in our previous work, which we continue to use here. For a more detailed description of the model we refer readers to Paper I. Following MT03, a star-forming core is defined as a region of a molecular cloud that forms a single star or a close binary via gravitational collapse. The core is assumed to be spherical, self-gravitating and in near virial equilibrium. In the fiducial case presented here, we study a core with an initial mass of . The core is surrounded by a larger, self-gravitating clump with a mean surface density g cm, which corresponds to pressure at the surface of the core of dyn cm. Such a core has a radius pc ( AU in the fiducial case). We also consider two variants of this model in which the surface density is smaller or larger than the fiducial value by a factor of 0.5 dex. We assume that the initial density distribution in the core is a power-law in spherical radius, , and adopt . Such structures have recently been confirmed to be present in Infrared Dark Clouds (IRDCs) by Butler & Tan (2012). Here we study such a fiducial core at the particular moment when the growing protostar reaches 8 . Such a protostar has a radius and a nuclear luminosity (MT03).
The collapse of the core is described with an inside-out expansion wave solution (Shu 1977, McLaughlin & Pudritz 1997). The expansion wave front moves outward, causing the material inside of it to collapse; the outer part of the core remains static. We use the solution by Ulrich (1976, hereafter the Ulrich solution) to describe the effects of rotation on the velocity field and the density profile inside the sonic point, where the infall becomes supersonic. In this solution, angular momentum is conserved until the material falls onto the midplane and forms an accretion disk. Assuming the rotational energy of the core is 2% of its gravitational energy (), the radius of the disk is AU in the fiducial case (Eq. 13 in Paper I).
Given a sufficiently large effective viscosity, the disk transfers angular momentum outwards and enables the high accretion rates required to form a massive star. The disk around a massive protostar can have a high mass. In all our models with disks, we assume the mass ratio between the disk and the protostar is a constant (), motivated by an expected rise in effective viscosity due to disk self-gravity at about this value of (e.g. Adams et al. 1989, Shu et al. 1990, Gammie 2001). As in the previous paper, the disk structure is described with an “-disk” solution (Shakura & Sunyaev 1973), but we have generalized this to include the effects of the outflow and the accretion infall to the disk (Section 2.1).
The outflow sweeps up part of the envelope and thus controls the efficiency with which the star forms out of the core. We use the same basic setup as in Paper I, but we have reduced the total mass loading rate of the wind relative to the stellar accretion rate to (rather than 0.6 in Paper I), which is a typical value for disk wind models (Königl & Pudritz 2000, also see Section 2.2). For better comparison to Paper I, we keep the opening angle of the outflow cavity to be 51. These assumptions then lead to an increased star formation efficiency of (rather than 0.5 of Paper I). The major difference from Paper I, where we treated the outflow cavity as being completely optically thin, is that we have developed a theoretical wind solution to evaluate a more realistic distribution of gas and dust in the outflow cavity, which will be discussed in Section 2.2.
In Paper I, we considered a standard Shakura-Sunyaev disk with a constant accretion rate and a zero torque inner boundary condition. However, we will now consider several additional effects that may be important in protostellar accretion disks: 1) the accretion rate varies with radius because material in the outflow (or the inflow) leaves (or joins) the disk; 2) similarly, outflows and inflows alter the angular momentum distribution in the disk; 3) the disk itself is massive and contributes to the gravitational potential; and 4) the potential changes with time due to the growth of the disk and the central object. As in Paper I, we to neglect any transfer of angular momentum between the disk and the star via magnetic torques. Furthermore, we continue to ignore the effect of stellar irradiation on the structure of the disk, although we do include its effect on the emission by the disk.
From conservation of the mass in an annulus at , we have
where is the surface density of the disk at and is the radial velocity of the accretion flow (negative for inward motion). Here we have supplemented the standard equation of mass conservation for disks (e.g., Frank et al. 2002) with times , the mass loading rate per unit area of the disk from the (Ulrich) inflow, and times , the mass loss rate per unit area of the disk into the wind. Integration of this relation implies the accretion rate through the disk, , varies with radius :
where is the rate of change of the mass of the disk inside , and are the outflow and inflow rates inside , and is the accretion rate onto the star. The disk extends from the stellar surface at to . We use the Ulrich solution to determine the rate at which mass flows from the envelope onto the disk. For consistency with Paper I, we adopt the same stellar accretion rate, yr, for most of the models we consider.
Similarly, we supplement the standard equation of angular momentum conservation for disks (Frank et al. 2002) with terms representing the gain and loss of angular momentum due to inflow and outflow, respectively,
where is the viscous torque, is the specific angular momentum of the inflow joining the disk, is the specific angular momentum (both rotational and magnetic) of the wind, and is the angular frequency. We have assumed that the disk extends to the surface of the protostar, and we further assume that there is no viscous transport of angular momentum between the disk and protostar (i.e., at the inner edge of the disk). Combining Equations (1) and (3) we have
where is the specific angular momentum of the disk gas. We see that decreases (angular momentum is transferred outwards) (1) due to the viscous torque, provided ; (2) the material falls onto the disk with sub-Keplerian speed; and/or (3) the wind leaves the disk with a super-Keplerian speed due to magnetic torques. Note that the deepening of the gravitational potential well because of the growth of the star and the disk enters into .
As mentioned in the previous section, we assume the disk mass always a constant fraction of the protostellar mass,
where is the rate of change of the total mass in the disk. and are given by the wind solution (see Section 2.2), with the total mass loss rate to the wind also assumed to be a constant fraction of the stellar accretion rate,
Then the total mass loading rate from the inflow to the disk is
with the mass loading rate per unit area () and its specific angular momentum given by the Ulrich solution. With these boundary conditions, from Equations (1) and (3), , , and other profiles of midplane temperature, disk scale height, can be solved with the -disk model. A detailed formulation can be found in Appendix A. Here we assume is a constant throughout the disk that is self-consistently determined by these boundary conditions. In the final fiducial model of the series .
Although we do not explicitly use an equation of global energy conservation in constructing our model, we note that it can be written as
where is Newton’s gravitational constant, the l.h.s of this equation is the total potential energy released in the accretion process from infinity to the stellar surface, is the energy emitted as radiation when the accreting material hits the surface of the protostar, is the total mechanical energy carried away by the wind, is the energy radiated from the disk due to viscous dissipation (note that radiation from star or part of may be absorbed and remitted from the disk; they are not counted in ), and is the rate of change of the energy content (gravitational and kinetic) of the disk, which is actually much smaller (by 2 orders of magnitude) than and in our case. For the disk wind, in the fiducial model (see Section 2.2), which means the disk luminosity is now smaller by a factor of 0.45 compared with that in Paper I, since more than half the disk energy is extracted by the wind. The accretion luminosity as it transitions from the disk to the star ( in the fiducial case) is assumed to be emitted from the stellar surface homogeneously, which indicates an effective surface temperature of the protostar of K, including the internal luminosity () from the protostar.
2.1.1 Thin Disk
The disk profile is calculated numerically as described in Appendix A with the opacities used in our model. Unlike in Paper I, where we then fit the density and scale height with power-law functions of radius for implementation in the radiative-transfer code, here we use the numerical solutions directly. The accretion luminosity is emitted from the disk, obeying the profile , where is the optical depth from the midplane to the surface (see Appendix A). The midplane temperature, , is calculated from the -disk model. We find that it is also consistent with the temperature from the final radiative transfer simulation, since the heating from the star has negligible penetration to the disk midplane.
The boundary between the inflow and the outflow is described by an Ulrich (1976) streamline, which is the trajectory of an infalling particle with an initial angular momentum but with no pressure forces:
Here is the radial distance to the point at which the trajectory crosses the disk plane (); the disk radius is determined by trajectories in the plane (), which gives . Therefore, under the condition of a thin disk, the boundary between the wind launching region and the inflow joining region of the disk is , which is 270 AU in the fiducial case.
Figure 1 compares the radial profiles of the thin disk in the various models. There are jumps in the profiles of density, scale height and temperature due to the drop of the opacity from dusty conditions at 1400 K to gas-only conditions at 1600 K and then a rapid increase of the gas opacity with temperature to K. The material is piled up just before this transition area at K. The red dashed curves are numerically calculated with a standard disk with a constant accretion rate. In Paper I (Model 8) and in Model 9 of this paper (see Sec. 2.5), we used the fitted power laws, shown as the red dotted lines. The dashed blue curves (Model 10) show the effects of the radial variation in the accretion rate due to the inflow and the outflow, the effects of the angular momentum of the inflow (but not of the outflow), the effects of a time-dependent disk, and the effects of a non-Keplerian potential. The first panel shows the accretion rate through the disk (Eq. 2). Inside 100 AU, increases due to the mass loss to the wind (). From 100 AU to 270 AU the negative (see the discussion in Appendix A) becomes dominant, so starts to decrease. Outside 270 AU, drops rapidly due to the inflow (). One effect of this spatially varying accretion rate is that the surface density of the inner disk (inside 200 AU) is less than that of the outer disk. The black curves show the effects of including the wind torque. Panels 2-5 show that this brings the dust/gas transition inward, while the last panel shows that it reduces the viscous dissipation and the corresponding luminosity radiated from the disk by a factor of 0.45.
2.1.2 Effect of Disk Thickness
The finite thickness of the disk has a substantial effect on where the boundary between the outflow and inflow intersects the disk surface, so we develop a simple approximation to allow for this. The vertical distribution of the density in the disk () has a Gaussian profile determined by the temperature at the disk mid plane (see Eq. 16 in Paper I). In the absence of outflows, the height of the disk at was defined as the height of the innermost Ulrich streamline. This sets a volume density at the disk surface. In Paper I, the disk surface at interior radii was set by having this same volume density. In the presence of an outflow, the disk height in the wind-launching region is defined by the point at which the density in the disk equals the density at the base of the wind, which we have somewhat arbitrarily defined as the point where the wind velocity is 1% of the Keplerian velocity (see Appendix B2). For disks of finite thickness, the innermost Ulrich trajectory intersects the disk surface at . Our method of determining the disk height in the wind-launching region becomes inconsistent with the condition that the wind extend to the innermost Ulrich streamline if that streamline would intersect our model disk beyond its outer edge at . Now, we model the disk as being thin, so that the radiative flux in the disk is vertical, and correspondingly the outer boundary of the disk is a vertical wall at . The problematic case–which includes the fiducial case we consider–corresponds to having the outermost wind trajectory exit our model disk through the vertical wall instead of the disk surface. In fact, the disk is likely to be thinner near the outer edge, since radiation can escape radially as well as vertically. Correspondingly, the density at the base of the wind is likely to be larger than the value we estimate in Appendix B via a self-similar model. We adopt the simple approximation for the problematic case that the height of the disk at the outer radius, , is the same as that of the outermost wind streamline at that radius. Since we do not adjust the wind density, this leads to a density discontinuity in the outer disk such that the disk density is somewhat higher than the density at the base of the wind. We assume that the density at the top of the disk remains constant as the radius decreases, until the density at the base of the wind becomes large enough to match it; inside that point, the density at the top of the disk is defined to be the density at the base of the wind. This approximation has a relatively small effect: it needs to be applied only for AU, i.e., a relatively small fraction of the disk, and it leads to a density discontinuity of only 20% between the disk and the wind in this region.
Bipolar outflows are ubiquitous around protostars (e.g. Königl & Pudritz 2000). These outflows (or winds: we will use these two terms interchangeably in this paper) are believed to be driven by the magnetic fields threading the rotating disk (or the star) and powered by the accretion. Current wind models can be categorized into two classes: disk winds, in which the wind flows along the magnetic field lines threading the disk (Königl & Pudritz 2000), and X-winds, in which the wind is launched from a narrow interaction region between the stellar magnetic field and the inner edge of a truncated accretion disk (Shu et al. 2000). Although the properties of these two winds in the near-disk region are quite different, their momentum distributions become similar at regions far from the disk, with (Matzner & McKee 1999). In Paper I, this property, together with an assumed star formation efficiency from the core, was used to determine the opening angle of the outflow cavities, which however were left empty. Here, we now include a detailed density distribution for such outflows.
In this paper we consider only disk winds. We base our analysis on the wind calculated by Blandford & Payne (1982, we hereafter refer to this as the BP wind) and refined by Contopoulos & Lovelace (1994), but we alter this so that the outermost streamline of the wind, at cylindrical radius , coincides with the innermost streamline of the accretion flow. In the Ulrich solution for the accretion flow, this streamline intersects the disk midplane at . We assume that there is a central cavity in the wind due possibly to the poloidal field of the protostar, as in the X-wind model (Shu et al. 1995), or to a line-driven wind from the protostar; the radius of this cavity is . For X-winds, one can show that is qualitatively similar to a BP streamline; since we are considering a BP-like wind here, we assume that does in fact have the shape of a BP streamline. Wind streamlines begin at a cylindrical radius and at a height . We assume that the innermost wind streamline emanates from the surface of the star, . BP winds have a density at the base of the outflow that varies as and a velocity that scales as the Keplerian velocity, . As a result, the mass loss scales logarithmically with the radius of the outermost streamline,
where we have assumed that the disk is thin, so that the mass flux is determined by the -component of the velocity. The density in the outflow at a height above the disk surface is given in Appendix B, in terms of the cylindrical radius of the footpoint, .
As discussed in the last section, the outflow extracts angular momentum and accretion energy from the disk. The specific angular momentum of a self-similar disk wind is (e.g. Königl & Pudritz 2000)
where is the cylindrical radius of the Alfvén surface (where the wind becomes super-Alfvénic). The ratio between the Alfvén radius and the footpoint radius of a streamline, , can be considered as a constant and has typical value of (e.g. Königl & Pudritz 2000). For the BP wind we approximate here (see Appendix B), . For such a wind, the specific energy of the wind, , is given by (Blandford & Payne 1982)
2.3. Dust and Gas Opacities
Dust shapes the SED of the protostar by transforming visible and UV radiation into the infrared. In Paper I, the dust was assumed to be completely destroyed by sublimation at a temperature K, which made the transition between dust and gas opacities relatively sharp. In order to improve the temperature convergence we make a smoother transition between the dust and gas opacities by linearly reducing the fraction of the dust content from 100% (gas-to-dust ratio 100) at 1400 K to 0% (no dust content) at 1600 K, approximating the results of Semenov et al. (2003). We identify 1600 K as a reference sublimation temperature.
We assume that the wind is dusty if it is launched outside the radius at which dust sublimates, neglecting other dust destruction processes such as shocks. Even though the temperature in the the wind launched from the inner dust-free disk soon drops below the dust sublimation temperature, this part of the wind is still likely to remain dust-free since dust grains are not expected to form in such a short time at these densities. We therefore define the boundary between the dust-free and dusty outflow by a streamline starting from the disk surface where the surface temperature is at the sublimation temperature (1600 K). In our fiducial model, this radius is 2 AU. Figure 2 shows the streamlines (white dotted curves) and the density profile (color scale in the upper panels). Each interval of the streamline contains 10% of . The dusty and dust-free parts of the wind can be recognized in the temperature profile. About 66% of the wind is dusty. Note because of the finite disk thickness, the outermost wind streamline emerges from the disk surface at , as discussed in Section 2.1.2. On the other hand, we used the thin-disk equations in calculating the dynamics and vertical profile of the disk, and in that approximation the wind emerges from the midplane at (see Fig. 1).
|, , , pc, , AU|
|8||gas ( K)||standard disk model||empty|
|dust ( K)||with constant accretion rate|
|9||smoother transition||as above||empty|
|(1400 K 1600 K)|
|effects of inflow / massive growing disk on|
|10||as above||accretion rate / angular momentum distribution;||empty|
|effects of outflow on accretion rate|
|11||as above||including effects of outflow in extracting angular||empty|
|momentum and accretion energy from disk|
|12||as above||as above||dust free|
|13||as above||as above||dust + gas|
|Model 13l||Model 13||Model 13h|
|initial core mass ()||60||60||60|
|mean surface density (g/cm)||0.316||1||3.16|
|core radius (pc)||0.10||0.057||0.032|
|outer radius of disk (AU)||801.4||449.4||253.4|
|stellar accretion rate (/yr)|
|stellar radius ()||11.3||12.0||5.93|
|isolated stellar surface temperature (K)|
|isolated stellar + accretion temperature (K)|
|isolated stellar luminosity ()|
|stellar accretion luminosity ( )|
|disk luminosity ( )|
|total luminosity ()|
We divide the dusty part of the simulation into four different regions: (1) the envelope, (2) the low density regions of the disk (), (3) the high density regions of the disk, and (4) the dusty part of the outflow. Opacity models for the first three regions were described in Paper I. The opacity model for the outflow has smaller grains than that in the envelope and has no ice mantles (Kim et al. 1994). These dust opacity models are default options in Whitney et al.’s (2003a) radiation transfer code.
Gas opacities are used in the dust-free, hot inner disk and in the outflow launched from it. The mean gas opacity at 1500 K is about 3 orders of magnitude lower than the opacity of the dust at K, but becomes comparable or even higher when the temperature reaches K, where the gas becomes collisionally ionized (note that we are not including photoionization from the protostellar radiation field). For temperatures higher than K, we adopt gas opacities from the OP Project (e.g. Seaton 2005). For regions with K, the opacity model is from Alexander & Ferguson (1994). However, the opacity given by Alexander & Ferguson (1994) has a minimum at K and increases when drops below that because of molecule and dust formation. Since we expect minimal dust formation to occur in the rapidly expanding outflow, in its dust-free region where K we adopt the opacity at K. More details about these gas opacities are given in Paper I.
We use the latest version of the Monte Carlo radiation transfer code by Whitney et al. (2003b, 2012 in prep.) to perform our calculations. To calculate the equilibrium temperature, we use the algorithm by Lucy (1999), which is implemented in the new version of the code. It sums the pathlengths of all the photons passing through the cell rather than counting only those that are absorbed (e.g., Bjorkman & Wood 2001, used in Paper I), and it thus can quickly reach temperature convergence with fewer photon packets. Starting from a uniform and cold ( 0.1 K) state, the temperature profile becomes stable within iterations with a number of photon packets equal to the number of grid cells (). There are still some oscillations at the transition regions between gas-only and dusty opacities because of the U-shaped opacity curve as a function of temperature, but they occur mostly in small, localized regions close to the disk midplane. The temperature profile on the disk surface and in the outflow is quite stable. We iterate 20 times, average the temperature profiles of the last 10 iterations, and then use this averaged temperature profile as the equilibrium temperature for the final iteration.
Corrections made by adiabatic (expansion) cooling and (compressional) heating in the outflow and in the accretion flow are also considered. We assume that the temperature of the gas varies only slowly, so that it is approximately in thermal equilibrium. The thermal energy equation for cells in these regions is then
where is the Planck mean opacity, is the velocity field, is the pressure of the gas, is the internal energy of the gas, is the total luminosity, is the total number of photon packets, is the volume of the cell, and sums all the optical depths that photon packets have traveled across the cell (Lucy 1999). From left to right, the four terms in the equation are the emission, the adiabatic heating or cooling, advection, and the absorption respectively. Adiabatic cooling becomes important in the outer region () of the outflow, where the temperature drops to K, since the adiabatic cooling rate is proportional to whereas the thermal emission term is proportional to . In these regions adiabatic expansion reduces the temperature by a factor of 2 to 3. Adiabatic heating in the infalling envelope is unimportant (the emission term is larger than the adiabatic heating term). Note that Equation (13) is valid for both the gas and the dust only in regions where they are well coupled, which happens for densities in dense molecular clouds (Doty & Neufeld 1997). This density threshold is at a radius 1000 AU in the outflow. Beyond this radius, the calculated value of the dust temperature should remain valid, although the gas will become colder. The temperature profiles are shown in Figure 2. The low-temperature region in the outflow that extends to a cylindrical radius of 100 AU at height of 600 AU is due to the lower opacity of the dust-free gas. The black region along the axis is the inner cavity due to the stellar magnetic field or the line-driven stellar wind.
A 1000 1500 grid is used to resolve the azimuthally-symmetric space. In space, 750 cells are used to resolve the disk with a finer grid at smaller ( 400 cells are used to cover the region inside 5 AU). In space, an especially fine grid is made to resolve the midplane as well as the inner and the outer boundary of the outflow. About 700 cells are used between above and below the disk midplane. About 100 cells are used within , and 300 cells between and (the opening angle of the outflow cavity is 51). For each run, a high S/N SED and images at different wavelengths and instrument bands (Spitzer IRAC and MIPS, 2MASS J, H, K bands, GTC CanariCam, Herschel PACS and SPIRE, and SOFIA FORCAST) are produced for viewing angles of and . For most of the models shown here, photon packets are used to obtain SEDs and photon packets are used to produce images.
2.5. Model Series
Following the same methodology of Paper I, we develop a series of models to show the effects of the improvements discussed above. We start this series with Model 8, the last model in Paper I, but now calculated with the new version of the code. In Model 9, we include the smoother transition between the gas opacities and dust opacities in the disk. In Model 10, we improve the disk by adding the effects of varying accretion rate in the disk due to the inflow and the outflow, angular momentum of the inflow, disk growth, and a non-Keplerian potential. The effect of the outflow in extracting angular momentum and accretion energy from the disk is included in Model 11. In Model 12, we apply the wind solution to the outflow cavity but assume it to be dust-free (i.e., only gas opacities in the outflow cavity). We also start to consider the flow terms in the energy balance equation from here. Model 13 is the fiducial model, which has all the elements discussed above: a smooth transition between dusty and dust-free gas, an accretion disk structure that includes the effects of time dependence and outflows, and an outflow that is dusty when launched outside the sublimation radius. The density and temperature structure of Model 13 is shown in Figure (2). The properties of each model are summarized in Table 1.
We also construct two variants of Model 13, one with g cm (Model 13l) and another one with g cm (Model 13h). The size of the core scales as , the density as , the pressure as , and the accretion rate as (MT03). The parameters we use here correspond to a change in the clump surface density by a factor of 0.5 dex and the surface pressure by a factor of 1 dex. Table 2 compares the parameters of Model 13l, 13, and 13h.
3. Results and Discussion
Figure 3 shows the SEDs of the models at inclination angles between the line of sight and the rotational axis of and . At an inclination of , the line of sight passes through the dense envelope, which obscures the protostar and the disk. At an inclination of , the line of sight passes through the outflow cavity; because the density in the cavity is lower than in the envelope, and because only part of the outflow is dusty, we can see much more deeply into the cavity than into the envelope. Some important features in the SEDs include the water ice feature at 3 m and the silicate feature at 10 m. As discussed in Paper I, the original output SEDs show some significant emission features such as the H and Paschen lines, but we subtract these features and smooth the SEDs at wavelengths m, since the line profiles may not be exactly correct given the spectral resolution of the radiative transfer simulation and since this line emission is not the focus of our study.
The differences of the SEDs among Models 8-12 are relatively small. The angle-dependent SEDs of Model 8 are the same as in Paper I, showing that the new method of calculating the temperature produces the same results as the method used in the previous paper. Model 9 has a smoother transition between the gas and dust opacities, which is both physically more realistic and achieves better temperature convergence. Model 10 and Model 11 update the disk profile by including the effects of the inflow, the massive disk, and the outflow. Compared with Models 8 and 9, they have less dense but thicker inner disks, so that starlight can heat the inner parts of the disk more easily, producing warmer disk surfaces, while the base of the envelope is more shielded and becomes cooler, thereby making the near- and mid-IR emission higher and the far-IR peaks lower. In Model 11, half of the accretion energy drives the outflow, reducing the emission compared to that in Model 10, particularly between 1 m and 10 m. In Model 12, gas opacities are included in the outflow cavity. However, most of the outflow is very cool so that the gas opacity is low, making the SEDs very similar to those of Model 11.
The dust in the outflow cavity affects the SED more significantly. As shown by the SED of Model 13 at a viewing angle of 30 in Figure 3, the emission at short wavelengths (m) is now suppressed. The reprocessed emission comes out as extra near- and mid-infrared radiation. At larger inclinations, the fluxes at short wavelengths become higher than in the previous models due to scattering by dust grains. Figure 4 compares the SEDs of Models 12 and 13 at an inclination of in more detail. Flux components of different origins are also shown. The existence of dust in the outflow cavity makes the disk warmer, producing higher mid- and far-IR fluxes. The outflow also becomes hotter and emits much more infrared radiation, especially around 10 to 20 m, where the outflow emission becomes dominant. Depending on the sensitivity of the observations, some of the differences among these models may not be detectable, such as from Model 8 - 12. However, the dust in the disk wind (Model 13) can produce factors of several differences to the SED in near- and mid-IR, which can be constrained by the observations. As we show below, dust in the outflow cavity also has significant effects on the infrared images, which are likely to be more useful for observationally testing these components of the theoretical model.
The SEDs of the fiducial model (Model 13) at four inclinations are shown in Figure 5.111We note that in our highly idealized model, in which the axis of the outflow is exactly straight and the BP streamlines persist to large distances, it is possible to see the protostar through the dust-free part of the outflow only for viewing angles that are very small. For Model 13, in which the dusty outflow commences at AU and the core size is AU, the central dust-free cavity has a size of about 600 AU at , corresponding to an opening angle of about . At greater distances from the protostar, this angle becomes even smaller. In fact, however, protostellar outflows are not observed to remain exactly straight to large distances, so it is quite possible that there is no inclination at which the protostar is completely unobscured by outflow dust. However, in the fiducial model considered here, the amount of obscuration by outflow dust is negligible for moderate angles of inclination (). Here we extend the wind to . If the wind extended to , would increase by 1% at an inclination of 10, and 0.1% at 30. Since the density of the wind drops as (see Appendix B), we conclude that the SEDs are not affected by truncating the wind at . As a result of dust in the outflow, the optical emission is suppressed by a factor of 2–3 when the line of sight passes through the outflow cavity at a viewing angle of 30). However, at a viewing angle of 10, the extinction is negligible and the SED closely follows the stellar spectrum at short wavelengths. The SEDs (in ) at these viewing angles are relatively flat from the near-IR to about 60 m, which may be used as an identifier of near face-on massive protostars. The red curves show the SEDs in scattered light. If the line of sight passes through the envelope, all the light at wavelengths shorter than m has been scattered. Direct emission dominates for models with lines of sight passing through the outflow cavity, even at short wavelengths.
As the extinction varies with inclination of the viewing angle, the IR (or bolometric) luminosity inferred from the observed flux varies, which is known as the “flashlight effect” ( Nakano et al. 1995, Yorke & Bodenheimer 1999). Figure 6 shows how the inferred IR luminosities from the near-IR to the far-IR change with viewing angle. In the near-IR, the flux for a face-on view can be higher than that for an edge-on view by 3 orders of magnitude, and for wavelengths m, it can be larger than the value averaged over the sky by a factor . This contrast decreases as the wavelength increases; at 100 m and longer wavelengths, it becomes very small. When the line of sight passes through the outflow cavity close to face-on, the inferred IR luminosities (from the observed ) in the (1 - 8) m band give a good estimate of the total bolometric luminosity, being low by a factor of 2. Figure 7 shows how the bolometric luminosities inferred from the SED depend on the viewing angle and how they compare with the true bolometric luminosities. When looking down the outflow cavity, the bolometric luminosity inferred from the SED is higher than the true value by a factor of 3 The inferred luminosity gradually decreases with increasing inclination angle and drops rapidly when the angle is large enough that the line of sight passes through the envelope. At high inclinations, the inferred bolometric luminosity is lower than its true value by a factor of 6.
Figure 8 compares the SEDs of Model 13, Model 13l and Model 13h, which have cores with the same mass but which are embedded in clumps of different surface densities, . As noted above, the density in the core scales as , which has two effects on the SED: First, since the accretion rate scales inversely with the free-fall time, which varies as , the accretion luminosity, which scales with the accretion rate, varies as ; this affects the SED at all wavelengths. Second, since the surface density of the core scales linearly with that of the clump in which it is embedded, the extinction due to the core is proportional to ; this affects the SED mainly at mid-IR and shorter wavelengths. Note, that we are not including the additional extinction that would result from the clump itself. So, as shown in the SEDs, Model 13h generally has higher fluxes at all wavelengths than Model 13 and Model 13l, due to the much higher luminosity. At lower inclination, the higher extinction due to the outflow makes the FUV fluxes in Models 13 lower than those in Model 13l (although these fluxes will in general be difficult to observe due to additional foreground extinction from dust in the clump and the diffuse Galactic ISM). At higher inclination, the lower extinction of the core makes the mid-IR flux of Model 13l higher than that of Model 13. In this case, the mid-IR flux from the disk in Model 13l is higher and less buried by the flux from the envelope, producing a more prominent 10m silicate absorption feature than in Model 13. This again implies that the depth of the silicate feature depends not only on the optical depth but also on the ratio of the mid-IR emission from the disk to that from the envelope.
An additional factor is the dependence of the protostellar evolution (especially size of the protostar) on the accretion rate. The protostar in Model 13h is much smaller (by a factor of about 0.5) in radius than the other models: it has not yet reached the luminosity wave stage that swells the star (MT03). Thus, its accretion luminosity is significantly higher than that expected from just its enhanced accretion rate. The dependence of these results on the use a simplified one zone protostellar evolution model (MT03) compared to a more detailed treatment, will be examined in a future paper.
Figures 9-12 present the images for the fiducial model at inclinations of and . We show both resolved images and those after convolving with the PSFs of 16 different instrument bands. A distance of 1 kpc is assumed and no foreground extinction is applied. Here, photons are used to produce images, but one still can see the noise due to the Monte Carlo method in the extended outflow, especially in the bands with low fluxes.
These results confirm the conclusion of our previous study: The outflow cavity has a significant influence on the images up to m, especially in the mid-IR, where the cavity wall (the inner part of the accretion flow that is heated by the protostar) dominates the morphology of the image. At around 20 m, the thermal emission from the dusty outflow itself dominates the morphology. At an inclination of , the receding side of the outflow is much dimmer (by 1 to 2 orders of magnitude at 10 - 20 m) than the facing side. The source becomes more symmetric at longer wavelengths. At m, the receding outflow and cavity appear with a flux within a factor of 10 of that from the approaching side. The envelope dominates the image at wavelengths m. In addition to the wide angle outflow, a much narrower jet-like structure due to the denser region in the outflow close to the axis is also seen at an inclination of . There is also a trend that from near-IR to m, the outflow appears to be narrower. Recall that we have assumed that the outflow extends to a radius of . This extended outflow material is brightest in the near-IR because it is dominated by the light from the hot central region that is scattered towards the observer. Note that our model does not include the material of the self-gravitating clump that is pressurizing the core. A real source would most likely be embedded in such a clump, and the appearance of the outflow outside the core would most likely be affected by additional foreground extinction.
At an inclination of , the star and disk could be seen directly in the previous models (up to Model 11) because the outflow cavity was treated as being optically thin. How are the locations of the surfaces changed by the existence of the dust in the outflow cavity? Figure 13 shows the surfaces at different wavelengths seen from outside of the source at an inclination of for Model 13. Even though there is dust in the outflow, the observer may still see down to the very base of the disk wind close to the protostar and disk even in the optical band. The surface basically lies on the surface of the accretion disk in the near infrared bands, and reaches deeper as the wavelength increases. At inclinations large enough that the line of sight passes through the envelope, the source is relatively unobscured at wavelengths m. With the angular resolution achieved by large sub-mm interferometers such as ALMA, such a disk (with a size of 1000 AU) could be readily resolved at a distance of several kpc.
Figure 14 shows the intensity profiles along the outflow axis at inclinations of and . The profiles are convolved with a beam with FHWM of 4 ″ and normalized by the mean intensity at that wavelength of a wide strip along the axis across the core. This strip extends to where the intensity drops to of the maximum intensity. As discussed above, at an inclination of , as one goes from shorter to longer wavelengths, the profile becomes more symmetric and the brightest point moves closer to the protostar. The extended facing outflow () is brightest in the near-IR; the normalized brightness at 70 m is similar to that at 10 m. The opposite side of the outflow can be best seen in the far-IR. At an inclination of , the profiles are much more symmetric. The degree of symmetry in the brightness profile appears to be a promising indicator of the inclination of observed massive protostars.
Figure 15 compares the images of Model 12 (dust-free outflow) and Model 13 (dusty outflow), showing that the dust makes both sides of the outflow cavities brighter at all the wavelengths we considered (which are between 2 and 70 m) except at 18 m, where the emission of the dusty outflow dominates and the flux is more concentrated to the central region. The dust in the disk wind also causes the outflow to appear narrower at and 18 m than in the near-IR. Figure 16 shows similar intensity profiles as Figure 14, but perpendicular to the axis at an inclination of 60. At 10 and 18 m, the profiles of the bright component of the emission are narrower in Model 13. The difference between the two models becomes small for m.
The effects of the surface density of the surrounding clump, , are shown in Figure 17. The surface brightnesses are normalized to the maximum value of Model 13h at each wavelength. Model 13h has an envelope with higher extinction, which causes a higher contrast between the two sides of the outflow than in Model 13 or Model 13l. The sides becomes more symmetric when is lower (Model 13l). Also, we can see that the outflow cavity is brighter at mid-IR wavelengths (8 and 10.3 m) in the two cases with higher extinction, as opposed to the model with the lowest extinction, for which the emission is concentrated at the center.
We have constructed radiation transfer models for individual massive star formation, with a 60 initial core and an 8 protostar forming at the center.
1. Compared with Paper I, our model is improved in two important aspects: the effects of the outflow, inflow and time dependence on the disk surface-density profile are considered, and a detailed distribution of dust and gas in the outflow cavity has been studied in a self-consistent way. While most of the improvements in the disk model (except for the reduction in the accretion energy, and hence the luminosity, due to driving the wind) have only minor influence on the results, the inclusion of dust in the outflow cavity can strongly affect both the SEDs and the images.
2. We have presented a simple analytic model for the density distribution in a disk wind that is based on the Blandford & Payne (1982) solution, generalized to allow for an arbitrary outer boundary.
3. SEDs of the series of models at inclinations of 30 and 60 are presented. Including dust in the outflow cavity significantly influences the SEDs at both inclinations. The dust in the outflow suppresses emission at short wavelengths, while contributing a significant amount of near- and mid-IR emission. The disk also becomes warmer due to the re-radiated light from the outflow material, producing higher mid- and far-IR emission. Close to a face-on view, at inclinations from to , the SED from near-IR to about 60 m is relatively flat, which may be used to identify face-on MYSOs. The 10 m silicate feature appears only if the mid-IR flux from the disk is high enough (e.g. due to lower extinction of the envelope) so that it is not buried in the flux from the envelope, which indicates that the depth of this feature is determined both by the total optical depth at 10m and by the relative brightness of the disk and the envelope.
4. We have presented images of the fiducial model for different observational bands, both resolved and convolved with the resolution beams of the instruments. As discussed in Paper I, the outflow cavity is the most significant feature on the images up to 70 m. The outflow cavity facing the observer may dominate the mid-IR flux and determine the morphology, while the opposite side can be seen only at wavelengths 30 m. Even though the outflow cavity is wide, at 10 to 20 m, the dust distribution in the disk wind can make the outflow appear narrower in the observed images. The distribution of the dust in the outflow affects the intensity profiles both along the axis and perpendicular to the axis. Comparison of such profiles with observation may potentially determine the properties of the outflow. At a lower inclination (more face-on), the disk becomes visible even in the optical bands, although foreground extinction from the natal clump and the larger-scale intervening ISM may limit detection of such emission in practice.
5. We studied the effects of the surface density ( and hence the ambient pressure, of the environment in which the core is embedded. Higher makes the core denser, leading to a higher accretion rate. It affects the SED in two ways: the higher accretion rate leads to higher luminosity, and the denser envelope (and outflow) leads to higher extinction, which can suppress emission at mid-IR and shorter wavelengths. The higher extinction also affects the images, making higher contrast between the two sides of the outflow, while in case of lower extinction, the two sides are more symmetric and the central region is more prominent.
6. We studied how the inferred luminosity depends on the viewing angle (the flashlight effect) The luminosity inferred from the flux from a nearly face-on view can be higher than that from an edge-on view by 2 orders of magnitude in the near-IR and mid-IR. The inferred luminosity from a nearly face-on view is typically higher than the true IR luminosity by a factor of about 2. In the mid- and far-IR up to 70 m, the face-on fluxes, , can also give a good prediction of the bolometric luminosity. When integrated from a SED, the inferred bolometric luminosity is higher than the true value by a factor of 3 if looking down the outflow cavity, i.e. most of the radiation escapes from the polar direction.
Appendix A A. An Improved -Disk Solution
Here we summarize the formulae we use to include the effects of outflow/inflow and of a massive, growing disk in an Shakura-Sunyaev -disk model. Following the approach in standard textbooks (e.g. Frank et al. 2002), we write the conservation of the mass and angular momentum for an annulus at radius of as
where and are the surface density of the disk at and its rate of growth, and are the surface mass loss/loading rates into the wind/onto the disk at , is the radial velocity of the accretion flow, is the viscous torque, and are the specific angular momenta of the wind and the inflow, and is the angular frequency. Integration from gives
where subscript * indicates values on the surface of the star (i.e. the inner boundary of the disk), and
We can write Equation (A4) as
compared with a disk with constant accretion rate through :
The equation set for a Shakura-Sunyaev disk is:
The radiation pressure term is small for the current fiducial model, and thus neglected here.
With the boundary conditions , , , , and , , and given by the solutions of the wind and inflow, at any moment , once we have a profile of , we can calculate and other disk profiles (, , etc.) by combining Eq. (A3), Eq. (A9) and Eqs. (A11). In fact, we set at first, and calculate at two close moments in time ( and ), and use these two profiles to re-estimate , and iterate until the two converge. Note the size of the disk can also vary at these two moments, i.e. the effect of an increasing disk size is also consistently included. We also use the surface density of the inner adjacent annulus to estimate in the r.h.s of Equation (A9) to make the calculation easier.
The effect of the increasing size of the disk on the disk surface density can be estimated as below. The total derivative of the disk mass can be written as
The size of the disk (Eq. 13 in Paper I), so that . Then we can write
where . We can see that because the size of the disk is increasing, the mass of the disk inside may either increase or decrease depending on the ratio of the surface density at the boundary of the disk to the average surface density of the disk. In our fiducial case, , then , i.e. the disk mass inside is decreasing.
Appendix B B. General Formulation of the Density Distribution of the Outflow
We wish to develop an approximate model for a disk wind that will fill the cavity described in Paper I. Protostellar disk winds are believed to be hydromagnetic in origin (Königl & Pudritz 2000), but we shall not explicitly consider the magnetic field here. Our goal is to determine the density distribution in the wind, and we shall make use of the results for the disk wind evaluated by Blandford & Payne (1982, BP wind) and refined by Contopoulos & Lovelace (1994). Our results asymptotically approach those of Ostriker (1997) for self-similar cylindrical winds. We assume that the disk is thin and Keplerian, that the wind is steady and that it is launched from the surface of the disk in a self-similar fashion.
We define the base of the wind as being at the point at which the vertical velocity reaches some fixed, small fraction of the local Keplerian velocity ; we label this velocity . The precise value of the ratio is irrelevant, but the ratio must be sufficiently small that the base of the wind is close to the surface of the disk. We allow for the finite thickness of the disk, but we assume that it is small so that the velocity normal to the surface of the disk is approximately in the direction. Let be the cylindrical radius of a streamline and denote the value of at the base of the wind. The corresponding values of the innermost and outermost radii at the base of the wind are and , respectively. The height of the base of the wind, i.e. the surface of the disk, is denoted by . For a streamline starting from [, ], we define the effective height to be . We assume that there is a central cavity in the wind created by the magnetic field threading the star (Shu et al. 1995); let be the radius of this cavity. We further assume that the wind is confined by the ambient medium, which sets , the radius of the outer boundary of the wind; we shall specifiy this below. We introduce the normalized radius at the base of the wind,
We then have and , where and are measured on the innermost streamline. Let be the density at the base of the wind for the streamline originating at . We assume that disk wind is launched in a self-similar manner, so that varies as a power-law in radius,
where is the density at the base of the wind along the innermost streamline.
The mass-loss rate from a ring of width on both sides of the surface of the disk is
where we have made use of assumed power-law scalings of the density and velocity. The total mass-loss rate from the disk is then
BP winds have , and in this case .
The mass loss rate through rings of width at a height above and below the surface of the disk is approximately the same as that in the plane of the disk,
be the ratio of the vertical velocity to the Keplerian velocity at the footpoint, and let
This is a generalization of the results of Ostriker (1997) for cylindrical winds, which have . Note that if the variation of is small (compared to that of ), then Equation (B10) implies that (Matzner & McKee 1999). In particular, this is the case for X-winds (Shu et al. 1994, 1995), for which .
b.1. B1. Approximate Streamlines
In order to complete our estimate of the density in the wind, we need to specify the streamlines. We describe the streamlines in terms of the dimensionless radius,
and the dimensionless height,
The innermost and outermost streamlines are and , respectively. We assume that the streamlines emanating from points in the disk are a weighted mean of the innermost and outermost streamlines:
Since this relation applies at the surface of the disk (), one finds that
Thus, for the innermost streamline (), Equation (B13) gives , and for the outermost one () it gives , as it should. The ansatz we have adopted thus provides a smooth transition from the innermost to the outermost streamline as the footpoint moves from to . In terms of the dimensionless radius, we have
It is important to note that in both this expression and in Equation (B13), the radius of the innermost and outermost streamlines are evaluated at the same value of the normalized height, , not at the same height above the disk, .