Magnetic fields in the infrared dark cloud G34.43+0.24

TL;DR

This study maps magnetic fields in the infrared dark cloud G34.43+0.24 using polarized dust emission, revealing ordered field geometries aligned with the filament and estimating field strengths that suggest magnetic influence on cloud dynamics.

Contribution

First detailed magnetic field mapping of G34.43+0.24 across multiple scales using various instruments, with new estimates of magnetic field strength and criticality parameters.

Findings

Magnetic fields are ordered and aligned perpendicular to the filament in key regions.

Magnetic field strengths are estimated at 470, 100, and 60 microGauss in different regions.

G34.43+0.24 is marginally critical with sub-Alfvénic turbulence, indicating magnetic fields influence cloud contraction.

Abstract

We present the B-fields mapped in IRDC G34.43+0.24 using 850\,$\mu$m polarized dust emission observed with the POL-2 instrument at JCMT. We examine the magnetic field geometries and strengths in the northern, central, and southern regions of the filament. The overall field geometry is ordered and aligned closely perpendicular to the filament's main axis, particularly in regions containing the central clumps MM1 and MM2, whereas MM3 in the north has field orientations aligned with its major axis. The overall field orientations are uniform at large (POL-2 at 14$\arcsec$ and SHARP at 10$\arcsec$) to small scales (TADPOL at 2.5$\arcsec$ and SMA at 1.5$\arcsec$) in the MM1 and MM2 regions. SHARP/CSO observations in MM3 at 350\,$\mu$m from Tang et al. show a similar trend as seen in our POL-2 observations. TADPOL observations demonstrate a well-defined field geometry in MM1/MM2 consistent with MHD simulations of accreting filaments. We obtained a plane-of-sky magnetic field strength of 470$\pm$190\,$\mu$G, 100$\pm$40\,$\mu$G, and 60$\pm$34\,$\mu$G in the central, northern and southern regions of G34, respectively, using the updated Davis-Chandrasekhar-Fermi relation. The estimated value of field strength, combined with column density and velocity dispersion values available in the literature, suggests G34 to be marginally critical with criticality parameter $\rm \lambda$ values 0.8$\pm$0.4, 1.1$\pm$0.8, and 0.9$\pm$0.5 in the central, northern, and southern regions, respectively. The turbulent motions in G34 are sub-Alfv\'{e}nic with Alfv\'{e}nic Mach numbers of 0.34$\pm$0.13, 0.53$\pm$0.30, and 0.49$\pm$0.26 in the three regions. The observed aligned B-fields in G34.43+0.24 are consistent with theoretical models suggesting that B-fields play an important role in guiding the contraction of the cloud driven by gravity.

Tables

Table 1. Table 1: Values calculated in G34 Center, North, and South regions.

Region	${\sigma }_{\theta }$	$\Delta v$	$n_{H_2}$	$B_{\mathrm{pos}}$	${\lambda }_{\mathrm{corr}}$	$M_A$
	(∘)	($\mathrm{km}{s}^{-1}$)	(${\mathrm{cm}}^{-3}$)	($\mu$G)
Center	11$\pm$4	1.1$\pm$0.1	$1.8\cdot {10}^5$	470$\pm$190	0.8$\pm$0.4	0.34$\pm$0.13
North	16$\pm$9	0.8$\pm$0.2	$0.6\cdot {10}^5$	100$\pm$40	1.1$\pm$0.8	0.53$\pm$0.30
South	15$\pm$8	0.6$\pm$0.2	$0.2\cdot {10}^5$	60 $\pm$34	0.9$\pm$0.5	0.49$\pm$0.26

Table 2. Table 2: Parameters derived from modified DCF methods without and with correction for beam integration .

		Without correction		With correction
		for beam integration		for beam integration
Paramters	Description	SF	ACF	SF	ACF
$\Delta \theta {(}^o)$	Angular dispersion	$26.0\pm 0.3$	$23.2\pm 2.7$	$52.0\pm 2.8$	$54.0\pm 5.1$
$⟨\delta B^2⟩/⟨B_0^{}{}_{}{}^2⟩$	Turbulent-to-ordered magnetic field energy ratio	$0.21\pm 0.02$	$0.16\pm 0.01$	$0.90\pm 0.04$	$0.81\pm 0.01$
$B_{pos}$ ($\mu$G)	Plane-of-sky magnetic field strength	$300\pm 120$	$200\pm 70$	$150\pm 90$	$90\pm 50$

Equations

$$P=\frac{1}{I}\sqrt{Q^2+U^2-\frac{1}{2}(\delta Q^2+\delta U^2)},$$

$$\delta P=\sqrt{\frac{(Q^2\delta Q^2+U^2\delta U^2)}{I^2(Q^2+U^2)}+\frac{\delta I^2(Q^2+U^2)}{I^4}}.$$

$$\theta =\frac{1}{2}\mathrm{ta}{\mathrm{n}}^{-1}(U/Q).$$

$$\delta \theta =\frac{1}{2}\frac{\sqrt{Q^2\delta U^2+U^2\delta Q^2}}{(Q^2+U^2)}\times \frac{180^{\circ }}{\pi }.$$

$$I_{\nu }=B_{\nu }(T_{\mathrm{d}})(1-{\mathrm{e}}^{-{\tau }_{\nu }}),$$

$$B_{\nu }(T_{\mathrm{d}})=\frac{2h{\nu }^3}{c^2}\frac{1}{{\mathrm{e}}^{h\nu /k_{\mathrm{B}}{T}_{\mathrm{d}}}-1},$$

$${\tau }_{\nu }={\mu }_{{\mathrm{H}}_2}{m}_{\mathrm{H}}{\kappa }_{\nu }{N}_{{\mathrm{H}}_2},$$

$$B_{pos}=Q\sqrt{4\pi \rho }\frac{{\sigma }_v}{{\sigma }_{\theta }},$$

$$B_{pos}\approx 9.3\sqrt{n({\mathrm{H}}_2)}\frac{\Delta v}{{\sigma }_{\theta }}\mu \mathrm{G},$$

$$\frac{\delta B_{pos}}{B_{pos}}=\frac{1}{2}\frac{\delta n({\mathrm{H}}_2)}{n({\mathrm{H}}_2)}+\frac{\delta \Delta v}{\Delta v}+\frac{\delta {\sigma }_{\theta }}{{\sigma }_{\theta }},$$

$${\lambda }_{obs}=\frac{(M/\varphi {)}_{obs}}{(M/\varphi {)}_{crit}},$$

