© ESO 2022

1 Introduction

In order to derive stringent constraints on the cosmological parameters from the study of the cosmic microwave background (CMB) radiation, the Planck satellite mapped the whole sky at 353 GHz in both intensity and polarization (e.g., Planck Collaboration I 2020). At this frequency, the polarized sky is dominated by the polarized thermal emission from dust grains of the magnetized interstellar medium (ISM) of the Galaxy. Because most of the cosmological information can be obtained from the CMB radiation through the analysis of its T, E, and B angular power spectra (intensity and linear polarization) and their correlation (e.g., Kamionkowski et al. 1997b,a; Zaldarriaga & Seljak 1997), the dust polarization sky has been characterized using that metric. The characterization of the dust polarization auto- and cross-angular power spectra has led to several properties that intrigued the community (Planck Collaboration XXX 2016; Planck Collaboration XI 2020). Among those properties, the most significant ones are the observed values for the E∕B asymmetry and the TE correlation; that is to say, that the power in the E-mode polarization is about twice than that of the B modes and that the E modes show a positive cross-correlation with the total intensity (T).

It was quickly hypothesized that these features in the polarization power spectra should be related to the physics of the ISM, at least to some extent and on certain angular scales. In particular, it has been argued that they might reflect the correlation between the orientations of anisotropic density structures and magnetic field lines, which is expected from magneto-hydrodynamic (MHD) physics in some regimes. This is an interpretation that is supported by observational evidence of such alignments (e.g., Clark et al. 2014; Planck Collaboration XXXVIII 2016).

Several authors have invoked the turbulent properties of the ISM as the physical origin of these statistical features. Caldwell et al. (2017) and Kandel et al. (2017) developed theoretical methods to derive polarization power spectra of thermal dust emission from MHD turbulence theory. Both works suggest that the observed E∕B power asymmetry and TE correlation could be predictive of the MHD turbulence parameters in the ISM. In fact, Kritsuk et al. (2018) demonstrate that dust polarization maps synthesized from multiphase MHD turbulent simulations of the local ISM (Kritsuk et al. 2017) can naturally lead to the observed E∕B asymmetry. Along the same lines, Kim et al. (2019) analyzed the power spectra of synthetic dust polarization maps from simulations of multiphase, supernova-driven MHD turbulence in a kiloparsec-sized stratified shearing box (Kim & Ostriker 2017). They report an E∕B asymmetry and a TE correlation roughly consistent with observations, although at somewhat lower values. Interestingly, they also report fluctuations of the power spectrum characteristics with the observer’s galactic environment (simulated by varying the differential rotation rate of the box) and with temporal variations of the ISM properties due to bursts of star formation. Brandenburg et al. (2019) show that helical turbulence induces strongly asymmetric distributions of E-mode contributions along the line of sight whereas B-mode contributions are distributed more symmetrically. In this view, the E∕B asymmetry naturally stems from the fact that, on average, for a given line-of-sight, there is less cancellation of E- compared to B-modes. This observation remains true on large scales if the Sun is embedded into a large-scale helical field (Bracco et al. 2019a). According to the evidence gathered so far, the dust polarization power spectra contain information about the properties of the magnetized ISM. However, the variance induced by the observer’s viewpoint (also referred to as cosmic variance) might be strongly diluting this information.

In this paper, we quantify the effect of cosmic variance on the E∕B power asymmetry and the TE correlation. For this purpose, we use a set of global MHD simulations of Milky-Way-sized galaxies to synthesize full-sky observationsof polarized thermal dust emission at 353 GHz from different observer positions and study the corresponding power spectra. With these calculations we are not aiming to replicate (or model) the Planck observations. Instead, we are looking for observable trends based on physical processes, which we can then use to interpret polarization power spectra of Galactic observations.

Our input galaxy models are the six simulations presented in (Ntormousi 2018, hereafter N18). They include different initial magnetic fields (both in topology and strength) and are realized with or without feedback. Although the selected models are not exact replicas of the Milky Way, they cover a range of parameters wide enough to address our objective: different initial conditions for the magnetic field, which represent different realizations of the ordered component of the magnetic field; the presence or not of feedback, which allows us to explore the impact of the turbulent component of the magnetic field on the observables; and a complex three-dimensional environment in which to move the hypothetical observer. Among the (few) galactic models with magnetic fields available to the community, these models have the additional benefit of having a divergence-free magnetic field.

In Sect. 2 we present the galaxy simulations, the polarized thermal dust emission model, and our method to synthesize dust polarization maps from the simulation grid. In Sect. 3 we present the computation of the polarization power spectra and define the quantities of interest. Section 4 contains our main results: the measurements of the E∕B asymmetry and TE correlation in the whole sample and per galaxy setup. Their we also explore the viewpoint-induced variance in more details. Finally, Sect. 5 contains summarizing comments and conclusions.

This document contains appendices: Appendix A illustrates the convention used in this work to synthesize maps from a simulation performed on an Adaptive Cartesian grid. Appendix B explains details of the ray-tracing algorithm we developed. Finally, in Appendix C we show that creating synthetic maps using interpolation schemes can lead to artifacts in the resulting power spectra.

2 Synthetic dust polarization maps

2.1 MHD galaxies

2.1.1 Numerical code

The galaxy simulations, presented in N18, are performed with the publicly available MHD code RAMSES (Teyssier 2002), which solves the ideal MHD equations on a Cartesian grid and has Adaptive Mesh Refinement (AMR) capabilities. RAMSES uses a constrained transport scheme to evolve the magnetic field B, which guarantees ∇⋅ B = 0 always (Fromang et al. 2006). This provides a significant advantage as compared to codes that rely on divergence-cleaning algorithms which are known to create spurious effects in studies of turbulent environments (Balsara & Kim 2004).

2.1.2 Setup

The initial conditions for the MHD fluid, the stellar and dark matter particles were created using DICE (Perret et al. 2014; Perret 2016). They represent a Milky-Way-like galaxy (total mass M_tot = 2 × 10¹² M_⊙) at redshift zerowith different initial morphologies and strengths of the magnetic field. The virial velocity of the galaxy is 200 km s⁻¹, the mass fraction in stars is about 4.5% (including a stellar bulge with a mass fraction of 0.5%) and the mass fraction in gas is about 1%. The dark matter halo follows a NFW profile (Navarro et al. 1996), while the gas and stars are initially placed in exponential disks with a scale length of 9 kpc. The galactic disk contains two spiral arms, starting at 2 kpc and ending at 12 kpc. The mean gas temperature is set at 10⁴ K and subsonic turbulence, with an rms value of 8 km s⁻¹, is introduced throughout the disk.

The base (coarsest) simulation grid has a resolution of 256³ cells filling the box of 60 kpc³. Four AMR levels are used to capture the disk dynamics, and an additional AMR level in regions of star formation. In physical units, the highest resolution corresponds to 29 pc and the lowest to 234 pc.

Table 1 contains the initial parameters of the different galaxy models. The parameters varied are the strength and the initial topology of the magnetic field. The initial field is either toroidal (labels start with T), with a scale height and scale length of 1 kpc, or a poloidal field (labels start with P), with a scale height of 1 kpc and a scale length of 2 kpc. The poloidal magnetic field is model C from Ferrière & Terral (2014). Star formation is simulated in all models by forming sink particles when the density exceeds 1000 cm⁻³. Models whose label contain (“nf”) “f” (do not) include stellar feedback from supernovae, resulting from previously formed sink particles with a time delay of 3 Myr. Supernovae are implemented by injecting thermal energy into the cells around the sink particle according to the number of supernovae estimated for the predicted size of the formed stellar cluster. Table 1 also contains the physical time of the outputs we use in this work and the mass-weighed median of the magnetic field strength at that time. Unfortunately, evolving these galaxies at this resolution is very computationally demanding, so we only have one stellar generation (30 Myr).

The power spectra of magnetic and kinetic energy of these galaxies are given in N18. Histograms and cumulative histograms of the mass-weighed magnetic field strength are shown in Fig. 1 for the six galaxies.

2.2 Synthetic dust polarization maps

2.2.1 Synthesis

Mapping 3D data to 2D through line-of-sight integration generally leads to artifacts and systematics that depend on the integration method. In the specific case of mapping a (nonuniform) Cartesian onto a spherical grid, certain regions are over- or under-sampled and spurious effects appear.

One approach for removing these artifacts is to interpolate the values of neighboring cubic cells at the positions of the spherical grid. However, this process can lead to a mixing of angular scales, which in turn has a measurable effect on the power spectra (see Appendix C for a detailed explanation). Therefore, instead of interpolating, here we sample the Cartesian grid with a spherical grid sufficiently dense to include all the mass in the simulation. Going one step further in this direction, we develop a ray-tracing algorithm that computes, for each line of sight, the path length through each intervening cell of the Cartesian grid. Our algorithm, presented in more detail in Appendix B, maps 3D Cartesian grids on to 2D maps that follow an HEALPix tessellation with a resolution parameter N_side (Górski et al. 2005). N_side is the only free parameter of our map-making process.

To produce the synthetic maps used in this work, we place an observer in the xy plane of the simulated galaxies at a radial distance of 8 kpc from the center, that is roughly the distance from the Sun to the Galactic center. From each observer location (see Sect. 2.2.3) we integrate the simulations to produce 2D maps of the required observables using the algorithm described above.

