Determining the Parameters of Massive Protostellar Clouds via Radiative Transfer Modeling

Determining the Parameters of Massive Protostellar Clouds via Radiative Transfer Modeling

Ya. N. Pavlyuchenkov Institute of Astronomy, Russian Academy of Sciences, Moscow, Russia    D. S. Wiebe Institute of Astronomy, Russian Academy of Sciences, Moscow, Russia    A. M. Fateeva Institute of Astronomy, Russian Academy of Sciences, Moscow, Russia    T. S. Vasyunina Max-Planck Institute for Astronomy, Heidelberg, Germany

A one-dimensional method for reconstructing the structure of prestellar and protostellar clouds is presented. The method is based on radiative transfer computations and a comparison of theoretical and observed intensity distributions at both millimeter and infrared wavelengths. The radiative transfer of dust emission is modeled for specified parameters of the density distribution, central star, and external background, and the theoretical distribution of the dust temperature inside the cloud is determined. The intensity distributions at millimeter and IR wavelengths are computed and quantitatively compared with observational data. The best-fit model parameters are determined using a genetic minimization algorithm, which makes it possible to reveal the ranges of parameter degeneracy as well. The method is illustrated by modeling the structure of the two infrared dark clouds IRDC-320.27+029 (P2) and IRDC-321.73+005 (P2). The derived density and temperature distributions can be used to model the chemical structure and spectral maps in molecular lines.



Accepted for publication in Astronomy Reports

I Introduction

Massive stars are one of the most important components of our Galaxy and other galaxies, since they determine the energy balance in the interstellar medium and are the principal source of the heavy elements. However, the formation of massive stars is understood less well than the formation of low-mass stars (). It is not yet known whether massive stars form like low-mass solar-type stars or stars with different masses form in different ways (see Zinnecker ()).

The general scenario of solar-type star formation has been developed with a fairly high degree of confidence. It starts with the formation of a dense collapsing core in a molecular cloud. During the core collapse, a central protostellar object and circumstellar torus appear, which are later transformed into a young star and a gas–dust disk, respectively. The durations of the various stages of this process and the details of the transitions between them remain open problems. However, objects at different stages of star formation have long been observed: prestellar (starless) cores, class 0 and I protostars, classical T Tauri stars (class II protostars) and weak-lined T Tauri stars (class III protostars).

No similar sequence of formation stages has been developed for massive stars, as a result of both observational and theoretical difficulties. All regions of massive star formation are far from us, and so require observations with high angular resolution (the distance of the closest, in Orion, is approximately 500 pc; typical distances are of the order of 2 kpc or more). The high gas densities in these regions lead to high extinctions, making optical and even IR observations difficult. As the energy of massive stars influences the surrounding gas starting from the earliest stages of their existence, regions of massive star formation have very complex structure, which requires the application of three-dimensional models in theoretical studies.

The detection of massive prestellar cores — the counterparts of low-mass prestellar cores — could be an important step toward understanding the nature of massive star formation. We expect them to have the simplest structure, while still providing important information concerning the initial conditions for massive star formation. The most promising candidate objects are Infrared Dark Clouds (IRDCs). These clouds are seen in absorption against the Galactic IR background over wavelengths from several m to several tens of m; they were detected as a result of surveys with the ISO ref1 () and MSX ref2 () space telescopes. Their cores, i.e., the densest IRDC regions, can also be detected in emission at millimeter and submillimeter wavelengths. Some of these demonstrate signs of star formation ref3 (), but the crucial test will be the detection of massive starless cores, i.e., without embedded compact sources and with collapse signature in their spectra.

This makes searches for molecules whose spectral transitions could trace the kinematics of massive prestellar cores topical. Finding such molecules and transitions requires modeling the chemical evolution and radiative transfer based on available information about the distributions of the gas and dust densities and temperatures in these objects. In the case of low-mass cores, information about their structure is obtained from data on either their emission in the submillimeter and millimeter launhardt (), or the absorption of the background stellar emission alves () in the optical and NIR.

Spectra of infrared dark clouds can be studied in more detail. First, due to the presence of the infrared background, they can be observed in absorption over a much broader wavelength range than typical clouds in regions of low-mass star formation. Second, like other gas–dust clouds, IRDC cores can be observed both in absorption and emission, making it possible to construct their spectral energy distributions (SEDs) in more detail, and, consequently, to reproduce their physical structure more reliably. However, this requires detailed numerical modeling and is a fairly resource-intensive task.