$$(M/\varphi {)}_{obs}=\frac{\mu m_H{N}_{{\mathrm{H}}_2}}{B_{pos}},$$

$$(M/\varphi {)}_{crit}=\frac{1}{2\pi \sqrt{G}}.$$

$$M_A=\frac{\sqrt{3}{\sigma }_v}{{\sigma }_A},$$

$${\sigma }_A=\frac{B_{tot}}{\sqrt{4\pi \rho }},$$

$${⟨\Delta {\varphi }^2(l)⟩}_{tot}\simeq b^2+m^2{l}^2+{{\sigma }_M}^2(l),$$

$$\frac{{⟨\delta B^2⟩}^{1/2}}{B_0}=\frac{b}{\sqrt{2-b^2}},$$

$$B_0\simeq \sqrt{(2-b^2)4\pi \mu m_H{n}_{H_2}}\frac{{\sigma }_v}{b}.$$

$$B_{pos}=Q{B}_0.$$

$$1-⟨cos[\Delta \varphi (l)]⟩\simeq \frac{1}{N}\frac{⟨\delta B^2⟩}{⟨B_0^2⟩}\times [1-e^{-l^2/2({\delta }^2+2W^2)}]+{a_2}^{\prime }{l}^2,$$

$$N=\frac{({\delta }^2+2W^2){\Delta }^{\prime }}{\sqrt{2\pi }{\delta }^3},$$

$$B_0\simeq \sqrt{4\pi \mu m_H{n}_{H_2}}{\sigma }_v[\frac{⟨\delta B^2⟩}{⟨B_0^2⟩}{]}^{-1/2}.$$

Full text

stars: formation — ISM: kinematics and dynamics — ISM: magnetic fields

1 Introduction

Filamentary structures exist in molecular clouds, with sizes ranging from a few to tens of parsecs (André et al., 2014; Wang et al., 2016). Recent magnetohydrodynamic (MHD) simulations (Klassen et al., 2017; Li et al., 2018; Gómez et al., 2018) probing the formation of large-scale filamentary clouds suggest a complex evolutionary process involving the interaction and fragmentation of dense, velocity-coherent, fibers into chains of cores, resembling observations in nearby clouds (e.g., L1495/B213 and Musca cloud; Hacar et al., 2013, 2016). The simulations show that global magnetic fields are expected to be roughly perpendicular to the longer axes of dense filamentary clouds. Several velocity coherent fibers are identified inside the clouds and appear to be supportable along the main filament. In 3D MHD simulations of cluster-forming turbulent molecular cloud clumps, Klassen et al. (2017) found that B-fields are oriented parallel to sub-virial clouds and perpendicular to denser gravitationally-bound clouds.

Recent ideal MHD simulations by Li & Klein (2019) found that the magnetic field helps in shaping the long filamentary structures with field orientation perpendicular to the long axis of the filaments. Their simulation produces fibre-like substructures similar to those observed in L1495 (Hacar et al., 2013). There are some other MHD simulations available which include magnetic fields in filaments. Federrath et al. (2016) presented MHD simulations studying the effect of magnetic fields, gravity and turbulence on the formation of filaments finding that filament width does not depend on the orientation of filament with respect to the magnetic fields in G0.253+0.016 region. A statistical analysis of nearby clouds such as Taurus, Musca, Ophiuchus, and Chameleon has revealed that B-field lines tend to become parallel to the filament long axes at low densities (or “diffuse” with a few $\rm cm^{-3}$) (e.g., Chapman et al., 2011; Cox et al., 2016; Planck Collaboration et al., 2015, 2016) but are perpendicular to the denser filamentary structures. Koch et al. (2014) presented a statistical analysis of 50 sources (from 4000 independent measurements observed with the SMA and the CSO) on the scales of 0.1 pc to 0.01 pc with densities $\rm\gtrsim 10^{5}~{}cm^{-3}$. Their analysis of B-fields and intensity gradients reveals that the field orientation is perpendicular to the sources’ major axes.

Polarized thermal dust emission at submillimeter wavelengths probes the magnetic field structure in high-density regions. The Radiative Torque Alignment (RAT) theory of grain alignment is currently one of the most promising models to explain the polarization of light towards clouds and cores (Dolginov & Mitrofanov, 1976; Lazarian et al., 1997). This model predicts the asymmetric non-spherical dust grains rotate due to radiative torque and align with their long axes perpendicular to ambient magnetic field. Due to low angular resolution (e.g., $\sim 5{\arcmin}$ with Planck) or high dust extinction (optical or near-infrared polarimetery), previous studies of magnetic fields in filamentary clouds have been mostly limited to nearby clouds. So far, magnetic fields have only been investigated in a few infrared dark clouds (Pillai et al., 2015; Liu et al., 2018a; Juvela et al., 2018; Liu et al., 2018b). Additional observations with higher angular resolution towards filamentary clouds and cores are still needed.

To this end, we are conducting a series of dust polarization observations toward the brightest filaments identified in the JCMT legacy survey of $\sim$1000 Planck Galactic Cold Clumps (PGCCs), called SCOPE (SCUBA-2 Continuum Observations of Pre-protostellar Evolution; Liu et al., 2018c; Eden et al., 2019), with the POL-2 polarimeter at the JCMT. The observational results of two PGCCs, G35.49-0.31 (hereafter G35) and G9.62+0.19 (hereafter G9) are published in Liu et al. (2018a) and Liu et al. (2018b), respectively.

In this work, we report POL-2 observational results toward a more evolved filament, G34.43+0.24 (hereafter G34). At a distance of $\sim$3.7 kpc (Sanhueza et al., 2012; Foster et al., 2014; Xu et al., 2016), G34 is an active high-mass star-forming filamentary cloud (Molinari et al., 1998; Rathborne et al., 2011; Sakai et al., 2018). G34 harbors multiple cores, including G34-MM1 through MM9, that are likely at different evolutionary stages (Chen et al., 2011). Figure 1 shows the locations of these MM sources. G34-MM2 was found to be the most evolved core (Rathborne et al., 2006) associated with the ultra-compact H II (UCH II) region IRAS 18507+0121 of spectral type B0.5 (Molinari et al., 1998; Shepherd et al., 2004, 2007). Investigating the cores in G34, Rathborne et al. (2008) found that the brightest millimeter core, G34-MM1, exhibits a typical chemical signature of a high mass core. On the other hand, the clump MM3 hosts a hot-corino (Yanagida et al., 2014; Sakai et al., 2015). Chambers et al. (2009) found G34-MM1, MM3, MM4, MM5, and MM8 associated with extended Spitzer 4.5 $\mu$m emissions, indicating possible outflow activities. Sanhueza et al. (2010) also observed these cores and found molecular outflows associated with cores G34-MM1, MM2, MM3, and MM4. G35 is a filament similar to G34 with several embedded low-luminosity massive protostars (Nguyen Luong et al., 2011) and massive starless clumps (Liu et al., 2018a). A network of filaments covering a broad range of densities is also revealed in G35. The magnetic field lines in G35 tend to be perpendicular to the densest part of the most massive filament, whereas they tend to be parallel in the low-density regions as well as in other less dense filaments (Liu et al., 2018a). The magnetic fields together with turbulence, however, do not appear able to support against the gravitational collapse of the densest clumps in G35. The northern region of G34 with MM3 is associated with the PGCC G34.50+0.27. G34 has a mass per unit length of $\sim$1600 M*⊙* pc*-1* (Xu et al., 2016) which is about four times larger than that ($\sim 410$ M*⊙* pc*-1*) of G35 (Liu et al., 2018a). By comparing G34 with G35, we can determine which of the three mechanisms, B-fields, gravity, or turbulence is dominant in filament evolution and dense core formation.