We fix the targeted map resolution setting the HEALPix parameter N_side = 128. Technically, because the minimum side length of the cells populating the galactic disk is 29 pc and the maximum distance along the line-of-sight is about 50 kpc, ensuring that every single cell is crossed by at least one line of sight would require an N_side of 2048. However, by construction, the smallest cells are found toward high density regions and, therefore, mostly in the galactic disk. Since we are about to disregard the disk by masking out the bright (mostly equatorial) regions, an N_side of 128 is sufficient for the purposes of our analysis. Indeed, an N_side of 128 allows us to ray-trace the smallest cells (29pc) up to a distance of about 3.5 kpc (i.e., at least one sightline goes through each cell). Toward the poles, such a distance is about an order of magnitude larger than the scale height of the disk inferred in our simulations. Of course, going at lower latitudes, material at a lower distance from the disk (lower |Z|), could be missed. We checked that synthetic maps computed with an N_side of 512 (effectively “missing” no material up to 15 kpc distance, but also tracing each voxel with 16 times more sightlines) do not show difference in the sky regions of interest. Subtle differences were spotted only in the disk and toward the galactic center, both of which we disregard in this paper.

Fig. 1 Normalized distribution (top) and cumulative distribution (bottom) of the mass-weighed magnetic field strength (log10 (|B|)) in the disk region for each simulation as indicated by the legend. We define the disk region as a vertical cylinder centered on the center of the simulation box, with radius of 15 kpc and extending up to 500 pc above and below the plan Z = 0.

2.2.2 Polarized dust emission model

Each cell of the grid contains information on the matter density and the magnetic field. We use these physical quantities to produce synthetic maps of the thermal dust polarized emission as seen at 353 GHz.

We start from the integral equations for the Stokes I, Q, and U similar to those given by Lee & Draine (1985), Planck Collaboration Int. XX (2015) (Appendix B) and reviewed in Pelgrims et al. (2020). For optically thin emission at frequency ν and following the HEALPix convention for the polarization position angle¹, we write: $$I={\displaystyle \int S_{\nu }}\left[1-p_0\left({\cos }^2\gamma -\frac{2}{3}\right)\right]n_{\text{H}}{\sigma }_{\text{H}}\text{d}s,$$(1) $$Q={\displaystyle \int p_0}{S}_{\nu }\cos \left(2\varphi \right){\cos }^2\gamma n_{\text{H}}{\sigma }_{\text{H}}\text{d}s,$$(2) $$U={\displaystyle \int p_0}{S}_{\nu }\sin \left(2\varphi \right){\cos }^2\gamma n_{\text{H}}{\sigma }_{\text{H}}\text{d}s,$$(3)

where the integrals are computed along the line of sight over the emitting region (here, the full simulation box); S_ν is the source function, p₀ a parameter related to dust polarization properties combining grain asymmetric cross sections and the degree of alignment with the magnetic field, n_H the gas density, σ_H the dust cross-section per hydrogen atom averaged over angles, γ the angle of the local magnetic field to the plane of the sky, and ϕ the local polarization angle (see Fig. 14 in Planck Collaboration Int. XX 2015). Following Planck Collaboration Int. XLIV (2016) and Kim et al. (2019), we adopt a universal value for p₀ = 0.2, the black-body source function B_ν with a constant dust temperature of 18 K and a dust opacity at 353 GHz of σ_d,353 = 1.2 × 10⁻²⁶ cm⁻². A constant dust-to-gas ratio is assumed.

In this work, we want to understand the influence of the ISM conditions on the observables. Therefore, we neglect any variations in the emissivity of the dust grains throughout the Galaxy (despite known observational evidence for slight variations in the Milky Way; Finkbeiner et al. 1999; Planck Collaboration XI 2014; Pelgrims et al. 2021). We also assume that the dust grain alignment physics is the same throughout the Galaxy, an assumption that is valid on Galactic scales (Reissl et al. 2020; Seifried et al. 2019; Vandenbroucke et al. 2021).

It is important to notice that, unlike in works using parametric modeling, such as Planck Collaboration Int. XLIV (2016) and Pelgrims et al. (2020), the parameter p₀ is not intended to account for variations in the magnetic field orientation along the line of sight (in the integral step). The reason is that, unlike parametric modeling, direct numerical simulations follow the evolution of fluctuations self-consistently. Additional fluctuations at the subgrid level would interfere with the physics encoded in the simulations, and they would alter the power spectra and presumably their synthesis characteristics (Reissl et al. 2019) in a way that depends on the type of fluctuations and on the implementation (see, e.g., Sect. 4.3 of Wang et al. 2020).

2.2.3 Map making

For each of the six simulation setups discussed above, we synthesize 72 sets of polarization maps (I, Q, and U) by allowing the observer to move in the xy plane and make a complete circular excursion around the galactic center, at 8 kpc distance and with an angular step of 5°. In total we generate 432 sets of polarization maps. The used conventions for the rotation of the observer are explained and illustrated in Appendix A.

We use our ray-tracing algorithm, discussed in Appendix B, to generate the synthetic dust polarization maps at the resolution of N_side = 128. This choice of resolution fixes a set of 196 608 lines of sight, corresponding to pixels with angular size of about 0.46 degrees. To compute the integrated observables (Eqs. (1)−(3)) we split the integrals into the contribution from each cell. For each cell the algorithm computes the lengths of all the line-of-sight segments (corresponding to a map at fixed N_side) that cross it. The lengths are used to weigh the observables using the midpoint rule. After rotation and projection of the AMR grid on the observer spherical coordinate system (see Appendix A), the contribution from a given cell, labeled i, to the integrated Stokes parameters for the line of sight subtended by e_r is given by $$\begin{array}{lll}{I}_{e_r}^i\hfill & =\hfill & {\lambda }_{e_r}^i{n}_{\text{d}}^i\left(1-p_0\left(\frac{B_{\theta }^i{}^2+B_{\varphi }^i{}^2}{|B^i{|}^2}-2/3\right)\right)\hfill \\ Q_{e_r}^i\hfill & =\hfill & {\lambda }_{e_r}^i{n}_{\text{d}}^i{p}_0\frac{\left(B_{\theta }^i{}^2-B_{\varphi }^i{}^2\right)}{|B^i{|}^2}\hfill \\ U_{e_r}^i\hfill & =\hfill & -2{\lambda }_{e_r}^i{n}_{\text{d}}^i{p}_0\frac{\left(B_{\theta }^i{B}_{\varphi }^i\right)}{|B^i{|}^2}\hfill \end{array}$$

where $n_{\text{d}}^i$ is the dust density in cell i and ${\lambda }_{e_r}^i$ is the path length of the given line of sight through that cell i. The path lengths ${\lambda }_{e_r}^i$ are computed through our ray-tracing algorithm as explained in Appendix B. No further post-processing is applied to the synthetic maps for the power spectrum analysis. As a result, projected cubic patterns can be spotted on some sets ofmaps. Interpolation would produce visually better maps but not without affecting the power spectra (Appendix C). An example set of synthetic polarization maps is shown in Fig. 2 for the Tsf galaxy.

3 Power spectrum analysis

In this section we follow the same statistical characterization as that of the Galactic thermal dust emission performed in the context of CMB foregrounds by producing polarization power spectra of the synthetic dust maps. We also rename the intensity map T ≡ I to conform with the relevant literature.

3.1 Formalism

To the Stokes Q and U maps correspond rotation-invariant quantities, the E (gradient) modes and B (curl) modes which are the even and odd parts of the polarization vector field under parity transformation (Kamionkowski et al. 1997b; Zaldarriaga & Seljak 1997; Hu & White 1997). The statistical description of these three scalar and pseudo-scalar quantities (T, E, and B) is commonly based on their auto- and cross-angular power spectra as a function of multipole, $C_{\mathcal{l}}^{XY}$, where X and Y refer either to T, E, and or B and ℓ is the multipole number. A given ℓ roughly corresponds to an angular scale α ≈ 180°∕ℓ. Following the commonly adopted formalism, we carry out our analysis using the pseudo power spectra: ${\mathcal{D}}_{\mathcal{l}}^{XY}=\mathcal{l}(\mathcal{l}+1)C_{\mathcal{l}}^{XY}/(2\pi )$. (See e.g., Bracco et al. 2019b for a recent review of the formalism.)

The E∕B power asymmetry is measured through the ${\mathcal{R}}_{EB}$ ratio, which is obtained by averaging the ratio of the auto-power spectra ${\mathcal{D}}_{\mathcal{l}}^{EE}$ and ${\mathcal{D}}_{\mathcal{l}}^{BB}$ over a specified multipole range $$\text{begin{align}}{\mathcal{R}}_{EB}\equiv \text{ left}<\frac{{\mathcal{D}}_{\mathcal{l}}^{EE}}{{\mathcal{D}}_{\mathcal{l}}^{BB}}\text{ right}>\text{nnend{align}}$$(7)

where $\langle \cdot \rangle$ stands for the mean over multipole bins.

To quantify the correlation between power spectra we use the normalized parameter $r_{\mathcal{l}}^{XY}$ introduced by Caldwell et al. (2017). $r_{\mathcal{l}}^{XY}$ takes values 1, − 1 and 0 in case of perfect correlation, perfect anticorrelation and absence of correlation, respectively, and is defined as $$r_{\mathcal{l}}^{XY}=\frac{C_{\mathcal{l}}^{XY}}{\sqrt{C_{\mathcal{l}}^{XX}{C}_{\mathcal{l}}^{YY}}}.$$(8)

Even though we used all three $r^{XY}=\langle r_{\mathcal{l}}^{XY}\rangle$ (with XY = {TE, TB, EB}) to generate theidealized power spectra for the mask validation (see Sect. 3.3.3), we only consider the r^TE parameter to characterize our power spectra.

Fig. 2 Example of synthetic polarization maps of thermal dust emission. The observer is immersed in the Tsf galaxy at 8 kpc from the galactic center and at angular coordinate ϕ⊙ = 300°.

3.2 Computation from synthetic maps

