© ESO 2018

1 Introduction

Magnetic fields are dynamically important constituents of galaxies, playing a major role, for example, in the star formation process and controlling the density and propagation of cosmic rays (see e.g. Beck 2016, and references therein). Observing magnetic fields, however, is non-trivial, as indirect observations are required, based primarily on dust polarization, Zeeman effect, synchrotron radiation, its polarization, and Faraday rotation of the polarization angle, which is referred to as rotation measure (hereafter RM; see e.g. Klein & Fletcher 2015, and references therein). Because all such methods have strong limitations, interpretation of the data is difficult, especially for the Milky Way, inside of which we reside. This is where radiative transfer simulations combined with numerical models may become useful, bridging the differences between physical models and indirect observations (e.g. Ostriker et al. 2001; Falceta-Gonçalves et al. 2008; Pelkonen et al. 2009; Planck Collaboration Int. XX 2015, hereafter PlanckXX).

Planck is a space mission that mapped the anisotropies of the cosmic microwave background (CMB; see e.g. Planck Collaboration I 2011, 2016). This mission also has the capability to measure thermal emission and its polarization from dust grains in all bands up to 353 GHz. Particularly, the foreground dust can be studied with the High Frequency Instrument in the frequency range 100–857 GHz (HFI; see e.g. Lamarre et al. 2010). Polarized dust emission and its spatial variations have been mapped with high resolution and sensitivity in a series of papers. Planck Collaboration Int. XIX (2015, PlanckXIX hereafter) studied the all-sky dust emission at 353 GHz, where polarized emission is most significant, and PlanckXX computed the statistics of polarization fractions and angles outside the Galactic plane. In Planck Collaboration Int. XXI (2015) thermal dust emission was compared with optical starlight polarization. Planck Collaboration Int. XXII (2015) presented a study of the variation of dust emission as a function of frequency inthe range 70–353 GHz. The results relating to the polarized thermal dust emissions are summarized in Planck Collaboration I (2016, Sect. 11.2). For this paper, the most relevant study in this series is the all-sky study of PlanckXIX (and subsequent updates reported in Planck Collaboration X 2016), as our modelling efforts concentrate on kiloparsec (kpc) scale magnetohydrodynamic (MHD) models including all the three phases (cold, warm, and hot) of the interstellar medium (ISM), regulated by supernova (SN) activity and subject to large- and small-scale dynamo instabilities. Here cold would correspond to cold neutral medium (CNM) with typical temperatures of 100 K and number densities of ~100 cm⁻³, warm to 10⁴ K and 0.1 cm⁻³, and hot to 10⁶ K and 0.001 cm⁻³. The major findings of PlanckXIX include the discovery of anti-correlation between the polarization fraction, p, and polarization angle dispersion, , and the decrease in the maximum polarization fraction, p_max, as column density increases. These major findings generally hold in zoomed-in regions near molecular cloud complexes (PlanckXX).

Efforts to compare polarization measurements with synthetic observations from MHD models have concentrated on the scale of molecularclouds. For such comparisons, the relevant MHD models normally include cold and warm phases of the ISM, describe the magnetic field as a uniform background field, and may include artificial flows to enhance the formation of shocks (e.g. Ostriker et al. 2001; Padoan et al. 2001; Bethell et al. 2007; Falceta-Gonçalves et al. 2008; Hennebelle et al. 2008; Soler et al. 2013; Planck Collaboration Int. XX 2015; Chen et al. 2016). Appropriately the Planck data provide high enoughsensitivity and resolution for studies at this length scale. In general, a satisfactory agreement between the synthetic and observed dust polarization properties has been found, and the anisotropic and turbulent character of the magnetic field has been identified as the most decisive factor, particularly PlanckXX demonstrated the connection of p and with turbulent magnetic fluctuations.

The large-scale dynamics of the ISM in the star-forming parts of spiral galaxies can be described with a three-phase medium regulated by stellar energy input (McKee & Ostriker 1977). By far the dominating energy source to power turbulence at the 100 pc scale (see Abbott 1982) originates from SNe. These SNe bring significant input of thermal and kinetic energy to the ISM. In the solar neighbourhood Tammann et al. (1994) estimated that SNe inject approximately 3 × 10⁵² erg kpc⁻² Myr⁻¹ thermal energy, which is dissipated mainly as heat into the ISM, but with some 10% converted into kinetic energy (Chevalier 1977; Lozinskaya 1992). In addition to driving expanding shock fronts that interact with each other, the SNe generate bubbles of hot gas near the galactic disk; these SNe, along with the cold ISM and molecular clouds, are embedded within the diffuse warm component. The filling factor of the hot gas is small near the galactic midplane but approaches unity in the halo (Ferrière 2001).

In addition to the re-structuring and mixing of the ISM, SN-forcing powers the galactic dynamo in the rotating anisotropic galactic disk. Anisotropic turbulence and a non-uniform rotation profile combine to provide the ingredients for large-scale dynamo instability, leading to the generation and maintenance of magnetic fields that are dynamically significant on a galactic scale. Along with the mean magnetic field, a strong fluctuating field is also generated. The dynamo processes are intrinsically connected to the three-phase structure of the ISM, such that both the large- and small-scale filling factor and topology are different in various phases and locations of the galactic disk. Recent numerical MHD models have attained sufficient realism to model these processes self-consistently (Gent et al. 2013b; Hennebelle & Iffrig 2014; Kim & Ostriker 2015; Bendre et al. 2015; Evirgen et al. 2017; Hollins et al. 2017). These developments enable us to study the influence of the three-phase medium, SN shock fronts, and dynamo-generated magnetic fields on the observable properties of dust polarization at large scales.

In this work we study the influence of all the afore mentioned physical effects on the polarization properties of the galactic ISM. Apart from PlanckXX and Planck Collaboration Int. XXXV (2016), our approach is different from most of the previous modelling efforts related to Planck observations in that we are not building models to explain the observations specifically. However, we aim to synthesize a set of independently built turbulence models to test the relevance of the physical effects they contain to explain real observations. Our strategy is also designed to be contrasted with studies that build statistical or phenomenological models to the observed all-sky polarization properties (e.g. Planck Collaboration Int. XXXII 2016; Planck Collaboration Int. XLII 2016; Planck Collaboration Int. XLIV 2016).

To be able to resolve, on the one hand, the SN-generated turbulence and, on the other hand, to allow for self-consistent dynamo action on kiloparsec scales, the simulation set-up is limited by these requirements, and we are therefore not free to arbitrarily choose the resolution optimally suited to the observations. In this study, in particular, we are limited to 4 pc grid resolution with the kiloparsec box size chosen.

The chosen resolution in practice means that we can model the cold cloud components only up to certain densities, therefore preventing these clouds from becoming gravitationally bound. The main regulatory mechanism creating the cold and warm phases is, however, properly included, namely the thermally unstable cooling function. Therefore, we consider that we provide a realistic description of the three-phase medium. At the box scale, the dimensions are yet too small for the large-scale effects such as spiral arms or central bulge to be realistically included. We note that in contrast Planck Collaboration Int. XLII (2016), who have built a phenomenological whole-sky model of the Milky Way, do model the spiral arms, but they do not include a physically self-consistent model of turbulence. Therefore, we restrict ourselves to consider only large-scale effects arising from rotation and its non-uniformities (differential rotation), both of which are needed to enable and sustain dynamo action in the system.

The paper is organized as follows. In Sect. 2 we describe the tools and methods used in this study. In Sect. 3 we present the simulated polarization results and compare these results to the observations of PlanckXIX. In Sect. 4 we consider how the shock and magnetic field affect the interpretation of observations. In Sect. 5, concluding the paper, we discuss the implications of our results and potential for further studies.

2 Methods

2.1 Numerical MHD simulations

In PlanckXX, polarization statistics were compared to MHD simulations, which include cold and warm phases of the ISM. These employ adaptive mesh refinement in a computational cube 50 pc across (Hennebelle et al. 2008) from which an 18³ pc³ subset is selected for analysis. We added comparisons to MHD simulations of the ISM, in which the turbulence is driven by SNe (Gent et al. 2013b; Gent 2012, Chapters 8 and 9). In this model, the cold and warm phases are produced, as in the two-phase models, through regulation by thermally unstable radiative cooling, but with the addition of a hot phase generated by SN heating.

