SOFIA Observations of 30 Doradus. II. Magnetic Fields and Large-scale Gas Kinematics

In Figure 3, the integrated intensity maps of [C ii] and ¹²CO(2-1) are compared with the dust continuum maps from HAWC+ at 214 μm. This figure shows that the [C ii] peak intensities are closely correlated with the dust continuum peak and that it also traces the extended emission. We find that ¹²CO(2-1) is more sensitive to the gas located in the high column density regions, which can be expected for this lower metallicity region. As a result, the [C ii] emission will be more sensitive to the kinematics of the ambient gas, whereas CO is more sensitive to the dense clumps in the cloud.

Zoom In Zoom Out Reset image size

In Okada et al. (2019) and Figure 3, several [C ii] and CO spectra of 30 Dor are presented. These observations show that most of the emission is found between 220 and 280 km s⁻¹. In this velocity range, multiple velocity components and high-velocity wings are observed. In Figure 3, the wings are most clear in spectra (1)–(4) and specifically in the velocity range between 220 and 240 km s⁻¹ and between 255 and 280 km s⁻¹. Such wing-like structures were recently also found in spectrally resolved [C ii] observations toward galactic H ii regions (e.g., Schneider et al. 2018, 2020; Tiwari et al. 2021; Bonne et al. 2022a). The 30 Dor region has a complex dynamical structure. This is illustrated in more detail for [C ii] and ¹²CO(2-1) in Figures 4 and 15 (in Appendix A). Despite the complexity of the kinematics, Figure 4 shows that both lines have a broadly similar velocity structure on large scales. Inspecting the velocity structure in Figure 4, it can be noted, in particular with [C ii], that the blueshifted and redshifted gas appear to create intersecting axes toward the middle of the SOFIA map. This intersection is shown in Figure 4 by the northeast-to-southwest axis for blueshifted velocities and the east-to-west axis for redshifted velocities. This intersecting structure formed by the blue and red axis is also observed in the channel maps of the region that are presented in Figure 15 of Appendix A. The channel maps also confirm that the southern part of the maps consists of two subregions, which are discontinuous in velocity. This is consistent with the SOFIA/HAWC+ integrated maps, which show that this western region contains two clumps. Furthermore, these two subregions have a noticeably different B-field structure in Figure 9. In the northern region of the map, it is also observed that there is a noteworthy change in B-field orientation at the transition from blueshifted to redshifted gas, which confirms that the B-field morphology is directly related to the large-scale kinematics of the 30 Dor region; see Figure 4.

Zoom In Zoom Out Reset image size

From the data cubes, the moment maps can be calculated. However, before presenting these moment maps, it has to be noted that there are limitations to the moment maps due to the inherent complexity (i.e., multiple components and high-velocity wings) in the spectra. The resulting moment maps are presented in Figure 5. The first-moment map (or the velocity field) has an organized gradient that looks similar in [C ii] and ¹²CO(2-1) and fits with the velocity field obtained at smaller scales with Atacama Large Millimeter/submillimeter Array observations (Indebetouw et al. 2013). Even though both lines show a similar morphology, it is found that the [C ii] velocity field covers a significantly larger velocity range. The CO kinematics thus tend to follow the [C ii] kinematics, but less drastically.

Zoom In Zoom Out Reset image size

Both lines are expected to trace different regions in the cloud because of the critical density ¹⁷ and the fact that far-UV radiation propagates farther into the cloud in low-metallicity regions. From the second-moment maps, presented in Figure 5, it is observed that the line width significantly depends on the observed line. For this, it has to be taken into account that there are indications of multiple components, also reported in Melnick et al. (2021), which are particularly clear in [C ii]. For example, two velocity components are prominent in the most northeastern part of the map (see Figure 3), which is translated in a large second moment in that region. From the maps, it is also observed that the northern region has a higher second moment than that of the southern region, which could be related to the smaller velocity range covered by the velocity field in that region. Both lines do show a common behavior for the second moment over the map. Specifically, they show a reduction in the line width toward the densest gas of the northern and southern regions.

The B-field observations of 30 Dor, presented in Figure 9, unveil a relatively organized structure in most regions of the cloud. As a result, it defines the main direction for both the northern clump and the southern regions of the map. This direction of the magnetic field in both regions will be used to study the [C ii] kinematics of 30 Dor using position–velocity (PV) diagrams. These two directions, which are aligned with the mean field lines, are presented in Figure 3. Additionally, cuts through 4 regions in the RA direction were made for PV diagrams of [C ii]. The location of these cuts is presented in Figure 5. All the resulting PV diagrams are presented in Figure 6.

Zoom In Zoom Out Reset image size

