A survey of magnetic field properties in Bok Globules

ABSTRACT

Bok globules are small, dense clouds that act as isolated precursors for the formation of single or binary stars. Although recent dust polarization surveys, primarily with Planck, have shown that molecular clouds are strongly magnetized, the significance of magnetic fields in Bok globules has largely been limited to individual case studies, lacking a broader statistical understanding. In this work, we introduce a comprehensive optical polarimetric survey of 21 Bok globules. Using Gaia and near-infrared (IR) photometric data, we produce extinction maps for each target. Using the radiative torque alignment model customized to the physical properties of the Bok globule, we characterize the polarization efficiency of one representative globule as a function of its visual extinction. We thus find our optical polarimetric data to be a good probe of the globule’s magnetic field. Our statistical analysis of the orientation of elongated extinction structures relative to the plane-of-sky magnetic field orientations shows they do not align strictly parallel or perpendicular. Instead, the data is best explained by a bimodal distribution, with structures oriented at projected angles that are either parallel or perpendicular. The plane-of-sky magnetic field strengths on the scales probed by optical polarimetric data are measured using the Davis–Chandrasekhar–Fermi technique. We then derive magnetic properties such as Alfvén Mach numbers and mass-to-magnetic flux ratios. Our findings statistically place the large-scale (A|$_{\mathrm{V}} < 7 \, \text{mag}$|⁠) magnetic properties of Bok globules in a dynamically important domain.

1 INTRODUCTION

Bok globules, named after the astronomer Bart Bok (Bok & Reilly 1947; Bok 1948), are small, isolated dark clouds that play a role in low-mass star formation. They have representative masses of |$\sim$|1–10 |$\mathrm{ M}_\odot$| within |${<}1$| pc and hosts only one or two dense cores (Kandori et al. 2005; Reipurth 2008; Launhardt et al. 2010). Compared to giant molecular clouds with multiple sites of star formation and considerable feedback, Bok globules therefore serve as simpler environments to constrain the physics of star formation. Multiwavelength studies have revealed thermal emission from dust (Clemens & Barvainis 1988; Clemens, Yun & Heyer 1991; Bourke, Hyland & Robinson 1995; Launhardt & Henning 1997; Moreira et al. 1997; Henning & Launhardt 1998; Moreira et al. 1999), submillimeter emission from protostars (Huard, Sandell & Weintraub 1999; Sadavoy et al. 2018; Pattle et al. 2022), and near-infrared observations identifying young stellar objects (YSOs; Alves & Yun 1994; Yun & Clemens 1995; Racca, Vilas-Boas & de la Reza 2009) in some Bok Globules. Molecular line observations provide insights into physical conditions, identifying collapsing clouds and molecular outflows (Yun & Clemens 1992, 1994; Wang et al. 1995; Marka et al. 2012).

The study of magnetic field strength in Bok globules is crucial for understanding the star formation process. Dust grains in the interstellar medium (ISM) are responsible for the polarization of starlight (Hall 1949; Hiltner 1949) and polarized thermal emission (e.g. Planck Collaboration XXXV 2016) caused by such aligned dust grains is the most widely used method for probing the projected magnetic fields on the plane of the sky (Andersson 2015; Tram & Hoang 2022). A growing body of theoretical and observational evidence now supports the radiative alignment torque (RAT) mechanism as the most likely explanation for the observed optical/near-infrared (O/IR) and far-infrared (FIR)/sub-millimetre (submm) polarization in the ISM (Draine & Weingartner 1997; Lazarian & Hoang 2007a), wherein magnetic fields align non-spherical dust grains through radiative torques. O/IR starlight polarimetry efficiently probes the large-scale magnetic field structure of Bok Globules (Jones, Hyland & Bailey 1984; Klebe & Jones 1990; Kane et al. 1995; Sen et al. 2000, 2005; Bertrang, Wolf & Das 2014; Das et al. 2016; Kandori et al. 2020c). FIR/submm polarization maps have revealed magnetic field distribution towards dense cores within Bok Globules (Davis et al. 2000; Henning et al. 2001; Vallée, Greaves & Fiege 2003; Wolf, Launhardt & Henning 2003; Ward-Thompson et al. 2009; Yen et al. 2020; Zielinski, Wolf & Brunngräber 2021; Pattle et al. 2022). These polarization patterns display significant diversity with respect to cloud structure, with the degree of polarization decreasing toward the dense cores. Such studies also applied the Davis–Chandrasekhar–Fermi (DCF; Davis 1951; Chandrasekhar & Fermi 1953) formalism to estimate magnetic field strengths in the range of a few hundred |$\mu$|G (Davis et al. 2000; Henning et al. 2001; Vallée et al. 2003; Wolf et al. 2003; Zielinski et al. 2021), suggesting a important role of magnetic fields in these environments. Polarization studies conducted so far have been largely case studies or involving at most a handful of objects. A dedicated polarization survey that can explore trends across various physical properties of Bok Globules is critically missing.

This paper is based on optical polarization data from the same sample of 21 Bok Globules analysed by Racca et al. (2009). A preliminary analysis of these data shows a tendency of globules having YSOs have less organized magnetic fields than the quiescent ones (Rodrigues et al. 2014). These globules are not associated with bright nebulae or molecular complexes. By systematically analysing for the first time a large sample of starlight polarimetric maps of Bok globules, this study seeks to address several unresolved questions: To what extent do magnetic fields dictate the morphology of Bok globules? Are there systematic differences in magnetic field structure and strength across different evolutionary stages of Bok Globules with and without star formation? These questions are explored by combining Gaia and Planck based dust column density data and optical starlight polarization data, allowing for a more detailed understanding of the interplay between magnetic fields and star formation in Bok globules.

This paper is organized as follows. In Section 2, we outline the observational data utilized in this study. Section 3 details the methodologies adopted, including the extraction of cloud orientations (Section 3.1), the generation of extinction maps (Section 3.2), the modelling of dust grain alignment efficiency (Section 3.4), and the derivation of magnetic field strengths (Section 3.5). Our findings based on these methods are presented in Section 4, followed by an interpretation of the results in Section 5. Finally, we summarize the key conclusions of this work in Section 6.

2 OBSERVATIONS

We have obtained optical polarimetric data of a sample of Southern Bok Globules at the Observatório do Pico dos Dias (OPD), using the 0.6-m Boller & Chivens telescope and the IAGPOL polarimeter (Magalhães et al. 1996), which was equipped with a rotating half-wave retarder plate. The charge-coupled device (CCD) used was an Ikon L, with 2048 |$\times$| 2048 pixels, configured with a 1 MHz readout frequency and a pre-amplifier gain of 4. This reading mode results in a gain of 0.9 electrons per analog-to-digital unit (ADU) and a readout noise of six electrons. The data were collected using a |$I_c$| filter (centred at around 840 nm with a FWHM of 150 nm) to minimize the effect of the interstellar extinction. The images have a field of view of 11 |$\times$| 11 arcmin|$^2$| without the focal reducer and twice this value with the reducer. The observations log is shown in Table 1. Polarimetric standard stars were also observed to correct the instrumental polarization angle to the equatorial system.

