© The Authors 2026

1 Introduction

The transition from the warm neutral medium (WNM) to the cold neutral medium (CNM) is a critical step in the interstellar medium (ISM) cycle, and pivotal for the formation of molecular clouds and stars. This transition, largely governed by the thermal instability (e.g., Field 1965; Wolfire et al. 2003; Audit & Hennebelle 2005; Saury et al. 2014; Godard et al. 2024), marks a key stage in the condensation of diffuse gas into denser, star-forming regions. Despite extensive theoretical and numerical investigations, the precise mechanisms regulating this phase change remain incompletely understood. Among the factors influencing this transition, magnetic fields have emerged as a possible player.

The polarization of thermal dust emission, whicH is extensively mapped by the Planck satellite (Planck Collaboration Int. XXXII 2016), has offered high-resolution insights into the magnetic-field morphology, particularly in relation to filamentary structures in the diffuse ISM. These studies revealed statistical alignments between magnetic fields and H i filaments, suggesting a strong coupling between magnetic-field orientations and the physical processes governing the ISM. Magnetic fields in the ISM have been increasingly recognized as a key factor in shaping its structure and dynamics, particularly through their interaction with neutral gas phases. Recent work by Clark et al. (2015) and Clark & Hensley (2019) used H i observations and Planck dust-polarization data to reveal alignment between the morphology of H i filaments and the orientation of magnetic fields. These studies provided evidence that filamentary structures in the CNM may trace the structure of the local magnetic field.

While starlight and dust-polarization studies have advanced our understanding of the magnetic field’s large-scale morphology, another proxy of the magnetic field has been less studied for the local ISM. Synchrotron emission, arising from relativistic electrons spiraling along magnetic-field lines, provides a complementary perspective. Synchrotron radiation traces the component of the magnetic-field perpendicular to the line of sight (LOS), offering insights into its strength and orientation (Rybicki & Lightman 1979; Longair 2011; Padovani et al. 2021). Observations at low radio frequencies are particularly valuable, as they capture Faraday rotation effects, wherein the plane of polarized light rotates as it propagates through a magnetized, ionized medium (Burn 1966; Sokoloff et al. 1998; Ferrière et al. 2021). This effect depends on the thermal electron density and the magnetic-field component along the LOS, making it a powerful tool for probing the magneto-ionic medium. The degree of rotation experienced by a polarized wave is quantified by the Faraday depth, which depends on the thermal electron density, n_e, the magnetic-field component along the LOS, B_∥, and s, whicH is the path length.

The technique of Faraday tomography (Brentjens & de Bruyn 2005) exploits these low-frequency radio observations to decompose polarized synchrotron emission by Faraday depth, enabling a non-spatial 3D reconstruction of the magneto-ionic structure of the ISM (Haverkorn & Heitsch 2004; Van Eck 2018; Erceg et al. 2022; Erceg et al. 2024). Recent advancements in radio facilities, such as the Low-Frequency Array (LOFAR; e.g. van Haarlem et al. 2013) have significantly improved the sensitivity and resolution of low-frequency synchrotron-polarization observations, bringing new questions concerning the interplay between magnetic fields and the thermal instability.

A growing body of research seeks to connect these radio observables with neutral gas structures (e.g., Van Eck et al. 2017). Observationally, Bracco et al. (2020) demonstrated that CNM structures, extracted via phase decomposition of 21 cm data with the Regularized Optimization for Hyper-Spectral Analysis (

ROHSA

) code (Marchal et al. 2019), exhibit stronger spatial correlations with Faraday tomography data compared to WNM structures. However, numerical simulations by Bracco et al. (2022) reported an inverse trend, with WNM showing greater correlation with synchrotron data. This discrepancy was attributed to the low CNM mass fraction (1 %) in the simulated environment and the absence of instrumental effects in the analysis, highlighting the need for a more systematic investigation. Besides this, depolarization canals, which are narrow regions of reduced polarized intensity observed in LOFAR data (e.g., Jelić et al. 2015; Erceg et al. 2022), have been shown to align with filamentary H i structures in the CNM (Jelić et al. 2018; Clark & Hensley 2019), further supporting a close link between small-scale magnetized structures in the neutral ISM and features seen in polarization observations (e.g., Haverkorn & Heitsch 2004; Jelić et al. 2015). In particular, this behavior was observed in the 3C196 field, which was used in this study for comparison. However, it only represents a single LOS through the Galaxy and therefore cannot fully capture the morphological and physical diversity of the diffuse ISM.

The study presented here revisits these questions using a suite of magneto-hydrodynamic (MHD) simulations (Rashidi et al. 2024) that span a range of CNM mass fractions (10 %–50 %), providing a more representative depiction of the diffuse ISM. In addition, we incorporated a simplified instrument model to simulate telescope beam effects and instrumental noise, ensuring a more realistic comparison between simulations and observations.

To achieve this, we developed a Julia-based (Bezanson et al. 2017) code, Mock Observation Of Synchrotron Emission (

MOOSE

), to produce synthetic observations of polarized synchrotron emission from MHD simulations, perform Fara-day tomography, and compute spatial correlations between H i phases and Faraday observables. This approacH is intended to test whether low-frequency synchrotron emission can provide new insights into magnetic-field-driven effects in ISM phase transitions, bridging the gap between theoretical predictions and observational constraints.

The paper is organized as follows. Sect. 2 features a review of the physical properties of the 3C196 field, which was used throughout the analysis for comparison with numerical simulations. In Sect. 3, we present the methods used to analyze the correlation between Faraday tomography observables and neutral gas structures. This includes the histogram of oriented gradients (HOG) method (Soler et al. 2019), the process of generating synthetic synchrotron polarization maps using the

MOOSE

code, and the impact of instrumental effects such as noise filtering and beam convolution. In Sect. 4, we present our findings, focusing on the dependence of the H i-Faraday correlation on CNM fraction, turbulence, and magnetic-field topology. We compare the results of different simulation setups and contrast them with the LOFAR Two-metre Sky Survey (LoTSS) observations of the 3C196 field (Erceg et al. 2022). In Sect. 5, we examine the limitations of our simulations and explore possible physical processes – such as electron-density prescription, cosmic-ray transport, and supernova-driven turbulence – that may influence Faraday tomography observables and their correlation with H i data. Finally, in Sect. 6, we summarize the key findings.

2 The 3C196 field of view

2.1 Observations

The 3C196 field is one of the primary fields of view of the LOFAR epoch-of-reionization key science project. It is located in a diffuse region at high Galactic latitude (l = 171^◦ , b = +33^◦) centered on the extragalactic source 3C196 (Jelić et al. 2015; Jelić et al. 2018). The angular size of the field is 64 deg². Observations were carried out between 115 and 189 MHz with 3.2 kHz spectral resolution. It is at the core of the analysis of Bracco et al. (2020), which revealed morphological correlations between low frequencies radio polarization and H i data from the Effelsberg– Bonn H i survey (EBHIS; Kerp et al. 2011; Winkel et al. 2016). In their study, Bracco et al. (2020) used an early single-pointing version of the 3C196 field LOFAR observations, which shows a sharp increase in noise at the edges of the field of view. Later on, Erceg et al. (2022) presented a complete Faraday tomographic mosaic of the LoTSS survey conducted with the LOFAR telescope, whicH includes the 3C196 field. Observations were conducted between 120 and 168 MHz with a resolution of 97.6 kHz. The mosaicking strategy, involving overlapping fields of view, results in a more uniform noise level across the observed area.

In this work, we reproduced the data-analysis methodology of Bracco et al. (2020). To do so, we used the LOFAR data product of Erceg et al. (2022) and, as Bracco et al. (2020), the 21 cm data from EBHIS. The data products have a respective angular resolution of 4′ for LOFAR and 10′8 for EBHIS. To compare the two datasets, the LOFAR polarimetric data cubes were convolved at the EBHIS resolution and then projected onto the EBHIS grid.

2.2 Physical properties of the 3C196 field

In this section, we describe the general physical properties of the 3C196 field, which we compared to the simulations. We are mostly interested in the multiscale properties of the H i data and the intensity of the magnetic field.

Based on EBHIS 21 cm data (Winkel et al. 2016), the average H i column density in the 3C196 field is N_{H i} = 4.2 × 10²⁰ cm⁻². The H i phase decomposition done in Bracco et al. (2020) using