In this paper, we present a method for studying the density and temperature distributions in prestellar cores, based on modeling the radiative transfer in them taking into account both their own radiation and the absorption of the background radiation at NIR to millimeter wavelengths. As an example, we use two dense IRDC cores from the sample of Vasyunina et al. vasyuninaetal2009 (), IRDC 320.23+0.32 and IRDC 321.71+0.07. The observational data at 1.2 mm described in vasyuninaetal2009 ()111Available at the Strasbourg astronomical Data Center (CDS) at are available for these cores, as well as Spitzer maps at wavelengths from 3.5 to 70 m, obtained as a result of the GLIMPSE glimpse () and MIPSGAL mipsgal () surveys.

We selected the sources P2 in IRDC 320 and P2 in IRDC 321 for our study. Both these sources appear round in millimeter maps (Fig. 1), suggesting that a spherically symmetric approximation may be applicable. There is an important difference between the sources: 70 m emission is associated with IRDS 321, whereas IRDS 320 is not detected in either emission or absorption at 70 m. If we are dealing with massive prestellar or protostellar cores, the presence of emission at 70 m may indicate a later evolutionary stage. For brevity, we further call these cores IRDC 320 and IRDC 321. Table 1 presents their coordinates, distances, and the masses and column densities derived from millimeter observations.


5mm \onelinecaptionsfalse\captionstylenormal \includegraphics[scale=0.35,clip=]./Pictures/fig01a.eps \includegraphics[scale=0.35,clip=]./Pictures/fig01b.eps


[scale=0.3,clip=]./Pictures/fig01c.eps \includegraphics[scale=0.3,clip=]./Pictures/fig01d.eps

Figure 1: Intensity distributions in the IRDC 320 (left column) and IRDC 321 (right column) at 8 m (top row) and 70 m (bottom row). The grey scale shows the IR emission and the contours show the emission at 1.2 mm.

The special feature of the method used here is that we reproduce the observational data at millimeter, near-IR, and mid-IR wavelengths in the framework of a single model. Similar modeling had been performed earlier, but using smaller numbers of wavebands. In particular, observations at 450 and 850 m were used to determine the thermal structure of an IRDC core seen toward the giant molecular cloud W51 ormel (). However, as we will show below, the use of only long wavelength data may be insufficient for this problem.

Parameters of the modeled IRDC cores from vasyuninaetal2009 (). Source (J2000.0) (J2000.0) Distance, Mass (1.2 mm), Column density, kpc cm IRDC 320.23+0.32 (P2) 15 07 56.7 –57 54 27 1.97 50 1.5 IRDC 321.71+0.07 (P2) 15 18 26.7 –57 21 56 2.14 110 3.2

Table 1:

Ii Search for Best-Fit Cloud Parameters

The reconstruction of a source structure based on observational data is an inverse problem. Such problems are usually ill-defined and cannot be solved without some initial assumptions about the structure of the studied region. Usually, a specified model with several free parameters is used, which can be found by fitting the model to the observational data. When the number of free parameters is large, some optimization algorithm must be used to find the best fit. Such algorithms often implement a search for the minimum of a function of several variables, where some criterion for the agreement between the model and observations serves as this function and the relevant free parameters are the variables. Here, we searched for the best-fit parameters applying the PIKAIA numerical code pikaia (), which uses a genetic algorithm. This code yields both the localization of the functional minimum in the space of several parameters and the degeneracy of the parameters. Note that PIKAIA has already been used in some astrophysical studies (see., e.g.pikaia_agb (); pikaia_proplyd ()).

ii.1 Example of Applying the Method: A Two–Component Model

