© ESO 2021

1 Introduction

The role of magnetic fields in the formation and evolution of molecular clouds (MCs) is the object of active research, and represents one of the keys to understanding the earliest phases of star formation (see e.g., the reviews by Li et al. 2014; Pattle & Fissel 2019; Hennebelle & Inutsuka 2019). It has been established that MCs have a filamentary nature (e.g., Ungerechts & Thaddeus 1987; Bally et al. 1987; Ward-Thompson et al. 2010; Arzoumanian et al. 2011), and indeed filaments are observed to develop in numerical simulations (e.g., de Avillez & Breitschwerdt 2005; Padoan et al. 2007; Federrath & Klessen 2013). However, the exact mechanisms leading to the formation and evolution of such elongated structures, and in particular the relative contributions of turbulence, gravity and magnetic fields, are not yet fully understood.

The study of the relative orientation between the filamentary column density structures and the projection of the magnetic field on the plane of the sky can provide insights into the formation mechanisms of the filaments and, ultimately, into the evolutionary stage of the star-forming region hosting them. Different methods have been used so far to identify linear structures in astronomical data and quantify their orientation. For instance, the DisPerSE method, originally developed by Sousbie (2011) to recover density skeletons in cosmic web data, was used by Peretto et al. (2012) and Palmeirim et al. (2013) to study the formation of filamentary structures in the Pipe nebula and in the Taurus MC, respectively. The inertia matrix has been used to investigate the physical origin of filaments in magnetized media (Hennebelle 2013), while the Hessian matrix has been used to analyze the relative orientation between density structures and the magnetic field traced by dust polarization in the all-sky Planck data (Planck Collaboration XXXII 2016).

In this paper we focus on the comparison between two other methods that have been widely used to investigate the relative orientation between column density structures and the magnetic field. The first method is the Rolling Hough Transform (RHT; Clark et al. 2014), and the second is the gradient technique (GRAD; Soler et al. 2013; Planck Collaboration XXXV 2016). Both methods recover the same quantities, namely, the angle between the direction of the density structures and the direction of the magnetic field projected on the plane of the sky. While this latter is inferred in both examples from dust polarization, the former is derived in two different ways.

The RHT method has been used for instance by Koch & Rosolowsky (2015), Malinen et al. (2016), Panopoulou et al. (2016), and Alina et al. (2019) as well as in numerical simulations by Inoue & Inutsuka (2016). GRAD has been employed to analyze avariety of observational data (Planck Collaboration XXXV 2016; Soler et al. 2017; Jow et al. 2018; Soler 2019) and in simulations (Chen et al. 2016; Soler & Hennebelle 2017).

One of the main results that emerged from both observational and theoretical studies of the alignment between density structures and magnetic field is that the relative orientation has a bimodal distribution, regardless of the method used to determine it. Low-density column density structures are preferentially aligned parallel to the magnetic field, while high-density regions are mostly orthogonal to the field or do not show a favored orientation.

In terms of this result, RHT and GRAD give results which are consistent in most of the cases. In consideration of the intrinsic differences between the two methods, this could be seen as a bit surprising. To our best knowledge, however, no systematic comparison between the two methods has been performed. Such a comparison is therefore the focus of our paper, to elucidate which factors determine the differences and similarities between the output from the two techniques. This is necessary to properly interpret the data and investigate the physical processes underpinning the observed interplay between gas column density structures and magnetic field.

We perform our analysis on density structures generated via state-of-the-art magnetohydrodynamic (MHD) simulations and traced via their dust emission at 353 GHz (850 μm). Because dust grains have a preferred orientation relative to the magnetic field, their emission is polarized. Dust is therefore a tracer of gas column density (through its emission) and of magnetic fields (through its polarized emission).

This paper is organized as follows. In Sect. 2 and Sect. 3 we describe the MHD simulation and the radiative transfer calculations, respectively. Section 4 outlines the methods used to quantify the orientation of the local magnetic field with respect to the density structures. Our results are illustrated in Sect. 5 and discussed in Sect. 6 where we also present our conclusions.

2 MHD simulation

We use an MHD simulation describing the ISM in a region with a volume of (250 pc)³, where the turbulence is driven by supernova (SN) explosions. The simulation, carried out with the Ramses adaptive-mesh-refinement (AMR) code (Teyssier 2002), was first presented in Padoan et al. (2016a), and further analyzed in Pan et al. (2016); Padoan et al. (2016b); we only describe it briefly here and refer the reader to Padoan et al. (2016b) for further details. The MHD runs used a minimum cell size (maximum spatial resolution) of d x = 0.24 pc, periodic boundary conditions, a mean density of 5 cm⁻³ (corresponding to a total mass of 1.9 × 10⁶ M_⊙) and a SN rate of 6.25 Myr⁻¹, with the SN explosions randomly distributed in space and time. Individual SN explosions are implemented with an instantaneous addition of 10⁵¹ erg of thermal energy and 15 M_⊙ of gas, distributed according to an exponential profile on a spherical region of radius r_SN = 3dx = 0.73 pc, which guarantees numerical convergence of the SN remnant evolution (Kim & Ostriker 2015). The energy equation includes the pd V work, the thermal energy from the SN explosions, photoelectric heating up to a critical density of 200 cm⁻³, and parametrized cooling functions from Gnedin & Hollon (2012).

The simulation starts with zero velocity, uniform density, n_H,0 = 5 cm⁻³, uniform magnetic field, and uniform temperature, T₀ = 10⁴ K. The SN-driven turbulence brings the mean thermal, magnetic and kinetic energy to an approximate steady-state, where the amplified magnetic field has an rms value of 7.2 μG and an average of |B| of 6.0 μG, consistent with the observations. The simulation was integrated for 45 Myr without self-gravity and then continued with self-gravity for 11 Myr. The snapshot used in this work is taken from the end of this second part of the simulation including self-gravity.

The analysis of this simulation has shown that SN-driven turbulence can explain the formation and evolution of MCs, their internal turbulence, their lifetimes, their mass and size distribution. With a higher-resolution continuation of the simulation, including sink particles, it was shown that SN-driven turbulence may also explain the low value of the star-formation rate in MCs and the even lower global star-formation rate in the Galaxy (Padoan et al. 2017). Because of their realistic properties, MCs selected from this simulation are an ideal tool to study properties that are difficult to infer directly from the observations, such as the spatial structure of the magnetic field, and its relation to the morphology of density structures, such as filaments and clumps.