We use the Xpol code (Tristram et al. 2005), utilized in Planck analyses (e.g., Planck Collaboration XXX 2016; Planck Collaboration XI 2020), to compute the pseudo power spectra of the dust polarized sky and account for incomplete sky coverage. Since our synthetic maps do not contain noise, we consider only the sampling variance as a source of uncertainty in the power spectrum estimates. The uncertainty is thus linked to the number of unmasked pixels and their spatial arrangement on the sky. Under Gaussian approximation these uncertainties estimated for a defined multipole bin (ℓ_bin) and retained sky fraction (f_sky) can be computed analytically using $${\left({\sigma }_{{\mathcal{D}}_{\mathcal{l}}^{XX}}\right)}^2=\frac{2}{(2{\mathcal{l}}_{\text{bin}})f_{\text{sky}}\Delta {\mathcal{l}}_{\text{bin}}}{\left({\mathcal{D}}_{{\mathcal{l}}_{\text{bin}}}^{XX}\right)}^2$$(9)

where X = {T;E;B} and Δℓ_bin is the width of the multipole bins (Planck Collaboration XXX 2016).

For each set of polarization maps we consider six nested regions at high galactic latitudes with retained sky fractions ranging from 20 to 70% by step of ten per cent. Those regions are determined independently for each set of polarization maps according tothe masking procedure described below. Then, for each of the 2592 sets of masked polarization maps we evaluate thepolarization power spectra using Xpol. The power spectra are estimated in the multipole range of $\mathcal{l}\in \left[40,370\right]$ which are split in bins with a width of 20. However, as discussed below, a shorter range of ℓ is determined and used in our analysis so as to allow for a systematic and homogeneous study of our sample and to easily perform any comparison.

To determine the high-galactic latitude sky regions on which to estimate the power spectra, we proceed as follows. For each setof maps, we start with the intensity map smoothed to a five degree resolution (FWHM) using a symmetric Gaussian beam. We then successively mask all pixels brighter than threshold values defined such that the retained skyfractions range from 0.2 to 0.7 by step of 0.1. To avoid power leakage we apodize the masks using a 2.5° FWHM beam. Each mask is then applied to the corresponding set of synthetic polarization maps and the polarization power spectra computed through the use of Xpol.

3.3 Full sample of polarization power spectra

3.3.1 Power-law power spectra

Upon visual inspection, most of the EE and BB auto-power spectra present a steep power-law dependence on ℓ, at least within a certain multipole range (see Fig. 3 for an example). Such a power-law dependence is also observed in power spectra of the real sky (e.g., Planck Collaboration XXX 2016) and of MHD simulations (e.g. Kim et al. 2019). We make use of least-square fits to the power spectra with a power-law function of the form² $${\widehat{\mathcal{D}}}_{\mathcal{l}}^{XY}=A_{80}^{XY}{\left(\frac{\mathcal{l}}{80}\right)}^{{\alpha }_{XY}+2}$$(10)

where XY can be TE, EE, and BB. $A_{80}^{XY}$ and α_XY are the amplitudes of ${\mathcal{D}}_{\mathcal{l}}^{XY}$ at ℓ = 80, and the spectral index of the spectrum (C_ℓ), respectively.

We quantify the quality of the fits by computing the reduced χ² over the same multipole range as $${\chi }_{XX}^2\equiv \frac{1}{N_{\text{bin}}}{\sum }_{{\mathcal{l}}_b}\frac{{\left({\mathcal{D}}_{{\mathcal{l}}_b}^{XX}-{\widehat{\mathcal{D}}}_{{\mathcal{l}}_b}^{XX}\right)}^2}{{\sigma }_{{\mathcal{D}}_{{\mathcal{l}}_b}^{XX}}{}^2},$$(11)

in which the uncertainties are from Eq. (9) and N_bin is the number of multipole bins involved.

Fig. 3 Examples of polarization (pseudo) power spectra (${\mathcal{D}}_{\mathcal{l}}$) obtained for the polarization maps shown in Fig. 2 and for the six values of the retained sky fraction (fsky) in the analysis; from 0.2 to 0.7 (top to bottom). The gray-shaded areas are excluded from our analysis as explained in the text. The straight lines are power-law power spectra fitted in the white area, i.e. the range 60 ≤ ℓ ≤140.

3.3.2 Choice of ℓ range

In order to compare polarization power spectra from different sets of maps and masks which may look very different, it is necessary todefine a common range of multipoles within which to characterize the power spectra.

By visually inspecting the power spectra, we notice that most of them show power-law behavior up to large multipole bins (ℓ ~ 360), especially when the maps contain a significant fraction of the sky. However, some power spectra show a sharp drop of power or a flattening in the high-ℓ range of the spectrum. These features either originate from the limited resolution of the 3D AMR grid (e.g., in situations where only a few low-resolution voxels contribute to the simulated unmasked sky) or from sharp edges in the maps created by the projection. We find that cutting the ℓ > 140 mitigates both effects homogeneously over our entire sample of power spectra.

At low-ℓ values, our reliability test discussed below indicates that it is safer to discard the bin centered at ℓ = 50 because some masks may lead to biased estimates of the power spectra. Therefore, we decide to constrain the characterization of all the power spectra in the multipole range 60 ≤ ℓ ≤ 140. This low-ℓ cut mitigates the spurious effects of the projected cubic voxels (visible in some maps) on the power spectrum estimate, because it avoids the range of angular scales corresponding to the projected voxels at high and intermediate latitudes.

For this multipole range, the EE and BB power spectra can be described by power laws. Out of the 2592 sets of power spectra, about 82% show reduced ${\chi }_{XX}^2$ lower than 10, using Eq. (9) as uncertainties, for both EE and BB simultaneously.

3.3.3 Reliability of the power spectra

In order to check the estimates of the uncertainties, but more importantly, to make sure that the masks used do not introduce biases in our polarization power spectra, we rely on a simulation-based validation for a subset of observer locations.

For each galaxy model we produce 12 sets of polarization maps from positions 30° apart. We apply the six masks to each set and calculate the power spectra, like above. We fit the resulting TT, EE, and BB auto-power spectra with power-laws in the multipole range 40 ≤ ℓ ≤ 140 (see below). Finally, we compute the three correlation coefficients between the spectra using Eq. (8). The fitted power-law auto-power spectra and the correlation coefficients determine a set of (idealized) polarization power spectra. Then, using the synfast functionality of the HEALPix Python package, we use the idealized power spectra to create 300 sets of full-sky polarization maps through Gaussian (Monte Carlo) realizations. We then apply the sky-fraction masks to these simulated maps and compute their power spectra with Xpol. To detect and quantify possible biases, we compare the value of the idealized spectra to the distribution (mean and standard deviation) of the Monte Carlo (MC) power spectra at each multipole bin.

For all the tested cases, we find that the MC and analytic uncertainties generally agree in the range 40 ≤ ℓ ≤ 140. Therefore, to avoid the additional computational cost of running MC simulations, in what follows we only consider the analytic uncertainties. We then verify that our power spectra do not suffer from biases. We find that several polarization maps lead to a significant bias for the multipole bin centered at ℓ = 50. However, no significant bias is detected in the range 60 ≤ ℓ ≤ 140. Therefore, we focus the comparison on the multipole range 60 ≤ ℓ ≤ 140.

The polarization power spectra corresponding to the maps shown in Fig. 2 are presented in Fig. 3 for the six f_sky values. On those plots we are also presenting results from fitted power-law power spectra on TT, EE, BB and TE in the range of 60 ≤ ℓ ≤ 140.

Fig. 4 Histograms of REB and rTE as measured from all galaxy setups, observer positions and sky fraction together. The red cross marks the location $(r^{TE},{\mathcal{R}}_{EB})=(0.4,2.0)$, typical values deduced from Planck data at 353 GHz.

4 Results

In this section we first present the characteristics of the full sample of power spectra. Then we examine the differences arising fromthe underlying galaxy model. We characterize the properties of the power spectra through the ${\mathcal{R}}_{EB}$ ratio, the correlation coefficient r^TE and power-law fits of the EE and BB spectra. We quantify the variance induced by the observer’s location and since we find it to be significant, we then explore our sample and simulations in more detail in the attempt to identify the dominant factors that determine the shape of the power spectra and their correlations.

4.1 Full sample inspection

We use the parameters $A_{80}^{XX}$, α_XX, and ${\chi }_{XX}^2$ in addition to ${\mathcal{R}}_{EB}$ and r^TE to characterize the synthetic polarization power spectra. Figure 4 presents histograms of ${\mathcal{R}}_{EB}$ and r^TE as measured from all galaxy setups, observer positions and sky fractions.

The ${\log }_{10}({\mathcal{R}}_{EB})$ distribution is slightly shifted toward positive values, which means that E modes dominate most of the polarization signal in our synthetic maps. However, there is a large scatter about the mean and several sets of polarization maps are dominated by B modes. The 16, 50, and 84 percentiles of ${\log }_{10}({\mathcal{R}}_{EB})$ are − 0.17, 0.15, and 0.43 respectively (the median of ${\mathcal{R}}_{EB}$ is 1.42). The correlation coefficient r^TE is mostly positive, with a second minor peak at negative values, revealing a general positive correlation between T and E modes. Interestingly, the pair of typical values reported for the real sky by Planck (marked by the red cross in Fig. 4) is covered by the 2D distribution obtained from the synthetic maps.

Figure 5 (panel a) presents 2D histograms for some combinations of parameters from the power-law fits as measured in the multipole range 60 ≤ ℓ ≤ 140 from all galaxy setups, observer positions and sky fractions.

The EE and BB auto-power spectra are very steep (α_XX ≈−4.5); much steeper than those observed in the Planck data (α_XX ≈−2.5) or in kpc-scaled MHD simulations (Kim et al. 2019) (α_XX ≈−3.6). The reason for this behavior is that most of the high-resolution areas of the AMR grid are excluded from the analysis by masking. This prevents our synthetic skies from being populated by small angular-scale variations. The limited resolution of the retained sky also reduces the variance in the maps, and is therefore responsible for the low values of the $A_{80}^{XX}$ parameters compared to the real sky.