The used polarimetric technique splits the incident light in two beams of orthogonal polarizations, producing the so-called ordinary and extraordinary images of a given source. The ratio between the difference and the sum of the counts of those images depends on the polarization of the source and also on instrumental configuration. Specifically, this quantity is a function of the Q and U Stokes parameters and modulates as a function of the retarder position. Therefore, the observed modulation can be used to determine the source polarization. This technique naturally removes the sky polarization from the estimated Stokes parameters of the source. The interested reader can see more on this dual-beam polarimetric technique in Magalhães et al. (1996) and Magalhães, Benedetti & Roland (1984).

The reduction was performed using standard iraf routines (Tody 1986, 1993) to perform bias and dome flat-fields corrections as well as aperture photometry of the ordinary and extraordinary images. The polarization estimates were performed using the pccdpack (Pereyra 2000; Pereyra et al. 2018) and pcckpack_inpe¹  iraf packages. The polarimetric data was originally presented in Magalhães (2012), where the reduction process is described in detail.

3 METHODS

We primarily need two physical quantities for each target – the filament structure (to extract relative orientation w.r.t. the ambient magnetic field) and the average extinctions (that gives us an estimate of the density and mass in the region).

For structures, the ideal resolution and sensitivity is provided by Herschel, but due to unavailability of these for all of our 21 targets, we use Gaia (Gaia Collaboration 2023) stellar density maps in the region (which also correlates well with the Herschel structures, where available).

For extinctions, we can estimate them from near-infrared reddening with 2MASS (M. Skrutskie 2006), or use emission-modelling derived values with Planck (Planck Collaboration 2020). Planck, however, has a low resolution and is also contaminated by background galactic emission. We thus use 2MASS-derived values for our targets, also finding a good agreement with Planck-derived values where computable.

The procedures of creating maps and calculating these quantities is detailed in the following subsections.

3.1 Structure determination

Herschel maps of infrared dust emission at wavelengths of 250, 350, and 500 |$\mu$|m from the SPIRE instrument (Griffin et al. 2010) provided ideal images with fine spatial resolution to identify the dominant structure of each Bok globule. However, Herschel data were available only for six of our 21 globules: BHR 16, 34, 53, 74, 111 and 140. Since this would reduce our sample size drastically, we opted to use an alternative technique based on star counts to ascertain the cloud structure.

We used the optical catalogue of stars from Gaia DR3 and queried around each Bok globule within a square of size 18 arcmin |$\times {}$| 18 arcmin centred on the location of the globule obtained from the BHR catalogue (Bourke et al. 1995). The distances to the globules were taken from the same reference. We found almost no stars present along the line-of-sight closer than the distance (inverse of Gaia parallax) of the globule itself. We also noticed that each of the globules appears as a ‘hole’ devoid of any stars in optical visibility, in the RA–Dec space (see bottom-left panel of Fig. 1 for an example). Thus we attempted to define the boundary of the large-scale structure of the globule using the extent of this ‘hole’ following the procedure outlined below:

• All stars within a square of size |$0.3\times {}0.3{\ }\mathrm{deg}^2$| centred at the globule were queried for from Gaia DR3.

• The corresponding RA–Dec space was divided into a |$60\times 60$| grid of ‘pixels’, each pixel thus measuring |$18\,\mathrm{ arcsec}$|⁠.

• Each pixel was assigned a number equal to the number of stars lying within that pixel. The hole that defined the Bok globule outline had all pixels with values equal to 0. The pixel values were then inverted, i.e. the pixels with zero stars were assigned the maximum value of the 2D array and vice-versa.

• To reduce noise, the pixel values were smoothed by convolving with a 2D Gaussian kernel with |$\sigma =1.5~\rm {}pixels$|⁠.

• Finally, the RA–Dec grid was sectioned into a finer grid of |$300\times {}300$|⁠, each new pixel now being |$3{_{.}^{\prime\prime}} {}6$| in side, and intermediate values were assigned by a 2D cubic spline interpolation

This method worked well enough to return not only the large-scale structure of each cloud, but also several surrounding auxiliary filament-like structures of lower column densities with an accuracy comparable to the 2MASS maps in Racca et al. (2009). We do note, however, that these Gaia stellar density maps in no way can be a proxy for visual extinction maps, and shall only be used for defining cloud structure and boundary. These maps were compared against the Herschel maps for the six clouds where SPIRE data were available and had an excellent match for both the main cloud as well as the less dense filament structures around it as revealed by visual inspection (see, for instance, Fig. 1 for BHR 140). These are hence useful to demarcate the low and high-extinction regions qualitatively.

We then tried to define the cloud’s structure (by broadly taking each cloud to have an elliptical projection on the plane of sky) using the major and minor axis and its position angle (PA). The Gaia maps have a finer resolution where curved structures cannot be fitted easily to an ellipse, so we needed to reduce the resolution to get an elliptical shape. We convolved our original Gaia maps with a Gaussian kernel of |$\sigma =0.08~\rm {}pc$| size, with the angular pixel scale adjusted using Racca et al. (2009) distances.

These low resolution maps were then fed into the filfinder algorithm on python, designed by Koch & Rosolowsky (2015), which effectively isolated the dominant structure as a single filament. This procedure is depicted in Fig. 2. To find the PA, we applied the Rolling Hough Transform (RHT) as described in Clark, Peek & Putman (2014) on the filaments returned by filfinder. The RHT algorithm gives the power of each possible angle along the filament, described in detail at https://seclark.github.io/RHT/. A power-weighted mean PA of all points along the entire filament structure was taken to define the cloud’s PA.

Once the PA is obtained, we found the major and minor axis of the fitting ellipse as follows: the interior (and hence the boundary) of the ellipse was defined to be encompassing all pixels whose value exceeded a certain threshold. In most cases this threshold was set as 0.88–0.94 times the maximum pixel value. Two clouds (BHR 148 and 149) had slightly different thresholds owing to background substructures unrelated to the actual globule. The major axis is then calculated as the length of a straight line of pixels through the centre of the cloud, at the angle given by RHT. The minor axis is found the same way at an angular perpendicular to the PA.

Fig. 2 illustrates the steps of the procedure described above. The obtained PAs for the primary filament are listed in Table 2.

3.2 Extinction maps and mass determination

We use a simplified version of the method outlined in Racca et al. (2009) using 2MASS magnitudes to create extinction maps. Having high spatial resolution was not a primary objective since structure is determined well with Gaia as described above. Our procedure went as follows.

• All stars lying within a |$20\,\mathrm{ arcmin} \times ~20\,\mathrm{ arcmin}$| square centred at the globule’s coordinates are queried from 2MASS.

• The RA–Dec space is divided into a grid of pixels of |$0{_{.}^{\prime}} {}5$| size. This, however, left several pixels empty, so we could not create a map stars only within each individual pixel.

• For the same grid with RA/Dec spacing of |$0{_{.}^{\prime}} {}5$|⁠, we select all stars in a |$1{_{.}^{\prime}} {}5$| square of each grid point. This leads to a repetition in the star sample between neighbouring grid points, but does not return an empty set for any grid point. The average colours of these stars, |$\overline{H-K}$|⁠, where H and K are the standard near-infrared 2MASS bands at 1.65 and |$2.19~\mu \rm {}m$|⁠, was taken as the representative value of reddening.