These PV diagrams confirm that there are several organized VGs in the region. These gradients cover a velocity interval of 5–15 km s⁻¹ in most PV diagrams and come in the form of curves/half-elliptical features that have been associated with expanding-shells (Pabst et al. 2019, 2020; Luisi et al. 2021; Tiwari et al. 2021; Beuther et al. 2022; Bonne et al. 2022a). Based on the presence of these curved features, we identify 4 expanding-shell candidates over the full map using the PV diagrams. The curved velocity structures associated with the candidate expanding-shells are indicated in Figure 6. With the channel maps in Figure 15, we confirm that these candidates have a spatially curved/ring-like morphology as expected for expanding-shell features. We name these candidate expanding-shells: north, east, west, and south after their locations on the map. The location and velocity range of these expanding-shell candidates are presented in Table 1. The location of these shell candidates is also indicated over their velocity range in Figure 15. Two expanding-shell candidates are found in the range of 230–250 km s⁻¹, and two are found in the range of 250–265 km s⁻¹. As some of these shells appear to be only partly covered, larger maps will be required to study their full extent, better characterize their properties, and investigate whether some of the identified shells have the same expansion origin or not. This will also help to study, in combination with data at different wavelengths, whether radiation or stellar winds drive these expanding-shell motions. As the expanding-shells are directly identified by the two blue and red regions in the composite (RGB) image in Figure 4, this shows that the B-field curvature is associated with these expanding [C ii] shell candidates. The highest-density region in the northern region is located at the intersection of two expanding-shells. It might be possible that this region is formed by the collision of the two expanding features. This could result in additional magnetic field bending through oblique shocks (e.g., Inoue & Fukui 2013; Bonne et al. 2020b). However, in the complex dynamical structure, and with the limited spectral [C ii] resolution, it is not easy to establish indications of such B-field bending with confidence. Additionally, the increasing role of gravity might also alter the B-field morphology in this high-density region.

In this section, we computed the VGs of both [C ii] and CO(2-1) from their channel maps adopted from Hu et al. (2021), which is referred to as VChGs (Lazarian & Yuen 2018). Here we recall the principle of this method. The methodology follows these steps: (1) the initial position–position–velocity datacube was preprocessed using principle component analysis by splitting a very large number of velocity channels along the line of sight (LOS) as long as the thin channel map ¹⁸ criteria is satisfied (Lazarian & Pogosyan 2000); (2) the product was convoluted with a 3 × 3 Sobel Kernel to create a raw gradient map (pixels are blanked out if their intensity is less than 3 times root-mean-square level); (3) the gradient angle at each pixel was statistically computed from an adaptive subblock ¹⁹ average method (Yuen & Lazarian 2017a), in which all single gradient orientation within a rectangle subblock in the raw gradient map was taken into account; (4) the pseudo-Stokes Q and U of the gradient were created (Lu et al. 2020), which allows computing the VG morphology. As the principle of VGT, the B-fields are inferred from this technique by rotating the VG orientation by 90° (González-Casanova & Lazarian 2017). However, please note that the zero-angle of the VGs is defined along the east–west direction, which is an offset angle of the polarization angle defined by the SOFIA/HAWC+. Therefore, the VG orientation naturally infers the B-fields in the frame of the SOFIA/HAWC+ thermal dust polarization.

Figure 7 shows the B-fields inferred by the VGTs (red vectors) from both [C ii] (left panel) and CO(2-1) (right panel). To measure the correlation between B-fields inferred from VGTs and SOFIA/HAWC+ (black vectors), we computed the local alignment measurement (AM) as (González-Casanova & Lazarian 2017)

$$\begin{eqnarray}&&\mathrm{AM}=2\left({\cos }^{2}{\theta }_{r}-\displaystyle \frac{1}{2}\right),\end{eqnarray} \tag{ 1 }$$

with θ_r = ∣θ_VG − θ_HAWC+∣ the angle differences between VGT and HAWC+. These two vectors are parallel when AM = 1, while they become perpendicular once AM = −1. AM < 0 indicates where these two vectors are misaligned (i.e., VGs are parallel to B-field lines). One can see that the VGT agrees well with the dust polarization in the south (i.e., VGs are perpendicular to B-field lines), but only fairly in the north. The VGT computed by [C ii] shows that the misalignment occurs mostly at the edge of the regions where the gas density is not peaked (see Paper I). The VGT from CO(2-1) further shows that the misalignment occurs at the denser regions. The reason is that [C ii] traces for more diffuse gas than CO so that it is less affected by gravity.

Zoom In Zoom Out Reset image size

Assuming grains are perfectly aligned with the B-fields, ²⁰ Figure 9 shows the geometry of the fields for all three bands (top to bottom), which is inferred from the polarization orientations. The background color is the original total intensity (Stokes I). One can see that there are two main regions in 30 Dor, these are the north and south regions relative to the massive star cluster R136, whose center is located by the red star. The B-fields are complex but highly ordered, and show a curved feature around the peak intensity (or the peak gas density) in both the north and the south regions. The field lines are curved toward the central cluster. In addition, the field structure at the east side of the peaked intensity in the north region looks like an hourglass. These behaviors are observed similarly in all three bands. The discrepancy among the three panels shown in Figure 9 is mainly due to the fact that the shorter wavelength has higher spatial resolution than those of the longer ones.

The distortion of the B-field morphology could result from the local conditions. For example, this could possibly be explained by the compression of the B-fields by turbulence, shock (Inoue & Fukui 2013), or gravitational contraction (Ewertowski & Basu 2013). To address which is a likely cause of this interesting situation and the local effect of the B-fields, we thus estimate the map of the B-field strength.

We determine the spatial distribution of the B-field strength in 30 Dor according to the strategy developed by Guerra et al. (2021):