The AMR method provides a high spatial resolution in regions of interest in the large computational volume, namely MCs and their substructures. The resolution of our simulations (0.24 pc/cell) is not sufficient to resolve physical processes at the spatial scales probed by the Planck data (0.2 pc for a cloud at 150 pc). However, the resolution is sufficient to resolve the structure at larger distances. For the purpose of the analysis of this work, we should stress the caveat that low-density regions are not spatially refined in these simulations, so their spatial resolution corresponds to that of the 128³ root grid, that is ~2 pc. This resolution is clearly insufficient to capture the filamentary structure often observed in low-density regions surrounding MCs. Thus, the study of the relative orientation of the magnetic and thin, low-density filaments is beyond the scope of this work. However, the simulations allow us to study the orientation of the magnetic field with respect to dense structures, and to compare different methods of pursuing such a study.

3 Radiative transfer calculations

We use radiative transfer (RT) modeling to predict the total dust emission and the polarization observed from the model clouds. The calculations were performed with the SOC program, which is a Monte Carlo radiative transfer program for the calculations of dust emission and scattering. The SOC program has been used in some previous publications (Gordon et al. 2017; Juvela et al. 2018)and is described in detail in Juvela (2019).

The SOC program could directly use the octtree grid of the MHD simulations. However, in this paper the cloud data were resampled onto a regular 512³ grid with a 0.49 pc cell size. For example at a distance of 550 pc, this would thus correspond to a resolution of 3.1′.

We adopted the dust model from Compiègne et al. (2011). The model cloud was illuminated from the outside by an isotropic radiation field with the intensities given in Mathis et al. (1983). We did not include any discreteradiation sources inside the model volume. The average visual extinction through the model is A_V = 1.7 mag but, because of the inhomogeneity of the density field, the large-scale radiation field remains relatively constant throughout the volume. Therefore, the individual dense regions within the models are illuminated by an intensity that is only slightly weaker than the external field.

RT calculations were used to solve the dust temperature in each model cell and, based on this information, to produce surface brightness images at 353 GHz. The maps were calculated for observers in the directions of the main axes. In this paper, we do not model in detail the grain alignment processes and instead assume a constant efficiency of dust alignment, with an intrinsic polarization fraction of 20%. Thus, in addition to the total intensity I, we calculated maps for the Stokes parameters Q and U from the following equations:

where the Rayleigh polarization reduction factor R is constant. Here j_ν is the emissivity that depends on the local density and dust temperature. At 353 GHz, the optical depth between the point of emission and the observer, τ_ν(s), is very small and the exponential terms are practically equal to 1. The equations include two angles that depend on the direction of the magnetic field B: Ψ is the position angle of B projected onto the plane of the sky (with respect to a chosen reference direction) and γ is the angle between B and the plane of the sky.

Given maps of Q and U, one can estimate the polarization angle χ from (3)

and the polarization fraction p from (4)

Our model is by definition free of noise. One may have noise in dust temperature but once those values are fixed, the resulting (I, Q, U) values are internally consistent to a very high precision. Furthermore, because the Monte Carlo errors are contained in the emissivity j_ν, these leave both χ and p essentially unchanged, apart from small variations in the relative weighting of different cells along the line of sight.

The calculated surface brightness maps at 353 GHz are shown in Fig. 1. Each panel refers to one observing direction (x-left, y-center, z-right). We consider specifically the snapshot labeled as “315”, taken 11 Myr after self-gravity was switched on in the simulation. This ensures that the initial density field had time to collapse into denser structures. We use physical units (parsecs) for the lengths, so that it is not necessary to assume a distance to the clouds. The size of each pixel is 0.49 pc, and our synthetic observations correspond to the maps that have been convolved with a Gaussian beam with full width at half maximum (FWHM) equal to 2.5 pix = 1.22 pc. This ensures that the data are well sampled and that the end result is independent of the details of the original pixelization. The surface brightness maps are overlaid with a red “fingerprint” pattern. The pattern itself indicates the direction of the magnetic field, observed from the Stokes parameters (Eqs. (1) and (2)), while the intensity of the red color increases with the polarization fraction.

The lack of polarization when the cloud is viewed along the x-axis reflects the initial orientation of the magnetic field along this direction in the MHD simulation. Because magnetic flux is conserved in the computational volume, the mean magnetic field (averaged over the whole volume) remains in the x direction at any time. While observing in the y and z directions, we see a clear large-scale orientation of the magnetic field in the x-direction on the more diffuse regions, like the mean field in the simulation, while the polarization fraction is strongly reduced in the denser structures. To perform our analysis, we selected two regions whose physical size is comparable to those of the MCs analyzed in Planck Collaboration XXXV (2016). Our selected regions are highlighted in Fig. 1: Region 1 is almost entirely occupied by a dense and thick filamentary structure with some branches, while Region 2 includes a thinner filament with denser knots surrounded by a more tenuous medium. The blue dashed box identifies a low-brightness area selected as a reference region (see Sect. 4.2).

Fig. 1 Surface brightness maps at 353 GHz (850 μm – shades of gray) resulting from radiative transfer calculations performed on density structure maps from MHD simulations. Each frame shows the maps corresponding to one observing direction (x-left, y-center, z-right). The number “315” labeling each panel refers to a snapshot taken 11 Myr after self-gravity has been switched on in the simulation. The maps have been convolved using a Gaussian beam with FWHM = 2.5 pix = 1.22 pc. The “fingerprint” pattern shows the direction of the magnetic field derived from the Stokes parameter, and the intensity of the red color increases with the polarization fraction. The red boxes indicate the regions where individual analysis was performed while the blue dashed box identifies a low surface brightness area selected as a reference region (see Sect. 4.2 for details).

4 Relative orientation between magnetic field and density structures

4.1 The Rolling Hough Transform method