• The average intrinsic colour of a star, taken as 0.13 from Lada et al. (1994), was subtracted from this |$\overline{H-K}$|⁠. The average reddening is thus |$\overline{H-K}-0.13$|⁠.

• The larger area for selecting stars (⁠|$1{_{.}^{\prime}} {}5$|⁠) ensured a smooth map so separate Gaussian smoothing was unnecessary. However, the pixels were more finely divided into |$3{_{.}^{\prime\prime}} {}6$| each, with cubic interpolation, similar to our procedure for Gaia maps.

• Each pixel number corresponds to |$E(H-K)=\overline{H-K}-0.13$|⁠, where |$E(.)$| denotes the excess due to reddening. We convert it to visual extinction as |$A_\mathrm{ V} = 15.9\times (\overline{H-K}-0.13)$|⁠, as given in Lada et al. (1994).

Racca et al. (2009) however noted that the minimum measurable visual extinction |$A_\mathrm{ V}$| varied from 2 to 4 depending on the error in |$(H-K)$| measurements. We adopt a value of 3 – all pixels with |$A_\mathrm{ V} > 3$| and lying within a circle of a specific radius around the centre represent the Bok globule. The centre coordinates and the radii were tweaked slightly for different globules to accommodate for their sizes and orientations.

To translate the visual extinction to column density and subsequently to masses, we use the relation

as given by Bohlin, Savage & Drake (1978). This gives us the column density for each pixel. Distances from Racca et al. (2009) are used to convert the angular scale to length scale for each cloud, which subsequently gives us the number density of H|$_2$| as |$n_\mathrm{H_2} = N_\mathrm{H_2}/l$|⁠, where l is taken as the projection of |$10\,\mathrm{ arcmin}$| at the cloud distance, since our polarization measurements probe a region of |$10\,\mathrm{ arcmin} \times 10\,\mathrm{ arcmin}$| around the cloud. The mass density is then calculated as |$\rho =1.36\, m_{\rm {}H_2}\, n_{\rm {}H_2}$|⁠, where |$m_{\rm {}H_2}\approx {}2~{\rm {}amu}$| and the factor 1.36 is to account for the distribution of heavier elements including Helium (Kauffmann et al. 2008). Mass is summed over all pixels with |$A_\mathrm{ V} > 3$| and inside the circle.

To test the accuracy of this method, it was first applied on Barnard 68, a nearby and well-studied Bok globule, whose distance and mass are known to be 125 pc and about |$2\, \mathrm{ M}_\odot$| respectively, from Alves, Lada & Lada (2001). Our method with the above adopted distance yielded a mass of |$2.12\,\mathrm{ M}_\odot$|⁠, sufficiently close to the known value (shown in Fig. 3). We proceeded to determine the densities of all the clouds in our sample using this method, with values reported in Table 3.

3.3 Crossmatching with Planck dust emission

The Planck mission measured the total emission intensity across the entire sky. For the Galactic plane and for molecular clouds, the emission intensity primarily comes from interstellar dust. This emission is therefore a useful tracer of column density and was used for the analysis of giant molecular clouds and their magnetic fields in Planck Collaboration XXXV (2016). The emission map along with the all-sky temperature profiles was also converted to a corresponding visual extinction map (Planck Collaboration XXXV 2016) which can be used for both mass estimation as well as average extinction. For Bok globules, however, the primary problem lies with Planck’s resolution; its beam size being as large as |$5\,\mathrm{ arcmin}$|⁠. This results in most of our clouds occupying fewer than 10 pixels in the Planck all-sky map.

We extracted a resampled map of resolution 0.05 deg (3 arcmin) using the reproject_healpix function in python, and then increased the image resolution via cubic interpolation. Most maps include both the Bok globule itself along with extinction from the Galactic background, which should be subtracted to isolate the foreground cloud. This is done by extracting a |$2{_{.}^{\circ}} {}5\times {}2{_{.}^{\circ}} {}5$| region around the Bok globule and subtracting the median emission of this image from each pixel. Wherever this region showed presence of additional high-emission background structures that could potentially distort the median, the region was suitably shrunk until the background could be cleanly seen.

Masses and average surrounding extinctions were calculated with the Gaia stellar density-based structures defining the cloud outline, similar to the procedure for 2MASS extinction maps. Four of our clouds (BHR 59, 133, 138, and 149) lie very close to the Galactic plane in latitude, and also at a longitude near to 0, causing the Galactic disc dust emission to dominate over the cloud in the foreground. The median is also dominated by the Galactic background as a result and becomes very high irrespective of the surrounding region size and median subtraction causes a significant number of pixels to assume negative values of extinction. For these four cases Planck results (both mass and surrounding extinction) are not used any more.

We test the mass determination method again on Barnard 68, and find an excellent agreement of the literature mass of |$2.1\, \mathrm{ M}_\odot$| and our Planck derived mass of |$2.16\, \mathrm{ M}_\odot$|⁠. Additionally, we also find a good agreement to within a factor of 2 in almost all cases between Planck and 2MASS derived values of masses and surrounding extinctions (see Table 6).

The extinction estimated using Planck maps are shown in Table 6. The distribution of ratios of 2MASS and Planck derived extinctions is given in Fig. 4. As can be seen the ratios concentrate around 1 (with a mean of 0.96 and standard deviation of 0.39), implying good agreement between the values derived by two completely different methods (extinction versus emission) in two different wavelength regimes and different physical processes.

We choose to use the 2MASS values of masses and extinctions in our analysis for uniformity, as the high background emission in the case of few clouds in the Planck maps renders the median subtracted extinction maps from Planck unusable for those cases. 2MASS does not have a similar problem as its lower limiting magnitude and near-infrared wavelength regimes means that it only probes the much nearer Bok globules themselves and not the larger background structures farther away that appear in Planck.

3.4 Polarization efficiency

The efficiency with which grains align deep into clouds has been studied by analysing the fractional polarization per unit visual extinction (known as polarization efficiency) as a function of extinction (Goodman et al. 1992).

We use Starhorse derived |$A_\mathrm{ V}$| (extinction) values (Anders et al. 2022). All stars in the polarization data set having SNR below 3 i.e. |$P/\delta P \le 3$| were removed, where P and |$\delta P$| are the polarization percentage and its corresponding uncertainty as measured at the OPD, respectively. Crossmatching with Gaia/Starhorse stars was done within a 1|$\,\mathrm{ arcsec}$| radius of each star with polarization data. Since we only searched for Starhorse stars with a known |$A_\mathrm{ V}$| value, there were no cases where there were |${>}1$| crossmatches in the 1 arcsec  radius. Typical uncertainties in |$A_\mathrm{ V}$| were about |$0.1{\!-\!}0.2$| mag. Any star closer than the distance of the Bok globule itself i.e. part of foreground were removed, although in most cases there were no such stars.

Unlike the global extinction maps in the region around the globule, we do not use 2MASS colours for this analysis. The difference in both analyses is justified by two main reasons.