We illustrate the method used to search for the best fit using the example of a simple two-component model of a protostellar cloud; various modifications of this procedure are fairly popular when analyzing observational data 2layer (); robitaille (). We suppose that the cloud consists of two components, each with its own dust temperature and dust surface density . Physically these two components may be, for example, the cloud core and an envelope, but their physical interpretation does not matter from the computational point of view. For definiteness, we consider the first component to be located behind the second. The radiation intensities of the first and second components propagating toward the observer are then given by {align} &I_1(ν) = (1-e^-κ_νΣ_1)B_ν(T_1)
&I_2(ν) = e^-κ_νΣ_2 I_1(ν) +(1-e^-κ_νΣ_2)B_ν(T_2), where is the dust absorption coefficient per unit mass of dust and is the intensity received by the observer. These formulas neglect scattering of the light, and take into account only the intrinsic thermal emission of the medium. Thus, the radiation intensity in the two-component model depends on the four parameters , , and . For convenience, the surface densities and can be converted into molecular-hydrogen column densities and assuming the dust-to-gas mass ratio to be 0.01 bochkarev (). The model will be compared with the SED observed over six wavelengths (1.2 mm, 70, 24, 8, 5.8, and 3.6 m), using the observed intensities toward the assumed center of the cloud, determined as the peak of the 1.2-mm emission. The chosen criterion for agreement between the model and observations takes the form


where and are the observed and theoretical radiation intensities for the th frequency channel and is the dispersion of the logarithm of the observed intensity.


5mm \onelinecaptionsfalse\captionstylenormal \includegraphics[scale=0.62]./Pictures/fig02a.eps \includegraphics[scale=0.55]./Pictures/fig02b.eps \includegraphics[scale=0.55]./Pictures/fig02c.eps

Figure 2: Results of searching for best-fit parameters using the two-component model for the protostellar cloud IRDC 321. The upper plot shows the history of the model’s convergence. The lower left plot demonstrates the localizations of the model parameters in space, and the lower right plot shows the same in space. The black squares show models with and the dots show other models.

Figure 2 shows the results of the minimization based on the PIKAIA algorithm. The computation begins with arbitrary parameter values (only their limits are specified); therefore, there is initially poor agreement with the observations, i.e., the values are very large (Fig. 2, upper plot). After calculating about a thousand models, the algorithm finds a minimum with ; however, the search for the solution continues, since the algorithm checks whether the minimum is local or not. Further, there are no changes, even after computing 10 000 models. The lower plots in Fig. 2 show the locations of the model parameters in the and spaces. The black squares represent the models with , i.e., models that correspond to the minimum in the upper plot of Fig. 2. Figure 3 shows a comparison of the best-fit model with the observational data for IRDC 321. The model parameters are sharply localized and provide a good agreement between the theoretical SED and the observed points.


5mm \onelinecaptionsfalse\captionstylenormal \includegraphics[scale=0.6]./Pictures/fig03.eps

Figure 3: Comparison of the SED of the best-fit two-component model with observations of IRDC 321. The thin curves demonstrate the contributions of the first and second components, and the thick curve the total SED. The circles show the observed intensities.

Formally, the comparison of the model with the observations shows that the real cloud can be represented by a cold core with temperature  Ðš and gas column density  cm, surrounded by a warm envelope with  Ðš and  cm. However, a comparably good agreement with the observed SED is given by a model with the opposite order of the components, i.e., with a hot core and cold envelope, and it is impossible to choose between these two models using this approach.

Another shortcoming of the two-component model is its phenomenological character: we specify layer temperatures and column densities without considering how they could be realized physically, if they can be realized at all. Therefore, the results of such modeling should be interpreted with caution. For instance, the following questions arise in connection with the two-component model. First, if a cold core is surrounded by a warm envelope, how was the envelope warmed? Suppose that the cloud is warmed by some external radiation field. If so, is it possible to choose parameters of this field that can reproduce the observed spectrum of the cloud? Second, if we suppose in a model with a hot core and cold envelope that the inner parts of the cloud are heated by some source, e.g., a protostar, is it possible to find protostar parameters that yield the necessary thermal structure of the cloud and reproduce the observed SED?

We must bear in mind that this modeling took into account only the spectrum toward the center. However, we know both the SED toward the cloud center (or the integrated spectra) and the distribution of the radiation intensity, i.e., maps of the object at the different wavelengths, from the observations. In particular, the NIR intensity distributions for the IRDC cores show that the intensities decrease toward the cloud centers; i.e., these objects are seen in absorption, in contradiction with the initial assumptions of our two-component model. Thus, although such simplified models are formally able to describe the observed spectrum, they may not help us make progress in reconstructing the physical structures of protostellar clouds.

Therefore, it seems logical to develop more physically feasible models of protostellar clouds and compare them with observations in more detail. A model of this type is described in the next subsection.

