Stochastic grain heating and mid-infrared emission in protostellar cores

6.8 MB·English

by  Ya. N. Pavlyuchenkov

Mon.Not.R.Astron.Soc.000,000–000 (0000) Printed4January2012 (MNLATEXstylefilev2.2)
Stochastic grain heating and mid-infrared emission in
protostellar cores
Ya. N. Pavlyuchenkov1⋆, D. S. Wiebe1, V. V. Akimkin1, M. S. Khramtsova1,
2
and Th. Henning2
1
0 1Institute of Astronomy, Russian Academy of Sciences, Pyatnitskaya str. 48, Moscow 119017, Russia
2 2Max Planck Institute for Astronomy, K¨onigstuhl 17, Heidelberg D69117, Germany

n
a
J 4January2012

3

] ABSTRACT
A Stochastic heating of small grains is often mentioned as a primary cause of large in-
G frared (IR) fluxes from star-forming galaxies, e.g. at 24µm. If the mechanism does
work at a galaxy-wide scale, it should show up at smaller scales as well. We calcu-
.
h latetemperatureprobabilitydensitydistributionswithinamodelprotostellarcorefor
p four dust components: large silicate and graphite grains, small graphite grains, and
-
o polycyclic aromatic hydrocarbon particles. The corresponding spectral energy distri-
r butions are calculated and compared with observations of a representative infrared
t dark cloud core. We show that stochastic heating, induced by the standard interstel-
s
a larradiationfield,cannotexplainhighmid-IRemissiontowardthe centre ofthe core.
[ In orderto reproduce the observedemissionfromthe coreprojectedcentre,in partic-

ular,at24µm,weneedtoincreasethe ambientradiationfieldbyafactorofabout70.
1
v However, the model with enhanced radiation field predicts even higher intensities at
2 the core periphery, giving it a ring-like appearance, that is not observed. We discuss
4 possibleimplications ofthis findingandalsodiscussaroleofothernon-radiativedust
6 heating processes.
0
. Key words: stars:formation– Stars,infrared:ISM – Sources as a function ofwave-
1
length, radiative transfer – Physical Data and Processes
0
2
1
:
v
i 1 INTRODUCTION bletothelatter,thetemperatureevolutionofasinglegrain
X is stochastic and consists of short temperature spikes upon
ThankstoIRAS,ISO,and,in particular, theSpitzerSpace
r a photon absorption, followed by prolonged “cold” periods,
a Telescope,mid-infraredemissionnowattractssignificantat-
when grain temperature is lower than the equilibrium tem-
tention both at small scale, as an indicator of the pro-
perature it would have if heating were continuous (Duley
tostar formation (e.g., Ragan et al. 2007), and at large
1973). Because of the energy relation, mentioned above,
scale, asanindicatoroftheglobal starformation rate(e.g.,
stochasticheatingdependsonthegrainsizeandisonlyim-
Calzetti et al. 2007, 2010). In both cases making reliable
portant for very small grains (VSG), polycyclic aromatic
conclusions is only possible with a radiation transfer (RT)
hydrocarbons (PAH), or other similar particles. The over-
model and a detailed account of thedust thermal balance.
all result for aparticular grain is excess emission at shorter
Numerous models have been developed to describe
wavelengths(duetotemperaturespikes)andloweremission
emergent spectra for a dusty medium. Most of them treat
at longerwavelengths (duetocold intervals).Whenagrain
dust radiative heating as a continuous process that leads
ensemble with an MRN-like size distribution (Mathis et al.
to a common equilibrium temperature for dust grains
1977)isconsidered,onlyexcessmid-IRemission isseen,as
(Egan et al. 1988; Ivezic& Elitzur 1997; Wolf et al. 1999;
emissionatlongerwavelengthsisdominatedbylargegrains
Dullemond & Dominik 2004; Robitaille et al. 2006). How-
that are not susceptible to stochastic heating.
ever,ithasbeenrecognizedlongagothatgrainsofdifferent
sizes respond differently to a photon absorption, depending
Theimportanceofstochasticheatingiswellrecognized
ontheratioofthephotonenergyandthegrainthermalen-
by the community dealing with the galaxy-wide star for-
ergy (Greenberg1968).Iftheformer is greater orcompara-
mation. It now becomes a standard component of galactic
models to explain the emergent spectral energy distribu-
tion(SED)(e.g., Draine & Li2007;Compi`egne et al.2011;
⋆ E-mail:[email protected] Popescu et al. 2011; Baes et al. 2011). There are also mod-
(cid:13)c 0000RAS
2 Ya. N. Pavlyuchenkov et al.
els of protostellar objects and circumstellar discs, which in- cess 24µm emission, supposedly found in two IRDC cores
cludethestochastic heating(e.g., Manske & Henning1998; by Pavlyuchenkovet al. (2011). They attempted to model
Ferland et al. 1998; Wood et al. 2008). However, in studies SEDs of infrared dark cloud cores IRDC 320.27+029(P2)
of individual protostellar objects this process is often ig- and IRDC 321.73+005(P2) (Vasyuninaet al. 2009) in the
nored, and 24µm emission is rather assumed to be an in- range from 8µm to 1.2mm. Pavlyuchenkovet al. (2011)
dication of the presence of an internal heating source (e.g., foundthatamodelofasphericallysymmetriccore,thatonly
Beuther& Steinacker 2007; Rathborneet al. 2010). This is accountsforequilibriumdustthermalemission,allowstore-
probably justified in cases when a compact source is seen produce absorption of background radiation at 8µm, emis-
on a 24µm image, but the interpretation of diffuse 24µm sion(orthelackofemission)at70µmand1.2mm,butfails
emission can beless straightforward. at 24µm. Absorption minima at this wavelength, predicted
A common approach to the SED modelling of star- by the model, turned out to be much deeper than is ac-
formingregionsisanapplicationofthemodifiedblack-body tuallyobserved.Pavlyuchenkovet al.(2011)suggestedthat
law (Hildebrand 1983). This law is extensively used in the an excess emission at 24µm can be generated by stochasti-
analysisoffar-infrared,submillimetre,andmillimetreobser- callyheatedVSGs,thatarewell-knownpotentialsourcesof
vations.However,itsapplicationislimitedonlytotheemis- emission in the mid-IRrange (Kruegel 2003).
sion ofdustdistributedalongthelineofsight.Accordingly, In this paper we use an extended version of the RT
it only allows to find the dust column density and density- model from Pavlyuchenkovet al. (2011) to simulate emis-
weightedtemperature,whichisnotenoughforthechemical sion of a dense interstellar clump heated from outside by
modelling needed tointerpret lineobservations. Also, infer- the interstellar UV radiation and possibly from inside by
ences, based on this law, can beinaccurate andambiguous. a protostellar object. Emission from stochastically heated
Forexample,observedvariationsinthedustemissivityindex VSGsand PAHsis taken into account. Ourinitial goal was
β can be caused either by real changes in dust properties, toreproduceshallow 24µm shadows in IRDCcores studied
byobservationalnoise,orbythetemperaturegradientalong byPavlyuchenkovet al.(2011).Whilepursuingthisgoal,we
theline of sight (Shettyet al. 2009a,b). foundthatastochasticdustheatingmodelofakind,thatis
Infrared Dark Clouds (IRDC), that are believed to be usedingalactic SEDmodelling, beingapplied toindividual
massivecounterpartsoflow-massprestellarcores,aresimul- clumps, predicts their distinct morphological features. We
taneously seen in emission (millimetre and submillimetre) discuss possible implications of this finding and also con-
andabsorption(near-IRandmid-IR).Thus,theirproperRT sider a role of other non-radiativedust heating processes.
modellinghopefullymakesthederiveddensityandtempera- Thestructureofthepaperisthefollowing.InSection2
turestructurelessambiguousandmorereliable.Densityand we discuss the problem of foreground emission and show
temperaturedistributionsinacoreaffectboththeshapeand how its level of uncertainty affects the derived core param-
the depth of the intensity minimum at shorter wavelengths eters.Theprotostellar coremodelandtheRTmethodwith
andtheintensitymaximumatlongerwavelengths.Thus,to stochasticheatingalgorithmaredescribedinSection3.Also,
reproducetheIRDCshadows and theirmillimetre emission in this section results of protostellar core simulations with
counterparts, one needs a spatially resolved, at least, 1D stochasticallyheateddustgrainsarepresented.InSection4
model of theRTin thecloud. we consider possible ways to solve the excess mid-IR emis-
ItneedstobetakenintoaccountthatextendedIRemis- sion problem. Our conclusions are summarized in the last
sion from an infrared dark core consists of three compo- section.
nents,namely,attenuatedbackgroundemission,propercore
emission, and foreground emission. Proper core emission is
caused by dust grains heated both by an internal source
2 CONTINUOUS DUST HEATING AND THE
and by the external radiation and can be estimated using
FOREGROUND EMISSION
the core RT model. A major problem of infrared studies is
theproblemofforegroundsubtraction.Therearetwopossi- Aswementionedin theintroduction,therearetwowaysto
ble approaches to this problem. The first is to assume that tackle the foreground problem. One of them is to assume
theobservedintensity at thedarkest spot of theconsidered that the optical depth at the darkest spot of the studied
region is entirely caused by foreground emission of any na- region is so large that all the emission from this location is
ture, like zodiacal light, foreground interstellar matter, or causedbyforegroundsources.Theotherapproachistoesti-
instrument noise (Stutzet al. 2009). The second is to take matetheforeground emission asasumofthezodiacal light
thezodiacal light contributionfrom theinterplanetarydust (Kelsall et al.1998)andtheinterstellarmattercontribution
cloud model (Kelsall et al. 1998), then to subtract it from (Butler & Tan2009).Thetwoapproachesmayleadtoquite
the signal, and to assume that all the remaining intensity different inferences.
comes from thesource itself. LetustaketheIRDC321.73+005(P2)core(IRDC321
Thefirst approach issafer asin thiscaseonedefinitely for short) from Pavlyuchenkovet al. (2011) as an example.
knows that derived results represent some limiting values. The background intensity at 24µm in its vicinity is about
On the other hand, its usage implies an assumption that 35MJy/ster,whiletheintensityatthepeakofthemillimetre
the core is optically thick in the considered range and does emission islessthan30MJy/ster.Thezodiacal lightcontri-
not produce any proper emission. These are two additional bution at 24µm, as given by the SPOT software, is about
and maybeunwanted constraints to the model. The second 18MJy/ster. As this object is relatively nearby (∼ 2 kpc,
approach seems tobemore honest,butit bringsupaques- Vasyuninaet al. 2009), the smoothed Galactic foreground
tion of a foreground subtraction. contribution, estimated as in Butler & Tan (2009), is less
The foreground problem can be related to the ex- than 10%, so we neglect it here. Subtracting the zodia-
(cid:13)c 0000RAS,MNRAS000,000–000
Mid-infrared emission in protostellar cores 3
cal light only, we get ∼ 17MJy/ster for the background The absorption and scattering coefficients as functions of
emission (plus any unaccounted foreground emission) and frequency for amorphous silicate grains of a given radius
10MJy/ster for the emission at the millimetre peak. These are computed using the Mie theory. To calculate the total
are intensities that havebeen used by Pavlyuchenkovet al. absorptionandscatteringcoefficientsforadustensemblewe
(2011). Alternatively, if we would assume that all the ex- assumethatthesizedistributionofdustgrainsisdescribed
tra emission at the bottom of theshadow is the foreground by a power law, f(a) ∝ a−3.5 (Mathis et al. 1977), with
emission, we would end up with zero emission for the core minimumandmaximumgrainradiiof0.001and10µm. For
and ∼5−7MJy/ster for thesurrounding region. eachparametercombination ofthecoremodel,wesimulate
To illustrate the significance of the choice, we repeat theradiativetransfer,computetheintensitydistributionsat
here calculations from Pavlyuchenkovet al. (2011) using eachconsideredwavelength,andquantitativelycomparethe
thesetwoapproachestotheforegroundestimation.Priorto computedandobserveddistributionsusingthestandardχ2
giving results we recall the basic equations of the adopted criterion. The theoretical distributions are convolved with
RTmodelwithcontinuousdustheating.Thecoreisassumed therelevanttelescopebeamsforeachwavelength.Thesearch
to be spherically symmetric, with the density distribution forbest-fitparametervaluesisperformedwiththePIKAIA
givenbythesameexpressionasisusedforlow-massprestel- genetic algorithm (Charbonneau 1995). More details about
lar cores (Tafalla et al. 2002): the fitting procedure can be found in Pavlyuchenkovet al.
n(H2)= 1+(nr0/r0)p, (1) (2011R)e.sultsof thefit areshown in Figure 1. Intheleft col-
umn we present data for the case when only DIRBE-based
where n0 is the central density, r0 is the radius of an inner zodiacal light estimate is subtracted from the signal (here-
plateau, p describes the density fall-off in the envelope. To
inafter DZF case). Practically, we subtract zodiacal light
takeintoaccounttheheatingbyanembeddedprotostar(or
from 8µm and 24µm data. Its contribution at 70µm and
agroupofprotostars)weassumethatthereisablack-body
1.2mm is negligible. The right column in Figure 1 corre-
source in the core centre with the temperature T∗ and a sponds to thecase when the foreground level is assumed to
fixed radius of 5R⊙. The inner core radius is set to 50AU be equal to the lowest intensity within the IRDC 321 core
(which is the upper limit for the dust sublimation radius),
(hereinafter complete foreground subtraction, CFS case).
and the outer core radius is assumed to be equal to 1pc.
Obviously,intheCFScasewedonothaveproblemsei-
The core is illuminated by the diffuse ambient interstellar
therwith 24µm emission orwithotherwavelengths.Atthe
field that is characterized by a certain colour temperature
sametime,somecoreparameters,thatarederivedunderthis
T and dilution D .
bg bg assumption, are quitedifferent from the results obtained in
We use the Accelerated Lambda Iteration (ALI)
the DZF case. While temperature profiles and central den-
method for the radiative transfer modelling, which is simi-
sities n0 are almost identical, slopes of the power-law enve-
lar to that described in Pavlyuchenkovet al. (2004) and in
lope p are markedly different (Figure 2). In the CFS model
Hogerheijde & van derTak (2000), but with modifications
p = 2.5, which is typical for prestellar objects. In the DZF
for thermal radiation (see also Hubeny 2003, for general
modelthedensityfall-offismuchsteeper,withp=3.7.Ap-
principles of theALI).The mean radiation intensity
parently,intheDZFcasethealgorithmtriedtoreconcilethe
needtohavemorematerialtoaccountforemissionatlonger
Jν =(4π)−1Z Iν(~n)dΩ (2) wavelengths and the need to have less material to account
4π for low absorption at shorter wavelengths. As a result, the
isdeterminedbyintegratingtheradiativetransferequation core mass is 310M⊙ in the DZF model and 750M⊙ in the
CFSmodel.Onemayarguethatthefactoroftwodifference
(~n∇)Iν =κν(Sν −Iν) (3) in the core mass is not that significant. However, different
alongrepresentativedirections.HereIν isthespecificradia- slopes of thedensityradial profilein theenvelopemay lead
tionintensity,Sν isthesourcefunction,κν =αν+σν isthe to a greater effect when applied to chemical modelling and
extinctioncoefficient,αν istheabsorptioncoefficient,andσν molecular line RTmodelling.
isthescatteringcoefficient.Thetemperatureofthemedium WeshouldnotethatingeneraltheCFSmodeldoeslook
T is found from theequation of radiative equilibrium more attractive. It produces more physical density profile,
is consistent with the 24µm intensity radial distribution.
∞ ∞
Also, analysis of the algorithm convergence shows that the
ανJνdν = ανBν(T)dν, (4)
Z Z solution in the CFS case is well-defined, while in the DZF
0 0 case an alternative solution exists with higher central den-
whereBν isthePlanckfunction. Theadoptedconvergence sity(∼109cm−3)andsmallerplateau(Pavlyuchenkovet al.
criterion is that the relative difference in temperature be- 2011)whichisnearlyasgood(byformalχ2criterion)asthe
tweensubsequentiterationsissmallerthan10−3 atallradii. low density solution presented here.
Although this formalism is multidimensional, in this study So,afterall,toreproduceallthebands,itseemsthatit
the RT equation is solved in 1D, under the assumption of isonlynecessarytoassumethatalltheextraemissioncomes
spherical symmetry. from the foreground, i.e. to rely on the CFS approach. But
In the radiative transfer equation, scattering is taken how justified is this assumption? Let usconsider some gen-
intoaccountintheapproximationofisotropiccoherentscat- eral arguments. The distance to the IRDC 321 core that
tering, so that thesource function Sν takes theform we use as an example is about 2kpc. To account for the
Sν = ανBνκ+νσνJν. (5) 2te4rµstmellaemr ifsosrieognr,ouasnddefsrcarcitbioend aasbohvieg,hwaesw70o%ul.dHnoewedevaenr,iann-
(cid:13)c 0000RAS,MNRAS000,000–000
4 Ya. N. Pavlyuchenkov et al.
30 30
IRDC 321 data IRDC 321 data
8 µm model 8 µm model
20 20
Sr Sr
y/ y/
J J
M M
10 10
0 0
0 20 40 60 80 100 0 20 40 60 80 100
Offset, arcsec Offset, arcsec
30 10
IRDC 321 data IRDC 321 data
9
25 24 µm model 24 µm model
8
7
20
Sr Sr 6
y/ 15 y/ 5
J J
M M 4
10
3
2
5
1
0 0
0 20 40 60 80 100 0 20 40 60 80 100
Offset, arcsec Offset, arcsec
400 400
IRDC 321 data IRDC 321 data
70 µm model 70 µm model
350 350
Sr 300 Sr 300
y/ y/
J J
M 250 M 250
200 200
150 150
0 20 40 60 80 100 0 20 40 60 80 100
Offset, arcsec Offset, arcsec
IRDC 321 data IRDC 321 data
0.6 1.2 mm model 0.6 1.2 mm model
m 0.4 m 0.4
a a
e e
b b
Jy/ 0.2 Jy/ 0.2
0 0
0 20 40 60 80 100 0 20 40 60 80 100
Offset, arcsec Offset, arcsec
Figure1.TheoreticalandobservedintegratedintensitiesintheIRDC321corecomputedwithdifferentassumptionsontheforeground.
In the left column all the foreground emission is caused by the zodiacal light estimated using the DIRBE model (DZF). In the right
column it is assumed that the foreground emission intensity is equal to the intensity at the darkest region of the core (CFS). Dots
representindividualobservationalpixels.
(cid:13)c 0000RAS,MNRAS000,000–000
Mid-infrared emission in protostellar cores 5
107
106 F+E+Be–t
E
105 F B
(cid:150)3m 104
, cH 103 G
n
102
101 DIRBE-based foreground
Completely subtracted foreground
100 Figure3.Schemeexplainingvarioussymbolsenteringequations
(6)and(7).
DIRBE-based foreground
Completely subtracted foreground
K
e, 102 0.6
r
u
at
r
e
mp 0.5
e
st t η
Du st, 0.4
a
r
nt
o
101 c 0.3
n
o
102 103 104 105 si
mis 0.2
R, AU
E
Figure2.Densityandtemperatureradialprofilesforthemodels 0.1
presentedinFigure1.
0.0
approximatemethoddescribedinButler & Tan(2009)gives
0 10 20 30 40 50
theexpectedforegroundcontributionofonlyabout10% for
Total 24 µm background, G
this object. Also, according to the Galaxy structure map,
based on the GLIMPSE survey (Churchwell et al. 2009),
there should not be much intervening material in the di- Figure 4. Emission contrast η at 24µm for cores from samples
rection of theIRDC321 core. byVasyuninaetal.(2009)andButler&Tan(2009)asafunction
of the total IR background intensity. Zodiacal light contribution
Let us assume that all IRDC cores are identical, that
issubtracted.
is, they all have the same optical depth τ and emissivity
E (at a certain wavelength). Also, let us assume that the
observedemissionGoutsideofthecoreprojectionisthesum
thenstaysconstant.Ifthepropercoreemissionwereabsent,
of the ‘true’ background emission B (that comes from the
thatis,E =0,wewouldhaveηindependentofG(ofcourse,
material behind the core) and the foreground emission F,
under the assumption that f is independent of G). On the
thatconstitutesthefractionf ofthetotalemission,F =fG
otherhand,if E is not zero, η would bean increasing func-
(Figure 3).Then, theobserved emission contrast
tionofGaslongasGisnottoobigtomakethesecondterm
G−(Be−τ +E+F) small. This is really observed for G < 30MJy ster−1. The
η= . (6) increaseofη by0.4fromG=10MJyster−1 toG=30MJy
G
ster−1 corresponds to E ≈ 5MJy ster−1. At larger values
After some algebra we obtain
of G thesecond term in eq. (7) is small, and η stays nearly
η=(1−f)(1−e−τ)− E. (7) constant with some scatter caused by scatter in f (that is,
G in distances and Galactic structure features).
InFigure4weshowtheobservationally inferred η forcores Solid linein Figure 4isdrawn usingeq.(7) forf =0.4
fromVasyuninaet al.(2009)andButler & Tan(2009)sam- and E = 5MJy ster−1 (τ24 = ∞). Thus, observed η val-
ples at 24µm (after zodiacal light subtraction). As we do ues indicate that there exists a typical foreground frac-
not intend to make any precise conclusions, values are sim- tion f ∼ 0.4, that is close to the value of 0.54 inferred
ply read off by eye from the MIPSGAL cutouts. The plot by Peretto & Fuller (2009). This is a statistical value, of
apparently shows that η grows for G < 30MJy ster−1 and course,thatisnotdirectlyapplicabletospecificobjects.For
(cid:13)c 0000RAS,MNRAS000,000–000
6 Ya. N. Pavlyuchenkov et al.
example, in Figure 4 a point corresponding to IRDC 321 component, Bν(T) is the Planck function. Distributions
lies above the solid line, which is consistent with f being Pi(T) are supposed to be known from the previous itera-
somewhat smaller than 0.4. An analogous plot for an 8µm tion.
contrast does not show any clear correlation with G, with SimilarlytotheRTmodelwithcontinuousheating,out-
an average η8 value of 0.6 and significant scatter. This also lined in the previous section, the method with stochastic
implies f ∼0.4 (assuming E8=0 and τ8 =∞). dustheatingisalsobasedontheALItechniqueandconsists
So,theobservedIRDC24µmcontrastsarebroadlycon- oftwosteps.Atthefirststep,themeanradiationintensityis
sistent with the omnipresence of the intrinsic mid-IR emis- computedateachcellbyintegratingtheRTequation(3)in
sion in IRDCcores. In other words, the DZF approach can 1D along representative directions. At the second step, the
actuallybevalid,indicating,though,thatwemisssomeim- meanintensityJν isusedtoupdatePi(T)foralldustcom-
portant mechanism(s) responsible for the intrinsic mid-IR ponents. Evaluation of Pi(T) is based on the Monte Carlo
emission of the cores. method.Thetemperatureevolutionofanisolateddustgrain
Itisimportanttorememberthatdifferentassumptions issimulatedandthenconvertedintotheprobabilitydensity
on the foreground lead to different inferences on the core distribution. A single grain temperature evolution is com-
structure,mass,chemical composition, and,byimplication, putedin threesteps:
on its evolutionary status. If we intend to use the informa-
• generationoftimeandfrequencysequencesofabsorbed
tion, deduced from the continuum emission modelling, to
photons;
simulatemolecularlineprofiles,wedoneedamoresubstan-
• calculation of grain temperature spikes due to absorp-
tiatedapproachtotheforegroundproblem.Inthefollowing
tion of photons;
sectionwecheckwhetherthediscrepancyinthe24µmemis-
• calculation of continuous grain cooling due to thermal
sion can bealleviated bytaking intoaccount thestochastic
emission between two consecutive absorption events.
heating of small dust grains that is known to affect a SED
in themid-IRrange. More details about thestochastic heating algorithm can be
found in Appendix.
Inthisstudyweassumethatdustconsistsoffourcom-
3 STOCHASTIC DUST HEATING ponents: large silicate grains, large graphite grains, small
graphite grains, and PAH particles. Their parameters are
To include effects of stochastic dust heating, the RT model listedinTable1.Absorptionandscatteringefficiencyfactors
needs to be upgraded. First, it is no longer possible to as- Qabs and Qsca for silicate and graphite particles are taken
ν ν
sume that all dust grains have the same temperature irre- from Laor & Draine(1993).Opticalpropertiesof PAHsare
spective of their size. Second, a small grain is not in ther- calculatedfollowingDraine & Li(2007),where30PAHfea-
malequilibriumwiththeradiationfield,anditstemperature tures (described by the Drude profile) are taken into ac-
varieswithtime.So,ingeneral,dusttemperatureofasingle count. Heat capacities CV(T) are taken from Draine & Li
particle is a function of grain size and time. (2001). SimulationsofIRDCsareperformed using100 log-
To keeptheproblem tractable, wedrop thecontinuous arithmically spaced radial cells. AsweusetheMonte Carlo
sizedistributionandassumethatdustconsistsoffinitenum- method to compute P(T), the resultant distributions are
ber of components, N, so absorption, scattering and emis- quitenoisy.Soweadoptamuchsofterapproachhere,assum-
sion coefficients are given by ing 5% difference in P(T) between subsequentiterations as
astopping criterion. Each modelrequires about 5 Lambda-
N N N
αν = αiν; σν = σνi; jν = jνi. (8) iterations for convergence. The total computation time is
Xi=1 Xi=1 Xi=1 about one day per model on an average PC.
The code was tested against a number of problems. In
Individual absorption and scattering coefficients are given
the case of large grains (when stochastic heating is not im-
by
portant)resultsareconsistentwithtemperaturescomputed
αiν =niπa2iQaνb,is (9) usingthecodewithcontinuousdustheating.Thestochastic
σνi =niπa2iQsνc,ia, (10) hsuecacteinssgfualllgyorreitphrmoduwcaesStEesDtedansdepPair(aTte)lyd.isItnribpuatritoicnuslashr,owwne
whereni,ai,Qaνb,is,andQsνc,ia arethenumberdensity,radius, in Figures 2 and 6 from Popescu et al. (2011).
absorption and scattering efficiency factors for dust grains
of ith type.
3.1 Stochastic heating by interstellar UV
Since an isolated dust particle is subject to tempera-
radiation
ture fluctuations, we consider a large ensemble of identical
particles and describe their thermal state by temperature Forthecomputationsbelow,densitydistributionandthein-
probability density distribution P(T). The P(T)dT is the nersourceparametersarethosefrom thebest-fitIRDC321
fraction of particles with temperatures lying in the interval model from Pavlyuchenkovet al. (2011) with the DIRBE-
(T,T +dT). The emission coefficient for each dust compo- based zodiacal background subtraction. The diffuse inter-
nent is given by stellar UV field is represented by a Planck spectrum with
∞ Tbg =20000K and dilution of 10−16.
jνi =αiν Pi(T)Bν(T)dT, (11) In Figure 5 (left) we present radial dependence of the
Z normalized probability density distribution P(T) for PAH
0 particles(P(T)forVSGslookssimilar).Thespreadoftem-
where Pi(T) is the probability density distribution for ith peratures in the core envelope is very broad. While at any
(cid:13)c 0000RAS,MNRAS000,000–000
Mid-infrared emission in protostellar cores 7
Table 1.Dustparameters.
Component Radius, Heatcapacitypa- Density, Atomic Mass
cm rameterT∗,K gcm−3 mass fraction
d
Largesilicategrains 3·10−5 175 3.5 28 0.70
Largegraphitegrains 2·10−5 450 1.81 12 0.15
Smallgraphitegrains 3·10−7 450 1.81 12 0.05
PAHs 7·10−8 450 2.24 12 0.10
∗ SeeAppendixforthedefinitionofTd.
giventimemostofPAHparticlesarecold(20K),theirsmall
fraction is heated up to 1000K. Maximum temperatures large sil
104 large gra
near the protostar also exceed 1000K. In the rest of the small gra
core P(T) distribution is relatively narrow, corresponding PAH
total
to theequilibrium temperature. 102
TheSEDsforthemodelcoretowardtheposition,offset
by 1′′ from the core centre (to exclude the direct emission
from the central source), are shown in Figure 5 (right). No τ 100
convolution with a telescope beam is applied in this and
subsequent SEDs. It is hard to isolate the relative contri- 10-2
bution from various dustcomponents on asingle spectrum.
To show their roles, after simulating the thermal structure
of the object with all four components simultaneously, we 10-4
computed the emergent spectra separately for each compo-
nent. More precisely, at the ray-tracing stage we only took 0.1 1 10 100 1000
into account emission, absorption, and scattering by a sin- λ, µm
gledustcomponent.Thesecontributionsareshowninright
Figure 6. Optical depths through the entire core toward the
panel of Figure 5 with coloured lines. With the solid black position,offsetby1′′ fromthecorecentre,forvariousdustcom-
line we show the combined spectrum computed for all the
ponents.
four dust populations.
Physically, the spectrum can be divided into three in-
tervals,accordingtotheopticaldepthshowninFigure6.At
IRspectrumofanIRDCcore.Thecontributioncanbemade
mmandsub-mmwavelengthsthermalemissionfromsilicate
more significant by increasing the ambient UV irradiation.
and graphite grains dominates the spectrum. The optical
InFigure7weshowsameplotsasinFigure5,butcomputed
depth at these frequencies is low (Figure 6), and intensities
with theenhanced ambient UVfield,havingT =26000K
bg
depend significantly on details of the thermal structure in and dilution of 2×10−15. The relative fraction of hot PAH
the interior of the core, in particular, on the presence or
particlesismuchhigherinthiscase,and,asaconsequence,
absence of a protostar. At the IR range (1–24µm) thermal
the IR radiation intensity toward the core projected centre
emission from stochastically heatedPAHsismostly respon- (again with 1′′ offset) also increases, this time matching all
sible for the emergent radiation. The total optical depth at
theobservedvalues.Computedintensitiesforthemillimetre
thesewavelengthsisrelativelyhigh,andthespectrumisnot
emission are higher than observed ones because we do not
sensitive to the presence or absence of the internal heating
convolvethemwith thebeamsizethatisquitelarge atthis
source.Atevenhigherfrequencies(λ<1µm)theemergent
wavelength.
spectrum is formed by large silicate grains via scattering of
However, having solved the problem for the central
background radiation. The optical depth here is very high,
spectrum,wehavesimultaneously created anevenmorese-
andthespectrumis notaffected bytheinterior ofthecore.
vere problem for offset positions. This is illustrated in Fig-
InFigure 5(right)intensities observedtoward thecen- ure 8, where we show distributions of 24µm intensity for
treoftheIRDC321 coreareindicated byfilled circles. Ob- themodelcorewiththestandard(bluelines)andenhanced
viously, we do not solve the mid-IR intensity problem even (red lines) ambient UV fields. With solid lines we show the
by adding the stochastic heating effects. More accurately, ‘true’emergentspectrumthatincludesboththeattenuated
if there were only PAHs in the core, combined emission of background and the core emission. The proper core contri-
stochastically heated particles in the core envelope and in butionisshownseparatelywithdashedlines.Obviously,for
the vicinity of the protostar would provide the needed in- the standard UV field stochastic heating makes negligible
tensity at 24µm. However, in the model with all four dust contributiontotheoverall coreemission (bluedashed line),
components the emission of PAHs from the core centre is and the computed intensity at the projected core centre is
completely absorbed byinterveninglarge grains.Atshorter much smaller than the observed one. As we make the UV
wavelengths even combined (core+envelope) PAH emission irradiation stronger, the intensity in the central part does
is not sufficient to explain observed intensity. grow uptotheobserved value,butat thesametime it gets
Apparently, stochastic heating by the standard diffuse even higher at the core periphery. The radial intensity pro-
UVfielddoesnotmakeanoticeablecontributiontothemid- file has a central flat region of about 5′′ in size and reaches
(cid:13)c 0000RAS,MNRAS000,000–000
8 Ya. N. Pavlyuchenkov et al.
background
104 large sil
large gra
102 small gra
PAH
total
100
Sr
Jy/ 10-2
M
, ν
I 10-4
10-6
10-8
0.1 1 10 100 1000
λ, µm
Figure 5.ResultsofRTsimulationsforthemodelwiththestandardinterstellarfield.Left:shownwithcolourislogarithmofthenor-
malizedtemperatureprobabilitydistributionP∗(T)forPAHparticles(P∗(T)=P(T)/max{P(T)}).Right:spectralenergydistribution
towardthepositionoffsetby1′′fromthecorecentre.Spectraarecomputedseparatelyforlargesilicates(redlines),largegraphitegrains
(thickbluelines),smallgraphitegrains(thinbluelines),andPAHparticles(redlines).Theemergentspectrumwithalldustpopulations
taken into account is shown with solid black line. The background interstellar radiation is shown with black dashed line. Filled circles
indicateobservedvaluesfortheIRDC321core.
background
104 large sil
large gra
102 small gra
PAH
total
100
Sr
Jy/ 10-2
M
, ν
I 10-4
10-6
10-8
0.1 1 10 100 1000
λ, µm
Figure 7.ResultsofRTsimulationsforthemodelofIRDC321withenhanced interstellarradiationfield.
themaximumvalueat35′′.Thisring-likeintensitydistribu- in IRDCs,so we mustadmit that stochastic heating donot
tion is formed by stochastically heated small grains in the seemtobeaviableexplanationfortheallegedexcessmid-IR
envelope as schematically shown in Figure 9. The thickness emission in IRDC cores. But the problem would be solved
of the heated layer is low, so the projection effect plays a if we can find some other dust heating mechanism which
role in the formation of the intensity radial distribution. In would work not only in the core envelope but in its whole
other words, column density of heated VSGs and PAHs to- volume. Among various non-radiative dust heating mecha-
ward the centre of the core is lower than toward the core nisms,thereisonethatisoftenincludedinchemicalmodels
periphery, which results in the intensity difference. The in- of prestellar cores. Specifically, these models include man-
tensitydistributionat70µmoverthecoresurfaceissimilar tle evaporation due to stochastic heating of dust grains by
to the one for 24µm except for the strong emission peak cosmicrayparticles(L´eger et al.1985).Inthefollowingsub-
towardthelocationoftheprotostar.Thispeakappearsdue section,wecheckifsamedustheatingthatisresponsiblefor
to lower optical depth at 70µm. mantle destruction can change appreciably thecore SED.
Suchring-likestructures,apparently,arenotubiquitous
(cid:13)c 0000RAS,MNRAS000,000–000
Mid-infrared emission in protostellar cores 9
60 processes is a challenging problem that is further restricted
24mm EnhancedUV,withBG
EnhancedUV, no BG byuncertainties of involved physicaldata.
NormalUV, withBG Hereweuseasimplifiedapproach.Themostsignificant
50 NormalUV, no BG
contributiontostochastic heatingofdustgrainsisassumed
to come from free electrons formed due to propagation of
40 cosmic rays. The mean energy of free electrons produced
Sr] by cosmic ray ionization is 20-35eV (Glassgold & Langer
y/ 1973).Weassumethatthisenergyentirelygoestothether-
MJ 30 malenergyofadustgraininasinglecollision.Weadoptthe
[n standard ionization rate 10−17s−1 and assume that these
I
electronsdonotinteractwithgasandonlycollidewithdust
20
grains.Sincethelastassumptionisunrealistic(thetotalef-
fective cross-section of molecules is much higher than the
10 total effective cross-section of dust grains), we significantly
overestimate the effect.
InFigure10wepresentP(T)distributionforPAHpar-
0
ticles in the model with CR stochastic heating. To isolate
0 20 40 60 80 100 120
the CR effect, we treat heating by external UV radiation
Offset, arcsec
in the integral non-stochastic way which provides the min-
Figure 8. Intensity distributions at 24µm over the model core imum temperature for grains. While PAHs are sometimes
for cases with and without the background radiation as well as heated up to 1000K, theprobability to find a particle with
with and without UV stochastic heating. Observed intensity for such a temperature is low. In other words, heating events
theIRDC321coreisschematicallyshownwithgrayband (DZF areveryrare,andPAHsspendmostofthetimeat thelow-
case). est temperature. The contribution of grains, stochastically
heated by CR, to the SED is even smaller than in the case
of standard UV (right panel of Figure 10). The computed
24µm intensity is only a few hundredth MJy/ster, while
theobserved intensityis about twentyMJy/ster. Given the
significant overestimation of thecollision ratebetween elec-
trons and grains in our model, we conclude that stochastic
heating of dust grains due to CR cannot solve the problem
of high mid-IRintensity toward IRDCcores.
4 DISCUSSION
TheanalysisofIRDCmid-IRpropertiesindicatesthatthese
clouds may possess intrinsic emission in this range. On a
larger scale, it is typically assumed that emission in the
24µm Spitzer band is produced by stochastically heated
small grains. Our study shows, however, that on a smaller
scale stochastic dust heatingby UVphotonsand CR parti-
clescannotsolvetheproblemofexcess24µmemission,that
Pavlyuchenkovet al. (2011) have inferred for some IRDC
Figure9.Schemeexplainingtheformationofthering-likeinten-
cores.
sityprofileinmid-IRforthecoreheatedbyinterstellarradiation.
For a standard interstellar UV field the contribution
of stochastically heated grains to the core mid-IR emission
isnegligible. Ifwetry toovercomethisproblem byenhanc-
3.2 Stochastic heating by cosmic rays
ingtheambientUVfield,theintensitydistributionbecomes
The interaction of cosmic rays (CR) with interstellar ring-like. The reason behind this is quite simple and qual-
medium is a complex process. The collision of highly en- itatively similar to the “limb brightening” effect considered
ergetic CR particles with interstellar molecules results in byLeung et al. (1989)forlonger wavelengths(140–300µm)
formation of secondary particles (pions, mesons, positrons, whicharenotaffectedbystochasticheatingofsmallgrains.
gamma-rays, etc.) that also interact with interstellar gas. The UV irradiation of the core creates the overheated dust
Both primary and secondary non-thermal particles dissoci- layer on its surface. Geometrically, this layer is narrow due
ate and ionize interstellar molecules as well as collide with to strong UV absorption. On the other hand, it is nearly
dust grains and contribute to their stochastic heating. Ide- transparent in the mid-IR range so that the hot surface of
ally, to calculate stochastic heating of dust grains induced the dense core should be observed as a bright rim. Such
by cosmic rays one should know 1) the flux/energy distri- structures are not seen in any of the cores we have checked
bution of all non-thermal particles, and 2) the fraction of in this study. While 24µm rings are actually observed in
the particle kinetic energy that goes into the grain thermal other sources, most of them are classified as planetary neb-
energy upon the collision. The proper modelling of these ulae and circumstellar shells around massive (proto)stars
(cid:13)c 0000RAS,MNRAS000,000–000
10 Ya. N. Pavlyuchenkov et al.
background
104 large sil
large gra
102 small gra
PAH
total
100
Sr
Jy/ 10-2
M
, ν
I 10-4
10-6
10-8
0.1 1 10 100 1000
λ, µm
Figure 10.SameasinFigure5,butforthemodelwithCRstochasticheating.
(Mizuno et al.2010;Wachter et al.2010).Weconcludethat scatteringisisotropic.Ifourassumptionisrelaxedinfavour
thestochastic heatingofsmall particles, atleast, initstyp- of forward scattering we would have some additional con-
icalimplementation,doesnotseemtobeaviablecandidate tribution from thecore centralhot region. However, thein-
for the source of mid-IR emission in IRDC cores. Other ferred core structure implies that it is opaque at 24µm, so
known potential stochastic heating factor, cosmic rays, is evenwithforwardscatteringthecontributionfromtheinner
alsonoteffective,evenifsomewhat extremeparametersare region should not be important.
adopted.
Scattering of the ambient infrared emission can be
Our inability to fit mid-IR observations may at least significant if we assume some grain evolution in the
in part be related to the deficiencies of the model. For ex- cores. Extended mid-IR emission in the 3.6 and 4.5µm
ample, we assume spherical symmetry for our cores, with bands, that was termed MIR cloudshine and coreshine by
smooth density distributions, resulting in large UV optical Steinackeret al. (2010) and Pagani et al. (2010), observed
depths in the core outer parts. The problem with 24µm in some molecular clouds and cores, was interpreted as a
emission would probably be alleviated if UV radiation is result of scattering of the ambient infrared light by large
able to penetrate deeper into the core due to its irregular grains. In principle, non-zero E due to scattering at 24µm
structure. In this case instead of the narrow hot layer we is also possible. In this study we consider only four rep-
wouldhaveamoreextendedhotenvelope.Also,weconsider resentative grain types and do not vary their parameters
only UV and CR heating, but another bulk emission gen- along the core radius. If we would include very large grains
eration mechanisms are possible. Duley& Williams (2011) (a > 10µm) in the model, scattering at 24µm would be
suggested that surface chemical reactions can be an energy greater, making mid-IR shadows less prominent. However,
sourcethatexcitesinfraredemissionbandseveningrainsof these large grains would need to be present in sufficient
‘classical’ sizes. Some chemically induced energy deposition amountsuchthatscatteringdominatedovertheabsorption
industgrainsisalsoimpliedbythenon-thermalmantledes- by small grains. It is also possible that scattering on large
orption mechanism suggested by Garrod et al. (2007). The porousdustaggregates,likethoseconsideredbyOrmel et al.
presence of some unaccounted energy source indirectly fol- (2011), may contribute intoobserved mid-IRflux.
lows from the low sticking probability (of about 0.3) found
Obviously, by varying dust optical properties and size
byPavlyuchenkovet al.(2006)intheirdetailedstudyofthe
distribution along the core radius, we may achieve a better
CB17 core. They suggested that this low value may be an
fit.However,theseparameterscannotbevariedinauncon-
indicationofsomedesorptionmechanismwhichispresentin
strainedway.Grainsintheinterstellarmediumcangrowby
these objects and is different from thermal desorption and
accretion of refractory elements (Voshchinnikov& Henning
cosmic ray induced desorption. Such a desorption mecha-
2010) or volatile species but the potential for this growth
nismwouldalsomeananenergydepositionintodustgrains.
is limited bytheabundanceof heavyelements. So, thepre-
If the core irregular structure is related to turbulence, ferredwayforgrainstogrowtolargesizesiscoagulation.In
thenturbulentdissipation mayprovidesomeheating. How- otherwords,thenumberoflargegrainscanonlyincreaseat
ever,itseemslikelythatintheabsenceofstrongshocksthis theexpenseofsmallgrainsandPAHparticles.Thetwo-fold
mechanism would rather somewhat raise themean temper- outcome of this process is that scattering becomes more ef-
ature of grains, enhancing the core emission in the far-IR fective,butatthesametimecontributionofstochasticheat-
and submillimetre bands, butnot in thenear- and mid-IR. ingdiminishes.Toestimatethephysicallymotivatedbalance
Yet another factor that is capable to affect the mid-IR between these two processes, one needs to include a grain
core emission is scattering. In this study, we assume that growth in the model. In general, the search for dust pa-
(cid:13)c 0000RAS,MNRAS000,000–000