• For polarization efficiency calculations, we need the extinction values of individual stars. Using 2MASS colours in such cases increases the uncertainty in |$A_\mathrm{ V}$| because the |$A_\mathrm{ V} = 15.9 \times (H-K-0.13)$| assumes the star’s intrinsic colour to be 0.13, which can vary between 0 and 0.25, and lead to up to 2 mag of error in |$A_\mathrm{ V}$|⁠. Starhorse  |$A_\mathrm{ V}$| values are calculated using the star’s physical properties, parallax, and observed colours in its three bands and has a much lower error.

• For extinction maps, the high error in individual stellar extinctions from 2MASS is not as large a problem since the colour is averaged over several stars in a pixel. Further, any analysis (like DCF and mass/column density estimation using average extinction) uses the bulk average of several or all of those pixels, further reducing the overall error. 2MASS, however, offers a key advantage over Gaia/Starhorse of being able to probe the higher extinction regions too, being in near-IR. In fact, the highest extinction central regions in 2MASS correspond to empty regions in Gaia and consequently in Starhorse, rendering its |$A_\mathrm{ V}$| values unusable for determining extinctions in the denser regions of the core.

After the |$A_\mathrm{ V}$| is determined from Starhorse, we find the polarization efficiency as |$P/A_\mathrm{ V}$| where P denotes the polarization percentage. We remove stars with |$A_\mathrm{ V} < 0.2$| as they have high relative errors, which propagate into high errors in polarization efficiency.

We then pick a particular cloud, BHR 121, which has a relatively large number of data points with Starhorse |$A_\mathrm{ V}$| crossmatches, and attempt to see if dust grain alignment modelling can explain the observed values (see Section 4.2). Since this modelling effort requires us to acquire meaningful constraints on the cloud structure, temperature structure and radiation field and thus customize these models to individual sources, we refrain from conducting a complete analysis on all sources. A dedicated effort on this topic would follow.

3.5 Quantifying magnetic field strength

We use the DCF method, as detailed in Chandrasekhar & Fermi (1953) to find the magnetic field strength from the projected plane-of-sky polarization vectors as follows:

Here, |$\rho$| is the mass density of the region we are probing the clouds in g cm|$^{-3}$|⁠, |$\delta v$| is the velocity dispersion in cm s|$^{-1}$|⁠, and |$\delta \phi$| is the polarization angle dispersion in radians. The factor |$f=0.5$| is taken from studies with synthetic polarization maps (Heitsch et al. 2001; Ostriker, Stone & Gammie 2001).

We obtained the |$\delta v$| values from the observation method as described in Otrupcek, Hartley & Wang (2000), using line widths from the |$^{13}$|CO |$(1 \rightarrow {}0)$| transition line, the values given by Tyler Bourke (via private communication) and listed in Table 4. The spectral data used for the calculation of |$\delta v$| is obtained from a 43-arcsec beam centred at the Bok globule. |$\delta \phi$| values are obtained from fitting a Gaussian to the distribution of polarization PAs and taking its standard deviation.

We extrapolate the velocity dispersion using Larson’s line–width size relation (Larson 1981), |$\delta v \propto L^{0.38}$| where L is the size of the region (10 arcmin converted to a length).

The densities |$\rho$| can in principle be obtained directly by using the masses obtained from 2MASS and assuming an ellipsoidal structure and using the structural parameters as obtained from Gaia. However, this approach delivers a density estimate for the dense inner regions of the Bok globules, whereas the polarization data probes the outer, less dense regions. We therefore use a length scale equivalent to |$10\,\mathrm{ arcmin}$| in our density estimates, as described below.

The values of |$n_\mathrm{H_2}$| i.e. the number density obtained by using masses and structures directly were on the order of |$n_\mathrm{H_2} = 10^4{\!-\!}10^6~\rm {}\, cm^{-3}$|⁠, clearly very high values where optical polarization data cannot be obtained due to high extinction. To infer the actual density |$n_\mathrm{H_2}$| of the less dense outer regions where the polarization data comes from, we used the 2MASS maps weighed with the Gaia star densities in the |$10\,\mathrm{ arcmin} \times 10\,\mathrm{ arcmin}$| square centred on the Bok globule. To convert the column densities and masses to volume densities, we assume the third axis (along line-of-sight) to also be equal to the linear length equivalent of |$10\,\mathrm{ arcmin}$| projected at a distance of the respective cloud from Earth, taken from Racca et al. (2009). This implies that the dense central regions where the star density in optical catalogue is 0, is excluded from the density calculations and the outer regions are probed. These returned densities |$n_\mathrm{H_2} = 10^2{\!-\!}10^4$| cm|$^{-3}$|⁠, in agreement of what is expected of the regions (see for example Kandori et al. 2020a, 2020d). The corresponding column density values, |$N_\mathrm{H_2}$| range between |$10^{20}$| and |$10^{22}$| cm|$^{-2}$|⁠. The length scale of 10 arcmin provides a reasonable lower limit for the densities. To ensure that our globules are not significantly larger than the mapped area in our polarization maps, we examined |$1{}^{\circ }\times 1{}^{\circ }$|  Planck-based dust column density maps (Section 3.3). These visual inspections reveal that Bok globules are compact features largely unassociated with any extensive large-scale structure. The absence of discernible structures on larger scales suggests it is unlikely that our globules extend significantly beyond the 10 arcmin along the line of sight. Even so, an extent twice as large along the line of sight would proportionally increase the average density by a factor of 1.2–2. We then have values of |$\rho , \delta v$| and |$\delta \phi$| all coming from the same region around the Bok globule, which allows us to use the DCF method to find the plane-of-sky magnetic fields and the total magnetic field using the relations presented above.

We then further use the values and the same equations as outlined in Pillai et al. (2015) to find the relative strengths of the magnetic field w.r.t. the turbulent motion of the gas and the self-gravity of the cloud using the following two metrics.

• The Alfvén Mach number |$\mathcal {M}_\mathrm{ A}$| that finds the relative strength of magnetic field and turbulence, given by

  $$\begin{eqnarray} \mathcal {M}_\mathrm{ A} = \sqrt{3} \sigma _v/v_\mathrm{ A}, \end{eqnarray}$$

  (2)

  where |$v_\mathrm{ A}$| is the Alfvén velocity given by |$B_{\rm {}tot}/\sqrt{4\pi \rho }$|⁠. Substituting gives |$\mathcal {M}_\mathrm{ A }\propto \delta \phi$|⁠. |$\mathcal {M}_\mathrm{ A}$| less than 1 signifies that the force due to magnetic field dominates over that of turbulent force and is therefore in sub-Alfvénic regime. Similarly, |$\mathcal {M}_\mathrm{ A}>1$| is the super-Alfvénic regime.