ii.2 Spatial Model of a Protostellar Cloud

We suppose that the cloud is spherically symmetric, and its density distribution can be described, like those in low-mass prestellar cores Tafalla2002 (), using the expression


where is the central number density of hydrogen nuclei, is the radius of an inner region with approximately constant density , and is the index for the density decrease in the envelope. The outer radius of the cloud is taken to be 1 pc. In the radiative transfer modeling, we must also specify the inner radius of the cloud, which is always taken to be 50 AU.

To take into account possible differences in the evolutionary status, e.g., the existence of a protostar, we assume that there is a source with radius at the center of the cloud, which emits as an ideal black body with temperature . The exact nature of this source is unimportant: this may be the emission of either the protostar or accreting matter. The central source determines the temperature of the inner parts of the cloud. The cloud is irradiated from the outside by isotropic interstellar background emission with color temperature and dilution , which determines the temperature in the outer parts of the cloud.

In addition, we specify the isotropic IR background at the wavelengths used for comparison with the observations. The intensity of this background is equal to the observed intensity at the edge of the cloud. At wavelengths shorter than 10 m, the IR background is probably due to the emission of small interstellar dust grains, polycyclic aromatic hydrocarbons (PAHs) draineli (). Note that the intensity of the NIR background is higher than the intensity of blackbody emission with temperature and dilution , and significantly affects the cloud’s appearance in the observations.

The main agent determining the thermal structure of the cloud is dust, which absorbs, scatters, and reemits the continuum radiation. Therefore, we must solve for the continuum radiative transfer in order to find the dust temperature distribution inside the cloud and compute the distribution of the radiation intensity.

We used the Accelerated Lambda Iteration (ALI) method for the radiative transfer modeling, in which the mean radiation intensity is determined by integrating the radiative transfer equation along randomly chosen directions. The algorithm used is similar to that described in Pavlyuchenkov2004 (), but applies modifications necessary for thermal radiation. The temperature of the medium is found from the equation of radiative equilibrium:


where is the absorption coefficient, is the mean radiation intensity, is the spectral intensity of the radiation, and is the Planck function. In integrating the radiative transfer equation,


scattering is taken into account in the approximation of isotropic coherent scattering, when the source function takes the form


In these equations, is the extinction coefficient and the scattering coefficient. The absorption and scattering coefficients as functions of frequency were computed according to the Mie theory for amorphous silicate grains using a code presented by D. Semenov (Max-Planck Institute for Astronomy, Heidelberg, Germany). We assumed that the size distribution of the dust grains is described by a power law, mrn (), with minimum and maximum grain radii of 0.001 and 10 m.

When the temperature distribution is found, we can compute the theoretical distribution of the radiation intensity. For this, we used the temperature of the medium and the mean radiation intensity, which were determined when modeling the thermal structure of the cloud. To compare the theoretical and observational distributions of the radiation intensity, we convolved the theoretical distributions with the relevant telescope beams for each wavelength.


5mm \onelinecaptionsfalse\captionstylenormal \includegraphics[scale=0.6]./Pictures/fig04a.eps \includegraphics[scale=0.6]./Pictures/fig04b.eps

Figure 4: Results of radiative transfer modeling both with and without scattering for a representative protostellar cloud model. The left plot shows the temperature distribution over the cloud, and the right plot shows the intensity distribution at 24 m.

We will now demonstrate that we must include scattering in the radiative-transfer modeling. Figure 4 shows the model temperatures and intensity distributions at 24 m obtained for a representative cloud with a central density  cm and a central source temperature  Ðš, both with and without scattering. The temperature distribution in the cloud essentially does not depend on the scattering, but scattering significantly affects the distribution of the emergent radiation. Without scattering, the intensity distribution has deep minima, with the minimum depth determined by the column density toward the cloud center. After taking scattering into account, the decrease in the intensity toward the cloud center grows weaker, with the minimum intensity being determined by the ratio of the absorption coefficient to the scattering coefficient.

The method used to compute the temperature and intensity distributions was tested using a number of model problems. In particular, our computational results on the dust temperature agree with the solution obtained with the TRANSPHERE-1D numerical code, which was developed by C. P. Dullemond for modeling radiative transfer in dust envelopes taking into account the thermal emission of the dust222

