The JCMT Gould Belt Survey: radiative heating by OB stars

Radiative feedback can influence subsequent star formation. We quantify the heating from OB stars in the local star-forming regions in the JCMT Gould Belt survey. Dust temperatures are calculated from 450/850 micron flux ratios from SCUBA-2 observations at the JCMT assuming a fixed dust opacity spectral index $\beta=1.8$. Mean dust temperatures are calculated for each submillimetre clump along with projected distances from the main OB star in the region. Temperature vs. distance is fit with a simple model of dust heating by the OB star radiation plus the interstellar radiation field and dust cooling through optically thin radiation. Classifying the heating sources by spectral type, O-type stars produce the greatest clump average temperature rises and largest heating extent, with temperatures over 40 K and significant heating out to at least 2.4 pc. Early-type B stars (B4 and above) produce temperatures of over 20 K and significant heating over 0.4 pc. Late-type B stars show a marginal heating effect within 0.2 pc. For a given projected distance, there is a significant scatter in clump temperatures that is due to local heating by other luminous stars in the region, projection effects, or shadowing effects. Even in these local, `low-mass' star-forming regions, radiative feedback is having an effect on parsec scales, with 24% of the clumps heated to at least 3 K above the 15 K base temperature expected from heating by only the interstellar radiation field, and a mean dust temperature for heated clumps of 24 K.


INTRODUCTION
Radiative feedback from young, higher-mass stars is thought to affect subsequent star formation by heating the gas and reducing fragmentation (Bate 2009;Howard et al. 2016).
We have been exploring the observational evidence for this effect in a series of papers looking at heating in local star-forming regions observed as part of the JCMT Gould Belt Survey (Ward-Thompson et al. 2007), specifically Perseus NGC1333 (Hatchell et al. 2013), Serpens MWC297 (Rumble et al. 2015) and the Aquila W40 complex (Rumble et al. 2016)). Using dust colour temperatures derived from submillimetre ratios, we have found evidence that OB stars heat the dust in nearby filaments and cores and potentially increase their stability against gravitational collapse through raising the Jeans mass.
Although these nearby regions are thought of as low-mass star-forming regions (with the exception of Orion A), there are O or B stars associated with several of them. Because only one or two OB stars are associated with each region, it is simpler to consider the influence of those stars on the surrounding cloud material than in high-mass star-forming regions, where the heating effect of several stars is combined.
In this study, we extend our investigation of radiative feedback to include all clouds in the local star-forming regions within the James Clerk Maxwell Telescope Gould Belt Survey (JCMT GBS), considering the dust heating as a function of stellar spectral type and projected distance from the OB star.

METHODS
Maps of dust colour temperatures based on SCUBA-2 450 µm and 850 µm maps from the JCMT Gould Belt Survey Data Release 1 (DR1) were created for the 30 GBS subregions listed in Table 1. For details of the observations (2012 to 2015), data reduction, and maps, see Kirk et al. (2018) and earlier survey publications (Rumble et al. 2015;Pattle et al. 2015;Mairs et al. 2015a;Chen et al. 2016;Kirk et al. 2016;Rumble et al. 2016). SCUBA-2 maps contain spatial scales up to 5 (0.75 pc at 500 pc) so they trace the dense cores and filaments most directly associated with star formation without the large-scale cloud background (Ward-Thompson et al. 2016).
Dust temperatures T d were calculated pixel-by-pixel from   Abt & Levato (1977); 9 - Strom et al. (1974);10 -Herbst (2008);11 -Walawender et al. (2008);12 -Drew et al. (1997);13 -Houk & Swift (1999);14 -Smith et al. (1985) the ratio of SCUBA-2 450 µm (S450) and 850 µm (S850) fluxes following Reid & Wilson (2005) in using where β is the dust opacity spectral index (Beckwith et al. 1990). The 450 µm map was first convolved to match the 850 µm resolution using the Aniano et al. (2011) kernel convolution algorithm adapted for SCUBA-2 images by Pattle et al. (2015) and Rumble et al. (2016) with an analytical model-beam kernel based on the two component beam model of Dempsey et al. (2013); see Rumble et al. (2016) for details. The resulting dust temperature maps have an effective resolution of 14.8 , comparable to the JCMT 850 µm effective beam FWHM of 14.6 . Temperature estimates are limited to the cores and filaments detected at 450 µm signal-to-noise ratios (SNR) greater than 5 and pixels with temperature uncertainties < 34.1% (Rumble et al. 2016). This is a choice in favour of reliability over area coverage: relaxing these requirements mostly adds hot, high-uncertainty pixels around the periphery of the clumps.
An example temperature map for Ophiuchus L1688 is shown in Fig. 1; the full set of dust temperature maps are publicly available from https://dx.doi.org/10.11570/20. 0007 1 . A constant β of 1.8 was assumed, consistent with most dust opacities measured in Perseus from joint Herschel/SCUBA-2 fits (Chen et al. 2016;Sadavoy et al. 2013), Planck observations (Juvela et al. 2015) and that used in other GBS papers (Hatchell et al. 2005, Salji et al. 2015, Rumble et al. 2015and Pattle et al. 2015. Using a fixed β avoids the inherent anticorrelation between T d and β when fitting fluxes (Chen et al. 2016;Shetty et al. 2016) at the expense of temperature errors when the assumed β is wrong (for example in protostellar discs where β ∼ 1 is a better estimate). The choice of β = 1.8 is also consistent with the Ossenkopf & Henning (1994) (Lynds 1962). The OB stars S1 and HD 147889 are marked as cyan stars. Also marked are the young stellar object candidates from Dunham et al. (2015) (green circle = Class 0, green plus = Class I, red cross = Class II, red star = Class III). The 850 µm SNR=3 contour is shown in black (0.26 mJy/arcsec 2 ) solute dust temperature scale throughout this work depends on the choice of β as shown in Fig. 2, with lower β values giving higher temperatures for the same observed flux ratio -for example, the median clump temperature derived in Sect. 3 would increase to 18.2 K for a lower assumed β = 1.6 or reduce to 12.9 K for β = 2.0, and this spread increases at higher temperatures. Simultaneous determinations of β and temperature are now available for a few local clouds based on joint SCUBA-2 and Herschel analysis (Sadavoy et al. 2013;Chen et al. 2016;Howard et al. 2019Howard et al. , 2021; but until these are consistently available for many clouds, this simplifying assumption is necessary. The dust colour temperature obtained by solving Equation 1 is an average over the line-of-sight components weighted by the emission at 450 µm which is longward of the peak of the spectral energy distribution even for dust as cold as 10 K. In the 450/850 ratio, there is a relatively high contribution from the colder dust (< 20 K) that makes up the bulk of the mass of (starless) cores and filaments. In addition, SCUBA-2 maps spatially filter out scales greater than 5 (Mairs et al. 2015b;Kirk et al. 2018), so it is insensitive to extended emission, for example, emission from warm diffuse dust heated by clusters or low-extinction cloud layers. Nonetheless, emission-weighted line-of-sight average dust temperatures are always higher than the equivalent massweighted temperatures.
The SCUBA-2 850 µm regions were split into clumps using the Starlink fellwalker algorithm (Berry 2015) following Rumble et al. (2016) with an 850 µm SNR threshold of 5. Clumps smaller than the beam were rejected, eliminating many Class II discs (which we wish to exclude as their temperatures are dominated by internal heating and for which our assumed β is incorrect) as well as artefacts. In total, 1134 clumps were detected across the 30 subregions listed in Table 1 (743 excluding Orion A). Clump sizes along each of the pixel axes were calculated by fellwalker as the root-meansquare flux-weighted offset of each pixel from the clump centroid, equivalent to the Gaussian width for a Gaussian clump, from which we take the geometric mean (for details, see Berry 2015). The beam-deconvolved flux weighted clump sizes range between 0.01 pc and 0.37 pc with a median of 0.08 pc, larger than a typical protostellar core (0.05 pc, Rygl et al. 2013), confirming that the fellwalker clumps represent more extended, filamentary structures. Average clump temperatures were calculated for the subset of 465/743 clumps for which at least 50% of the pixels within the clump had a reliable temperature estimate (i.e. a 450 µm SNR > 5 and a temperature uncertainty of less than 34.1%).
Although it is the best known massive star forming region in the Gould Belt, Orion A is excluded from our analysis for two reasons. Firstly, temperatures around OMC-1 are so high that the SCUBA-2 450 µm fluxes lie on the Rayleigh-Jeans tail and so Equation 1 cannot be used to calculate reliable temperatures. Secondly, there are so many OB stars in Orion that identifying the one(s) heating each clump and assigning a projected distance is not trivial.
The OB (and A) star(s) potentially heating each subregion were compiled based on the Reed (2003) Catalogue of Galactic OB Stars, rejecting line-of-sight contaminants based on proper motions from the SIMBAD database (Wenger et al. 2000). For regions with multiple OB stars a major OB star was identified (see Table 1): these were selected primarily on spectral class, favouring the most massive stars, but also proximity to any clumps, as in the cases of Ophiuchus L1688, Perseus NGC1333, and NGC2023/24 where lower mass OB stars are directly heating their immediate environment in addition to more distant, higher mass OB stars (Pattle et al. 2015;Hatchell et al. 2005;Meyer et al. 2008). We calculated the projected distance (in parsecs) of each SCUBA-2 clump to its nearest major OB star. Figure 3 shows that the 16 subregions with OB stars contain larger areas of warm dust (20-50 K) than those without OB stars.

RESULTS AND ANALYSIS
In regions with OB stars, the distribution is shifted to higher temperatures with a median temperature of 21 K compared to 15 K for regions without OB stars. Regions with OB stars have a factor 3.3 more pixels above 20 K than regions without OB stars.
We further subdivide the 16 regions with OB stars into three sets -O-type, Early B-type (B0-B4) and Late B-type (B5-A0) -based on the classification of the major OB star in the region (see Table 1 for assignments). A0 is included in the Late B class in order to assess the influence of V892Tau in Taurus L1495. The distributions of pixel temperatures for each subset, and for the regions individually, are shown in Fig. 3 (second panel from top and lower panels). Median temperatures are 33 K, 21 K, 18 K and 15 K for regions with O, early B, late B-type, and no OB stars, respectively. Regions with O, early B, or late B stars show greater fractions of warm dust (20-50 K) than regions with no OB stars, with the greatest fraction in the two regions with O stars (Orion B NGC2024 and W40).
The variation of clump temperature with projected distance from the region's major OB star is shown in Figure 4. As expected, more of the clumps in regions with O stars are warm (20-50 K), consistent with the temperature distributions in Figure 3. All subsets show the anticorrelation of clump temperature with projected distance identified in the Aquila W40 complex (Rumble et al. 2016). Nevertheless, there is some scatter around this trend: projected distances that are shorter than distances along the line-of-sight will produce clumps that are colder than expected and clumps that are locally heated by stars other than the major regional ones listed in Table 1 will be warmer.
Following Draine (2011), we model the variation in dust temperature by considering a dust parcel of mass M that is heated by the incident power from an OB star, ΓOB, plus the interstellar radiation field (ISRF), ΓISRF, and cooled by the power radiated by the dust, Λ cool . In thermal equilibrium, Λ cool = ΓOB + ΓISRF. (2) The OB star heating varies as the inverse square of distance d and consists of the integrated flux absorbed by the parcel over all wavelengths, where L * (ν) is the OB star luminosity, κν is the dust emissivity responsible for the absorption and M the mass. A layer of extinction between the star and clump is factored as e −τν (Av ) where τν (Av) is the optical depth as a function of visual extinction Av.
Additionally, an incident interstellar radiation field (ISRF) with mean intensity JISRF(ν) heats the dust mass from all directions. Assuming no local attenuation of the ISRF, the absorbed power is given by The dust cools radiatively in all directions following the Stefan-Boltzmann law, where T d is the dust temperature and 4 κν M is the effective surface area of the dust mass. κν is the emission-averaged opacity weighted by a blackbody at temperature T d .   . Clump temperature as function of projected distance from the major OB star in the region (red circles: clumps in regions with O stars; blue triangles: early-type B stars; cyan squares: late-type B stars). The solid lines are OB + ISRF heating fits as described in the text (colours as data). The dotted lines show the expected temperatures vs. distance for clumps heated by stars at the spectral type boundaries with extinction Av = 0.5 (B0: red dots, B5: blue dots; A0: cyan dots); the heated clumps in each spectral type class should lie above these. Low temperature points excluded from the fits are shown as small '+' . See Sect. 3 for details of models and fit.
In this formulation, the mass M appears on both sides of the equation and cancels out. The integrals are carried out over a range in frequency from the equivalent of centimetre wavelengths (for which dust opacity is negligible) to the Lyman limit (91.2 nm), assuming all the shorter-wavelength UV emission is absorbed by hydrogen atoms.
We assume Ossenkopf & Henning (1994) model 5 dust as appropriate for dense molecular clouds. The long-wavelength spectral index of β = 1.8 gives a Planck-averaged dust opacity that depends on temperature as κ T d 10 cm 2 g −1 × (T d /10 K) β . We use the same dust properties for the heating/cooling balance and for the extinction layer. For the ISRF intensity spectrum, we use the tabulated version from Evans et al. (2001) which has contributions from the microwave background (Wright et al. 1991), IR/visible/UV (Mathis et al. 1983) and short UV (Draine 1978). With these choices, the power absorbed from the ISRF and the cooling factor per unit mass are ΓISRF = 2.33 × 10 −2 W and Λ cool = CT 4+β with C = 3.59 × 10 −9 W K −5.8 .
The heating/cooling balance gives a temperature of 15 K for clumps only heated by the ISRF. This value is within the range of steady-state dust temperatures calculated by Mathis et al. (1983). It is also consistent with the median temperature of clumps in regions without OB stars (see Fig. 3), where the emission-weighted dust temperatures are dominated by the ISRF-heated outer layers of clumps. Gas temperatures are expected to be a few Kelvin lower than dust temperatures in these outer layers as the gas is heated primarily by collisions with the dust (Forbrich et al. 2014). The strength of the ISRF varies from region to region depending on the local stellar distribution and extinction (Bernard et al. 2010), and regions with significant numbers of clumps and cloud area below the 15 K median indicate a weaker ISRF (see Fig. 3 and Fig. 4).
We calculate the OB heating power ΓOB for a range of stellar spectral types and extinctions Av. Table 2 gives the values of ΓOB for unextincted OB stars and OB stars with an extinction layer Av = 0.5. Combining the heating from OB stars with that from the ISRF as above, ΓOB + ΓISRF, the dust temperature was calculated using Equations 2 and 5. Results for B0, B5 and A0 spectral types with an extinction layer Av = 0.5 are shown for comparison with the data in Figure 4. Unextincted OB stars produce temperatures that are too high by many Kelvin. The addition of a layer of extinction Av ∼ 0.5 between the OB star and the clump, however, reduces the dust temperature to moderate levels consistent with those observed. The impact of this assumption is similar to that of the AV = 1 layer assumed by Galli et al. (2002) when modelling prestellar cores. Alternatively, these results suggest that the emission-weighted dust temperatures are typical of dust that is at an extinction of AV = 0.5 within the clumps, consistent with results for isolated dense cores such as CB68 (Nielbock et al. 2012;Lippok et al. 2016). This subtlety could be explored further with a more sophisticated radiative transfer model, though a diverse range of structures would have to be considered, as the clumps here are not specifically designed to represent individual starless/prestellar/protostellar cores.
Projection effects can also reduce the dust temperature: a clump always has a higher actual distance than its projected distance from the OB star, reducing the 1/d 2 term in the OB heating (Equation 3). Projection effects are most signif- Table 2. OB spectral types, power absorbed per unit mass of dust from ISRF, power absorbed per unit mass of dust from OB star without extinction, power absorbed at 1 pc with extinction Av = 0.5, least squares fit power at 1 pc and heating range for dust temperatures above 18 K. See Sect. 3 for details of models and fit.

Type
Γ ISRF Γ OB Γ OB (Av = 0.5) Γ OB (fit) range W / kg W / kg W / kg W / kg pc O9.5 -08 2.33 × 10 −2 5.2 -9.2 (1.2 -1.9) × 10 −1 2.5 ± 0.3 × 10 −1 2.4 B4 -B0 2.33 × 10 −2 (4.3 -450) × 10 −2 (3.0 -100) × 10 −3 4.2 ± 0.7 × 10 −3 0.4 A0 -B5 2.33 × 10 −2 (0.9 -31) × 10 −3 (1.7 -24) × 10 −4 8.5 ± 2.4 × 10 −4 0.2 icant for clumps at small projected distance but large lineof-sight distance from the OB star, such as the front/back of shells, and are unimportant if the cloud is filamentary in the plane-of-the-sky and clumps are offset along the filament. The overall effect is to increase the scatter in temperatures at every distance, and particularly at small projected distance. As cloud geometry is varied and uncertain, and to estimate even an average correction to the temperatures requires assumptions about the distribution of clumps, in our modelling, we take account of the projected distance simply by excluding a small fraction of clumps with particularly low temperatures for a given distance (see below). Lastly, the dust temperatures produced by the model with combined OB star and ISRF heating were fitted to the measured dust temperatures as a function of distance (as shown in Figure 4) to determine a characteristic value of the OB heating parameter ΓOB for each of the spectral type categories. The fitting was carried out using a nonlinear least-χ 2 fit (MPFIT 2 ). The fit results for the three stellar spectral type ranges are given in Table 2 and overplotted on Figure 4. We excluded clumps within a JCMT beam radius of the OB star from the fit, because their fluxes are affected by freefree emission and discs, if present, and iteratively excluded clumps with temperatures well below the model to avoid fitting those clumps that likely have significantly greater actual distances than their projected distance (i.e., those that were > 3σ below the fit after 10 iterations), ultimately excluding 9% of clumps.

DISCUSSION AND CONCLUSIONS
The amount of heating and range of the heating varies significantly with the spectral type of the major OB star in the region. Taking the range of influence of the stars as the distance at which the temperature drops below 18 K, the O-type stars can be seen to have an impact over more than 2 parsecs, essentially the entire region in which they appear, whereas the B stars' influence is less than 0.5 pc (0.4 pc for early and 0.2 pc for late-type stars). The impact of O stars on dust heating is large, with clump dust temperatures raised to 40 K close to the OB star and above 20 K to the furthest extent that clumps are detected in the map (for example, the 2.5 pc transverse extent of the W40 and Orion B NGC2024 maps). Although the absolute values of these temperatures and thresholds will vary if dust opacities vary (where we have assumed a fixed β) the trend is robust. The impact of the heating is also apparent as bright emission in the midinfrared (e.g. Herschel 70 µm in Aquila; Bontemps et al. 2010 or Spitzer 24µm in Orion B; Megeath et al. 2012 ). Midinfrared emission can be produced by small masses of hot dust (above around 50 K) whereas the submillimetre-based temperatures suggest that the heating also penetrates deeper into the clumps to at least Av = 0.5. The O-star heating range is larger than the radio Hii region extent (6 ∼ 0.8 pc for W40 Rodney & Reipurth 2008;Vallee & MacLeod 1991, 0.98 pc for NGC2024 Gordon 1969, which indicates that the penetration of the O-star radiation varies with direction, consistent with a clumpy medium, and supported by the larger size of the associated infrared nebulosity than the radio Hii region. By contrast, the heating range of B-type stars is typically less than half a parsec. At the luminous end of the range, the power absorbed by dust at 1 pc from a B0 star (with extinction of Av = 0.5) is increased by a factor of 5.3 over that from the interstellar radiation field and results in a factor 1.3 increase in the dust temperature (as Γ cool ∝ T 4+β d ), raising temperatures from 15 K to 20 K (see red dashed line on Fig. 4)), but this reduces rapidly with spectral type. For late B types, the heating range of 0.2 pc is comparable with the typical 0.1 pc size of a protostellar core, suggesting that any radiative feedback will be limited to only to the immediate environment. This range is consistent with the location of heated clumps observed in Serpens MWC 297 (Rumble et al. 2015), Perseus NGC1333 (Hatchell et al. 2005, Chen et al. 2016) and L 1495 , for example.
Splitting the sample into 179 'heated' and 564 'unheated' clumps (those that lie within the heating ranges listed in Table 2), the 'heated' clumps have a mean temperature of 24.0 K, significantly larger than 15.5 K for the 'unheated' clumps.
Is stellar radiative heating influencing the subsequent star formation? We expect the main effect of heating on core stability is to increase the thermal support in the inner, densest regions of cores (n > 10 4.5 cm −3 ; Goldsmith 2001) where the dust and gas are well coupled. Models of externally-heated cores show that gas temperatures interior to the core rise by a similar factor due to FIR emission from small grains penetrating the interior where the dust and gas temperatures are well coupled (Galli et al. 2002). This coupling is supported by temperatures measured in the dense gas by the Green Bank Ammonia Survey NH3, which rise above 20 K near the B stars in NGC1333 (Perseus) and L1688 (Ophiuchus) (Friesen et al. 2017). Although the dust temperatures are generally higher than the gas temperatures (as expected from theory and supporting our choice of β), the difference falls to less than 2 K towards the dense pre/protostellar cores, with gas and dust kinetic temperatures both approaching 20 K in the dense cores within 0.4 pc of the B2II star S1 in Ophiuchus (in NGC1333, there are no dense cores within heating range of the B5 star). Temperatures of around 20 K for these cores are also found from the joint Herschel-SCUBA-2 analysis of Ophiuchus which solves for T and β independently (Howard et al. 2021). In both Ophiuchus and Perseus, gas and dust temperatures rise towards the heated edge of the cloud, with the dust temperature rising more rapidly than the gas temperature as opacities drop and dust cooling reduces. With clump temperatures enhanced by 55% on average, depending on distance from the star, then the thermal support term in the virial equation is increased proportionally. As most cores lie within a factor two of virial equilibrium, this extra heating could stabilise cores in the 'heated' clumps against gravitational collapse.
Radiative feedback can also have a substantial effect on the stellar initial mass function (IMF), suppressing fragmentation and reducing the number of low-mass stars (Bate 2009;Krumholz et al. 2011;Mathew & Federrath 2020). This impact can result from even moderate temperature rises as the critical mass (Jeans mass MJ or Bonnor-Ebert mass MBE) is proportional to T 3/2 . Namely, the 55% mean increase in temperature from 15 K to 24 K that we see here (typical of B0 stars at 0.5 pc) raises the critical mass by a factor of two.
A secondary effect of OB stellar heating can be to raise the external pressure in the gas surrounding the clumps. This pressure is not well traced by the dust temperature partly because the dust and gas temperatures are decoupled in these low density regions but mostly because the main contribution to the pressure comes from turbulence. External pressure counteracts thermal pressure and acts alongside selfgravity to bind cores. From a virial analysis in Ophiuchus using linewidths of C 18 O as a tracer, Pattle et al. (2015) found that external pressure is responsible for binding several starless cores, particularly those around the B star S1. Indeed, the Green Bank Ammonia Survey is finding similar results in other regions Kerr et al. (2019). Further studies are needed to understand if the increase in external pressure counteracts the increase in internal support in the heated regions and mitigates the consequences of stellar heating for the star formation efficiency and IMF.
OB stars can produce much more dramatic changes to star formation regions than simply heating dust. For example, molecular material can be ionised, photodissociated, collected, disrupted, shocked and moved around. The overall positive and negative effects of OB star feedback are hard to disentangle (Dale et al. 2015). Nonetheless, given that we find that 24% of clumps are heated externally, even in the nearer star-forming regions with few massive stars, simulations which take into account the effects of radiative feedback, and not just ionising UV radiation, from nearby highand intermediate-mass stars are needed.