Unveiling polarized emission from interstellar dust of the Large Magellanic Cloud with Planck

D. Alina^1,2, J.-Ph. Bernard², K. H. Yuen³, A. Lazarian^3,4, A. Hughes², M. Iskakova⁵, A. Akimkhan⁶, A. Mukanova¹
¹Department of Physics, School of Sciences and Humanities, Nazarbayev University, Kabanbay batyr ave, 53, Nur-Sultan 010000, Kazakhstan
²IRAP, Université de Toulouse CNRS, UPS, CNES, F-31400 Toulouse, France
³Department of Physics, University of Wisconsin-Madison, Madison, WI 53706, USA
⁴Centro de Investigación en Astronomía, Universidad Bernardo O’Higgins, Santiago, General Gana 1760, 8370993, Chile
⁵Washington University in St Louis, One Brookings Drive, St. Louis, MO 63130
⁶Nazarbayev Intellectual School, Khussein Ben Talal, 21, Nur-Sultan 010000, Kazakhstan dana.alina@nu.edu.kz

(Accepted XXX. Received YYY; in original form ZZZ)

Abstract

Polarization of interstellar dust emission is a powerful probe of dust properties and magnetic field structure. Yet studies of external galaxies are hampered by foreground dust contribution. The study aims at separating the polarised signal from the Large Magellanic Cloud (LMC) from that of the Milky Way (MW) to construct a wide-field, spatially complete map of dust polarization using the Planck 353 GHz data. To estimate the foreground polarization direction, we used velocity gradients in HI spectral line data and assessed the performance of the output by comparing to starlight extinction polarization. We estimate the foreground intensity using the dust-to-gas correlation and the average intensity around the LMC and we assume the foreground polarization to be uniform and equal to the average of the MW around the galaxy to derive foreground I, Q, and U parameters. After foreground removal, the geometry of the plane-of-the-sky magnetic field tends to follow the structure of the atomic gas. This is notably the case along the molecular ridges extending south and south-east of the 30 Doradus star-forming complex and along the more diffuse southern arm extending towards the Small Magellanic Cloud. There is also an alignment between the magnetic field and the outer arm in the western part. The median polarization fraction in the LMC is slightly lower than that observed for the MW as well as the anti-correlation between the polarization angle dispersion function and the polarization fraction. Overall, polarization fraction distribution is similar to that observed in the MW.

keywords:

Galaxies: general, ISM, Magellanic Clouds, magnetic fields.

^†^†pubyear: 2022

1 Introduction

Along with turbulence and gravity, the magnetic field is thought to play a crucial role in star formation and in shaping the interstellar medium (ISM). It is therefore essential to consider magnetic fields in the studies of the ISM, and, first of all, to map the geometry of the magnetic fields in galaxies. The Large Magellanic Cloud (LMC) is an irregular galaxy (Wielebinski, 1993) with low metallicity (Olszewski et al., 1991; Matsuura et al., 2009; Grady et al., 2021). It is the closest gas-rich galaxy and is viewed almost face-on. Polarimetric observations of interstellar dust reveal information about the magnetic field and the dust physics. The direction of the observed polarization is related to the plane-of-the-sky (POS) projection of the magnetic field. The amount of polarization depends in principle on dust properties and the intensity of the ambient radiation (Lazarian & Hoang, 2007), but for diffuse interstellar media it is shaped primarily by the 3D structure of the magnetic field (Planck Collaboration Int. XIX, 2015).

The alignment necessary to produce a net polarisation is thought to occur through phenomena such as Radiative torques (RATs) which lead to align the shortest grain axis parallel to the ambient magnetic field. Consequently, the dichroic absorption of the incident starlight by these elongated grains yields the observed visible linear polarisation to be parallel to the POS magnetic field direction, while in emission, the polarization direction is perpendicular to POS magnetic field direction.

The LMC magnetic field was first studied using visible extinction measurements in polarization (Visvanathan, 1966). Spatially incomplete measurements showed that the LMC has an ordered large-scale magnetic field over some regions, and that the polarization is maximum in the bright 30 Doradus region (Visvanathan, 1966; Mathewson & Ford, 1970; Schmidt, 1976). More recently Wisniewski et al. (2007) mapped the optical polarization with imaging polarimetry around the NGC2100 supershell and showed regular polarization direction over the $15\arcmin$ field of view. Near-infrared polarimetric observations complement optical observations which are limited to very low interstellar extinctions. For instance, Kim et al. (2016) examined a $\simeq 1$ square degree region in the Northern part of the LMC (in Galactic coordinates) in $J,H,K$ bands and inferred, using the DCF method (Davis-Chandrasekhar-Fermi, Chandrasekhar & Fermi (1953); Davis (1951)), a mean field strength of around $3-25\,\mu G$ to be aligned with molecular clouds within that region.

At larger scales, synchrotron emission radio observations and Faraday rotation measures (RM) have been used to study the magnetic field of the LMC. Haynes et al. (1990) showed that the polarization direction to be roughly aligned with filamentary structures seen in HI near the 30 Doradus region. Gaensler et al. (2005) inferred the magnetic field strength to be larger than $1\mu G$ and to have a spiral shape with a random component that dominates the ordered component by a factor of $\simeq$ 3, based on RM. Mao et al. (2012) also analysed RM data and suggested the large-scale magnetic field to be spiral axisymmetric with the vertical component which reverses orientation across the disk and to be aligned with HI filaments. Still, it is likely that this structure applies to the large scale-height sampled by cosmic ray particles and does not probe the geometry of the field in the thin disk where most of the molecular material forming stars is located.

The large-scale structure of the magnetic field in the LMC associated with dense regions is still largely unexplored. Optical and near-infrared observations of dust polarization in the LMC are sparse, and radio observations bring information about the line-of-sight (LOS) magnetic field component located mostly beyond the thin disk. Here, we aim to complement the knowledge of the LMC magnetic field structure, taking advantage of the Planck satellite measurements at sub-millimeter wavelengths.

Thermal dust in the LMC has been shown to have striking differences from what is observed in the Milky Way (MW). In particular, Bernard et al. (2008) has shown the presence of excess emission near $70\micron$ and in the FIR. This was confirmed by the analysis of the Planck data towards the LMC and the SMC in Planck Collaboration et al. (2011), which showed much flatter thermal dust submm spectral energy distributions (SEDs) in those galaxies than in the MW. Many low metalicity galaxies seem to share a similar behaviour (e.g., Rémy-Ruyer et al., 2013). The origin of these differences are not clearly understood today. Some of the possible explanations invoke magnetic nanoparticles (e.g. Draine & Hensley, 2012) with specific predictions in polarization. It is therefore of interest to see if those differences also reflect in polarization, in order to gain knowledge on their possible origin.

The purpose of this study is to derive the dust polarization structure of the LMC using the Planck data. To achieve this goal, a key challenge is to estimate and subtract the foreground polarization. The emission of the interstellar dust grains is maximum at sub-millimetre wavelengths where the radiation is optically thin (Draine, 2003). This allows us to study distant objects, including those located outside of our Galaxy. However, a particular attention should be paid to quantifying the contribution of the foreground material to the total radiation even on the LOS out of the Galactic disk, because the polarised emission of the interstellar dust is observed not only in the disk of the Galaxy, where the dust is mixed with the molecular gas, but also at high Galactic latitudes ( $|\textit{b}_{II}|>30^{\circ}$ ). There, the emission of the diffuse dust and atomic gas, traced by far-infrared continuum and HI emission respectively, are observed to be correlated (Boulanger et al., 1996; Kalberla et al., 2021). Moreover, the polarisation fraction in the diffuse ISM is high compared to what is observed in dense regions: Planck Collaboration Int. XIX (2015) showed that the average polarization fraction at the column densities $10^{21}$ cm^-2 is around $7\%$ and can be up to $\simeq 20\%$ for lower column densities. Thus, it is necessary to account for foreground contamination in polarization data.