As illustrated in Fig. 5, the amplitudes at ℓ = 80 of the EE and BB auto-power spectra span several orders of magnitude but appear to be correlated. Generally, $A_{80}^{EE}$ is found to be larger than $A_{80}^{BB}$. The distribution of ${\log }_{10}(A_{80}^{EE}/A_{80}^{BB})$ is roughly symmetric with percentile 16, 50, and 86 of − 0.12, 0.16, and 0.42 (the median of $A_{80}^{EE}/A_{80}^{BB}$ is 1.45). The spectral indices, α_XX, of the EE and BB power spectra are also positively correlated and have similar values for the entire set of maps. The two bottom panels of the figure show that very steep power-law fits to the spectra have large ${\chi }_{XX}^2$ values. These poorly-fitted spectra with a lack of power at small scales likely correspond to quiescent sky realizations.

Figure 5 (panel b) shows that ${\log }_{10}({\mathcal{R}}_{EB})$ and ${\log }_{10}(A_{80}^{EE}/A_{80}^{BB})$ are very well correlated. This correlation implies that, despite the small differences in the spectral indices of the EE and BB auto-spectra, the ${\mathcal{D}}_{\mathcal{l}}^{EE}/{\mathcal{D}}_{\mathcal{l}}^{BB}$ ratios are not strongly scale-dependent. Therefore, the ${\mathcal{R}}_{EB}$ parameter can be used to derive meaningful information.

4.2 Dependence of ${\mathcal{R}}_{\text{EB}}$ and r^TE on the galaxy model

In this section we infer a possible dependence of the ${\mathcal{R}}_{EB}$ and r^TE parameters on the simulation setup and on the retained sky fraction. We find that different galaxy models lead to different distributions of ${\mathcal{R}}_{EB}$ and r^TE and a small but systematic increase of ${\mathcal{R}}_{EB}$ and r^TE values with increasing sky fraction.

Figure 6 summarizes the distributions of ${\log }_{10}({\mathcal{R}}_{EB})$ (left) and r^TE (right) per galaxy (Pwf, Pwnf, Twf, Twnf, Tsf and Tsnf, from top to bottom panel respectively) as a function of retained sky fraction (f_sky). The centralpoint is the median and the asymmetric error bars show the 1σ spread of the 72 data points we have per galaxy and f_sky.

4.2.1 The ${\mathcal{R}}_{\text{EB}}$ parameter

For all six galaxy setups, ${\log }_{10}({\mathcal{R}}_{EB})$ is found to be mainly positive with medians in the range from ${\mathcal{R}}_{EB}~1$ to ${\mathcal{R}}_{EB}~2$. Although the spread of the distributions are larger for low f_sky values, a striking difference between galaxies with and without feedback appears: Galaxies for which the feedback is disabled show a much larger dispersion of the ${\mathcal{R}}_{EB}$ ratio.

Figure 6 shows a small but systematic increase of ${\log }_{10}({\mathcal{R}}_{EB})$ with f_sky. This trend reflects the fact that, as f_sky increases, more of the bright filamentary patterns clearly seen above and below the disk in Fig. 2 are included in the spectra. These features are generally well-aligned with the magnetic field (see Fig. 7 where the integrated and projected magnetic field orientation is superimposed over the intensity map close to one of these features), and therefore produce more E-modes than B-modes (e.g. Zaldarriaga 2001; Clark et al. 2021; Konstantinou et al. 2021).

The increase of ${\log }_{10}({\mathcal{R}}_{EB})$ with f_sky, however, is not significant with respect to the variance induced by the observer’s location for each f_sky, shown as error bars in the same figure. Therefore, for each galaxy model, we can group measurements from different f_sky values together to generate the histograms shown in Fig. 8. The large difference of ${\log }_{10}({\mathcal{R}}_{EB})$ distributions between galaxies that include feedback (top panel) from those that do not (bottom panel) is clearly demonstrated. The presence of feedback reduces the range of ${\mathcal{R}}_{EB}$ values, inducing a peak between $1\lesssim {\mathcal{R}}_{EB}\lesssim 2$. Table 2 contains summarizing statistics of these distributions.

We quantify the similarity of ${\mathcal{R}}_{EB}$ distributions obtained from the different galaxies using a 2-sample Kolmogorov- Smirnov (KS2S) test, computing the probabilities that distributions in pairs are drawn from the same parent distributions. These probabilitiesare reported in the lower left triangle of Table 3.

The distributions of ${\log }_{10}({\mathcal{R}}_{EB})$ in model Pwf and Twf are very similar to each other. ${\mathcal{R}}_{EB}$ does not seem to be related to the initial topology of the magnetic field, at least when the field strength is low (models Twf, Twnf, Pwf, Pwnfb). Comparison of the ${\log }_{10}({\mathcal{R}}_{EB})$ distributions of Twf and Tsf reveals a possible effect of the magnetic field strength. The stronger the field, the narrower the distribution and the higher the central value.

The ${\mathcal{R}}_{EB}$ distributions are more difficult to compare between the galaxies without feedback. However, the ${\log }_{10}({\mathcal{R}}_{EB})$ distributions of Twnf and Tsnf appear more similar to one another than to that of Pwnf, which extends to larger values.

Fig. 5 2D histograms of combinations of parameters from power-law fits used to characterize the shape of the power spectra for the full sample; all galaxy setups, observer position and sky fractions merged together. Panel a: from top left to bottom right we have $({\log }_{10}(A_{80}^{EE}),{\log }_{10}(A_{80}^{BB}))$, (αEE, αBB), $({\alpha }^{EE},{\chi }_{EE}^2)$ and $({\alpha }^{BB},{\chi }_{BB}^2)$. The main diagonals are shown in top panels for reference. In the bottom panels the horizontal lines mark values of ${\chi }_{XX}^2=1$ and the vertical lines show the median of the spectral indices. Panel b: ${\log }_{10}({\mathcal{R}}_{EB})$ versus ${\log }_{10}(A_{80}^{EE}/A_{80}^{BB})$. The horizontal and vertical lines mark ratio values of 1 and 2 and are given for visual reference. The main diagonal is also shown.

4.2.2 The r^TE parameter

Figure 6 (right) shows that the synthetic polarization maps from all the simulations except Pwnf have mainly positive r^TE values, independently of the f_sky. For these five models, r^TE increases with f_sky. This trend stems from the inclusion of more structures as f_sky increases. Nevertheless, the r^TE values at f_sky = 0.2 and 0.7 agree within uncertainties and, as for ${\mathcal{R}}_{EB}$ we can combine measurements obtained with different f_sky values to compare characteristics of the power spectra between different galaxies. The right panel of Fig. 8 shows the histograms of r^TE for different galaxies. Summarizing statistics are provided in Table 2. We quantify the similarity of the r^TE distributions obtained from the different galaxies using a KS2S test, computing the probabilities that distributions in pairs are drawn from the same parent distributions. Those probabilities are reported in the upper right triangle of Table 3.

In general (as seen from Fig. 8, right panels), the r^TE histograms indicate that a toroidal magnetic field topology (i.e. field lines parallel to the disk) leads to positive TE correlations, and that a stronger magnetic field strengthens the correlation. As for the ${\mathcal{R}}_{EB}$ ratio, the presence of feedback in the simulations appears to shrink the distribution of possible r^TE values. However, here the effect is milder.

The r^TE values from the Pwnf galaxy are consistently negative for all f_sky values and present no significant trend with f_sky. This anticorrelation of T and E modes is caused by the peculiar magnetic field topology of this model. In Pwnf, the field lines are by construction perpendicular to the disk, and they remain so because of the absence of feedback and the short integration time of the simulation. As Konstantinou et al. (2021) demonstrate, a configuration where density structures and magnetic field are perpendicular to each other leads to a negative r^TE. In the cases without feedback and a toroidal field geometry, the field lines remain parallel to the density structures of the ISM, a configuration that leads to positive r^TE. Interestingly, the rare occurrences with positive r^TE in model Pwnf correspond to observers located at the edges of the spiral arms. In these cases magnetic field and density structures locally follow coherent orientations and cover a large area of the sky. In other cases, the polarization sky is very quiescent.

Fig. 6 Summarizing statistics of ${\mathcal{R}}_{EB}$ (left) and rTE (right) per galaxy and per sky fraction fsky. Asymmetric error bars around the median encompass 68% of the data points. For visual reference, the horizontal dashed lines mark the values of ${\mathcal{R}}_{EB}=1$ and 2 to the left and rTE = 0 and 0.4 to the right.

4.3 The viewpoint-induced variance

Both ${\mathcal{R}}_{EB}$ and r^TE show an increase with f_sky. As f_sky increases, lines of sight with higher column density are included in the power spectrum computation. The positive correlations of ${\mathcal{R}}_{EB}$ and r^TE with f_sky may thus reflect a correlation with density and or line-of-sight complexity of the probed ISM. The probed ISM density can be quantified through the weighted average of the column density or the weighted average of the emission intensity (Ī), where the weights come from the mask.

For a given simulated galaxy the weighted average of emission intensity depends on the observer’s location and on the retained sky fraction in the analysis. Figure 9 illustrates this complex dependence. The dependence of the mean intensity, ${\mathcal{R}}_{EB}$ and r^TE with azimuth angle ϕ_⊙ (the angular coordinate of the observer in the galactic disk as introduced in Appendix A) is shown on the left, middle, and right panels, respectively. The difference between galaxies with and without feedback is again striking: the presence of feedback creates much more constrained distributions of ${\mathcal{R}}_{EB}$ and r^TE. However, based on these figures it is difficult to conclude whether or not the ${\mathcal{R}}_{EB}$ and r^TE azimuthal modulations are correlated with those of Ī.