The RHT technique has been introduced by Clark et al. (2014) and it is based on the Hough transform from Hough (1962) as implemented by Duda & Hart (1972). The Hough transform was originally devised to track the motion of particles in high-energy physics experiments. Clark et al. (2014) developed a rolling version of the Hough transform to quantify the linearity and spatial coherence of density structures in HI data. The basic steps of the procedure can be summarized as follows (see Fig. 2 in Clark et al. 2014). The image to analyze is high-pass filtered with a circular filter of diameter D_K to enhance the filamentary structure, and the image is then thresholded to form a bitmask. A circular patch of diameter D_W is centered on each pixel and rolls across the data. The pixels inside the extraction window that pass an intensity threshold Z and that are aligned with a local direction θ are counted. The adopted values for the RHT parameters are the following: D_K = 10 pix = 4.9 pc, D_W = 25 pix = 12.3 pc, Z = 70%. The output of the procedure is the RHT intensity R(θ, x, y), where (x, y) are the coordinates of each pixel.

The RHT output provides the intensity as a function of orientation for every image pixel. For our purposes, we need to obtain a single estimate for the orientation of column density structures at the location of each pixel. To do so, we use the point estimator described in Clark et al. (2014), which quantifies with a single angle the direction of a given region. In our case, the point estimator for each pixel is the RHT angle expectation value , given by the equivalent value on the interval θ ∈ [0, 180°) of the angle (see Eqs. (7) and (8) in Clark et al. 2014): (5) (6)

The angle is elected asbeing the local filament direction. With this, for each pixel we can derive directly the angle ϕ_RHT between the direction of the magnetic field and the filament: (7)

where is the position angle of B on the plane of the sky derived from the synthetic observations (using the Stokes parameters). The angle ϕ_RHT is defined in the interval [−90°, 90°]. A value close to zero indicates that the magnetic field and the filamentary structures have a similar orientation.

To visualize the linear structures identified by the RHT technique, we use the backprojection R(x, y) from Clark et al. (2014). This is calculated via integration of the RHT intensity R(θ, x, y) over the angle θ: (8)

4.2 The GRAD method

The GRAD method consists in finding the local orientation of the magnetic field with respect to the column-density isocontour. This is equivalent to calculating the local orientation of the unit polarization pseudo-vector (orthogonal to the magnetic field) with respect to the local gradient, which is by definition perpendicular to the column density isocontour. In practice, the gradient derives directly from the finite differentiation of the intensity field on a mesh, using a Gaussian convolution kernel. For our implementation, we use a kernel with FWHM = 1.2 pc.

The unnormalized gradient of the column density N_H is given by (9)

with a = (i, j). The gradient is unnormalized as there should be a scaling factor to account for the grid mesh spacing in the finite differentiation, 0.49 pc in our case, However, such a scaling is unimportant as we are only interested in the unit vector giving the direction (10)

The polarization angle χ (Eq. (3)) characterizes the direction of (11)

from which we derive the unit pseudo-vector : (12)

From the above equations, we obtain the relative angle between the normal to the column-density isocontour and the polarization direction (see Eq. (2) in Planck Collaboration XXXV 2016): (13)

In the above, we used û_z, the directly oriented unit vector normal to the plane of observation. From Eq. (13) it is clear that the absolute amplitudes of the two vectors do not matter as only the value of the ratio of the amplitudes is used.

The pixels over which we perform our analysis are selected applying the following basic criterion on the above-mentioned set of coordinates : (14)

with being the set of pixels included in our selected Region 1 and Region 2 (Fig. 1). Following Planck Collaboration XXXV (2016), we define the threshold for the gradient as (15)

with being the reference region and its cardinality. For this we chose a smooth and low-brightness area in the z-axis snapshot where we scaled down the intensity by a factor of ten to get an average column density of ~2 × 10¹⁹ cm⁻². Such a value is smaller than the one of our lowest-column-density bin (see Sect. 4.3). Our reference region is identified bythe blue dashed rectangle in Fig. 1.

4.3 Histograms of relative orientations and the shape parameter ξ

We use the calculated values of ϕ_RHT(a) and ϕ_grad(a) to build the histograms of relative orientations (HROs) for the angles between the magnetic field and the density structures (e.g., Soler et al. 2013; Planck Collaboration XXXV 2016; Planck Collaboration XXXVIII 2016; Malinen et al. 2016) for each of the 15 bins in which we grouped our column density values. The bins have been designed to contain the same number of pixels. Planck Collaboration XXXVIII (2016) and Malinen et al. (2016) use the same term to refer to the histograms, regardless of the method used to determine the orientation of the structures. We adopt the same naming convention and use the term “HROs” to indicate the histograms, irrespective of the method used to calculate the orientations.

The HROs are here presented (Sect. 5.2) as a probability density P(ϕ) obtained via normalization of Ã_i based on a set of a coordinates: (16)

where (17)

with h_i = − 90° + 180° × i∕12 and i ∈ [0,12]. A histogram peaking at around ϕ = 0° means that the magnetic field is preferentially oriented parallel to the density structures, while a histogram having a minimum at the same angle indicates that the magnetic field is mostly orthogonal to the density structures.

From the HROs we compute the histogram shape parameter ξ (Planck Collaboration XXXV 2016; Soler et al. 2017) to quantify the variations occurring in the HROs as a function of column density. We calculate A₀ and A₉₀ as the area of the central part of the histogram (−22.°5 < ϕ < 22.°5) and the area of the far ends of the histogram (−90.°0 < ϕ < −67.°5 and 67.°5 < ϕ < 90.°0), respectively. Finally, we define ξ, as in Planck Collaboration XXXV (2016): (18)

To get an approximation of the error bar due to Poisson count we use the Gaussian limit of the Poisson law. For a Poisson intensity λ the standard deviation is . The Gaussian error on ξ is (19)

For each column density bin, ξ > 0 means that the magnetic field and the density structures are mostly parallel (concave histogram), ξ < 0 indicates a mostly perpendicular orientation (convex histogram) and ξ ≈ 0 reveals that there is no preferred orientation (flat histogram).

4.4 The projected Rayleigh statistic