2 Observations, data acquisition, reduction and validation

The POL-2 observations were conducted in August 2018 (M18AP041; PI: Soam A.) in Band-2 weather conditions using the POL-2 daisy map mode of JCMT (Holland et al., 2013; Friberg et al., 2016; Bastien et. al. in prep., 2019) at 850 $\mu$m. In this mode of observations, a fully sampled circular region of 11$\arcmin$ diameter is produced with a high signal-to-noise coverage over the central 3$\arcmin$ of the map. This observing mode is based on the SCUBA-2 CV daisy scan pattern (Holland et al., 2013) but modified to have a slower scan speed (i.e. 8$\arcsec$/s compared to 155$\arcsec$/s) to obtain sufficient on-sky data for good Stokes Q and U values. Coverage decreases, with a consequent significant increase in the rms noise, toward the edges of the map. The POL-2 polarimeter with a rotating half-wave plate at a frequency of 2 Hz and a fixed polarizer is placed in the optical path of the SCUBA-2 camera. The total on-source integration time was $\sim$3 hours with $\tau_{225}$ ranging from 0.05 to 0.08, where $\tau_{225}$ is atmospheric opacity at 225 GHz. We adopted the same observational strategy as described by Ward-Thompson et al. (2017). POL-2 simultaneously collects the data at 450 $\mu$m and 850 $\mu$m wavelengths with full-width half maximum (FWHM) of 9.6$\arcsec$ and 14.1$\arcsec$, respectively (Dempsey et al., 2013). We have not reported 450 $\mu$m data in this work since the instrumental polarization (IP) model for 450 $\mu$m data was not commissioned when this project started.

The data were acquired from the Canadian Astronomy Data Center (CADC) and reduced using the STARLINK/SMURF package pol2map (Chapin et al., 2013; Currie et al., 2014) specifically developed for reducing sub-millimeter data obtained from JCMT. The details of the data reduction steps and procedure are described in Wang et al. (2019). In the first run of pol2map, the raw bolometer time-streams for each observation are converted into separate Stokes Q, U, and I time-streams using the process calcqu. Then a Stokes I map is created from all I time-streams using a routine makemap which is an iterative map-making process in SMURF package. Individual I maps corresponding to each observations were coadded to produce the initial I map of the region. The details of this step can be seen in Chapin et al. (2013). The final I, Q, and U maps were obtained by running pol2map a third time. The initial I map described in a previous step is used to generate a fixed SNR-based mask for all further iterations of makemap. The pointing corrections determined in the previous step were applied during the map-making process. During the final process, we invoked an additional parameter called skyloop111http://starlink.eao.hawaii.edu/docs/sc22.pdf in pol2map and corrected for the loss of synchronization between data values and pointing information in the data reduction process. Skyloop improves the recovery of faint, extended structures by performing one iteration of the mapmaker on all of the observations, co-adding the result, and testing each successive iteration for convergence (see Wang et al., 2019). This is in contrast to the traditional map-making method of deriving an iterative solution for each observation individually. The final polarization values obtained here are de-biased by using the mean of Q and U variances to remove statistical bias in regions of low signal-to-noise ratio (SNR).

The calibrated I, Q, and U maps were obtained in $\rm Jy\,beam^{-1}$ units using a Flux Calibration Factor (FCF) of 537 Jy $\rm pW^{-1}$ given for 850 $\mu$m. The output maps are multiplied by 1.35 to account for additional losses due to POL-2 (Dempsey et al., 2013; Friberg et al., 2016). The final co-added total intensity map has an rms noise222This value was measured using SCUBA2-MAPSTATS recipe under PICARD package in STARLINK of $\rm\sim\,7.0\,mJy\,beam^{-1}$. The rms noise in Q and U maps were found to be $\rm\sim\,7.9\,mJy\,beam^{-1}$ and $\rm\sim\,6.8\,mJy\,beam^{-1}$, respectively.

After the final step of running pol2map, we obtain a polarization vector catalogue which is produced by co-adding Stokes I, Q, and U maps. The data were reduced with a 4$\arcsec$ pixel size but to improve the sensitivity, we binned the co-added Stokes I, Q and U maps to 12$\arcsec$ pixel size using binning over 3$\times$3 pixels.

The debiased (Wardle & Kronberg, 1974) polarization fraction values were estimated (see Soam et al. (2018); Wang et al. (2019)) as

[TABLE]

where P is the debiased polarization fraction and I is the total intensity. Q, U, $\delta Q$, and $\delta U$ are the Stokes parameters and their uncertainties. The uncertainty in polarization fraction is estimated using

[TABLE]

The polarization position angles were measured increasing towards the east from the north in the sky projection using relation

[TABLE]

The corresponding uncertainties in $\theta$ were calculated using

[TABLE]

The plane-of-sky B-field orientation is inferred by rotating polarization angles by 90∘ (assuming that the polarization is caused by elongated dust grains aligned perpendicular to the magnetic field). We have used only the data points where the observed uncertainties in position angles are less than 20∘. The large-scale B-fields are examined using Planck 850 $\mu$m (353 GHz) dust emission polarization maps (Planck Collaboration et al., 2015, 2016). The image is smoothed to the 7$\arcmin$ resolution to ensure good SNR data. The vectors are drawn at 3.5$\arcmin$ (half-resolution) steps.