ROHSA

(Marchal et al. 2019) indicates a WNM column density of N_WNM = 2.4 × 10²⁰ cm⁻². The fraction of the H i in the CNM phase, defined as f_CNM = N_CNM/N_{H i}, is 38%.

The

ROHSA

decomposition also allows us to estimate the H i velocity dispersion along the LOS. Here, we took the WNM as the reference as it is volume-filling. The WNM velocity dispersion σ_u in the 3C196 field is 16.3 km s⁻¹. This velocity dispersion results from the combination of thermal and turbulent contributions and can be expressed as ${\sigma }_u^2={\sigma }_{\text{turb}}^2+{\sigma }_{\text{therm}}^2$. Assuming a reference temperature of 7000 K for the WNM, the thermal velocity dispersion is σ_therm ≈ 7.6 km s⁻¹, leading to a turbulent velocity dispersion of σ_turb ≈ 14.4 km s⁻¹. Following Kolmogorov (1941), the turbulent velocity dispersion is related to scale l as $${\sigma }_{\text{turb}}={\sigma }_{u,1\text{pc}}{\left(\frac{l}{1\text{pc}}\right)}^{1/3},$$(1)

where σ_u,1 pc is the 1D turbulent velocity dispersion along the LOS at a scale of 1 pc.

According to the 3D-dust map of Edenhofer et al. (2024), the Galactic interstellar matter observed in the 3C196 field is distributed over a depth of ∼350 pc and located approximately 200 pc from the Sun, corresponding to the wall of the Local Bubble. Assuming the WNM is filling the whole depth of 350 pc, we obtain σ_{u,1 pc} = 2.0 km s⁻¹. At the 50 pc scale (the simulation scale), the extrapolated 1D velocity dispersion is σ_{u,50 pc} = 7.4 kms⁻¹, yielding a 3D turbulent velocity dispersion of σ_3D = 12.8 km s⁻¹.

Moments 0, 1, and 2 of the Faraday spectra of the 3C196 field are displayed in Fig. 1. Previous analyses of this region (Jelić et al. 2015; Turić et al. 2021; Erceg et al. 2022) indicate that the magnetic field is largely dominated by its plane-of-the-sky component, with only a modest contribution along the LOS. Therefore, the observed Faraday depths do not necessarily imply locally enhanced LOS magnetic-field strengths, but rather result from the interplay of electron density and path length along this particular LOS. According to Heiles & Troland (2005), the magnetic-field strength at high Galactic latitude is estimated as 6 ± 2 µG based on 21 cm Zeeman-splitting measurements.

3 Methods

3.1 Description of numerical simulations

In this study, we used ideal-MHD simulations designed with the

RAMSES 

code (Teyssier 2002; Fromang et al. 2006), which uses adaptive mesh refinement (AMR) up to 1024³. We used a version of these simulations computed on a Cartesian grid of 256³ pixels. The box’s physical size is L = 50 pc aside, leading to a resolution of ∆x ≈ 0.20 pc.

The physical conditions we want to reproduce are the ones of atomic gas clouds or those transitioning to molecular gas as expected on the walls of the Local Bubble (Inoue & Inutsuka 2016). In this diffuse medium (0.1 to 100 cm⁻³, Ferrière (2001)), the self-gravity energy density is significantly lower (two-to-three orders of magnitude) than turbulent and magnetic energy densities; hence, self-gravity is excluded from the simulations. The cooling function is based on Audit & Hennebelle (2005) and closely follows Wolfire et al. (2003). It includes photo-electric heating of dust grains and cooling by electron recombination onto dust grains, as well as collisional excitation of Ly α, C II, and OI emissions.

To maintain turbulence in the simulation, a Fourier space forcing was modeled using an Ornstein–Ulhenbeck stochastic process (Eswaran & Pope 1988; Saury et al. 2014). A turbulent velocity field was added on a large scale (1 < k < 3), peaking at k = 2 (half the box size, L). The forcing was parametrized by its amplitude, A_f,RMS, and the compressive-to-solenoidal mode decomposition parameter, χ, where χ = 0 only contains solenoidal modes and χ = 1 only compressible modes.

In this study, we used the set of 14 simulations described in Table 1. The mean B-field intensity of all simulations is 7.64 µG. The ratio of compressible to solenoidal modes of the forcing was varied from χ = 0 to χ = 1, and its amplitude, A_f,RMS, was varied between 9000 and 72 000, producing 3D velocity dispersions from 4.3 to 14.5 km s⁻¹.

Each simulation began with WNM conditions in place, a density of n_H = 1 cm⁻³, and a temperature of T = 8000 K. The magnetic field was initiated uniformly along the x-axis, which gives a prevailing direction for the B field in the simulation box; ∇ · B = 0 was guaranteed by the use of the constrained transport method. Due to the nature of the periodic box, the coherence of the magnetic field along the x-axis was maintained at each time step.

As the simulations evolve, shear flows and shock waves emerge in the volume, causing local pressure variations in the WNM leading to the formation of CNM structures. We let simulations proceed until a steady state was achieved when velocity dispersion, gas-phase mass fractions, density, and temperature PDFs reached a stationary state, within at least three dynamical timescales, t_dyn = L/σ_u,3D, where σ_u,3D is the velocity dispersion in 3D. An illustration of one simulation is given in Fig. 2, where two slices of the density field are shown, one along and one perpendicular to the mean magnetic field, represented here using the line integral convolution (LIC) technique (Cabral & Leedom 1993).

For each simulation, we decomposed the gas in the three main H i components: CNM, lukewarm neutral medium (LNM), and WNM. As in Saury et al. (2014), we notice that CNM structures form preferentially in over-pressured and less turbulent WNM regions. The overall CNM mass fraction in our set of simulations tends to increase with the rate of compressible modes, and even more so for less turbulent ones. This is illustrated in Fig. 3, where the CNM mass fraction is plotted for the last time step of each simulation. Table 1 also provides the CNM mass and volume fractions for the last time step of each simulation, together with the root-mean-squared (RMS) values for B_RMS, u_RMS, and the velocity dispersion, σ_u,1pc, of the gas at 1 pc.

In this paper, we regularly refer to simulation 7, which we used as our case study throughout our analysis. This simulation has a CNM fraction of 29% and σ_u,1pc = 1.1 km s⁻¹, both of which are slightly lower than the values of the 3C196 field, but still representative of the diffuse ISM. Moreover, the Fara-day depth distributions of the simulations and of the 3C196 field, shown in Fig. 4, exhibit comparable orders of magnitude, namely a few rad m⁻².

3.2 Synthetic observations

In order to reproduce the analysis of Bracco et al. (2020), which compared polarized synchrotron and 21 cm emission data of the 3C196 field, we had to create corresponding synthetic observations for each numerical simulation. This section describes how these were created, including the effects of Faraday rotation andthe LOFAR beam on the polarized synchrotron emission.

3.2.1 Synchrotron emission

In order to produce synthetic observations of synchrotron emission and its polarization, whicH is the basis of the Faraday tomography, we used Eqs. (1)–(9) of Padovani et al. (2021) derived from Ginzburg & Syrovatskii (1965). As cosmic-ray electrons (CRe) travel through the ISM, they lose energy due to interactions with matter, magnetic fields, and radiation (Longair 2011). These interactions modify the cosmic-ray energy spectrum, j_e(E), where E represents the CRe energy.

For the CRe energy spectrum, j_e(E), we used the same power-law definition as in Padovani & Galli (2018) and Orlando (2018): $$j_e(E)=j_0\frac{E^a}{{\left(E+E_0\right)}^b},$$(2)

where j₀ = 2.1 × 10¹⁸ eV⁻¹s⁻¹cm⁻²sr⁻¹, E₀ = 710 MeV, a = −1.3, and b = 1.9. In our case, j_e(E) was considered spatially constant in the simulated box.

Noting m_e as the electron mass, e the electron charge, and c the speed of light, the synchrotron-specific emissivity at a position, s, can be divided into two linearly polarized components: one parallel and one perpendicular to the magnetic-field component orthogonal to the LOS, B_⊥: $${\epsilon }_v^{\Vert }(s)=\frac{\sqrt{3}{e}^3}{2m_e{c}^2}{{\displaystyle \int }}_{m_e{c}^2}^{+\infty }\frac{j_e(E)}{v_e(E)}{B}_{\perp }(s)(F(x)-G(x))dE,$$(3)