Besides the HRO statistics described above, we also applied to our simulated data the so-called projected Rayleigh statistic (PRS; Jow et al. 2018). The PRS is based on the classic Rayleigh test typically used in circular statistics to determine whether theangles of a specific set are uniformly distributed (e.g., Batschelet 1981; Glimm 1996; Mardia & Jupp 1999) and can be considered a specific case of the V statistics (Durand & Greenwood 1958; Mardia & Jupp 1999).The V statistics allows us to test for uniformity against a specific mean direction characterized by an angle θ. The V statistics for the case θ = 0, corresponding to a parallel orientation, has been renamed PRS. Following Jow et al. (2018) we use our angles ϕ(a) (from both RHT and GRAD) to define the PRS for the same set . The parameter used to quantify the relative orientations is Z_Jow defined as (20)

As from Jow et al. (2018), the typical uncertainty on this quantity is given by: (21)

As for ξ, positive values of Z_Jow indicate a parallel orientation, negative values a perpendicular orientation and values compatible with zero show that there is no preferred orientation.

5 Results

5.1 Identified structures

The structures identified by RHT (visualized by the backprojection R(x, y)) and by GRAD are shown in Figs. 2 and 3 for Region 1 and Region 2, respectively. The top-left panel shows the filamentary structures identified by RHT applied to our surface brightness maps after convolution with a Gaussian beam with FWHM = 2.5 pix = 1.22 pc. The set of parameters for RHT is the following: (D_K, D_W, Z) = (4.9 pc, 12.3 pc, 70%). We identify this case as RHT-H.

The top-right panel shows the structures selected by GRAD applied using the same Gaussian convolution kernel with FWHM = 1.22 pc, imposing that the module of the gradient must be greater than the average over a reference smooth region (Sect. 4.2). We refer to this case as GRAD-H.

To perform a meaningful comparison between the two methods, however, we need to make sure that they are both working at the same resolution, that is, on the same scale. This means that the gradient must be calculated over a region that matches the size of the extraction region in RHT, where this region has diameter D_W. To derive the size of the new Gaussian convolution kernel required to calculate the gradient, we proceeded as follows.

The relation between the FWHM and the standard deviation, σ, of the Gaussian distribution is given by the following expression: (22)

We make the reasonable assumption that the extraction performed in the RHT method resembles a top-hat filtering (compare with Fig. 2 in Clark et al. 2014). For Gaussian and top-hat to behave in the same way at large scales, we impose that σ = R∕2 where R is the radius of the top-hat filter. In our case R = D_W∕2, therefore: (23)

For D_W = 25 pix = 12.3 pc we obtain FWHM = 14.7 pix = 7.2 pc. We refer to this case as GRAD-L, to indicate that GRAD is applied with a lower resolution. The identified structures are shown in the bottom-left panel of Figs. 2 and 3. Because Region 1 and Region 2 are part of a much larger map, there are no edge effects and the convolution is well defined over the whole area shown in the figures.

The comparison between the RHT-H and GRAD-H/GRAD-L selections reveals that RHT clearly picks narrow and elongated features while GRAD selects extended regions. This can be ascribed to the fact that in the GRAD technique there is no interdependence between the detected pixels, combined with the selected degree of thresholding. Indeed, GRAD is sensitive to contours, which makes it more adapted to perform orientation analysis of extended regions. In the GRAD-H case the excluded areas are typically small and scattered, while in the GRAD-L case they appear bigger and locally concentrated. This is due to the larger size of the smoothing kernel used to calculate the gradient in the GRAD-L case.

In the RHT-H case, the selected structures cover 45.7 (Reg. 1) and 37.7% (Reg. 2) of the total area of the regions, while the GRAD method selects much more extended areas. In the GRAD-H case the coverage is 99.5 (Reg. 1) and 95.8% (Reg. 2), and in the GRAD-L case the values are 96.5 (Reg. 1) and 83.6% (Reg. 2). This is due to the fact that the threshold for the gradient is low with respect to the value of the gradient in our considered regions, therefore most of the pixels pass the selection. We repeated our analysis adopting a reference region with higher average column density (3 × 10²⁰ cm⁻²) which provides a higher threshold for the gradient (see the Appendix for the corresponding figures). As expected, the selected regions occupy a much smaller fraction of the total area: 79.7% (Reg. 1) and 43.9% (Reg. 2) in the GRAD-H case, and 51.7% (Reg. 1) and 19.6% (Reg. 2) in the GRAD-L case. However, we find that the behaviors of the HROs and of the ξ and Z_Jow parameters remain practically unchanged.

Fig. 2 Comparison between the selected structures in Region 1, obtained from the RHT and GRAD methods. The recovered structures, indicated in blue, are overlaid on the brightness map (shades of gray). The RHT-H case (top-left) illustrates the results from RHT applied to surface brightness maps convolved using a Gaussian beam with FWHM = 1.2 pc. The GRAD-H case (top-right) results from GRAD using a convolution kernel of the same size, with the constraint that the module of the gradient must be greater than the average over a reference (smooth) region. The same constraint has been imposed in the GRAD-L case (bottom-left) with the difference that GRAD is applied using a Gaussian convolution kernel with FWHM = 7.2 pc. The bottom-right panel shows the superposition of the results for the RHT-H case (green) and the GRAD-L case (red). The yellow color indicates the regions in common between the two cases.

Fig. 3 Same as Fig. 2 but for Region 2.

Fig. 4 Histogram of relative orientation between the RHT angle expectation value (RHT-H case) and the direction angle derived from the GRAD method (GRAD-L case) for Region 1 (top) and Region 2 (bottom). The angles were calculated for the pixels that are common to both methods (yellow areas in Figs. 2 and 3). The curves show the HROs for three column-density bins: the highest and lowest ones (red and black, respectively) and the intermediate bin between these two (blue). A ± 1 σ error is represented by the width of the shaded areas.

5.2 Relative orientations

To make a direct comparison between the orientations of the density structures recovered from RHT and from GRAD, we constructed the HROs (Fig. 4) between the RHT angle expectation value (RHT-H case) and the direction angle derived from the GRAD method (GRAD-L case), and we plotted the value of the relative orientation on top of the intensity map (Fig. 5), for both Region 1 and Region 2. The angles have been calculated for the pixels that are commonto both methods, corresponding to the yellow areas in Figs. 2 and 3.