We checked the quality of our data used for analysis by examining the different SNR values derived from polarization intensity (PI) and its uncertainty ($\sigma_{PI}$). In panel (b) of Figure 2, the B-fields inferred from SNR$>$ 2 (PI/$\sigma_{PI}>2$; 211 red vectors) and SNR$>$3 (PI/$\sigma_{PI}>3$; 146 purple vectors) are generally consistent in the regions where both are available. The other two panels (c) and (d) of Figure 2 show comparisons of the distributions of B-field position angles and polarization fraction of the two subsamples. The aim is to test the validity of the data with $\rm 2<SNR<3$ (which is generally used in such studies) and to decide whether or not data with SNR$\geq$2 could be used for studying B-field morphology and strength. The very similar trends in distributions of position angles and polarization percentages reassures us that we can use the $\rm 2<SNR<3$ data for further analysis.

3 Results and discussion

The Stokes I map of the G34 filament in 850 $\mu$m continuum emission with inferred B-field geometry is shown in panel (b) of Figure 2. The elongated shape of the filament is clearly visible and three regions of interest (‘N’, ‘C’, and ‘S’) are indicated by labeled dashed black rectangles. The overall filament appears to have a small head to the North (containing MM3), a dense clump (consisting of MM1 and MM2) in the center, and a diffuse tail-like structure to the south.

3.1 Magnetic field morphology

Panel (a) in Figure 2 represents the large-scale B-fields inferred from Planck measurements in the region containing G34 (Planck Collaboration et al., 2015). There is a clear indication of field lines aligned in the south-west to north-east directions. Panel (b) of the figure shows zoomed-in B-fields in G34 from our POL-2 observations at sub-parsec scales. The lengths of line-segment are proportional to the fractional polarization. Magnetic field geometry and properties are studied individually in the regions center (C), north (N), and south (S) labelled in panel (b). The northern part containing MM3, has field orientations almost along the elongated clump. The central region, however, has field lines perpendicular to the long axis of the filament with an apparent smooth change in orientation when seen from west to east. The southern diffuse region has most field lines closely perpendicular to the tail. The large-scale field in the northern region from Planck observations is also nearly parallel to the filament (see left panel (a)), which is similar to the fields seen in the region ‘N’ from POL-2 observations. This suggests that the B-field is connected from parsec to sub-parsec scales, despite orders of magnitude difference in density and the physical scales. However, it should also be noted that compared to G34, the region measured by Planck next to it is mostly background and foreground. Hence, it is not evident a priori that the field orientations are identical. We found that the background subtraction would tend to make Planck polarization vectors more perpendicular to the filament but details depend on the selection of the reference regions chosen to represent the background and the filament remains unresolved in the Planck data.

Figure 3 shows a better view of magnetic field morphology in the G34 filament where we use the normalized vectors with their lengths independent of the polarization fraction. The smooth change in field lines from being perpendicular to almost parallel from the center to north regions can be clearly seen in this figure.

Figure 4 shows the Gaussian fitted distributions of B-field position angles in the center, north, and south regions. The distributions in all regions peak around $75^{\circ}$, which is close to an east-west orientation.

3.2 Dust temperature and column density

We estimated the dust temperature ($T_{\mathrm{d}}$) and $\mathrm{H_{2}}$ column density ($N_{\mathrm{H_{2}}}$) of the G34 filament using archival Herschel PACS/SPIRE (70, 160, 250, 350, and 500 µm) and JCMT 850 µm data fitted with a modified black-body function. In this procedure, the different resolution Herschel images and JCMT 850 µm image were smoothed to the SPIRE 500 µm FWHM beam size of 35″and reprojected on the same grid. The G34 filament is found embedded in a large-scale molecular cloud in Herschel images causing additional emission from surrounding material in the line-of-sight. In order to obtain an accurate column density for G34, this background was subtracted. Then the spectral energy distribution (SED) was fitted to the fluxes obtained in Herschel and JCMT maps for each pixel position using the formulae (Kauffmann et al., 2008).

[TABLE]

[TABLE]

[TABLE]

where $B_{\nu}(T_{\mathrm{d}})$ is the Planck function at a given dust temperature ($T_{\mathrm{d}}$), $\tau_{\nu}$ is the optical depth, $\mu_{\mathrm{H_{2}}}$ is the mean molecular weight per hydrogen molecule, $m_{\mathrm{H}}$ is the hydrogen atom mass, $\kappa_{\nu}$ is the dust opacity (absorption coefficient), and $N_{\mathrm{H_{2}}}$ is $\mathrm{H_{2}}$ column density. The value of $\mu_{\mathrm{H_{2}}}$ is 2.8 and $\kappa_{\nu}$ for each used frequency are 1.76, 0.4, 0.195, 0.1, 0.05, and 0.0197 cm2 g*-1*, respectively, adopted from Ossenkopf & Henning (1994) for a dust-to-gas ratio of 0.01. The temperature and column density maps of G34 made using this procedure are shown in Figure 5. The temperatures (left panel of Figure 5) throughout the filament vary from $\sim$10$-$25 K with hot dust present in the central region containing MM1/MM2 and colder in the north and southern regions. The column density values (right panel of Figure 5) are found peaking at $\sim 10^{23}\,cm^{-2}$ in the central region.

We also estimated the $\rm H_{2}$ volume densities of the three regions of G34 north, center, and south assuming them to have cylindrical geometry and adopting the procedure explained in section 3.2 of Liu et al. (2018a). The projected lengths (L) of the cylinders corresponding to ‘N’, ‘C’, and ‘S’ regions of G34 shown in middle panel of Figure 2 are 1.9, 3.0, and 2.6 pc, respectively. The mean values of the projected radius (r) of circular ends of these cylinders are measured to be 1.8, 2.2, and 1.1 pc, respectively. We used these values to estimate volumes of the cylinders and their number densities. The estimated values of volume densities are shown in Table 1.

The estimated values of column and volume densities are used in sections 3.3.1 and 3.3.2 for further calculations.

3.3 Magnetic Field, Gravity, and Turbulence in G34

3.3.1 Magnetic field strength

We estimate the plane-of-sky B-field ($\rm B_{pos}$) strengths in the central dense region ‘C’, north ‘N’, and south ‘S’ regions of G34 using the Davis-Chandrasekhar-Fermi relation (DCF; Davis 1951; Chandrasekhar & Fermi 1953). The DCF relation assumes a regular field geometry with dispersion indicating a measure of the distortion in the field geometry caused by turbulence. The vector distribution is considered to be Gaussian with a well-characterized standard deviation. The DCF method is represented by the expression

[TABLE]