• The mass to magnetic flux ratio can be compared to a critical value. The formula is given by

  $$\begin{eqnarray} \frac{(M/\Phi _\mathrm{ B})}{(M/\Phi _\mathrm{ B})_{\mathrm{cr}}} = 0.76 \left(\frac{\mathinner {\langle {N_\mathrm{H_2}}\rangle }}{10^{23} \mathrm{cm}^{-2}} \right) \left(\frac{B_{\rm {}tot}}{1000 \mu G} \right)^{-1}, \end{eqnarray}$$

  (3)

  where |$\mathinner {\langle {N_\mathrm{H_2}}\rangle }$| is the average column density of the region we are probing. In this case again the stellar density maps are used to do a weighted average over the less dense regions of the clouds. If this ratio is less than 1, then magnetic field dominates over the gravity of the system (sub-critical). A value |${>}1$| signifies that gravity is dominant, and is thus super-critical.

3.6 Quantifying error measures

Having metrics for errors was essential to arrive at reliable conclusions for the cloud properties. We needed to quantify errors for the following quantities.

• Distances: Racca et al. (2009) does not have error margins on the distance, however their distance bins are at most 25 pc wide in all cases, which we henceforth adopted as the error on the distances of the individual clouds.

• Velocity dispersions: The data obtained from Tyler Bourke did not have error measures for the velocity dispersions; however their channel spacing for the CO velocity dispersions was |$0.1~\rm {}km\, s^{-1}$| along the line-of-sight. We henceforth adopt an error upper limit of |$0.2~\rm {}km\, s^{-1}$| on the velocity dispersions of individual clouds.

• Mean and sigma of polarization angle (PA): By normal statistics, the error margin on each of these was |$\delta \phi /\sqrt{N}$|⁠, where |$\delta \phi$| is the Gaussian-fit standard deviation on the star sample and N is the sample size.

• Angle of cloud filament: Since this was obtained purely from images themselves generated via star counts and heuristic convolution and binning, an unambiguous quantification is not possible. This quantity is important as it is also a measure of how close to being highly elliptical (filamentary) a cloud is, or how close to being nearly spherical.

  For each cloud with the ellipse identified, we draw a line through the centre on both sides till it touches the boundary of the elliptical structure. The angle from the major axis of the ellipse at which this line acquires a length equal to 95 per cent of length of major axis, is taken to be the error in the dominant angle of the ellipse. The value of 95 per cent was also chosen by inspection; any lower values gave absurdly large values for the angle error and any higher values did not distinguish well between highly elongated versus spherical clouds due to small non-uniformities on the structure boundary (no cloud boundary was shaped like a smooth ellipse) which was used to find the angle error. This is shown in Figure 5 for BHR 75.

  This angle error is also the dominant source of error for the relative angle between the cloud and mean polarization which is detailed in the results section.

• Mass of cloud (M): While not used directly anywhere in further calculations, determination of mass using the structure identification and extinction determination methods detailed above was an important outcome of this work. Hence as an additional check, we found the masses from both 2MASS maps as well as Planck emission (converted to extinction) maps. We find a good agreement of both the measures to within a factor of 2 in all cases except the four outlined earlier, and proceed by assuming the average of the two values to be the true mass of the cloud. Further, we adopt a standard value of error as a factor of 2, which in most cases overestimates the actual error but does not significantly alter the analysis.

• Average extinction in low-density surroundings (weighted |$A_\mathrm{ V}$|⁠): This is again derived from images binned from 2MASS photometric data but devoid of any heuristic cuts (like |$A_\mathrm{ V} > 3z$|⁠). The weighted average extinction can incur errors from mainly two sources: stellar densities of Gaia used for weighing and the extinction derived from stellar passband magnitudes from 2MASS, which have average measurement errors from 0.15–0.2 mag, which with the formula given earlier in Racca et al. (2009) translates to an error in |$A_\mathrm{ V}$| of about 4. We however note that the RA–Dec measurements of stars in Gaia astrometry are extremely accurate compared to this and hence stellar density errors can effectively be neglected. However, the stellar densities by themselves cannot be claimed to act as good weights for the extinction maps, hence we do not have a strictly quantifiable error measure on the average extinction from 2MASS. Hence the same extinctions are then also derived with the Planck maps, but that individually suffers from the strong galactic disc background contamination talked about in the previous section which is estimated by subtracting the median extinction from the entire map. The same stellar density-based weighing and averaging is used for these maps too, and we find a good agreement between the 2MASS and Planck measures within a factor of 2, excepting those four clouds with strong galactic background. We adopt a factor of 2 error in the average extinction as well and the actual average extinction to be the mean of the two values obtained from Planck and 2MASS where both are available, and only 2MASS else. Our results remain robust even to this relatively large error on weighted |$A_\mathrm{ V}$|⁠.

All other quantities are computed using the above main quantities, so errors for the derived quantities was computed using standard Gaussian error propagation methods.

3.7 Magnetic field position angle distribution

In Figs C1–C2, we present the PA (polarization angle) distribution for all the clouds in our sample. In the DCF method, the dispersion of the polarization PAs is the key parameter used to estimate the magnetic field strength. A strongly ordered magnetic field distribution would appear as a Gaussian distribution with a well-defined mean and small angle dispersion, while a weak or randomly oriented field would have a significant dispersion. Several numerical MHD simulations have investigated the uncertainty of the DCF approach in determining field strength (Heitsch et al. 2001; Ostriker et al. 2001) and derived correction factors based on the angle dispersion. The correction factor of 0.5 used in this work (Ostriker et al. 2001) is valid only for B-field angle dispersions, |${\delta \phi < 25^\circ }$|⁠. In our sample of 21 Bok globules, we find 15 to be consistent with having a well-ordered magnetic field, as shown in Figs C1–C2. Two of the clouds, BHR 44 and BHR 145, appear to have randomly distributed fields. Upon visual inspection, the polarization vectors in these clouds appear to trace shell-like structures around the clouds, and attempts to fit a Gaussian result in a poor fit, with dispersions |${\ge} 50^\circ$|⁠, leading to their exclusion from the analysis.

Additionally, four more globules (BHR 53, 59, 148, and 149) show two distinct peaks in the PA distribution histogram. We fit a double Gaussian to account for the bimodal distribution. The two distinct modes are likely influenced by intervening dust from two separate clouds. To confirm this, determining the distances to the stars belonging to each of the two modes is necessary.

We use our Gaia crossmatched stars for sources belonging to each of the two modes to see if they show a clear distinction in distance. For three of the four cases, we find the average distances of the two groups (with two different polarization angles) to have a statistically significant (⁠|${>} 1\sigma$|⁠) difference, as shown for an example in Fig. 6. As any stars whose Gaia-derived distances placed them in front of the target cloud were removed already, both the distance modes are still behind the Bok globule in all cases. The means and standard deviations of each polarization component and the corresponding distances given in Table 5. Since all the Bok globules in our sample lie at distances of 200 to 700 pc, we conclude that the polarization angle of the nearer set of stars behind the cloud is representative of the true polarization in the vicinity of the Bok globule. The farther component is likely affected by other intervening clouds in the line-of-sight, unrelated to our target globule and hence discarded for further analysis. One of the four clouds, BHR 53, did not have a clear |$1\sigma$| difference in the distance of the stars belonging to the two groups of polarization angles, but we still take the one with the nearer distance for further analysis for uniformity.

Additionally, we notice that in all of the four clouds, the stars belonging to the two different PA components have mean polarization fractions that are similar at the |$1\sigma$| level.