The curves shown in Fig. 4 represent the HROs for three of the column-density bins which we grouped our pixels into: the highest and lowest ones (red and black, respectively) and the intermediate bin between these two (blue). The width of the shaded areas represents a ± 1σ error. In Region 1 the angle agreement, indicated bypeak at around 0° is modest, while in Region 2 it is practically nonexistent.

The modest angle agreement is confirmed by Fig. 5, where the pixels in common to both methods have been color-coded as a function of the relative orientation between the density structures. While in Region 1 (left panel) there is some balance between green areas (parallel orientation) and magenta areas (orthogonal orientation), in Region 2 the magenta and blue areas dominate.

Figures 6–9 present our results in terms of quantification of the relative orientation between the magnetic field and the column density structures in our simulated regions. In each figure, the left column shows the intensity map of the region under consideration overlaid with the red drapery pattern and with a set of solid and dashed contours. The areas between the solid and dashed contours (or inside them, if only one contour is present) identify three specific column density bins: red is the highest one, black the lowest one and blue the intermediate one. The central column shows the HROs for each of the column density bins described above. The width of the shaded areas represents a ± 1σ error. Finally, the right column presents the parameters ξ and Z_Jow defined in Sect. 4, plotted as a function of the column density, N_H, which has been grouped into 15 bins.

We recovered the relative orientations applying the HRO analysis and evaluating the ξ and Z_Jow parameters for the three structure selections described in Sect. 5.1 and highlighted in blue/purple in Figs. 2 and 3: RHT-H, GRAD-L and GRAD-H. We point out that we use these three terms to indicate both the selected regions and the method/conditions used to calculate the relative angle between density structures and magnetic fields. These can be decoupled, as we show afterward. Specifically, in Fig. 6 (for Region 1) and in Fig. 7 (for Region 2) we have chosen the RHT-H sample and applied to it the RHT-H, GRAD-L and GRAD-H cases. This allows us to evaluate the three cases having for all of them the same bias on the selected regions. The three panels in the first column are equal, to illustrate that we are considering the same sample. The top, middle and bottom rows show the results for the RHT-H, GRAD-L and GRAD-H cases, respectively.

In Region 1 (Fig. 6), the global preferred orientation for the density structures is orthogonal to the magnetic field, as can be deduced from the HROs and the ξ and Z_Jow parameters plots. In these latter, the values of the two parameters remain negative or close to zero for most of the column density bins, revealing the perpendicular privileged orientation. Within each case, the ξ and Z_Jow curves are similar between each other, with the errors on Z_Jow slightly smaller than those on ξ. However, there are some differences between the three considered cases. In the RHT-H case the curves show some ample oscillations in the low-density bins followed by a sharp decrease. In the GRAD-L case, the oscillations are smaller and concentrated in the midle-density bins, the curves have a global descending trend except for the highest-column-density bin. In the GRAD-H case, the curves show a flat behavior with some oscillations, with the highest-column-density bin pointing toward less negative values for both ξ and Z_Jow.

In Region 2 (Fig. 7), the preferred orientation of the density structures is again perpendicular to the direction of the magnetic field. All curves show some moderate oscillations. The global trend for ξ and Z_Jow as a function of N_H is flat in all three cases.

Figures 8 and 9 (for Region 1 and Region 2, respectively) are the analogs of Figs. 6 and 7 but in this case the relative orientations have been derived by applying each method to its own native selected area, that is: RHT-H method on RHT-H selection (top row, unchanged) GRAD-L and GRAD-H cases applied to GRAD-L, and GRAD-H samples, respectively (middle and bottom rows). The use of different samples is illustrated by the three different panels in the left column.

In Region 1, the preferred orientation of the density structures is again perpendicular to the magnetic field. In the GRAD-H case, ξ and Z_Jow show the same flat behavior as a function of the column density recovered in the previous case, but with the amplitude of the oscillations and the error bars strongly reduced. In the GRAD-L case, the two parameters show a tendentially flat behavior with some oscillations, and a value close to zero in the highest-column-density bin.

In Region 2, the GRAD-H case shows again a flat trend with few oscillations. In the GRAD-L case, the behaviors of ξ and Z_Jow show a remarkable difference with respect to the previous cases. The two curves, which are essentially coincident and present a few peaks, start at a positive value for the lowest-column- density bin and show a clear decrease toward the negative region of the plot, although the highest-column-density bin reverts the trend moving toward less negative values. We recover in this case the expected transition of the preferred orientation of the density structures, going from mostly parallel to the magnetic field in tenuous media to mostly perpendicular in denser regions.

Figures 10 and 11 show the surface brightness maps (shades of gray) for Region 1 and Region 2 respectively, overlaid with our analyzed selections, color-coded according to the relative orientation between the magnetic field and the density field. The color goes from blue for fully parallel orientation (relative angle equal to zero) to red for fully perpendicular orientation (relative angle equal to 90°). The top rowshows the results when the different methods are performed on the same sample, selected by RHT. In the bottom row, each case has been applied to its native selection.

In both Region 1 and Region 2, the global shape of the selected structures (those color-coded in blue-red) is very similar, which confirms the fact that the analysis has been indeed performed over the same pixel selection. However, some interesting differences can be noted. For both RHT-H and GRAD-L, the colored regions are fairly smooth and have similar thicknesses, indicating that the two methods have been performed at the same resolution. On the other hand, the colored regions for GRAD-H (right column) exhibit a finer granularity, confirming that the method has worked at a higher resolution than in the other two cases.

In terms of the specific orientations (the color of the regions) all three cases exhibit a similar level of mixing, with blue and red often interwoven in the same region. The kind of recovered orientation, preferentially parallel (blue) or orthogonal (red), is generally consistent in the three considered cases, with some notable exceptions. In Region 1 for instance, the vertical feature in the top-right corner of the maps is mostly red for RHT-H and GRAD-H and mostly blue for GRAD-L. In Region 2, the elliptical ring-shaped feature in the top part of the maps is blue for RHT-H, red for GRAD-L and mixed for GRAD-H.