Iii Modeling Results for Irdc 320 and Irdc 321

The method used to model the structures of protostellar clouds and the intensity distributions was applied to IRDC 320 and IRDC 321. We selected five variable model parameters (Table 2); three of these (, and ) characterize the density distribution, and the other two ( and ) the intrinsic and external radiation fields. The other parameters are fixed, and are also presented in Table 2. The inner radius of the cloud was specified in order to approximately take into account the absence of dust near the internal source (due to sublimation). In principle, the dust sublimation radius depends on the source parameters, but we chose a fixed radius of 50 AU, which in our calculations corresponds to the upper limit for the size of the sublimation zones around hot stars. If there is no star, than the inner radius of the cloud essentially does not affect the distribution of emitted radiation. The outer radius of the cloud, 1 pc, is equal to the observed radii of the studied cores. The radius of the central source (star), 5  whitney (), is chosen fairly arbitrarily, assuming that, within the explored range of stellar parameters, the heating will depend on the star’s luminosity, and an incorrect choice of the star’s radius can be compensated for through a suitable correction of its temperature. The temperature of the interstellar radiation field was taken to be  K.

For each parameter combination, we modeled the radiative transfer, computed the intensity distributions, and quantitatively compared the computed and observed distributions. The search for the best-fit values of the variable parameters was done with the PIKAIA genetic algorithm. The models were compared with observations at 1.2 mm and 70, 24 and 8 m. As the clouds were assumed to be spherically symmetric, the radiation intensity depends only on the angular distance from the cloud center. When constructing the one-dimensional observed intensity distribution, we took the center to coincide with the peak of the 1.2 mm emission.

The agreement between the observed and model intensity distributions at each wavelength can be characterized using the standard criterion. The generalized agreement parameter that must be minimized by the genetic algorithm is equal to the sum of the values for the individual wavelengths. We aim to determine the physical structure of the object using observations at both IR and millimeter wavelengths. The need to use both these wave bands is demonstrated by the ambiguity of the results obtained using only long-wavelength data. The modeling of IRDC 320 using only data at 1.2 mm formally yields the following best-fit parameters of the cloud: the central H density is  cm, the plateau radius AU, the density-profile index 3.2, and the mass 580 . However, this minimum is very flat. Figure 5 exhibits the model distribution in the plane. We can see that the algorithm cannot distinguish between the models with and without an internal source. As an illustration, we selected two models denoted in Fig. 5 by the black squares. These are located at the edges of the range of allowed models: a less dense core with an internal source and a denser core without such a source. Figure 6 displays the distributions of the temperature and intensity for these two models. Even with the very different temperature distributions, the difference between the 1.2 mm intensity profiles is negligible. Thus, it is undesirable to use only millimeter wavelengths when determining the temperatures and densities of prestellar and protostellar cores.


5mm \onelinecaptionsfalse\captionstylenormal \includegraphics[scale=0.62]./Pictures/fig05.eps

Figure 5: Allowed models (), obtained as a result of fitting the theoretical intensity distribution for IRDC 320 to the observational data at 1.2 mm only.

5mm \onelinecaptionsfalse\captionstylenormal \includegraphics[scale=0.62]./Pictures/fig06a.eps \includegraphics[scale=0.62]./Pictures/fig06b.eps

Figure 6: Comparison of the radial temperature (left) and intensity (right) distributions for the two IRDC 320 models denoted by squares in Fig. 5. The dashed curve represents the model with central density  cm and central source temperature 7300 K, and the solid curve the model with central density  cm and without a central source. The spectra were fitted to the observational data at 1.2 mm only.

Further, we present our modeling results obtained using four wavelengths. The localization of the allowed parameters for both sources is shown in Fig. 7, which displays all models with for IRDC 320 and for IRDC 321. The corresponding best-fit parameters for IRDC 320 and IRDC 321, are given in Table 2, together with the masses and hydrogen column densities toward the cloud centers. The temperature and column-density distributions for the best-fit models are shown in Fig. 8.


5mm \onelinecaptionsfalse\captionstylenormal \includegraphics[scale=0.55]./Pictures/fig07a.eps \includegraphics[scale=0.55]./Pictures/fig07b.eps

Figure 7: Localization of the allowed parameters for the IRDC 320 (circles) and IRDC 321 (triangles) models in space (left) and space (right)