A first attempt to estimate foreground polarized extinction in the optical was made by Visvanathan (1966). Mathewson & Ford (1970) compiled different star catalogs and provided catalogs of starlight polarization towards the LMC. It was shown that most of the foreground stars show polarization oriented in a similar direction, which was later confirmed by the study of Schmidt (1976). There are two main difficulties when estimating the contribution of foreground material to the observed dust polarisation using stellar extinction measurements. First, the point-like stellar extinction measurements must be robustly extrapolated across the entire field with dust polarisation information. Second, the target stars are located at various depths within the dusty material that is responsible for the extinction, and hence sample different columns of interstellar material. In emission however, we would like to construct a fully sampled map of the foreground polarized emission, and that this map includes the contribution of all MW material along the LOS. Therefore, we test two techniques to infer maps of the polarization direction of the foreground emission which we compare to the measured stellar polarization directions, and we use the observed polarized emission around the LMC as a proxy for the polarization fraction of the foreground towards the LMC.

In order to explore the foreground polarization structure in this paper, we use two recent approaches that have established a link between the magnetic field direction and the properties of interstellar matter. The attractive feature of these approaches is that they use widely available data for the atomic gas distribution and kinematics to infer the magnetic field structure. Both techniques use fluctuations of intensity in channel maps, but analyze them differently.

The first technique is based on the notion of the magnetically aligned density filaments in cold HI data that can be detected using the Rolling Hough Transform algorithm (Clark et al., 2014). The latter is an algorithm designed to detect linear structures in two-dimensional images. Using this approach, Clark & Hensley (2019) constructed three-dimensional maps of the linear Stokes parameters using HI data and showed them to be statistically consistent with filaments in the Planck data, for which parallel alignment with the magnetic field in the diffuse ISM was previously reported (Planck Collaboration Int. XXXV, 2016; Alina et al., 2019).

The second technique employs a prediction of MHD turbulence theory (Goldreich & Sridhar, 1995a; Beresnyak & Lazarian, 2019) and the theory of turbulent reconnection (Lazarian & Vishniac, 1999) that rotation axes of magnetic eddies are aligned with the local direction of the magnetic field (Lazarian & Vishniac, 1999; Cho & Vishniac, 2000). Such aligned eddies induce the maximal variance of the velocity in the direction perpendicular to the eddy rotation axis. Therefore, the gradients of the velocity amplitude in MHD turbulence get maximal in the direction perpendicular to the local magnetic field, i.e. the gradient direction acts as a proxy of polarization. The corresponding Velocity Gradients Technique (VGT) has been applied not only to HI data, but also to CO, CHN and other molecular and atomic transitions, i.e. to the settings where the existence of filaments in real space has not been reported. These studies have shown a good correspondence between the magnetic field directions inferred using VGT and that obtained from polarimetry, e.g. with the Planck, BLASTPOL, and HAWC+ data (Hu et al., 2019; Alina et al., 2022; Hu et al., 2020, 2021; Liu et al., 2022). For the VGT, the contribution arising from the actual HI density structures provides an unwelcome interference. To minimize its effects, we employ a recently developed Velocity Decomposition Algorithm (VDA) approach (Yuen et al., 2021) that separates the velocity and density fluctuations within the spectral data cubes. The property of intensity gradients being related to the magnetic field has also been used in other studies, e.g., Koch et al. (2012) to infer magnetic field strength in molecular clouds.

In this paper, we use both techniques to estimate the contribution of the Milky Way foreground to the observed polarization signal. We apply the techniques to HI data obtained by the GASS survey McClure-Griffiths et al. (2009). In order to assess the accuracy of the magnetic field directions inferred using these techniques, we compare their results to polarization directions measured via optical polarization of stars by Mathewson & Ford (1970) and Schmidt (1976).

The paper is organized as follows. In Sect. 2 we describe data sets used in the study. In Sect. 3 we briefly present the two existing techniques for estimating the magnetic field direction and describe our method for foreground removal. We present and discuss our results in Sect. 4. In Sect. 5 we summarize our key conclusions and discuss future perspectives.

2 Data

2.1 HI data

Refer to caption — Figure 1: HI integrated intensity map, in K km s^-1, towards the LMC. It is overlaid with a drapery pattern showing the POS magnetic field direction derived from the Planck data at $16\arcmin$ , computed using the LIC technique. The blue, red, green and pink contours correspond to the HI integrated intensity levels of $1000,1400,2400$ , and $2800$ K km s^-1, respectively.

We use the third data release of the Galactic All-Sky Survey (GASS), which observed the 21-cm HI emission with the Parkes 64-m Radio Telescope between January 2005 and October 2006 (Kalberla & Haud, 2015; McClure-Griffiths et al., 2009). This data release has been corrected for stray radiation and cleaned for radio frequency interference. The data has an effective spatial resolution of 16.2’. The survey covers local standard of rest (LSR) velocities ranging from -495 km/s to 495 km/s, with a channel separation of 0.82 km/s. The total integrated intensity map towards the LMC is shown in Figure 1 with the magnetic field direction inferred from polarization in emission in the Planck data shown using the line integral convolution (LIC) technique (Cabral & Leedom, 1993).

To separate the foreground Galactic emission and the emission from the LMC, we use the velocity ranges proposed by Staveley-Smith et al. (2003) and use velocities from $-50$ km s^-1 to $90$ km s^-1 for the Milky Way, and from $182$ km s^-1 to $366$ km s^-1 for the LMC. Hereafter, the HI maps computed over the adapted Milky Way velocities range above will be referred to as the foreground maps. The foreground HI map is shown in Figure 2.

The total gas column density, used to provide the average column density further in Table 1 was calculated from Parkes HI and NANTEN CO (Fukui et al., 2008) data, assuming optically thin HI and CO conversion factor of $\mathrm{X_{CO}}=2\,10^{20}$ cm^-2 K km s^-1.

2.2 Planck data

