CHARACTERIZING COSMIC-RAY PROPAGATION IN MASSIVE STAR-FORMING REGIONS: THE CASE OF 30 DORADUS AND THE LARGE MAGELLANIC CLOUD

Spitzer imaging of the LMC was carried out as part of the Surveying the Agents of a Galaxy's Evolution (SAGE; Meixner et al. 2006) legacy program. The observations covered an area of ∼7° × 7° over the LMC in all IRAC (3.6, 4.5, 5.8, and 8 μm) and MIPS (24, 70, and 160 μm) bandpasses. The individual observations were calibrated, combined, and mosaicked by the SAGE team, using a custom pipeline (see Meixner et al. 2006 for further details). Because we are only interested in diffuse emission, we use point-source-subtracted versions of these images prepared by the SAGE team (see right panel of Figure 1). For the present study, we only make use of the MIPS 24, 70, and 160 μm data, which have native resolutions of ≈57, 17'', and 38'', respectively. The photometric uncertainties are conservatively taken to be 5%, 10%, and 15% at 24, 70, and 160 μm, respectively (see the MIPS Instrument Handbook⁷).

Zoom In Zoom Out Reset image size

The majority of the emission at 24 μm arises from dust heated in the vicinity of massive stars, which are the progenitors of core-collapse supernovae. Consequently, the 24 μm data are used as a proxy for the CR source distribution in our phenomenological modeling. We use the full infrared data (i.e., 24, 70, and 160 μm imaging) to estimate the radiation field energy density over the LMC and 30 Dor (see Section 2.4).

We use combined 1.4 GHz single-dish (Parkes 64 m) and interferometric (ATCA) data presented in Hughes et al. (2006), where a much more detailed description of the data preparation can be found. We employ a version of the radio map that has point sources removed (see left panel of Figure 1) from the map using a median filter technique (see Hughes et al. 2006) because more than 90% of such sources are expected to be background active galactic nuclei (AGNs; Marx et al. 1997). The angular resolution of the final combined mosaic is 40'', and sets the resolution for our analysis between the radio and infrared data. At our assumed distance for the LMC, this projects to a linear scale of ≈10 pc.

Because we are interested in the non-thermal component of the radio data, we subtract the contribution from thermal (free–free) emission using a scaled version of the 24 μm image following Equation (2) of Murphy et al. (2006b). This relation is based on the empirical correlation found between the 24 μm and extinction-corrected Paα luminosities from star-forming regions within NGC 5194 by Calzetti et al. (2005). Although not universal (e.g., Pérez-González et al. 2006; Calzetti et al. 2007), it has proved to be a good first-order estimate of the thermal radio emission in nearby galaxies when compared to thermal fractions derived from multifrequency radio data (e.g., Murphy et al. 2008).

To check the reliability of our free–free emission estimate, we compare it with the thermal fraction derived from standard radio spectral decompositions using single-dish radio data from the literature: 1.4 GHz (Klein et al. 1989), 2.3 GHz (Mountfort et al. 1987), and 2.45 GHz (Haynes et al. 1991). There is a known discrepancy between the total flux from the 1.4 GHz single-dish map and the combined Parkes+ATCA map being used here, as Hughes et al. (2007) report a global 1.4 GHz flux density that is 1.3 times smaller than the single-dish estimate from Klein et al. (1989). This difference is largely attributed to an underestimation of the beam width for the single-dish data. Because the single-dish data were reduced and analyzed self-consistently by Haynes et al. (1991), we assume that each map is affected in a similar way.

Using these multifrequency radio data, and assuming a constant non-thermal spectral index of α_NT ≈ 0.7 (Haynes et al. 1991), we derive a thermal fraction at 1 GHz that is ≈1.2 times larger than the value of 25% found using the 24 μm map over our region of interest (RoI) surrounding 30 Dor shown in Figure 1. This difference is consistent with results from a high-resolution radio continuum survey of M 33 (Tabatabaei et al. 2007), where a standard radio spectral decomposition was found to yield thermal fractions for 11 H  ii complexes that were overestimated by an average factor of 1.5 at 8.3 GHz. However, even if we scale the 24 μm derived free–free map by a factor of 1.2 to match the total thermal fraction measured by the radio spectral decomposition in our RoI, our results are not significantly affected (see Section 3.2).

2.2.1. CR Electron Energy Estimate