To address this question we carry out Spearman rank order correlation tests to quantify whether or not ${\mathcal{R}}_{EB}$ and or r^TE are correlated with Ī either when all galaxies and f_sky are considered together or when galaxies are considered separately. The results are reported in Table 4.

The correlation of ${\mathcal{R}}_{EB}$ with Ī is strong whenmeasurements from all sky fraction values and all galaxies are considered together. There is a relatively strong correlation between ${\mathcal{R}}_{EB}$ and Ī for Twf andTwnf, and a weak correlation for Pwf and Pwnf. There is no correlation between ${\mathcal{R}}_{EB}$ and Ī for Tsf and Tsnf. The absence of correlation for Tsf is likely due to the smaller range of ${\mathcal{R}}_{EB}$ values encountered in this galaxy.

The correlation of r^TE with Ī is strong and significant when gathering all sky fraction and all galaxies into a unique sample. It also persists when examining each model individually. The strongest correlations are found for galaxies with feedback and low value of the initial magnetic field strength.

In general, ${\mathcal{R}}_{EB}$ and r^TE are well correlated with Ī for the Twf galaxy whereas Tsf (an identical model with stronger field) shows almost no correlation. Again, the loss of correlation is likely due to the smaller range of ${\mathcal{R}}_{EB}$ and r^TE observed forTsf as compared to Twf.

We conclude that ${\mathcal{R}}_{EB}$ and r^TE depend strongly on the observer’s position. This dependence is generally reflected by a correlation of both ${\mathcal{R}}_{EB}$ and r^TE with the weighted averaged of the emission intensity.

Fig. 7 Zoom in the set of maps presented in Fig. 2 to illustrate the alignment between the projected plane-of-the-sky magnetic field orientation and projected density structure orientation. The background color shows the intensity map while the streamlines show the magnetic field orientation as inferred from the Q and U Stokes parameters.

4.4 The observer environment

Given that the summarizing characteristics of the polarization power spectra appear to depend sensitively on the observer’s location, we want to check whether placing the observer in similar environments yields similar ${\mathcal{R}}_{EB}$ and r^TE values independently of the galaxy model. In the following, we start by investigating possible differences in the power spectra obtained by observers located within or outside spiral arms (Sect. 4.4.1). Since we do find substantial differences between these two subgroups, we proceed to search for any correlations of the power spectrum characteristics with the mean dust density and magnetic field strength in the observer’s vicinity (Sect. 4.4.2). We cannot find any, but caution against the fact that, due to masking, those local estimates are not fully representative of the ISM regions that are imprinted in the power spectra. Finally, we also investigate the effect of placing the observer within supernova-driven bubbles (Sect. 4.4.3). We find that power spectra taken from within bubbles are peculiar within each galaxy model, even if the scatter remains large. However, the effect of a bubble environment generally supersedes the imprint of the underlying galaxy model.

4.4.1 Arm vs. inter-arm regions

We first want to investigate the possible effects of having the observer in an arm or in an inter-arm region. Figure 9 (top right) shows the normalized angular profile (ñ_d(ϕ_⊙)) of the mean density ${\overline{n}}_{\text{d}}$ in a sphere of given radius (R_sph) surroundingthe observer which we normalize by the median of the profile. The location of the two arms are well spotted by eye. Within the arms they are several substructures with local maxima and minima, possibly due to the presence of large aggregates of clumps and large under-densities or bubbles. To automate the determination of the angular coordinates of the edges of the arms at the galacto-centric radius of 8 kpc, we consider the mean-density curves (ñ_d(ϕ_⊙)) obtained with R_sph = 200 pc measured at an angular step of 1°. We smooth those curves with a wide Gaussian kernel (FWHM of 10°) to eliminate small-scale features and consider the successive derivatives of ñ_d (ϕ_⊙). The edges of the arms are defined as the locii where the third derivatives (${\partial }_{{\varphi }_{\odot }}^3({\tilde{n}}_{\text{d}})$) vanish on either sides of the two main maxima of ñ_d(ϕ_⊙). For each galaxy, the range of angular coordinates span by the arms are marked by gray bands on the panels of Fig. 9.

It is interesting to note that on the right side of most of the arms one can observe sharp dips that reveal the trails of the density waves that propagate through the disk toward decreasing ϕ_⊙. This feature is best seen with R_sph = 200 pc (and with angular sampling of 1°, not shown here).

In order to study whether or not the fact that the observer is in an arm influences the characteristics of the polarization power spectra, we create subsamples of measurements corresponding (i) to observers within the arms (i.e. with ϕ_⊙ falling in the gray bands of Fig. 9) and (ii) to observers located away from the arms. The latter are obtained by shifting the limits of the arms by 90°. This ensures that, for a given galaxy, the sub-samples of spectra taken within and away from the arms have the same size. The number of observer locations thus obtained ranges from 15 to 17 per model. We visually check (on column density maps like the one in Fig. A.1) that our arm/inter-arm determination is effective. In Table 5 we report descriptive statistics of ${\mathcal{R}}_{EB}$ and r^TE distributions obtained for the six galaxies, merging measurements obtained with all f_sky.

For most of the cases, we see that both ${\mathcal{R}}_{EB}$ and r^TE are larger when observed from within the arms. Only models Pwnf and Tsf do not show this trend.

We test the null hypothesis that distributions of ${\mathcal{R}}_{EB}$ (r^TE) as measured from inside or away from the arms are drawn from the same parent distribution using a KS2S test. The results are reported in Table 6. Generally, the arm and inter-arm distributions differ. This effect is stronger for galaxies without feedback and is not observed for the Tsf galaxy. The latter is consistent with the fact that the variance of both ${\mathcal{R}}_{EB}$ and r^TE distributions are smaller for the Tsf galaxy. The distributions of ${\mathcal{R}}_{EB}$ and r^TE within and away from the arms still differ according to the underlying galaxy simulations.

According to the latest observational evidence, the Milky Way probably has four spiral arms (e.g., Reid et al. 2019), while our input models only include two. Doubling the number of arms would certainly affect our results. According to our findings, the general trend is for both ${\mathcal{R}}_{EB}$ and r^TE to increase for maps obtained from within the arms compared to other locations. Therefore, we may speculate that we would observe more frequently large values of ${\mathcal{R}}_{EB}$ and r^TE in histograms like the ones in Fig. 8 for a galaxy with four instead of two spiral arms. However, the difference seen between arm and inter-arm regions should persist. In future work, this conjecture needs to be tested against new simulations because other non-trivial effects could impact the power spectra.

Fig. 8 Histograms of ${\log }_{10}({\mathcal{R}}_{EB})$ (left) and rTE (right) as measured on polarization maps synthesized from the three models with feedback (top) and without feedback (bottom). Power spectra corresponding to all values of fsky are put together, making an ensemble of 432 measurements per galaxy setup. The vertical gray dot-dashed (red dashed) lines mark the values of ${\mathcal{R}}_{EB}=1$ (2) on the left and of rTE = 0 (0.4) on the right, for reference.

Fig. 9 Variations of the averaged weighted emission intensity (Ī) (top left), ${\mathcal{R}}_{EB}$ (bottom left) and rTE (bottom right) as a function of the azimuthal position of the observer (ϕ⊙ [°]) for the different galaxies (indicated in the captions) and for different values of sky fraction (fsky, indicated by the colors in the legend). Top right: evolution of the averaged dust density (normalized to the sample median) in spheres of radii (100 pc, 200 pc and 500 pc) around an observer as a function its azimuthal coordinate (ϕ⊙). The colors represent different values of the radius of the sphere surrounding the observer as informed by the legend. Spheres are independent if their ϕ⊙ differ by at least 1.43°, 2.86° and 7.16° for radius Rsph of 100 pc, 200 pc, and 500 pc, respectively. The vertical gray bands indicate the range of ϕ⊙ values in which we consider the observer to be located in the arms (see Sect. 4.4.1).

4.4.2 Local density and magnetic field strength

In this section we search for correlation of ${\mathcal{R}}_{EB}$ and r^TE with the local density and the local strength of the magnetic field. We determine the local values for each observer location in measuring the mean values within a sphere of a given radius (R_sph) around the observer. We consider two values for the radius: 200 and 500 pc. These values are chosen such that the spheres encompass most of the dusty ISM matter that contribute to the simulated polarization skies at high galactic latitudes. However, because the brightest regions are masked, these local values may not be fully representative of the part of the ISM that is included in the power spectrum estimates.

The ranges of local density values encountered in all the six galaxy models are similar and follow the variation about the coordinate ϕ_⊙ shown in Fig. 9 (top right). We are not able to find convincing correlation of ${\mathcal{R}}_{EB}$ or r^TE with the local density. Only the no-feedback galaxies do show weak correlations (with p-values at the order of 1%) between r^TE and the local density. This is because both r^TE and the density increase within the arms (see also previous section) as it can be inferred by inspecting the bottom right panel of Fig. 9.

The range of values of the strength of the magnetic field that is local to the observer ($\text{left}<\left|B\right|\text{ right}{>}_{R_{\text{sph}}}$) encountered in all the six galaxy models spans several orders of magnitude. We concentrate on the Twf and Tsf galaxies in order to infer the effect of the strength of the magnetic field while keeping as much similar as possible the other parameters. In this case, $\text{left}<\left|B\right|\text{ right}{>}_{R_{\text{sph}}}$ spans four orders of magnitude, as shown in Fig. 10. Despite this large coverage we find no correlation of ${\mathcal{R}}_{EB}$ or r^TE with $\text{left}<\left|B\right|\text{ right}{>}_{R_{\text{sph}}}$ as also illustrated in that figure and as confirmed through Spearman rank order tests that we performed for all values of the sky fraction and for both R_sph = 200 and 500 pc. For R_sph = 500, the probability that any observed correlation happens by chance are 55 and 12% for the pairs $(\text{left}<|B|\text{ right}{>}_{R_{\text{sph}}}\text{},{\mathcal{R}}_{EB})$ and $(\text{left}<|B|\text{ right}{>}_{R_{\text{sph}}}\text{},r^{TE})$, respectively. Therefore, based on this analysis we cannot conclude that the strength of the local field to the observer is an important factor in shaping the polarization power spectra. We emphasize that this estimate of $\text{left}<\left|B\right|\text{ right}{>}_{R_{\text{sph}}}$ does not necessarily represent adequately the strength of the field permeating the medium that is imprinted in the polarization maps. The test carried in Sect. 4.5 addresses this point in particular.