To capture all the relevant dynamics of the three-phase model, the simulation domain size has to be increased. The grid is 256 × 256 × 560 and spans horizontally 1.024 kpc and vertically ± 1.12 kpc about the galactic midplane. The SN supersonic forcing naturally generates a highly shocked turbulent flow, so no artificial forcing is applied. Moreover, the interaction of rotation and anisotropic turbulence with the galactic shear flow induces a natural magnetic field through dynamo action. To model dynamo action, we solved non-ideal MHD, including viscous, thermal, and magnetic diffusivities. With temperatures spanning 7–8 orders of magnitude, and high Mach numbers, it is not possible to apply physically motivated values for diffusivity. To resolve the flows in the hot gas, while obtaining optimal small-scale structure to the turbulence in the cold and warm ISM, we set the viscosity proportional to the sound speed (or ), which may be significant for analysis in Sect. 3.1. This viscosity scheme is not chosen to model the actual turbulent diffusivities, as little is known about their dependence on the key physical quantities in the ISM. The analytical estimates obtained in the first-order smoothing approximation framework (Steenbeck et al. 1966), however, imply that the magnitude of turbulent diffusion is orders of magnitude larger than the Spitzer molecular values. The motivation for the used viscosity scheme is to resolve much finer structures in the cool and warm ISM than would be possible with constant diffusivity.

The resolution of 4 pc along each side in these simulations without magnetic fields (Gent 2012, Table 5.1) yields a maximum gas number density for the ISM of about 100 cm⁻³ and the fractional volume of the cold gas is 0.4%. The fractional volume of the warm gas is 60%; the hot gas fractional volume is about 28% near the SN active midplane and increases to 41.5% elsewhere. With the magnetic field included (Evirgen et al. in prep.) the model fractional volume of the warm gas increases to about 80% because the hot gas is pushed further into the halo, and the cold gas is confined to an even smaller fractional volume. A snapshot of the thermodynamic profile of the model ISM is shown in Fig. 1 (left), from which the distribution of the different gas phases is evident. Observed densities and those arising from the MHD simulations of Hennebelle et al. (2008) extend to much higher densities and increased fractional volume of cold ISM. The characteristic properties of the MHD simulation data are listed in Table 1. Temperatures, velocities, and magnetic field strengths are better representatives of the observed ISM, but the smallest scales of their fluctuations are limited by the grid resolution. The saturation of the magnetic field has the effect of restricting the flow and increasing the homogeneity of the ISM, such that the maximum density reduces to about 10 cm⁻³, which must be taken into account when making comparisons with the Planck observations and the earlier MHD molecular cloud model. The Mach number in the simulations reaches as high as 25.

Both large- and small-scale dynamo instabilities are present in the system. A 3D rendition of the magnetic fieldlines embedded in this atmosphere is illustrated in Fig. 1 (right). As reported by Evirgen et al. (2017), the strength of the generated mean magnetic field at the midplane in the MHD model is in good agreement with the global observational estimates of 1.6 μG by Rand & Kulkarni (1989), using pulsar RM. The random field in the model, however, is much weaker than observed, as the random field is only 20–50% of the 5 μG estimate of Rand & Kulkarni (1989) or 5.5 μG of Haverkorn (2015), combining RM with thermal electron density measures. These estimates are supported by the reviews of galactic magnetic fields (Beck et al. 1996; Beck 2016). Some preliminary results applying synthetic RM measurements to the MHD modelconsidered here are reported in Hollins et al. (2017), and further such application shall be considered in future work. The model and characteristics of the multiphase structure of the simulated ISM are described in detail in Gent et al. (2013a), and a summary is included in Appendix B. The coherence and fluctuations of the magnetic field are important to the polarization measurements, so it is useful to decompose the field into the mean field and fluctuations b, where (1)

In this model, the entire field is subject to stirring due to the action of the many SN remnants, so the convention of separation into large-scale and turbulent components needs more careful consideration. We separate B by volume averaging with a Gaussian kernel (at l = 50 pc scale; see Gent et al. 2013b). In Fig. 2 (left) fieldlines and (right) b fieldlines are plotted over density isosurfaces on background slices of temperature. With this treatment the mean field, which is what we shall refer to as the large scale has strength and orientation that varies with space and time, but is coherent, i.e. similarly oriented, over scales above about 100 pc. The small-scale field is derived by subtracting the large-scale field from the total field for any given snapshot. Hence the large-scale field still exhibits spatial structure, which would influence polarization fractions (see Sect. 2.3). The random field is clearly incoherent and would contribute to depolarization. Some indicative values are listed in Table 1 for the ranges across the MHD snapshots of gas number density, temperature, magnetic field, and speed.

The gas density ρ from each snapshot is used with the radiative transfer code SOC, as described in Sect. 2.3, to model the dust temperature distribution, and the magnetic field B is then employed when simulating the dust polarization observations. From the velocity field we can compute directly its divergence, ∇ × u, which is also used in the analysis of the polarization to determine the influence of the shock structure of the ISM on the observations.

The simulation used for this analysis applies parameters for gas density, stellar, and dark halo gravitation, galactic rotation, and SN rates and distribution matching estimates for the solar neighbourhood. The magnetic field is amplified for a period exceeding a Gyr by dynamo action from a seed field of a few nG, which then saturates with an average field strength of a few μG. The strength of the mean field portion is consistent with observational estimates, while the random field strength is 2–5 times weaker than estimated. This has some influence on the dispersion of polarization angles in the simulated observations, discussed in Sect. 4. We used 12 snapshots from the saturated statistical steady stage of the model, each separated by 25 Myr, commencing at 1.4 Gyr. For simulated maps presented in this paper, we used the snapshot at 1.7 Gyr. The system is in a statistical steady state, so individual snapshot characteristics are representative, but strong temporal and spatial differences are also present. For some of the analysis we consider ensemble averages of the snapshots to identify the most persistent structures.

2.2 Stellar radiation field

Dust in the ISM is illuminated by the stellar radiation field. We inverted the measurements of the average stellar radiation field in the solar neighbourhood from Mathis et al. (1983) into a distribution of radiation sources to model the stellar radiation for our radiative transfer model. We then modelled this radiation with a horizontal density profile (see Gent et al. 2013a, excluding the dark matter component), and emissivity j_ν reflecting the vertical distribution of stars as (2)

Here a₁ = 4.4 × 10⁻¹⁴ km s ⁻², z_a = 200.0 pc and j_ν is the emissivity over the spectrum of frequencies ν. The normalization coefficients, j_ν,0, were determined to return the expected total intensity J_ν from (3)

where , and τ is the optical depth along the line-of-sight (LOS). The distribution was generated using the Monte Carlo method and was scaled to match J_ν,Mathis from Mathis et al. (1983). Using this inversion, we obtained an approximate z-dependent radiation field, which produces reasonable dust temperatures with our radiative transfer simulations.

Fig. 1 Left: representative slice in the xz-plane of gas number density n [1∕cm−3] from a singlesnapshot overplotted with contours of temperature T [K], illustrates the multiphase ISM. Right: snapshot of the model ISM with temperature slices in the background and isosurfaces of the density in the foreground; fieldlines of the magnetic field are overplotted.

2.3 Radiative transfer calculations

To calculate the dust emission, we used the program SOC, which is a new Monte Carlo code for continuum radiative transfer tested against the CRT program (Juvela 2005) and also against several other codes participating in the Transport of Radiation through a DUSTy Medium (TRUST)¹ benchmark project on 3D continuum radiative transfer codes (Gordon et al. 2017). The density distribution of the models can be defined withregular or modified Cartesian grids or hierarchical octree grids. In this paper, all calculations employ regular Cartesian grids. In addition to the density field, the program requires a description of the dust grains (i.e. absorption and scattering properties) and radiation sources. The program employs a fixed frequency grid to simulate the radiation transport atdiscrete frequencies, one frequency at a time. The information of the absorbed energy is used to solve the grain temperatures for each cell of the model. The dust model could also include small stochastically heated grains. However, to predict dust emission at submillimetre wavelengths, the calculations are limited to large grains that are assumed to remain at a constant temperature in equilibrium with the local radiation field. Once the dust temperatures have been solved, synthetic images of dust emission can be calculated towards selected directions or, as in the case of the present study, over the whole sky as seen by an observer located inside the model volume.

We assumed a constant gas-to-dust ratio, where we followed the dust model BARE-GR-S of Zubko et al. (2004), which has been created to match the observations of typical Milky Way dust with the extinction factor R_v =3.1. For simplicity, we kept the properties of the dust similar throughout the whole computational domain. Dust properties may vary between different ISM phases, but the investigation of these effects is beyond the scope of this initial study.

During the simulation of the internal radiation field, the radiation transport was calculated without taking the polarization into account. However, SOC includes tools to produce synthetic polarization maps. The local grain alignment efficiency could be calculated following the predictions of the radiative torques theory (Draine & Weingartner 1996, Pelkonen et al. 2009, PlanckXX). Because our study concentrates on emission from low-density medium, we omitted this step and weessentially assumed a constant dust grain alignment efficiency. The polarization reduction factor p₀ was simply set to a constant value that results in a maximum polarization fraction that is consistent with observations. Maps were calculated separately for the Stokes I, Q, and U components, taking into account the local total emission, local magnetic field direction, and value of p₀. These data were then finally converted to maps of polarization fraction and polarization angle dispersion. A representative set of synthetic Stokes I, Q, and U maps is presented in Appendix A.