where $\rho$ is the gas density, $\sigma_{v}$ is the observed velocity dispersion of the gas, and $\sigma_{\theta}$ is the dispersion in polarization angle. The DCF field model assumes that $Q$ is a factor of order unity that accounts for variations in the B-field on scales smaller than the beam. Ostriker et al. (2001) compared their mean values of the known plane-of-the-sky magnetic fields with DCF estimates and found $Q$ in the range of 0.46-0.51. They suggested that the DCF estimate, modified by a multiplicative factor of $\sim$0.5 to account for a more complex magnetic field and density structure, can provide an accurate value of B-field strength when polarization angles are quite uniform. Therefore, we adopted $Q$ as 0.5 for our calculations. Following the simplification introduced by Crutcher et al. (2004), eq. 8 can be written as

[TABLE]

where $n({\rm H}_{2})$ is the number density of molecular hydrogen in cm*-3*, $\Delta v$ = $\sigma_{v}\sqrt{8\ln 2}$ is FWHM in km s*-1* and $\sigma_{\theta}$ is in degrees.

The uncertainty in field strength is measured by combining uncertainties using the relation

[TABLE]

where $\delta n({\rm H}_{2})$, $\delta\Delta v$, and $\delta\sigma_{\theta}$ are the uncertainties in $n({\rm H}_{2})$, $\Delta v$, and $\sigma_{\theta}$, respectively.

Tang et al. (2019) has recently investigated the magnetic field strengths in MM1, MM2, and MM3 regions using CSO/SHARP polarization data at 350 $\mu$m wavelength and $\rm N_{2}H^{+}$(J=1-0) line observations. We used their velocity dispersion in $\rm N_{2}H^{+}$(J=1-0) line with FWHM $\Delta v$= 1.1$\pm$0.1 $\rm~{}km~{}s^{-1}$. We estimated the average volume densities as $\rm\sim 1.8\times 10^{5}~{}cm^{-3}$ in the central region containing MM1/MM2, $\rm\sim 0.6\times 10^{5}~{}cm^{-3}$, and $\rm\sim 0.2\times 10^{5}~{}cm^{-3}$ in the north and south regions, respectively. We subtracted the mean value of position angles from all measured position angles (with $\rm PI/\sigma PI>3$) in region ‘C’ giving residual angles ($\rm\delta\theta$). The Gaussian fit to the distribution of $\rm\delta\theta$ values provides a dispersion in position angle ($\rm\sigma_{\theta}$) of $\rm 12.9\pm 4.0^{\circ}$ (see Figure 6). We corrected the value of $\rm\sigma_{\theta}$ by mean value of observed position angle uncertainties which is measured to be $\rm 7.0^{\circ}$. The corrected value of $\rm\sigma_{\theta}$ is $\rm\sqrt{12.9^{2}-7.0^{2}}\approx 10.8^{\circ}$. This value of dispersion in position angles satisfies one of the assumptions of the DCF relation which limits the maximum value of $\rm\sigma_{\theta}$ to be $\rm\leq 25^{\circ}$ (Heitsch et al., 2001). Using the above-mentioned values of FWHM in $\rm N_{2}H^{+}$(J=1-0), $n({\rm H}_{2})$, and $\rm\sigma_{\theta}$, the strength of $\rm B_{pos}$ in G34 center is found to be 470$\pm$190 $\mu$G. This field strength is similar to that found in other IRDCs such as $\sim$270 $\mu$G in G11.11-0.12 (Pillai et al., 2015) and $\sim$100 $\mu$G in G035.39-00.33 (Liu et al., 2018a) but smaller than the value 790 $\mu$G in G9.62+0.19 (Liu et al., 2018b).

We have used a similar approach as described above to estimate the magnetic field strengths in the North and the South regions of G34. The values of velocity dispersion i.e. FWHM of $\rm N_{2}H^{+}$(J=1-0) in these regions are adopted from Tang et al. (2019) for estimating the field strength using eq. 9. Dispersion in position angles towards these regions is estimated similarly as described above for the central region. The uncertainties in $B_{pos}$ are estimated using eq. 10 which is derived from eq. 9 by propagating the errors in the quantities. We found the values of $\rm B_{pos}$ as 100$\pm$40 $\mu$G and 60$\pm$34 $\mu$G in ‘N’, and ‘S’ regions, respectively.

3.3.2 Mass-to-flux ratio

We will use our estimates of magnetic field strength to calculate the standard parameters of mass-to-flux ratio ($\rm M/\phi_{B}$) and Alfvénic Mach number ($\rm M_{A}$). These measure the relative importance of magnetic fields versus gravity and turbulence, respectively.

$\rm M/\phi_{B}$ is the ratio of the mass (M) of the object to the flux ($\rm\phi_{B}$) of the magnetic fields threading the object. Crutcher et al. (2004) (and references therein) discussed that the maximum mass that can be supported by a given magnetic flux is known as critical mass, $\rm M_{crit}=\frac{\phi_{B}}{2\pi\sqrt{G}}$. We tested the importance of the B-field in the context of gravity in all three regions of G34 where plane-of-sky B-field strength is estimated using eq. 9. This can be investigated by calculating the value of criticality parameter ($\rm\lambda_{obs}$) using the relation

[TABLE]

where the observed mass-to-flux ratio is estimated as

[TABLE]

and $\rm\mu$, $\rm m_{H}$ and $N_{\mathrm{H_{2}}}$ are the mean molecular weight per $\rm H_{2}$ molecule, mass of atomic hydrogen, and molecular hydrogen column density, respectively. The average values of column densities in the center, north, and south regions are found to be $\rm\sim 15\times 10^{22}$, $\rm 4.5\times 10^{22}$, and $\rm 2.5\times 10^{22}~{}cm^{-2}$, respectively.

The clouds that are not collapsing due to the support by magnetic fields are called magnetically “subcritical” ($\rm\lambda<1$), whereas those with gravity which overcomes the support of the magnetic field are referred as magnetically “supercritical” ($\rm\lambda>1$).

[TABLE]

Using the column density in $\rm cm^{-2}$ and measured B-field strength in $\rm\mu G$, we estimated the value of $\lambda_{corr}$ after applying a geometric correction to $\lambda_{obs}$ following Crutcher (2004). The value of $\lambda_{obs}$ can be overestimated by a factor of 3 due to geometrical effects suggesting $\rm\lambda_{corr}=\lambda_{obs}/3$. The errors in $\lambda_{obs}$ come mainly from the uncertainty in B-field strength. We calculated the $\lambda_{corr}$ in all three regions of G34 and the results are given in Table 1. The values of $\rm\lambda_{corr}$ obtained for the center, north, and south show that these regions are transcritical. All the values are close to criticality (i.e. $\lambda=1$) suggesting that gravity and magnetic fields are equally important in these regions.