• 1.  
  A map of the strength of the B_POS can be determined using the DCF (Davis 1951; Chandrasekhar & Fermi 1953), based on the equipartition of the turbulent magnetic energy and the turbulent kinematic energy approximation:

  $$\begin{eqnarray}&&{B}_{\mathrm{POS}}\simeq \sqrt{4\pi \rho }\displaystyle \frac{\delta {V}_{{\rm{t}}}}{\delta \phi },\end{eqnarray} \tag{ 2 }$$

  where δ V_t is a map of the velocity dispersion of the gas, ρ is the map of the gas density, and δ ϕ is a map of the angle dispersion of the polarization orientation. Following a modification in Houde et al. (2009), the angle dispersion is replaced by the angular dispersion (or structure-function) proportional to the ratio of large-to-small-scale magnetic energies as $\delta \phi \simeq {[\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle ]}^{1/2}$. This improvement is able to avoid the spatial change in the B-field morphology with $\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle$.
• 2.  
  The structure-function dispersion analysis (Houde et al. 2009) was applied to every map of B-field angles (ϕ) in a pixel-by-pixel fashion. This means, for a given pixel i, a circular kernel of radius w-pixel is defined around it. All the pixels inside this kernel are then used to calculate the dispersion function $1-\langle \cos [{\rm{\Delta }}\phi ({\ell })]\rangle$ for pixel i. This dispersion function can be described by the two-scale model:

  $$\begin{eqnarray}&&1-\langle \cos [{\rm{\Delta }}\phi ({\ell })]\rangle =\displaystyle \frac{1-{e}^{-{{\ell }}^{2}/2({\delta }^{2}+2{W}^{2})}}{1+{ \mathcal N }{\left[\tfrac{\langle {B}_{t}^{2}\rangle }{\langle {B}_{0}^{2}\rangle }\right]}^{-1}}+{a}_{2}{{\ell }}^{2},\end{eqnarray} \tag{ 3 }$$

  where Δϕ(ℓ) is the difference between a pair of ϕ separated by the distance ℓ. The first term in Equation (3) describes the small-scale, turbulent B contribution to the dispersion. In this term, δ is the turbulence correlation length, ${ \mathcal N }(={{\rm{\Delta }}}^{{\prime} }({\delta }^{2}+2{W}^{2})/\sqrt{2\pi {\delta }^{3}})$ is the number of turbulent cells along the LOS, and W is the beam's radius (e.g., the σ value) of the polarimetric observations. ${{\rm{\Delta }}}^{{\prime} }$—the cloud's effective depth—is used as a proxy for the depth along the LOS (see Houde et al. 2009 for definition). The second term in Equation (3) quantifies the contribution from the large-scale, ordered magnetic field to the dispersion function.
• 3.  
  The dispersion function in every pixel is fitted with Equation (3) using a Markov Chain Monte Carlo (MCMC) solver. This fitting process can determine values of $\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle$, a₂, and δ in every pixel. The map of $\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle$ is required to estimate the B_POS strength.
• 4.  
  Before combining all three maps through Equation (2), they must have the same angular resolution. In order to do that, the lowest angular resolution (target resolution; σ_T ) among δ V, ρ, and $\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle$ is chosen, and the other maps are smoothed using a Gaussian kernel with width ${\sigma }^{2}\to {\sigma }_{T}^{2}-{\sigma }^{2}$.

The choice of kernel size, w, is important for the fitted values of $\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle$, a₂, δ (if fitted), and therefore the B_POS maps. According to Guerra et al. (2021), the kernel size needs to be large enough to result in statistically significant dispersion functions but small enough for any resulting map to have meaningful angular resolution and avoid smoothing them over large areas. The optimal kernel size, w_opt, can be found by monitoring the distribution of ρ ²¹ values—between the constructed dispersion function and its fit in every pixel. We performed tests using the 214 μm map and values of w = 7, 9, 11 pixels, and found that w = 9 pixels produced the largest fraction of pixels with ρ ∼ 1. Table 2 shows the values of w_opt, global δ, and correlated beam size ($\sqrt{2}W$) for all three bands. From these values, it is clear that the 9 pixel-radius circular kernel contains 7.5 correlated beams and 2–3 turbulence correlation lengths.

Table 2. Values of Turbulence Correlation Length (δ), Correlated Beam Size ($\sqrt{2}W$), and Optimal Kernel Size (w_opt) for the Map-making Procedure (Top)

30 Dor	λ	δ	$\sqrt{2}W$	wopt
	(μm)	('')	('')	('')
	89	16.00	4.68	17.55
	154	21.56	8.17	30.51
	214	23.68	10.93	40.95

North	λ	${{\rm{\Delta }}}^{{\prime} }$	${{\rm{\Delta }}}^{{\prime} }$
	(μm)	(arcmin)	(pc)
	89	0.55	7.84
	154	0.74	10.53
	214	0.76	10.82

South	λ	${{\rm{\Delta }}}^{{\prime} }$	${{\rm{\Delta }}}^{{\prime} }$
	(μm)	(arcmin)	(pc)
	89	0.33	4.70
	154	0.40	5.70
	214	0.38	5.41