and $${\epsilon }_v^{\perp }(s)=\frac{\sqrt{3}{e}^3}{2m_e{c}^2}{{\displaystyle \int }}_{m_e{c}^2}^{+\infty }\frac{j_e(E)}{v_e(E)}{B}_{\perp }(s)(F(x)+G(x))dE\text{,}$$(4)

where $x=\frac{v}{v_c}$ and ν_c is defined as the frequency at which the synchrotron intensity of CRe at energy, E, is the greatest: $$v_c(B(r),E)=\frac{3e}{4\pi m_{e}c}{B}_{\perp }(r)\left(\frac{E}{m_e{c}^2}\right).$$(5)

Synchrotron emission is thus $$I_v(s)={{\displaystyle \int }}_0^L{\epsilon }_v^{\Vert }(s)+{\epsilon }_v^{\perp }(s)d(s).$$(6)

3.2.2 Faraday rotation

In this paper, we address the scenario where Faraday rotation is completely intertwined with synchrotron emission, leading to differential Faraday rotation (e.g., Sokoloff et al. 1998) as in Bracco et al. (2022). Each layer within the simulated cubes contributes to both synchrotron emission and Faraday rotation.

Stokes Q_ν and U_ν for synchrotron emission are given by the following equations: $$U_v={{\displaystyle \int }}_0^L\left({\epsilon }_v^{\Vert }(s)-{\epsilon }_v^{\perp }(s)\right)\sin \left(2\left({\text{Ψ}}_{\text{loc}}(s)+{\text{ΔΨ}}_v(s)\right)\right)ds$$(7)

and $$Q_v={{\displaystyle \int }}_0^L\left({\epsilon }_v^{\Vert }(s)-{\epsilon }_v^{\perp }(s)\right)\cos \left(2\left({\text{Ψ}}_{\text{loc}}(s)+{\text{ΔΨ}}_v(s)\right)\right)ds,$$(8)

where Ψ_loc is the polarization angle at the position s defined as the orientation of B_⊥ rotated by 90^◦. ∆Ψ_ν describes the Fara-day rotation, whicH is the integral along the LOS between the observer and s. The Faraday rotation depends on the square of the observation wavelength, λ; the electron density, n_e; and the magnetic field along the LOS. It is also expressed as $${\text{ΔΨ}}_v(s)=\varphi (s){\lambda }^2,$$(9)

where $\lambda =\frac{c}{v}$ and ϕ is the Faraday depth, a quantity independent of wavelength that describes the physical quantities responsible for the angle rotation from the emitting source toward the observer, corresponding to the cumulative Faraday rotation experienced by the polarized emission along the LOS: $$\frac{\varphi (s)}{\left[{\text{rad\hspace{0.17em}mm}}^{-2}\right]}=0.81{{\displaystyle \int }}_0^s\frac{n_e}{\left[{\text{cm}}^{-3}\right]}\frac{B}{[\mu \text{G}]}\cdot \frac{ds}{[\text{pc}]},$$(10)

where ds = ∆x in the discrete case of our numerical simulations. From Equations (7) and (8), we define the polarized intensity, P, as $$P=Q_v+i{U}_v.$$(11)

P is commonly computed as an analog of the Fourier transform into the Faraday spectrum: $$F(\varphi )=\frac{1}{\pi }{{\displaystyle \int }}_{{\lambda }_{\min }^2}^{{\lambda }_{\max }^2}P\left({\lambda }^2\right)e^{-2i\varphi {\lambda }^2}d\left({\lambda }^2\right).$$(12)

For the limit case, where the integration is made all along the LOS, i.e., s = 0 and s + ds = L, ϕ becomes the rotation measure (RM).

The Faraday spectrum is represented as a polarized brightness temperature by taking the modulus of F(ϕ), making it a non-negative quantity. Hereafter, we use F(ϕ) as representative of the modulus of Eq. (12).

3.2.3 Electron density

In order to produce synthetic synchrotron observations including Faraday rotation, we needed to estimate the density of thermal electrons, n_e, at each position in the cube. The numerical simulations were thus post-processed to estimate n_e in units of cm⁻³ following Eq. (6) of Bracco et al. (2022), which derives the following from Wolfire et al. (2003) and Bellomi et al. (2020): $$n_e\approx 2.4\times {10}^{-3}{\zeta }_{16}^{0.5}{T}_2^{0.25}\frac{G_{\text{eff}}^{0.5}}{{\varphi }_{\text{PAH}}}+n_{\text{H}}{X}_{C^{+}},$$(13)

where ζ₁₆ represents the total ionization rate per hydrogen atom due to energetic photons (including extreme UV and soft X-rays) and CRs expressed in units of 10⁻¹⁶ s⁻¹, T₂ is the temperature of the gas expressed in units of 100 K, ϕ_PAH denotes the recombination rate of electrons onto small dust grains such as polycyclic aromatic hydrocarbons (PAHs). Furthermore, $X_{{\text{C}}^{+}}$ is the relative abundance of ionized carbon compared to the number of hydrogen atoms (n_H). Following Bracco et al. (2022), we used $X_{{\text{C}}^{+}}=1.4\times {10}^{-4}$ and ϕ_PAH = 0.5, G_eff = 1, and ζ = 2.5 × 10⁻¹⁶ s⁻¹ (Padovani et al. 2022) as representative values for the diffuse ISM. The impacts of these choices on our results are discussed in Sect. 5.

3.2.4 Rotation measure

Before computing mock LOFAR observations, it is useful to evaluate the range of Faraday depth that we expect from our simulations. For example, Fig. 5 shows the distributions of RM for simulation 7 (computed using Eq. (10) and integrated over the full LOS) for cubes produced with the perpendicular, parallel, or 45^◦ LOS with the mean direction of B. For each orientation, the range of RM values is centered on specific values, illustrating the fact that the strength of B_∥ is different in each case. When the LOS is perpendicular to B, the RM values are very close to zero, as expected. For the case where the mean magnetic field is aligned with the LOS, ⟨RM_∥⟩ = 5.5 rad m⁻². In the intermediate case, where the simulation cubes are produced with a LOS at 45^◦ with the mean direction of B, the RM distribution is intermediate, with a mean RM of 4.8 rad m⁻². In the three projections, the dispersion of RM is only about 1–2 rad m⁻². These RM values fall within the range observed with LOFAR in the 3C196 field (Jelić et al. 2018; Bracco et al. 2020; Erceg et al. 2022). They are also similar to what was reported by Hutschenreuter et al. (2022) in the 3C196 field, with most RM values ranging from −5 to 10 rad m⁻². This is interesting as the RM values from Faraday measurements on extragalactic point sources of Hutschenreuter et al. (2022) are probing greater physical depth, as they estimated.

3.2.5 Mock LOFAR observation

In order to compute synchrotron Stokes parameters, we developed a Julia package called

MOOSE

¹, whicH is described in Appendix A. For each of the 14 simulations, we used the B-field and n_e cubes to compute I_ν, Q_ν, and U_ν based on Equations (3)–(9).

MOOSE

is able to treat any frequency range, but in this work we focused on LoTSS frequencies used for the study of the 3C196-field (Erceg et al. 2022); i.e., ν is in the range of [120, 167] MHz, with a frequency resolution of δν = 98 kHz. As LOFAR interferometric observations are subject to large-scale filtering, we wanted to simulate this effect. To do this, we used a Gaussian kernel, G_ker, of standard deviation σ_LS. We applied a convolution on a given image, A_i,j, as follows: $${\tilde{A}}_{i,j}=A_{i,j}-A_{i,j}⊛G_{\ker }\left({\sigma }_{\text{LS}}\right),$$(14)

where ${\tilde{A}}_{i,j}$ is the image with the large scales are removed.

With Q_ν and U_ν calculated, we computed F(ϕ) using a Julia version of the rotation-measure synthesis (

RM-Synthesis

) (Brentjens & de Bruyn 2005) implemented in

MOOSE