3.3.3 B-fields and turbulence

The nature of turbulent motions in the G34 clump can be studied by estimating the value of $\rm M_{A}$ which describes the relative importance of magnetic fields and turbulence in molecular clouds (Padoan et al., 2001; Nakamura & Li, 2008). When the fields are uniform and strong, the turbulence is regulated by the magnetic fields, yielding a sub-Alfvénic scenario (with $\rm M_{A}\leqslant 1$). On the other hand, if the cloud is super-Alfvénic (i.e., $\rm M_{A}>1$), the magnetic field is not strong enough to resist scrambling by turbulent motions. The value of $\rm M_{A}$, using molecular line and polarization observations, can be estimated as

[TABLE]

where $\rm\sigma_{v}$ is the mean non-thermal velocity dispersion, measured from the FWHM (i.e. $\rm\sigma_{v}=FWHM/\sqrt{8ln2}$) of $\rm N_{2}H^{+}$(J=1-0) line observations Tang et al. (2019) which we used in estimating magnetic field strength. $\rm\sigma_{A}$ is the Alfvénic velocity calculated as

[TABLE]

Crutcher et al. (2004) found from a statistical study that the total magnetic field strength ($\rm B_{tot}$) is 1.3 times the plane-of-sky field strength. In the absence of knowledge of the 3D geometry of G34, this is a reasonable correction to apply. The value of $\rm\sigma_{A}$ in three regions of G34 is calculated using different magnetic field strengths and volume densities in these regions. Using the values of $\rm\sigma_{v}$ and $\rm\sigma_{A}$, we calculate $\rm M_{A}$ in all three regions of G34.

The value of dispersion in position angle, FWHM of $\rm N_{2}H^{+}$(J=1-0) line, volume density, estimated plane-of-the-sky magnetic field strength ($\rm B_{pos}$), projection corrected mass-to-flux ratio ($\lambda_{corr}$), and Alfvénic Mach number ($\rm M_{A}$) in regions ‘C’, ‘N’, and ‘S’ are given in Table 1. The values of $\lambda_{corr}$ in all three regions of G34 suggest it to be marginally critical. The values of Alfvénic Mach number suggest the sub-Alfvénic nature of turbulence in G34 filament.

Some analytical studies investigated the stability and fragmentation of filaments in context of turbulent motions (Ostriker, 1964; Inutsuka & Miyama, 1992; Toalá et al., 2012; Heitsch, 2013) and B-fields (Heitsch, 2013) and found B-fields are important in filament formation. Soler et al. (2013) studied the dependence of B-fields on the initial magnetization of filament using combination of synthetic polarization maps and numerical simulations of magnetized clouds and concluded that strong compression is caused by super-Alfvénic turbulence. Whereas, the sub-Alfvénic turbulence allows the gravitationally collapsing material to move along the B-field lines (Nakamura & Li, 2008). In case of G34, turbulence is found to be sub-Alfvénic in all three regions of G34 and the field lines found to be perpendicular to the elongated axes. This is mostly true in the central and southern regions whereas in the northern region the field appears to change from perpendicular to parallel. However, the field orientation seem mostly perpendicular in the lower part of the northern region. Tang et al. (2019) found the similar geometry in G34 north part near MM3 (see Figure 1 for the location of MM3) and suggest that B-fields must be playing different role here than in central region. Liu et al. (2018a) notices field lines getting parallel from perpendicular in northern region of an IRDC G35 suggesting that fields in that region are likely to be poloidal. Similar trend in north of G34 around MM3 agrees with the finding of Tang et al. (2019) and can also indicates that magnetic fields in this region may be poloidal.

3.3.4 Structure function and auto-correlation function analysis

We attempted to separate the large-scale and the turbulent scale B-fields in the cloud. In the structure function (SF) method of Hildebrand et al. (2009), the magnetic field consists of large-scale structure, $B_{0}$ and a turbulent component, $\delta B$. The SF analysis provides the variation of angular dispersion of position angles obtained from polarization observations as a function of separation length $\ell$. The turbulent component $\delta B$ reaches the maximum at some scale larger than the turbulent-scale $\delta$. At scales smaller than d (where d is the correlation length scale which characterizes the variation in $B_{0}$ (Hildebrand et al., 2009)), the higher-order terms in a Taylor expansion of regular component $B_{0}$ can be ignored. In case of $\rm\delta<\ell<<\textit{d}$, the angular dispersion function can be written as:

[TABLE]

where ${\langle\Delta\phi^{2}(l)\rangle}_{tot}$ is the square of the total measured dispersion function, where $b^{2}$ is a constant turbulent contribution, $m^{2}l^{2}$ is the contribution from the large-scale field structure, and ${\sigma_{M}}^{2}(l)$ is the contribution of the measured uncertainty. The ratio of the turbulent to large-scale magnetic field components given by

[TABLE]

and $B_{0}$ is estimated as

[TABLE]

$B_{pos}$ is corrected by using a correction factor Q as

[TABLE]

The value of Q is taken as 0.5. The angular dispersion function (ADF) corrected by uncertainty $({\langle\Delta\phi^{2}(l)\rangle}_{tot}-{\sigma_{M}}^{2}(l))$ is shown in Figure 7 plotted as a function of distance measured in polarization map. We followed Hildebrand et al. (2009) and divided data into separate distance bins with separations corresponding to the pixel size. At the scales of $0{\arcsec}-25{\arcsec}$, the ADF increases steeply, probably due to the contribution from the turbulent field. After $25{\arcsec}$ length, the dispersion function increases with shallower slope which may be a contribution from the large-scale regular magnetic fields. It reaches the maximum at $\sim 90{\arcsec}$, the maximum ADF value seen here is less than $52{{}^{\circ}}$ the one expected for random field structure (Poidevin et al., 2010). The structure function is fitted over $25{\arcsec}<l<90{\arcsec}$. The calculated parameters are given in Table 2.

The autocorrelation function (ACF) method (Houde et al., 2009) is the expansion of structure function analysis with the inclusion of the effects of signal integration along the line-of-sight and within the beam. The ADF by Houde et al. (2009) is written as

[TABLE]

where $\Delta\phi(l)$ is the difference between the position angles at a separation of $\ell$, W is the beam radius (6${\arcsec}$ in case of JCMT which is the FWHM beam divided by $\sqrt{8ln2}$), ${a_{2}}^{\prime}$ is the slope of second-order term in the Taylor expansion, and $\delta$ is the turbulent correlation length. N is the number of turbulent cells in telescope beam which is given by