With the final components taken from these, and after discarding BHR 44 and 145, we have a clean set of polarization PA distributions with all of them having a dispersion less than |$25^\circ$|⁠.

4 RESULTS

4.1 Polarization efficiency and extinction

A power law of the standard form |$(P/A_\mathrm{ V}) = k(A_\mathrm{ V})^{-p}$| was fitted to the data of polarization efficiency w.r.t. the extinction |$A_\mathrm{ V}$| as obtained in Section 3.4. For the four clouds which have a bimodal distribution of PAs (BHR 53, 59, 148, 149) and indicated two separate cloud structures at different distances influencing the polarization (discussed in Section 3.7), we only used the stars belonging to the nearby component for the power law fitting. We also observed that even if we perform the fitting with all stars instead of just the ones in the nearby component, the best-fitting index (p) changes in value by less than 0.15 in all cases.

We took the index p of the star-forming and starless cores separately and found their variance-weighted mean, where the variance is returned by the curve_fit function itself.

The star-forming clouds had a mean power law index of |$-0.70$| and standard deviation |$\sigma = 0.32$|⁠. The starless cores have a mean |$= -0.39$| and |$\sigma _\mathrm{ p} = 0.24$|⁠.

There is thus only a marginally statistically significant difference in the power-law indices for the star-forming cores and the starless cores. A complete effort on individual modelling of the clouds with physical constraints would be needed to explain this trend.

4.2 Modelling grain alignment efficiency

To place our power law fits the polarization efficiency above in context, we utilize the dustpol-py model² for an isolated cloud (without an internal radiation source, see Tram et al. 2025 for details regarding the model basis) to analyse BHR 121 – a cloud with extensive data on starlight polarization and corresponding Gaia-derived extinctions. The model is fundamentally based on the radiative torque theory for grain alignment (RAT-A; see Dolginov & Mitrofanov 1976; Draine & Weingartner 1996; Lazarian & Hoang 2007a; Hoang & Lazarian 2008; Andersson, Lazarian & Vaillancourt 2015 for further details) and presumes that grains are entirely aligned with the magnetic fields. The most important parameters in this model are the gas volume density and temperature, the radiation field and the dust temperature, and the grain composition and shape.

We utilize a cylindrical configuration incorporating a Plummer density profile across the core of the globule. The gas volume density (⁠|$n_{\rm H}=10^{5}\, \rm cm^{-3}$|⁠) and changes as a function of the distance (r) from the centre following |$n_{\rm H} \sim r^{-2}$|⁠. Our cloud is embedded in the standard ISRF with the mean wavelength of 1.3|$\, \mu$|m, and the surrounding envelope temperature of 15 K. The local radiation field, mean wavelength and dust temperature within the cloud are computed through the radiative process. We make use of the oblate astrodust grain composition (Hensley & Draine 2023) with the axial ratio of 1.4. The grain size is varied to investigate the impact of the largest grain size on polarization efficiency across different wavelengths. Fig. 7 presents the best models that well cover our data.

As seen in Fig. 7, 95 per cent of our observed data lies on or between the four model curves for different grain sizes. For |$A_\mathrm{ V} \gtrsim 0.7$|⁠, the observed polarization efficiencies can be well explained by the typical radiation fields and dust grain sizes found in such environments. At lower extinctions, the polarization efficiencies are systematically higher than the model predictions, implying that they are still aligned with, and hence are a probe of, the local magnetic field, but additional physics like a stronger radiation field or more elongated grains might be needed to explain the higher polarization efficiencies observed.

In Fig. 7, we also show the ‘small grain’ model (e.g. |$a_\mathrm{max}=0.1\, \mu$|m to |$a_\mathrm{max}=0.15\, \mu$|m where |$a_\mathrm{max}$| is the maximum grain size) that exhibits a slope of |$-1$| at the upper end of the visual extinction range we probe (⁠|$A_\mathrm{ V}$| of 2 to 6). To provide a reference, we include an arbitrarily scaled line representing |$P/A_\mathrm{ V} \propto A_\mathrm{ V}^{-1}$|⁠, plotted over the same |$A_\mathrm{ V}$| range and positioned below the model predictions. This addition demonstrates that, within the majority of the observed |$A_\mathrm{ V}$| range for BHR 121, grain alignment is not lost at optical wavelength. Furthermore, it supports the conclusion that the maximum grain size in this region exceeds that of typical ISM grains.

4.3 Relative angles between cloud and field

To determine the relative orientation of clouds and magnetic fields, we use the RHT algorithm (Section 3.1) to calculate the dominant orientation of each cloud. For the magnetic field, we calculate the mean angle of polarization vectors by fitting a Gaussian distribution to their PA data. The relative angle is then the difference between these two values, computed modulo |$180^\circ$|⁠.

Since there is an inherent uncertainty of approximately |$10^\circ$| (corresponding to our 3 |$\sigma$| cutoff on polarization angles) in each angle measurement, we construct a histogram of the relative angles using |$10^\circ$| bins. As shown in the upper panel of Fig. 8, this allows us to analyse the distribution of relative orientations more robustly.

Additionally, we examine potential correlations between the relative angles and the cloud mass, as derived from 2MASS observations. This step helps us investigate whether the relative orientation depends on cloud mass or follows a specific trend. The results are shown in the lower panel of Fig. 8.

To ensure the reliability of the analysis, we exclude the following clouds: BHR 44 and BHR 145. These clouds exhibit a polarization PA dispersion greater than |$25^\circ$|⁠, which means the Gaussian fitting does not yield a clear peak. The resulting mean PA is unreliable and likely random. For BHR 138, the Gaia stellar density structure map shows a nearly circular shape, which lacks a distinct filamentary structure necessary for defining a dominant orientation.

After these exclusions, we proceed with the analysis using a final sample of 18 clouds.

Following the approach in Stephens et al. (2017), we generate the cumulative distribution (CDF) of the observed projected angles between our starlight derived B-field orientations and the filament orientation. The observational results are shown as the blue step curve in Fig. 9. To find the actual distribution of 3D angles between the magnetic field vectors and filament structure vectors, we performed Monte Carlo simulations to generate |$5\times 10^4$| 3D vector pairs, representing the orientations of the filament axis and the magnetic field respectively. The angle between each vector pair was taken uniformly from |$(0^\circ , 20^\circ)$| with probability |$\alpha$| (cases where the filament and magnetic fields are parallel to each other), and uniformly from |$(70^\circ , 90^\circ)$| with probability |$1-\alpha$| (cases where they are perpendicular to each other). The value of |$\alpha$| is varied from 0 to 1 in steps of 0.05. Viewing angles are drawn uniformly from a 3D sphere, and the distribution of the projected plane-of-sky angle between each vector pair is calculated. We then apply the Anderson–Darling (ad) test (Anderson & Darling 1954) to determine which distribution, from a set of simulated distributions ranging from purely parallel, a mix of parallel and perpendicular, to purely perpendicular relative angles, most accurately fits our data. We also perform the ad test using simulated vector pairs for a random orientation, where their relative angles in 3D are uniformly drawn from the range |$(0^\circ , 90^\circ)$|⁠. The p-value in the ad test quantifies the probability of observing the given data if it were truly drawn from a hypothesized distribution. In the case of a random alignment of true 3D vectors, with a p-value of 0.055, we can reject, at the 10 per cent confidence level, the hypothesis that the vectors are randomly distributed.