We estimate the typical energies of the CR electrons emitting the observed non-thermal 1.4 GHz emission. For electrons propagating with a pitch angle θ in a magnetic field of strength B with isotropically distributed velocities, such that 〈sin ²θ〉 = 2/3 leading to B_⊥ ≈ 0.82B, then a CR electron emitting at a critical frequency ν_c will have an energy

$$\begin{equation} \left(\frac{E_e}{\rm GeV}\right) = 8.8 \left(\frac{\nu _c}{\rm \,GHz}\right)^{1/2} \left(\frac{B}{\rm \mu G}\right)^{-1/2}. \end{equation} \tag{ 1 }$$

Integrating the free–free-corrected 1.4 GHz flux over the region shown in Figure 1 (S^NT_1.4 GHz = 102.26 Jy; see Section 2.4), we use the revised minimum energy calculation of Beck & Krause (2005), assuming a proton-to-electron number density ratio of K₀ ≈ 100  ±  50, a non-thermal spectral index of α_NT ≈ 0.7  ±  0.1, and a path length of l ≈ 1  ±  0.5 kpc, to calculate a minimum energy magnetic field strength of ∼11 ± 4 μG. Using this range of values in Equation (1), our 1.4 GHz maps are most sensitive to ∼3 ± 0.5 GeV CR electrons. If the minimum energy condition does not hold, and we instead assume an extreme range in B, such as 3–50 μG, the corresponding CR electron energies will range from 6 to 1.5 GeV. However, the (in)validity of the minimum energy assumption is uncertain and often assumed (as we do in this paper). We briefly discuss the impact of this on the derived properties of the diffusion coefficient in Section 4.2.1.

The Fermi-LAT instrument, event reconstruction, and response are described in Atwood et al. (2009). In this paper, we use events and instrument response functions (IRFs) for the standard low-background "Clean" events corresponding to the Pass 7 event selections.⁸ To minimize the contribution from the very bright Earth limb, we restrict the event selection and exposure calculation to zenith angles <100°. We selected all front-converting events, for which the point-spread function (PSF) is narrowest, within a 20° square RoI centered on R.A. 05^h38^m42^s, decl. −69°06'03'' (J2000) for ≈32 months of sky survey data from 2008 August 4 until 2011 March 24. Exposure maps and the PSF for the pointing history of the observations were generated using the standard Fermi-LAT ScienceTools package available from the Fermi Science Support Center.⁹ The exposure of the instrument over the RoI for the data taking period used in this analysis is very uniform.

We constructed foreground-subtracted γ-ray maps over the same region using the standard Fermi-LAT diffuse emission model¹⁰ together with the appropriate isotropic background model and γ-ray point sources from the Fermi-LAT Second Source Catalog (Nolan et al. 2012), forward folding the combination with the exposure and PSF to obtain the total counts expected for these foregrounds. These were subtracted from the data and the residual emission for the 1–3 GeV energy interval is shown in Figure 1 overlaid on the free–free-corrected 1.4 GHz and 24 μm maps. The effective 68% containment radius for the Fermi-LAT PSF over this energy bin for the event selection used is ≈050 (≈440 pc), assuming an E^−2.1_γ spectrum for the γ-ray emission, which is a sufficient approximation to the spectral shape determined by Abdo et al. (2010b) for the LMC γ-ray emission over this energy range.¹¹ Note that the effective FWHM of the Fermi-LAT PSF over this energy bin is ≈038 (≈330 pc), significantly smaller than the 68% containment radius quoted above.

Similarly, we constructed a map for the 3–10 GeV residual emission. This is shown in Figure 2 where we overlay the 1–3 and 3–10 GeV contours on the H i column density data (Staveley-Smith et al. 2003; Kim et al. 2003). The effective 68% containment radius for this map is ≈024 (≈210 pc), with an effective FWHM of ≈019 (≈170 pc), again assuming an E^−2.1_γ spectrum for the γ-ray emission over this energy bin. However, the 3–10 GeV map has much lower statistics, even near 30 Dor, contrasting with the 1–3 GeV γ-ray map, which has ≳ 3 times more events. Therefore, we base the majority of our spatial analysis on the 1–3 GeV map, and only use the higher energy γ-ray data for testing the energy dependence of the CR propagation (see Section 4.2.2).

Zoom In Zoom Out Reset image size

2.3.1. CR Nuclei Energy Estimate