Fig. 10 Scatter plots of strength of the local magnetic field (in μG) versus ${\mathcal{R}}_{EB}$ (left) and rTE, measured from the Twf (green) and Tsf (purple) galaxies for fsky = 0.7. The strength of the field is the volume-weighed average of the field of all cells having their centers within a sphere radius of 500 pc centered on the observer. Gray and red vertical lines are the same as in Fig. 8. The gray tine solid line show the medians of the data points and the shacked black lines show running means along |B | considering 10 data points at a time.

4.4.3 Within bubbles

There is ample observational evidence that the Sun resides in a special environment called the Local Bubble. The Local Bubble is a cavity of hot plasma created by supernova explosions, surrounded by a magnetized shell of cold, dusty gas (e.g., Lallement et al. 2018; Pelgrims et al. 2020). Guided by this observational fact, we create maps as before, this time placing observers inside supernova-driven bubbles, self-consistently generated in the simulations. We visually identify 14, 14, and 11 bubbles within the Pwf, Twf, and Tsf galaxies in various locations (center, arm, and inter-arm regions) of the galactic midplane. We do not impose any other selection criteria on our bubble sample. Their sizes range from ~ 100 pc to ~1.3 kpc. Interestingly, we find that the strength of the magnetic field within the bubbles ranges from 0.5 to 1 μG independently of their environment (center, arm, inter-arm) and their underlying galaxy model. This coincidence may be related to the physical conditions necessary to lead to the formation and explosion of stars. However, a detailed study of this feature is outside the scope of this paper. In what follows we compare the distributions of ${\mathcal{R}}_{EB}$ and r^TE as obtained from the inside of bubbles (shown in Fig. 11) to those from the main samples (see top row of Fig. 8).

First, for each galaxy model, we investigate whether the ${\mathcal{R}}_{EB}$ and r^TE distributions from the bubbles can be considered random occurrences within the main samples. We use a KS2S test to address this question. The probabilities that the distributions corresponding to the bubble locations and those of the main sample are drawn from the same parent distribution – for the same galaxy model – are found to be low, as reported in Table 7. Therefore, it is unlikely that the distributions of ${\mathcal{R}}_{EB}$ and r^TE measured from the inside of bubbles are random realizations of the main sample. The polarization skies and power spectra, from inside bubbles are peculiar realizations in their respective galaxy: Both the ${\mathcal{R}}_{EB}$ and r^TE distributions from within bubbles are shifted toward lower values. The effect is more pronounced for Twf and Tsf.

Second, we study whether the distributions corresponding to the bubbles are statistically similar irrespective of their underlying galaxy. We compute the KS statistics (ks^⋆) and the corresponding probabilities ($p_{KS}^{\star }$) for each comparison in pairs. The values are reported in the upper blocks of Table 8. We see that the distributions are generally more similar to each other than the main samples are to one another (see Table 3). The distributions of ${\mathcal{R}}_{EB}$-bubbles and r^TE -bubbles from Twf and Tsf are those that keep deviating the most with a minimum probability for ${\mathcal{R}}_{EB}$ of ≲ 1% of being drawnfrom the same parent distribution.

To assess whether this similarity between the distributions is not simply an effect of decreased sample size we proceed as follows. We generate 10.000 bootstrap sub-samples from the main sample, each with the size of the bubble samples. The bootstrap selection is performed by means of observer location. Then, we compare bootstrapped samples from two different galaxy setups and for each compute the KS2S statistics ks. This leads to distributions of ks for the comparison of two sub-samples drawn from, say, Pwf and Twf. We then compare the same statistics but between Pwf-bubble and Twf-bubble (ks^⋆) to the distribution obtained above in computing the number of draws with ks ≤ ks^⋆ and divide by the total number of draws. This gives the one-sided p-value that the bubbles samples look similar by chance, i.e. that random sub-samples from the main samples lead to at least this level of similarity due to the loss of statistic. These probabilities are reportedin the bottom block of Table 8.

We see that the probability of randomly observing a ${\mathcal{R}}_{EB}$ from the main sample that is similar to one from the bubble sample is high. However, that is not true of the r^TE distributions. There, the probability reaches a minimum ≲1% for the comparison of Pwf-bubble and Tsf-bubble.

This observation is consistent with our earlier finding that, for the toroidal field geometry, taking power spectra from within a bubble reduces the r^TE values. This effect is more pronounced the stronger the magnetic field and negligible for the poloidal field topology, where the observer’s position induces a larger scatter.

In conclusion, we observe that given an underlying galaxy model, the polarization power spectra taken from within bubbles are peculiar. The distribution of ${\mathcal{R}}_{EB}$ and r^TE taken from within bubbles are generally similar, independently of the galaxy setup, but the scatters remain large. Furthermore, the influence of the underlying galaxy model seems small. More specifically, observing from inside a bubble seems to decouple the observed sky from the underlying galaxy. The information about the topology of the magnetic field at large scales is lost (comparing Pwf-bubbles and Twf-bubbles) while a mild effect from the strength of the field may survive (comparing Twf-bubbles and Tsf-bubbles). However, we find no correlation between ${\mathcal{R}}_{EB}$ or r^TE with the (local) field strength in bubbles, the size of the bubble, or with the environment of the bubbles within the galaxy. Additionally, for the bubble samples, we do not observe a correlation between the averaged weighted intensity and the ${\mathcal{R}}_{EB}$ and r^TE values. Combining all bubbles and all values of f_sky in one sample, the Spearman rank order coefficients are 0.04 and 0.16, to which p-values of 56 and 1.2% correspond, respectively.

Fig. 11 Normalized histograms of ${\log }_{10}({\mathcal{R}}_{EB})$ (left) and rTE (right) as measured on polarization maps synthesized from the inside of bubbles detected in the three models with feedback. Power spectra corresponding to all values of fsky are put together. Vertical lines are the same as in Fig. 8.

4.5 Strength of B in the cone of observation

An important result from the previous subsections is that the distributions of ${\mathcal{R}}_{EB}$ and r^TE for Twf and Tsf differ significantly. This suggests that the overall strength of the magnetic field in the galaxy may play a role in shaping the polarization power spectra. This observation holds true for the spectra taken from within bubbles, even though we have shown that the bubble environment has the general effect of regularizing the distributions.

Here, we quantify the dependence of the ${\mathcal{R}}_{EB}$ and r^TE on the magnetic field strength in the retained regions of the ISM. For this purpose, for each observer position, we consider a double (two hemispheres) vertical cone with an apex angle of 120°, a height of 1 kpc, and with the observer at the tip (the polar caps defined by the intersection of such cone and the celestial sphere cover half of the sky). Then, we calculate the mass-weighed average strength of the field using all the simulation cells locatedin the cone. This process is repeated for all observer positions in the Twf and Tsf galaxies and for f_sky = 0.5³ Fig. 12 shows ${\mathcal{R}}_{EB}$ and r^TE as a function of these magnetic field strength estimates for each observer. Although the field strength measurements span several orders of magnitude, we do not find any correlation with ${\mathcal{R}}_{EB}$ or r^TE (a Spearman correlation test shows that the probability that these scatter plots are random are 5% and 3% respectively). Therefore, we conclude that the magnetic field strength in the retained ISM regions is not directly responsible for the differences of distributions seen in Figs. 8 and 11.

Fig. 12 Scatter plots of strength of the magnetic field (in μG) in the cone of observation versus ${\mathcal{R}}_{EB}$ (left) and rTE (right), measured from the Twf (green) and Tsf (purple) galaxies for fsky = 0.5. The strength of the field is the mass-weighed average of the field for all cells having their centers in the cone centered on the observer, with an opening angle of 60° around the z axis and with height lower than 1 kpc. Vertical and shacked lines are built the same way as in Fig. 10.

5 Summary and conclusion

In this work, we synthesized full-sky observations of polarized dust emission at 353 GHz from a set of global, Milky-Way-sized galaxy simulations. We studied the resulting T, E, and B power spectra in terms of the E∕B power asymmetry (${\mathcal{R}}_{EB}$) and the correlation between T and E modes (r^TE), aiming to identify the physical parameters that affect these correlations. Specifically, we investigated the dependence of ${\mathcal{R}}_{EB}$ and r^TE on (i) the physical conditions in the galaxy simulations (presence of feedback, topology, and strength of the large-scale magnetic field) and (ii) the particular location of the observer in the galaxy. Our intention is not to reproduce or model actual observationsof the Galactic magnetized ISM but rather to identify the dominant factors shaping the polarization power spectra.

We found that the distributions of ${\mathcal{R}}_{EB}$ and r^TE from different galaxy models differ significantly: The presence of feedback, the topology, and the strength of the magnetic field at large scales all impact the shape of the distributions. However, the cosmic variance dominates: the ${\mathcal{R}}_{EB}$ and r^TE parameters measured from polarization power spectra depend sensitively on the observer’s location.