5mm \onelinecaptionsfalse\captionstylenormal \includegraphics[scale=0.62]./Pictures/fig08a.eps \includegraphics[scale=0.62]./Pictures/fig08b.eps

Figure 8: Distributions of the density (left) and temperature (right) in the models for IRDC 320 (solid) and IRDC 321 (dashed).

In the case of IRDC 321, the allowed models are located in the region of higher densities, with a concentration of the density toward the cloud center (demonstrated by the high values of ). IRDC 320 is characterized by flatter matter distributions and lower central densities. However, one must recognize that there remains some ambiguity in the solutions, even after fitting the model SEDs at four wavelengths.

The formal best fit for IRDC 321 (when the computation was terminated) has the central density  cm, central plateau radius 5000 AU, , and central source temperature  K. However, there is another solution having nearly the same , with a higher central density of  cm, central plateau radius of 2200 AU, and central source temperature of 6700 K. Thus, there remains some ambiguity for this source, which prevents us from distinguishing between models with denser, compact cores and colder stars and those with less dense, extended cores and hotter stars. Nevertheless, neither of the best-fit models admits the absence of a central source. In other words, the SED features of this source, especially the existence of the 70 m emission, cannot be explained if it is a starless source.

The situation with IRDC 320 is slightly more definite. At the termination of the computations, the best agreement with the observations was achieved for the model with central density  cm, central plateau radius 4000 AU, , and central source temperature  K. Other models with similar values have approximately the same parameters; in contrast to the situation for IRDC 321, there is no alternative group of allowed models. Note that our fitting of the SED of this source over four wavelengths also indicates the existence of a central source.

To test this conclusion, we determined the best-fit parameters for IRDC 320 excluding the heating by the central source, and obtained a worse agreement with the observed SED. In this case, we found that 70 m is the crucial wavelength for determining whether the central source is present or not. The profiles of the source integrated intensity at other wavelengths can be reproduced equally well (or poorly) as in the model with the central source.

Parameters of the IRDC 320 and IRDC 321 models Parameter Notation IRDC 320 IRDC 321 Adjustable parameters Central density of H, cm Plateau radius, AU Index of density profile 3.1 3.8 Star temperature, K 7300 9300 Background dilution Derived parameters Cloud mass, 170 230 Column density, cm Star luminosity, 60 160 Fixed parameters Internal cavity radius, AU 50 Cloud radius, pc 1 Star radius, 5 Background temperature, K

Table 2:

5mm \onelinecaptionsfalse\captionstylenormal \includegraphics[scale=0.6]./Pictures/fig09a.eps \includegraphics[scale=0.6]./Pictures/fig09b.eps \includegraphics[scale=0.6]./Pictures/fig09c.eps \includegraphics[scale=0.6]./Pictures/fig09d.eps \includegraphics[scale=0.6]./Pictures/fig09e.eps \includegraphics[scale=0.6]./Pictures/fig09f.eps \includegraphics[scale=0.6]./Pictures/fig09g.eps \includegraphics[scale=0.6]./Pictures/fig09h.eps

Figure 9: Comparison of the model and observed intensity distributions for IRDC 320 (left) and IRDC 321 (right). The solid curves show the results for the best-fit models, and the dashed curves show the results for the same models without the central sources.

Figure 9 compares the observed and theoretical intensity distributions for the best-fit models. On the whole, a good agreement was achieved for both sources at all wavelengths. The models for IRDC 320 can explain the 1.2 mm emission, the absorption profiles at 8 and 24 m, and the flat intensity distribution at 70 m. The models for IRDC 321 yield emission at both 1.2 mm and 70 m, whereas absorption profiles appear at both 8 and 24 m. The cloud masses are comparable, but the central density and column density is appreciably higher in the IRDC 321 model than the IRDC 320 model. The higher central density, column density, and stellar temperature in the IRDC 321 model compared to those in the IRDC 320 model lead to a higher 1.2-mm intensity and an increase in the emission at 70 m. The dashed curves in Fig. 9 show the intensity distributions calculated for clouds with the same parameters but without the central sources. The presence of a central star alters the intensities at 1.2 mm and 70 m, whereas the 8 m and 24 m intensities are not sensitive to the presence of this source.

Iv Discussion