As for electrons and the 1.4 GHz emission, we estimate the typical energies of the CR nuclei emitting γ rays for the energy ranges that we consider in this paper. For non-AGN-dominated star-forming galaxies like the Milky Way and the LMC, the assumption that the CR nuclei are the predominant relativistic particle population producing the γ-ray emission in the Fermi-LAT energy range is reasonable (Strong et al. 2010).

We define the function

$$\begin{eqnarray} G(E_\gamma ,E_p) &=& \left(\int _{E_p}^\infty dE_p^{\prime } \, \beta \frac{dJ_p}{dE_p^{\prime }} \frac{d\sigma _{pp}(E_\gamma ,E_p^{\prime })}{dE_\gamma }\right) \nonumber \\ &&\bigg / \left(\! \int _{E_p^{\rm min}(E_\gamma)}^\infty dE_p^{\prime } \, \beta \frac{dJ_p}{dE_p^{\prime }} \frac{d\sigma _{pp}(E_\gamma ,E_p^{\prime })}{dE_\gamma } \!\right) - \frac{1}{2}, \nonumber\\ \end{eqnarray} \tag{ 2 }$$

where E_p and E_γ are the proton and γ-ray energies, β = v/c, dJ_p/dE'_p∝E'^−α_p is the CR proton flux, dσ_pp/dE_γ is the differential cross section for γ-ray production calculated as in Moskalenko & Strong (1998), and E^min_p is the minimum proton energy required to produce a photon of energy E_γ. The effective proton energy, E^eff_p, for a given E_γ is given by the root of G(E_γ, E_p) so that contributions to the γ-ray flux at E_γ below and above E^eff_p are equal. The effective proton energy is a function of the proton spectral index E^eff_p = E_p^eff(α), which we show for different α for the 1–3 GeV and 3–10 GeV energy ranges in Table 1. We also show the effective (average) γ-ray energy for the 1–3 GeV and 3–10 GeV energy ranges, calculated assuming the same spectral index for the γ rays as for the protons. This assumption is valid for thin-target interactions by the CR nuclei with gas, which applies for the ISM surrounding 30 Dor and the greater LMC.

Table 1. Effective γ-ray and Proton Energies for Different γ-ray Energy Ranges

Proton Energy	γ-ray Energy Range
Index	1–3 GeV		3–10 GeV
	Eeffγ	Eeffp	Eeffγ	Eeffp
	(GeV)	(GeV)	(GeV)	(GeV)
2.1	1.63	22.4	5.10	69.2
2.2	1.61	19.5	5.04	60.4
2.3	1.60	17.2	4.99	53.6
2.4	1.58	15.4	4.93	48.1
2.5	1.60	14.1	4.87	43.6
2.6	1.55	12.5	4.82	39.8
2.7	1.54	11.5	4.77	36.6
2.8	1.53	10.5	4.72	33.9

Download table as:  ASCIITypeset image

The effective proton energy changes by a factor of ∼2 when the power-law index of the proton spectrum changes from 2.1 to 2.8, while the effective energy of γ rays does not change appreciably. For the γ-ray energies, this is not surprising because the effective energy is calculated over a relatively narrow energy bin. The effective proton energy change is a combination of the fairly steep rise in the inclusive pion production cross section as the momentum of the interacting proton increases, approximately ∝E^0.5 between 10 and 100 GeV (see, e.g., Figure 2(a) from Dermer 1986). When the proton spectrum is steep (e.g., with an index ∼2.8), this rise is not large enough to compensate for the steeply falling number of high-energy protons. When the proton spectrum is flat (e.g., with an index ∼2.1) the rise in the inclusive cross section increases the contribution of high-energy protons, making the E^eff_p significantly larger.

All images were cropped to a common field-of-view and regridded to a common pixel scale. To properly compare each image at the same resolution, maps were convolved with an appropriate PSF. For the radio and infrared (24, 70, and 160 μm) data, we convolve all images to the resolution of the 1.4 GHz image, which has an FWHM of 40'', using a Gaussian beam having an FWHM equal to the quadrature difference between the final and original FWHM values of each image.

To match the resolution of the 24 μm data to that of the γ-ray maps, we convolve the 24 μm image (at its native resolution) with the in-flight Fermi-LAT PSF, described in Section 2.3. Because the FWHM of the 24 μm PSF is <1% of the effective FWHM over the γ-ray energy bins used in our analysis, we do not attempt to correct for this additional broadening.