When each method is applied to its own native pixel selection (bottom row of Figs. 10 and 11), the results are much more diverse. The high-coverage and fine-grain structures analyzed in the GRAD-H case (right panel) show a complex pattern where preferentially parallel and preferentially perpendicular orientations are mixed together at small scales. The RHT-H case (left panel) shows the capability of the method to recover narrow elongated structures. The blue filaments develop in the horizontal direction while the red ones develop in the vertical direction. This is consistent with the fact that in Regions 1 and 2, the large-scale magnetic field is oriented in the x-direction, corresponding to the horizontal direction in the z-axis snapshot where the regions have been identified (right panel in Fig. 1).

The general trend shown by the RHT-H case, with primarily horizontal blue filaments and primarily vertical red filaments, is not shown in the GRAD-L and GRAD-H cases. The red and the blue filaments from these methods can both be either horizontal or vertical, with no clear preference for one orientation versus the other. In the GRAD-L case, the colored patches are both extended and filamentary, showing that the method is able to follow the elongated higher-density structures present in both regions. The parallel (blue) and perpendicular (red) orientations alternate across the regions. In Region 1, a blue vertical filament seems to trace the backbone of the main density structure, while the horizontal branch is traced by a red filament. The main density structure in Region 2 seems also to be traced by blue strips.

Fig. 5 Relative orientation between the RHT-H angle expectation value and the GRAD-L direction angle (magenta-blue-green) overlaid on the surface brightness map (shades of gray) for Region 1 (left) and Region 2 (right). The relative orientation goes from fully parallel (green, 0°) to fully perpendicular (magenta, 90°). As in Fig. 4, the angles have been calculated for the pixels that are in common to the RHT-H and GRAD-L cases.

Fig. 6 Analysis of Region 1 from Fig. 1, where the RHT and GRAD methods have both been applied to the structures selected by RHT (RHT-H case, top-left panel in Fig. 2). The columns show the surface brightness map overlaid with the red fingerprint pattern (left); the histograms of relative orientation (HROs – middle) between the magnetic field direction and the density structures; the ξ and ZJow parameters plotted as a function of column density, NH (right). For both parameters, a positive value indicates a preferred parallel orientation, a negative value a preferred perpendicular orientation and a value compatible with zero indicates that there is no preferred orientation. The areas between the solid and dashed contours on the maps (or inside them, if only one contour is present) highlight the regions corresponding to the column density bins shown in the HRO plots. Red and black identify the highest and lowest bins respectively, while blue indicates the intermediate bin between the two above. In the HRO plots, a ± 1σ error is represented by the width of the shaded areas. The top row shows the results for RHT-H, the middle row for GRAD-L and the bottom row for GRAD-H.

Fig. 7 Same as Fig. 6 but for Region 2.

Fig. 8 Same as Fig. 6 but in this case the RHT and GRAD techniques have been applied each to their own native selections, as shown by the three different panels in the left column.

Fig. 9 Same as Fig. 8 but for Region 2.

Fig. 10 Relative orientation between the magnetic field and the density structures (blue-red) overlaid on the surface brightness map (shades of gray) in Region 1. The relative orientation goes from fully parallel (blue, 0°) to fully perpendicular (red, 90°). The top row shows the results from the analysis performed on the same pixel subsample for the RHT-H case (left), the GRAD-L case (middle) and the GRAD-H case (right). In the bottom row, each method is performed on its own native selection. Because the subsample chosen for the same-pixel analysis is the one from RHT (RHT-H case), the two panels on the left are the same.

Fig. 11 Same as Fig. 10 but for Region 2.

6 Discussion and conclusions

In this paper we have studied the structure of magnetic fields as outlined by dust polarization in synthetic maps of MCs generated via MHD simulations. We focused on the comparison between two methods which are widely used to analyze the relative orientations between the magnetic field, traced by polarized dust emission, and the density structures, for which dust emission properly propagated through radiative transfer is used as a proxy. The techniques that we have investigated are the RHT method and the GRAD method and our goal was to establish under which conditions the comparison between the two methods is meaningful. We considered three specific cases: RHT applied to data convolved with a Gaussian beam with FWHM = 1.22 pc (RHT-H case), GRAD applied using the same convolution kernel (GRAD-H) and, finally, GRAD applied using a convolution beam with FWHM = 7.2 pc (GRAD-L), to ensure that GRAD works on the same scale as RHT.

We selected two regions in our simulated maps, named Region 1 and Region 2, and performed the relative-orientation analysis in the three aforementioned cases. We first considered the same region selection in all cases, and then we applied each method to its own native selection, as it is done when the analysis is performed on observational data. To quantify the behavior of the relative orientation as a function of the column density N_H we used two statistics: the HRO providing the ξ parameter with uncertainty given by the Gaussian limit of the Poisson law, and the PRS giving the parameter Z_Jow with its characteristic error.

Our first finding is that the sample selection in the three cases is very different, both in terms of the number and location of the selected pixels, although the number is of course dependent on the selected thresholds (see below; Appendix A). This implies that each method traces different regions and is sensitive to different scales. When the methods are applied to the same pixel selection, the comparison between the ξ and Z_Jow curves as a function of column density in our three cases shows that there are some noticeable differences in Region 1, especially in terms of the amplitude of the oscillations and the global trend (decreasing in the RHT-H and GRAD-L cases and flat in the GRAD-H case) while the curves are more consistent between each other in Region 2. This implies that, when the analysis is performed on the same data sample, the output in terms of relative orientations between the density structures and the magnetic field mostly depends on the resolution at which the methods are applied and on the characteristics of the analyzed regions (Sect. 5.2). Region 1 and Region 2 are substantially different, with the former dominated by a thick filamentary structure with some branches, and the latter showing a thinner filament dotted with denser clumps.