. Due to the frequency range, which we chose as that of the LoTSS survey (Erceg et al. 2022), the Faraday resolution is 1.14 rad m⁻² (see Eq. (61) in Brentjens & de Bruyn 2005). We computed F(ϕ) from ϕ = −10 rad m⁻² to 10 rad m⁻² with steps of 0.25 rad m⁻² to include the range of values of RM found in our simulations (see Fig. 5).

Figure 6 (left) shows an example of polarized intensity at a specific Faraday depth, F(ϕ), for simulation 7 in the case where the mean B field is along the LOS. For all our simulations, F(ϕ) displays complex structures that reflect the interplay between the electron density and magnetic-field morphologies.

3.2.6 Mock 21 cm line observation

To simulate 21 cm line observations of neutral hydrogen, we computed the brightness temperature, T_B(u), as a function of velocity, u, along each LOS using the following radiative transfer equation: $$T_B(u)={\displaystyle \sum }_{k}T(k)\left(1-e^{-\tau (k,u)}\right)\cdot \exp \left(-{\displaystyle \sum }_{i<k}\tau (i,u)\right).$$(15)

Here, the l index runs over the cells along the LOS, τ(l, u) is the 21 cm optical depth at position l and velocity u, and T(l) is the gas temperature (taken as a proxy for the spin temperature). In this formulation, the term $T(l)\left(1-e^{-\tau (l,u)}\right)$ represents the emission from cell l and velocity u, while the exponential factor exp $\left(-\underset{}{{\displaystyle {\sum }_{i<l}}}\tau (i,u)\right)$ describes the attenuation of this emission by the gas located between position k and the observer.

We assumed that the opacity τ(l, u) is $$\tau (l,u)=\frac{N_H(l,u)}{C{T}_s(l)},$$(16)

where C = 1.82243 × 10¹⁸ cm⁻² K⁻¹ (km/s)⁻¹ and N_H(k, u) is the column density in cell l expressed as a Gaussian with thermal broadening centered on the local velocity, u = v(l); and normalized by the local density, n_H(l): $$N_H(l,u)=-\frac{n_H(l)}{\sqrt{2\pi }\text{Δ}}\times \exp \left(-\frac{{(u-v(l))}^2}{2{(\text{Δ}(l))}^2}\right)dl,$$(17)

with $$\text{Δ}(l)=\sqrt{\frac{kT(l)}{\mu m}}.$$(18)

In theory, ∆ should also include the turbulent velocity dispersion at scales smaller than the cell scale. We made this simplification as the typical turbulent velocity dispersion at the scale of a cell (0.2 pc) is ∼0.6 km s⁻¹, assuming a typical value of σ_u,1pc = 1 km s⁻¹ and using Eq. (1). This would become significant only for gas at T < 45 K, whicH is very rare in the simulations used here.

We produced synthetic 21 cm observations for each simulation using this formalism. An example of a 21 cm channel map is shown in Fig. 6, right. As for F(ϕ), the 21 cm channel maps show a complex multiscale structure with several narrow and twisty filaments. On the other hand, contrary to F(ϕ), where the structure is very different depending on the orientation of B with respect to the LOS, the structure in T_B(u) is morphologically very similar in all directions.

3.2.7 H i phases cubes

Bracco et al. (2020) used

ROHSA

to decompose the EBHIS 21 cm observation of the 3C196 field into phases of the neutral hydrogen. They defined the CNM as Gaussian components with σ ≤ 3 km s⁻¹, and the WNM with components with σ ≥ 6 km s⁻¹. The LNM is everything in between. To avoid having to run

ROHSA

on every 21 cm synthetic observation of our simulation set, we chose to segment the H i phases based on the 3D temperature field directly. Assuming a Mach number of 1 and 2 for the WNM and CNM, respectively (Marchal et al. 2024; Heiles & Troland 2003), these σ_u thresholds correspond to the following temperature ranges: T_CNM ≤ 500 K, 500 < T_LNM < 3500 K, and T_WNM ≥ 3500 K.

In practice, for a given phase (CNM, LNM, WNM) we identified the voxels in 3D corresponding to its temperature range; then, we produced 21 cm synthetic observations using only those voxels. Therefore, for each simulation setup we produced four 21 cm cubes for the total, H i, CNM, LNM, and WNM.

4 Results

4.1 Global properties of T_B(u) and F(ϕ)

The H i and Faraday structures of simulation 7 are illustrated in Figures 7 and 8 with views parallel and perpendicular to the mean B field. On the left, these figures show the column-density maps of the H i components (total, CNM, and WNM) highlighting the striking filaments of the CNM as well as the diffuse structure of the more volume filling WNM. On the right, we show the moments (M₀, M₁, and M₂) of the Faraday cube; see Appendix B. These can be directly compared to the moments of the 3C196-field LoTSS observations (Fig. 1).

The Faraday moments of the simulations (Figures 7 and 8) are strongly influenced by the orientation of the magnetic field. In the case perpendicular to the mean field, the values of M₁ remain close to zero due to the weak field component along the LOS. On the other hand, when the magnetic field is parallel to the LOS we observe values of M₁ that are consistently positive and fall within the positive range observed in the M₁ map of Erceg et al. (2022).

Regarding the spread of the Faraday spectra estimated using M₂, the values are very small in the case where B is perpendicular to the LOS, very close to theoretical RMSF of LoTSS, whicH is 1.03 ∼ rad m⁻² estimated for the theoretical RMSF of LoTSS. On the other hand, the case where B is parallel to the LOS shows M₂ values comparable to the ones estimated in the LoTSS 3C196 field.

The histograms of the first-order moment M₁ are shown in Fig. 9 (top panel). For the 3C196 field, M₁ has a mean value of approximately 0.2 rad m⁻² and a standard deviation of 0.9 rad m⁻². Among the simulations, the perpendicular case provides the closest matcH in terms of mean value, which reaches ⟨M₁⟩ = −0.2 rad m⁻², while the parallel and 45^◦ cases have ⟨M₁⟩ = 2.4 rad m⁻² and ⟨M₁⟩ = 1.7 rad m⁻², respectively.

The Faraday spread M₂ (de-biased from noise) measured in the 3C196 field has a median value of approximately 1.7 rad m⁻² (see right panel of Fig. 9). Observed values, even after de-biasing from noise, are generally slightly higher than what was recovered in the synthetic observations constructed without noise. The best agreement was found when the emission was integrated along the direction rotated by 45^◦ to the mean magnetic field, yielding a mean of ⟨M₂⟩ = 1.5 rad m⁻². This is close to the the case where B is parallel to the LOS (⟨M₂⟩ = 1.4 rad m⁻²), but with a narrower distribution. The perpendicular configuration shows lower values, as expected because the B component along this direction is very low, with ⟨M₂⟩ = 1.0 rad m⁻².

One striking feature of the LoTSS polarization data shown in Erceg et al. (2022) is the presence of long and thin depolarization canals; i.e., filamentary regions with the maximum F(ϕ), P_max, close to zero (for more details, see Sokoloff et al. 1998; Haverkorn & Heitsch 2004). Maps of P_max (LOS perpendicular and parallel to B) are shown in Fig. 10. Structures with almost zero intensity are seen in the simulations, especially in the case where the mean B field is aligned with the LOS. However, the long, straight canals detected in the LoTSS survey are not reproduced.

4.2 Correlation between LoTSS and EBHIS data

Following the same methodology of Bracco et al. (2022), we compared cubes of H i brightness temperature with Faraday cubes using the HOG method (Soler et al. 2019). The core principle of HOG is to measure the morphological alignment between two images, A_i,j and B_i,j, where i and j denote the image coordinates. It relies on the assumption that the local appearance and shape of an image can be fully characterized by the distribution of its local intensity gradients or edge directions. From this assumption, the gradients of the two images can be computed using Gaussian derivatives, yielding ∇A_i,j and ∇B_i,j. Following the definition in Soler et al. (2019), the angle between the two gradient vectors, α, can then be determined as follows: $${\alpha }_{i,j}=\arctan \left(\frac{\left(\nabla A_{i,j}\times \nabla B_{i,j}\right)\cdot \widehat{z}}{\nabla A_{i,j}\cdot \nabla B_{i,j}}\right),$$(19)

