Astronomy&Astrophysicsmanuscriptno.Juvela (cid:13)c ESO2015 January29,2015 Galactic cold cores V. Dust opacity (cid:63) (cid:63)(cid:63) (cid:63)(cid:63)(cid:63) M.Juvela1,I.Ristorcelli2,3,D.J.Marshall4,J.Montillaud1,5,V.-M.Pelkonen1,6,N.Ysard7,P.McGehee8,R.Paladini8, L.Pagani9,J.Malinen1,A.Rivera-Ingraham2,C.Lefe`vre9,L.V.To´th10,L.A.Montier2,3,J.-P.Bernard2,3,P.Martin11 1 DepartmentofPhysics,P.O.Box64,FI-00014,UniversityofHelsinki,Finland,mika.juvela@helsinki.fi 2 Universite´deToulouse,UPS-OMP,IRAP,F-31028Toulousecedex4,France 3 CNRS,IRAP,9Av.colonelRoche,BP44346,F-31028Toulousecedex4,France 4 LaboratoireAIM,IRFU/ServicedAstrophysique-CEA/DSM-CNRS-UniversitParisDiderot,B¢t.709,CEA-Saclay,F-91191, Gif-sur-YvetteCedex,France 5 Institut UTINAM, CNRS UMR 6213, OSU THETA, Universite´ de Franche-Comte´, 41 bis avenue de l’Observatoire, 25000 Besançon,France 5 6 FinnishCentreforAstronomywithESO(FINCA),UniversityofTurku,Vaïsa¨la¨ntie20,FI-21500Piikkio¨,Finland 1 7 IAS,Universite´Paris-Sud,91405Orsaycedex,France 0 8 IPAC,Caltech,Pasadena,USA 2 9 LERMA,CNRSUMR8112,ObservatoiredeParis,61avenuedel’observatoire75014Paris,France n 10 LorańdEo¨tvo¨sUniversity,DepartmentofAstronomy,Pa´zmańyP.s.1/a,H-1117Budapest,Hungary(OTKAK62304) a 11 CITA,UniversityofToronto,60St.GeorgeSt.,Toronto,ONM5S3H8,Canada J ReceivedSeptember15,1996;acceptedMarch16,1997 8 2 ABSTRACT ] A Context. TheprojectGalacticColdCoreshascarriedoutHerschelphotometricobservationsofinterstellarcloudswherethePlanck satellitesurveyhaslocatedcoldandcompactclumps.Thesourcesrepresentdifferentstagesofcloudevolutionfromstarlessclumps G toprotostellarcoresandarelocatedindifferentGalacticenvironments. . Aims. Weexaminethissampleof116Herschelfieldstoestimatethesubmillimetredustopacityandtosearchforvariationsthat h mightbeattributedtotheevolutionarystageofthesourcesortoenvironmentalfactors,includingthelocationwithintheGalaxy. p Methods. ThesubmillimetredustopacitywasderivedfromHerscheldata,andnear-infraredobservationsofthereddeningofback- - groundstarsareconvertedintonear-infraredopticaldepth.Weinvestigatedthesystematicerrorsaffectingtheseparametersandused o modellingtocorrectfortheexpectedbiases.Theratioof250µmandJbandopacitiesiscorrelatedwiththeGalacticlocationand r t thestarformationactivity.Wesearchedforlocalvariationsintheratioτ(250µm)/τ(J)usingthecorrelationplotsandopacityratio s maps. a [ Results. Wefindamedianratioofτ(250µm)/τ(J) = (1.6±0.2)×10−3,whichismorethanthreetimesthemeanvaluereported for the diffuse medium. Assuming an opacity spectral index β = 1.8 instead of β = 2.0, the value would be lower by ∼ 30%. No 1 significantsystematicvariationisdetectedwithGalactocentricdistanceorwithGalacticheight.Examinationoftheτ(250µm)/τ(J) v mapsrevealssixfieldswithclearindicationsofalocalincreaseofsubmillimetreopacityofuptoτ(250µm)/τ(J)∼4×10−3towards 2 thedensestclumps.Theseareallnearbyfieldswithspatiallyresolvedclumpsofhighcolumndensity. 9 Conclusions. We interpret the increase in the far-infrared opacity as a sign of grain growth in the densest and coldest regions of 0 interstellarclouds. 7 0 Keywords. ISM:clouds–Infrared:ISM–Submillimeter:ISM–dust,extinction–Stars:formation–Stars:protostars . 1 0 1. Introduction detectingandclassifyingofalargenumberofcoldandcompact 5 1 Galactic sources. These probably represent different phases in Theall-skysurveyofthePlancksatellite(Tauberetal.2010)has : the evolution of dense interstellar clouds that leads to the for- v enabled a new approach to studying the earliest stages of star mation of stars. Careful analysis of the Planck data has led to i formation. The sub-millimetre measurements, with high sensi- X a list of more than 10000 objects that form the Cold Clump tivity and an angular resolution down to ∼ 4.5(cid:48), have enabled Catalogue of Planck Objects (C3PO, see Planck Collaboration r a XXIII 2011). At Planck resolution, it is not possible to resolve (cid:63) Planck (http://www.esa.int/Planck) is a project of the European gravitationallyboundcoreseveninthenearestmolecularclouds. Space Agency – ESA – with instruments provided by two scientific The low colour temperature of most of the sources (T < 14K) consortia funded by ESA member states (in particular the lead coun- strongly suggests that the Planck clumps must have high col- tries:FranceandItaly)withcontributionsfromNASA(USA),andtele- umn densities, possibly at scales not resolved by Planck, and scopereflectorsprovidedinacollaborationbetweenESAandascien- they probably contain even dense cores. A significant fraction tificConsortiumledandfundedbyDenmark. (cid:63)(cid:63) Herschel is an ESA space observatory with science instruments oftheclumpsmaybetransientstructuresproducedbyturbulent providedbyEuropean-ledPrincipalInvestigatorconsortiaandwithim- flows,however. portantparticipationfromNASA. (cid:63)(cid:63)(cid:63) Tables 1 and E.1 are only available in electronic form at the Within the Herschel Open Time Key Programme Galactic CDS via anonymous ftp to cdsarc.u-strasbg.fr (130.79.128.5) or via Cold Cores, we have carried out dust continuum emission ob- http://cdsweb.u-strasbg.fr/cgi-bin/qcat?J/A+A/ servations of 116 fields that were selected based on Planck de- 1 M.Juvelaetal.:GalacticcoldcoresV.Dustopacity tections listed in C3PO. The fields, which are typically ∼40(cid:48) shape,structure,materialcomposition,andaccretionofmantles in size, were mapped with Herschel PACS and SPIRE instru- canalsoallcontribute. ments (Pilbratt et al. 2010; Poglitsch et al. 2010; Griffin et al. Moreover, Malinen et al. (2011); Juvela & Ysard (2012); 2010) at wavelengths 100–500µm. The higher angular resolu- Ysard et al. (2012) have investigated the impact of radiative tion of Herschel (Poglitsch et al. 2010; Griffin et al. 2010) en- transfer on the results derived from observations under the as- ables studying the internal structure of the Planck clumps, de- sumption of a single average colour temperature. They showed tecting individual cores, and, in conjunction with mid-infrared thatthemixingofdifferenttemperaturesalongthelineofsight data, studying protostellar sources (Montillaud 2014 submit- produces a tendency that is opposite to the observed one. They ted). The inclusion of far-infrared wavelengths helps to deter- concluded that in the absence of internal heating sources, the mine the physical characteristics of the regions and, in particu- observedemissivityincreasetowarddensecloudscannotbeex- lar, to study the properties of dust emission. First results have plainedbyradiativetransfereffects.Itmustoriginateinintrinsic been presented in Planck Collaboration XXIII (2011); Planck variationsoftheopticalpropertiesofthegrains. Collaboration XXII (2011) and in Juvela et al. (2010, 2011, Itis,however,importanttonotethatthedustemissivityin- 2012)(papersI,II,andIII,respectively).Montillaud(2014sub- creaseisnotsystematicallyobservedintheinterstellarmedium mitted) presented the analysis of all clumps and cores found in (ISM,seeNutteretal.(2008);Juvelaetal.(2009);Paradisetal. ourHerschelfields,includingacomparisonwiththepopulation (2009)) . These intriguing results call for a broader investiga- ofyoungstellarobjects(YSOs).Furtherstudiesconcentratedes- tion,makinguseofthelargeobservationsstatisticsprovidedby pecially on high latitude clouds (Malinen et al. (2014), Rivera- Herschel and Planck, probing different Galactic environments. Ingrahametal.andRistorcellietal.,inpreparation). The key questions are still open today: when, where, and how In this paper we concentrate on dust properties and es- dustevolvesbetweendiffuseanddenseregions,whatthephys- pecially on the submillimetre dust opacity. Variations of dust ical conditions enhancing (or preventing) the efficiency of the emission properties have been investigated with far-infrared coagulation process are, what the time scales are, and whether (FIR) and submillimetre observations of diffuse and molec- the process is directly related to specific stages in the cloud or ular clouds, using data from IRAS, COBE, ISOPHOT, the core evolution. Understanding these questions is critical since PRONAOS balloon-borne experiment, and ground-based tele- knowing the dust opacity has a direct impact on many key pa- scopes.Itwasshownthatthelowtemperaturesfoundinasample rametersderivedfromdustemission,suchasthecolumndensi- of molecular clouds (Laureijs et al. 1991; Abergel et al. 1994, ties,masses,andvolumedensitiesoftheclouds.Forthisreason, 1996) and the translucent Polaris Flare cirrus cloud (Bernard itisalsonecessarytoinvestigatethepossiblesystematiceffects et al. 1999) cannot be explained by the extinction of the radi- on the emission and extinction measurements that could cause ation field. An increase of the dust emissivity by a factor of 3 errorsintheopacityestimates. comparedtothestandarddiffusevaluewasneededtoreproduce Thestructureofthepaperisasfollows:Theobservationsare the cold temperatures observed in the Taurus filament L1506 described in Sect. 2. The main results are presented in Sect. 3, (Stepnik et al. 2003). In dense regions, several studies have includingtheestimatesofsubmillimetreandnear-infrared(NIR) shownanopacityincreasebyafactorbetween2to4(Cambre´sy optical depths, the correlations between these variables, and et al. 2001; Kramer et al. 2003; Bianchi et al. 2003; del Burgo thecorrelationswithenvironmentalfactors.Theresultsaredis- et al. 2003; Kiss et al. 2006; Ridderstad et al. 2006; Lehtinen cussedinSect.4beforewelistthefinalconclusionsinSect.5. et al. 2007). More recently, similar results have been obtained with Herschel and Planck in molecular clouds and cold cores 2. Observationsandbasicdataanalysis (Juvela et al. 2011; Planck Collaboration XXV 2011; Martin et al. 2012; Roy et al. 2013). In their detailed modelling of 2.1. Targetselection HerschelobservationsoftheL1506filament,Ysardetal.(2013) characterisedthedustevolutiontowardthedensepartofthefil- TheselectionoftheHerschelfieldsisdescribedinJuvelaetal. ament. The dust emissivity in the outer layers of the filament (2012) and an overview of all the maps is given in Montillaud was found to be consistent with standard grains from the dif- (2014submitted).Weonlyrepeatthemainpointshere. fusemedium,whereastheemissivityincreasesbyafactorof∼2 Planck sub-millimetre observations, together with IRAS above gas densities of a few times 103cm−3. This change has 100µm data, enabled the detection of over 10000 compact beenattributedtotheformationoffluffyaggregatesinthedense sourcesinwhichthedustissignificantlycolderthaninthesur- medium,resultingfromgraincoagulation,assuggestedbypre- rounding regions (Planck Collaboration XXIII 2011). The de- viousstudies(Cambre´syetal.2001;Stepniketal.2003;Bernard tectionprocedureisbasedonthiscolourtemperaturedifference et al. 1999; del Burgo et al. 2003; Kiss et al. 2006; Ridderstad and,furthermore,limitsthesizeofthedetectedclumpstovalues et al. 2006). The average size of aggregates required to fit the below∼12(cid:48)(Montieretal.2010).Thesourcesarebelievedtobe FIR,submillimetre,andextinctionprofilesintheL1506filament Galacticcoldclumpsor,atlargerdistance,entireclouds(Planck is about 0.4 µm (Ysard et al. 2013). This value is close to the CollaborationXXIII2011). smallestgrainsizeneededtoscatterlightefficientlyinthemid- The fields for Herschel follow-up observations were se- IR,whichproducesthe’coreshine’observedtowardanumberof lected using a binning of Planck cold clumps with respect to densecores,whichhasalsobeeninterpretedasaresultofgrain the Galactic longitude and latitude, the estimated dust colour growth(Paganietal.2010;Steinackeretal.2010). temperature, and the clump mass. At the time of source selec- Onthe theoreticalside, variousdust opticalproperty calcu- tion, distance estimates existed for approximately one third of lationshavepredictedasignificantincreaseintheemissivityof the sources in C3PO and, therefore, some sources of unknown aggregates at long wavelengths compared to compact spherical mass were also included. The binning ensured full coverage grains (Wright 1987; Bazell & Dwek 1990; Ossenkopf 1993; of the clump parameter space, especially of the high Galactic ◦ Ossenkopf & Henning 1994; Stognienko et al. 1995; Ko¨hler latitudes and of the outer Galaxy. Galactic latitudes |b| < 1 etal.2011,2012).Thisvariationisshowntobemainlyduetothe wereexcludedbecausethatareaiscoveredbytheHi-GALpro- increase in the porosity fraction with aggregate growth, but the gramme (Molinari et al. 2010). Similarly, regions included in 2 M.Juvelaetal.:GalacticcoldcoresV.Dustopacity other Herschel key programmes such as the Gould Belt sur- polated Planck and IRIS data. The zero points were calculated vey (Andre´ et al. 2010) and HOBYS (Motte et al. 2010) were iteratively, including colour corrections calculated using colour avoided. temperaturesthatwereestimatedfromSPIREdatawithafixed A total of 116 separate fields were observed. The SPIRE spectralindexofβ=2.0. maps are on average over 40(cid:48) in linear size, with an average area of ∼ 1800(arcmin)2. The PACS maps are smaller, with 2.2.3. Calculatingsubmillimetreopticaldepth an average area of ∼ 660(arcmin)2. Most fields contain more thanonePlanckclump,themapsaltogethercovering∼350indi- TheHerschelmapswereconvertedintoestimatesofdustoptical vidualPlanckdetections.Therangeofprobedcolumndensities depthat250µm.Thesurfacebrightnessmapswereconvolvedto extends from diffuse fields with N(H2) ∼ 1021cm−2 to cores acommonresolutionof40(cid:48)(cid:48),andcolourtemperatureswerecal- withN(H2)>1023cm−2.ThefieldsarelistedinTable1andthe culated by fitting the spectral energy distributions (SEDs) with HerschelobservationnumbersareincludedinTableE.1. modified blackbody curves with a constant opacity spectral in- dexofβ=2.0.The250µmopticaldepthwasobtainedfrom 2.2. Herscheldata I (250µm) τ(250µm)= ν , (1) 2.2.1. Herscheldatareduction Bν(T) The fields were mapped with the SPIRE instrument at wave- usingthefitted250µmintensityIν(250µm)andthecolourtem- lengths 250, 350, and 500µm and with the PACS instrument perature T. The calculations were made with 250–500µm data at wavelengths 100 and 160µm. One field was observed with and160–500µmdata.Thefitswereweightedaccordingto15% SPIRE alone (G206.33-25.94, part of the Witch Head Nebula, and7%errorestimatesforPACSandSPIREsurfacebrightness IC 2118). The Herschel observations are discussed in detail measurements,respectively(seeSect.2.2.1). in (Juvela et al. 2012) and (Montillaud 2014 submitted). The The assumed value of β=2.0 may be appropriate for dense SPIRE observations at 250µm, 350µm, and 500µm were re- clumps, although at lower column densities the average value duced with the Herschel Interactive Processing Environment is lower, β ∼ 1.8 (Boulanger et al. 1996; Planck Collaboration HIPE v.10.0.0, using the official pipeline with the iterative XXV 2011), and the value of β may further depend on destriper and the extended emission calibration options. The the Galactic location and the wavelength range (e.g. Planck maps were produced with the naive map-making routine. The Collaboration Int. XIV 2014). If the true value of β were 1.8 PACS data at 100µm and 160µm were processed with HIPE insteadof2.0,ourcolourtemperatureestimateswouldbehigher v. 10.0.0 up to level 1, and the maps were then produced with by∼1Kandtheτ(250µm)valueslowerby∼30%.Furthermore, Scanamorphos v20 (Roussel 2013). In the order of increasing ifthevaluesofβwerecorrelatedwithcolumndensity,theslope wavelength, the resolution of the maps is 7(cid:48)(cid:48), 12(cid:48)(cid:48), 18(cid:48)(cid:48), 25(cid:48)(cid:48), of τ(250µm) vs. τ would be similarly affected. We return to J and 36(cid:48)(cid:48) for the five bands. The raw and pipeline-reduced data theseeffectsinSects.3.7and4. areavailableviatheHerschelScienceArchive,theuser-reduced To estimate the statistical uncertainty of τ(250µm) values, mapsareavailableviaESAsite1. we used Markov chain Monte Carlo (MCMC) runs. The prior The accuracy of the absolute calibration of the SPIRE ob- distribution of temperature values is flat but limited between servations is expected to be better than 7%2. For PACS we as- 5.0K and 35K. In addition to the relative errors quoted above, sume a calibration uncertainty of 15%. This is a conservative we included the uncertainty of the intensity zero points. These estimate that is compatible with the differences of PACS and aretypicallymuchsmallerthantheassumedrelativeerrors,but SpitzerMIPSmeasurementsofextendedemission3. maybeimportantatlowcolumndensities,especiallyat160µm. Thezero-pointerrorsaresystematicbutareincludedsimplyas another component of statistical noise. Their effect is thus re- 2.2.2. Estimatingintensityzeropoints flectedintheerrorestimatesofindividualpixels.Theerrordis- tributionofτ(250µm)isnearlyGaussian,andweusedthestan- Todeterminethezeropointoftheintensityscale,wecompared darddeviationoftheMCMCτ(250µm)samplesastheerroresti- the Herschel maps with Planck data complemented with the mates.Theseestimateswerecalculatedseparatelyforeachpixel IRIS version of the IRAS 100µm data (Miville-Descheˆnes & oftheτ(250µm)maps. Lagache 2005). The Planck and IRIS measurements were in- Because of line-of-sight temperature variations, the derived terpolatedtoHerschelwavelengthsusingfittedmodifiedblack- body curves, B (T )νβ, with a fixed value of the spectral in- τ(250µm) estimates probably systematically underestimate the ν dust dex,β = 2.0.ThelinearcorrelationsbetweenHerschelandthe truevalues(Shettyetal.2009;Malinenetal.2011).Wecannot reference data were extrapolated to zero Planck (+IRIS) sur- directlydeterminethemagnitudeoftheseerrorsbut,withsome face brightness to determine offsets for the Herschel maps. For assumptions,wecanuseradiativetransfermodellingtoestimate SPIRE the uncertainties of these fits are typically ∼1MJysr−1 the magnitude of the bias. The simulations, described in detail inAppendixC,wereusedtoderivebiasmapsthataretakeninto at 250µm and at longer wavelengths smaller in absolute value. accountwhenthedatawerecorrelatedwithτ values. ThederivedintensityzeropointsareindependentofPlanckcal- J ibrationandofanymultiplicativeerrorsinthecomparison.For PACS the correlations are often less well defined, and the zero 2.3. Near-infrareddata points were set directly based on the comparison of the average values of the Herschel maps and the corresponding inter- We used the Two Micron All Sky Survey (2MASS, Skrutskie et al. 2006) to derive estimates of dust column density that 1 http://herschel.esac.esa.int/UserReducedData.shtml are independent of dust emission. We used the method NICER 2 SPIREObserver’smanual, (Lombardi & Alves 2001) and the standard extinction curve http://herschel.esac.esa.int/Documentation.shtml (Cardellietal.1989)toconvertthereddeningofthebackground 3 http://herschel.esac.esa.int/twiki/bin/view/Public/PacsCalibrationWebstars to estimates of J-band optical depth, τJ. Because the cal- 3 M.Juvelaetal.:GalacticcoldcoresV.Dustopacity culations involve only near-infrared bands, the results are ex- Forτ theerrorestimateswereprovidedbytheNICERrou- J pected to be insensitive to the value of the ratio of total to se- tine. For τ(250µm) these were obtained from MCMC calcula- lective extinction, R (Cardelli et al. 1989). The shape of the tions (see Sect. 2.2.3). The comparison between the different V NIR extinction curve is believed to be relatively stable, even cases (for example, regarding the use of 160µm data, back- at high extinctions (e.g. Draine 2003a; Indebetouw et al. 2005; ground subtraction, or gradient corrections) provides informa- Lombardietal.2006;Romań-Zuń˜igaetal.2007;Ascensoetal. tion on the uncertainty caused by some sources of systematic 2013; Wang & Jiang 2014). Some variations are observed with errors. Galacticlocationand/ordensity,butgenerallyonlyatalevelof The τ(250µm) vs. τ(J) points of individual fields and sam- 5%oftheNIRpower-lawindex(e.g.Stead&Hoare2009;Fritz ples of fields were fitted with a linear model to derive the ratio etal.2011).ThisquestionisdiscussedinmoredetailinSect.4.2. τ(250µm)/τ(J). These total least-squares fits take into account Theτ(J)valuesarederivedusingboththeJ-HandH-Kcolours theuncertaintiesinbothvariables.Thefitsweremadeusingei- but,withtheextinctioncurveused,wehavethecorrespondence ther all data points or only data below or above a given τ(J) of E = 0.65τ(J). Flags in the 2MASS catalogue were used limit.Toreducethebiascausedbythesecuts,thedataweredi- J−K to avoid galaxies (ext key not null or gal contam not zero) and vided with the help of a preliminary linear fit to all data points sourceswithuncertainphotometry(ph qualworsethanC). (see Sect. 3.1 for details). The limiting value of τ(J) thus cor- FiveofourfieldsarefullycoveredintheVISTAHemisphere responds to a position on this line, and the cut was performed Survey, VHS (McMahon et al. 2013), which has more than ten using a line that is perpendicular in a coordinate system where times the sensitivity of 2MASS (in H band VHS has a 5-σ de- theaverageuncertaintiesofthetwovariablesareequal. tection threshold of 19.0, compared to a 2MASS point source The τ(250µm)/τ(J) ratios were also calculated for alterna- catalogue completeness limit of ∼16 mag). One of these fields tiveversionsoftheτ(250µm)data,usinglocalbackgroundsub- istoodistanttoobtainareliableextinctionmap,butthedatafor tractionorusingancillarydatainanattempttocorrectforpossi- thefourotherfieldswereanalysedandtheresultscomparedwith blelarge-scaleerrorsinthesurfacebrightnessdata.Thesealter- those obtained with 2MASS data. The fields are G4.18+35.79 nativedataarediscussedinAppendixA. (LDN134),G21.26+12.11,G24.40+4.68,andG358.96+36.75. For the fields G21.26+12.11 and G24.40+4.68, only J- and Ks-band data exist. The data are available in VISTA Survey 3. Results Archive4 , and the VISTA Data Flow System pipeline process- 3.1. Apparentτ(250µm)/τ(J)values ingandsciencearchivearedescribedinIrwinetal.(2004)and Hamblyetal.(2008). Wecalculatedτ(J)andτ(250µm)mapsofthe116fieldsasde- Extinction maps are produced by averaging extinction es- scribedinSects.2.2and2.3.Thecorrelationsbetweenτ(J)and timates of individual stars with a Gaussian weighting func- τ(250µm) were calculated at a resolution of 180(cid:48)(cid:48). In addition tion with FWHM=180(cid:48)(cid:48). We also tested a higher resolution of tothefullrangeofcolumndensities,therelationshipswereex- FWHM=120(cid:48)(cid:48). For distant sources, the extinction of the target amined separately below and above the limit of τ(J)=0.6 (see clouds cannot be reliably reproduced because of the poor res- Sect.2.4).Thiscorrespondstovisualextinctions A ∼ 2.3mag V olution and the increasing number of foreground stars. This is andA ∼ 2.0magfortheR valuesof3.1and5.0,respectively V V the main factor that limits the number of fields where the ratio (Cardellietal.1989).Insteadofahigherlimit,weselectedthe τ(250µm)/τ(J) can be reliably estimated. The extinction mea- relativelylownumberofτ(J) = 0.6tomaximisethenumberof surements can be significantly biased even in nearby fields if fieldswherealinearfitcouldalsobemadeabovetheτ(J)thresh- these contain steep column density variations. No special steps old.Theτ(250µm)valueswerederivedfromHerscheldatawith were taken to eliminate the contamination by foreground stars either250–500µmor160–500µm(seeAppendixAforanalysis (see,e.g.,Schneideretal.2011),apartfromthesigmaclipping withadditionalalternativedatasets). procedurethatispartoftheNICERmethodandwasperformed Inagivenfield,thenumberofpointseitherbeloworabove at 3σ level. The reliability of the extinction maps and the bias theτ(J)limitisofteninsufficienttodetermineanyreliablevalue caused by sampling problems and the presence of foreground for the slope k = ∆τ(250µm)/∆τ(J). In a few fields no reliable stars was examined with simulations (see Sect. B). The results value of k can be determined at all, mainly because of the low ofthesesimulationsareusedtoderivemapsoftheexpectedun- quality of the τ(J) data. This especially affects the most dis- certaintyandthebiasoftheτ(J)valuesforeachfield. tantfieldsbecauseofthecontaminationbyforegroundstarsand because the structures are too small to be resolved with the 3(cid:48) beam.Theformalerrorsofthekparameterwereusedtoexclude 2.4. Correlationsbetweensub-millimetreandNIRopacity theclearlyunreliablefits.Thecriterionδk/k < 0.1leavesinthe The ratio k of sub-millimetre opacity τ(250µm) and the NIR defaultcase106fitstoalldatainafield,103fitsbelowτ(J)=0.6, opacity τ was estimated for all 116 fields. The τ(250µm) and 38 fits above τ(J)=0.6. These fits appear relatively reliable J maps were convolved to the 3(cid:48) resolution of the τ maps. The alsobasedonvisualinspection. J τ(250µm) and the τ(J) data were read at 90(cid:48)(cid:48) steps (half-beam Figure1showsanexampleoftherecovereddependencebe- sampling), excluding the map borders where the result of the tweenτ(J)andτ(250µm)values,includinglinearfitstothethree convolution to 3(cid:48) resolution is poorly defined. For local back- τ(J)ranges.Inthisexample,theslopeappearstobecomesteeper ground subtraction, only areas where the signal was more than as τ(J) increases.This mightbe anindication ofan increasein 2σ above the average value of the reference area were used thedustsubmillimetreopacity,whichinturnmightbeattributed (see2.2.2).Hereσisthestandarddeviationofthevaluesinthe tograingrowth(e.g.,Ossenkopf&Henning1994;Stepniketal. referenceregion.Thisisaconservativelimitbecausepartofthe 2003; Ormel et al. 2011; Ysard et al. 2013). However, before fluctuations is caused by real surface brightness variations and drawinganysuchconclusions,wemustconsiderthesystematic notbynoisealone. effectsthataffectthetwoparameters.Figure2showsasummary of all the ∆τ(250µm)/∆τ(J) values where τ(250µm) values are 4 http://horus.roe.ac.uk/vsa/index.html basedonHerschel250-500µmdata.Beforeanybiascorrections 4 M.Juvelaetal.:GalacticcoldcoresV.Dustopacity Fig.1. Relation between τ(250µm) and τ(J) in the field G95.76+8.17. The black solid line is a linear weighted total least-squares fit to all data points. The blue and red points and linesofthecorrespondingcolourshowthedataandthefitsbe- low and above the threshold of τ(J) =0.6, the dashed line in- dicates the division. The values of the slopes k are given in the plot.Errorbarsareshownforasetofrandomdatapoints. Fig.3. Comparison of ∆τ(250µm)/∆τ(J) values in three τ(J) (seebelow),thevaluesareseentoclusteraround∼ 2.0×10−3, intervals (three frames), obtained using either three of four Herschel bands in deriving the τ(250µm) values. The two his- withsometendencyforhighervaluesinthehigherτ(J)range. togramswithouthatching(thickoutlines)showallfieldswhere Figure 3 summarises the statistics of the dust opacity mea- the estimated uncertainty is below 10%. The two hatched his- surementsashistograms,includingallfitswheretheformalerror tograms contain only the intersection with better than 10% ac- of the slope of the least-squares fit τ(250µm) vs. τ(J) is below curacywithboththreeandfourbands(79,83,and14fieldsfor 10%. Weneedtoobserveasufficientnumberofbackgroundstars thethreepanels,respectively). foreachresolutionelement,evenathighcolumndensities.This means that τ(J) estimates and the comparison with submillimetre emission can be made only at a low resolution (2–3(cid:48)). Averaged over such large areas, the statistical uncertainty of Herschel data is very small. In Sect. 3.2 we show that the bias isprobablyalsodominatedbytheerrorsinτ(J).Inthefollow- pecially towards column density peaks. We estimated the ex- ing,werelymainlyontheHerscheldatasetthatconsistsofob- tent of the problem with simulations using the stellar statistics servations 250–500µm (the “default” data set). There are three in low-extinction areas near each field. The contamination by reasons. First, in theory, the inclusion of the 160µm data re- foreground stars was evaluated with the help of the Besançon ducesstatisticaluncertaintyofthecolourtemperatureestimates modeloftheGalacticstellardistribution(Robinetal.2003).We butincreasessystematicerrorscausedbyline-of-sighttempera- used the Herschel column density maps as a model of the col- turevariations(Shettyetal.2009;Malinenetal.2011).Second, because of the smaller (and, for parallel mode, different) area umn density structure, simulated the distribution of foreground andbackgroundstars,analysedthesimulatedobservationswith covered by the PACS observations, the use of the 160µm band NICERroutine,andcomparedtheresultswiththeknowninput significantly reduces the area where correlations with τ(J) can becalculated.Third,160µmdatamaybeaffectedbyadditional τ(J) map. The procedure is described in detail in Appendix B. systematic effects related to the relative calibration of the two Weobtainedforeachfieldamapoftheexpectedsystematicrel- ative error in τ(J) that gives a multiplicative correction factor instruments,uncertaintiesinthezero-pointdetermination(inter- fortheτ(J).Thesimulationsdonotconsidertheeffectofcloud polation between IRAS and Planck channels and the contribu- structures at scales below 18(cid:48)(cid:48) but the procedure probably pro- tion of stochastically heated grains in the IRAS 100µm band) videsareasonableestimateofthemagnitudeoftheeffect. andtoimperfectionsinthemapmakingthatcouldbeincreased by the smaller size of the PACS maps (see Sect. A.2). We are Werepeatedtheanalysisoftheprevioussectionandreplaced particularlyinterestedinthecoldestregionswhere250–500µm the original τ(J) maps with bias-corrected estimates. Figure 4 dataprovideadequateconstraintsonthedusttemperature. comparesthe∆τ(250µm)/∆τ(J)distributionsforthedefaultcase with and without bias correction. The statistics include all fits forwhichtheformalerroroftheslopeisbelow10%.Thiscor- 3.2. Biasinτ(J)values respondstoFig.3,butthenumberofpointsisdifferent.Because Bias in τ(J) values is very likely a significant problem, espe- thebiascorrectionsdependontheclouddistance,fieldswithout cially for distant fields in which all high column density struc- distanceestimateshadtobedropped.However,intheτ(J)>0.6 tures are not resolved and the results begin to be affected by interval the number of fields fulfilling the δk/k < 0.1 criterion foreground stars. Both effects decrease the τ(J) estimates, es- hasdoubled. 5 M.Juvelaetal.:GalacticcoldcoresV.Dustopacity Fig.2. Slopesk = ∆τ(250µm)/∆τ(J)forallcaseswithuncertaintiesδk/k < 0.1.Theblack,blue,andredsymbolscorrespondto valuesderivedforthefullτ(J)rangeandfordatabelowandabovethelimitofτ(J)=0.6.Thevaluesofτ(250µm)havebeenderived fromSPIREdatawithoutthesubtractionofthelocalbackground.Neitherτ(250µm)norτ(J)hasbeencorrectedfortheexpected bias. temperatureand,subsequently,tounderestimationofτ(250µm) (e.g.Ysardetal.2012). We constructed for each field a three-dimensional radiative transfer model that covered a projected area of 30(cid:48) × 30(cid:48) witha10(cid:48)(cid:48)pixelsize.Themodellingassumedspatiallyconstant dustproperties,thedustmodel(Draine2003b)correspondingto R =5.5 (see Appendix C for details). The density distribution V and external heating were adjusted until the model exactly reproduced the observed 350µm surface brightness and, for the area above median column density, the average 250µm/500µm ratio. The model-predicted surface brightness maps were anal- ysed as in the case of the actual observations, to produce maps ofτ(250µm).Toestimatethebias,thesevalueswerecompared totheactualτ(250µm)valuesofthemodeltoderivemultiplica- tivecorrectionfactors. The results depend on the assumed cloud structure in the line-of-sight direction (Juvela et al. 2013). In our models, the line-of-sight density distribution only has one peak. This en- hances temperature contrasts and increases our bias estimates. Ontheotherhand,thedensestobservedcoresareprobablyeven morecompact,andintheircasewemaybeunderestimatingthe bias.Ifthecloudscontainembeddedsources,theactualbiasmay againlocallybeverydifferentandoftenlowerthanpredictedby our models. Although the bias estimation is more difficult than Fig.4. Comparison of ∆τ(250µm)/∆τ(J) distributions without forτ(J),themodelsshouldagainprovideareasonableestimate bias corrections (“Default”) and with bias corrections applied ofthemagnitudeoftheeffect.Therelativesystematicerrorsare either to τ(J) or to both τ(J) and τ(250µm). The three frames smallerinτ(250µm)thaninτ(J)sothattheireffectonthefinal correspondtodifferentrangesofτ(J)values. resultislessstrong. Thek=∆τ(250µm)/∆τ(J)valueswerere-calculatedinclud- ingbiascorrections inbothvariables.Theresulting histograms are included in Fig. 4. The bias correction makes the distribu- 3.3. Biasinτ(250µm)values tionssignificantlynarrower.Figure4alsoshowsthatthecorrec- tionsaremuchstrongerforτ(J)thanforτ(250µm).Asaresult, The systematic errors in τ(250µm) values were estimated with the average value of ∆τ(250µm)/∆τ(J) now decreases with in- radiativetransfermodelling.Theline-of-sighttemperaturevari- creasingτ(J).Thenumberoffieldsfulfillingtheδk/k <0.1cri- ations are expected to be the main source of error that, in stan- terionhasdoubledto76fields(usingSPIREbands).Compared dardanalysis,leadstooverestimationofthemass-averageddust totheoriginaldata,themedianvaluesofk×104havedecreased 6 M.Juvelaetal.:GalacticcoldcoresV.Dustopacity Fig.5. Slopevaluesk = τ(250µm)/τ(J)ofFig.2asafunctionofestimateddistance(framesa,b),galactocentricdistance(frames c,d),andgalacticheight(framese,f).Theleftframesshowtheoriginalslopes,therightframestheslopesafterthebiascorrections of τ(250µm) and τ(J). The black, blue, and red colours correspond to the full τ(J) range and to data below and above τ(J)=0.6, respectively.Thesolidcurveswiththesamecoloursaretheweightedmovingaverages(windowsizes30%indistance,800pcin Galactocentricdistance,and100pcinGalacticheight).Allτ(250µm)valuesarecalculatedwithSPIREbandsalone. from 20.2, 21.1, and 23.4 to 15.3, 16.0, and 12.2, the numbers in the corrected data there is a slight decrease in k values as a correspondingtothefullτ(J)range,databelowτ(J) = 0.6,and function of distance. This could point to some over-correction dataaboveτ(J) = 0.6,respectively.Thestrongchangeinthek of the τ(J) estimates, although the bias correction should not valuessuggeststhattheuncertaintyofkisprobablyoftenseveral onlydependondistance,butevenmainlyonthecloudstructure. tensofpercent,especiallyintheτ(J)>0.6interval.Therefore, However,thedecreaseofkcanbeanindicationofselectionef- Fig.4doesnotexcludeasystematicincreaseofkasthefunction fects or direct resolution effects. For example, higher k values ofτ(J),ifthatbecomesvisibleonlyathighcolumndensities. mightbefoundinindividualdenseclumpsthatareonlyresolved atshortdistances. Figure5displaystheslopesk = ∆τ(250µm)/∆τ(J)asfunc- There is little difference between the k values found in the tionofdistanceandGalacticlocation.Theleftframesshowthe threeτ(J)intervals.Inthenextsectionweexamineinmorede- relations without and the right frames with the bias corrections tailtheglobalτ(J)dependenceoftheτ(250µm)/τ(J)ratios,es- appliedtoτ(J)andτ(250µm).Onlyfitswithδk/k < 0.1arein- peciallyregardingthehighestobservedcolumndensities. cluded.Theoriginaldatashowedsometrends,includinganin- creaseinkasafunctionofdistanceandgalactocentricdistance. The first is visible especially in the high τ(J) interval, but was 3.4. Globalrelationτ(250µm)vs.τ(J) expectedbecauseτ(J)valuesofdistantsourcescanbeseverely underestimated.Afterbiascorrections,thescatterofthekvalues To further test the hypothesis that τ(250µm)/τ(J) ratios may issignificantlyreduced.Thissuggeststhatthecorrectionsareof change systematically as a function of column density, we car- correctmagnitude.Thedistancedependencehaschangedsothat ried out non-linear fits τ(250µm)= A+B×τ(J)+C×τ(J)2. 7 M.Juvelaetal.:GalacticcoldcoresV.Dustopacity The fits were first performed using the combined data of all fields. To reduce the mismatch in the zero levels of individual fields,wesubtractedfromeachτ(J)andτ(250µm)mapthelo- calbackgroundusingtheoffregionslistedinTable1.Theoffre- gionsarenotcompletelyvoidofemissionbutprovideacommon referencepointforthequantities.Thus,therelationisexpected todevelopviatheoriginforeachfieldseparately,theparameter Abeingclosetozeroforthecombineddataaswell.Thesignof thefittedparameterCindicatesthepossibleincreaseordecrease ofk=∆τ(250µm)/∆τ(J)asthefunctionofcolumndensity. Figure 6 shows the results obtained with the bias-corrected data. We included data from all fields in which individual lin- Fig.6. Fitofτ(250µm)= A+B×τ(J)+C×τ(J)2 tothecom- ear fits had δk/k < 0.2, thus relaxing the previous constraint bined data of all fields in which individual linear fits showed a of δk/k < 0.1. The second-order polynomial was fitted to all strongcorrelationwithδk/k <0.2.Theblueandthegreenlines dataandseparatelytodatapointswithτ(J) > 1.0.Intheprevi- correspond to fits to the full column density range and to data oussectionathresholdvalueofτ(J) = 0.6wasused.However, pointsτ(J) > 1.0alone,respectively.Inthesecondframe,only some70%ofallpointsarebelowτ(J)=0.6and,whenincluded, datawithcolourtemperaturesbelow14Kareused. theydominatethefitsthatsystematicallyunderestimatethedata above τ(J) ∼ 5. With the combined data set, there are enough high column density data points so that the lower limit can be movedupwards.Theuseoftheτ(J) = 1.0thresholdenablesan adequatefittoalldatawithhigherτ(J)values.Theτ(J)calculations employ a different off region for each field. These may contain different amounts of extinction, which leads to small relative shifts along the τ(J) axis. Based on dust emission, the extinction in the off regions is typically τ(J)=0.2–0.4. The uncertainty of the relative zero points contributes to the scatter in Fig.6,buttheeffectisweakerthanthetotaldispersionandthe non-linearityseenathighextinctions. To prevent the τ(J) cut itself from biasing the fits, the data wereselectedusinglinesperpendiculartoalinearleast-squares linefittedtoalldata(cf.Sect.2.4).Thus,thequotedτ(J)limits correspondtoapointonthefittedline,andthecutitselfisper- pendiculartothefittedline.Allfitstakeintoaccounttheuncer- taintiesinbothvariables,whichareassumedtobeuncorrelated. The error distributions of the parameters A-C were calculated withanMCMC. Fig.7. Distributions of the parameters of the fit τ(250µm) = The sign of the parameter C depends on the range of τ(J) A+B×τ(J)+C×τ(J)2.Thefitislimitedtodatawithτ(J)>1. values but is less dependent on the field selection, for example regardingtheδk/k limitthatwasusedtoselectthefields.Most fields with high τ(J) values (and thus with a wide dynamical range) have low values of δk/k. When all points are included, column densities the relation is linear, the curvature increases thevalueoftheparameterC isnegative,butthefitisverypoor onlybeyondτ(J)∼5,andthelowestdatapointsτ(J)<1.0are at high τ(J). When pixels τ(J) < 1.0 are excluded,C becomes notpartofthefit. positive (see Fig. 6a), which points to an increase in the submillimetre opacity above τ(J) ∼ 1. Systematic additive errors ineitherparametermightalsoexplainthedifferentbehaviourat very low τ(J). When the fit is made using all data τ(J) > 0.6 (notshown),theparameterC ismarginallypositive,butbeyond τ(J)=10thefittedlineisbelowallthedatapoints.Thesecond Figure7showstheerrordistributionsoftheparametersA–C. frameofFig.6showsthefitswhenpixelswithcolourtempera- Thefitwasmadetodataτ(J)>1.0forallfieldswithalinearfit turesabove14Kareexcluded.ThevaluesofC arenowhigher accuracy δk/k < 0.2. In most fields, the formal error estimates and positive even when data τ(J) < 1.0 are not excluded. The of δτ(J) and δτ(250µm) are smaller than the actual scatter of best fit to the high τ(J) end of the relation is still obtained by points. Therefore we used the residuals of the linear fits before excludingdatawithτ(J)<1.0,thisresultsintherelation the MCMC calculation to determine a scaling factor, typically 2.0–3.0, that makes the error estimates in each field consistent τ(250µm)=0.73×10−3+1.25×10−3τ(J)+0.11×10−3τ(J)2.(2) withtheactualscatter.Evenafterthisincreaseofuncertainties, TheformalerrorestimatesoftheparametersA-Careoftheorder MCMC gives a 100% probability for a positive value of C. In of5%,probablylowerthanthesystematicuncertainties.Alldata reality,theresultisnotthatstrongbecausetheuncertaintymay beyondτ(J)∼5areaffectedbylargebiascorrectionsand,con- be dominated by systematic errors. The sign of C was already sequently,thevalueofC alsodependsontheaccuracyofthese seentochangedependingontherangeofτ(J)valuesfitted.The corrections.Thus,Fig.6stronglysuggestsbutdoesnotyetpro- result also depends on a relatively small number of fields with videafinalproofofthevariationsoftheratioτ(250µm)/τ(J). dataaboveτ(J)>5.Therefore,wemustconsidertheτ(250µm) ThepositiveoffsetA=0.7×10−3resultsfromthefactsthatatlow vs.τ(J)relationinindividualfieldsinmoredetail. 8 M.Juvelaetal.:GalacticcoldcoresV.Dustopacity 3.5. Correlationsinselectedfields Figure6showedhintsofanincreaseoftheτ(250µm)/τ(J)val- uesasthecolumndensityincreases,buttheglobalstatisticsmay beconfusedbythemixofdifferentfields.Furthermore,thesign oftheC parameterisdeterminedbythehighestτ(J)pointsthat originateinasmallnumberofindividualfields.Eachfieldmay beaffectedbydifferentsystematiceffectsrelatedtothesurface brightnesszeropoints,distanceuncertainty(viabiascorrection), anddifferencesinthelocalradiationfield.Severaldiffusefields are even entirely below τ(J) ∼ 1.0. Therefore, we also need to examinethefieldsindividually. Three criteria were used to select a subset of fields. We re- quiredthat(1)theuncertaintyofthefittedparameterCisbelow 0.3×10−4, (2) there are at least ten data points (selected from themapsat90(cid:48)(cid:48) steps)withτ(J)above0.6,and(3)thebiascor- rectionschangetheslopeofthelinearfitofτ(250µm)vs.τ(J)by Fig.10. Linear slope k (upper frame) and the parameters C lessthan30%.Thefirsttwocriteriaensurethatthereareenough (lower frame) for the 23 selected fields. The values obtained datapointsatlargeτ(J)withasmallscattertogainsomeinsight withoutbiascorrectionsareshownwithblacksymbols.Theval- aboutthecolumndensitydependenceofthe∆τ(250µm)/∆τ(J). ues obtained with corrected τ(J) and τ(250µm) data are shown The third criterion excludes distant fields for which the uncer- with red symbols, the shaded area corresponding to the uncer- taintyofthebiascorrectionofτ(J)renderstheresultsuncertain, taintyofthebiascorrectionthatisduetotheuncertaintyofthe even for the apparently well-defined relation between τ(J) and distance estimates. The dashed lines show the median values τ(250µm). The selection leaves 23 fields with distances mostly corresponding to the black and red symbols. The fields are ar- intherange100-500pc.Therearetwoexceptions,G216.76-2.58 rangedinorderofincreasingdistance,andtheτ(250µm)values andG111.41-2.95,forwhichtheestimateddistancesare2.4kpc arebasedonSPIREdataalone. and3.0kpc.Wekeptthetwofieldsinthesampleeventhoughthe resultsareknowntobeunreliablebecauseofthelargedistance. To avoid underestimating the fit errors, we scaled the error did not change appreciably (by less than 0.1×10−3). The me- estimates of τ(J) and τ(250µm) up to correspond to the actual dian value of C increased to 2.0 × 10−4, which in that case is scatterofpoints(seeSect.3.4).Allobservations,sampledat90(cid:48)(cid:48) still lower than the scatter. The differences between the fields steps,arefittedwithalinearmodelτ(250µm)=b+k×τ(J)us- arelargerthantheestimatedformaluncertainties(includingthe ing total least-squares. We continued to use as the default data statisticalerrorsofτ(J)givenbyNICERandτ(250µm)derived setonewithτ(250µm)derivedfromSPIREbandsalone,includ- from the uncertainty of the surface brightness measurements). ing bias corrections in τ(250µm) and τ(J). However, for com- Theerrorbarsonlyreflectthestatisticalerrorsofthefitsanddo parison, we also examined results obtained with four Herschel not include the uncertainty of the bias correction, for example. bands(160–500µm;includingthebiascorrections)and,finally, Figure 8 demonstrates that the uncertainty of the distances can withthreeHerschelbandsbutwithoutanybiascorrections. beasignificantsourceoferror.InFig.10,theshadedareasshow The non-linear fits were made using MCMC (with 2×105 thedifferencebetweenthevaluesobtainedwithdistancesd−δd samplesperfield)andbootstrapsampling(2000realisationsper andd+δd,aslistedinTable1.Thesewereestimateddirectlyby field).Bothfitsusedtotalleast-squares5 andtheerrorestimates repeatingtheanalysisusingthesetwodistancevalues.Asmaller oftheindividualpoints.Theparametervaluesandtheuncertain- distancecorrespondstosmallerbiascorrectioninτ(J)and,thus, ties derived with these two methods should be similar, except to a steeper slope k and typically a higher value ofC. In some forrarecasesinwhichtheresultdependsonasmallnumberof cases,thevalueofCobtainedwiththedefaultdistancedisout- influential points, which are always present in the MCMC cal- sidetheshadedregion,showingthattheeffectisnotalwaysthis culation,butnotinallbootstrapsamples. simple.Thedistanceuncertaintyisnotyetenoughtoexplainall Figures8-9showtheresultsofthesefits.Eachframeshows thescatterinkandespeciallyinC.Achangeinthedistancees- the values of the linear slopes for the three cases discussed timateresultsatfirstapproximationinanearlylinearscalingof above.Theparametersofthenon-linearfitsareshowntogether τ(J)values(seeFig.8,comparisonofthedashedmagentalines). with the error estimates calculated with the bootstrap method. Inreality,thesituationmaybemorecomplex.Inparticular,ifa Inadditiontothefittothedefaultdataset(solidredcurve),the field contains cloud structures at different distances, this might dashedmagentalinesshowtheeffectofthedistanceuncertainty. resultinlargeerrorsinbothkandC.Figure10alsoshowsthat Using the distance uncertainties δd listed in Table 1, we also in spite of the small formal errors of the least-squares fits, we calculatedthebiascorrectionsforτ(J)fordistancesd−δd and cannot constrain the opacity values in the last two fields with d+δd.Thus,theupperdashedlinecorrespondstodistanced−δd distancesexceeding2kpc. andasmallerbiascorrection. Figure 10 shows the linear slopes k and the values of pa- rameterC for this sample of fields. The fits were performed to In the sample of Fig. 10, the median value ofC is close to all data, without a τ(J) threshold. If the data below τ(J) = 0.6 zerowithanumberoffieldswithnegativevalues. Thepositive were removed, the median slope τ(250µm)/τ(J) = 1.6×10−3 values of C in Fig. 6 are due to a small number of fields, and 5 Distancebetweena(τ(J),τ(250µm))pointandthemodelcurveis the increase of τ(250µm)/τ(J) values was only visible above measured in a coordinate system where the error region of the point τ(J) ∼ 5. There are only 20 fields with any data points above is circular. We used the smallest distance to the curve, ignoring the τ(J)=5.Onlysixfieldshavetenormoredatapointsabovethis marginaleffectresultingfromthecurvatureofthemodelcurve. limit: G6.03+36.73, G70.10-1.69, G82.65-2.00 G92.04+3.93 9 M.Juvelaetal.:GalacticcoldcoresV.Dustopacity G107.20+5.52,andG202.02+2.85.Ofthese,onlyG6.03+36.73 isincludedinthesampleofFig.10.Alltheotherswereexcluded because the bias correction changed the slope k by more than 30%. Thus, a clear steepening of the relation τ(250µm vs. τ(J) is seen exclusively in fields with the highest column densities, which for the same reason also have the largest uncertainty re- gardingthebiascorrections. 3.6. Mapsofτ(250µm)/τ(J)ratio Wealsoexaminedtheratiosτ(250µm)/τ(J)intheformofmaps. Thisisusefulifkchangesinsmallregionsthathavelittleeffect when all data of a field are fitted. Unlike in Fig. 8, where the offset between τ(250µm) and τ(J) is a free parameter, the ap- pearance of the ratio maps depends on the consistency of the τ(250µm)andτ(J)zeropoints.Becausewedonothaveanab- Fig.11. FieldG4.18+35.79(LDN134).Theupperframesshow solute zero point for τ(J), we used the reference areas listed τ(250µm)(framea)andtheratioτ(250µm)/τ(J)(frameb).The in Table 1 and subtracted from τ(250µm) and τ(J) the average lowerframesshowthelower(framec)andupper(framed)limits valuefoundinthereferencearea.Thislimitstheregionwherea ofτ(250µm)/τ(J)calculatedas(τ(250µm)+δτ(250µm))/(τ(J)− reliableratiocanbecalculated,excludingregionsoflowcolumn δτ(J))and(τ(250µm)−δτ(250µm))/(τ(J)+δτ(J)).Theareasnot density.ThedetailsofthecalculationsaregiveninAppendixE, coveredbyHerschelobservationsandregionswithaSNbelow where we also show the figures of selected fields. As an ex- 0.5 have been masked. In frame a, the solid black contour and ample,Fig.11showsthefieldG4.18+35.79(LDN134),where thedashedwhitecontourcorrespondtoτ(250µm = δτ(250µm) theratioτ(250µm)/τ(J)isstronglycorrelatedwithcolumnden- andτ(J) = δτ(J).Themapshavearesolutionof180(cid:48)(cid:48) andτ(J) sity. In Fig. 11, the increase of k remains clear even in maps isderivedusing2MASSdata. of (τ(250µm)±δτ(250µm)/(τ(J)±δτ(J)). The error estimates δτ(250µm) include the statistical errors due to Herschel pho- tometryanduncertaintyofthesurfacebrightnesszeropoint(see in the extinction maps. The analysis was repeated for the four Sects.2.2.3and2.2.2).Theparameterδτ(J)correspondstothe fieldswithVISTAdatawithresolutionsof180(cid:48)(cid:48) and120(cid:48)(cid:48) .The uncertaintyoftheτ(J)zeropoint(seeAppendixE).Forτ(J)the resultsaresummarisedinFig.12.Theresolutionhasnostrong formalerrorestimatescalculatedwithNICERarebelow10%. systematic effect on the parameters. The largest differences ap- For high-opacity sources like G4.18+35.79 the bias correc- pearwhenparameterC isestimatedexcludinglowcolumnden- tionofτ(J)isveryimportant.Ifnobackgroundstarsarevisible sity points. However, even in that case the difference between through some part of the core, the values of k naturally remain the results with a resolution of 120(cid:48)(cid:48) and 180(cid:48)(cid:48) is smaller than more uncertain (see, however, Sect. 3.7). Conversely, the zero- theeffectofexcludinglowτ(J)datafromthefits.WhenVISTA pointuncertaintyonlybecomesimportantindiffuseregionsbut data are available, the results are close to those obtained with mightevenreversethecorrelationwithcolumndensity.Theratio 2MASS.BecauseweusedthesameHerscheldataandbiascor- mapsarealsoaffectedbytheassumptionofaconstantvalueof rectionsderivedinthesameway,theresultsarenotindependent. βandofpotentialerrorsinthebiascorrections.However,wear- However,becausetheuncertaintyofτ(J)isexpectedtobeone gueinSect.3.7thatthesearemainlymultiplicativeerrors(that ofthemostsignificantsourcesoferror,thisgivesussomecon- do not affect the morphology of the maps) and/or tend to de- fidencethattheobserveddifferencesbetweenthefieldsarereal. crease the variations seen in the ratio maps. Therefore, we are The τ(250µm)/τ(J) ratio in the field G4.18+35.79 was confidentthattheincreaseofsubmillimetreopacitythatisseen shown in Fig. 11. Figure 13 shows the corresponding figure insomeofthemapsisreal. obtained with VISTA NIR data and a spatial resolution of Based on the maps, the submillimetre opacity is correlated 120(cid:48)(cid:48). The highest value of τ(250µm)/τ(J) has increased from with column density in the fields G4.18+35.79, G6.03+36.73, 3.6×10−3 inFig.11to6.7×10−3.Thisismainlyattributedto G111.41-2.95,G161.55-9.30,G151.45+3.95,andG300.86-9.00 theincreasedspatialresolution,althoughalsoatthe180(cid:48)(cid:48) reso- (seeAppendixE).InG4.18+35.79andG6.03+36.73thevalues lution the peak value is ∼25% higher than in Fig. 11. Even in rise to close to 4 × 10−3. In some fields the background sub- VISTAdatathereareonly∼10starswithinthe2(cid:48)×2(cid:48)areacen- tractionreducestheavailablemapareatosuchanextentthatno tredontheτ(250µmmaximum,andthereforethepeakvalueof conclusionscanbedrawn. τ(250µm)/τ(J)issubjecttosomeuncertainty. 3.7. Potentialsystematicerrors Sofar,allNIRextinctionmapswerecalculatedat180(cid:48)(cid:48) res- Because of the significance of the bias corrections (Sects. 3.2 olution. Depending on the number of background stars, extinc- and 3.3), we tried to characterise the effects that systematic er- tion map could be derived at a higher resolution and possibly rors in these corrections could have on the τ(250µm)/τ(J) ra- with smaller bias. This especially applies to the four fields for tios. Furthermore, the assumption of a constant dust emissivity which VISTA observations are available. We recalculated the spectralindexmaybeincorrect.Belowweexaminethepossible extinctionmapsat120(cid:48)(cid:48) resolution,repeatingMonteCarlosim- systematiceffectscausedbythesefactors. ulationstoestimatethebiasofτ(J).Thesmallerbeamincreases The bias correction made to the τ(250µm) values is in it- the noise per resolution element, but does not yet cause holes selfsmall(seeFig.4)and,consequently,theerrorsmadeinthat 10

Galactic cold cores V. Dust opacity

Similar Galactic cold cores V. Dust opacity