Note. Cloud's effective thickness (${{\rm{\Delta }}}^{{\prime} }$) for the north (middle) and south (bottom) regions of 30 Dor. These values were calculated as the half-width half-max of the 1D autocorrrelation of the polarized flux at each wavelength. Linear depths are calculated from angular depths assuming a distance of 49 kpc.

Download table as:  ASCIITypeset image

Using kernels with w > w_opt will not result in better-determined fitting parameters since the addition of more pairs of vectors will not modify the small-scale portion of the dispersion function. Instead, it will result in the contribution of scales larger than those described by the term ∝ℓ², which are not described by the model in Equation (3).

In order to evaluate Equation (2), we used the N(H₂) map from Paper I, and for δ V we used the second-moment map of either CO or [C ii] (Figure 5). To transform column density to mass density, the cloud's depth must be used. In this work, we used the ${{\rm{\Delta }}}^{{\prime} }$ defined above, which is calculated as the width of the autocorrelation function of the polarized intensity. However, taking into consideration that the 30 Dor complex is distinctly separated into two clumps, assuming one single value for ${{\rm{\Delta }}}^{{\prime} }$ can result in less accurate values of B_POS. Therefore, we created a map of ${{\rm{\Delta }}}^{{\prime} }$ that consists of two areas (north and south) each with a constant, different value. These values are displayed in Table 2. On the other hand, the second moment (δ V), calculated from CO and [C ii] observations, contains contributions from thermal motions of the molecules as well as turbulent motions. Therefore, the velocity dispersion values were corrected as ${(\delta {V}_{{\rm{t}}})}^{2}={(\delta V)}^{2}-{k}_{{\rm{B}}}{T}_{\mathrm{gas}}/m$, where k_B, T_gas, and m are the Boltzmann constant, the gas temperature, and the molecule's mass. In this work, we assume 30 Dor is in LTE condition and adopt T_gas = T_ex. We computed the excitation temperature T_ex at each pixel using its maximum main-beam temperature as in Pineda et al. (2008).

According to Guerra et al. (2021), there are two alternatives to fitting the dispersion function (see their Appendix A): (1) One can solve for all three parameters (${a}_{2},\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle ,\delta$) simultaneously. If this provides the majority of pixels in the map with values of δ greater than $\sqrt{2}W$, then the turbulent contribution is well resolved by the polarimetric observations, and the values of $\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle$ are well defined. (2) If the condition $\delta \gt \sqrt{2}W$ is not met for the majority of pixels, the parameter δ can be fixed to a prescribed value, and the MCMC solver computes only two parameters, ${a}_{2},\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle$. For this investigation, the first approach was tried with all polarimetric maps, but only the 214 μm satisfied the condition mentioned above. For the results here presented, we used the second approach, with a prescribed δ value for each pixel equal to that from the global analysis (e.g., when only one dispersion function is calculated for the entire observation). Using the δ-fixed approach results in underestimated values of $\langle {B}_{t}^{2}\rangle /\langle {B}_{0}^{2}\rangle$ and overestimated values of B_POS. However, a test performed with the 214 μm data (the only map for which $\delta \gt \sqrt{2}W$ was met in the majority of pixels) showed that differences in B_POS values using these two approaches were small, with an average value of ∼16%.

The resulting maps of B_POS strength, calculated using the [C ii] data, are presented in Figure 10 for all three bands 89, 154, 214 μm (from the top to bottom). These maps have angular resolutions of 33'', 58'', 77'', respectively. The estimated strengths show a dependency with the wavelength of the polarimetric data. The values range between few and ∼350, ∼450, and ∼600 μG for 214, 154, and 89 μm. This wavelength discrepancy should be mainly because of the angular resolution of the map—shorter wavelengths have smaller resolution. However, there is also a possibility that such differences are the result of each individual wavelength tracing different layers of the cloud. At the same wavelength, the strength structure using CO is similar to the one using [C ii], but the amplitudes are slightly lower (see Figure 16,  Appendix B). The uncertainty of our calculation is discussed in Section 4.5. The black contours in Figure 10 indicate the spatial area in which I/σ_I ≥ 100, and p/σ_p ≥ 3 (e.g., the polarization measurement is statistically significant and intrinsically associated with the source as discussed in Paper I).

An important factor indicating the influence of the B-field is the mass-to-flux ratio. We adopt the ratio as in Crutcher et al. (2004) as

$$\begin{eqnarray}&&\lambda =7.6\times {10}^{-21}\displaystyle \frac{N({{\rm{H}}}_{2})}{{B}_{\mathrm{tot}}},\end{eqnarray} \tag{ 4 }$$