where $\widehat{z}$ is the unit vector perpendicular to the plane of sky. HOG tests whether the distribution of angle, α, is uniformly distributed or has a preferential alignment at α = 0, using the projected Rayleigh statistics defined as $$V=\frac{{{\displaystyle \sum }}_{i,j}{\omega }_{i,j}\cos \left(2{\alpha }_{i,j}\right)}{{{\displaystyle \sum }}_{i,j}{\left({\omega }_{i,j}^2/2\right)}^{1/2}},$$(20)

where ω_i,j is the statistical weight corresponding to each α_i,j. The Rayleigh statistics parameter, V, is thus the estimate of how well the two images are morphologically correlated. For more details about this technique, see Soler et al. (2019).

The observational maps do not have the same dimensions as those from the simulations. To address this, we introduced a normalization factor, V_max, defined as the maximum value of V. The maximum of V is attained when the gradients of the two maps are perfectly aligned (i.e., α_i,j = 0) and supposing the weights are all the same (i.e., ω_i,j = ω) (also detailed in Mininni et al. 2025). From Eq. (20), this yields V_max = (2N)^1/2.

Specifically, let F_i,j,ϕ represent the Faraday cube and $T_{b_{i,j,u}}$ the 21 cm brightness temperature cube, where i, j are pixel indices on the plane of the sky. For each pair (ϕ, u), we applied HOG using a local kernel around F_i,j and $T_{b_{i,j}}$. This allowed us to construct maps parameterized by (ϕ, u), illustrating the spatial correlation between specific F(ϕ) and corresponding T_b(u).

The correlation in the 3C196 field reported in Bracco et al. (2020) has since been updated using LoTSS data and the same phase decomposition of the EBHIS data done with

ROHSA

. In Fig. 11, the ratio $\tilde{V}=V/V_{\max }$ is represented on a map with axes ϕ and v_LSR. In the H i 21 cm data, two distinct velocity regions can be identified: one around v_LSR = 0 km s⁻¹ and another around v_LSR = 20 km s⁻¹, similar to the findings of Bracco et al. (2020). Along the ϕ-axis, a broad structure is observed, spanning approximately 5 rad m⁻², with a central peak around ϕ = 0 rad m⁻².

This structure appears bright across all velocities in H i and WNM. For the CNM, velocities extend from v_LSR = −10 km s⁻¹ to v_LSR = 20 km s⁻¹, with a strong correlation at 20 km s⁻¹. In contrast, for the LNM phase, this structure spans a more limited velocity range, from v_LSR = −20 km s⁻¹ to v_LSR = 0 km s⁻¹.

The highest value is $\tilde{V}\approx 0.06$, but as the maximum is sensitive to outliers, we show the 98^th percentile of $\tilde{V}$, ${\tilde{V}}^{98\%}$ for each ϕ in the lower panels of Fig. 11. Based on this HOG metric, the strongest correlation between T_B(u) and F(ϕ) is found for the CNM (also seen in the total H i) compared with the WNM and LNM components. This is compatible with the results of Bracco et al. (2020).

In order to further quantify the level of correlation of the CNM and WNM phases, we used the maximum of ${\tilde{V}}^{98\%}$ for each phase and computed the following ratio: $$\eta =\frac{\max \left({\tilde{V}}_{\text{WNM}}^{98\text{\%}}\right)}{\max \left({\tilde{V}}_{\text{CNM}}^{98\text{\%}}\right)}.$$(21)

This value reflects the strength of the dominant correlation among the two: η < 1 is indicative of a stronger correlation of the Faraday data with the CNM, while η > 1 shows a stronger correlation with the WNM. In the 3C196 field of view, we find η = 0.73, which supports the visual observation in Fig. 11 that the CNM exhibits the strongest correlation with the Faraday data.

4.3 Correlation between Faraday tomography and H i emission with numerical simulations

We applied the same HOG methodology to the synthetic observations of our simulation set. Figures 12 and 13 present the HOG results for the correlation of the H i phases with the Fara-day cube for simulation 7 for two orientations (mean B parallel and perpendicular to the LOS). For both orientations the WNM and LNM exhibit higher correlation with Faraday structures than the CNM, in contrast with the results for the 3C196 field of view.

For the case where the mean B field is aligned with the LOS (Fig. 12), the WNM and LNM have $\tilde{V}$ values peaking at 0.12. The CNM, although less correlated, reaches values ($\tilde{V}\sim 0.07$) very close to the maximum $\tilde{V}$ observed in the 3C196 field. The curve of ${\tilde{V}}^{98\%}\left(\varphi \right)$ confirms this trend, with the WNM displaying the broadest and strongest peak, followed by the LNM and CNM. The strongest correlation is observed around ϕ ≈ 3 rad m⁻², whicH is consistent with the expected Faraday depth range for a magnetic field aligned with the LOS.

For the case where the magnetic field lies in the plane of the sky (Fig. 13), the values of $\tilde{V}$ are similar, but the structure in the map of $\tilde{V}$ is much narrower in ϕ space compared to the case B parallel to the LOS. The WNM and LNM still exhibit the highest correlation, but here the CNM correlation is even lower. In this configuration, the correlation is concentrated near ϕ ≈ 0 rad m⁻², whicH is expected since a predominantly perpendicular field results in minimal Faraday rotation.

Overall, the case-study simulation produces ${\tilde{V}}^{98\%}$ values of total H i, WNM, and LNM that are approximately twice as high as those observed in the 3C196 field. The CNM phase has $\tilde{V}$ values lower than the other phases, but they reach levels consistent with observational values ($\tilde{V}\approx 0.06$).

Interestingly, for the case where B is perpendicular to the LOS (Fig. 13), the total H i displays a ${\tilde{V}}^{98\%}\left(\mu \right)$ curve with a distinctive double-peaked, U-shaped profile. As shown in Fig. 2, right, in this orientation some dense structures associated with the CNM appear misaligned or even orthogonal to the projected B field. This morphological decoupling reduces the coherence between cold H i structures and the Faraday tomographic data. The result is a decrease of $\tilde{V}$ in velocity channels dominated by CNM emission, inducing this U-shaped form in ${\tilde{V}}^{98\%}\left(\mu \right)$.

4.4 Noise effects on HOG results

To assess the impact of noise on the correlation between the H i phases and the Faraday structures, we added white noise to the mock H i and Faraday data. For H i, noise was added to each velocity channel at various levels with respect to the σ_EBHIS = 90 mK noise (Winkel et al. 2016): σ_EBHIS/5, σ_EBHIS/2, σ_EBHIS, 2 × σ_EBHIS, and 10 × σ_EBHIS. For Faraday, the noise was added to Q_ν and U_ν, with a noise amplitude tuned to reproduce specific S/Ns in Faraday depth space, including the S/N of the observational data – whicH is five – as well as both lower and higher S/N regimes.

As shown in Fig. 14, the addition of noise to the H i data predominantly affects the WNM structures. Specifically, when increasing the noise level from σ_EBHIS/5 to 10 × σ_EBHIS, the maximum correlation value for the CNM decreases by a factor of approximately 1.5, whereas for the WNM, the corresponding decrease reaches a factor of 2.5.

Two effects are at play here. First, the WNM has a broad emission line due to its larger thermal broadening. Therefore, a given column density of WNM is more spread out in velocity with lower values of T_B(u) than the same column density of CNM that will be narrow and brighter. Second, the WNM has an intrinsic diffuseness and lower spatial contrast compared to the CNM, which forms narrow, high-contrast filaments. As a result, the WNM large-scale structures with smoother gradients are more easily suppressed by noise, particularly at high spatial frequencies where gradient-based methods such as HOG are most sensitive.

Contrary to T_B(u), adding noise to Q_ν and U_ν affects both the WNM and CNM in a similar manner. As shown in Fig. 15, increasing the S/N in the Faraday space (S/N_ϕ) from two to six leads to an increase in the maximum correlation value by a factor of approximately 1.5 for both phases. In the case where the signal amplitude equals that of the noise (i.e., S/N_ϕ = 1), the resulting correlation between the Faraday cube and the H i data is effectively random (${\tilde{V}}^{98\%}\approx 0$), as expected. The case perpendicular to the mean field gives similar results.

In conclusion, noise in T_B(u) tends to decrease η (i.e., the noise decreases the correlation of WNM compared to CNM), while noise in F(ϕ) does not affect η as it reduces the absolute value of $\tilde{V}$ for all phases.

4.5 Kernel size effect on HOG results