To calculate the polarization within a single cell, we applied the following method. We used the Planck/HEALPix² convention for the polarization angle ψ (Planck Collaboration Int. XIX 2015; Planck Collaboration Int. XXI 2015). Within a cell we normalized the direction of the magnetic field, (4)

and calculated ψ with (5)

where the and are the directional vectors of the HEALPix coordinate directions. The angle between the magnetic field and the plane of the sky (POS) is (6)

where is the direction of the LOS. Based on the non-polarized emitted intensity of the cell, I_orig, we get the Stokes components I, Q, and U with

To match the polarization degrees observed in the PlanckXIX and PlanckXX we set p₀ = 0.2.

From the integrated Q and U we can calculate the polarization fractions p and the polarization angle dispersion functions over the whole sky. Simulated p and yield familiar measures of system properties and allow comparison to previous studies. The polarization fraction is defined as (10)

and the polarization angle dispersion function (Falceta-Gonçalves et al. 2008; Hildebrand et al. 2009; Planck Collaboration Int. XIX 2015) as (11)

where ψ(r) is the polarization angle in the given position in the sky r and ψ(r + δ_i) the polarization angle at a position displaced from the centre by the vector δ_i. The sum extends over pixels whose distances from the central pixel in the location are between δ∕2 and 3δ∕2. As in PlanckXIX, we calculated the dispersion function with (12)

In addition, to be comparable with the analysis presented in PlanckXIX we applied 1° Gaussian smoothing to our simulated observations and set δ = 30′ for all calculations of . This sets the widest distance between points in the annulus to be 1.5 times the size of the FWHM of the Gaussian filter. The size of the annulus is small because, according to PlanckXIX the angle dispersion function gradually loses its coherence with increasing δ. PlanckXIX also noted that their measurements of are not an artefact of either the choice of δ itself or instrumental bias. However, the small annulus includes some influence due to spatial correlations induced by the 1° beam. The choice of δ is independent of our radiative transfer model, and its primary justification is to compare our results with PlanckXIX.

We calculated our results with several integration distances R_max, namely R_max = 0.25, 0.5, 1, 2, 4 kpc, to explore the effect of depth. A part of the emission observed in PlanckXIX is received from distances larger than what our model is capable of exploring. However, with this in mind, our motivation is to inform, with an independently and physically generated ISM and magnetic field, how the impact of near and distant features might contribute to the observed images. We assumed periodic boundary conditions in the horizontal direction of the computational domain. However, if the boundary is crossed in the vertical direction, the integration stops at that point, as happens in the cases R_max = 2 and 4 kpc. For our analysis, we masked latitudes |b| > 45° from our maps. This more closely reflects the region of the sky where the Planck observations are included in the analysis of PlanckXIX. This is an arbitrary choice, as we are not limited by the signal-to-noise issues of PlanckXIX. However, this has also other benefits. It avoids the problem of calculating near the poles and some of the issues of asymmetry caused by the limited vertical extent of the domain, which would affect the comparison of R_max ≤ 1 kpc and R_max ≥ 2 kpc. To save space, we omitted from the figures related to the case with R_max = 0.5 kpc, but the results are similar to R_max = 1 kpc.