where ${B}_{\mathrm{tot}}=\sqrt{{B}_{\mathrm{POS}}^{2}+{B}_{\mathrm{LOS}}^{2}}$ in μG is the total B-field strength in 3D, which is missing in this work because we do not have the LOS component of the field, and ${N}_{{{\rm{H}}}_{2}}$ in cm⁻² is the gas column density. Statistically, ${\overline{B}}_{\mathrm{tot}}$ could be approximated as $4/\pi \times {\overline{B}}_{\mathrm{POS}}$ (Crutcher et al. 2004). Using the measured Zeeman LOS component of the field, Guerra et al. (2021) showed that B_tot ≃ B_POS in the specific case of Orion Molecular Cloud (OMC-1; see their Equation (11)). Lacking information on the LOS component, we simply adopt B_tot = B_POS with the same caveat as many other studies (e.g., Soam et al. 2018b; Eswaraiah et al. 2021; Hoang et al. 2021b; Ngoc et al. 2021). The gas column density is derived by a modified blackbody fitting from Herschel data at 100, 160, 250, 350, and 500 μm (Meixner et al. 2013) as in Paper I. The cloud is called supercitical for λ > 1, in which the B-field is insufficient to provide support against the gravitational potential. Otherwise, the cloud is called subcritical for λ < 1 in which the B-field prevents the cloud from collapsing.

Figure 11 (left panels) shows the maps of the mass-to-flux ratio for all three bands (from top to bottom) overlaid with the field lines using [C ii]. The similar maps estimated by CO are shown in Figure 16, (right panels.). Generally, the cloud is subcritical (λ < 1) for most of the region, except at the peak intensity in both the north and south (λ ≥ 1), and the regions where λ ∼ 1 spreads larger in area at longer wavelengths. Since λ ∼ N(H₂), and ${B}_{\mathrm{pos}}\sim \sqrt{N({{\rm{H}}}_{2})}$, the gas density is the main uncertainty factor, which is discussed in Section 4.5. Our derivation of N(H₂) agrees quite well with the one from photodissociation region (PDR) models (Chevance et al. 2020). As the maximum of gas column density is about 10²² cm⁻², a hundred μG B-field can make λ < 1 (see Equation (4)).

The interplay between the B-field and turbulence can be defined through the Alfvénic Mach number. Because the DCF assumes that the turbulent magnetic energy balances the turbulent kinematic energy, the total magnetic energy estimated by the DCF method is always greater than the turbulent kinematic energy, which leads to an Alfvénic Mach number that is always lower than 1. Therefore, we use the 3D Alfvénic Mach number to assess the relative importance between the B-field and turbulence as

$$\begin{eqnarray}&&{{ \mathcal M }}_{{\rm{A}}}=\sqrt{3}\displaystyle \frac{{\sigma }_{\mathrm{NT}}}{{{ \mathcal V }}_{{\rm{A}}}},\end{eqnarray} \tag{ 5 }$$

where σ_NT is the nonthermal velocity dispersion, and ${{ \mathcal V }}_{{\rm{A}}}={B}_{\mathrm{tot}}/\sqrt{4\pi \rho }$ is the Alfvénic velocity. ${{ \mathcal M }}_{{ \mathcal A }}\lt 1$ stands for the sub-Alfvénic, the impact of the gas turbulence is minimal, and the B-fields is able to regulate the gas motion. ${{ \mathcal M }}_{{ \mathcal A }}\gt 1$ stands for the super-Alfvénic; the gas turbulence drives the B-field to be random.

The right panel in Figure 11 shows the maps of the Alfvénic Mach number in three bands (from top to bottom) using [C ii], while the ones using CO are shown in the left panels of Figure 17. These maps are color scaled such that the blue color represents ${{ \mathcal M }}_{{\rm{A}}}\lt 1$, and the red color is for ${{ \mathcal M }}_{{\rm{A}}}\gt 1$. One can see that the entire cloud is made up of sub-Alfvénic motions, while it is trans- or super-Alfvénic at around the peaked intensity both in the north and south, except at 214 μm where the the south peak becomes sub-Alfvénic. Similar to mass-to-flux ratio maps, the trans-Alfvénic area is more widespread at longer wavelength.

The magnetic and isotropic turbulent pressures in units of dyn per square centimeter are computed as

$$\begin{eqnarray}\begin{array}{ccl}{P}_{\mathrm{mag}} & = & \frac{{B}_{\mathrm{tot}}^{2}}{8\pi };\\ {P}_{\mathrm{turb}} & = & \frac{3}{2}\rho {\left({\sigma }_{\mathrm{NT}}\right)}^{2}.\end{array}\end{eqnarray} \tag{ 6 }$$

If P_turb/P_mag is greater than 1, the turbulence pressure is higher than the magnetic pressure; otherwise the magnetic pressure is dominant. The maps in all three bands are shown in the left panels of Figure 12, indicating that the turbulent pressure is comparable and higher than that of the B-field at which ${{ \mathcal M }}_{{\rm{A}}}\gtrsim 1$; elsewhere it becomes lower (see Figure 11, right panels).  The similar trends are seen using CO as shown in Figure 17, right panels.

The relative importance between the magnetic and thermal pressure could be quantified through a so-called plasma beta parameter, which is a ratio of these two pressures as

$$\begin{eqnarray}&&\beta =\displaystyle \frac{{P}_{\mathrm{thermal}}}{{P}_{\mathrm{mag}}}=\displaystyle \frac{2{c}_{s}^{2}}{{{ \mathcal V }}_{{\rm{A}}}^{2}},\end{eqnarray} \tag{ 7 }$$