An important parameter of the HOG method is the derivative kernel size, σ_HOG, that dictates the scale over which the spatial gradient is calculated. This parameter might have a different effect on $\tilde{V}$ for the WNM and CNM as they are not structured at the same scale; CNM structures exhibit a more small-scale filamentary morphology, whereas the WNM is present across a broad range of scales with limited small-scale contrast. In this section, we detail how we evaluated the effect of σ_HOG on η. We did so for different values of the noise level on T_B(u) as it has an effect on η (see Sect. 4.4). Figure 16 shows the variation of η as a function of σ_HOG (in pixels) for different H i noise levels added to simulation 7, for B parallel (left) and perpendicular (right) to the LOS. The observational reference field 3C196 is also shown in this figure for comparison. The values of η for the 3C196 field remain consistently below unity across all kernel sizes, σ_HOG, indicating that the CNM is systematically more correlated with the Faraday structures than the WNM.

The results on the mock data of simulation 7 show a significant variation of η as a function of σ_HOG, H i noise level, and orientation of B with respect to the LOS. For case B parallel to the LOS, the value of η is almost independent of σ_HOG, while η increases significantly with σ_HOG in the case where B is perpendicular to the LOS. Strong noise in the H i data produces an increase of η with σ_HOG as the effect of the noise affecting the WNM is gradually reduced. These variations for a single simulation show that observational results on η can depend significantly on the observational conditions (noise level, physical resolution of the observations), but also on the mean orientation of B with respect to the LOS. That said, our analysis shows that the only way to reproduce the values of η seen in the 3C196 field with simulation 7 is to add a high noise level in the H i mock data in order to lower the WNM correlation to the benefit of the CNM. This is discussed further in Sect. 5.

4.6 Variation with physical parameters

4.6.1 Turbulence strength

The HOG analysis described in the previous sections and illustrated in simulation 7 was also performed on our full simulation set. First, in order to understand the nature of the signature observed in Figures 12 and 13, we show the ${\tilde{V}}^{98\%}\left(u\right)$ curves for all simulations in Figures C.1 and C.2. We chose not to apply any instrument model (LOFAR beam, noise on Faraday and H i) to estimate the effect of turbulence level alone.

We observed that simulations with lower turbulence strength have higher $\tilde{V}$ values than simulations with stronger turbulent driving; see Fig. 17. This suggests that lower turbulence, and therefore weaker magnetic fluctuations, enhances the correlation between magneto-ionized observables and the H i phases. However, we also note that η is always higher than what is seen in observational data whatever the turbulence level (i.e., the WNM remains more correlated to F(ϕ) than the CNM). We recall that this experiment was done without any noise in the mock data and, as shown in Fig. 16, noise in the H i data tend to lower η.

4.6.2 CNM fraction

Finally, one goal of the current study was to evaluate if the fraction of the H i in the CNM phase has an impact on the measured correlation between T_B(u) and F(ϕ). To do so, we examined how η varies as a function of f_CNM. The aim is to determine whether a higher proportion of CNM in the simulated H i data systematically enhances the correlation with Faraday structures and whether this trend depends on the orientation of the magnetic field with respect to the LOS.

Figure 18 displays η as a function of f_CNM and u_rms. A horizontal dashed line at η = 1 indicates the threshold separating WNM-dominated (η > 1) from CNM-dominated (η < 1) regimes. The value observed in the 3C196 field, η_obs = 0.72, is also shown. In the perpendicular configuration, η reaches high values, particularly for simulations with low CNM fractions – exceeding η = 7 in some cases – whicH indicates a strong WNM–Faraday alignment. As f_CNM increases, η decreases but generally remains above one, confirming the persistent dominance of the WNM. In contrast, the parallel configuration shows η values consistently close to or below one across the entire range of CNM fractions with a flat trend. Only at high CNM fractions (f_CNM ≳ 35%) does the correlation approach the observational regime where CNM dominates.

To further disentangle the contributions from each phase, the two bottom panels show the respective alignments of WNM and CNM structures with the Faraday map. These panels confirm that the WNM largely dominates the correlation in the perpendicular case, while the CNM only begins to compete in the parallel configuration at high CNM fractions. This suggests that while the presence of CNM is a necessary ingredient for achieving a CNM-dominated regime, it is not the primary driver of the correlation, which depends more critically on the relative orientation of the magnetic field.

In the configuration where the LOS is aligned with the mean magnetic field, the maximum values of ${\tilde{V}}^{98\%}$ for the CNM and WNM phases follow remarkably similar trends. This suggests that the morphological structures of CNM and WNM are comparable in this case and may both trace the same Faraday structures. In contrast, the perpendicular configuration reveals a clear divergence between the CNM and WNM η values, indicating that their structures differ significantly. As also supported by Fig. 8, this behavior suggests that the WNM structures are more closely aligned with the Faraday morphology, whereas the CNM structures deviate more substantially. This asymmetry highlights the dominant role of the WNM in shaping the observed Faraday features, particularly when the magnetic field lies in the plane of the sky.

5 Discussion

5.1 Well-reproduced observation features

The numerical experiment presented here, representing the multiphase H i post-processed using a classical prescription to estimate the local ionization fraction, is able to reproduce several features seen in observations of the diffuse ISM at high Galactic latitudes. These features are listed below.

• The velocity dispersion (σ_1pc = 0.7–2.3 km s⁻¹) and CNM fractions (f_CNM = 17–46%) are representative of what is measured at 21 cm in the diffuse ISM of the solar neighborhood (Marchal & Miville-Deschênes 2021; Marchal et al. 2024).

• The range of F(ϕ) values is comparable to the ones measured with LOFAR in the 3C196 field with an average Faraday depth of −2 < M₁ < 3 rad m⁻² and spread of the Faraday spectra of M₁ ∼ 1–2.5 rad m⁻².

• The simulations of our experiment produce a range of peak polarization intensity between 0 and 12 K (Figures 10 and D.1). The highest values might not be representative as they were obtained when B is exactly perpendicular to the LOS, but our experiment indicates that such relatively low H i column density (⟨N_HI⟩ = 1.55 × 10²⁰ cm⁻²) can easily produce P_max values of a few kelvin, similarly to what is observed in the 3C196 field (P_max ≤ 7 K Jelić et al. 2015).

• The synthetic Faraday cubes produced, including the effect of a proxy of the LoTSS uv minimum, show elongated structures with no polarization, reminiscent of the observed depolarization canals (e.g., Erceg et al. 2022).

• The HOG analysis used to evaluate the correlation between H i and Faraday structures shows a significant correlation of the CNM features with F(ϕ), as revealed in the 3C196 field (Bracco et al. 2020).

Our results are indicative that the physics implemented in such numerical experiment is capturing fundamental aspects of the physics of the diffuse ISM. Even though these experiments were designed to reproduce the diffuse neutral gas (H i), this study indicates that a simple prescription for computing the ionization fraction (and therefore n_e) reproduces several properties of the LOFAR data, including a correlation between CNM structures and the morphology of F(ϕ) maps. We recall that before the results of Van Eck et al. (2017), Jelić et al. (2018), and Bracco et al. (2020), the RM values and Faraday structures observed in low-frequency radio data were thought to originate mostly from thermal electrons located in the warm ionized medium (WIM). Our numerical study indicates that even witH its low ionization fraction (∼10⁻² in the WNM and ∼10⁻⁴ in the CNM), the H i gas contributes significantly to the Faraday rotation signal, at levels strong enough to capture a correlation between F(ϕ) and the CNM 21 cm emission. This is compatible with the results of Boulanger et al. (2024) which, through a totally different approach based on UV-absorption data, also concluded that electrons in the WNM are a major contributor to the Faraday depth.

5.2 CNM-Faraday tomographic data correlation

Even though the numerical experiment presented here reveals a significant correlation between CNM and Faraday structures, when we looked closely we could not find a single example in our simulation set that would perfectly reproduce the relative correlation of CNM and WNM with the Faraday tomographic data (i.e., η).

The Faraday depth values observed in the 3C196 field are almost centered on zero (see M₁ in Fig. 9, left). As mentioned in Jelić et al. (2018), this implies that the magnetic field is likely close to perpendicular to the LOS in this area of the sky. Indeed, we reproduced this feature (⟨M₁⟩ ∼ 0) via synthetic observation with B perpendicular to the LOS, but the distribution of M₁ and the range of Faraday spread (M₂ ∼ 1–2.5 rad m⁻²) for this configuration are narrower than what is observed, even after de-biasing for noise.