We have described our method for studying the structure of dense molecular clouds based on modeling their SEDs at both millimeter and IR wavelengths. A number of studies reconstructing the density and estimating the masses of prestellar cores have used observations of dust emission only in the millimeter launhardt (). Since the intensity of this radiation in the optically thin approximation is proportional to the product of the dust column density and temperature, the dust temperature must be known to estimate the mass and density.

Including shorter-wavelength (IR) data makes it possible to obtain a self-consistent reconstruction of the density and temperature distributions. If the studied object is seen in absorption, i.e., the contribution of the dust thermal emission is negligible at these wavelengths, the absorption intensity does not depend on the dust temperature, and is determined only by the dust surface density. On the other hand, if the cloud temperature is high enough to generate the intrinsic IR emission, the intensity of this radiation places strong constraints on the dust temperature. Supplementing these data with millimeter wavelength data significantly narrows the range of allowed model parameters for the cloud.

In other words, the millimeter-wavelength data must be supplemented with shorter-wavelength observations. However, modeling of the IR spectra encounters difficulties. First, the shorter the wavelength, the greater the importance of scattering. This does not seriously affect the temperature distribution, but significantly alters the shape of the SED and the distribution of the integrated intensity (as was noted above). Including the effect of scattering appreciably increases the complexity of the models and makes calculations more resource-intensive. Further, the contribution of PAHs may be important at wavelengths of several m. To take this into account correctly, we must make use of the parameters of both the external UV radiation and the PAH particles, which also increases the model complexity.

In the current study, we included scattering in our modeling, but neglected possible PAH emission. Even in this form, the model makes us able to successfully reproduce a significant part of the observations at the four wavelength ranges used. Among shortcomings of the models, we note the poor fit quality at 24 m: the central intensity minima in the models of both sources are deeper than the observed minima. The higher observed intensity may be related to scattering between the cloud and observer, which is not taken into account in our model. However, efficient scattering at 24 m requires large dust grains, which are unlikely to exist in the interstellar medium. Another possible explanation is PAH emission in the cloud envelope (due to the influence of interstellar UV radiation). However, this explanation is also unlikely: the strongest PAH bands are in the 7-8 m band, whereas our model successfully reproduces the intensity of the 8 m emission. An important factor that could affect the spectrum of the outer parts of the cloud in the NIR is the stochastic heating of small dust grains draineli_st (). Although including this effect should not significantly affect the internal structure of the cloud, it may be useful for determining the UV radiation intensity at the cloud periphery. Taking the foreground into account is also very important. In particular, strictly speaking, our conclusion that massive prestellar cores at 70 m should be seen in absorption is valid for the cores considered only if all the background emission arises behind them. This assumption may be approximately correct for the relatively nearby clouds IRDC 320 and IRDC 321, but the contribution of foreground emission may be a critical component of the model in the more general case, and should be determined together with other parameters.

Some of the computed models successfully reproduce the intensity profile at 24 m, but not the distributions at other wavelengths. These are models with either extended, low-density cores (densities of the order of  cm), or very compact cores with rarefied envelopes. This may suggest a more complex morphology of the studied clouds, in particular, clumpiness and deviations from spherical symmetry. Another reason for the discrepancy between the model and the observations may be the failure of the assumption that the central source is compact. In reality, the cloud interiors may harbor a group of protostars, or may undergo volume heating as a result of gravitational collapse. We must also bear in mind the significant uncertainty in the distances of IRDCs, which affect the reconstruction of their geometric parameters. The method used here can be used to find the dust temperature distribution. However, modeling of the gas-phase chemical processes requires knowledge of the gas temperature rather than the dust temperature. The gas temperature can be directly determined from analysis of molecular line intensities; however, in this case, we encounter an ambiguity related to the complex chemical structures of the studied clouds and non-LTE conditions of the line formation. In particular, the temperatures obtained using ammonia lines are relevant only to those parts of the clouds where the ammonia molecules exist. Modeling of the dust temperature allows us to determine this quantity over the whole cloud. For the purpose of chemical modeling, we can assume that the gas and the dust temperatures are the same. This assumption is valid for the dense regions of molecular cloud cores goldsmith2001 (), but fails at their periphery. Note, however, that the chemical composition at the periphery is essentially determined by photo-processes, whose rates do not depend on the temperature.