where ${c}_{s}=\sqrt{{k}_{B}{T}_{\mathrm{gas}}/{\mu }_{{{\rm{H}}}_{2}}{m}_{{\rm{H}}}}$ is the thermal sound speed. Figure 12, right panels, shows the maps of the thermal to magnetic pressure ratio. The magnetic pressure is much stronger than the thermal one and varies across the cloud, which is confirmed by the constraints observed by Lee et al. (2019), Chevance et al. (2020) using the PDR Meudon code. Indeed, the authors constrained that P_thermal ∼ 10⁴–10⁶ K cm⁻³, while we showed that P_mag ∼ 10⁷–10⁸ K cm⁻³, with B = 200–500 μG.

The interplay between the gravity, magnetic fields, and turbulence in the star formation has been raised and discussed in numerous studies. The probability distribution function (PDF) of the gas column density was demonstrated to be a key to study the dynamics of molecular clouds. The PDF density of the isothermal and non-self-gravity gas follows closely a log-normal distribution; while this PDF develops a power-law tail once the self-gravity dominates, and the collapse is significant (e.g., Klessen 2000; Federrath & Klessen 2013; Schneider et al. 2013, 2015, 2022; Girichidis et al. 2014; Kainulainen et al.2014).

However, the shape of the PDF strongly depends on the way observations are performed as studied in Körtgen et al. (2019). The authors demonstrated that a small or insufficient field of view will result in a cut-off in the log-normal part of the distribution, and argued that a sufficiently large field of view is needed for the PDF to be built completely. Thus, we used the whole column density map to build the PDF in Figure 13. Then, we fit a log-normal to the histogram as shown by the solid lines. One can see that a log-normal could not fit entirely the PDFs of both the north's and south's gas density, but up to S_t ≃ 0.657, or $N\geqslant 1.8\times {\bar{N}}_{\mathrm{North}}=3\times {10}^{21}\,\,{\mathrm{cm}}^{-2}$ in the north (upper panel), and S_t ≃ 1.239, or $N\geqslant 3.5\times {\bar{N}}_{\mathrm{South}}=3.8\times {10}^{21}\,\,{\mathrm{cm}}^{-2}$ in the south (lower panel). Above these values, power-law tails can be seen in both regions. This reveals that the gravitational instability likely sets in relatively low gas density in 30 Dor cloud in general, and the north experiences it at a slightly lower density than that of the south because its power-law tail pops up earlier. The gravitational collapse in parsec scale has been proposed in various studies (e.g., Hartmann & Burkert 2007; Peretto et al. 2006, 2013; Schneider et al. 2010; Jackson et al. 2019; Bonne et al. 2022b). Another note is that the gravitational collapse in the south is much more compact than that in the north.

Figure 8 shows the location of the protostars candidates adopted from Indebetouw et al. (2009). These locations span from low to the peak gas column density. Interestingly, the protostars' locations are coincident with the misaligned vectors between the rotated VGs by [C ii] and CO and the B-fields, or at which these gas motions are parallel to the field lines.

Zoom In Zoom Out Reset image size

The power-law tails have been seen to develop even in low column density gas, and we suspect that the turbulence could be the key gradient to make that happened. There are three main turbulent types, i.e., compressive (curl-free), solenoidal (divergence-free), and the mixing, in which the compressive turbulence results in stronger compression and thus larger deviation of the PDF (Federrath et al. 2010). These three types of turbulence could be identified through the turbulence driving parameter b defined as

$$\begin{eqnarray}&&{b}^{2}=\displaystyle \frac{{e}^{{\sigma }^{2}}-1}{{M}_{S}^{2}}\displaystyle \frac{\beta +1}{\beta },\end{eqnarray} \tag{ 8 }$$

where M_S is the sonic Mach number, σ is the standard deviation of the log-normal distribution in Figure 13, $\beta =2{M}_{{\rm{A}}}^{2}/{M}_{S}^{2}$ is the compressibility (or the plasma beta parameter in the right panel of Figure 12) with M_A as the Alfvenic Mach number. Purely solenoidal driving has b = 1/3, the compressive driving turbulence yields b = 1, and a stochastic forcing mixture has b ∼ 0.4 (Federrath et al. 2010). Table 3 shows the mean values in the north and south regions. One can see that b ≃ 1 in both regions, which indicates that the driving force of turbulence in 30 Dor is mainly compressive. This is consistent with the existence of the supersonic expanding-shells within the cloud discussed in Section 2.3, and similar to other star-forming regions; such as Serpens cloud (Hu et al. 2021).

Table 3. Mean Values of $\bar{b}\simeq 1$ in Both North and South Regions Indicates That the Turbulence in 30 Dor Is Likely the Compressive Driven

Band ${\bar{T}}_{\mathrm{gas}}$${\bar{{ \mathcal M }}}_{{\rm{A}}}$${\bar{\sigma }}_{v}$σ$\bar{b}$

Turbulent Driving Force
North Region
C	31.68	0.41	8.46	0.608	1.15
D	31.68	0.46	8.46	0.608	1.03
E	31.68	0.61	8.46	0.608	0.78
South Region
C	30.68	0.56	7.60	0.535	0.73
D	30.68	0.52	7.60	0.535	0.78
E	30.68	0.62	7.60	0.535	0.66

Download table as:  ASCIITypeset image