[TABLE]

where $\Delta^{\prime}$ is the effective thickness of the cloud derived from the distance corresponding to the half-maximum of polarized flux of the cloud (Houde et al., 2009). The ordered magnetic field strength can be estimated using

[TABLE]

The upper panel in Figure 8 shows the ADF of polarization segments in G34 “C” region. The lower panel of the figure shows the correlated component of the dispersion function. The function is fitted at $l<90{\arcsec}$ distance. The reduced $\chi^{2}$ of the fitting is 5.3. The turbulent correlation length $\delta$ is $\rm 7.4\pm 0.9{\arcsec}$ (0.13$\pm$0.02 pc). As mentioned above, the turbulent correlation length characterizes the turbulent component of magnetic fields. This is typical the size of a turbulent magnetized cell. Some previous studies have reported values of the turbulent correlation lengths as $\sim$16 mpc and $\sim$10 mpc in the high-mass star forming regions OMC1 (Houde et al., 2009) and Orion KL (Houde et al., 2011). In the starless core Oph-C, Liu et al. (2019) reported correlation length of $\sim$4.3 mpc. All these regions are much closer compared to G34. The turbulent correlation length in G34 is larger compared to nearby regions due to insufficient power to resolve it at 3.7 kpc. The number of turbulent cells in G34 is derived as 5.5$\pm$0.3. Other calculated parameters are given in Table 2. The uncertainties in derived parameters are statistical uncertainties from the dispersion function method. The uncertainty in $B_{pos}$ is taken as a factor of two, as seen in other measurements in several studies. We did not perform this analysis on the northern and southern regions as we do not have enough vectors (20-25 vectors only) for dispersion function analysis and the ADF is too scattered to fit the function. A detailed investigation of change in estimated parameters of SF and ACF analysis on correction with and without beam integration methods can be found in Liu et al. (2019).

A detailed comparison of magnetic fields, gravity, and turbulence on filament to core scale in G34 has been presented by Tang et al. (2019) at 350 $\micron$. In this work, we are comparing these quantities in three different regions of the filament using our POL-2 measurements at 850 $\micron$.

3.4 Comparison to other studies

There have been several attempts to investigate the B-fields in G34 filament in various wavelengths using dust and line emission polarization measurements. Figure 9 shows the field morphologies mapped by JCMT/POL-2 (this work), TADPOL/CARMA333Combined Array for Research in Millimeter-wave Astronomy (Hull et al., 2014), and the SMA444Submillimeter Array (Zhang et al., 2014) observations towards the G34 center containing MM1 and MM2. The field orientation seems to be similar from large to small scales when seen from POL-2 (red vectors) and TADPOL observations (green vectors). But the field geometry changes on even smaller scales seen with the SMA (white vectors). The difference in the POL-2 and SMA field geometries can be seen in zoomed-in lower right panel of Figure 9 where white line segments are misaligned with red lines. A quantitative comparison using histograms of B-field position angles from JCMT, CSO, SMA, and CARMA is shown in Figure 10. The details of these other investigations of B-fields in G34 are given below.

Tang et al. (2019) studied the details of magnetic fields in the regions of G34 containing MM1, MM2, and MM3 using high-resolution (i.e. 10$\arcsec$) 350 $\mu$m CSO/SHARP polarization observations and kinematics using $\rm N_{2}H^{+}$(1-0) line observations. The B-field orientation found perpendicular to the main axis of filament as also seen in this work suggests that field lines are guiding material towards the filament. They found a close alignment between local velocity gradients derived from $\rm N_{2}H^{+}$(1-0) line and local B-field orientation. Since our 850 $\mu$m polarization measurements are consistent with 350 $\mu$m polarization results of Tang et al. (2019), we expect the similar correlation of local velocity gradients and B-field lines at 850 $\mu$m. This kind of correlation suggests a coupling of B-fields and gas motion in G34 filament. Tang et al. (2019) also propose varying relative importances of B-fields, gravity, and turbulence in MM1/MM2 and MM3 resulting in different patterns of small scale fragmentation in the clumps at 0.2 pc scale. The clump containing MM1 shows no fragmentation at all. They found that clumps containing MM2 show an aligned fragmentation and the other clumps with MM3 show a clustered fragmentation. We refer to Tang et al. (2019) for detailed explanation of these findings.

Cortes et al. (2008) presented interferometric observations of polarized continuum in 3 mm wavelength (with 16${\arcsec}$ resolution) and line emission using CO (J=1-0) from G34 filament using BIMA555Berkeley-Illinois-Maryland Association array. They found a very uniform polarization pattern from both dust and line emission as seen in present work in 850 $\mu$m and at 350 $\mu$m by Tang et al. (2019). This is a remarkable consistency of polarization measurements in different wavelengths tracing different dust grains.

Hull et al. (2014) have studied the B-fields in G34 central region using $\lambda$ 1.3mm TADPOL/CARMA observations of dust polarization with 2.5${\arcsec}$ resolution. The observations from the present work at 850 $\mu$m, from the CSO at 350 $\mu$m and BIMA at 3 mm wavelength show uniform and ordered field geometry in G34 central region but the results of Hull et al. (2014) reveal a much more complex polarization pattern with a dragged B-field geometry. They even see a hint of an hourglass morphology in the densest part of the core.

Zhang et al. (2014) investigated the small scale ($\leq 0.1$ pc) B-field structure in G34 center high-density region using SMA at 870 $\mu$m wavelength with 1.5${\arcsec}$ resolution. Their findings also suggest that the magnetic fields are roughly perpendicular to the major axis of the filament and consistent with those of Hull et al. (2014) in MM1. The SMA polarization measurements are uniform but deviate from our 850 $\mu$m B-field orientations.

The zoomed-in panels of Figure 9 show a deviation in field lines. The sub-parsec scale fields are misaligned and even become perpendicular to the large-scale field lines. Similarly in the MM2 region, sub-parsec field lines probed with SMA observations are almost, if not exactly, perpendicular to the large-scale fields seen with JCMT. The field might be strong enough on the clump scale to guide the material along the field lines which eventually get concentrated into cores. The concentration can pinch the B-field lines inside the cores, but does not necessarily lead to complete misalignment with the large-scale field lines. This may be a potential explanation of the change in field geometry from large clump to small core scales. The MHD simulations by Li & Klein (2019) also revealed the deviation of core scale magnetic fields from large-scale average field orientations with a deviation as strong as 90${}^{{}^{\circ}}$. They suggested that change may be caused by the gravitational collapse, enhanced turbulence, and the gas flow along the cloud’s long axis.

