Statistical analysis of the interplay between magnetic fields and filaments hosting Planck Galactic Cold Clumps

We present a statistical study of the relative orientation in the plane of the sky between interstellar magnetic fields and filaments hosting cold clumps. For the first time, we consider both the density of the environment and the density contrast between the filaments and their environment. Moreover, we geometrically distinguish between the clumps and the remaining portions of the filaments. We infer the magnetic field orientations in the filaments and in their environment from the Stokes parameters, assuming optically thin conditions. Thus, we analyze the relative orientations between filaments, embedded clumps, and internal and background magnetic fields, depending on the filament environment and evolutionary stages. We recover the previously observed trend for filaments in low column density environments to be aligned parallel to the background magnetic field; however, we find that this trend is significant only for low contrast filaments, whereas high contrast filaments tend to be randomly orientated with respect to the background magnetic field. Filaments in high column density environments do not globally show any preferential orientation, although low contrast filaments alone tend to have perpendicular relative orientation with respect to the background magnetic field. For a subsample of nearby filaments, for which volume densities can be derived, we find a clear transition in the relative orientation with increasing density, at $n_{\rm H} \sim 10^{3}~{\rm cm}^{-3}$, changing from mostly parallel to mostly perpendicular in the off-clump portions of filaments and from even to bimodal in the clumps. Our results confirm a strong interplay between interstellar magnetic fields and filaments during their formation and evolution.