When each method is applied to its own native pixel selection, the differences are more evident and the role of the specific method used for the analysis becomes much more apparent. In particular, in both regions, the ξ and Z_Jow curves in the RHT-H and GRAD-L cases show opposite trends for increasing values of N_H. In Region 1, the RHT-H ξ and Z_Jow curves start to steadily decrease just after N_H = ~3.5 × 10²¹ cm⁻². In the GRAD-L case, the behavior is almost specular and the curves tend to rise in the highest-column-density bins, with the transition occurring at a similar value for N_H. In Region 2 the behavior is almost the opposite: a slightly rising trend for RHT-H and a clearly decreasing trend for GRAD-L. In both cases, the highest-density bin behaves opposite with respect to the global trend. In the GRAD-H case, the trend is flat in both Region 1 and Region 2. In all our considered cases, the values of ξ and Z_Jow are predominantly negative (perpendicular orientation), with a small fraction compatible with zero (no preferred orientation). The values of the parameters are positive (parallel orientation) only for the GRAD-L case in Region 2, in the two lowest-density bins (Fig. 9, middle row, right panel).

The Z_Jow curves follow the shape of the ξ curves, although the absolute values of the parameters for a given column-density bin can be quite different in some cases. The error bars from the PRS are typically smaller than for the HRO statistics.

We repeated our analysis of the GRAD method adopting a higher threshold for the gradient (see Appendix A). While the selected regions are less extended, we retrieved trends of ξ and Z_Jow consistent with the case based on the lower threshold. We fully recovered the flat behavior of the GRAD-H case, while in the GRAD-L case the trends are more extreme. In Region 1, the curves have a very steep rise, reaching a positive value in the highest-density beam. In Region 2, a sharp decrease is followed by a steep rise.

As already mentioned, one of the main finding from the analysis of Planck data is that the alignment of the density structures in MCs tends to be parallel to the local magnetic field, or without a preferred orientation, for low column densities (N_H ≲ 10^21.7 cm⁻²) and orthogonal to the direction of the field for high column densities (Planck Collaboration XXXV 2016). In the ten Gould Belt nearby MCs analyzed in Planck Collaboration XXXV (2016), the values of the histogram shape parameter ξ tend to decrease from positive to negative for increasing column density, indicating that the preferred relative orientation, calculated using the GRAD technique, goes from parallel to perpendicular when moving toward denser regions. The steepness of the slope changes in the different regions. This trend has been also recovered using HROs calculated from simulations (Chen et al. 2016). The analysis of BLAST-Pol data of the Vela C MC (Soler et al. 2017) has led to the same conclusions, and the observed trend has been confirmed by Jow et al. (2018) where both the Planck Gould Belt and BLAST-Pol Vela C data have been re-analyzed applying the PRS.

Malinen et al. (2016) analyzed Herschel data of the high-latitude MC L1642 using the RHT technique. They also found that the preferred relative orientation transitions from parallel to perpendicular to the magnetic field for increasing column density, and the transitions occurs at a relatively low threshold of ~N_H = 1.6 × 10²¹ cm⁻², which is equivalent to the value reported in Planck Collaboration XXXV (2016) if the same dust opacity value is assumed in both cases.

Alina et al. (2019) performed a statistical analysis of the relative orientation in the filaments hosting the Planck Galactic Cold Clumps (Planck Collaboration XXVIII 2016), using an improved version of the RHT method called supRHT. Their analysis revealed a change in the relative orientation between the magnetic field and the filaments, which depends on both the column density of the environment and the density contrast of the filaments. The combination of low-column-density environment and low-column-density contrast results in the preferential alignment of the filaments with the background magnetic field, while low-column-density contrast filaments embedded in high-column-density environments tend to be orthogonal to the magnetic field. When the clumps inside filaments are considered, the observed alignments are both parallel and perpendicular.

When we compare our findings with the results discussed above, it appears that they are consistent in some cases (Sect. 5.2). In Region 1, this happens in the RHT-H case while in Region 2 it occurs in the GRAD-L case. When the other method is applied (or at a different resolution) the results are different. In our MHD simulation low-density regions are not spatially refined, so their spatial resolution corresponds to that of the 128³ root grid, that is ~2 pc. Because this resolution is insufficient to capture the filaments often observed in low-density regions surrounding MCs, it is to be expected that positive values of ξ and Z_Jow, indicating alignment of B and low-column-density filaments, are mostly missing from our results. Thus, although a comparison between the results from the two methods is still appropriate, the comparison of our analysis with the observations should be limited to high-density regions. In such regions, we do find primarily negative values of ξ and Z_Jow, in agreement with the results reported in the literature based on Planck data.

Our findings reveal the complex relationships between magnetic fields and the density structures produced by turbulence and gravity, which change depending on, for example, the view direction and the general density of the region. The results also depend on the spatial scales that the analysis is targeting. However, the outcomes of the RHT and GRAD analysis were also seen to differ significantly, even when applied to the same area and focusing on similar scales.

The differences are not surprising, given that the methods try to recognize different topological features. RHT is searching for ridges, where the gradient is by definition close to zero. Conversely, a linear slope gives a strong signal in the GRAD method but none in the RHT analysis. The techniques are thus complementary. Only in the case of clear filaments at the chosen scale, the methods give consistent results, although still based on different parts of the structure.

The methods differ also regarding their sensitivity to different spatial scales. In the GRAD method, the size of the smoothing kernel of the Gaussian-derivatives technique (Soler et al. 2013) sets a lower limit for the size of the detected structures. There is no strict upper limit, although the method is still more sensitive to structures close to the smoothing scale, where the absolute values of the gradients are likely to be larger. In comparison, RHT is more clearly tied to a specific range of scales. That is first limited by the data resolution and the parameter D_K of the high-pass filter. The selected diameter D_W finally sets the preference for linear structures of at least comparable size. It is important that one can choose the scale precisely. This allows for a separate analysis of filament-scale and striation-scale objects or the generalization to multi-resolution analysis.

We conclude that RHT and GRAD can serve a similar purpose but are partly complementary in the investigation of structure orientations. The fact that the methods concentrate on different parts of the structure can mean that they perform differently depending on the resolution of the observations, for example, whether the gradients around the filaments are spatially resolved or not. The results of the two methods could also be combined into a single map of structure orientation. Although the interpretation of such a map may be more challenging, this can be useful when the methods provide significant measurements for clearly disjoint sets of pixels.

Acknowledgements