Fig. 2 Snapshot of the model ISM, with temperature slices in the background and isosurfaces of the density in the foreground, as in Fig 1, but with the total magnetic B fieldlines replaced by fieldlines (left) and b fieldlines (right). Further visualization of the MHD, including video representations, is available online. (the full set is provided at http://fagent.wikidot.com/astro).

2.4 Combining data into 2D histograms

Many of theresults in this paper are represented using joint 2D histograms. This allows us both to compare our results directly to the analysis in PlanckXIX and PlanckXX, and separate the snapshot and point-of-view (POV) specific details from the general large-scale behaviour.

For 2D histograms, we combined radiative transfer simulations from all 12 snapshots and for 5 different POV positions in the approximate midplane, specifically from the points (0.5, 0.5, 1.0) kpc, (0.25, 0.5, 1.0) kpc, (0.5, 0.25, 1.0) kpc, (0.75, 0.5, 1.0) kpc, and (0.5, 0.75, 1.0) kpc, under the assumption that the use of periodic boundary conditions in x- and y-directions should give a reasonable representation of the ISM in the neighbourhood of the solar system.

The data were combined as follows. From each SOC-generated map, i = 1–60, separate 2D histograms were constructed over a common set of 512 × 512 bins for eachcomparison, specifically C_bin,i(N_H, p), , . Here C_bin,i represents the number of counts in a single bin calculated from a single synthetic observation. Then, assuming C_bin,i to be lognormally distributed between separate 2D histograms, we calculated the mean and standard deviation over i, (13)

We included only the bins C_bin for which C_bin,i > 0 for all i. Second, to reduce additional noise, we masked out the elements where (14)

yielding the composite 2D histograms, as presented in Fig. 3 and subsequent 2D-histogram figures. Combined 2D histograms, summing counts over many snapshots and POV, contain the patterns most consistent across the individual 2D histograms. Despite masking low count noise, many low-frequency outliers persist.

The various simulated observations may each exhibit distinct features worthy of further analysis in the future. However, in this study we focus on the general behaviour of the system. For example, the 2D histograms only include pixel-by-pixel numerical values. Therefore, such analysis cannot regard the filament-like shapes visible in -maps (see Figs. 5, 9 and 14).

3 Properties of synthetic polarization observations

We present our results and analysis based on the synthetic observations obtained using SOC with polarization tools on top of our MHD simulations. For our analysis, we use the maps of column densities N_H, polarization fractions p, and polarization angle dispersions . As discussed in Sect. 2.1, we take (in the figures) as the reference case the 1.7 Gyr snapshot, where the observer is situated in the centre of the computational domain. We opt for projection onto a sphere as the most appropriate observational reference frame. With the domain spanning over 1 kpc in each direction, it is too large to be treated as a single observation along a Cartesian coordinate. Thus, projection onto a sphere is the most realistic observational frame, and is especially so when nearby regions are concerned. The all-sky point of view affords the best comparison with Planck on the scales relevant to the MHD model. All of the maps of synthetic observations over the whole sky are presented using Mollweide projection.

3.1 Simulated column density

The maps of column density resulting from integrating over each distance, up to R_max = 0.25, 2 and 4 kpc from the observer, are shown in Fig. 3, right panels. The dependence of N_H on R_max is immediatelyapparent. The densest regions are situated near the midplane, as would be expected from horizontal averaging of the gas density in the MHD model represented in Fig. 1. However, this is not visible on the map from R_max = 0.25 kpc, but clear for the higher R_max. The vertical anisotropy in the latter reflects the temporal upward shift of the disk centre of mass, which is evident from Fig. 1. Also, density should be near isotropic in the latitudinal direction because we do not model either the central bulge or spiral arms of the Galaxy. In a single snapshot, however, local SN bubbles or superbubbles (merging multiple SN remnants) may inject significant anisotropy into the overall density profile. This is clearly seen only for the R_max < 1 kpc, where bubbles are apparent on several locations. For the higher lengths of integration the influence of features near to the observer is almost negligible. As the combination of nearby and far-away emission is present in both real and synthetic observations, exploring the effect of distance to the observations is called for.

Figure 3, left panels, show joint 2D histograms of polarization angle dispersion and column density N_H for R_max  =  0.25, 2 and 4 kpc. At 0.25 kpc, N_H is clustered around only 1.5 × 10²⁰ cm⁻² and , while its corresponding N_H map, dominated by local gas structure, is without dense clouds. This is inconsistent with PlanckXIX and indicates that the integration length is too short. For higher R_max, N_H is typically above 10²¹ cm⁻² and is consistent with the observed column density featured in PlanckXIX, Fig. 24 top panel. The value also increases above 1, but not to the level obtained by PlanckXIX. The high R_max maps of N_H also capture the presence of these clouds, particularly evident near the midplane.

On these joint 2D histograms are also plotted lines of maximum C_bin (green) and weighted mean (blue) polarization fractions as functions of N_H. For R_max ≳ 1 kpc the trend has a positive correlation with N_H, in disagreement with PlanckXIX Fig. 24. There is no record in PlanckXX on how this relationship applies to their MHD model. Comparing our synthetic observations with PlanckXIX Fig. 24, there is a peak in the weighted mean profile for high column density dust. This corresponds to the molecular clouds, which are not resolved in the MHD model. In Fig. 3 (R_max > 1 kpc) the high columndensity range, representing the warm and cold ISM, is a good match for the PlanckXIX data. In the PlanckXIX data, however, the weighted mean remains as high or even increases at lower densities, while our simulations show a correlation in this range of column densities, which are in obvious disagreement.

Let us try to understand this discrepancy by assuming that the small-scale scatter seen through is a result of the underlying turbulent diffusion, taking the action of smoothing the flow. If this was the case, one would expect the higher viscosity regions to smooth out structures, along with those seen in , more efficiently; that is a relation would be expected,where ν_t is a turbulent diffusion coefficient related to turbulent mixing, but not necessarily to any dependence assumed in viscosity scheme used in the model. As explained in Sect. 2.1, gas viscosity is set as ν ∝ T^0.5 so as to resolve flows in the hot gas of the MHD model. The ISM is modelled as an ideal gas, in which all phases are in approximate statistical pressure balance, so this is statistically equivalent to ν ∝ ρ^−0.5. For the radiative transfer calculations, the dust density is assumed proportional to the ISM gas density. If higher viscosity in the hot gas tends more proportionally to smooth small-scale structures in the flow, we might expect the dispersion , which we indeed see for the low-density hot gas. This implies that the diffusion in the hot gas reflects the dependence input through the viscosity scheme, while that for the cold and warm gas does not do so, but is somewhat steeper than expected from this simple argumentation.

We have no reliable method to determine observationally the true turbulent viscosity in the ISM. Some elaborate and approximate methods make this possible within numerical experiments (see e.g. Käpylä et al. 2018). Our simple hypothesis presented above could be tested by applying the same analysis in this paper to MHD simulations with alternative prescriptions for viscosity to exclude a relationship between the weighted mean profile and the viscosity, as appears to exist here. However, if MHD models consistently display such a dependence on viscosity for the weighted mean , then we might be able to infer something about the actual turbulent viscosity in the ISM from the PlanckXIX profiles. Using PlanckXIX with the weak negative correlation in between the weighted mean values of and N_H would imply ν ∝ T^−λ for the ionized ISM, where 0 < λ ≪ 1. This is unlike, for example Spitzer molecular viscosity ν ∝ ρ⁻¹T^2.5. The scale of such simulations place that investigation beyond the scope of the present study, and of course there are many complicating factors that would permit alternative explanations of this trend.

At the lowest R_max the visible cloud structure in the N_H map traces our MHD model. With increasing R_max smaller details appear in the maps, which are the effects of projection from more remote sources. Hence, the observed results are inevitably a combination of both the MHD model volume and its projection through a mixture of features along the LOS. By assuming periodicity in the horizontal direction, we can reach higher R_max than the physical scale of our MHD model. However, with too high R_max the periodicity creates a mirror gallery effect of repeating patterns near the x- and y-axes of the MHD model. To avoid this, it is reasonable to restrict the highest integration distance to R_max = 2 kpc, or for the general case, R_max ≲ 2L_x, where L_x is the horizontal extent of the MHD model.

Fig. 3 Left: joint 2D histograms of polarization angle dispersion and column density NH. The green lines follow the maximum Cbin as a functionof NH, and the blue lines represent the weighted mean of . Right: sample maps of column density. For the case Rmax = 1 kpc, see Fig. 13.

3.2 Polarization fraction across the galactic plane

For each R_max we construct joint 2D histograms of polarization fraction p and N_H, presented in the left panels of Fig. 4, and maps of p in the right panels. With R_max = 0.25 kpc there is a broad range for p values. The distribution is independent of the N_H values, which in this case are much lower than in the PlanckXIX observations. The mapping of p can be seen to correlate over large smooth regions. There is negligible depolarization and the highest C_bin appears on the scale of p₀ = 20%.

With R_max = 2 kpc, the depolarization is stronger, particularly near the midplane; the map of p is an excellent analogue for PlanckXIX, Fig. 6. Its 2D histogram of p and N_H is also a better match with PlanckXIX, Fig. 19, although N_H > 10²² cm⁻² are absent and there are few values with p > 15%. The high column density may relate to features not modelled with the MHD, such as the central bulge, spiral arms, and molecular cloud properties requiring higher resolution to be resolved. Some of the absence of higher p is due to the masking of bins that lack signal across all synthetic maps combined with the 2D histogram in question. In addition, the fraction of the PlanckXIX data, which have p > 15%, is concentrated in the molecular cloud complexes in which high polarization fractions can be found, which are among the features not resolved in our MHD models.

The increased frequency of p > 5% with R_max = 4 kpc is more likethe PlanckXIX results, but the high number of points where p < 5% and N_H ≃ 5 × 10²¹ cm⁻² is not consistent with PlanckXIX; this is an indication that oversampling the same periodic domain near the midplane is distorting the distribution. To improve the range of column densities in a physically consistent manner, one should increase the horizontal extent of the MHD models. Also, increasing the MHD model resolution would improve both the proportion of high column densities and high polarization fractions, but this would be computationally demanding in the context of SN-driven MHD turbulence that is capable of producing self-consistent small- and large-scale dynamos.

Regardless of these current limitations, it is interesting to consider how the distances influence polarization and depolarization. In the framework of the radiative transfer calculations, the mechanism is easy to understand. Along the LOS, individual cells in the grid produce positive and negative contributions in Q and U, depending on the local magnetic field direction. For an incoherent magnetic field, the sign of Q and U fluctuates, with similar magnitude, such that over very long distances in an optically thin medium their averages approach zero. This is comparable to the analysis of Houde et al. (2009), who noted that multiple independent turbulent cells along the LOS make the apparent fluctuation approach zero. However, the mean field also varies in direction when exploring kiloparsec scales, which also contributes to the perceived effect. In the presence of a directionally coherent magnetic field, its direction dominates the Q and U over long distances. However, the magnetic fields are highly turbulent throughout the MHD domain, so we expect depolarization to increase with R_max.

Near the midplane the depolarization effect is strongly proportional to integration length with R_max ≥ 2 kpc (Fig. 4, bottom right panels), corresponding high N_H. This is also consistent with the above explanation. Near the midplane with large R_max it is possible to integrate polarization over a longer distance, unlimited by the upper or lower boundaries of the computational domain. Therefore, there is more influence on p from the depolarizing effect of the fluctuations that are most evident near the midplane. The role of magnetic field fluctuations inducing a broad spread of p has also been suggested by PlanckXIX and PlanckXX and our results support this idea.

Because of the averaging and masking method presented in Sect. 2.4, some of the outlying values persist as halos in the 2D histograms of Fig. 4 (also Fig. 13), but with higher sampling rates the masked points between the bulk data and halo could be restored. Also, the gradual loss of alignment of the dust grains within radiation shielded dense clouds can decrease polarization, which is also considered by PlanckXIX to be an effect influencing their observed depolarization with N_H ≥ 10²² cm⁻². We set the strength of alignment proportional to p₀ in Eqs. (7)–(9), neglecting any effects of such shielding.

As our methods resemble those of PlanckXX, some discussion is called for, although we cannot make a direct comparison of their MHD model with our results. PlanckXX have considered scales well below our 4 pc grid resolution. What is relevant is that the structure of the flow and the magnetic field in our MHD model is naturally driven by the forces on the scales of SN remnants cascading to the lowest turbulent eddies we can resolve. It is reasonable to expect that the turbulent structure would extend to lower scales, until new physical processes become active.

PlanckXX compared their MHD results with observations of the Chamaeleon-Musca and Ophiuchus fields. For their MHD simulation, they viewed the domain with respect to the mean field as POS, LOS, and 50/50, confirming that a high polarization fraction is indicative of a strong POS component to the field. This is evidence that the magnetic field has a strong POS component in Chamaeleon-Musca, while in Ophiuchus the field is morealigned along the LOS; this interpretation is consistent with Planck Collaboration Int. XXXV (2016), their Fig. 3, which has shown a more consistently ordered POS field for Chamaeleon-Musca compared to that of Ophiuchus.

For all three cases examined in PlanckXX, the values in the polarization fraction distribution are low compared to either Chamaeleon-Musca or Ophiuchus. It may be that the random component, and hence depolarization, in the MHD models was too strong. Alternatively, the forcing mechanism they used induced a Gaussian distribution to the magnetic field, while our analysis, with SN-driven turbulence, indicates that the magnetic field has a more exponential distribution, which could influence the efficiency of the depolarization (see Sect. 3.3). With respect to the mean field component, Gent (2012, Ch. 9 Fig. 9.12) found that the magnetic fields in the cold filamentary regions formed by SN driven turbulence are more regular than the ISM as a whole and that they are strongly aligned with the ambient warm ISM in which they are embedded. This is in contrast to the observations of Planck Collaboration Int. XXXII (2016), Planck Collaboration Int. XXXIII (2016), who found the orientation of the filamentary structure of the most dense molecular clouds to be perpendicular to the mean magnetic field. In the orientation of magnetic fields along dense molecular cloud filaments, gravitational collapse and runaway thermal instabilities leading to the formation of such high-density structures may be critical, which are absent on the scales considered in our MHD model.

Fig. 4 Left: joint 2D histograms of polarization fraction p and column density NH. Red dashed lines show the location of pmax = 19.8% from PlanckXIX. The green lines follow the maximum Cbin as a functionof NH, and the blue lines represent the weighted mean of p. Right: sample maps of p. The polarization fractions are weakened in the galactic midplane, which is seen in the 2D histograms as a growing distribution of low p at high NH. For Rmax = 1 kpc, see Figs. 13 and 14.

3.3 Polarization angle dispersion

Looking at the maps of polarization angle dispersion presented in Fig. 5 (right panels), we observe filamentary structures (-filaments) similar in appearance to PlanckXIX Fig. 12. However, for R_max = 0.25 kpc these are very large-scale structures, which span the full range of examined galactic latitudes in some locations and are much thicker than the PlanckXIX observations. With increasing R_max the -filament structure resembles very well PlanckXIX Fig. 12, becoming tangled and fragmented, increasing in number, and having more wiggles. PlanckXX report filament-like maps of from their synthetic observations. However, because of the differences in scale, resolution and forcing mechanism behind the MHD model used in PlanckXX and this study, an effective qualitative comparison is not reasonable.

In Fig. 5 (left panels) the joint 2D histograms of and p show some agreement with Fig. 23 of PlanckXIX. Angular dispersion is inversely proportional to polarization fraction and may be approximated by (15)

The 2D histograms show our best fits (black lines) for Eq. (15) and the fit from PlanckXIX (red lines). The parameters are summarized in Table 2. For R_max at 1 kpc and above, the slope α is near to the PlanckXIX fit, but by having smaller intercept β our 2D histograms are shifted to lower . There may be a case for inferring that the optimal integration length is about R_max = 1.5 kpc, if we are looking to match α for the PlanckXIX data. Repeating the analysis for an MHD domain with increased horizontal extent and/or increased resolution might indicate whether this is a robust physical relationship between the simulated and observed measurements. The ratio is averaged and listed in Table 2 for each R_max.

The dispersion values in our simulations increase towards PlanckXIX as R_max increases, but plateau at a ratio of 35% for R_max ≳ 1 kpc. The increase to 39%, when artificially strengthening the amplitudes of the small-scale fluctuations of the magnetic field, supports the view that the low simulated value for β is in part due to an insufficient small-scale field. The low values of in our MHD model may be attributed to truncated resolution at the grid scale of 4 pc.

To understand the distribution of velocity and magnetic fields in our MHD model data, we present probability density distribution functions (PDFs) of the both variables in Figs. 6 and 7. The PDFs are calculated over all MHD model cells from all snapshots where each cell, being of same volume, has an equal weight. In the top panel of Fig. 6, the PDF is exponential for B_x (B_y) and B_z, while B_y (B_x) is skewed by the global shear. The bottom panel of Fig. 6 shows the PDFs for the velocity components, where within the velocity range ± 400 km s⁻¹, the PDFs reflect ISM physics as resolved by the model. We evolve SN remnants only from the latter stages of the Sedov-Taylor phase, however, so this is exhibited in the truncated PDF for |u| ≳ 400 km s⁻¹.

The spiked PDF shown in Fig. 6 for the magnetic and velocity field components arise from the physics of repetitive shock-driven turbulence, independent of the model and resolution. Few authors have discussed this property, but a similar distribution for the velocity profile is illustrated in Gressel (2010, Fig. 3.13). In their Fig. 9 Mac Low et al. (2005), using approximately 1.5 pc resolution, show similar PDF profiles for the divergence of the velocity field, also supporting this physical interpretation of the turbulent structure of the ISM.

The multiphase medium plays a role in this distribution. In Fig. 7 we divide PDFs into three components, corresponding with cold, warm, and hot medium (T < 100 K, 100 ≤ T < 10⁵ K and T > 10⁵ K, respectively). The Fig. 6 velocity profile mainly exhibits the warm phase PDF depicted in Fig. 7, apart from at high velocities where the hot phase is more visible. The PDF is approximately Gaussian for the cold phase. Based on this component separation, the sharp PDF profile of magnetic field in the hot phase, and less so for the warm, is likely connected with the large-scale compressive forcing in the warm and hot phases and the subsequent turbulent cascade. The cold phase contains a magnetic field with a distribution that is between an exponential and a Gaussian, resembling that of PlanckXX, Fig. 11, apart from having a weaker magnetic field strength. It is possible that the higher densities in the cold phase may act as a sponge for these compressions and rarefactions, but we cannot assume that the effects of the turbulent cascade driven by SN are absent even at the scale of molecular clouds.

Therefore, some caution should be attached to the velocity PDF for the cold phase for at least two reasons. In their Fig. 15 c Gent et al. (2013a) showed that the hot phase flows are mainly subsonic, that is the warm transonic and cold supersonic. However, the cold clusters are typically entrained within the bulk flows of their ambient warm gas. If the bulk velocity of the ambientwarm gas were subtracted, then the Mach number of the cold phase would likely be reduced, and the residual flow might retain more of the PDF structure of the hot and warm phases. Also, in this model the scale of the cold structures tend to be only a few grid spaces across, that is they are near the limit of the model resolution, and hence much of the substructure of the magnetic and velocity fields in this phase are truncated. Therefore, although the MHD model here is truncated at 4 pc, it is our contention that the physics that drive the structure of the magnetic field in the hot and warm phases are still relevant to the flow driving the dynamics at smaller scales. This would require comparison with higher resolution multiphase turbulence simulations. We would only expect to introduce changes to the structure of the turbulence from new physical processes such as self-gravity.

Fig. 5 Left: joint 2D histograms of angle dispersion and polarization fraction with the effect of increasing Rmax. The red line depicts the fit to the observations as presented in PlanckXIX, . The black lines are the best fits to our simulated data, to which α and β are given in the Table 2. Right: sample maps of polarization angle dispersion with corresponding Rmax.

4 Shock and magnetic structure interpretation

We now consider how -filament structures seen in the polarization angle dispersion measurements are related to physical properties of the ISM. These are difficult to measure directly through observations, but can be measured easily in the MHD models. In the analysis that follows, we mostly refer to integration along the LOS with R_max = 1 kpc. This range, within the properties and horizontal extent of the MHD model, is sufficient to adequately capture the key features present in the PlanckXIX observations. For more a demanding analysis, it would be recommended to integrate R_max ≃ 2L_x.

Fig. 6 Probability density function combined from all 11 snapshots for the components of B and u.

4.1 -filaments compared with shocks

Changes in the direction of polarization angle ψ and therefore are related to changes in the magnetic field, and these are driven and generated by SN-driven turbulence. Generally, follows a lognormal distribution (see PlanckXIX Fig. 14). The lognormal nature of the distribution in the observed and simulated ISM appears consistent with the effect shocked turbulence has on the statistics of the gas density (as noted in Vazquez-Semadeni 1994; Elmegreen & Scalo 2004; Gent et al. 2013a), which encourages us to look into the connection between distribution and shocks present in our simulation data.

To investigate the effect of the shocks, we first compute a proxy of their magnitude, , where only the negative divergence contributes. This corresponds to regions where the flow is convergent, where therefore the shocks created by SNe are compressing the surrounding ISM. The values of C_{shock, local} are calculated within the numerical grid of the MHD model. All C_shock maps show C_{shock, local} averaged over the LOS up to the defined R_max, or (16)

In Fig. 8, maps are shown for average C_shock, , and the POS magnetic field, |B_POS|, averaged over the LOS up to R_max = 1 kpc. A zoom-in area is indicated on each map, for which the local maps are shown in Fig. 9. The joint 2D histograms of and C_shock are presented in the Fig. 10.

Apart from the energy input to the turbulence being stronger in the midplane because of the general distribution of SNe, there is no visible correlation between the SN shocks and the average POS magnetic field nor polarization angle dispersion. The 2D histograms presented in the Fig. 10 do not show any clear dependence between and C_shock, apart from the effect of midplane-weighted SN-distribution being more pronounced. However, there is a large-scale pattern in the magnetic field, which is discussed in the Sect. 4.2.

Examining the zoomed-in area shown in Fig. 9, the filamentary structures in (top right panel) overlap very sharply with the areas of low polarization fraction (top left panel). There is likely a connected phenomenon, which links these effects. Directly relating this to specific physical features in the model is a challenge because locally the alignment of polarization is a combination of effects layered on the top of each other. One approach to understanding this is to perform a series of calculations over a range of discrete integration lengths R_max and to analyse in detail how the maps change as certain features of the model are included or excluded. It is also possible to explore the relationship between magnetic fields and dispersion over large scales.

We stress that in advocating such an approach, we do not aim to explain directly particular observational features from Planck, but to inform how the inclusion or absence of physical phenomena along the LOS, which we can identify precisely in the MHD models, might be expected to contribute to the simulated observations. Such results can also be used to deepen the physical interpretation of the Planck findings.

When comparing the shock profile in Fig. 9 (bottom right panel) with the polarization angle dispersion (top right panel), the strongest filamentary structures correspond to locations where the shocks are negligible. In the top half of Fig. 9 (bottom left panel) the strength of the POS magnetic field seems to correlate very well with the polarization fraction (top right panel), however in the bottom right quadrant of the map a strong field is anti-correlated to p, so the relationship is not at all straightforward. In principle, all of these relationships should be explored further by varying R_max as described above, but it appears likely that we can exclude that the -filaments are indicative of the shock structure of the ISM.

Fig. 7 Probability density function combined from all 11 snapshots for the components of B and u with cold, warm and hot components separated.

4.2 Dispersion and the large-scale magnetic field

Heiles (1996), for example, already used starlight polarization measurement to estimate the direction and curvature of the Galactic magnetic field. Their results showed signs of spiral curvature, although the fluctuating component of the field is also significant. More recent studies (Planck Collaboration Int. XLII 2016; Planck Collaboration Int. XLIV 2016) have used the new Planck results to estimate the structural properties of the Galactic magnetic field by fitting magnetic field models to observations.

In our synthetic observations, the distribution of polarization angle dispersion exhibits a significant dependence on galactic latitude and longitude. This becomes more apparent when taking averages from all simulated observations with the same observer location, but for different snapshots, as is presented in Fig. 11 (bottom right panel). The measurements of Q reverse twice in one full rotation in latitude (top left panel), while U switches sign across the galactic plane and also exhibits the same latitudinal sign reversals as Q (top right panel). The polarization fractions are minimized where the brightest filamentary structures are most pronounced. The general nature of this pattern may be expected. As outlined in Sect. 2.1, the magnetic mean field is strongly aligned in the direction of the differential rotation of the model.

The dependence of polarization properties on the orientation of the mean magnetic field is a known relation. Already, in the context of turbulent molecular clouds, Ostriker et al. (2001), Soler et al. (2013), and PlanckXX have shown that distribution of fluctuations in polarization is connected with the direction of the magnetic mean fields of molecular clouds. Such analysis has also been used with Planck observation in relation to the molecular clouds (Planck Collaboration Int. XXXII 2016; Soler et al. 2016). This urges us to look into this phenomenon with our modelling results. However, as we look into large-scale patterns, Planck Collaboration Int. XLII (2016); Planck Collaboration Int. XLIV (2016) provide the most fruitful points of comparison.

The top panels of Fig. 11 are remarkably consistent with Planck Collaboration Int. XLII (2016, Fig. 13), who examine a set of Galactic magnetic field models without dynamo nor SN-driven turbulence to generate a realistic field, but including the Galactic centre and spiral arms, which they fit to observational data. Their synthetic maps do not have the small-scale features associated with the turbulence, but on larger scales our maps have very similar structure; the only characteristic difference is longitudinal variation caused by the included spiral arms. So, apparently, local structures in the ISM may be predominant in the observed data, although further investigation along the lines of Planck Collaboration Int. XLII (2016) and inclusion of spiral arms in a model of MHD turbulence would need to be pursued.

The large-scale variation follows from the presence of the mean field in our MHD model. We can illustrate it with a simple analytical example starting from the first principles. Let us assume a simple uniform y-directional mean field with random fluctuations at the smallest scales of the grid (17)

where |B₀ŷ|≫|Δb|. This configuration generates a large-scale structure of the polarization (see Fig. 12, and also Planck Collaboration Int. XLIV 2016, their Fig. 4, top panels). In addition to this, the direction of the magnetic field affects the sensitivity of the observed polarization to the small magnetic fluctuations. If we apply Eq. (17) to Eqs. (4)–(6), we notice that near the HEALPix coordinates ϕ ≈±π∕2 and θ ≈ π∕2 the influence of the magnetic field approaches the values ψ ≈ π∕2 + Δψ and γ ≈ π∕2 + Δγ, where we divided the contribution from the mean and fluctuating field. Therefore, when calculating the polarization components, with Eqs. (8) and (9), we get,

This signifies that, when the LOS approaches the direction of the consistent mean magnetic field, the POS field is highly sensitive to small, local fluctuations caused by turbulence. This, in turn, shows up as variations in polarization angle and therefore relatively high . To summarize, when the strong mean field is perpendicular to the LOS, its direction dominates the polarization angles, but when the field is parallel to the LOS, the observed polarization angles are more sensitive to the small fluctuations in the field. However, the polarization fraction p is weak in the mean field aligned with the LOS, as the small fluctuations themselves produce less strong Q and U. Thus, we have a similar interpretation to PlanckXX. In their study, is strongest when the POV faces towards the mean field direction, along with a weaker polarization fraction. In contrast, they observe a higher polarization fraction and a coherent polarization angle when the direction of the mean field follows the POS.

We note that B₀ŷ is distinct from as defined in Eq. (1). The parameter is defined by local averaging and includes varying x and y components, although it is most strongly aligned along y. In the MHD model, the fluctuations b have the same order of magnitude as , which complicates visualizing the large-scale mean field even with the simulated observations. For observations, the structure of the mean field is even more opaque, as the mean field direction is subject to large-scale diversions when interacting with spiral arms and the central bulge of the Galaxy. Therefore, using this interpretation to understand the observed by PlanckXIX is not trivial. In the case of our MHD simulation, the mean field is clearly stronger than the fluctuating field. In the case of our Galaxy, the general structure of the large-scale field is more complicated and the fluctuations are stronger (e.g. Rand & Kulkarni 1989; Haverkorn 2015). However, the effect of a large-scale magnetic field on could encourage further study, as long as measurement error considerations are taken into account (Montier et al. 2015; Alina et al. 2016).

Fig. 8 Top: distribution of the average values of Cshock in the LOS. Middle: map of polarization angle dispersion . Bottom: projected average magnetic field strength in POS. The green rectangles refer to the area featured in Fig. 9. Rmax = 1 kpc.

Fig. 9 Local comparison of p, , |BPOS |, and shocks within the area indicated in Figs. 14 and 8.

Fig. 10 Joint2D histograms of polarization angle dispersion and average LOS shock compression Cshock.

Fig. 11 Averaged maps of Q, U, p, over all 12 snapshots.

Fig. 12 Distribution of polarization assuming a simple analytical model of the magnetic field as described by the Eq. (17).

4.3 Effect of magnetic fluctuations

To assess the relevance of the mean and fluctuating field to the synthetic observations, we explore the impact of increasing the amplitude of the fluctuations relative to the mean field. Using the decomposition of the field as illustrated in Fig. 2, we double the amplitudes of the fluctuating component of the magnetic field and then sum it with the mean field to obtain a physically generated field with stronger fluctuations. This means that as for the original dataset we get and with the doubled amplitudes .

Doubling the amplitudes of the fluctuating field somewhat reduces p and increases , as can be seen by comparing the left and right joint 2D histograms in the top two rows of Fig. 13. The column density is not affected, which is understandable as the thermodynamic properties of the model are unchanged. The trends in the maximum count and weighted mean, traced onto the 2D histograms by the greenand blue lines, respectively, do not change, but they shift correspondingly with the general shift in the distribution. Comparing the joint 2D histograms of and p in the last row of Fig. 13, the log fit relating to p in the simulation shifts marginally closer to the PlanckXIX fit, but the slope of the line steepens. In this sense the original MHD model still better explains the relationship between and p. The increases in for the enhanced perturbations is associated mainly with low values of p, while values associated with high p ≳ 10% remain less affected by the enhanced perturbations.

From the maps of p in Fig. 14 (top row) the polarization fraction is in general not only dampened, but also exhibits increasingly fine structure. This is also evident for (bottom row), in which the filamentary structure for the perturbed field with doubled amplitudes (right) is highly fractured, compared to the map for the original field.

Increasing the relative strength of the fluctuations in the magnetic field does not bring us visibly closer to the values observed by PlanckXIX. Nor can a longer integration distance address this, although this does help with increasing column densities. The distribution of polarization fraction in the MHD simulations is spread slightly higher than in PlanckXIX, but doubling the strength of the random field component makes this a very good fit to PlanckXIX. Part of the contribution to comes from the mean field, which in the MHD model is far from uniform. In the Galaxy, spiral arms and the central bulge add to the variations in the mean field. Inclusion of such features in an MHD model would serve to enhance dispersion angles, but most likely the strongest factor is the limiting scale of the magnetic field fluctuations and gas density concentrations. These are truncated at the grid scale of 4 pc.

Given the long computation times and the large scales necessary to model a SN driven dynamo, increased resolution in the near future is likely to be modest. Nevertheless, the trends and characteristics exhibited with these simulations have much in common with the PlanckXIX results. Exploration of these simulated observations has helped to reveal how different aspects of the magnetic field contribute to the observations. A sweep of integration lengths R_max, comparing observations from many viewing angles and results across even a limited range of MHD model resolution may offer further valuable insights.

Fig. 13 Joint 2D histograms of (top) polarization fraction and column density, (middle) polarization angle dispersion and column density, and (bottom) polarization angle dispersion and fraction. The red line depicts the fit to the observations as presented in Planck Collaboration Int. XIX (2015), . The black lines are the best fits to our simulated data, to which α and β are given in the Table 2. The green lines follow the maximum counts in the 2D histograms as a function of NH, and the bluelines represent a weighted mean of the polarization fractions. Left: normal case. Right: doubled perturbations case. In all cases Rmax = 1 kpc.

5 Discussion and conclusions

In this paper, using MHD models supplemented with radiative transfer computations, we set out to study the effect on the polarization of dust in the ISM of the following physical ingredients:

• SN regulated multiphase ISM with a hot component and with longitudinal and latitudinal anisotropy;

• the presence of ubiquitous shock fronts driven by SNe;

• the presence of self-consistently generated inhomogeneous and anisotropic magnetic fields both by large- and small-scale dynamos.

Previous investigations have been limited to the two-phase ISM, including only the cold molecular and diffuse warm gas components, either without shocks or with artificially induced shock fronts, and with imposed magnetic field configurations.

We find a very good correspondence with the simulated maps, exhibiting astrongly filamentary structure, and the all-sky observations of PlanckXIX implying that our MHD models capture some essential features related to the formation of -filaments. In accordance with the observations, we find a good match to the anti-correlation between the polarization fraction p and polarization angle dispersion . The power law relation is very accurately reflected in the simulation, although differs by a factor 1/3 because it is sensitive to small-scale fluctuations and cold dense clouds, which our MHD model cannot sufficiently resolve.

The mean magnetic field has both a systematic orientation in the direction of the galactic shear and a non-uniform structure. This significantly affects the observed polarization properties. A strong POS mean field is found to dominate over the contribution of the small-scale component to the polarization angles in such a way that the observed is reduced. Conversely, when the mean field is parallel to the LOS, the observed polarization angles are more sensitive to the small-scale fluctuations. Because of its varying orientation, the mean field also partially contributes to the generation of -filaments. In the light of these general findings, our key results can be listed as follows:

• We have demonstrated a means of probing the relationship between the observations and physical features along the LOS by varying the integration length (R_max) in radiative transfer calculations.

• The value p is correlated with the strength of the mean field in the POS and is correlated positively with the fluctuating field b and inversely with p. We confirm the inverse correlation of the form , observed by PlanckXIX. Our results support a view that a high polarization fraction indicates a strong POS coherent magnetic field, while a low polarization fraction is consistent with a strong LOS mean field that is supposed to be approximately isotropic. It may be possible to apply this to the measurements of and p to make inferences about the strength and orientation of the mean field.

• The filamentary structure of becomes smaller scale, brighter, and more fragmentary with increasing R_max as depolarization accumulates along the LOS. The general occurrence of brightest -filaments are well correlated with the large-scale shifts in POS orientation of the magnetic field.

• The -filaments do not correlate with the column density nor the location of SN shocks, but can be attributed to the distribution of the small-scale magnetic field. This is because the small-scale magnetic fields are the result of a small-scale dynamo also present in the MHD model, enabling field generation throughout the domain.

• The fluctuations in the magnetic field are primarily driven by SN turbulence and follow an exponential distribution in the hot and warm phase medium, while the cold phase medium follows a more Gaussian distribution. This indicates that the methods of model fitting assuming a Gaussian random magnetic field might not be the most sensible choice at least for the diffuse ISM, where the warm and hot components play a part.

• Comparing joint 2D histograms of and column density N_H, we probethe relationship applying in real ISM between turbulent viscosity and . We tentatively assert that turbulent viscosity in the fully ionized ISM reduces slightly with increasing temperature as T^−λ, for some small positive λ. Further MHD simulations with varying models of viscosity are required to test this interpretation.

• We compare simulated polarization observations here and in Planck Collaboration Int. XLII (2016) with Planck observations. It appears likely that the strong variations in PlanckXIX observations are attributable to the physical structure of the ISM in the solar neighbourhood. Alternative models of the spiral structure and combining spiral arms with SN turbulence MHD models are required to explore this further.

• Increasing R_max above 1 kpc increases N_H and reduces p in line with observed distributions, suggesting that there is a minimum value for R_max that is needed for simulated observations. However, to exclude artificial artefacts, we have to limit R_max ≲ 2L_x, where L_x is the model horizontal extent.

For future work, MHD models with increased resolution and/or horizontal extent to probe the effects of longer integration ranges and higher densities and fluctuations on and N_H would be helpful. As these are essential for many scientific priorities, these opportunities will undoubtedly be fulfilled. Even with the current simulations the role of the fluctuating magnetic field can be further investigated by retaining the existing MHD field and adding to this an appropriate isotropic fluctuating field of the correct magnitude with disturbances at various smaller scales. Another way forward would be to use the saturated stages of the MHD models, re-mesh these stages to include denser grids and thereby finer scales, and re-run only up to a new saturated stage, eliminating the long dynamo evolutionary stage.

Conducting a series of experiments probing a range of R_max to analyse how physical features are being captured by the simulated observations in relation to the POV would shed more light on how polarization features are connected with the magnetic field structure. Including spiral arms in MHD simulations and exploring how this impacts on the structure of the magnetic field and anisotropy in the synthetic polarization observations could improve the interpretation of model results in terms of the Milky Way Galaxy. Yet another possibility is to focus on zoomed-in features of the MHD models with similarities to the observed features from Planck, and this way to explore what we can learn about the 3D structure of the magnetic field at that location.

The often used Davis-Chandrasekhar-Fermi method (Davis 1951; Chandrasekhar & Fermi 1953) allows for the determination of the POS magnetic field strength if the dispersion of polarization angle is known. It is especially useful for estimating magnetic field strengths in regions, such as molecular clouds, where the Zeeman effect is weak. An important direction of future work is to test the predictions of this method against self-consistently generated large- and small-scale fields, ideally when the models can be refined to reach the limit where both quantities resemble their observed counterparts.

As we have shown in this paper, the interplay of the large- and small-scale magnetic fields can cause systematic effects in the polarization measures that may well be used to map the mean magnetic field of the Milky Way. The ratio of and p, reacting to the presence of different levels and orientations of the magnetic field components, may be used as a tracer of the orientation of the mean field and, on the other hand, of the ratio of the strengths of large- and small-scale magnetic fields.

Fig. 14 Top: sample maps of polarization fraction p. Bottom: sample maps of polarization angle dispersion . Left: normal case. Right: doubled perturbations case. In all cases Rmax = 1 kpc.

Acknowledgements

MSV, FAG, and MJK acknowledge support of the Grand Challenge project SNDYN, CSC – IT Center for Science Ltd.(Finland) and the Academy of Finland Centre of Excellence ReSoLVE (project number 272157). The simulations were performed using the supercomputers hosted by the CSC – IT Center for Science Ltd. in Espoo, Finland, which is administered by the Finnish Ministry of Education. MJ acknowledges the support of the Academy of Finland Grant No. 250741. MSV thanks University of Helsinki and the Jenny and Antti Wihuri Foundation for financial support. Some of the results in this paper have been derived using the HEALPix package. This research has made use of NASA’s Astrophysics Data System.

Appendix A Maps of polarized emission

We feature maps of Stokes I, Q, and U used for calculating p and shown by the maps in Figs. 4, 5, 8, and 14. See Figs. A.1–A.3 for R_max = 0.25, 1.0, 2.0, and 4.0, respectively.

Fig. A.1 Stokes I, Q, and U with Rmax = 0.25 kpc.

Fig. A.2 Stokes I, Q, and U with Rmax = 1 kpc.

Fig. A.3 Stokes I, Q, and U with Rmax = 2 kpc.

Fig. A.4 Stokes I, Q, and U with Rmax = 4 kpc.

Appendix B The MHD model

The focus of this paper is to explore the unique insights possible by using physically driven multiphase ISM and dynamo generated simulation magnetic field data to derive synthetic observations. These can then be compared with astronomical observationsand those of simulations based on imposed flows and magnetic field configurations. We stress that the MHD simulations are not designed to model features in the ISM that have been observed or modelled by Planck or other observers, but independently to examine the multiphase structure and dynamo action of SN-driven turbulence in the ISM. These simulations have been described in detail elsewhere, but without breaking the focus on the core purpose of this paper, we clarify more of the basic assumptions and techniques relating to the MHD models.

In this particular model SNe are modelled by injecting thermal and kinetic energy of 10⁵¹ erg into a spherical region, from which the pressure gradient and momentum drive shock fronts and heating into the ambient ISM. The thermal component has been found to be critical for efficient generation of vorticity through baroclinicity (Käpylä et al. 2018), and this is beneficial for dynamo action, unlike models applying momentum forcing, which are dominated by potential flows (Iffrig & Hennebelle 2017). The timing and location are motivated by observational estimates of SN distribution in the solar neighbourhood and placement is not dependent on the local ISM density or temperature, avoiding excessive damping in dense material or the accumulation of immovable overdense regions. The environment applies an initial vertical density profile based on estimates summarized in Ferrière (2001), subject to gravitational acceleration following Kuijken & Gilmore (1989) and differentially rotating with angular speed and rotation consistent with the solar neighbourhood. The resulting turbulence, magnetic field structure, and the multiphase structure of the ISM are thus determined by the solution of the physical equations and random SN events, and not imposed. The system of equations being solvedinclude continuity, momentum, energy, and induction treating the ISM as a monatomic ideal gas, i.e.

with the symbols having their usual meaning. For a full explanation refer to the earlier papers.

Modelling a galactic dynamo required an integration time in excess of 1 Gyr with a time step below 100 yr, and on occasions less than 1 yr; our chosen resolution is 4 pc along each edge. A lower resolution would permit faster progress with reduced resources, but careful experiments with various resolutions established 4 pc provided the minimum resolution with our numerical scheme to faithfully evolve SN remnants with the relevant essential physics (Gent 2012, Appendix A). There is a seed magnetic field of a few nG, which is amplified by the SN driven turbulence to saturate with a mean field strength a few μG. Unlike most other similar models with imposed magnetic fields and unphysical forcing, the field structure is derived from the physical processes modelled. The motivation for the 1 kpc horizontal extent was to ensure no artificial self-feedback from SN superbubbles, which can expand as far as 500 pc, and to include all estimated dynamo modes. Some similar models apply AMR with local resolution near 1 pc, although elsewhere within the models the resolution can be an order of magnitude weaker, but they are not seeking dynamo solutions. Improving the resolution fourfold would multiply the size of the calculation at each iteration by 4³ and increase the number of iterations per Myr by at least 4². A detailed analysis of the multiphase structure and properties of the ISM obtained with this model is reported in Gent et al. (2013a), and these are very consistent with results reported by other authors adopting higher resolution and employing AMR (Joung & Mac Low 2006; de Avillez & Breitschwerdt 2007; Hill et al. 2012). Other than obtaining higher maximal values for density and higher proportion of cold gas due to the enhanced resolution, the phases obtained are still similarly located and distributed, and exhibit the same temperature and velocity characteristics.

The reason for this is that the difference between 1 and 4 pc does not permit the inclusion of qualitatively distinct physics. Even a tenfold increase in maximal densities would not bring the model within the threshold demanding self-gravity. The Jeans mass for gas at 100 cm⁻³ at 100 K is 4.4 × 10³³ kg. The total in a volume element is less than 10³¹ kg, hence none of the regions within the model are close to being susceptible to gravitational instability, and this would be the case even for models with ten times more dense regions. The typical cooling times in our MHD model are about 100 yr, but occasionally as low as 1 yr; this is never sufficiently dense to introduce unresolved thermal instabilities, which lead to cloud collapse, and these processes are also absent for models with 1 pc resolution. Even where isolated regions within such models may approach the critical Jeans mass, they are unlikely to have enough significance to the global dynamics to justify the numerical overheads.

Increased resolution is, of course, desirable, but comes with a numerical cost. Higher resolution would increase the Reynolds numbers applicable in the model, permitting smaller magnetic and velocity structures to be resolved and perhaps more efficient dynamo action. Further illustrations, including video representations, of the simulation data are available online³.

References

Movies

This movie traces the density of the ISM. (Access here)

This movie shows the temperature. (Access here)

This movie depicts the magnetic energy distribution. (Access here)

This movie illustrates the azimuthal component of the magnetic field. (Access here)

Images below show a snap shot of the ISM simulation, where the system has not yet grown fully turbulent.

The animation traces x-planes through the computational domain.

Density. (Access here)

Temperature. (Access here)

Vertical velocity. (Access here)

All Tables

All Figures

	Fig. 1 Left: representative slice in the xz-plane of gas number density n [1∕cm−3] from a singlesnapshot overplotted with contours of temperature T [K], illustrates the multiphase ISM. Right: snapshot of the model ISM with temperature slices in the background and isosurfaces of the density in the foreground; fieldlines of the magnetic field are overplotted.
In the text

	Fig. 2 Snapshot of the model ISM, with temperature slices in the background and isosurfaces of the density in the foreground, as in Fig 1, but with the total magnetic B fieldlines replaced by fieldlines (left) and b fieldlines (right). Further visualization of the MHD, including video representations, is available online. (the full set is provided at http://fagent.wikidot.com/astro).
In the text

	Fig. 3 Left: joint 2D histograms of polarization angle dispersion and column density NH. The green lines follow the maximum Cbin as a functionof NH, and the blue lines represent the weighted mean of . Right: sample maps of column density. For the case Rmax = 1 kpc, see Fig. 13.
In the text

Fig. 4 Left: joint 2D histograms of polarization fraction p and column density NH. Red dashed lines show the location of pmax = 19.8% from PlanckXIX. The green lines follow the maximum Cbin as a functionof NH, and the blue lines represent the weighted mean of p. Right: sample maps of p. The polarization fractions are weakened in the galactic midplane, which is seen in the 2D histograms as a growing distribution of low p at high NH. For Rmax = 1 kpc, see Figs. 13 and 14.

In the text

	Fig. 5 Left: joint 2D histograms of angle dispersion and polarization fraction with the effect of increasing Rmax. The red line depicts the fit to the observations as presented in PlanckXIX, . The black lines are the best fits to our simulated data, to which α and β are given in the Table 2. Right: sample maps of polarization angle dispersion with corresponding Rmax.
In the text

	Fig. 6 Probability density function combined from all 11 snapshots for the components of B and u.
In the text

	Fig. 7 Probability density function combined from all 11 snapshots for the components of B and u with cold, warm and hot components separated.
In the text

	Fig. 8 Top: distribution of the average values of Cshock in the LOS. Middle: map of polarization angle dispersion . Bottom: projected average magnetic field strength in POS. The green rectangles refer to the area featured in Fig. 9. Rmax = 1 kpc.
In the text

	Fig. 9 Local comparison of p, , |BPOS |, and shocks within the area indicated in Figs. 14 and 8.
In the text

	Fig. 10 Joint2D histograms of polarization angle dispersion and average LOS shock compression Cshock.
In the text

	Fig. 11 Averaged maps of Q, U, p, over all 12 snapshots.
In the text

	Fig. 12 Distribution of polarization assuming a simple analytical model of the magnetic field as described by the Eq. (17).
In the text

Fig. 13 Joint 2D histograms of (top) polarization fraction and column density, (middle) polarization angle dispersion and column density, and (bottom) polarization angle dispersion and fraction. The red line depicts the fit to the observations as presented in Planck Collaboration Int. XIX (2015), . The black lines are the best fits to our simulated data, to which α and β are given in the Table 2. The green lines follow the maximum counts in the 2D histograms as a function of NH, and the bluelines represent a weighted mean of the polarization fractions. Left: normal case. Right: doubled perturbations case. In all cases Rmax = 1 kpc.

In the text

	Fig. 14 Top: sample maps of polarization fraction p. Bottom: sample maps of polarization angle dispersion . Left: normal case. Right: doubled perturbations case. In all cases Rmax = 1 kpc.
In the text

	Fig. A.1 Stokes I, Q, and U with Rmax = 0.25 kpc.
In the text

	Fig. A.2 Stokes I, Q, and U with Rmax = 1 kpc.
In the text

	Fig. A.3 Stokes I, Q, and U with Rmax = 2 kpc.
In the text

	Fig. A.4 Stokes I, Q, and U with Rmax = 4 kpc.
In the text