Next, we apply our image smearing analysis (described below) to the resolution-matched images to determine the extension of the non-thermal radio and γ-ray morphologies relative to the 24 μm morphology. The procedure largely follows that described in Murphy et al. (2006b, 2008), where we refer the reader for a more detailed description. Our underlying assumption is that the only difference between these distributions is due to the diffusion and energy losses of the CRs because of the common origin for the CR leptons and nuclei. We focus on the 24 μm data in the smearing analysis rather than the 70 μm data, as was done by, e.g., Murphy et al. (2008), because it traces warmer dust and has been found to be more peaked near star-forming regions than the 70 μm morphology (Helou et al. 2004; Murphy et al. 2006a). For instance, typical 24 μm knots surrounding 30 Dor tend to be ∼10 pc in size. Thus, the 24 μm maps likely act as a better source function for recently accelerated CRs.

To summarize the procedure, we convolve the entire 24 μm image of the LMC by a parameterized kernel κ(r) and compute the residuals between the free–free-corrected 1.4 GHz and smoothed 24 μm maps. The smoothing kernel is a function of a two-dimensional position vector r, having a magnitude of r = (x² + y²)^1/2, where x and y are the right ascension and declination offsets on the sky, respectively. We investigate both exponential and Gaussian kernels because a preference of kernel type may suggest different CR transport effects. The exponential and Gaussian kernels have the forms of $\kappa (\textbf {\textit {r}})= e^{-\textbf {\textit {r}}/l}$ and $\kappa (\textbf {\textit {r}})= e^{-\textbf {\textit {r}}^{2}/l^{2}}$, respectively, where l is the e-folding scale length. Thus, for the Gaussian kernel, σ² = l²/2.

Gaussian kernels suggest a simple random walk diffusion scenario. On the other hand, exponential kernels can arise as the result of energy losses or escape on timescales comparable to diffusion, modifying the Gaussian kernels (e.g. Bicay & Helou 1990; Helou & Bicay 1993). Because the LMC lacks a well-defined disk, and has a low inclination (i ≈ 31°; Nikolaev et al. 2004), we assume spherical symmetry for 30 Dor. Hence, there is no need to modify the position vector by geometric factors.

We calculate the normalized-squared residuals between the radio and the infrared images

$$\begin{equation} \phi (l) = \frac{\Sigma [Q^{-1}\tilde{I}_{j}(l) - R_{j}]^2}{\Sigma R_{j}^{2}}, \end{equation} \tag{ 3 }$$

where R is the radio image with free–free emission removed, $\tilde{I}(l)$ is the infrared image smoothed by a kernel with scale length l, Q = ΣI_j/ΣR_j is used as a normalization factor (i.e., the global 24 μm/1.4 GHz ratio), and the subscript j indexes each pixel. Only pixels detected with a significance of >3σ of the rms noise in each map are used in the calculation of the residuals. The minimum in ϕ defines the best-fit scale length, which is taken to be the typical distance traveled by the CR electrons. Similarly, we perform the same procedure using the 24 μm and 1–3 GeV γ-ray maps to determine the typical distance traveled by the CR nuclei.

For this procedure, generically, if there is no broadening due to propagation effects, the minimum of Equation (3) will be at l = 0, because both the radio and γ-ray data have already been reduced to the same angular resolution. Below, we obtain kernel scale lengths indicating that significant smoothing of the 24 μm image is required to improve the match with the radio and γ-ray images. While the scale lengths obtained for the radio/24 μm residual analysis are much larger than the angular resolution of these maps, for the γ-ray/24 μm data, the derived scale lengths are comparable to the angular resolution after convolution with the Fermi-LAT PSF. The technique is capable of detecting scale lengths that are a small fraction of the resolution (i.e., ≪Δl = 50 pc, which is the scale length step size used in the analysis), and is more sensitive for higher signal-to-noise ratio maps. We emphasize that the kernel scale lengths obtained in this case are in addition to the smoothing of the 24 μm data to match the Fermi-LAT PSF. If there was no additional effect from propagation, the derived scale lengths would be zero, whereas we detect a meaningfully non-zero scale length.

The uncertainty in ϕ is estimated by numerically propagating the uncertainties in the input images as measured by the 1σ rms noise of each map, with the uncertainty for the best-fit scale length then estimated as the range in scale length corresponding to that from min (ϕ) to min (ϕ) + unc(ϕ) along the residual curve. We note that this places a lower limit on the uncertainty for the best-fit scale lengths since there is additional uncertainty on the 1.4 GHz thermal fraction estimation and foreground subtraction for the γ-ray maps. The calculation of the residuals and photometry was carried out within an aperture having a radius of 1° (≈875 pc), encompassing 30 Dor, centered at 05^h39^m07^s, −69°22'02'' (J2000, see Figure 1).