We are grateful to the referee for the careful reading and the insightful comments, and to S. Clark for useful discussions. We would like to thank D. Kodi Ramanah and G. Lavaux for providing plotting routines that were used as a basis for computing the “drapery pattern” in Fig. 1. E.R.M. and M.J. wish to acknowledgethe support from the Academy of Finland grant No. 285 769, P.P. acknowledges support by the Spanish MINECO under project AYA2017-88754-P, D.A. acknowledges fundings from the Ministry of Education and Science of the Republic of Kazakhstan state-targeted programme BR05236454 and from Nazarbayev University grant FCDRG110119FD4503, and J.M. acknowledges the support of ERC-2015-STG No. 679852 RADFEEDBAC.

Appendix A Results from the GRAD method applied with a high threshold

Fig. A.1 Same as Fig. 1 but showing the new region selected as reference, which provides a higher threshold for the gradient.

Fig. A.2 Analog of Fig. 2 (Region 1) but using a higher threshold for the gradient to make the pixel selection.

Fig. A.3 Analog of Fig. 3 (Region 2).

Fig. A.4 Analog of Fig. 8 (Region 1).

Fig. A.5 Analog of Fig. 9 (Region 2).

Fig. A.6 Analog of Fig. 10 (Region 1).

Fig. A.7 Analog of Fig. 11 (Region 2).

References

All Figures

Fig. 1 Surface brightness maps at 353 GHz (850 μm – shades of gray) resulting from radiative transfer calculations performed on density structure maps from MHD simulations. Each frame shows the maps corresponding to one observing direction (x-left, y-center, z-right). The number “315” labeling each panel refers to a snapshot taken 11 Myr after self-gravity has been switched on in the simulation. The maps have been convolved using a Gaussian beam with FWHM = 2.5 pix = 1.22 pc. The “fingerprint” pattern shows the direction of the magnetic field derived from the Stokes parameter, and the intensity of the red color increases with the polarization fraction. The red boxes indicate the regions where individual analysis was performed while the blue dashed box identifies a low surface brightness area selected as a reference region (see Sect. 4.2 for details).

In the text

Fig. 2 Comparison between the selected structures in Region 1, obtained from the RHT and GRAD methods. The recovered structures, indicated in blue, are overlaid on the brightness map (shades of gray). The RHT-H case (top-left) illustrates the results from RHT applied to surface brightness maps convolved using a Gaussian beam with FWHM = 1.2 pc. The GRAD-H case (top-right) results from GRAD using a convolution kernel of the same size, with the constraint that the module of the gradient must be greater than the average over a reference (smooth) region. The same constraint has been imposed in the GRAD-L case (bottom-left) with the difference that GRAD is applied using a Gaussian convolution kernel with FWHM = 7.2 pc. The bottom-right panel shows the superposition of the results for the RHT-H case (green) and the GRAD-L case (red). The yellow color indicates the regions in common between the two cases.

In the text

	Fig. 3 Same as Fig. 2 but for Region 2.
In the text

Fig. 4 Histogram of relative orientation between the RHT angle expectation value (RHT-H case) and the direction angle derived from the GRAD method (GRAD-L case) for Region 1 (top) and Region 2 (bottom). The angles were calculated for the pixels that are common to both methods (yellow areas in Figs. 2 and 3). The curves show the HROs for three column-density bins: the highest and lowest ones (red and black, respectively) and the intermediate bin between these two (blue). A ± 1 σ error is represented by the width of the shaded areas.

In the text

Fig. 5 Relative orientation between the RHT-H angle expectation value and the GRAD-L direction angle (magenta-blue-green) overlaid on the surface brightness map (shades of gray) for Region 1 (left) and Region 2 (right). The relative orientation goes from fully parallel (green, 0°) to fully perpendicular (magenta, 90°). As in Fig. 4, the angles have been calculated for the pixels that are in common to the RHT-H and GRAD-L cases.

In the text

Fig. 6 Analysis of Region 1 from Fig. 1, where the RHT and GRAD methods have both been applied to the structures selected by RHT (RHT-H case, top-left panel in Fig. 2). The columns show the surface brightness map overlaid with the red fingerprint pattern (left); the histograms of relative orientation (HROs – middle) between the magnetic field direction and the density structures; the ξ and ZJow parameters plotted as a function of column density, NH (right). For both parameters, a positive value indicates a preferred parallel orientation, a negative value a preferred perpendicular orientation and a value compatible with zero indicates that there is no preferred orientation. The areas between the solid and dashed contours on the maps (or inside them, if only one contour is present) highlight the regions corresponding to the column density bins shown in the HRO plots. Red and black identify the highest and lowest bins respectively, while blue indicates the intermediate bin between the two above. In the HRO plots, a ± 1σ error is represented by the width of the shaded areas. The top row shows the results for RHT-H, the middle row for GRAD-L and the bottom row for GRAD-H.

In the text

	Fig. 7 Same as Fig. 6 but for Region 2.
In the text

	Fig. 8 Same as Fig. 6 but in this case the RHT and GRAD techniques have been applied each to their own native selections, as shown by the three different panels in the left column.
In the text

	Fig. 9 Same as Fig. 8 but for Region 2.
In the text

Fig. 10 Relative orientation between the magnetic field and the density structures (blue-red) overlaid on the surface brightness map (shades of gray) in Region 1. The relative orientation goes from fully parallel (blue, 0°) to fully perpendicular (red, 90°). The top row shows the results from the analysis performed on the same pixel subsample for the RHT-H case (left), the GRAD-L case (middle) and the GRAD-H case (right). In the bottom row, each method is performed on its own native selection. Because the subsample chosen for the same-pixel analysis is the one from RHT (RHT-H case), the two panels on the left are the same.

In the text

	Fig. 11 Same as Fig. 10 but for Region 2.
In the text

	Fig. A.1 Same as Fig. 1 but showing the new region selected as reference, which provides a higher threshold for the gradient.
In the text

	Fig. A.2 Analog of Fig. 2 (Region 1) but using a higher threshold for the gradient to make the pixel selection.
In the text

	Fig. A.3 Analog of Fig. 3 (Region 2).
In the text

	Fig. A.4 Analog of Fig. 8 (Region 1).
In the text

	Fig. A.5 Analog of Fig. 9 (Region 2).
In the text

	Fig. A.6 Analog of Fig. 10 (Region 1).
In the text

	Fig. A.7 Analog of Fig. 11 (Region 2).
In the text