This can also be explained with numerical simulations (Ostriker et al., 2001; Nakamura & Li, 2008; Van Loo et al., 2014) showing less disturbed and organized field geometry when B-fields are stronger i.e. $\beta=(P_{th}/P_{B})<<1$, where $\beta$ is the ratio of thermal pressure ($P_{th}$) to magnetic pressure ($P_{B}$). To estimate the $\beta$ values from our observations in G34, we calculated the magnetic pressure ($P_{B}=B^{2}/8\pi$) and thermal pressure ($P_{th}=nkT$), where B is plane-of-the-sky B-field strength, $n$ is the volume density, and T is the gas temperature (Dirienzo et al., 2015). The values of $\beta$ are found to be 0.1, 0.2, and 0.6 in the north, central, and southern regions, respectively. Wareing et al. (2016) investigated 3D MHD simulations to understand the formation of clumps and filaments and to determine the driving processes responsible for filament formation and fragmentation. They explored the range of magnetic field strengths in clouds with $\beta$ varying from 0.1 to 1.0. They found that with no magnetic fields, clumps are found to be forming within the cloud whereas in the case of strong magnetic fields ($\beta$=0.1) these clumps start appearing as filaments. Our findings of $\beta$ values in G34 are consistent to these simulations and suggest that B-fields are playing important role in the formation of this filament.

Outflow patterns of in MM1, MM2, and MM3 of G34 are plotted with SMA polarization measurements in figure 1 of Zhang et al. (2014). The outflows in MM1 are compact and mostly aligned with the small scale B-fields (as seen in their figure) but the outflows in MM2 are highly complex with red and blue shifted lobes overlapping each other. Therefore, it is hard to check for any correlation in B-fields and outflows in the core MM2. Sanhueza et al. (2010) reported the discovery of outflows in MM3 using CO (J=3-2) line observations. The outflow mass and kinetic energy associated to outflows in MM3 suggest a high-intermediate mass star embedded in the core. The outflow orientation is not indicated in their work so it is not possible to relate the field orientations and outflow direction in MM3. The authors also report a possible association of outflows with the core MM4 in G34 central region. The highly ordered and the uniform field geometry of G34 seen in the above mentioned studies including the present work, suggests that feedbacks from these detected outflows associated to MM1, MM2, MM3, and MM4 are not significantly affecting the field geometry of the region. This may be further investigated on much smaller scales using ALMA polarization capabilities.

Among the cores embedded in G34, MM2 has an associated UCHII region (Shepherd et al., 2004). We did not see any prominent change in the B-field lines due to the compression by HII region in this core MM2 as seen by Liu et al. (2018b) in an actively high-mass star-forming region G9.62+0.19. The B-field strength in G34 is found to be less strong than in G9.62+0.19. To further investigate the effect of UCHII regions on B-fields in G34-MM2, we have to probe the fields and kinematics at much smaller scales using ALMA observations as done by Dall’Olio et al. (2019) in G9.62+0.19.

4 Conclusion

We present the plane-of-sky projected magnetic field in G34, obtained using 850$~{}\mu$m dust polarization observations. We investigated the relative importance of gravity, turbulence, and magnetic fields in G34 at sub-parsec scales. The main findings of the study are as follows:

The overall B-field structure in G34 is ordered and perpendicular to the long axis of the filament. The small-scale field geometry is found connected to the large-scale field lines seen with Planck dust polarization observations. The observed aligned B-fields in G34 are consistent with theoretical models suggesting that B-fields play an important role in guiding the contraction of the cloud driven by gravity. 2. 2.

Our measurements of field geometry in G34 using JCMT 850 $\mu$m wavelength are found consistent with previous studies which inferred field morphology at 350 $\mu$m (CSO) and 3 mm (BIMA) wavelengths. However, there is some deviation in the field lines seen at core scale at 870 $\mu$m (SMA) and 1.3 mm (CARMA) wavelengths. 3. 3.

The present study, combined with several similar studies of other IRDCs, suggests that field lines are mostly perpendicular to the filament major axes but change direction at sub-parsec scales in embedded cores which may be caused by relatively different roles of gravity and B-fields than that on clump scale. 4. 4.

We used an updated form of the Davis-Chandrasekhar-Fermi relation to estimate a plane-of-sky magnetic field strength of 470$\pm$190 $\mu$G, 100$\pm$40 $\mu$G, and 60$\pm$30 $\mu$G in the central, northern, and southern regions of G34, respectively. Our results are consistent with those found in several other observations of IRDCs and behavior predicted by theoretical simulations. 5. 5.

From the estimation of mass-to-flux ratio, G34 filament is found to be marginally critical with a criticality parameter $\rm\lambda_{corr}$ of 0.8$\pm$0.4, 1.1$\pm$0.8, and 0.9$\pm$0.5 in the central, northern, and the southern regions, respectively. 6. 6.

The values of Alfvénic Mach number in all three regions correspond to a sub-Alfvénic nature of turbulence in G34 filament.

The authors thank referee for a constructive report which considerably helped to improve the content of the manuscript. AS and B-GA are supported by National Science Foundation Grant-1715876. TL is supported by KASI and EACOA fellowships. AS acknowledges the support from Korea Astronomy & Space Science Institute (KASI) for our JCMT observations. MJ, ERM and VMP are supported by Academy of Finland Grant-285769. VMP is also supported by the European Research Council, Advanced Grant No. 320773, and the Spanish MINECO under projects MDM-2014-0369 and AYA2017-88754-P. CWL was supported by National Research Foundation of Korea (NRF-2019R1A2C1010851). AS thanks Simon Coud$\acute{e}$ for a useful discussion during the revision and Piyush Bhardwaj for a critical reading of the draft. This work was carried out in part at the Jet Propulsion Laboratory, operated for NASA by the CalTech. DJ is supported by National Research Council Canada and an NSERC Discovery Grant. WK was supported by Basic Science Research Program through the National Research Foundation of Korea (NRF-2016R1C1B2013642). JCMT is operated by the East Asian Observatory on behalf of National Astronomical Observatory of Japan; Academia Sinica Institute of Astronomy and Astrophysics; the Korea Astronomy and Space Science Institute; the Operation, Maintenance and Upgrading Fund for Astronomical Telescopes and Facility Instruments, budgeted from the Ministry of Finance of China and administrated by the Chinese Academy of Sciences and, the National Key R&D Program of China (No. 2017YFA0402700).

Facility: James Clerk Maxwell telescope (JCMT) Softwares: Starlink (Currie et al., 2014), Astropy (Astropy Collaboration et al., 2013).