In particular, we showed that measuring the ${\mathcal{R}}_{EB}$ and r^TE distributions from within a spiral arm as opposed to an inter-arm region leads to a measurable difference. However, we did not find any correlation between the local matter density or the strength of the magnetic field and the measured ${\mathcal{R}}_{EB}$ and r^TE values. We further demonstrated that taking power spectra from the inside of supernova-driven superbubbles minimizes the dependence on the underlying galaxy model. However, it does not eliminate it: the distributions of ${\mathcal{R}}_{EB}$ and r^TE from within superbubbles still differ between models.

We found evidence that the strength of the global magnetic field in the galaxy model (which is below equipartition values in all of the models) may play a role in shaping the distribution of power spectrum characteristics; it decreases the viewpoint-induced variances and increases both the E∕B asymmetry and the TE correlation. However, we found no correlation between the power spectra characteristics and the strength of the magnetic field in the regions of the ISM that enter the maps we analyzed.

In summary, we have shown that the statistics of the sky in polarization depend predominantly on the observer’s point of view. However, no correlation with the local mean values of physical quantities such as density or magnetic field strength was found. Therefore, any statistical properties of the polarized sky from intermediate to high-Galactic latitude regions may reflect the complexity of the line-of-sight structure of the magnetized ISM rather than mean local or global properties of the ISM. This result reinforces previous works that demonstrate the importance of local structures in the study and characterization of the observed sky in polarization (e.g., see Alves et al. 2018; Bracco et al. 2019a; Skalidis & Pelgrims 2019; Pelgrims et al. 2020). This evidence suggests that modeling the dynamical history of the Solar neighborhood and its 3D structure may be necessary to obtain a better description of the polarization sky at intermediate and high-Galactic latitudes. However, restricting our study to the local Solar neighborhood may not be sufficient: the statistical properties of the polarization sky also depend on more distant regions of the ISM probed by the lines of sight. This motivates modeling efforts to connect models of the local ISM to large-scale models of the magnetized Galaxy.

Acknowledgements

The authors thank Andrea Bracco, Fabio Del Sordo, and Vincent Guillet for helpful comments and discussions and their anonymous referee for providing a constructive and thorough report. This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme under grant agreements No. 771282 and 740120. E.N. acknowledges funding from a Marie Curie Action of the European Union (grant agreement No. 749073).

Appendix A Observer excursion: Angular convention

In this work we consider an hypothetical observer on the xy plane of the Cartesian grid of the simulation (z = 0) at a distance of R_⊙ = 8 kpc from the galaxy center. To quantify the effect of the observer’s location on our inferences of the dust polarization maps, we move the observer in a circle on the xy plane; as illustrated in Fig. A.1 on a top-down view of the Tsf galaxy model. We choose the galactic center to be placed at longitude zero (l = 0°) on the synthetic maps for any observer location. To achieve this, we rotate the AMR grid according to $x^{\prime }={\mathcal{R}}_{{\varphi }_{\odot }}{x}_{\text{AMR}}$ where ${\mathcal{R}}_{{\varphi }_{\odot }}$ is the clockwise rotation matrix by an angle ϕ_⊙ about the z axis. Then we translate the full grid to obtain the coordinates of each cell in the observer Cartesian coordinate system (x = x′ + R_⊙ (1, 0, 0)) where the conversion to spherical coordinates is straightforward. The spherical coordinate system centered on the observer is the natural reference frame in which we want to synthesize the observables. In the observer reference frame the galactic center is always at coordinates (x, y, z) = (8.0, 0.0, 0.0) kpc. The used conventions and coordinate transformation are sketched in Fig. A.2. We apply the same vector transformation to the magnetic vector field ($B(x^{\prime })={\mathcal{R}}_{{\varphi }_{\odot }}{B}_{\text{AMR}}(x_{\text{AMR}})$) before we convert it to spherical coordinates, thus, centered on the observer: $$B_r=B\cdot e_r;B_{\theta }=B\cdot e_{\theta };B_{\varphi }=B\cdot e_{\varphi }$$(A.1)

where (e_r, e_θ, e_ϕ) is the orthonormal vector basis of an observer-centered spherical coordinate system with e_r pointing away from the observer and e_θ pointing southward. We discuss the synthetic maps in terms of galactic latitudes and galactic longitudes which are related to the spherical angular coordinates as (b, l) = (90 − θ, ϕ) when θ and ϕ are expressed in degrees.

Fig. A.1 Illustration of the observer’s circular excursion around the galaxy center on a top-down view of the Tsf galaxy model considered in the paper. A black star is placed every 5° at 8 kpc from the galaxy center. The background color trace the gas density.

Fig. A.2 Convention for observer angular position and AMR grid rotation. The third dimension (z axis) is dropped as it is left unchanged during our coordinate transformation. The observer is at (x, y) = (0, 0) and the galaxy center at (xAMR, yAMR) = (0, 0). The AMR grid is rotated clockwise by ϕ⊙ about the galaxy center and then shifted at (x, y) = (8, 0) to make the circular excursion of the observer.

Appendix B Ray-tracing algorithm

We rely on vector formalism to compute the path length through a cubic cell of a line of sight. This choice is motivated first by the fact that we consider rotated AMR grids by arbitrary angles and second by the very large number of cells in our simulations, from about 16 × 10⁶ to 40 × 10⁶. Both call for an automated and homogeneous procedure that can be implemented in an efficient algorithm.

Fig. B.1 Computing the distance di from the observer to the intersecting point (Ii) of a sightlinewith the bounding plane i of a cubic cell.

We consider the equations of the six planes bounding the cell and the equation of the line of sight with its unit direction vector, e_r, oriented outward from the observer. As sketched in Fig. B.1, the distance from the observer at the origin ($\mathcal{O}$) to the intersection (in I_i) between the line of sight with a plane can be written as $$d_i=\frac{P_i\cdot n_i}{e_r\cdot n_i}$$(B.1)

where i denotes one of the six planes and takes values from one to six, n_i is the unit vector normal to the plane and points outward from the cell center and P_i is one arbitrary point on the plane. We specify the latter from the cell center (c) and by a translation along the corresponding median of the cell: P_i = c + δ n_i, where δ is half the cell size.

To compute the path length through a cell, we encounter two cases. The observer is either fully outside the cell, including the faces, (distant cells) or it is inside (or on face) of the cell (observer cells). In both cases we compute the six distances d_i. We sort them by increasing order, defining the sorted set of distance d^[j] (j ∈ [1, 6]). For a distant cell of the simulation the path length through the cell is always given by λ = d^[4] − d^[3]. For the case of an observer cell, the intersection point at distance d^[3] is antipodal to the observer (toward −e_r) or is zero if the observer stands on a face. In those cases, the path length of the line of sight through the cell is simply given by d^[3].

It has to be noted that due to the AMR nature of the Cartesian grid that is ray-traced, the number of observer cells can be anything from one to eight depending on whether the observer is inside a cell or is placed on a face, an edge or even at a corner of the cells, and that adjacent cells may correspond to different resolution levels, so have different sizes. In such cases particular care is needed not to count more than once the path along cell faces since given lines of sight can be attributed to several cells. If not accounted for; this may produce artifacts at constant longitude(s) and or latitude(s) depending on the precise observer location within the grid and the ϕ_⊙ value which specify the grid orientation with respect to the celestial coordinate system of the observer.

In practice, the algorithm to ray-trace the simulation is defined as follows⁴. First, we fix the resolution of the map. Ideally it has to be sufficiently good to guarantee that all simulation cells are crossed at least once. The resolution fixes a set of pixels, i.e. of lines of sight. From this choice the algorithm then goes as follows:

For each simulation cell:

• 1.

  from cell position and size; select sightlines that may intersect the cell through the selection criteria: $\arccos \left(e_r\cdot c/\left|c\right|\right)\le 2\delta \sqrt{3}/\left|c\right|$.

• 2.

  for each selected lines of sight:

  • (a)

    compute the six distances d_i,

  • (b)

    sort the d_i’s and generate the list of d^[i],

  • (c)

    compute the path length through the cell λ = d^[4] − d^[3],

  • (d)

    make sure the intersection are on cell faces,

  • (e)

    use λ to update observable maps (e.g., Eqs. 4–6) or store λ and pixel index attached to the cell.

To check the validity of our ray-tracing algorithm, and thus of our synthetic maps, we generate the maps corresponding to the full length integrated up to the simulation box boundaries and corresponding only to the observer cells, for all observer view points and all simulation snapshots considered in this work. Examples of such maps are presented in Fig. B.2.

Fig. B.2 Integrated path length through the full simulation box (top panels) and through only the observer cells (lower panels). The full box is of 60 kpc on a sideand is integrated from an observer placed at R⊙ = 8 kpc from thecenter. The observer is immersed in the Tsf galaxy at 8 kpc from the galactic center and his angular coordinate is is ϕ⊙ = 180° for the left panels and ϕ⊙ = 60° for the rightpanels.

The path length from the observer cells is very small compared to the path length through the full box. However, if we do not include it in the sightline integration, the highly nonuniform features induced by the observer cells may produce strong patterns at the synthetic-map level when, for example, the disk-scale height right above and or right below the observer is small. Accounting for the path spent in the observer cell mitigates such artifacts.

Appendix C Map making and interpolation details

In this work, we wanted to avoid any effect that could mix angular scales because we are interested in characterizing the angular power spectra. Therefore, we developed the ray-tracing algorithm presented in Appendix B to circumvent the use of interpolation techniques that could, in principle, lead to such unwanted artificial effects. In this appendix, we show that the synthetic polarization maps may indeed depend on the choice of interpolation parameters. It is not our intention to conduct a detailed analysis of the differences that may appear on a large sample of maps or to study the use of different interpolation schemes. Such a detailed analysis, while important, is beyond the scope of this paper. Therefore, we carry out our analysis for the case of one interpolation scheme and one observer position in one simulated galaxy (corresponding to the maps shown in Fig. 2).