Figures 8 and 14 show that the VGs of both C ii and CO misalign with the B-fields at certain locations in both north and south regions, which indicates that (1) gas moves parallel to the field lines, and (2) there is gravitational collapse. However, these positions distribute from low to maximum gas density, and mostly in magnetically subcritical, sub-Alfvénic regions where the magnetic pressure is higher than that of the turbulence and thermal. The PDF of gas density in Figure 13 illustrates the power-law tail, an indicator of star formation activities, at not very dense gas (column density few times higher than 10²¹ cm⁻²) supporting the local collapse scenario, where protostars candidates are likely to locate. How can we explain the ongoing star formations in strong B-fields?

The turbulence driving parameters b ≃ 1 (see Table 3) reveal that the turbulence is mainly driven by compressive forcing in 30 Dor. The plasma parameter β ≪ 1, (see Figure 12, right panel) shows that the sonic Mach number ${{ \mathcal M }}_{s}\gg 1$, indicating the presence of supersonic turbulence. The supersonic compressive turbulence could create overdense gas in which the gravity gets unstable, and stars can form by directly compressing the gas and/or accumulating the gas along the B-fields line. The latter will not be affected by the magnetic pressure because this pressure acts preferentially perpendicular to the field lines. As a result, the supersonic compressive turbulence could help to trigger the gravitational collapse. That turbulence could result from the cluster-wind from R136 (Melnick et al. 2021) or probably the interaction with much larger giant H i-bubbles in the LMC (Kim et al. 1999). After that, the B-fields get amplified and slow down the star formation activities, as we observe today.

30 Dor contains a spectacular feature with multiple large expanding-shells (see Figure 1). The kinematic of these structures is just about 10⁵⁰–10⁵¹ erg (Chu & Kennicutt 1994), which could be generated by a cluster-wind (Melnick et al. 2021) or a combination with supernovae (Chu & Kennicutt 1994) from R136. Melnick et al. (2021) estimated the virial mass within a distance of 25 pc from R136 as a few times of 10⁶ M_⊙, which is larger than the mass of 30 Dor gas (a few times of 10⁵ M_⊙). Pellegrini et al. (2011) demonstrated that the thermal pressure (in order of 10⁻⁹ dyn cm⁻²) is about 1.6 times lower than that of the radiation pressure (see their Figure 18), and is about 3 times higher than the hot gas pressure (i.e., the hot gas pressure is higher elsewhere; see their Figure 21). Thus, how can the 30 Dor cloud survive?

We suspect that B-fields play a crucial role here in holding the cloud integrity. The B-field morphology orients perpendicular to the radiation direction, so that the magnetic pressure could resist pressure coming from this direction. Indeed, the magnetic pressure is P_mag ≃ 1.6 × 10⁻⁹–10⁻⁸ dyn cm⁻² for B = 200–500 μG, which is absolutely higher than both thermal and radiation pressure and thus is capable to resist against the radiative impact. This calculation is in agreement with the right panel in Figure 12. The role of B-fields could be deduced from a kinematic point of view. Melnick et al. (2021) showed that the kinematic energy of the region within 25 pc from R136 is about 24% the total kinematic energy (which is equivalent to ≃5 × 10⁵⁰ erg ²² ). For B = 200–500 μG, the magnetic energy is orders of magnitude larger as ≃10⁵¹ − 2 × 10⁵² erg within the same spherical volume. In addition, Melnick et al. (2021) showed the turbulent kinematic energy could be up to 10⁵¹ erg; thus, the turbulence only dominates over the B-fields for B < 200 μG, which is consistent with the maps of ${{ \mathcal M }}_{{\rm{A}}}$ (right panel in Figure 11) and P_turb/P_mag (left panel in Figure 12).

The next question could be how large external shells developed around our cloud? 30 Dor opens to the eastern direction from R136, which helps hot gas (traced by X-ray; see Figure 14 in Townsley et al. 2006) easily escape and creates the large expanding-shell along this direction. The mechanism could otherwise be complicated along other directions, especially northern and western. There are multiple expanding-shells on a smaller scale within the cloud (see Figure 6) referring to the fact that the cloud is being carved by the cluster feedback. Melnick et al. (2021) proposed that this impact could lead the certain cloud's structure to break, which enables hot gas to go through and "inflate" the external expanding-shells along these directions. Interestingly, Figure 10 shows that the B-field strength is relatively weak (i.e., lower than 200 μG) at certain regions in 30 Dor along directions that point to or adjoin to these external structures. Along these directions, the field is insufficient to resist feedback; thus, there are nozzles where the hot gas can escape through the cloud.

The larger-scale B-fields have been measured using both near-infrared (NIR; e.g., Nakajima et al. 2007; Kim et al. 2016) and radio polarimetric observations (e.g., Klein et al. 1993; Mao et al. 2012).

The absorption polarization at NIR probed the B-field structure outside our region (see Figure 3 in Nakajima et al. 2007). These large field lines are seen to associate with the dust cloud in LMC and the large expanding-structures around 30 Dor. The B-field strength from this more diffuse region was estimated of 3–25 μG using DCF.