Another important limitation concerns the physical size of the simulation box. The comparison between the simulated and observed Faraday depth distributions implicitly assumes that all polarized emission and rotation detected by LOFAR originate within a volume comparable to that of the simulation (50 pc on one side). In reality, the observed emission in the 3C196 field likely arises from a more extended region along the LOS, possibly spanning several hundred parsecs through the local magneto-ionic medium. This difference in physical depth may partly account for the narrower range of simulated Faraday depths and could also influence other derived quantities such as rotation-measure dispersion and the H i–polarization correlations. Nevertheless, achieving realistic CNM structures requires sub-parsec resolution, which cannot be maintained in simulations spanning Galactic-scale path lengths.

Regarding the H i-Faraday correlation, the B ⊥ LOS case is the one for which the values of η depart the most from the observations. As shown in Figures 13 and 16, the CNM shows significantly smaller values of $\tilde{V}$ than the WNM (i.e., higher values of η) than in the case where B is parallel to the LOS. The only way to obtain η values close to the ones of the 3C196 field is by adding unrealistically large noise in the 21 cm synthetic data – i.e., ten times the one present in the real EBHIS data – in order to significantly affect the 21 cm WNM signal, but not as much the CNM. The case with B ∥ LOS produces values of η that are closer to the observed ones, but this configuration shows values of M₁ that are significantly larger than the observations. A larger CNM-F(ϕ) correlation (i.e., lower values of η) could be produced by an increase of n_e in the dense H i, producing a stronger Faraday signature associated with CNM structures. This could be adjusted by modifying the post-processing prescription used to estimate n_e locally. Another way would be a stronger morphological correlation between the CNM filamentary structures and the orientation of B. We noticed in our simulations that many CNM clouds seem to have a random orientation with respect to B (see Fig. 2). This might require another forcing scheme for the simulations, perhaps using stellar feedback (e.g. supernovae) that might introduce coherent high pressure fronts over larger scales.

Finally, another way to lower η would be to have a lower correlation between the WNM and the Faraday structure. This could be caused by the WIM itself, which we did not model in our setup. The WIM is likely to have a diffuse structure, similar to the WNM, with more large-scale fluctuations. If the WIM is not spatially correlated to the WNM, its contribution to the Faraday deptH is likely to affect the WNM-Faraday correlation more than the CNM one. The same is true if one considers the fact that the LOS in the 3C196 field is likely to be about four times larger than what is modeled here. Adding up WNM signal along the LOS could lower the morphological correlation between the WNM and Faraday. As the CNM occupies smaller volumes along the LOS, it is less affected by the projection effect of long LOSs. All these effects are beyond the scope of this work, but they could be tested with dedicated numerical setups. Finally, our comparison relied entirely on a single LOFAR field, and thus it sampled a limited volume of the local ISM. The magnetic and thermodynamic conditions along this LOS may not be representative of other regions of the sky. Consequently, any correspondence between simulations and observations should be interpreted as illustrative rather than universally representative of the Galactic ISM that could be tested with the whole LoTSS mosaic (e.g., Erceg et al. 2022; Erceg et al. 2024).

5.3 Depolarization canals

A second aspect that is not well reproduced by our numerical setup is that of depolarization canals. We observe them in our simulations, and some seem to be caused by beam depolarization. On the other hand, the morphology of the synthetic depolarization canals do not match the observations. Depolarization canals observed in the 3C196 field and over large areas at high Galactic latitude (Erceg et al. 2022) are often very elongated, with axis ratios from ten to 100. The depolarization canals produced in our instrumental experiment are more twisty and relatively short.

This calls for another numerical study with the properties of the magnetic field (we used a single value for the initial B) and those of the mechanism injecting kinetic energy in the system varied more. One might want to explore supernova feedback or, alternatively, Fourier driving with different timescales and/or with a more consistent orientation. In general, turbulence-in-a-box simulations have a hard time producing very elongated CNM features, even in cases where they are found to be aligned with B (Inoue & Inutsuka 2016).

6 Conclusion

The main goal of this study was to evaluate if current state-of-the-art MHD numerical experiments used to model the thermally bistable H i are able to reproduce properties of the Faraday sky and, in particular, the H i-Faraday correlation observed in the LOFAR data (e.g., Bracco et al. 2020). Our parametric study is based on 50 pc simulation boxes (MHD + cooling) with a range of turbulent-velocity forcing with different amplitudes and levels of compressibility. The forcing led to realistic ranges of velocity dispersion, CNM fraction, and magnetic-field fluctuations, reproducing values observed in the diffuse ISM. From this set, we produced synthetic 21 cm and synchrotron observations, including the effect of noise and the LoTSS uv minimum, exploring a range of orientations of the LOS from parallel to perpendicular to the mean B field.

Our study reveals that thermal electrons associated with the H i phase could contribute a significant fraction of the low-frequency (100–200 MHz) synchrotron polarized emission observed at high Galactic latitudes. Indeed, our synthetic Fara-day observations reproduce levels of polarization intensity and RM values that are commensurate with what was observed with LOFAR in the 3C196 field, our test region.

To study the contribution of CNM structures to the Faraday tomographic structures in more detail, we developed a metric, η, based on the HOG algorithm (Soler et al. 2019) that measures the ratio of WNM to CNM spatial correlation between 21 cm velocity channels and polarized intensity at specific Faraday depths. We show that this metric is quite robust with respect to the kernel size of the HOG algorithm. Our analysis reveals a significant correlation of the CNM with Faraday structures, similar to what is observed in the 3C196 field (Bracco et al. 2020). The level of CNM-F(ϕ) correlation specifically depends on several factors, which we list below:

• Simulations with low turbulent motions tend to have a stronger H i-Faraday correlation in general, and also smaller η;

• The CNM-F(ϕ) correlation level depends on the B − LOS orientation; B perpendicular (parallel) to the LOS shows the lowest (highest) CNM-F(ϕ) correlation;

• A larger CNM-F(ϕ) correlation can be measured in faint and noisy 21 cm data as higher 21 cm noise tends to affect the broad and low-brightness WNM more compared to the sharper and brighter CNM components;

• Noise in Faraday data reduces the level of correlation of all H i phases equivalently, while not affecting η;

• Finally, the level of CNM-F(ϕ) correlation does not depend strongly on f_CNM. We observe a very mild increase in CNM-F(ϕ) correlation with f_CNM, and only for f_CNM < 20% and in the case where the LOS is perpendicular to the mean B field.

Within the space we explored, the CNM is always correlated with the Faraday tomographic structures, but we were not able to precisely reproduce values of η as low as those found in the 3C196 field (η = 0.72). This indicates that CNM structures are more correlated to Faraday structures in the observations than what we were able to reproduce with our numerical experiment. This is specifically true if we consider simulations with physical conditions similar to the 3C196 field (B ⊥ LOS, σ_u,1pc = 2.0 km s⁻¹, and f_CNM = 38%). In these conditions, we find η ∼ 1.5, whicH indicates that WNM contributes more to Faraday structures than the CNM. Several elements could contribute to this result: the choice of the n_e prescription, the contribution of the WIM to the Faraday depths (not modeled), the longer length of the LOS in the data compared to the simulations, and the Fourier driving that tends to produce CNM structures not aligned with B.

Our results are indicative that low-frequency observations could serve as a probe of the magnetic-field intensity and morphology in the neutral edges of the Local Bubble. Moreover, the detailed comparison of simulations and observations in the 21 cm Faraday space revealed subtle discrepancies that call for an improvement of the way the multiphase, turbulent, magnetized, and partly ionized ISM is modeled today.

Finally, we note that our conclusions are drawn from the comparison with a single observational field. Future work should extend this analysis to additional LOFAR fields to test the robustness and generality of the trends identified here across different Galactic environments.

References

---

Appendix A MOOSE

Mock Observation of Synchrotron Emission (

MOOSE

) is an interactive Julia-based tool designed to process numerical simulations of synchrotron emission and Faraday rotation. It generates synthetic observations by integrating three-dimensional cubes of physical parameters, including density, temperature, and magnetic-field components. The electron density is computed using either empirical prescriptions or direct scaling with the neutral hydrogen density.