The Planck¹¹1Planck (http://www.esa.int/planck) is an ESA science mission designed and completed by a collaboration of institutes directly funded by ESA Member States, NASA, and Canada. instrument mapped the whole sky over a range of frequencies from 70 to 857 GHz to measure the Cosmic Microwave Background (CMB) anisotropies and Galactic foregrounds. Here, we use the Planck PR3²²2The Planck data is available at https://pla.esac.esa.int/#home data release. The CMB contribution is subtracted using the Generalized Needlet Internal Linear Combination algorithm Remazeilles et al. (2011); Collaboration (2020). We also subtracted the Cosmic Infrared Background (CIB) monopole contribution to the intensity map according to Planck Collaboration XI (2014). For polarization, we restrict the analysis to the $353$ GHz channel, which is the highest frequency channel of Planck with polarization capabilities.

The LMC was observed by Planck many times, since the Planck scanning strategy (Planck Collaboration I, 2011) repeatedly covered the region close to the ecliptic poles. In particular, the LMC is located at the boundary of the Planck deep field. Thus, the variation of the noise level across the LMC is one order of magnitude in intensity, as illustrated in the left panel of Figure 21, and the signal-to-noise ratio (SNR) of the intensity within the LMC is $\sim$ 1000 (see right panel of Figure 21). Given that these observations were repeated with many different analysis angles, the Planck data in polarization at 353 GHz in the direction of the LMC are in principle extremely reliable.

The polarization fraction $p$ and angle $\psi$ can be derived from the Stokes parameters $I,\,Q,$ and $U$ observed by Planck according to following equations:

p=\frac{\sqrt{Q^{2}+U^{2}}}{I}\,,

(1)

\psi=0.5\mathrm{atan}(-U,Q)\,.

(2)

The negative sign in front of $U$ is to pass from the COSMO (or Healpix) convention used for the Planck data to the IAU convention for polarization. $p$ and $\psi$ are biased non-linearly due to noise propagation in equations 1 and 2 (Serkowski, 1958; Simmons & Stewart, 1985; Vaillancourt, 2006; Quinn, 2012; Montier et al., 2015a). Estimates of $p$ and $\psi$ must be corrected for this bias.

We compute the Mean Posterior estimates of $p$ and $\psi$ using the Bayesian method described in Montier et al. (2015b), which takes into account the noise variance matrix of each map pixel. To do so, we performed Monte-Carlo simulations with $10^{7}$ realizations. The advantage of the Bayesian method is that it allows us to compute the polarization angle and polarization fraction, and their uncertainties from a common probability distribution function (PDF), even at low SNR (SNR( $p$ ) $\simeq 2$ ). The resulting polarization fraction is shown in Figure 3. Other existing methods only provide estimates of the polarization fraction or angle that are valid at intermediate to high SNR, SNR( $p$ ) $>3$ (Simmons & Stewart, 1985; Montier et al., 2015a).

We brought the Planck data to the resolution of the HI data. To do so, we degraded the HEALPix³³3https://healpix.sourceforge.io/ Planck maps at nominal angular resolution ( $\simeq 5\arcmin$ ) to $16\arcmin$ and pixel size corresponding to $N_{side}=1024$ , according to the procedure described in Planck Collaboration Int. XIX (2015). This method takes into account the rotation of the local polarization reference frame in the convolution. The angular size of $16\arcmin$ at the distance of the LMC corresponds to a physical size of 200 pc.

The local maps used in this study were computed using the Drizzlib⁴⁴4http://cade.irap.omp.eu/dokuwiki/doku.php?id=drizzlib approach, which is designed to calculate the intersection of pixel surfaces precisely and hence preserves the photometry accuracy.

In this paper, we only derive polarization parameters of the LMC for the data that satisfy the following criteria: the original intensity is larger than the median value of $I_{med}=0.276$ MJy/sr inside the region shown in the map in Figure 3 and the signal-to-noise ratio of polarization fraction is larger than three, $p/\sigma_{p}>3$ , where $\sigma_{p}$ is computed using the Bayesian method described above. The resulting mask is shown as the black contour in Figure 3.

2.3 Optical polarization data

The optical polarization data used are taken from following catalogues:

$-$

The Schmidt (1976) (hereafter S76) catalogue contains optical polarization data towards the Magellanic Clouds (MC) region. The polarization observations for 309 stars were obtained between 1971 and 1972, using the ESO 1m photometric telescope and the ESO two channel polarimeter at La Silla observatory in Chile. The catalogue also contains data from Schmidt (1970) obtained from the Boyden observatory. In total for the LMC, there are 149 galactic foreground stars. We manually copied the polarization data for foreground stars from Table 6 of Schmidt (1976) original paper.
$-$

Mathewson & Ford (1970) (hereafter MF70) catalogue contains optical polarization measurements taken during 1969 and 1970 using the 74-inch reflector at Mount Stromlo and the 24-inch telescope at Siding Spring Observatory. There are 31 foreground galactic stars toward the LMC.

The polarization angles in the above catalogues were originally given with respect to the equatorial coordinate system at epoch B1950, which we precessed to J2000.

In our star catalogs, distance estimates were missing or doubtful for a significant number of stars. As it is critical to control how far the Milky Way ISM is traced by the observations, we updated the distances and completed missing data using the SIMBAD database (Wenger et al., 2000). The new distances were taken from the following Gaia-based (Gaia Collaboration et al., 2016) catalogues: the Gaia 2nd release (Gaia Collaboration, 2018; Bailer-Jones et al., 2018) and the StarHorse catalog of astro-photometric distances (Anders et al., 2019). In total, we compiled distance data for 145 of the 149 stars in the S76 catalog, and 30 of the 31 stars in MF70 catalog. In case of discrepancies between the original polarization catalogues and recent distance estimated, we adopted the latter. This represents 97.3 $\%$ and 96.8 $\%$ of entries in the S76 and MF70 catalogs respectively. The compiled data is publicly available upon publication⁵⁵5https://github.com/danakz/starpol. We used background stars for which the distance is larger than $500$ pc in order to ensure that the corresponding polarization traces the foreground magnetic field direction over most of the corresponding emitting LOS. We also retained stars for which the uncertainty on the polarization angle is less than $20^{\circ}$ . These criteria yield a final sample of 56 foreground stars.

The magnetic field position angles inferred from the stellar extinction data are shown in Figure 2. The foreground magnetic field directions appear to be rather uniform, oriented almost vertically and close to the North-East to South-West direction. In Figure 4, we present histograms for all stars and for our final selection. We observe that restricting the selection to star distances larger than $500$ pc does not significantly change the pattern. However, the corresponding histogram is narrower. The direction of the visible polarization is globally aligned parallel to the magnetic field structure around the LMC, traced by Planck. Also, the HI background map shows a clear elongated structure crossing the LMC. This structure dominates the overall foreground, and actually extends further than shown in Figure 2, with a similar orientation. The magnetic field direction inferred from our stellar data appears to globally follow that direction, which is consistent with the expectation of a magnetic field oriented parallel to density filaments (Clark & Hensley, 2019). This motivates using filaments identification to constrain the magnetic field direction.

3 Methods

The separation between the polarization originating from the LMC and that from the MW foreground is best achieved through the Stokes parameters formalism. Given that both the MW foreground and the LMC emission is optically thin, we assume

$\displaystyle I_{LMC}$	$\displaystyle=$	$\displaystyle I_{obs}-I_{fg}$	(3)
$\displaystyle Q_{LMC}$	$\displaystyle=$	$\displaystyle Q_{obs}-Q_{fg}$	(4)
$\displaystyle U_{LMC}$	$\displaystyle=$	$\displaystyle U_{obs}-U_{fg}\,.$	(5)

where $I_{obs}$ , $Q_{obs}$ and $U_{obs}$ are the total Stokes parameters observed with Planck and the corresponding quantities with subscripts _LMC and _fg are for the LMC and MW foreground respectively. The foreground Stokes parameters can be written as following:

	$\displaystyle Q_{fg}=p_{fg}\,I_{fg}\,\cos{2\psi_{fg}}$		(6)
	$\displaystyle U_{fg}=p_{fg}\,I_{fg}\,\sin{2\psi_{fg}}\,.$		(7)

where $I_{fg}$ is the foreground total intensity and $p_{fg}$ and $\psi_{fg}$ are the assumed polarization fraction and angle of the foreground contribution respectively, which we estimate using the approaches described below.

3.1 Determination of the foreground polarization angle

To construct a map of the foreground polarization angle, we use two independent techniques that assume a relationship between the magnetic field direction and the intensity structures identified using spectroscopic data. The first uses the Rolling Hough Transform (RHT) to derive the direction of magnetically-aligned filaments. The second uses the Velocity Gradient Technique (VGT) to constrain the field direction based on MHD assumptions. We apply the methods to the HI data that accounts for the foreground contribution, as described in Sect. 2.1, that is from $-50$ to $90$ km s^-1. We assess the performance of each method to this study by comparing their predictions with the magnetic field direction derived from optical starlight polarization data, assumed to be a robust measurement of the projected field orientation on those LOS where stellar data is available.

3.1.1 The Rolling Hough Transform

Comparisons between the atomic hydrogen emission at 21 cm and magnetic field directions derived from polarization have shown that the diffuse HI gas is organized into elongated linear fibers, the orientation of which correlates with the orientation of the magnetic field (Clark et al., 2014). Based on this observation, Clark & Hensley (2019) suggested that it is justified to derive the magnetic field direction via the determination of the orientation of HI fibers. While a detailed description of the RHT can be found in the original article, here we briefly describe the main steps involved in the method.

RHT is a machine vision algorithm that identifies and quantifies linear structures in two-dimensional maps. The input two-dimensional image is smoothed using a top-hat function of a given diameter defined by the user. The smoothed image is then subtracted from the original map and a threshold is applied to the difference map to obtain a bitmask. Next, a disk of a given diameter is extracted around each point in the bitmask and the Hough Transform is applied to each disk. The Hough Transform measures the intensity of each rectangular kernel centered at the disk and stores the one above set threshold as $R(\theta,x_{0},y_{0})$ , RHT intensity, as a function of angle for each pixel. The RHT algorithm takes as input parameters a map, over which to operate, the size and width of the kernel window, smoothing radius, and the intensity threshold as a percentage and returns a map of identified structures, their position angle values ( $\theta_{\text{RHT}}$ ) and associated uncertainties ( $\sigma_{\theta,\text{RHT}}$ ). One of the main advantages of this technique is that it determines the linearity of structures that do not have sharp boundaries. A disadvantage of the method is that the size and the shape of the identified structures strongly depends on the user-defined kernel parameters.

Clark & Hensley (2019) provided a full sky data set⁶⁶6https://dataverse.harvard.edu/dataverse/ClarkHensley of $I,\,Q$ and $U$ in 3D (the third dimension being the velocity) of the filamentary structures identified using the RHT method applied to HI data from the GALFA (Peek et al., 2018) and the HI4PI (HI4PI Collaboration et al., 2016) surveys over the Galactic velocity range. We extract the data from the integrated RHT $Q$ and $U$ maps towards the LMC region. The left panel of Figure 5 shows the orientations of the structures derived from the RHT method as a LIC-drapery and the starlight polarization data with the black segments. Although those orientations seem to globally match structures that are observed in the HI integrated intensity map, the directions of the detected filaments vary a lot. The RHT filaments directions do not correspond to the starlight polarization data very well as can be seen from the left panel of Figure 6 where a histogram of the angle difference is represented with black solid line. This may be due to some of the following reasons. First, the RHT parameters are adjustable and those used in Clark & Hensley (2019) can be sub-optimal for the region. In fact, the size of the kernel was fixed to $75^{\prime}$ for the calculation for the whole Galaxy. Second, in the original work, the RHT procedure was applied in velocity channel bins, with the bin size fixed to $\Delta v=3.7$ km s^-1. The integration of the RHT Stokes $Q$ and $U$ over all the velocity channels, including channels with no obvious structures or less significant small features could yield many various directions to appear.

Even though Clark & Hensley (2019) found a good correlation between the orientations of RHT-derived structures and the magnetic field direction derived from the Planck data across the ISM in our Galaxy, the corresponding maps do not seem to us to be the most appropriate in the direction towards the LMC.

We tried to improve the HI-filaments output to correspond better to the observed polarization and to address possible issues regarding the variation of the direction of filaments detected in different velocity-channels and the choice of the RHT parameters. To do so, instead of detecting filaments in thin velocity channels and then integrate over the velocities, we ran the RHT method over the integrated foreground HI map. The kernel in the RHT technique is an input parameter, and we tried various kernel sizes. The best detection of HI intensity filaments that reflected visual perception of the HI map was achieved for the kernel with the length of $2^{\circ}$ and the width of $25\arcmin$ . It is worth noting that we did not introduce any matching of the filaments with the polarization data in our optimization of the RHT parameters. Since filaments are only detected in the brightest regions, there are parts on the map that have missing values. We therefore interpolated the results via the Stokes parameters $Q_{RHT}=\cos(2\theta_{RHT})$ and $U_{RHT}=\sin{2\theta_{RHT}}$ using a spatial exponential interpolation:

X^{\prime}_{i}=\frac{\sum\limits_{j}^{N}X_{RHT,j}\mathrm{e}^{-\frac{d_{ij}^{2}}{\sigma_{d}^{2}}}}{\sum\limits_{j}^{N}\mathrm{e}^{-\frac{d_{ij}^{2}}{\sigma_{d}^{2}}}}\,,

(8)

where $X^{\prime}_{i}$ is the interpolated $Q$ or $U$ parameter, $d_{ij}$ are the distances from a given pixel $j$ to pixels where RHT filaments were detected, $X_{RHT,j}$ with $j\in[1,N]$ are the Stokes parameters of these filaments and $\sigma_{d}$ is a characteristic length, taken to be $100\arcmin$ , N is the number of pixels. We also smoothed the integrated foreground HI map at 20 pixels scale prior to the application of the RHT kernel. This gives the same spatial averaging as the one required in VGT calculation shown in the next Section and allowed us to compare the results. The right panel of Figure 5 shows the directions of the derived filaments as well as the interpolated directions as a LIC-drapery pattern and the blue contours show the extent of the detected filaments. The resulting map of RHT-directions using this approach seem to better follow the orientation of the observed elongated structures in the foreground HI map. However, the histogram of the difference with respect to the starlight polarization in the left panel of Figure 6 shows that the peak is shifted to 25 degrees. Also, the distance between the detected filaments in some cases is as large as the largest filament structure, and the spatial interpolation may be doubtful. Thus, although there is some parallel alignment between filaments and starlight polarization, there are also large differences for both RHT-based results.

3.1.2 Velocity gradients technique

Interstellar gas is turbulent (Armstrong et al., 1995; Chepurnov et al., 2010) and the theory of magnetohydrodynamic (MHD) turbulence (Beresnyak & Lazarian, 2019) connects the turbulent motion of the plasma to the magnetic field structure. On this basis, the velocity gradients technique (VGT) was developed specifically to trace interstellar magnetic fields (González-Casanova & Lazarian, 2017; Yuen & Lazarian, 2017; Lazarian & Yuen, 2018).

The foundations to the technique are rooted in the turbulence anisotropy discovered in Goldreich & Sridhar (1995b). Turbulent reconnection enables mixing of the magnetic field lines by eddies that are elongated along the local direction of the magnetic field (Lazarian & Vishniac, 1999). Within these eddies, magnetic field lines reconnect and do not prevent eddies to evolve in a hydrodynamic fashion. Since the eddies are closely aligned with the direction of the magnetic field that they are embedded in, the changes of amplitudes of velocities are maximal in the direction perpendicular to the magnetic field. Therefore the gradients of the velocity amplitude are expected to be perpendicular to the magnetic field direction and can trace the magnetic field orientation, in a way similar to magnetic field orientation being traced by the polarization of the dust emission (see also a pictorial illustration in Pattle (2019)).

The velocity fluctuations are reflected in the fluctuations of intensity in spectroscopic data cubes and the latter are most pronounced in thin channel maps (Lazarian & Pogosyan, 2000). The analysis of gradients of intensity in channel maps are suggested to use in the version of the Velocity Gradient Technique (VGT) that we employ in the paper (Lazarian & Yuen, 2018). It is clear from the description that while the RHT is focused on structures arising from density variations, the VGT is focused on exploring the effects induced by velocity fluctuations.

It is worth noting that additional effects can influence velocity gradients. The propagation of shocks and the effect of gravitational collapse were considered in the literature as possible influencing factors. Shocks are most important for gradients of intensities (Yuen & Lazarian, 2017; Hu et al., 2019). For velocity gradients, the effects of shocks is subdominant as the change in velocity is very localized and the sub-block averaging employed in the technique, and described further in this Section, suppresses the effect of such localized velocity amplitude changes. The velocity gradients arising from gravity can be important for localized regions where gravity-induced acceleration is larger than that induced by turbulence (Hu et al., 2019) In practical terms, these are relatively small collapsing regions of molecular clouds and are not relevant for the present study.

The effects of gravity and shocks can be amplified if the employed measures of velocity are contaminated by density effect. This is more pronounced when using velocity centroids (Esquivel & Lazarian, 2005), rather than when using individual channel maps (Lazarian & Pogosyan, 2000). Thus, in the present paper we use the version of the VGT that employs velocity channels (Lazarian & Yuen, 2018) rather than maps of velocity centroids (González-Casanova & Lazarian, 2017; Yuen & Lazarian, 2017).

More recently, on the basis of the theory in Lazarian & Pogosyan (2000), Yuen et al. (2021) proposed a new technique, the Velocity Decomposition Algorithm (VDA), for filtering effects of density in spectroscopic data. This approach significantly reduces the effects of density fluctuations and was shown to reduced the area over which the measured gradients of intensities in velocity channels are not perpendicular to the magnetic field orientation. In our study we employ the VDA technique to increase the accuracy of the VGT method.

In practice, gradients are calculated in velocity-centroid or velocity channel maps for every pixel. Then, a histogram of the gradients is constructed for a sub-block of a given size. The size of the sub-blocks depends on the noise in the data and is chosen to get a distribution of gradient directions close to Gaussian. In this study, the conventional size of 20 $\times$ 20 pixels (Yuen & Lazarian, 2017) is employed, which corresponds to a total area of around 50 beams. The expectation value of the distribution gives the mean orientation of the gradients within a sub-block (Yuen & Lazarian, 2017), while the width of the distribution provides a measure of media magnetization, i.e. its Alfven Mach number, $M_{A}=V_{L}/V_{A}$ , where $V_{L}$ is the turbulent injection speed and $V_{A}$ is the Alfven speed (Lazarian et al. 2018). In this paper, we focus on the magnetic field directions derived from velocity gradients and do not study the $M_{A}$ distribution.

The fitting of a Gaussian into the distribution of gradient directions provides a practical advantage of the VGT, i.e. the size of the sub-block is being determined on the basis of the objective criterion. The obtained VGT maps do not require extrapolation either.

Our predicting of the galactic foreground polarization is similar to that in Lu et al. (2020), where the motivation was to improve the accuracy of distilling the Cosmic Microwave Background (CMB) polarization. The difference is that we use the VDA to remove the effects of densities and thus increase the accuracy of magnetic field tracing.

The magnetic field orientation of the foreground as derived using by VDA method is shown in Figure 7. The orientation is globally in the North-East to South-West direction, close to the vertical, and the starlight extinction polarization is generally aligned parallel to this orientation. As can be seen from the right panel of Figure 6, the histogram of the angles difference between the predicted magnetic field direction and the observed polarization in the visible is centered on zero, indicating a good match between the predictions of the method and the observed POS magnetic field orientation. In the following, we use the result from the VDA method to constrain the direction of the foreground polarization.

3.2 Determination of the foreground intensity

We estimate the Galactic foreground intensity $I_{fg}$ , assuming that the MW dust emission mostly originates from optically thin atomic gas and that dust properties remain uniform within the chosen field of view. This leads to the expectation that

I_{fg}=N_{H}\times\alpha\,\,

(9)

$\alpha$ is a constant and the column density $N_{H}$ is computed from the HI data using

N_{H}=1.823\times 10^{18}\int T_{B}dv\,,

(10)

with the integral carried over Galactic velocities. We correlated the gas column density from HI data with the original Planck intensity at 353 GHz in the outer region and found a slope of $\alpha=1.87$ ( $\times 10^{21}$ ). We used this value to estimate the foreground intensity over the LMC using Eq. 9. The resulting foreground intensity map is shown on Figure 8.

Figure 9 shows the histogram of the foreground intensity values towards the LMC at 353 GHz as well as the histogram of intensity values of the LMC after foreground removal. The average LMC brightness in this region is 0.68 MJy/sr, to be compared to the average foreground intensity of 0.29 MJy/sr. The foreground emission is therefore generally significantly lower than that of the LMC in total intensity. The resulting intensity of the LMC after foreground removal is shown in Figure 22, where the outer spiral density structures and central molecular regions can clearly be distinguished.

3.3 Determination of the foreground polarization fraction

The polarization fraction map derived from the Planck data at 353 GHz towards the LMC is shown in Figure 3. We measured the average value of $p$ around the LMC from this map, excluding the area delimited by the red contour around the LMC, shown on the figure. The red contour is built on the basis of a circle of $5^{\circ}$ radius around the center of the LMC. In addition, we excluded areas where SNR( $p$ ) $<\,3$ to avoid bias. The resulting average polarization fraction is $p_{fg}=7.2\%$ , which we assume constant over the whole LMC surface to apply to the foreground emission. This value is close to the value derived for the MW from the Planck data, on average, which is around $7\%$ (Planck Collaboration Int. XIX, 2015). Although the polarization fraction of the emission of the foreground material probably varies, we believe this approach to be justified as it does not privilege any particular region and the subtraction is done uniformly. This hypothesis is also supported by the globally uniform foreground magnetic field orientation, as derived from the gradients and observed from the foreground stars polarization.

4 Results and Discussion

We used the quantities derived in the previous section and Eqs. 6 and 7 to compute maps of the Stokes parameters for the foreground emission $Q_{fg}$ and $U_{fg}$ . These maps were then subtracted from the Planck 353 GHz data. We also compute the noise variance matrix on the derived Stokes parameters, as detailed in Appendix B. This allowed us to derive the SNR on the derived LMC polarization fraction. We represent in Figure 10 the variation of the median polarization fraction in the area defined earlier for the LMC as a function of SNR(p_LMC). We observe a drastic decrease of p for SNR(p_LMC) for values lower than $0.5$ . Thus, the new mask is taken above this threshold and is used for the remainder of the analysis.

Below we first describe the resulting polarization angle which gives the large-scale structure of the LMC. Second, we discuss the polarization fraction properties.

4.1 Large-scale magnetic field structure

Figure 11 shows the POS magnetic field direction derived from the Planck data at 353 GHZ after foreground removal. Although the initial and the foreground subtracted magnetic field directions are similar over the brightest LMC regions, as can be seen from Figure 22, the foreground removal produces a rotation of the magnetic field directions over faint regions the LMC, especially in the north-east quadrant. The histogram of the angles difference before and after the foreground removal is shown on Figure 12. We observe that the angles differ by at most $20^{\circ}$ at the $3\sigma$ level. The alignment between intensity structures of the LMC and the magnetic field direction appears tighter after foreground removal than in the initial data, in particular, in the southern parts of the LMC. We also observe that the magnetic field of the LMC after foreground removal tends to align with the outer faint arms. Moreover, we observe an alignment between the faint inner structures and the magnetic field, south-east from the center. On the contrary, in a part of the north-western region of the LMC, we still observe the regular north-west to south-east magnetic field orientation. It can be due to an incomplete foreground removal. However, such a direction may also be inherent to the LMC’s magnetic field. Schmidt (1976) retrieves similar polarization direction for stars located in the LMC, after their foreground removal performed for the visible polarization data (seen in the Figure 3 of Schmidt (1976)).

Figure 13 shows the map of the polarized intensity $P=pI=\sqrt{Q^{2}+U^{2}}$ after foreground removal. The figure shows a faint elongated structure in the outskirts of the southern part of the LMC, which is only clearly apparent in the polarized intensity map. It corresponds to Arm S introduced by Staveley-Smith et al. (2003). This feature is clearly seen in the integrated intensity from HI data, as shown by Mao et al. (2012). However, no peculiar structures in synchrotron polarized intensity are observed in that region. The polarized dust intensity highlights regions with intrinsically high polarized dust emission which allow to detect aligned dust in the low density Arm S of the LMC for the first time. There, the magnetic field shows a smooth spiral structure from the eastern to the southern part of the LMC.

We use the polarization angle dispersion function $\mathit{S}$ as a measure of the regularity of the magnetic field, given by

S=\sqrt{\frac{1}{N(l)}\sum_{i=1}^{N(l)}\left[\psi(\mathbf{x})-\psi(\mathbf{x}+\mathbf{l}_{i})\right]^{2}}\,,

(11)

where $\psi(\mathbf{x})$ and $\psi(\mathbf{x}+\mathbf{l}_{i})$ are the polarization angles at locations $\mathbf{x}$ and $\mathbf{x}+\mathbf{l}_{i}$ , respectively, and the sum is taken over all pixels contained within an annulus of radius $l=|\mathbf{l}|$ (called the lag) and width $\delta\,l$ from position $\mathbf{x}$ . $\mathit{S}$ measures the average amplitude of the local variations of the polarization angle at scale $\mathit{l}$ (Serkowski, 1958; Hildebrand et al., 2009; Alina et al., 2016).

The maps of $\mathit{S}$ at a lag of $16\arcmin$ before and after foreground removal are shown on Figure 14. We restrict the analysis to the data with SNR(S) $>3$ , where the uncertainty is computed using Monte Carlo simulations that take into account the whole noise covariance matrix. The difference between the maps reflects the difference between the initial and foreground corrected magnetic field directions shown on Figure 11. The foreground removal reveals noticeable rotation of the magnetic field especially correlated with locations of density structures. Planck Collaboration Int. XIX (2015) showed an anti-correlation between $p$ and $\mathit{S}$ in the Milky Way which indicated that, globally, on large scale the decrease of polarization is due to the rotation of the magnetic field angle. We represent on Figure 15 the scatter plot of $p$ versus $S$ for the LMC. We observe a clear anti-correlation trend, similar to what has been observed for the Milky Way in general in the Planck data (Planck Collaboration Int. XIX, 2015). Following the approach used in Planck Collaboration Int. XIX (2015), we fit a straight line to the data in log-log, using

\log_{10}(S)=\alpha\times\log_{10}(p)+\beta\,.

(12)

The fit to the original data leads to $\alpha=-0.789$ and $\beta=1.661$ while the same fit on the foreground subtracted data leads to $\alpha=-0.504$ and $\beta=1.287$ . The slope of the correlation before the foreground removal is therefore closer to the value derived for the ISM in the Milky Way, $\alpha=-0.834$ , in Planck Collaboration Int. XIX (2015) while it is shallower after the foreground removal.

4.2 Relative orientation between magnetic field and the LMC structure

The relative orientation between the magnetic field and the galaxy structure, which is almost seen face-on, provides a hint on the role of the magnetic field in the global dynamics of the LMC. Although the visual inspection of Figure 11 allows us to conclude that in vast regions such as the arms of the galaxy the magnetic field is aligned with density structures, we aim to quantify more precisely their relative orientation. To do so, the RHT algorithm was applied to the Planck 353 GHz intensity map, and the best-suited parameters of the kernel’s width $16\arcmin$ and length $50\arcmin$ were chosen based on the visual analysis. We represent the corresponding bitmask map on Figure 16. The absolute difference between the magnetic field angles and the position angles of the detected filaments in Figure 17 shows that there is a trend for alignment between the magnetic field direction and the density structures in the LMC.

4.3 Polarization fraction and polarized intensity

The map of polarization fraction $p$ after foreground polarization removal is shown in the left panel of Figure 18 and the corresponding SNR is shown in the right panel. Globally, $p$ varies between 2 $\%$ and 12 $\%$ , and in some regions reaches around 20 $\%$ . We also observe values that reach around 25 $\%$ , in particular at the location of the foreground structure detected in the HI map. Although these values may be real, they may be due to noise or to the imperfect foreground subtraction. Figure 19 shows the histogram of $p$ after foreground removal for the pixels within the mask defined in Sect. 3.3 and Figure 18. We observe that the foreground removal broadens the distribution. The median value is $p_{med}=$ $3.6$ $\%$ while the original median polarization is $4\%$ . Figure 20 shows how the median value varies with increasing column density. As expected, as the column density increases, corresponding to increasing intensity values, the polarization fraction decreases. More details can be found in Table 1 where we report statistics of the polarization fraction values. We follow the approach by (Planck Collaboration, 2018) and show in Table 1 the upper 90th, 95th and 99th percentiles. We take the 99th upper percentile as the conservative estimate of the lower limit of the maximum polarization fraction that is observed in the most favorable conditions, that is, when the magnetic field is oriented in the plane of the sky over most of the LOS. This value is $p_{\mathrm{max}}$ = $21.6$ $\%$ . As the LMC is indeed viewed almost face-on, we have here an indication on the maximum polarization observed by Planck which is limited by beam dilution and LOS averaging. This value is close to the corresponding values derived for the Milky Way of $19.8\%$ (Planck Collaboration Int. XIX, 2015) and the value of $22\%$ (Planck Collaboration, 2018) derived from the Planck data at the resolutions of $60\arcmin$ and $80\arcmin$ respectively. We note that the average gas column density of the LMC estimated inside the mask is $1.7\,10^{21}$ cm^-2 which is even higher than the average value of the MW, around $1\,10^{21}$ cm^-2 (Planck Collaboration Int. XIX, 2015). Planck evidenced a systematic decrease of the polarization fraction with column density (Planck Collaboration Int. XIX, 2015), and the difference in column density would lead to lower polarization fractions. However, this is not observed here for the LMC. The difference may indicate that the intrinsic polarization by dust particles is actually higher than the value derived from galactic regions of the solar neighborhood. It could also indicate that the magnetic field in the LMC is more ordered than that of the MW.

Vandenbroucke et al. (2021) simulated polarimetric observations of external galaxies, including dust properties, alignment, radiative transfer and spiral magnetic field geometries with different inclination angles, based on the MHD structure predicted by the simulations of Milky Way type galaxies. The latter simulations were performed in the frame of the Auriga project (Grand et al., 2017). The model, when observed from within the simulated galaxy, was constrained using the Planck data and was used to derive linear polarization fractions that could have been observed with the SPICA/B-BOP instrument. When observed from a distance of 10 Mpc with the size of 50 kpc, the angular resolution of the SPICA/B-BOP instrument ( $38"$ ) corresponds to the Planck resolution in the LMC. Their derived value of the lower limit for the maximum polarization fraction at this resolution is around $13\%$ and the average value of the polarisation fraction is around $6\%$ . As stated by Vandenbroucke et al. (2021), to determine the maximum polarization fraction of an external galaxy, it should be resolved with at least 1 kpc, and the LMC observations with Planck largely satisfy this requirement. The value we derive here for the polarization fraction of the LMC are therefore lower than the expectations for the average value and significantly higher for the maximum value.

In order to derive an approximate predicted value of polarization fraction that would be observed if the LMC was an unresolved source, we average original $I,\,Q$ and $U$ (Mean Bayesian estimates) inside a circle of a radius of $5^{\circ}$ and obtain $1.1\%$ while the average obtained using the foreground removed parameters, is $3.3\%$ .

p [ $\%$ ]	LMC	obs	MW
$\mathrm{max(PDF)}$ [ $\%$ ]	2.5	3.9
median $p$ [ $\%$ ]	$3.6$	4.0	5.3¹
mean p [ $\%$ ]	$4.9$	4.2
upper 90 perc. [ $\%$ ]	9.8	6.7	12.68²
upper 95 perc. [ $\%$ ]	13.2	7.8	15.0²
upper 99 perc. [ $\%$ ]	$21.6$	10.0	18.7²

Table 1: Statistics of polarization fraction

p

in percent after foreground removal, over the pixels enclosed in the black contour shown in Figure 3, corresponding the description in column 1. These are the polarization fraction at the maximum of the PDF, the median and mean values and the upper 90, 95 and 99 percentiles. Column 2 gives the quantities as measured in the foreground subtracted maps. Column 3 gives the corresponding values in the original Planck data. Column 4 gives the corresponding values for diffuse Milky Way ISM as published in the Planck papers.

1

Data taken from Planck Collaboration Int. XIX (2015).
2

Data taken from Planck Collaboration (2018), where resolution $80\arcmin$ was employed. Fiducial offset values are taken, with the bias removed as indicated in Planck Collaboration (2018).

5 Summary and perspectives

The aim of this work was to isolate the polarization signal from the Large Magellanic Cloud (LMC) in the Planck polarization data at 353 GHz. The main task consisted in estimating the foreground polarization signal from the Milky Way and subtracting it. This required an estimating the foreground total intensity, polarization direction and polarization fraction.

We tested two techniques to determine the foreground magnetic field direction using HI data. The first technique is based on magnetically-aligned filaments paradigm (Clark & Hensley, 2019). It is based on observations in diffuse media density filaments are aligned with plane-of-the-sky (POS) magnetic field, while in dense medium their relative orientation is mainly perpendicular (Clark et al., 2014; Planck Collaboration Int. XXXV, 2016; Alina et al., 2019; Xu et al., 2019). The interpretation of the striations observed in the velocity channel maps according to (Clark & Hensley, 2019) is due to cold HI filaments aligned by magnetic field.

The second technique is based on MHD-theory, according to which the turbulent velocities induce the striations in velocity channel maps even in the absence of the density filaments. These fluctuations reflect the gradients of velocity amplitude which are perpendicular to the local magnetic field. Therefore the gradients of intensities in channel maps, similar to polarization, represent the underlying magnetic field direction (Lazarian & Yuen, 2018).

We compared results from the two techniques to the starlight extinction polarization data, for which we updated distance estimates using the Gaia data. This allowed us to choose the most distant stars that trace foreground ISM volumes close to what is traced by sub-millimeter emission in polarization. A practical advantage of these techniques is that they provide a map, while starlight measurements are sparse, which complicates further pixel-to-pixel analysis.

The comparison showed a better agreement for the velocity gradients technique, which we adopted for the rest of the analysis. The correspondence between the magnetic field direction inferred from velocity gradients and starlight polarization makes us confident that the derived pattern indeed represents the magnetic field structure.

We estimated the foreground intensity using properties of Galactic dust in the vicinity to the line-of-sight to the LMC. We assume that Galactic dust properties in the foreground towards the LMC match the average properties in our Galaxy around the LMC. We derived the relationship between dust optical depth and gas column density empirically from the maps of dust emission and HI data. As for the foreground polarization fraction, we assumed it to be constant across the LMC and equal to the average value observed around the LMC.

We estimated and removed foreground contribution from polarization signal in the Planck data and derived polarization parameters of the LMC. Although the observed Planck signal is dominated by the LMC, the foreground removal allows us to derive global characteristics of magnetic field structure and dust properties of the LMC and compare them to the Milky Way. We note that there are caveats to the method that we propose. First, the VGT technique requires sub-block averaging which provides the coarse grading of magnetic field maps as compared to the resolution of the original HI maps. Second, for the determination of the foreground Stokes $Q$ and $U$ parameters, we assume the polarization fraction of the Milky Way’s dust in the direction of the LMC to be uniform, which in fact can vary. However, we do not expect large variations of $p$ because of the slight variations of the angular dispersion function due to the overall uniform magnetic field direction. Nevertheless, our method allows us to estimate the foreground polarization in a conservative manner.

We also stress that we focus here on large-scale LMC properties, due to the low resolution of data used. We observe that magnetic field and density structure elongation directions follow each other closer after foreground removal. In particular, magnetic field and the south-eastern spiral arm (arm S) are well aligned. We also observe a closer alignment of the LMC structures with the magnetic field in the eastern and western outer arms, but also in the inner regions. However, in a fraction of the observed galaxy, that is the central part of the northern region, the large-scale magnetic field direction persists after foreground removal and is oriented north-east to south-west. The polarization angle dispersion function shows an anti-correlation with the polarization fraction as it is observed in the Milky Way, with a shallower slope, which indicates that depolarization is mostly due to the magnetic field rotation inside the beam and along the LOS.

A conservative estimate of maximum observed polarization fraction is $21.6$ $\%$ which is similar to the value derived for the Milky Way. We obtained a median value of $p$ of $3.6$ $\%$ and an average of $4.9$ $\%$ . We use results from recent simulations by Vandenbroucke et al. (2021) to show that this is lower than the expected value for a galaxy similar to ours observed faced-on at the distance of the LMC. The maximum polarization value we observe is however significantly higher than the expected value from those simulations ( $13\%$ ). The difference could originate from a more ordered large-scale magnetic field structure or from more efficient dust alignment in the LMC compared to our Galaxy. Both alternatives appear plausible. The difference in the magnetic field structure could arise from the difference in the star formation in the LMC and the MW while the difference in the dust alignment efficiency may arise from the different dust composition that have been evidenced by observations of the dust emission (Clayton & Martin, 1985; Matsuura et al., 2009; Gordon et al., 2014; Planck Collaboration et al., 2011; De Vis et al., 2017). The relative role of these factors requires further investigation, which is out of the scope of this paper.

Acknowledgements

This work was supported by the Science Committee of the Ministry of Education and Science of the Republic of Kazakhstan, in the frame of the project "Study of the polarised emission of the interstellar medium of the Magellanic Clouds using big data analysis and machine learning", Grant No. AP08855858. DA and MI acknowledge the Nazarbayev University Faculty Development Competitive Research Grant Programme No110119FD4503. AL and KHY acknowledge the support of NASA ATP2303 80NSSC20K0542 and NASA TCAN 144AAG1967 as well as NSF AST 1715754 and 1816234.
We acknowledge the use of data provided by the Centre d’Analyse de Données Etendues (CADE), a service of IRAP-UPS/CNRS (http://cade.irap.omp.eu).
We thank the anonymous reviewer for the time and effort invested into the review of the manuscript and for the valuable comments and suggestions.

Data availability

Upon publication, maps of the foreground-subtracted $I,\,Q$ and $U$ at 353 GHz will be made available via the CADE⁷⁷7http://cade.irap.omp.eu database.

References

Alina et al. (2016) Alina D., Montier L., Ristorcelli I., Bernard J.-P., Levrier F., Abdikamalov E., 2016, A&A, 595, A57
Alina et al. (2019) Alina D., Ristorcelli I., Montier L., Abdikamalov E., Juvela M., Ferrière Bernard J.-P., Micelotta E., 2019, MNRAS, 485, 2825
Alina et al. (2022) Alina D., et al., 2022, A&A, 658, A90
Anders et al. (2019) Anders F., et al., 2019, VizieR Online Data Catalog, p. I/349
Armstrong et al. (1995) Armstrong J. W., Rickett B. J., Spangler S. R., 1995, ApJ, 443, 209
Bailer-Jones et al. (2018) Bailer-Jones C. A. L., Rybizki J., Fouesneau M., Mantelet G., Andrae R., 2018, AJ, 156, 58
Beresnyak & Lazarian (2019) Beresnyak A., Lazarian A., 2019, Turbulence in Magnetohydrodynamics
Bernard et al. (2008) Bernard J.-P., et al., 2008, AJ, 136, 919
Boulanger et al. (1996) Boulanger F., Abergel A., Bernard J. P., Burton W. B., Desert F. X., Hartmann D., Lagache G., Puget J. L., 1996, A&A, 312, 256
Cabral & Leedom (1993) Cabral B., Leedom L., 1993, in Special Interest Group on GRAPHics and Interactive Techniques Proceedings. pp 263–270
Chandrasekhar & Fermi (1953) Chandrasekhar S., Fermi E., 1953, ApJ, 118, 113
Chepurnov et al. (2010) Chepurnov A., Lazarian A., Stanimirović S., Heiles C., Peek J. E. G., 2010, ApJ, 714, 1398
Cho & Vishniac (2000) Cho J., Vishniac E. T., 2000, ApJ, 538, 217
Clark & Hensley (2019) Clark S., Hensley B. S., 2019, ApJ, 887, 136
Clark et al. (2014) Clark S., Peek J., Putman M., 2014, The Astrophysical Journal, 789, 82
Clayton & Martin (1985) Clayton G. C., Martin P. G., 1985, ApJ, 288, 558
Collaboration (2020) Collaboration P., 2020, A&A, 641, A4
Davis (1951) Davis Leverett J., 1951, Phys. Rev., 81, 890
De Vis et al. (2017) De Vis P., et al., 2017, MNRAS, 471, 1743
Draine (2003) Draine B. T., 2003, ARA&A, 41, 241
Draine & Hensley (2012) Draine B. T., Hensley B., 2012, ApJ, 757, 103
Esquivel & Lazarian (2005) Esquivel A., Lazarian A., 2005, ApJ, 631, 320
Fukui et al. (2008) Fukui Y., et al., 2008, ApJS, 178, 56
Gaensler et al. (2005) Gaensler B. M., Haverkorn M., Staveley-Smith L., Dickey J. M., McClure-Griffiths N. M., Dickel J. R., Wolleben M., 2005, Science, 307, 1610
Gaia Collaboration (2018) Gaia Collaboration 2018, VizieR Online Data Catalog, p. I/345
Gaia Collaboration et al. (2016) Gaia Collaboration et al., 2016, A&A, 595, A1
Goldreich & Sridhar (1995a) Goldreich P., Sridhar S., 1995a, ApJ, 438, 763
Goldreich & Sridhar (1995b) Goldreich P., Sridhar S., 1995b, ApJ, 438, 763
González-Casanova & Lazarian (2017) González-Casanova D. F., Lazarian A., 2017, ApJ, 835, 41
Gordon et al. (2014) Gordon K. D., et al., 2014, ApJ, 797, 85
Grady et al. (2021) Grady J., Belokurov V., Evans N. W., 2021, ApJ, 909, 150
Grand et al. (2017) Grand R. J. J., et al., 2017, MNRAS, 467, 179
HI4PI Collaboration et al. (2016) HI4PI Collaboration et al., 2016, A&A, 594, A116
Haynes et al. (1990) Haynes R. F., et al., 1990, Proceedings of the Astronomical Society of Australia, 8, 339
Hildebrand et al. (2009) Hildebrand R., Kirby L., Dotson J., Houde M., Vaillancourt J., 2009, The Astrophysical Journal, 696, 567
Hu et al. (2019) Hu Y., et al., 2019, Nature Astronomy, 3, 776
Hu et al. (2020) Hu Y., Lazarian A., Yuen K. H., 2020, arXiv e-prints, p. arXiv:2002.06754
Hu et al. (2021) Hu Y., Lazarian A., Stanimirovic S., 2021, arXiv e-prints, p. arXiv:2102.06225
Kalberla & Haud (2015) Kalberla P. M. W., Haud U., 2015, A&A, 578, A78
Kalberla et al. (2021) Kalberla P. M. W., Kerp J., Haud U., 2021, A&A, submitted
Kim et al. (2016) Kim J., Jeong W.-S., Pak S., Park W.-K., Tamura M., 2016, ApJS, 222, 2
Koch et al. (2012) Koch P. M., Tang Y.-W., Ho P. T. P., 2012, The Astrophysical Journal, 747, 79
Lazarian & Hoang (2007) Lazarian A., Hoang T., 2007, The Astrophysical Journal, 669, L77
Lazarian & Pogosyan (2000) Lazarian A., Pogosyan D., 2000, ApJ, 537, 720
Lazarian & Vishniac (1999) Lazarian A., Vishniac E. T., 1999, ApJ, 517, 700
Lazarian & Yuen (2018) Lazarian A., Yuen K. H., 2018, ApJ, 853, 96
Liu et al. (2022) Liu M., Hu Y., Lazarian A., 2022, MNRAS, 510, 4952
Lu et al. (2020) Lu Z., Lazarian A., Pogosyan D., 2020, MNRAS, 496, 2868
Mao et al. (2012) Mao S. A., et al., 2012, ApJ, 759, 25
Mathewson & Ford (1970) Mathewson D. S., Ford V. L., 1970, ApJ, 160, L43
Matsuura et al. (2009) Matsuura M., et al., 2009, MNRAS, 396, 918
McClure-Griffiths et al. (2009) McClure-Griffiths N. M., et al., 2009, ApJS, 181, 398
Montier et al. (2015a) Montier L., Plaszczynski S., Levrier F., Tristram M., Alina D., Ristorcelli I., Bernard J.-P., 2015a, A&A, 574, A135
Montier et al. (2015b) Montier L., Plaszczynski S., Levrier F., Tristram M., Ristorcelli I., Alina D., Bernard J.-P., Guillet V., 2015b, A&A, 574, A136
Olszewski et al. (1991) Olszewski E. W., Schommer R. A., Suntzeff N. B., Harris H. C., 1991, AJ, 101, 515
Pattle (2019) Pattle K., 2019, Nature Astronomy, 3, 692
Peek et al. (2018) Peek J. E. G., et al., 2018, ApJS, 234, 2
Planck Collaboration (2018) Planck Collaboration 2018, A&A, submitted
Planck Collaboration I (2011) Planck Collaboration I 2011, A&A, 536, A1
Planck Collaboration Int. XIX (2015) Planck Collaboration Int. XIX 2015, A&A, 576, A104
Planck Collaboration Int. XXXV (2016) Planck Collaboration Int. XXXV 2016, A&A, 586, A138
Planck Collaboration XI (2014) Planck Collaboration XI 2014, A&A, 571, A11
Planck Collaboration et al. (2011) Planck Collaboration et al., 2011, A&A, 536, A17
Quinn (2012) Quinn J., 2012, A&A, 538, A65
Remazeilles et al. (2011) Remazeilles M., Delabrouille J., Cardoso J.-F., 2011, MNRAS, 418, 467
Rémy-Ruyer et al. (2013) Rémy-Ruyer A., et al., 2013, A&A, 557, A95
Schmidt (1970) Schmidt T., 1970, A&A, 6, 294
Schmidt (1976) Schmidt T., 1976, A&AS, 24, 357
Serkowski (1958) Serkowski K., 1958, Acta Astronomica, 8, 135
Simmons & Stewart (1985) Simmons J., Stewart B., 1985, A&A, 142, 100
Staveley-Smith et al. (2003) Staveley-Smith L., Kim S., Calabretta M. R., Haynes R. F., Kesteven M. J., 2003, MNRAS, 339, 87
Vaillancourt (2006) Vaillancourt J. E., 2006, Publications of the Astronomical Society of the Pacific, 118, 1340
Vandenbroucke et al. (2021) Vandenbroucke B., Baes M., Camps P., Kapoor A. U., Barrientos D., Bernard J. P., 2021, A&A, 653, A34
Visvanathan (1966) Visvanathan N., 1966, MNRAS, 132, 423
Wardle & Kronberg (1974) Wardle J., Kronberg P., 1974, ApJ, 194, 249
Wenger et al. (2000) Wenger M., et al., 2000, A&AS, 143, 9
Wielebinski (1993) Wielebinski R., 1993, in Krause F., Radler K. H., Rudiger G., eds, Kluwer Academic Publishers Vol. 157, The Cosmic Dynamo. p. 317
Wisniewski et al. (2007) Wisniewski J. P., Bjorkman K. S., Magalhães A. M., Pereyra A., 2007, ApJ, 664, 296
Xu et al. (2019) Xu S., Ji S., Lazarian A., 2019, ApJ, 878, 157
Yuen & Lazarian (2017) Yuen K. H., Lazarian A., 2017, ApJ, 837, L24
Yuen et al. (2021) Yuen K. H., Ho K. W., Lazarian A., 2021, ApJ, 910, 161

Appendix A Additional figures

Figure 21 shows, from left to right, the maps of the signal-to-noise ratios of intensity ( $I/\sigma_{II}$ ) and polarization fraction ( $p/\sigma_{p}$ ) and the uncertainty of the polarization angle $\sigma_{\psi}$ . There, $\rm\sigma_{II}$ is given from the original smoothed Planck data while $\sigma_{p}$ and $\sigma_{\psi}$ are obtained using Bayesian analysis. Figure 22 shows the comparison between magnetic field of the LMC in the original Planck data and after foreground removal, represented by black and red segments respectively.

Appendix B Estimation of the noise covariance parameters after foreground removal

Below we describe the method we used to estimate the covariances between the Stokes I, Q, and U parameters after the foreground removal. First, we estimate the covariances in the region around the LMC, outside of the red contour in Figure 3: $C^{\prime}_{II},C^{\prime}_{QQ},C^{\prime}_{UU},C^{\prime}_{IQ},C_{IU},C_{QU}$ , where $C$ denotes the covariance, classically noted as $\sigma^{2}$ , to simplify the writing. Second, we estimate the average intensity in the same region: $I^{\prime}$ . Third, we estimate the variances of the LMC foreground using the following relationship:

C^{fg}_{XY}(i,j)=C^{\prime}_{XY}\frac{I_{fg}(i,j)}{I^{\prime}}\,,

(13)

where $X,Y$ are some Stokes parameters (I, Q or U), $(i,j)$ is a 2D position of a pixel, and $I_{fg}(i,j)$ is the foreground intensity at each pixel estimated in Sect. 3.2. Then, we add them to the initial covariances to obtain the covariances of the parameters after the foreground removal:

C^{LMC}_{XY}(i,j)=C^{fg}_{XY}(i,j)+t\sigma^{2}_{XY}(i,j)\,,

(14)

where $t$ is the factor that takes into account the dfference in the pixels size in HEALPix and WCS formats. Finally, the variance on the obtained polarization fraction of the LMC is computed using the classical approach (Wardle & Kronberg, 1974) extended to include non-diagonal terms of the covariance matrix (Montier et al., 2015b):

	$\displaystyle\sigma^{2}_{p}=\dfrac{1}{p^{2}I^{4}}\big{(}Q^{2}\sigma_{Q}^{2}$	$\displaystyle+$	$\displaystyle U^{2}\sigma_{U}^{2}+p^{4}I^{2}\sigma_{I}^{2}+2QU\sigma_{QU}$		(15)
		$\displaystyle-$	$\displaystyle 2QIp^{2}\sigma_{QI}-2UIp^{2}\sigma_{UI}\big{)}\,.$		(15)