To synthesize maps from a Cartesian grid using interpolation, we proceed as follows. We consider a spherical grid surroundingthe observer. This grid has an angular sampling determined by a HEALPix map with resolution N_side = 32, and a radialsampling of 30 pc (approximately the size of the smallest cells in N18 simulations). At each node of the grid, we evaluate the dust density and the magnetic field vector through an interpolation scheme (a multiquadric radial-basis interpolation, working in three dimensions) and seek the information from the n closest neighbors. The observables are then obtained through line-of-sight integration according to the mid-point rule.

We repeat the synthesis of the I, Q, and U maps for different values of n: 4, 7, 16, 27, and 36. First, we notice that the individual pixel values of the Stokes parameters change with the choice of n. Indications for convergence appear for large values of the closest neighbors involved in the interpolation. This is best illustrated in Fig. C.1 where we present the relative difference of the I, Q, and U maps as estimated for different n and arbitrarily compared to n = 36.

Fig. C.1 Summarizing statistics of the relative differences between maps obtained using n neighbors in the interpolation scheme and n =36. For each value of n we show themedian of the relative differences and the error bars are such that they contain 68 per cent of the data. The observables (X) are shown as indicated by the legend. We slightly shift horizontally the values for I and U for a better visualization.

In the following, we explore the effect of these differences on the polarization power spectra. We choose a sky fraction f_sky = 0.5 and generate the mask corresponding to each intensity map following the same procedure as in the core of the paper. We then define a common mask for all n by multiplying the individual masks. We apply this mask to the five sets of polarization maps, and compute the polarization power spectra using Xpol. Because the maps are at a lower resolution than the sample presented in the main text, we consider the multipole bins between 40 and 80 with a bin width of 5. Each set of maps leads to a set of auto- and cross-power spectra, for which we compute the ${\mathcal{R}}_{EB}$ ratio and the r^TE correlation coefficients, and we fit the spectra with power-law functions. We show these values as a function of n in Fig. C.2. The power spectrum characteristics vary depending on the specific choice of the used interpolation scheme. This effect justifies our choice not to use a map-making algorithm that involves interpolation.

Fig. C.2 Power spectrum characteristics and number of neighbors in the interpolation scheme. We show the slopes of the TT, EE, BB, and TE spectra (top), ${\mathcal{R}}_{EB}$ (middle) and rTE (bottom) as a function of n.

References

All Tables

All Figures

	Fig. 1 Normalized distribution (top) and cumulative distribution (bottom) of the mass-weighed magnetic field strength (log10 (|B|)) in the disk region for each simulation as indicated by the legend. We define the disk region as a vertical cylinder centered on the center of the simulation box, with radius of 15 kpc and extending up to 500 pc above and below the plan Z = 0.
In the text

	Fig. 2 Example of synthetic polarization maps of thermal dust emission. The observer is immersed in the Tsf galaxy at 8 kpc from the galactic center and at angular coordinate ϕ⊙ = 300°.
In the text

Fig. 3 Examples of polarization (pseudo) power spectra (${\mathcal{D}}_{\mathcal{l}}$) obtained for the polarization maps shown in Fig. 2 and for the six values of the retained sky fraction (fsky) in the analysis; from 0.2 to 0.7 (top to bottom). The gray-shaded areas are excluded from our analysis as explained in the text. The straight lines are power-law power spectra fitted in the white area, i.e. the range 60 ≤ ℓ ≤140.

In the text

	Fig. 4 Histograms of REB and rTE as measured from all galaxy setups, observer positions and sky fraction together. The red cross marks the location $(r^{TE},{\mathcal{R}}_{EB})=(0.4,2.0)$, typical values deduced from Planck data at 353 GHz.
In the text

Fig. 5 2D histograms of combinations of parameters from power-law fits used to characterize the shape of the power spectra for the full sample; all galaxy setups, observer position and sky fractions merged together. Panel a: from top left to bottom right we have $({\log }_{10}(A_{80}^{EE}),{\log }_{10}(A_{80}^{BB}))$, (αEE, αBB), $({\alpha }^{EE},{\chi }_{EE}^2)$ and $({\alpha }^{BB},{\chi }_{BB}^2)$. The main diagonals are shown in top panels for reference. In the bottom panels the horizontal lines mark values of ${\chi }_{XX}^2=1$ and the vertical lines show the median of the spectral indices. Panel b: ${\log }_{10}({\mathcal{R}}_{EB})$ versus ${\log }_{10}(A_{80}^{EE}/A_{80}^{BB})$. The horizontal and vertical lines mark ratio values of 1 and 2 and are given for visual reference. The main diagonal is also shown.

In the text

	Fig. 6 Summarizing statistics of ${\mathcal{R}}_{EB}$ (left) and rTE (right) per galaxy and per sky fraction fsky. Asymmetric error bars around the median encompass 68% of the data points. For visual reference, the horizontal dashed lines mark the values of ${\mathcal{R}}_{EB}=1$ and 2 to the left and rTE = 0 and 0.4 to the right.
In the text

	Fig. 7 Zoom in the set of maps presented in Fig. 2 to illustrate the alignment between the projected plane-of-the-sky magnetic field orientation and projected density structure orientation. The background color shows the intensity map while the streamlines show the magnetic field orientation as inferred from the Q and U Stokes parameters.
In the text

Fig. 8 Histograms of ${\log }_{10}({\mathcal{R}}_{EB})$ (left) and rTE (right) as measured on polarization maps synthesized from the three models with feedback (top) and without feedback (bottom). Power spectra corresponding to all values of fsky are put together, making an ensemble of 432 measurements per galaxy setup. The vertical gray dot-dashed (red dashed) lines mark the values of ${\mathcal{R}}_{EB}=1$ (2) on the left and of rTE = 0 (0.4) on the right, for reference.

In the text

Fig. 9 Variations of the averaged weighted emission intensity (Ī) (top left), ${\mathcal{R}}_{EB}$ (bottom left) and rTE (bottom right) as a function of the azimuthal position of the observer (ϕ⊙ [°]) for the different galaxies (indicated in the captions) and for different values of sky fraction (fsky, indicated by the colors in the legend). Top right: evolution of the averaged dust density (normalized to the sample median) in spheres of radii (100 pc, 200 pc and 500 pc) around an observer as a function its azimuthal coordinate (ϕ⊙). The colors represent different values of the radius of the sphere surrounding the observer as informed by the legend. Spheres are independent if their ϕ⊙ differ by at least 1.43°, 2.86° and 7.16° for radius Rsph of 100 pc, 200 pc, and 500 pc, respectively. The vertical gray bands indicate the range of ϕ⊙ values in which we consider the observer to be located in the arms (see Sect. 4.4.1).

In the text

Fig. 10 Scatter plots of strength of the local magnetic field (in μG) versus ${\mathcal{R}}_{EB}$ (left) and rTE, measured from the Twf (green) and Tsf (purple) galaxies for fsky = 0.7. The strength of the field is the volume-weighed average of the field of all cells having their centers within a sphere radius of 500 pc centered on the observer. Gray and red vertical lines are the same as in Fig. 8. The gray tine solid line show the medians of the data points and the shacked black lines show running means along |B | considering 10 data points at a time.

In the text

	Fig. 11 Normalized histograms of ${\log }_{10}({\mathcal{R}}_{EB})$ (left) and rTE (right) as measured on polarization maps synthesized from the inside of bubbles detected in the three models with feedback. Power spectra corresponding to all values of fsky are put together. Vertical lines are the same as in Fig. 8.
In the text

Fig. 12 Scatter plots of strength of the magnetic field (in μG) in the cone of observation versus ${\mathcal{R}}_{EB}$ (left) and rTE (right), measured from the Twf (green) and Tsf (purple) galaxies for fsky = 0.5. The strength of the field is the mass-weighed average of the field for all cells having their centers in the cone centered on the observer, with an opening angle of 60° around the z axis and with height lower than 1 kpc. Vertical and shacked lines are built the same way as in Fig. 10.

In the text

	Fig. A.1 Illustration of the observer’s circular excursion around the galaxy center on a top-down view of the Tsf galaxy model considered in the paper. A black star is placed every 5° at 8 kpc from the galaxy center. The background color trace the gas density.
In the text

Fig. A.2 Convention for observer angular position and AMR grid rotation. The third dimension (z axis) is dropped as it is left unchanged during our coordinate transformation. The observer is at (x, y) = (0, 0) and the galaxy center at (xAMR, yAMR) = (0, 0). The AMR grid is rotated clockwise by ϕ⊙ about the galaxy center and then shifted at (x, y) = (8, 0) to make the circular excursion of the observer.

In the text

	Fig. B.1 Computing the distance di from the observer to the intersecting point (Ii) of a sightlinewith the bounding plane i of a cubic cell.
In the text

Fig. B.2 Integrated path length through the full simulation box (top panels) and through only the observer cells (lower panels). The full box is of 60 kpc on a sideand is integrated from an observer placed at R⊙ = 8 kpc from thecenter. The observer is immersed in the Tsf galaxy at 8 kpc from the galactic center and his angular coordinate is is ϕ⊙ = 180° for the left panels and ϕ⊙ = 60° for the rightpanels.

In the text

Fig. C.1 Summarizing statistics of the relative differences between maps obtained using n neighbors in the interpolation scheme and n =36. For each value of n we show themedian of the relative differences and the error bars are such that they contain 68 per cent of the data. The observables (X) are shown as indicated by the legend. We slightly shift horizontally the values for I and U for a better visualization.

In the text

	Fig. C.2 Power spectrum characteristics and number of neighbors in the interpolation scheme. We show the slopes of the TT, EE, BB, and TE spectra (top), ${\mathcal{R}}_{EB}$ (middle) and rTE (bottom) as a function of n.
In the text