The code reconstructs polarization properties by computing synchrotron emission and applying Faraday rotation effects along the LOS. It accounts for observational constraints by filtering large-scale structures to matcH interferometric data, convolving signals with telescope beam responses, and incorporating noise models to replicate real observations.

MOOSE

employs 2D interpolation methods to efficiently estimate synchrotron emissivities based on precomputed lookup tables. The interpolation grid is currently designed to use LOFAR frequencies but can be adapted to any frequency range depending on the observational requirements. This technique allows for accurate integration of synchrotron radiation over a range of magnetic field strengths and frequencies while maintaining computational efficiency.

The code first reconstructs the Stokes parameters (I_ν, Q_ν, U_ν) as a function of frequency before applying

RM-Synthesis

to derive the Faraday spectrum. This approach enables the decomposition of polarized emission along different Faraday depths, making it possible to analyze the magneto-ionic properties of the ISM.

The tool assumes an optically thin medium, meaning absorption effects are neglected. While this assumption is generally valid for high Galactic latitude regions, it may introduce biases in denser environments.

The synthetic data can be analyzed to extract rotation measure maps, Faraday spectra, and polarization structures.

Figure A.1 represents an example of the polarized intensity as a function of the Faraday depth (in rad m²). Peaks in the spectrum indicate distinct polarized emission components along the LOS. The two most significant peaks in the spectrum correspond to Faraday depths of approximately ϕ₁ = 2.75 rad m⁻² and ϕ₂ = 5 rad m⁻², which likely indicate separate emitting regions along the LOS.

It provides a visual representation of how

MOOSE

reconstructs the Faraday spectrum from input simulation data. The presence of sidelobes in the spectrum is a common feature resulting from the limited frequency coverage and resolution in

RM-Synthesis

. These sidelobes can introduce ambiguities in identifying true polarized components.

Appendix B Faraday moments

M₀ is the polarized intensity integrated over the Faraday depth range. M₁ is the polarized intensity weighted mean of Faraday depth ϕ. M₂ is the width of the Faraday spectrum. These are analogous to moments in velocity space and are fully described in Dickey et al. (2019). These quantities are affected by noise; moments are often estimated by thresholding the data, for instance Erceg et al. (2022) used only values six times higher than the noise level. Here we propose alternatives to estimate the Faraday moments without rejecting any data.

Because F(ϕ) is a strictly positive quantity, noise in F(ϕ) has a specific effect. Noise fluctuations do not cancel out in calculating M₀ but they add up introducing a bias. For a given noise level σ_noise (computed here by taking the standard deviation at the first 30 and the last 30 Faraday depths ϕ) the noise bias on the first moment is equal to N σ_noise where N is the number of Faraday depths. To take this into account we compute the first moment as $$M_0={\displaystyle \sum }_{i=1}^N\left(F\left({\varphi }_i\right)-{\sigma }_{\text{noise}}\right)d\varphi$$(B.1)

where dϕ is the channel widtH in ϕ.

The second-order moment M₂ is also affected by noise but it is also very sensitive to the presence of sidelobes in Faraday spectra (see an example in Fig. A.1). These sidelobes arise from the limited range of λ² over which the Fourier transform is computed (see Eq. (12)). The spectral shape of F(ϕ) also depends significantly on the spectral shape of Q_ν and U_ν that can vary significantly due to Faraday rotation (Equations (7) and (8)). This results in spectral structures in P(λ²) that can vary significantly from one position to the next. These spectral features in P(λ²) are not instrumental as they depend on the physical distribution of emitting and rotating structures along the LOS (Takahashi 2023). In certain conditions the combination of Q_ν and U_ν can even cancels signal in P(λ²). As a consequence, the width of the Faraday spectrum can, in some cases, be smaller than the width of the RMSF, making its interpretation difficult.

As an alternative, we estimated the width of the Faraday spectrum using the effective width defined as $$W_{\text{eff}}=\frac{M_0}{\underset{\varphi }{\max }F(\varphi )}.$$(B.2)

This quantity reduces the influence of sidelobes by emphasizing the integrated power of the spectrum relative to its peak intensity. We found that W_eff provides a more robust and interpretable characterization of the spectral widtH in the presence of complex interference patterns. It is also not biased by noise like a standard second moment is. We computed W_eff using M₀ corrected for the noise bias. Then the second order moment is estimated as $$M_2=\frac{W_{\text{eff}}}{\sqrt{2\pi }}$$(B.3)

whicH is the exact solution for a Gaussian function.

Finally, the first moment M₁ defined as $$M_1=\frac{{{\displaystyle \sum }}_{i=1}^N{\varphi }_{i}F\left({\varphi }_i\right)}{M_0}$$(B.4)

is also biased in the presence of noise. As noise increases, the numerator tends to produce M₁ → ⟨ϕ_i⟩. If the ϕ vector is centered on 0, like it is the case in our study, lower S/N F(ϕ) will tend to have M₁ → 0. To circumvent this effect, we are using a matched filter approach. We leverage on the fact that W_eff is a non biased estimate of the width of the signal, especially if the complexity of the Faraday spectra is small (i.e., F(ϕ) is characterized by one dominating component, large or narrow). The first moment is found by convolving F(ϕ) by a Gaussian filter with a width given by W_eff. M₁ is then estimated at the value of ϕ where this convolution is maximum.

Appendix C HOG results for the whole set of simulations

The results of the HOG analysis for the whole 14 simulations set are presented in Fig. C.1 (B parallel to the LOS) and Fig. C.2 (B perpendicular to the LOS). Different color backgrounds are used to illustrate low (v_RMS < 8 km s⁻¹ – gray background) and high (v_RMS > 8 km s⁻¹ – light blue background) turbulence level.

Appendix D Rotation of the cubes

To analyze physical quantities along different lines of sight, it is necessary to rotate the simulation cubes containing both scalar and vector fields. The rotation is particularly important in the context of these simulations because the magnetic field has a preferred orientation with a strong mean component in a specific direction. In our simulations, the ratio of the mean field strength to the root-mean-square fluctuation (B/B_rms) is large, meaning that the field is strongly ordered. As a result, analyzing the simulations only along the three cartesian orientations provides a limited perspective on the system.

By rotating the cubes by arbitrary angles, it becomes possible to explore configurations that are not perfectly aligned with the dominant field direction. This is crucial for understanding how turbulence and ordered structures interact under different projection effects. It also allows for a more comprehensive analysis of how observables change with different viewing angles.

Rotation of scalar fields

For scalar fields such as electron density (n_e), temperature (T), and neutral hydrogen density (n_H), the rotation consists of a coordinate transformation from the original frame (x, y, z) to the new frame (x′, y′, z′). Given a rotation matrix R, the new coordinates are obtained as: $$\left[\begin{array}{l}{x}^{\prime }\hfill \\ y^{\prime }\hfill \\ z^{\prime }\hfill \end{array}\right]=R\left[\begin{array}{l}x\hfill \\ y\hfill \\ z\hfill \end{array}\right].$$(D.1)

The scalar field values remain unchanged at each transformed coordinate. In practice, the rotation is implemented by interpolating the values of the scalar field onto the rotated grid using bilinear interpolation.

Rotation of vector fields

For vector fields such as the magnetic field components (B_x, B_y, B_z) and velocity components (V_x, V_y, V_z), both the coordinate system and the vector components must be rotated. If the original vector field is represented as: $$U=\left[\begin{array}{l}{U}_x\hfill \\ U_y\hfill \\ U_z\hfill \end{array}\right],$$(D.2)

then after rotation, the new components are given by: $$U^{\prime }=RU,$$(D.3)

where R is the same 3 × 3 orthogonal rotation matrix applied to the coordinate system. However, to correctly rotate a vector field, it is necessary to first apply the inverse transformation to determine where in the original coordinate system the vector should be sampled, and then apply the rotation to the vector itself: $$U^{\prime }\left(r^{\prime }\right)=RU\left(R^{-1}{r}^{\prime }\right).$$(D.4)

These transformations ensure that the rotated vector components remain physically meaningful and aligned with the new coordinate system.

Figure D.1 shows maps of P_max and column density of the case-study simulation rotated around the z-axis by a 45^◦. Fig. D.2 shows the result of HOG in this configuration.

All Tables

All Figures