The B-fields measured by synchrotron polarized emission inferred the larger scale of LMC. The low spatial resolution (e.g., $14^{\prime}$ in Mao et al. 2012) was unable to resolve our region. The total B-field strength is measured <7 μG, assuming equipartition with cosmic-ray electrons, with an upper-limit B_POS of 11 μG in the H i southwestern filaments up to the 30 Dor region. The higher B-field strength in the filaments may be due to anisotropic turbulent B-fields mostly located in the plane of sky.

Our maximum B-field strengths of 350–600 μG at parsec scales using thermal dust polarization are much higher than those from the NIR and radio polarimetric observations in more diffuse media, which may be enhanced due to turbulent dynamo mechanism and star formation activity. The additional larger FIR and higher spatial resolution radio polarimetric observations of the LMC are required to better understand the connection between galactic B-fields and those B-fields responsible for the star formation activity and molecular cloud formation.

In this section, we discuss several possible uncertainties. Our first concern lies in the estimate of the column density, which is constructed by the modified spectral energy distribution fitting from Herschel data (see Paper I). Our derived column density agrees well in a range of 1.9 × 10²¹ − 1.3 × 10²² cm⁻² estimation from PDR modeling (Chevance et al. 2020), but 1.3 × 10²² cm⁻² is 1.12 times higher than our maximum derivation. Hence if ${B}_{\mathrm{pos}}\sim N{({{\rm{H}}}_{2})}^{0.5}$, the mass-to-flux ratio λ increases by a factor 1.06, and the Alfvénic Mach number correspondingly decreases by a factor of 0.94, the magnetic pressure increases by a factor of 1.12. The values are slightly changed but do not affect our conclusions.

Our second concern is related to the volume density derived from the column density. The depth of the cloud is fixed and estimated by the autocorrelation function from the polarized intensity. Thus, our derived volume density is insensitive to the local cloud depth, which could be either overestimated or underestimated.

The third concern is the nonthermal velocity dispersion. We computed this term by using [C ii] and ¹²CO emission lines. These lines are not necessarily optically thin and show multiple components; hence it is not ideal for probing the turbulence component of the gas motion even though it shows no evident signs of absorption. If these emissions are optically thick, the nonthermal velocity dispersion might be overestimated because (δ V_turb)² ∼ (δ V)² (e.g., Phillips et al. 1979; Hacar et al. 2016).

The fourth concern is the gas temperature, because (δ V_turb)² ∼ (δ V)² − k_B T_gas/m_gas. We assume the gas is in the LTE condition, so we adopt T_gas = T_ex across the entire cloud. This assumption will fail at the locations where the LTE is simply incorrect. However, we expect that T_gas has a little impact on the final result because (1) the thermal contribution is relatively small (≃2.6 km s⁻¹, for T_gas = 10⁴ K) compared to the total velocity dispersion (Figure 5), and (2) this velocity dispersion is supersonic for [C ii].

The fifth concern is the strength of the B-field in 3D space. In this work, we adopt the value of the POS component of B-field since its LOS component is unknown in 30 Dor. The higher value of the total B-field strength will result in lower values of λ, ${{ \mathcal M }}_{{\rm{A}}}$, while increased in P_mag. Chen et al. (2019) proposed a method to estimate the inclination of B-fields by assuming uniform B-fields (i.e., no variation in position angle and inclination of the field lines) and a constant polarization coefficient (i.e., p₀ in their Equation (10)) throughout a cloud. These assumptions are unlikely to be satisfied in the case of 30 Dor given a variation of the field lines and the polarization degree across the cloud.

The sixth concern is the method to estimate the B-field strength. In this work, we use an improved version of the DCF method from Houde et al. (2009), in which the structure-function of polarization angle is used rather than the dispersion function of this quantity. Skalidis & Tassis (2021), Skalidis et al. (2021) showed that the strength of the mean B-field component is proportional to $1/\sqrt{\delta \phi }$, which is different from Equation (2). However, this approach assumed an equipartition between turbulent motions and the parallel component of turbulent B-fields whose physical basics behind it is not clear (see Appendix A3 in Li et al. 2022; and Section 10.3 in Lazarian et al. 2022). Recently, at the time of this work, Lazarian et al. (2022) proposed a new method to estimate the field strength, namely DMA as $B\simeq f\sqrt{4\pi \rho }{D}_{v}^{1/2}/{D}_{\phi }^{1/2}$ where D_v , D_ϕ are the structure-functions of the velocity centroids and polarization angle, and f is the factor based on the anisotropic turbulence. The authors showed that the DMA method could estimate the field strength at the smallest scale (lower than the turbulence injection scale) and with anisotropic turbulence properties. These two modifications contrast the DCF method and could increase the estimation accuracy. A comparison between DCF and DMA is beyond the scope of this work, and it hence serves as an important improvement for future works (not only for 30 Dor but also other star formation regions).

To summarize, the quantitative value of B-field strengths could be different from what is presented in this work once using a proper value of the nonthermal velocity dispersion and the gas temperature for the DCF or different methods. However, it could not change by orders of magnitude, and we do not expect it will change the effect of the B-field. Indeed, when the strength is as high as hundreds of μG, the field is strong enough to support against the gravitational potential (Equation (4)) and is sufficient to resist the hot gas and radiation pressures (Section 4.3).