Fig. 9 shows the results of the simulated distributions as the set of yellow curves and random distribution as the black dashed line. Our CDF is inconsistent with a random alignment of the true 3D vectors (below a 5.6 per cent confidence level), as well as with a purely parallel or perpendicular alignment (below 0.1 per cent confidence level for both). The family of brown curves in Fig. 9, however, agrees well with our observed distribution where a bimodal distribution of a |$30~{{\ \rm per\ cent}}{\!-\!}70~{{\ \rm per\ cent}}$| to |$50~{{\ \rm per\ cent}}{\!-\!}50~{{\ \rm per\ cent}}$| mix of parallel versus perpendicular 3D vector alignments with a p-value greater than 25 per cent. (The caps of 0.1 per cent and 25 per cent on the p-value are imposed by the scipy function used for the test, and essentially represent thresholds where we can confidently accept or reject the respective hypotheses. The p-value for the range (⁠|$30~{{\ \rm per\ cent}}{\!-\!}50~{{\ \rm per\ cent}}$| of parallel alignments) capped at 0.25, thus indicates a significantly better match of the bimodal distribution with the observed distribution than a purely random distribution.

4.4 Magnetic field properties

Here, we present the key findings related to the distribution of magnetic field strengths, mass-to-flux ratios, and Alfvén Mach numbers, which provide insights into the magnetic support and turbulence within these Bok Globules, which is also summarized in Fig. 10, and presented in full in Table 6.

The total magnetic field strength, |$(B_{\text{pos}}$|⁠), derived using the DCF method (refer to equations 2–3), shows a range of values from 23 to 296 |$\mu$|G. As shown in Table 6, the mass-to-magnetic flux ratios, which quantify the balance between gravitational forces and magnetic support, are found to be in the 0.03–1.51 range All except one globule exhibit mass-to-magnetic values |${<}1$|⁠. Considering the observational uncertainties, this sample of Bok Globules appear to be close to being magnetically critical such that magnetic pressure may provide support on cloud scales against gravitational collapse. That does not necessarily prevent star formation occurring in dense cores within these clouds. The one cloud with a ratio exceeding 1 is BHR 133, which has an exceptionally high 2MASS-derived mass of |$393\, \mathrm{ M}_\odot$|⁠, found to be a starless globule (Racca et al. 2009). This also happens to be the cloud with the largest distance at 700 pc, suggesting that distance uncertainties could cause a significant overestimation of the total mass.

The Alfvén Mach number, |$(\mathcal {M}_\mathrm{ A}$|⁠), derived from the ratio of the turbulent velocity dispersion to the Alfvén speed, also exhibits a distribution in the range 0.26 to 0.92. All globules show sub-Alfvénic values |$(\mathcal {M}_\mathrm{ A} < 1$|⁠) except the two with high PA dispersion (BHR 44 and 145, where the Alfvén Mach number is anyways unreliable as DCF does not work for those clouds). This indicated that magnetic fields are strong enough to influence the cloud dynamics. This suggests that optical polarimetry probes a regime where magnetic forces dominate over those due to turbulence, as can be seen in Fig. 10.

5 DISCUSSION

Understanding the role of magnetic fields in the formation and evolution of dense molecular clouds is essential to advancing our knowledge of star formation Our optical polarimetric findings towards a sample of 21 Bok Globules reveal key insights into the magnetic field morphology, and their stability.

5.1 Grains are efficiently aligned

A characterization of grain alignment is necessary for accurately analysing magnetic field strengths and turbulence in the ISM through polarization studies. In one of our Bok Globules with a substantial number of detections (BHR 121), we observe that grain-alignment efficiency declines with increasing extinction, following a power-law index, |$< -1$|⁠. This finding indicates that, on scales probed by optical polarimetric data, dust grains are mostly aligned with the magnetic field of the cloud. Using a physical model of dust polarization driven by radiative torques toward BHR 121, we confirm that dust grains remain efficiently aligned at a few magnitudes of visual extinction. The polarization efficiency is consistent with previous optical polarimetric study in three Bok Globules (Chakraborty, Das & Paul 2014). Our results also align with previous studies on dense molecular clouds, which suggest that grain alignment can be maintained by radiative torques even in dense environments (Matthews & Wilson 2000; Henning et al. 2001; Alves et al. 2014; Kandori et al. 2018; Soam et al. 2018; Coudé et al. 2019; Pattle et al. 2019; Pillai et al. 2020).

A comprehensive modelling effort, similar to that conducted for BHR 121, across our full sample may provide deeper insight into the alignment of dust grains, dust grain distribution and the extent to which grain growth or other factors influence alignment efficiency in Bok Globules with varying environmental conditions and density structures.

5.2 Globules are preferentially extended parallel or perpendicular to the magnetic field

An effective way to use polarization measurements to understand the role of magnetic fields is to compare the inferred orientation of the magnetic field with the orientations of elongated molecular cloud structures (Tassis et al. 2009; Li et al. 2013; Planck Collaboration XXXV 2016). Bok Globules, as simple and isolated systems, are not significantly affected by large-scale turbulence or nearby star-forming events, making them ideal for probing this relationship. While several studies on magnetic fields in Bok Globules exist (Sen et al. 2000; Ward-Thompson et al. 2009; Bertrang et al. 2014; Choudhury et al. 2022), large-scale systematic studies on orientation of the Bok Globules versus B-field orientation remain limited. Ward-Thompson et al. (2009), found a 40 deg offset between the magnetic field orientation and the short axis in two Bok Globules, while Chakraborty & Das (2016) find the observed magnetic field in CB130 to be almost aligned perpendicular to its minor axis and Prokopjeva et al. (2014) found magnetic field to be oriented parallel to the Bok Globule’s filamentary short axis in CB67. Our survey results for 21 targets is now able to provide context to these case studies. The projected angles between large-scale B-fields and globule long axes do not align purely parallel or perpendicular, showing instead a distribution consistent with a bimodal alignment, where the alignment is parallel only 30–50 per cent of the time. Note that we cannot entirely rule out purely random orientation at |${\sim} 5$| per cent probability.

Chen & Ostriker (2014) suggested that gas flows along magnetic field lines, such that its orientation is perpendicular to the main filament. Subsequent theoretical work has further investigated the relative orientation between magnetic fields and column density structures (Chen, King & Li 2016; Soler & Hennebelle 2017; Seifried et al. 2020) finding a relation for strongly magnetized regimes, with Seifried et al. (2020) highlighting the impact that projection effects can have in assessing relationships between cloud structures and their magnetic fields, even when they exist.

Observational studies using Planck polarization maps have shown that, in most molecular clouds in nearby Gould Belt regions, B-fields are mostly parallel to filaments at low column density, and a transition occurs from parallel to perpendicular magnetic field alignment with increasing density (Planck Collaboration XXXV 2016). A similar relation is observed in nearby regions where column density structures are preferentially aligned with magnetic fields (Fissel et al. 2019; Lee et al. 2021; Bij et al. 2024). This transition from parallel to perpendicular alignment typically occurs at a visual extinction of 3–5 mag (Planck Collaboration XXXV 2016).

In our present work, we lack sufficient background stars at or beyond this visual extinction to directly detect such transitions. Deeper near-infrared and longer wavelength polarimetric data that can probe the cores of Bok Globules are needed to confirm the existence of this transition. However, our bimodal distribution results suggest that many of the Bok Globules in our sample exhibit at least moderately strong magnetic fields, indicating they are likely either trans- or sub-Alfvénic.

5.3 Magnetic forces dominate over turbulence and can balance self-gravity

Our DCF based estimates as shown in Table 6 further strengthens the case for an important role of magnetic fields. The plane-of-sky magnetic field strength in our sample ranges from 23 to 296 |$\mu$|G. Our analysis indicates that the magnetic field within the Bok globule is strongly sub-Alfvénic and sub-critical, which has significant implications for its evolutionary path and star formation potential. In the sub-Alfvénic regime, magnetic forces outweigh turbulent gas motions, potentially hindering large-scale collapse and stabilizing the globule against gravitational contraction. The sub-critical nature of the globule further suggests that magnetic support prevents gravitational collapse. The free-fall time-scale for globule collapse under gravity is approximately 0.1 to 1 Myr. However, certain globules in our sample are known to host YSOs (Racca et al. 2009), suggesting that while magnetic support exists on a global scale, it is not sufficient to counteract gravitational forces within the densest cores, enabling localized star formation.

A standard observational test distinguishing strong fields from that of weak fields is the scaling relation between the magnetic field strength and the density it probes: |$B \propto n^{\kappa }$| (Crutcher 2012). In the strong-field model (Mouschovias & Ciolek 1999), the collapse of a molecular cloud driven by ambipolar diffusion primarily compresses gas along magnetic field lines, increasing the density without significantly amplifying the magnetic field strength (⁠|$\kappa \le 0.5$|⁠). In contrast, weak-field models of cloud collapse, assuming flux freezing, suggest that the magnetic field is too weak to restrict gas flows in any particular direction, leading to a scaling of |$\kappa \approx 0.65$| (Mestel 1966). Observational data, based on a compilation of magnetic field strengths measured using Zeeman splitting and their associated densities, align more closely with the weak-field model predictions (Crutcher 2012). Since our magnetic field strength measurements are based on the DCF method, in Fig. 11, we compare our findings with an ensemble of DCF-based results from Pattle et al. (2023). The reference line represents the theoretically expected magnetic field-density relation (⁠|$B \propto n^{0.65}$|⁠) for weak fields (Crutcher 2012). The blue points in Fig. 11 show DCF-based Bok Globule observations compiled from the handful of previous studies (Vallée et al. 2003; Wolf et al. 2003; Bertrang et al. 2014; Chakraborty & Das 2016; Das et al. 2016; Choudhury et al. 2019; Kandori et al. 2020a, b, c; Myers & Basu 2021). Our results mostly lie above the n|$^{0.65}$| relation expected under pure flux-freezing, however it stays within the upper bounds of the literature compilation. They also broadly align with the limited observations of Bok globules available at similar densities. The few existing DCF-based B-field measurements at higher densities (from sub-mm polarimetric data) reveal magnetic fields well below the theoretical reference limit for those densities. The is consistent with gravitational collapse conditions in weak field models (Crutcher 2012). However a systematic DCF based study of star-forming regions and magnetic field-density relation by Myers & Basu (2021) provide a statistical explanation for the same trend when considering detailed physical properties of these cores under moderately strong magnetic fields.

Our finding makes a key assumption that the DCF method reliably estimates B-field strength. While these approaches are well-established, future research could benefit from multiwavelength polarimetric observations to more accurately capture B-field morphology across varying cloud depths. Additionally, high-resolution observations – potentially from instruments like ALMA – could offer new insights into small-scale magnetic field variations and further clarify how magnetic fields dissipate to allow star formation that is observed within some of these Bok Globules (Racca et al. 2009).

6 CONCLUSIONS

Our study underscores the role of magnetic fields in the formation and stability of Bok globules. Optical polarimetric data from 21 Bok globules probing their envelopes reveal that magnetic fields in these regions are typically strong and sub-Alfvénic, meaning magnetic forces largely counteract turbulence, potentially stabilizing the globules against at least large-scale gravitational collapse. Grain alignment modelling for BHR 121 confirms that radiative torques maintain dust grain alignment up to moderate extinction levels, suggesting that optical polarimetry can be a valuable tracer of magnetic field structure in these simple systems. Our findings also suggest that magnetic fields are not purely parallel or perpendicular to globule structures but are in best agreement with a bimodal alignment. The magnetic field strengths, ranging from 23 to 296 |$\mu$|G, align with past observations of dense clouds as well as limited studies of Bok Globules. Future studies with deeper infrared polarimetry and higher resolution observations could better probe these B-field and density interactions, providing more detailed insights into how magnetic support affects star formation within Bok globules.

ACKNOWLEDGEMENTS

TGSP gratefully acknowledges support by the National Science Foundation under grant no. AST-2009842 and AST-2108989 and by NASA award #09-0215 issued by USRA. We thank Dan Clemens and Phil Myers for helpful comments. CVR thanks the Brazilian Space Agency (AEB) by the support from PO 20VB.0009 and the Conselho Nacional de Desenvolvimento Científico e Tecnológico – CNPq (Proc:310930/2021-9). This paper is based on observations made at the Observatório do Pico dos Dias managed by Laboratório Nacional de Astrofísica (Brazil). This publication makes use of data products from the Two Micron All Sky Survey, which is a joint project of the University of Massachusetts and the Infrared Processing and Analysis Center/California Institute of Technology, funded by the National Aeronautics and Space Administration and the National Science Foundation. This work has also made use of data from the European Space Agency (ESA) mission Gaia, processed by the Gaia Data Processing and Analysis Consortium. Funding for the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.

DATA AVAILABILITY

The optical polarization data obtained from the Observatório do Pico dos Dias is available upon reasonable request. All other data (Gaia, Planck, Herschel and 2MASS) are obtained from publicly available surveys.

Footnotes

REFERENCES

APPENDIX A: GAIA OVERLAYED IMAGES

The Fig. A1–A2 show the stellar density maps obtained for each Bok globule, along with the polarization vectors of background starlight overlayed in blue.

APPENDIX B: 2MASS OVERLAYED IMAGES

The Fig. B1–B2 show the extinction maps obtained for each Bok globule using 2MASS reddening, along with the polarization vectors of background starlight overlayed in green. The lengths of the vectors are proportional to the polarization fraction.

APPENDIX C: DISTRIBUTION OF POLARIZATION ANGLES

The Fig. C1–C2 show the distribution of polarization PAs around each cloud, along with the best-fitting Gaussian profiles for each.

APPENDIX D: DISTRIBUTION OF POLARIZATION FRACTIONS

The Fig. D1–D2 show the distribution of the polarization fractions around each Bok globule.