We use the 24, 70, and 160 μm photometry to calculate the total infrared (IR, 8–1000 μm) luminosity over this region by fitting these data to the spectral energy distribution (SED) models of Dale & Helou (2002) and integrating the best-fit SED between 8 and 1000 μm. The individual 24, 70, and 160 μm flux densities are 0.42 ± 0.02, 4.1 ± 0.41, and 8.0 ± 1.7 × 10⁴ Jy, respectively. In addition to photometric uncertainties, the above errors include a term for the mean rms noise values in the convolved 24, 70, and 160 μm maps, which are measured to be ∼0.045, 0.685, and 13.2 MJy sr⁻¹, respectively. We obtain an IR luminosity of L_IR = 1.15 ± 0.12 × 10⁴² erg s⁻¹ (3.00 ± 0.31 × 10⁸ L_☉). Taking this value, we estimate the corresponding radiation field energy density

$$\begin{equation} U_{\rm rad} \approx \frac{2\pi }{c}I_{\rm bol} \gtrsim \frac{L_{\rm IR}}{2 A_{\rm IR}c}\left(1 + \sqrt{\frac{3.8\times 10^{42}}{L_{\rm IR}}}\right), \end{equation} \tag{ 4 }$$

where I_bol is the bolometric surface brightness, c is the speed of light, and A_IR ≈ 2.4 kpc² is the area over which the photometry was measured. All quantities in Equation (4) are in cgs units. This calculation is for radiation emitted near the surface of a semitransparent body, and the parenthetical term provides an empirically derived correction for non-absorbed UV emission (Bell 2003), resulting in a value of U_rad ≈ 2.39 ± 0.17 × 10⁻¹² erg cm⁻³ (1.49 ± 0.11 eV cm⁻³) averaged over the volume considered.

Figure 1 shows the residual 1–3 GeV γ-ray emission overplotted as contours on the free–free-corrected 1.4 GHz radio and 24 μm warm dust emission in the left and right panels, respectively. While there is a general correspondence between the γ-ray, radio, and infrared emission, the peak of the γ-ray emission is clearly offset from the peak of the radio and infrared emission, which appear to be nearly co-spatial. This lack of correspondence between the peak of the γ-ray emission with the radio and infrared emission leads one to question whether the observed γ-ray emission is indeed associated with 30 Dor.