INTRODUCTION
First observed in extinction by interstellar dust particles (Schneider & Elmergreen 1979) and then in dust and CO emission (Abergel et al. 1994;Falgarone et al. 2001), filamentary structures in molecular clouds recently became a subject of many studies. There is a growing evidence that filaments play a fundamental role in the onset of star formation for low and high-mass stars. The unprecedented angular resolution and sensitivity of the Herschel maps of dust emission in the FIR wavelength range enabled to discover a network of filamentary structures ubiquitous in a wide range of environments of the ISM (Molinari et al. 2010;Motte et al. 2010;Men'shchikov et al. 2010;Miville-Deschênes et al. 2010;Juvela et al. 2012). Moreover, it was shown that E-mail: dana.alina@nu.edu.kz prestellar cores and protostars form primarily in the densest and gravitationally-bound filaments (André et al. 2010;Polychroni et al. 2013;Könyves et al. 2015). Investigating the origin and evolution of filaments is then crucial to better understand the early stages of star formation.
MHD simulations show that a hierarchy of sheets and filaments can form as a result of shock-compression and shear flow in supersonic turbulence, at the same time as the collapse or fragmentation of gravitationally unstable structures (André et al. 2014;Li et al. 2014). Magnetic fields are believed to play a key role in the formation of structures, but their interplay with turbulence and self-gravity still needs to be understood better at different spatial scales and evolutionary stages.
Observations of the starlight polarization revealed a connection between the magnetic field geometry and the dense filamentary structures, finding predominantly perpen-2 D. Alina et al. dicular relative orientations between the mean field and the long axis of the filaments (Goodman et al. 1990;Pereyra & Magalhães 2004;Franco et al. 2010;Chapman et al. 2011;Alves et al. 2008). More recently, parallel relative orientation was also evidenced between the mean field and low-density sub-filaments (Sugitani et al. 2011) or the striation features revealed with Herschel in Taurus (Palmeirim et al. 2013), Musca (Cox et al. 2016), and the high-galactic latitude cloud L1642 (Malinen et al. 2016). A bimodal distribution of relative orientations was both observed in molecular clouds of the Gould Belt (Li et al. 2013) and predicted from MHD simulations (Soler et al. 2013). However, the studies of the magnetic field orientations derived from observations of polarization in extinction are limited to small regions and depend on the background star distribution. Alternatively, measurements of dust polarized thermal emission can also be used to trace the magnetic field geometry as projected onto the plane-of-the-sky (POS). Such observations performed from ground-based telescopes at high angular resolution are well suited to the scale of prestellar cores, but they are limited to bright and small regions of molecular clouds (Dotson et al. 2000;Matthews et al. 2009;Ward-Thompson & Kirk 2000;Tang et al. 2009;Tassis et al. 2009).
The Planck 1 satellite survey has recently provided the first all-sky map of dust polarization emission (Planck Collaboration Int. XIX 2015; Planck Collaboration 2016). With unprecedented sensitivity measurements at 353 GHz and a resolution of 5 , these data are uniquely suitable to probe the dust polarized emission over a large range of column densities, at intermediate scales between molecular clouds and cold cores. Studying the relative orientation between the magnetic field and column density structures at 15 resolution and over most of the sky, Planck Collaboration Int. XXXII (2016) found that most of the elongated structures in the diffuse ISM are preferentially aligned with the magnetic field, while perpendicular structures appear in molecular clouds. In a detailed study of ten nearby molecular clouds of the Gould Belt, Planck Collaboration Int. XXXV (2016) showed that, as gas column density increases, the relative orientation changes from mostly parallel or having no preferred orientation to mostly perpendicular at the high column densities. A similar trend was recently evidenced in the Vela C molecular cloud observed with BLAST-POL at higher angular resolution (3 ), with a sharp transition in N(H) characterizing the change of relative orientation (Soler & Hennebelle 2017). Simulations show that such a change could be related to the degree of magnetization of the cloud, particularly when magnetic energy is in equipartition with turbulent energy (Hennebelle 2013;Soler et al. 2013;Chen et al. 2016). In a detailed study of three nearby filaments, Planck Collaboration Int. XXXIII (2016) showed that the mean magnetic fields in the filaments are coherent along their length and do not have the same orientation as in the background. Moreover, the relative orientation between the filaments' and background magnetic fields has been shown to be different in the three cases studied. All these studies suggest that the magnetic field must play an important role in the assembly of interstellar matter and the formation of dense structures. Further statistical observations are required to better understand and characterize the range of scales and densities for which the magnetic field has a significant influence. To better understand early stages of star formation, it is particularly important to study if and how the properties of the ambient magnetic field are related to the filament evolution and prestellar core formation.
The Planck Catalog of Galactic Cold Clumps (PGCC) is built from the Planck full all-sky survey (Planck Collaboration Int. XXVIII 2016). The catalog contains more than 13,000 cold sources distributed across the whole sky, from the Galactic plane to high latitudes, and are mainly situated in the molecular clouds. The sources span a broad range in temperature, mean density, mass, and size. Because of the limited angular resolution, the most compact and nearby sources have a linear diameter of ∼ 0.1 pc, though at large distances, many sources have an intrinsic size of tens of parsec. As revealed by higher angular observations performed with Herschel follow-up on a representative sub-sample of 350 PGCC targets, these clumps are characterized by the presence of substructures. The clumps often contain colder and denser individual starless or prestellar cores and young protostellar objects still embedded in their colder surrounding cloud (Juvela et al. , 2011Juvela & Ysard 2012;Montillaud et al. 2015;Planck Collaboration XXII 2011). As shown with a statistical analysis based on molecular dense tracers, most of the Planck clumps correspond to an early stage of star formation (Wu et al. 2012;Liu et al. 2015;Tatematsu et al. 2017). Furthermore, it was found that a substantial fraction of them are embedded in filamentary structures . A statistical analysis of their polarization properties could bring new insights on the formation of filaments and clumps potentially related to early stages of star formation.
In the present work, our goal is to investigate the possible interplay between the magnetic field and the filaments hosting PGCCs. Having selected a sample of these filaments surrounded by an ordered magnetic field, we perform a statistical analysis of the relative orientations between the filamentary structures and that of the magnetic field projected onto the POS. This paper is organized as follows. The Planck data and the PGCC catalog are presented in Section 2. The description of the methods used for the filament detection and the separation of the polarization properties (I, Q, and U Stokes parameters) intrinsic to the filaments from those of the surrounding background are presented in Section 2.2. Section 3.1 presents the properties of the filaments selected for this statistical study. We analyze in Section 4 the relative orientation between the filamentary structures and the magnetic fields either in the background or within the filaments. Finally, we discuss the implications of the results and present our conclusions and perspectives in Section 5.

The Planck data used and the PGCC catalog
In our analysis, we use the Planck 353 GHz full mission data, which is a part of the Planck Product 2015 release. In order to improve the signal-to-noise ratio, the data are smoothed from the nominal resolution of the maps (5 ) to a resolution of 7 . The smoothing is performed using a convolution with a Gaussian kernel of 5 .
The Planck all-sky Galactic Cold Clump catalog (Planck Collaboration XXVIII 2015) contains 13,188 sources detected using a dedicated method (Montier et al. 2010;Planck Collaboration XXIII 2011). The method consists of a color-based detection of the excess emission in the cold residual maps after subtraction of the warm component traced by IRAS 100 µm observations. Clumps were fitted using elliptical Gaussian functions, and the catalog contains the resulting semi-major, semi-minor axes, and the position angles of the ellipses (noted as θ cat hereafter). With the exception of some sources in the Large Magellanic Cloud and in the Small Magellanic Cloud, the catalog excludes extragalactic sources (Planck Collaboration XXVIII 2015) and contains nearby cores, clumps, and distant molecular clouds. This makes the catalog particularly useful in the studies related to the early stages of star formation. The Planck data are in the HEALPix format, with a resolution parameter N side of 2048. We extracted small 2 • square maps centered on the positions of the clumps from the catalog. Such maps will be called minimaps hereafter.
We select cold clumps from the PGCC catalog using the following criteria: • clumps are within the galactic latitude range 2 • ≤ |b| ≤ 60 • to avoid confusion along the line-of-sight (LOS) in the Galactic plane; • clumps are detected with accurate flux density estimates in the Planck 545 and 857 GHz channels as well as in the IRIS 3 THz band. This corresponds to the PGCC catalog flag FLUX_QUALITY=1. We also include fainter sources: clumps for which there is no detection in the IRIS 3 THz band (flag FLUX_QUALITY=2).
• clumps are detected with the signal-to-noise ratio SNR ≥ 8 in the three highest Planck frequency bands (857, 545 and 353 GHz).
• clumps have an average signal-to-noise ratio (SNR) on the polarized intensity (P = Q 2 + U 2 ) larger than 2. We use SNR(P) as it reflects the joint SNR of the Stokes parameters Q and U used in the estimation of the polarization angle (see Equation 2). SNR(P) is estimated as P/σ P , where σ P = √ σ Q σ U . Such an estimate of the uncertainty on the polarized intensity has been discussed by Montier et al. (2015). SNR(P) is calculated in a disk of 10 radius around the position of the clumps and in the annulus with the inner and outer radii of 30 and 50 . We require P/σ P ≥ 2 in both regions.
This yields a selection of 3,743 clumps.

Filament detection method and detection parameters
Several methods exist to detect filamentary structures (Soler et al. 2013;Sousbie 2011;Planck Collaboration XI 2014). Among those, we opted for the SupRHT method, which is an improved version of the Rolling Hough Transform proposed by Clark et al. (2014). It is a machine vision method that allows one to determine whether each pixel in an image belongs to a curve of a given shape. The advantage of this method is that it is independent of the intensity of the structure with respect to the background intensity, as a top-hat filtering is performed before calculations. Another advantage is that it provides an estimate of the uncertainty on the angle. We use the 353 GHz Planck minimaps, which we convolved with a bar kernel of length k s and width k w as an input to the method. As a result, the SupRHT procedure gives the intensity I RHT , the angle θ, and the uncertainty on the angle σ θ . θ is counted positively in the direction from North to East and ranges between −90 • and 90 • (cf. Fig. 1). In the following, all the angles used in the analysis are brought to the same convention.
We test kernels of different sizes: k s = 11 , 21 , 31 noted as K11, K21, and K31, respectively, and observe their effect on the detection of the filaments in the intensity minimaps. These selected sizes are motivated by the sizes of the elongated clumps and by the resolution of the data. K11 yields positive detection for Bok globules. K11 and K21 yield positive detections for cometary shape clumps. K31 is the most restrictive in detecting linear filaments that span in two directions from the center of the minimaps. The kernel width k w = 6 is found to eliminate the round-shape isolated clumps. Larger widths do not allow us to detect the filaments at the position of the clumps because of the average size of the clumps, which is about 6 . Resulting maps of I RHT better reflect the variations of the structures compared to larger values of k w whereas smaller values are not appropriate as they do not allow us to take the clump width into account.
We run the intensity minimaps centered on the preselected PGCCs through the SupRHT procedure. The raw SupRHT maps contain several overlapping filamentary structures. We isolate only the filamentary structure associated with the PGCC in the center of each minimap. Hereafter, we call such filamentary structures as filaments. We build a mask to isolate the clump ellipses defined by the semi-major and semi-minor axes provided in the PGCC catalog and the position angles θ. We also build a mask that gives the filaments without the clump areas. The number of pixels in a filament ranges from ∼ 180 to ∼ 1200 with an average of ∼ 500. The number of pixels in a clump ranges from ∼ 50 to ∼ 500 with an average of ∼ 200, as shown in histograms in Figure 2.
For the clump in the center of each minimap, we fit the histograms of the angles θ inside the clump ellipse with a Gaussian function. Hereafter, we normalize the histograms to one and call them as the distribution functions (DFs). The fitting by a Gaussian allows us to derive the position of the peak, θ clump , and the standard deviation of the DFs, σ θ, clump . We use the Gaussian fit to consider the average angle of the majority of the pixels in a filament and not be affected by eventual separately distributed angles. We define δθ as the difference between θ clump from the Gaussian fit and the value of the SupRHT angle at the center of the clump θ 0,0 .
The parameters (σ θ , δθ) allow us to ensure that the angle θ does not change much over the area of the clump's ellipse so that the pixels belong to the same detected filamentary structure. We also calculate the average angles in the regions in the two directions given by the central angle  θ 0,0 delimited by two circles at the distances equal to 1.5 and 2 times the clump's major axis and the filament area (if detected), as shown in Figure 1. The absolute differences between θ 0,0 and the averages in the two regions are called ∆θ up and ∆θ down . These parameters allow us to control the apparent length of the filaments and to put constraints on their curvature. Moreover, as we deal only with the POS projection of the structures (and of the magnetic field), detectable filaments must lie in the POS, or within a small angle.
We apply the following criteria to the above parameters for the detection of the linear filaments in our study: • I RHT 0 in the center of each minimap; • the uncertainty of the θ at the center of each minimap satisfies σ θ < 1 • , This yields a selection of 160 fields, each centered on distinct PGCCs. Different PGCCs can belong to the same filamentary molecular cloud. To avoid statistical bias from considering the same filament several times, we withdraw the minimaps showing the same filament but centered on different PGCCs. In the rest of the analysis, we consider the contribution of all the PGCCs belonging to a filament Thus, we pre-selected 137 fields. Note that 38 filaments out of 137 fall into the regions studied previously using the Planck data by Planck Collaboration Int. XXXV (2016); Planck Collaboration Int. XXXIII (2016). However, most of the selected fields have not been subject to the previous statistical studies on the relative orientation of the magnetic field and matter structures in the Planck data.

Coherence of the background polarization angle
The polarization signal towards a filament comes from the filament itself and both the background and foreground interstellar medium. For convenience, we call the background and foreground contributions as the "background" in the following. We separate the background contribution from the observed signal using the following simple approach. We assume that the emission observed in the immediate surrounding of the filament is the same as the background emission at the location of the filament. Thus, the total signal towards the filament is where X = I, Q, U and the "fil" and "bkg" subscripts stand for "filament" and "background" respectively. Such an approach is only valid for cases where the background is the same on both sides of a filament and for optically thin media which is the case at the Planck 353 GHz channel. As our study is focused on the magnetic field and matter relative orientations, we seek filaments having uniform polarization angle in the background. We define three rectangular regions for each minimap: the "central" region containing the filament and two regions on both sides ("left" and "right") with their longest dimensions being parallel to each other. The position angles are given by θ at the position of the clump. The lengths of the rectangles are the same and they are defined by the extremities of the I RHT filament. The width of the central rectangle is equal to twice the width of the I RHT filament at the position of the clump. The width of the side rectangles is fixed to 30 , which corresponds to the average width of the central region. The two side regions together constitute the background region of each minimap. Figure 3 shows examples of how we divide a minimap into the three regions.
In the side rectangles, we calculate the average Stokes Q and U parameters (Q,Ū) from which we derive the corresponding average magnetic field orientation angleφ as pro-Statistical analysis of the interplay between magnetic fields and filaments 5 jected onto the POS: which we call magnetic field angle hereafter. Here ψ is the polarization angle. As the angles in the Planck data are counted positively in the direction from North to West, we bring it to the IAU convention adopted for all the angles in this study by taking the negative value of the Stokes U parameter. Left panels of Figure 4 show the averages and the dispersions of Q and U for each filament in each of the side background regions. We observe that, globally, the Stokes linear polarization parameters have similar values on both sides of the detected filaments. However, for some minimaps, the difference between them as well as the dispersion of the Stokes parameters on each side can be high. Dispersion of the Stokes parameters individually does not represent the angles dispersion. We calculate the standard deviation of φ in each region using circular statistics (Berens 2009), which we call σ φ . It allows us to characterize the coherence of the magnetic field in the background and to refine the field selection for our analysis.

Separation of the filament and the background polarization parameters
In this study, we focus on the relative orientation between the filaments and the magnetic field angles given by the sub-millimeter polarization angle measurements. In order to investigate the contribution from the filament and its background in the polarization angle maps, we separate their respective contributions. For each pixel in a filament, we estimate the average background Stokes polarization parameters. If the uncertainty for a given pixel is smaller than 10 • , we use cuts passing through the side rectangles perpendicular to the pixel position. Otherwise, the cut is done perpendicular to the angle of the central pixel, that is at the position of the corresponding PGCC. See Figure 3 for examples of cuts. Then, the average Q and U are calculated over the cuts and are assigned to the pixel i in the filament (Q bkg,i , U bkg,i ). As a first approximation, the polarization parameters of the filament at a pixel i are obtained by subtracting the background contribution from the total signal: The background magnetic field angle φ bkg and the magnetic field angle of the filament φ fil are then calculated using Equation (3) for each pixel in the filament. Figure B1 shows the average Stokes parameters as well as the standard deviations for each filament before and after the estimated background level has been subtracted. On average, the background subtraction does not critically change the linear polarization parameters and the overall values are within the amplitude of the initial Q and U parameters in the minimap. The above method is applied to the filaments selection described further in Section 3.1. These filaments are embedded in an environment where the magnetic field orientation is uniform so that the assumption of linear contribution is valid.

Column density
We determine the column density for each pixel using where we adopt the dust opacity κ ν =0.1 (ν/1 THz) β cm 2 g −1 (Beckwith et al. 1990), which is appropriate for the case of dense clouds at intermediate densities. µ is the mean molecular weight, which equals to 2.8 amu in our analysis. The color temperature T is estimated using the brightness emission in the Planck 857, 545, 353 GHz channels, and the IRIS 3 THz (100 µm) band. The maps were convolved to a 7 resolution and, for each pixel the SED was fitted with B ν (T)ν β using β = 2. This value is usually adopted for starless core studies and it is close to the mean value of β of the derived on Planck clumps (Planck Collaboration XXIII (2011), Planck Collaboration Int. XXVIII (2016)). I ν is derived from the flux density at 857 GHz, the closest Planck band to the 1 THz reference of Beckwith et al. (1990), which minimizes the impact of assuming a fixed emissivity spectral index. This frequency also corresponds to the Planck HFI band where the SNR is the strongest.

The "uniform background" selection
About 26% of the pre-selected 137 filaments have a magnetic field angle dispersion larger than 20 • in both background regions. Thus, most filaments are located in regions where the magnetic field seen by Planck is coherent. To define filaments that have a uniform magnetic field in the background, we require thatφ does not change more than 20 • between the two sides of the filament. This yields a selection of 92 objects distributed over the sky (see Figure 5). The rest of our analysis is focused on this selection of filaments. By comparing the left and the right panels of Figure 4, we note that the final selection excludes the fields with large dispersion in the Stokes parameters.

Characteristics of the selection
The temperatures are reported in the PGCC catalog for 80% of the selection and span between 9 K and 19 K with a mean value of 13.7 K. The distance estimates are provided for about 45% of the filaments and are mostly lower than 1 kpc. The corresponding histogram is presented in Figure B2 and shows a main peak around 0.2 kpc , and second at ∼ 0.4 kpc. The masses are determined for 38 clumps. The masses of the 5 most distant clumps exceed 100 solar masses and approach the masses of molecular clouds. However, 21 clumps, which represent 23% of the selection, have their mass determined to be smaller than 5M with a mean value of 2.02M . The gas column density values provided in the PGCC catalog are integrated over the detection ellipse. In our selection, their values span the range from 1.4 × 10 20 cm −2 to 3.1 × 10 22 cm −2 with a mean value of 3.2 × 10 21 cm −2 . The average gas column densities calculated according to Equation (5) range from 2.9 × 10 20 cm −2 to 2.5 × 10 22 cm −2 with a mean value of 3 × 10 21 cm −2 .   We estimate the average size of the clumps as follows: where a and b are the major and minor semi-axes. The average sizes range from 10 to 17 with an average of ∼ 13 . The corresponding histogram is presented in Figure B3. The lengths of the selected filaments that are computed as the sum of the two line segments connecting the center of the clump to the extremities of the filament, range from 34 to 104 with an average value of 66 . The widths of the filaments taken in the center of the clumps in directions perpendicular to the detected θ, spread from 3 to 20 with an average value of 8 . We emphasize that the length and width given here are not the physical dimensions, but are detection width and length. These parameters depend on the detection method and the kernel size used in the SupRHT method. A conservative estimate of the dispersion of the magnetic field angle in the filaments is calculated using circular statistics method. It ranges from 6 • to 73 • with a mean value of 40 • . The dispersions of φ for each filament including and excluding the clump areas are presented in Figure B4.

The low and high background column density subsamples
Following the observation that the preferential orientation of the matter in filaments and the background magnetic field changes with gas column density (Planck Collaboration Int. XXXV 2015), we divide our selection of filaments into two subsamples with high and low density environments. We separate these two using a threshold density value equal to the median background column density of 1.2 × 10 21 cm −2 . The corresponding subsamples will be called N high H,bkg and N low H,bkg . We will distinguish between the filaments that are more or less dense with respect to their environment also by calculating a rough estimate of the differential column density for each filament: whereN H is the total average column density over the filament area andN H, bkg is the total average column density in the background. Two of the filaments have negative column density excess. A detailed inspection of the corresponding intensity minimaps shows that these filaments are very faint and averaging over the background that have some prominent features yields negative ∆N H . These filaments are not taken into account in further analysis. We will denote ∆N low H and ∆N high H the filaments that have differential column densities lower and higher than 4×10 20 cm −2 , which is the median value over the ∆N H of the selection.
The histograms of the column densities normalized to one for various target selections are shown in Figure 6. The number of the detected filaments in each subsample is reported in Table 1. The distribution of number of pixels in clumps and filaments in the different subsamples are roughly the same and follow the general distribution shown in Figure  2. We emphasize that the PGCCs chosen for this study are outside the Galactic plane and the observed column densities are more probably due to the local structures rather than the averaging along large lines of sight. We observe in Figure  Figure 6. Distribution function of the column density in the pixels of the detected filaments. The black solid line represents all filaments, the green and purple lines show the low and high background column densities subsamples. The gray dotted line represents the background regions, while the turquoise line corresponds to the high background subsample with high differential column density. 6 that with our selection we probe environments with different densities, shown by the gray dotted line. We also note that the filaments embedded in high density environments have higher column densities. The same behavior was observed in Herschel data at smaller scales (Rivera-Ingraham et al. 2017), where dense filaments are embedded in dense environments.

Relative orientation between filaments and embedded clumps
In this section, we calculate the absolute difference | θ fil − θ clump | between the average orientation angles of the clumps and of the filaments (cf. Section 2.2). For this purpose, only elongated clumps with ellipticities larger than 1.5 have been used. As the orientations of both clumps and filaments do not have a preferential direction, we use the smallest angle between the two line segments, which is in the range from 0 to 90 • . Figure 7 shows the corresponding distribution function over the final selection. Most (∼ 98%) of the clumps are aligned with the filament within 20 • . The same is observed when the clumps of ellipticities smaller than 1.5 are included in the analysis. In order to check that this result is independent of the detection method, we also  compare the average clump position angles θ clump given by SupRHT to the position angles θ cat provided in the PGCC catalog, which is based on the cold residual maps. The corresponding histogram is shown in Figure 8. Both methods agree on the position angle of the clumps within 20 • for more than 70% of the selection. If we consider the clumps of all ellipticities, this fraction goes down to 62% but remains significant.

Relative orientation between the filaments and the background magnetic field
In this section, we investigate the statistics of the relative orientation between the background magnetic field and the filaments. For this purpose, we calculate the absolute difference between the average background magnetic field angle and the position angle |φ bkg − θ| for each pixel. This quantity is also defined in the range from 0 • to 90 • . We calculate the DFs over the pixels in all the filaments, in the N high H,bkg and N low H,bkg subsamples and over the pixels belonging to the clumps as well as over the pixels belonging to filaments without clump areas in Figures 9, 10 and 11. We do not calculate this quantity per filament because the size of our sample is not sufficient for significant statistics.
We see in the upper panel of Figure 9 that the DFs of both subsamples exhibit peaks at 0 • . However, the distribution for the filaments embedded in low density environment decreases strongly with increasing value of |φ bkg − θ| whereas the distribution for the filaments embedded in high density environment is almost flat. If we take into account only the filaments (middle panel of Figure 9), we observe that the filaments embedded in dense environment have no preferential alignment whereas those embedded into a low density environment have a clear tendency to be in parallel alignment with respect to the background magnetic field. When considering only the clump areas (lower panel), both DFs peak at 0 • , decrease towards intermediate angle values and increase again at large angles. The uncertainties of the DFs are estimated using the bootstrap technique consisting of random sampling with replacement (Efron & Tibshirani 1993). The details on the uncertainties calculation using independent pixels are described in Appendix A and the histograms of the relative uncertainties are given in the Appendix B. The errors shown in Figure B5 on the DFs represented in Figure  9 are around 10% for the DF on all pixels and 15% for the DFs only for filaments and clumps and for the DFs of the subsamples. For each subsample, we separate the contributions from the filamentary structures that are either denser or less dense than their environment. The resulting DFs are shown in Figures 10 and 11. In the N low H,bkg subsample, the filaments that are much denser than their environment tend to be much less aligned than those that are not prominent with respect to their background, as shown in the upper panel of Figure 10. Such a trend is more pronounced when considering only the clumps. The denser clumps are located in denser filaments with the lowest average density per clump of 1.6 × 10 21 cm −2 in the ∆N high H subsample of the N low H,bkg filaments. We continuously separate the clumps into two parts depending on their average column density and recover the same kind of figure as in the bottom panel of Figure 10 for the values from 1.4×10 21 cm −2 to 1.6×10 21 cm −2 . There is a clear dichotomy depending on the density: denser clumps show a perpendicular alignment whereas the lower density clumps show parallel alignment with respect to the background magnetic field. The corresponding relative uncertainties are presented in Figure B6. For illustration purpose, we report the error bars in the lower panel of Figure 10 in the bin where most difference between the histograms is observed. The rest of the uncertainties are of the same level of magnitude. In the N high H,bkg subsample, the tendency for the perpendicular relative ori- filaments show preferential alignment with the background magnetic field orientation in the histogram shown in Figure 11 with the corresponding uncertainties for each bin in Figure B7. In this case, we do not observe a dichotomy in the relative orientation of the clumps depending on their density (the DF for clumps only is not shown here).

Relative orientation between the filaments and the magnetic fields in the filaments
In order to analyze whether there is any coherence between the filament orientation and the magnetic field orientation  inside the filaments, we calculate the absolute difference between the respective angles (θ and φ fil ). The resulting distribution function and its clumps and filaments counterparts of the differential column density subsamples are presented in Figure 12 and their respective uncertainties are shown in Figure B8. All the DFs show strong peaks at 0 • . However, the filaments, which are denser than their environment (∆N high H ) regardless of their background density, show a peak between 70 • and 90 • . When comparing the contributions of the clumps and filaments in the middle and bottom panels of Figure 12, we observe that, for the ∆N high H subsample, the orientation of the clumps is mostly perpendicular or parallel, whereas the filaments themselves tend to be aligned to the inner magnetic field orientation. However, such a behavior is not observed in the ∆N low H subsample which show preferential parallel alignment with respect to the magnetic field in both filaments and clumps counterparts.
For the clumps counterpart in the ∆N high H subsample, we separate the contributions of the low and high density background subsamples. The corresponding histograms are shown in Figure 13 and the relative uncertainties for each bin in Figure B9. For illustration, we report the relative uncertainties for one bin with the largest difference between the histograms. The rest of the uncertainties are of the same level of magnitude. The perpendicular orientation between matter and magnetic field is clearly observed inside the filaments that are located in low density environment. One of the reasons why such a behavior is not observed in filaments located in dense environment could be the issue of the component separation method. The subtraction of high background level yields very low polarization in the filaments. For the filaments in the N high H,bkg subsample the subtracted background level and a proper determination of the inner magnetic field angles becomes impossible.
We also derive the relative orientation between the magnetic field in filaments and the orientation of clumps. The average magnetic field angle in the filaments ( φ fil ) is computed using the average Stokes Q and U parameters through Equation (3) after subtraction of the background contribution using Equation (4). The DF of the absolute difference between φ fil and θ clump (cf. Section 2.2) is presented in Figure 14. The statistics are low but a peak at low angles is observed suggesting that on average the clumps are aligned with the inner magnetic field in the filaments. There are less filaments showing the intermediate angle values between the orientation of clumps and the average magnetic field in the filaments. An increase in the histogram at larger angles, toward 70 • , is also observed. However, as the statistics are low and as there are much less clumps aligned perpendicular to the magnetic field, we can not conclude about the origin of this peak.

Relative orientation between the magnetic fields in the filaments and in the background
The difference between the magnetic field angles in the filament and in their background is calculated using the Stokes parameters of both regions (Planck Collaboration Int. XIX 2015) and is equal to the difference between the respective polarization angles: (8) Figure 15 shows the corresponding distribution functions over the whole selection and the N H, bkg subsamples. We observe that the fields are mostly parallel to each other. The relative orientation between the magnetic fields inside and outside the filaments appears to depend on the environment: the DF of the N high H,bkg subsample is flatter than the DF of the N low H,bkg subsample. For illustration purpose, repre-  sentative error bars are reported in dotted lines for the bin corresponding to the highest value. The rest of the uncertainties shown in Figure B10 and are of the same level of magnitude. Figure 16 shows the histograms of the absolute difference between the average magnetic field angles in the filaments and in the clumps calculated using average Stokes parameters. We observe that the shape of the histograms are similar to the shape of the DFs shown in Figure 15. Thus, Figure 15 describes mostly the relative angle variation within the two subsamples. However, we remind that the coherence of the magnetic field in the filaments varies case by case (cf. Figure B4), so that in Figure 15 there is also a counterpart which comes from the magnetic field variations within filaments. Also, the derived Stokes polarization parameters in the clump areas include both filament and clump contribution. We estimate that a proper determina-  tion of magnetic field angles in the clumps is not possible because of the low resolution of the data and low polarization signal.
Planck Collaboration Int. XXXIII (2016) estimated the average difference between the magnetic field orientation in the filaments and in their background for three filaments: Musca, Taurus L1506 and B211. They used polynomial fits to the Q and U parameters after subtraction of the background contribution. The first two filaments (Musca and L1506) are included in our final selection. For these two clouds, we obtain the average values of 8.6

SUMMARY AND DISCUSSION
We presented a statistical analysis of the relative orientation between the magnetic field angles and the filaments detectable in the 353 GHz Planck data and hosting the Planck Cold Clumps. We identified a selection of 90 such filaments all over the sky excluding the Galactic plane, that are embedded in an environment with a uniform POS magnetic field orientation. The selection of filaments used in this study span a range of background gas column density N H from 1.2 × 10 20 cm −2 to 3.8 × 10 21 cm −2 (averaged values at the angular resolution of 7 ). We separate the contribution to the observed polarization signal from the filaments and the background assuming optically thin emission and taking into account the two-dimensional geometry of each region. The advantage of such an approach is in its applicability to large samples. The relative angles derived using this method are consistent with the results from previous studies of a sample of Planck filaments, where a polynomial fit to the filament intensity profiles were used (Planck Collaboration Int. XXXIII 2016).
We confirm that the Planck Galactic Cold Clumps main axes are statistically oriented along the filaments main axes. We show that the relative orientation between the magnetic fields and the filaments strongly depends not only on the density, but also on the nature of the structures such as the clumps or the filaments excluding clumps. Inside the filaments as a whole, the perpendicular alignment that was previously observed by Planck Collaboration Int. XXXII (2016) in the high density clouds, seems to take place in the clumps whereas the filaments themselves are mostly parallel to the magnetic field. In the clump areas, both perpendicular and parallel alignment with respect to the magnetic field are observed.
The filaments embedded in dense environment that correspond to the densities above N H = 1.2 × 10 21 cm −2 in our study, do not show a preferential alignment with the background magnetic field whereas those situated in low density environment generally tend to be aligned. There is a clear alignment of the filaments in low column density regions, especially for the low contrast filaments seen against low density background. In the same subsample of the filaments embedded in low density background, we find that dense and low density clumps exhibit different behaviors: clumps with the average densities larger than 1.5×10 21 cm −2 in our study are mostly perpendicular to the background magnetic field whereas clumps with lower average densities are parallel. We also showed that the alignment between the magnetic fields in the filaments and in their background is weaker in the filaments embedded in dense environments. This could indicate a possible disconnection between the magnetic fields in dense regions.
In our analysis we did not account for projection effects. As pointed out by Planck Collaboration Int. XXXII (2016), the projection effects account for the broadening of the distribution functions at parallel alignment whereas the observed perpendicular trend indicates a true perpendicular relative orientation.
The results presented in our study indicate that particular processes between the matter and magnetic field occur in the clumps and are already detectable at the resolution of the Planck data. This study brings information on the role of the magnetic fields in the evolution of the dense ISM matter at the intermediate scale between molecular clouds and dense cores. A dedicated analysis on known molecular clouds could bring insights on the origin of the parallel or perpendicular alignment depending on the evolution of the clouds and/or efficiency of star formation in the clouds. We would like to emphasize that the contribution of the Cold Clumps should be considered in the future studies of the filamentary structures.
Both alignment regimes are observed in the Planck data. A more precise study of the interplay between matter and magnetic fields in the Cold Clumps at higher resolution could disentangle the dependence of the alignment on the geometry of the clump and of its evolutionary stage such as the Planck follow-ups of the PGCCs with the Scuba-2 polarimeter at JCMT (Tie et al. and Juvela et al. in preparation). Table 1: Characteristics of the filaments used in the study. Columns: number in the PGCC catalog, galactic longitude, galactic latitude, length, width, average position angle given by the SupRHT method using Gaussian fit, its uncertainty, average column density over the filament and its uncertainty, average column density over the filament excluding clumps areas and its uncertainty, average column density over the background region and its uncertainty, average magnetic field angle over the filament using Gaussian fit, average magnetic field angle over the filament using circular statistics method and its uncertainty.  Figure B1. Scatter plot between the average Stokes parameters for each filament before and after the background subtraction in the filament pixels (Q is shown in the right panel and U is shown in the left panel). The error bars show the standard deviation in each filament. The lines show the one-to-one correlation. Figure B2. Histogram of the distances provided in the PGCC catalog for a half of the filaments in our study.

APPENDIX A: ESTIMATION OF UNCERTAINTIES
In order to properly estimate the uncertainties, we proceed to bootstrap on independent pixels that corresponds to pixels situated at the minimum distances of 7' between each other, which is equal to the resolution of the data. This diminishes the number of the pixels used. Thus, there are only 3% of the pixels considered in each bootstrap round. The histograms of uncertainties are calculated as the ratio of standard deviations of the 5,000 realizations of bootstrap over the values in each bin of the corresponding histograms (σ N /N). We have checked that the average bootstrapped DFs are compatible with the DFs obtained directly from the data.

APPENDIX B: EXTRA FIGURES
This paper has been typeset from a T E X/L A T E X file prepared by the author. Figure B3. Histogram of the average clump sizes. Figure B4. Dispersion of the magnetic field angles in each filament excluding clump areas versus dispersion of the magnetic field angles in each filament including clump areas.
MNRAS 000, 1-18 (2017) Statistical analysis of the interplay between magnetic fields and filaments 19 Figure B5. Relative uncertainties to the DFs presented in Figure  9 calculated over 5,000 realizations of bootstrap. Figure B6. Relative uncertainties to the DFs presented in Figure  10 calculated over 5,000 realizations of bootstrap. Figure B7. Relative uncertainties to the DFs presented in Figure  11 calculated over 5,000 realizations of bootstrap.
MNRAS 000, 1-18 (2017) 20 D. Alina et al. Figure B8. Relative uncertainties to the DFs presented in Figure  12 calculated over 5,000 realizations of bootstrap. Figure B9. Relative uncertainties to the DFs presented in Figure  13 calculated over 5,000 realizations of bootstrap.