In spite of the difficulties noted above, the results obtained demonstrate the effectiveness of our proposed method for reconstructing the physical structures of massive protostellar clumps. Further, we will use this method to model the chemical structure of these clumps and construct spectral maps, making it possible to better establish their evolutionary status. For this purpose, we must take into account the decoupling of the gas and dust temperatures at the periphery of the clumps, where gas heating due to the photoelectric effect becomes important. The low densities in these regions decreases the efficiency of gas cooling due to collisions with dust grains, thereby increasing the importance of gas cooling via molecular-line and atomic-line radiation, when the gas equilibrium temperature is higher than the dust temperature.

The authors are grateful to Th. Henning, H. Linz, H. Beuther, B. M. Shustov, and M. S. Kirsanova for useful discussions, and to the anonymous referee for important comments. This work was supported by the Russian Foundation for Basic Research (project 10-02-00612) and a Grant of the President of the Russian Federation for the State Support of Young Russian PhD Scientists (MK-4713.2009.2). This work has made use of the NASA “Astrophysics Data System Abstract Service”, the data obtained in the GLIMPSE survey, a scientific program based on the SPITZER results and financially supported by NASA, and the NASA/IPAC Infrared Science Archive, which is supported by the Jet Propulsion Laboratory of the California Institute of Technology in the framework of a NASA contract.


  • (1) H. Zinnecker, H. M. Yorke, Ann. Rev. Astron. Astrophys. 45, 481 (2007)
  • (2) M. Perault, A. Omont, G. Simon, et al., Astron. Astrophys. 315, 165 (1996)
  • (3) M. P. Egan, R. F. Shipman, S. D. Price, et al., Astrophys. J. 494, 199 (1998)
  • (4) J. M. Rathborne, R. Simon, J. M. Jackson, Astrophys. J. 662, 1082 (2007)
  • (5) R. Launhardt, D. Nutter, D. Ward-Thompson, et al. Astrophys. J. Suppl. 188, 139 (2010)
  • (6) J.F. Alves, C. J. Lada, E.A. Lada, Nature 409, 159 (2001)
  • (7) T. Vasyunina, H. Linz, Th. Henning, et al., Astron. Astrophys. 499, 149 (2009)
  • (8) R. A. Benjamin, E. Churchwell, B. L. Babler, et al., Publ. Astron. Soc. Pacific 115, 953 (2003)
  • (9) S. J. Carey, A. Noriega-Crespo, D. R. Mizuno, et al., Publ. Astron. Soc. Pacific 121, 76 (2009)
  • (10) C. W. Ormel, R. F. Shipman, V. Ossenkopf, F. P. Helmich, Astron. Astrophys. 439, 613 (2005)
  • (11) P. Charbonneau, Astroph. J. Suppl. 101, 309 (1995)
  • (12) A. Baier, F. Kerschbaum, T. Lebzelter, Astron. Astrophys. 516, 45 (2010)
  • (13) C. Brinch, M. R. Hogerheijde, S. Richling, Astron. Astrophys. 489, 607 (2008)
  • (14) H. Beuther, J. Steinacker, Astrophys. J. 656, 85 (2007)
  • (15) Th. Robitaille, B.A. Whitney, R. Indebetouw, K. Wood, P. Denzmore, Astrophys. J. Suppl. 167, 256 (2006)
  • (16) N. G. Bochkarev, Fundamentals of Physics of Interstellar Media (URSS, Moscow, 2009) [in Russian].
  • (17) M. Tafalla, P. C. Myers, P.Caselli, et al., Astrophys. J. 569, 815 (2002)
  • (18) B. T. Draine, A. Li, Astrophys. J. 657, 810 (2007)
  • (19) Ya. N. Pavlyuchenkov, B. M. Shustov, Astron. Rep. 48, 315 (2004)
  • (20) J. S. Mathis, W. Rumpl, K. H. Nordsieck, Astrophys. J. 217, 425 (1977)
  • (21) B.A. Whitney, R. Indebetouw, J.E. Bjorkman, K.J.S. Wood, Astrophys. J. 617, 1177 (2004)
  • (22) B. T. Draine, A. Li, Astrophys. J. 551, 807 (2001)
  • (23) P. F. Goldsmith, Astrophys. J. 557, 736 (2001)
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description