In Abdo et al. (2010b), it was noted that γ-ray emission near 30 Dor may have a non-negligible contribution from two Crab-like pulsars, PSR J0540−6919 (05^h40^m112, − 69°19'54''; J2000) and PSR J0537−6910 (05^h37^m474, − 69°10'20''; J2000). In addition, there is a background X-ray source (RX J0536.9−6913: 05^h36^m578, − 69°13'26''; J2000) that is also coincident with the line of sight toward 30 Dor, which could also contribute to the observed emission. The locations of these point sources are shown in Figures 1 and 2, and are near, but not at, the peak of the 1–3 and 3–10 GeV γ-ray emission. It is also interesting to see that the 1–3 and 3–10 γ-ray contours appear to peak close to a region of high H i column density in Figure 2. It is possible to subtract fitted γ-ray point sources at the locations of the pulsars and AGN, resulting in a situation where 30 Dor does not exhibit γ-ray emission in excess of the subtracted foreground model. The implications of such a scenario are described in Section 4.1. However, we note that there is no detection of pulsed γ-ray emission so far reported from either of the pulsars. Also, no variability is detected from the region during our period of analysis. If RX J0536.9−6913 is a γ-ray-emitting AGN, the lack of variability makes it difficult to determine if this candidate source contributes significantly in our RoI.

In Figure 3, we plot the residuals (ϕ, see Section 2.4) between the free–free-corrected 1.4 GHz and smoothed 24 μm maps as a function of exponential and Gaussian kernel scale lengths. The residuals between the free–free-corrected 1.4 GHz and 24 μm maps are decreased by more than a factor of ∼2 after smoothing the 24 μm map using either exponential or Gaussian kernels. A slightly larger improvement is found by using an exponential kernel, suggesting energy losses and/or escape may be important for the CR electrons on timescales comparable to the diffusion timescale. The corresponding best-fit exponential and Gaussian kernel scale lengths are 100⁺¹⁰₋₁₀ and 200⁺¹⁰⁰₋₅₅ pc, respectively. We note that the estimated ∼3 GeV CR electron propagation length reported here assuming random walk diffusion (i.e., using a Gaussian kernel) is consistent within errors to the GeV CR proton confinement length reported by Abdo et al. (2010b), whose estimate assumed Gaussian profiles (see Figure 3). Abdo et al. (2010b) report the Gaussian σ to be σ = 170 ± 60 pc, so we multiply this number by $\sqrt{2}$ for proper comparison with the scale length definition for our Gaussian kernel.

Zoom In Zoom Out Reset image size

As discussed in Section 2.2, we test how our results are affected by assuming that the thermal fraction of the 1.4 GHz radio map is being underestimated by the 24 μm maps relative to what is derived using single-dish radio data at 1.4, 2.3, and 2.45 GHz. By scaling the 24 μm derived free–free maps by a factor of 1.2 to match the total thermal fraction measured by the radio spectral decomposition, and repeating the smearing analysis, we find that the best-fit exponential scale length is still 100 pc, while the best-fit Gaussian scale length is slightly increased to 250 pc. A larger improvement is again found using an exponential kernel relative to a Gaussian kernel. The results therefore do not appear to be significantly affected by increasing the thermal fraction estimate to match that from the radio spectral decomposition method.

Similarly, in Figure 4, we plot the residuals between the 1–3 GeV γ-ray and smoothed 24 μm maps as a function of exponential and Gaussian kernel scale lengths. The residuals between the 1–3 GeV and 24 μm maps are decreased by a factor of ∼2 after smoothing the 24 μm map using either exponential or Gaussian kernels. While a slightly larger improvement is found by using exponential kernels in the comparison between the 1.4 GHz and 24 μm maps, we do not find a preference in kernel type for the 1–3 GeV γ-ray map. The corresponding best-fit exponential and Gaussian kernel scale lengths are 200⁺⁴⁰₋₁₀ and 450⁺⁶⁵₋₅₀ pc, respectively. These best-fit scale lengths are significantly larger (i.e., a factor of ∼2) than what was measured for the 1.4 GHz maps. Because the 1.4 GHz and 1–3 GeV γ-ray maps probe different energy CR particle populations, the corresponding differences in their best-fit scale lengths may arise from different diffusion speeds. This scenario is discussed in Section 4.2.

Zoom In Zoom Out Reset image size

As stated in Section 3.1, there are both pulsars and a background X-ray-emitting source that may contribute to the observed γ-ray emission toward 30 Dor. For this case, the situation has 30 Dor emitting at 1.4 GHz, therefore containing diffusing CR electrons, but no observational signature for the CR nuclei. Explaining this would require the CR nuclei to escape the system without interacting with the interstellar gas.

Taking the above estimate for the bolometric luminosity from 30 Dor (i.e., the IR luminosity given in Section 2.4 corrected by the parenthetical term in Equation (4) for non-absorbed UV emission) together with the updated SFR calibrations given in Murphy et al. (2011), we estimate a corresponding SFR of ∼0.15 M_☉ yr⁻¹. This value is consistent with others in the literature (e.g., Hughes et al. 2007; Harris & Zaritsky 2009; Lawton et al. 2010). Assuming this value and the sensitivity of the Fermi-LAT data used in our analysis, the absence of γ-ray emission from 30 Dor clearly disagrees with the empirical correlations between SFR or the product of supernova rate and the gas mass with γ-ray luminosity found for the local group and nearby starbursts (e.g., Abdo et al. 2010a; Lenain & Walter 2011). Furthermore, theoretical expectations of such scaling relations and their implications for the contribution by normal star-forming galaxies to the diffuse γ-ray extragalactic background (e.g., Pavlidou & Fields 2001; Fields et al. 2010) would then be questionable. However, because there is currently no strong evidence linking the γ-ray emission from 30 Dor to the known point sources in the field, we discuss the physical interpretation of our analysis with the assumption that it arises from 30 Dor.

To interpret our findings on the kernel scale lengths, we must consider the dependence of gas density on the synchrotron and γ-ray emission. The diffuse γ rays in the energy range that we have used (1–3 GeV) are predominantly from CR nuclei interacting with the interstellar gas, and therefore S_γ∝n_pn_ISM where n_p and n_ISM are the densities of the CR nuclei and gas, respectively. The diffuse synchrotron emission is proportional to the CR electron density, n_e, and the square of the magnetic field strength, B, i.e., S_sync∝n_eB². Assuming flux-freezing scaling of $B \propto \sqrt{n_{\rm ISM}}$ (e.g., Ruzmaikin et al. 1988; Niklas & Beck 1997; Crutcher 1999), the synchrotron emission becomes proportional to the product of the electron density and the density of the interstellar gas, S_sync∝n_en_ISM. This simple scaling argument therefore suggests that the differences in appearance of the γ-ray and radio images are due to differences in the distributions of the CR nuclei and electrons. This was tested by normalizing the γ-ray and free–free-corrected radio images by the H i column density map before applying the smearing analysis, recovering results consistent with those described in Section 3.2.

4.2.1. Cosmic-ray Transport Properties

We assume that the transport of the CRs can be described by a simple random walk process with a rigidity-dependent spatial diffusion coefficient D_R (e.g., Ginzburg et al. 1980). The diffusion coefficient is defined as

$$\begin{equation} D_{R} = D_{0} \left(\frac{R}{\rm GV}\right)^{\delta } \approx \frac{l^{2}_{\rm diff}}{\tau _{\rm diff}}, \end{equation} \tag{ 5 }$$

where D₀ is the normalization constant and l_diff is the characteristic distance that CRs travel after a time τ_diff. For our Gaussian kernels, we can relate l_diff to the corresponding best-fit scale lengths such that l²_diff = σ² = l²/2.

If we assume that the CR electrons and nuclei are injected by the same source(s) and have been propagating through the ISM for the same length of time, we find that the ∼20 GeV CR nuclei diffusion coefficient is ∼4 times larger than that for the ∼3 GeV CR electrons. Solving Equation (6) for δ, we find that the diffusion coefficient scales as (R/GV)^{0.69 ± 0.15} from the exponential best-fit scale lengths and as (R/GV)^{0.81 ± 0.30} using the Gaussian best-fit scale lengths. Errors on δ were estimated by a standard Monte Carlo approach using the uncertainties in the best-fit scale lengths.¹² If the minimum energy assumption does not hold, and we instead consider the extreme range of B ≈ 3–50 μG described earlier in Section 2.2.1, we find corresponding δ values of ≈0.5–1.2.

These model-independent values of δ are consistent with the value of ∼0.6–0.7 used in empirical diffusion models to fit the observed secondary-to-primary ratios, typically boron-to-carbon (B/C), for the Milky Way. In addition, our derived values for δ generally exceed those for physically motivated turbulence theories, being larger than the nominal value of 1/3 for Kolomogorov (Kolmogorov 1941), and marginally consistent with the value of 1/2 for Iroshnikov–Kraichnan (Iroshnikov 1964; Kraichnan 1965) turbulence spectra. However, the values for the Milky Way are found for a large volume in a spiral galaxy, while the results of our analysis are for a single, highly active star-forming region: 30 Dor exhibits a complex network of kinematic features including slow (v ≲ 100 km s⁻¹) and fast (100 km s⁻¹ ≲ v ≲ 300 km s⁻¹) expanding shells powered by stellar winds from young massive stars and supernovae (Chu & Kennicutt 1994).

Because of strong radiative energy losses as they propagate, the corresponding diffusion lengths for the ∼3 GeV CR electrons may be underestimated by our best-fit scale lengths, which would result in an overestimate for δ. This is suggested by the fact that an exponential kernel resulted in a lower minimum residual between the 1.4 GHz and 24 μm maps than a Gaussian kernel. However, as discussed below, the additional cooling of electrons may not significantly affect our estimate for δ.

Taking the average derived value for δ (i.e., 0.75), together with the age of the current star formation activity responsible for supernova remnants accelerating the CRs, we can estimate the diffusion coefficient normalization factor D₀. The 30 Dor complex is known to contain many non-coeval stellar populations ranging in age from <1 to 10 Myr from Hubble Space Telescope spectroscopy (Walborn & Blades 1997). Recently, using a Bayesian analysis, Martínez-Galarza et al. (2011) fit Spitzer Infrared Spectrograph data for 30 Dor using mid-infrared SED models (continuum + lines), arriving at a luminosity-weighted age for the system of ∼3 Myr. This suggests that CRs associated with this star-forming event have been accelerated in the supernova remnants of very massive O-stars, and that the bulk of CRs have yet to be accelerated by the supernova remnants from the more numerous and less massive (i.e., ∼8 M_☉) stars with lifetimes of ∼30 Myr. Note that this age is less than the estimated cooling lifetime τ_cool ∼ 12 ± 2.3 Myr for the 1.4 GHz emitting CR electrons as they propagate through the ISM of 30 Dor.¹³

Using this mean age of 3 Myr in Equation (5), together with the best-fit Gaussian scale lengths corresponding to l_diff ≈ 140 pc and 320 pc for the 3 GeV and 20 GeV electrons and nuclei, respectively, we find a diffusion coefficient normalization constant of D₀ ≈ (0.9–1.0) × 10²⁷ cm² s⁻¹. This is more than an order of magnitude lower than those found in the Milky Way (e.g., ≈5 × 10²⁸ cm² s⁻¹; Trotta et al. 2011). However, analytical solutions for the perpendicular diffusion coefficient (i.e., diffusion across magnetic field lines) are found to be much smaller. For example, Shalchi et al. (2010) report a perpendicular diffusion coefficient of D_⊥ ≈ (3.0–30) × 10²⁶ cm² s⁻¹. Because 30 Dor is an active star-forming region, which likely has a strong turbulent magnetic field and lacks a large-scale regular field like the Galaxy, it may not be surprising that our derived value for the diffusion coefficient normalization factor is intermediate in the range for the perpendicular diffusion coefficients obtained by other authors.

4.2.2. Results Including 3–10 GeV Maps

4.2.3. Propagation Length versus Star Formation Activity

Having estimated the average distance traveled by CR electrons for an individual star-forming region, 30 Dor, it is interesting to see how these results compare with similar estimates of CR electron propagation distances for entire galaxies. Murphy et al. (2006b, 2008) reported a correlation between the CR electron propagation distance and galaxy surface brightness such that CR electrons are found to travel shorter distances, on average, in galaxies having higher star formation activity.

In Figure 5, we plot the results from Murphy et al. (2008, i.e., the best-fit exponential scale lengths projected in the plane of the sky versus U_rad) along with our estimates for 30 Dor. We note that various systematics and selection effects (e.g., distance, inclination) were rigorously investigated by Murphy et al. (2008). They additionally tested the differences between isotropic kernels and kernels projected in the plane of each galaxy disk, finding a preference, i.e., smaller residuals, for isotropic kernels.

Zoom In Zoom Out Reset image size

The U_rad estimate for 30 Dor was calculated from the surface brightness measured over an area of 2.4 kpc², the same region for which the residuals were estimated (see Section 2.4). For this case, the best-fit scale length was measured by comparing the free–free-corrected 1.4 GHz with smoothed 70 μm maps to allow a proper comparison with the results of Murphy et al. (2008). The best-fit exponential and Gaussian kernel scale lengths were 0 and 50 pc, respectively. Because Murphy et al. (2008) plot exponential kernel scale lengths, we plot the result from smoothing the 70 μm map of 30 Dor with a Gaussian kernel as an upper limit.

Included in Figure 5 are both spirals and star-forming irregulars. The latter have markedly small best-fit scale lengths for their values of U_rad relative to the sample of spirals. This is thought to be caused by increased escape of CR electrons in these systems (e.g., Cannon et al. 2005, 2006; Murphy et al. 2008). To see how the behavior of entire galaxies compares with the LMC as a whole, we repeat the image smearing analysis using the 70 μm maps after first projecting the LMC to the mean distance of the Murphy et al. (2008) sample irregulars, d = 3.9 Mpc to take possible resolution effects into account. We calculate the residuals and U_rad within an aperture having a radius of ≈3.1 kpc (275 at the distance of 3.9 Mpc) centered at 05^h18^m50^s, −68°50'49'' (J2000). This is approximately the area over which the 24 μm emission was detected at ≳ 3σ level. The best-fit exponential and Gaussian scale lengths are 200 and 400 pc, respectively. This places the LMC close to the main trend among the star-forming spirals.

Focusing only on the region containing the star-forming spirals, we re-plot the disk galaxies along with 30 Dor and the LMC in Figure 6. The dotted line is a least-squares fit to the Murphy et al. (2008) galaxy disks, excluding the LMC and 30 Dor. Clearly, the position of 30 Dor in the plot is consistent with the empirical trend describing spiral galaxies, but extrapolated to a surface brightness greater by almost an order of magnitude. Thus, this scaling relation appears to operate on the scales of entire galaxies all the way down to individual star-forming regions.

Zoom In Zoom Out Reset image size