© The Authors 2024

1 Introduction

Ultrahigh-energy cosmic rays (UHECRs) are cosmic rays with energies surpassing 10¹⁸ eV. Since they are charged, the Galactic magnetic field (GMF) deflects their paths, and this leads to a mismatch between the true and the observed arrival direction. It is a crucial challenge in the field of high-energy astrophysics to determine their origins. Successfully addressing this challenge could offer insights into astrophysical processes that generate UHECRs, and into their composition. Additionally, knowledge of UHECR sources would be a crucial ingredient in multimes-senger studies of high-energy systems (e.g., Fang & Murase 2018; Murase 2019).

Although numerous theoretical models have been proposed to explain the sources of UHECRs (e.g., Bhattacharjee & Sigl 2000; Torres & Anchordoqui 2004; Kotera & Olinto 2011), it has proven to be a complicated task to determine these sources concisely. The main challenge arises from the fact that UHECRs are charged particles, and they are deflected by both the GMF and the intergalactic magnetic field. As a result, even if multiple UHECRs were emitted from a single, intense, and proximate cosmic-ray source (di Matteo et al. 2023), their trajectories would be dispersed across the plane of the sky (POS). Consequently, any UHECR hotspot would not align with the source, but would be displaced away from it due to systematic deflections by the ordered component of the GMF, in addition to being spread out due to the random deflections by the turbulent component of the GMF. This situation is different from that of photons or neutrinos, for which it is more straightforward to establish a connection between observed events and their probable sources, even in the limit of low statistics and the poor angular resolution of their detectors. In addition to UHECRs, a better understanding of the 3D structure of the GMF would significantly contribute to our understanding of a multitude of astrophysical problems, such as endeavors to subtract the Galactic foreground from the cosmic microwave background foreground (Tassis & Pavlidou 2015; Rubiño-Martín et al. 2023), as well as the physics of star formation (Pattle et al. 2023; Doi et al. 2024).

The primary challenge in understanding the GMF lies in the difficulty of obtaining a 3D tomographic reconstruction of the intervening GMF because the majority of the currently accessible observations are integrated along the LOS, such as observations of the Faraday rotation induced by the GMF (Pandhi et al. 2022; Hutschenreuter et al. 2024). This limitation has guided the predominant approach in GMF modeling to rely on parametric models. This is typically achieved by fitting parameters to distinct analytic components, for instance, a toroidal component, a poloidal component, and a turbulent component. To model the latter, a Gaussian random field is employed (Sun et al. 2008; Sun & Reich 2010; Takami & Sato 2010; Jansson & Farrar 2012a; Jansson & Farrar 2012b).

However, direct insights into the 3D structure of the interstellar medium of the Milky Way are attainable. The Gaia mission has mapped the positions of over a billion stars in the Galaxy by accurately measuring stellar parallaxes (Gaia Collaboration 2016, 2021; Bailer-Jones et al. 2021). Combined with other spectroscopic data, these parallaxes have enabled the construction of 3D tomographic maps showing the dust density distribution in certain regions of the Galaxy (Lallement et al. 2018, 2019, 2022; Green et al. 2019; Leike & Enßlin 2019; Leike et al. 2020, 2022; Edenhofer et al. 2024b). Nevertheless, these are maps of the dust density and do not directly constrain the magnetic field.

Observational methods that probe the 3D structure of the GMF do exist, however. A notable example is the linear polarization of starlight. Typically, starlight originates from its source as unpolarized light, but it can become linearly polarized due to the dichroic absorption by interstellar dust particles, which align themselves with the surrounding magnetic field (Andersson et al. 2015).

Future optopolarimetric surveys such as the Polar-Areas Stellar Imaging in Polarization High-Accuracy Experiment (PASIPHAE) and SOUTH POL are poised to deliver high-quality stellar polarization measurements for millions of stars (Magalhães 2012; Tassis et al. 2018; Maharana et al. 2021; Maharana et al. 2022). When combined with the stellar distance data obtained from the Gaia survey, these measurements will enable direct tomographic measurements of the GMF POS component in regions with dust clouds (Davis 1951; Chandrasekhar & Fermi 1953; Panopoulou et al. 2017; Skalidis et al. 2021; Skalidis & Tassis 2021; Pelgrims et al. 2023). Notably, optical stellar polarimetry was recently used in order to probe the GMF in 3D in the Sagittarius spiral arm (Doi et al. 2021; Doi et al. 2024). Additionally, local information can be further constrained through the study of HI gas in different velocity bins, which provide GMF information on the number of clouds along a given LOS as well as local information through the use of 3D dust-reddening maps (Green et al. 2018; Tritsis et al. 2018, 2019; Clark & Hensley 2019). In conjunction with available LOS data (see, e.g., Tahani et al. 2022a,b), this information promises to provide localized and sparse GMF data in the future. All of these techniques will provide crucial information for creating 3D tomographic maps of specific areas of interest. With these maps, the paths of UHECRs can be traced back through these regions and the source localization on the sky can be improved (however, the contribution of the intergalactic magnetic field is still not accounted for). Specifically, there is an intense interest in mapping the GMF in the direction of UHECR hotspots, as well as in parts of the Galaxy that are likely to have been traversed by particles comprising these hotspots (Abbasi et al. 2014; Pierre Auger Collaboration 2017; Kawata et al. 2019).

This study is the second entry in our effort to reconstruct the GMF nonparametrically in 3D in a Bayesian setting. It directly follows Tsouros et al. 2024, hereafter Paper I. Essentially, we address an inverse problem within a Bayesian framework, where the goal is to sample the posterior distribution of GMF configurations in a specific part of the Galaxy using a combination of local and LOS-integrated information. In this work, local measurements only provide information for the POS component of the magnetic field. This corresponds to the information content of tomographic measurements of interstellar magnetized dust through optopolarimetry of starlight. On the other hand, LOS-integrated measurements provide information about the LOS component of the magnetic field as derived from Faraday rotation measurements, for instance. We approach this problem within the context of information field theory, which was developed specifically for Bayesian inference for fields and was successfully applied in various contexts (Enßlin et al. 2009; Enßlin 2019; Enßlin 2022). By reconstructing the posterior distribution of GMF realizations, we aim to accurately recover the true arrival directions of UHECRs given the observed arrival directions. We also account for the influence of the GMF.

In Section 2 we briefly describe the method, the forward models we used, and the sampling of the posterior. In Section 3 we present the main results of the algorithm for the considered scenarios, and in Section 4 we discuss the results further.

2 Method

In general, we are interested in inferring the configuration of the GMF, B(x) with x ∈ $\mathcal{V}$ over a domain $\mathcal{V}$ ⊂ ℝ³, given some observed dataset d. In the context of Bayesian inference for continuous signals, the task is to determine the posterior probability distribution of B(x) for d, $$P(\mathit{B}\mid d)=\frac{1}{Z}P(d\mid \mathit{B})P(\mathit{B}).$$(1)

Here, P(d|B) is the likelihood, representing the probability of observing magnetic field measurements d given a specific configuration B(x). The prior, P(B), encapsulates preexisting information about B(x), and Z = P(d) is the normalization factor.

In this work, the field that serves as a ground truth (the true field) was generated from a dynamo magnetohydrodynamic simulation discussed in Appendix A. The original simulation domain extended to ~ 1 kpc in the x-y direction and to ~2 kpc above the Galactic plane. The GMF was rescaled so that its root mean square (RMS) value was 5μG.

2.1 Likelihood

Tomography of the magnetized ISM from stellar polarization measurements is a highly nontrivial problem, and its full discussion is beyond the scope of this work (Pelgrims et al. 2023). However, we note that by combining Gaia data and stellar polarization data for stars with a known distance to the Sun, it is possible to acquire information on the Stokes parameters that each intervening dust cloud imposes on the observed starlight if enough stars have polarization measurements and known distances. This can then be translated into local information on the orientation of the POS component of the GMF at that cloud through the connection to grain alignment, as referenced briefly in the previous section and thoroughly used in Tassis et al. (2018). Information on the POS component of GMF in clouds can also be acquired by the use of 21 cm neutral hydrogen (HI) emission measurements (Clark & Hensley 2019). We assumed that the task of determining the locations to which the measurements correspond to has been carried out.

Thus, for the ith datapoint, we assumed a forward model of the form $${\mathbf{d}}_{\text{local }}^{(i)}=\int R_{\text{local }}(\mathbf{x},{\mathbf{x}}_i)\mathit{B}(\mathbf{x})d^{3}x+{\mathbf{n}}_{\text{local }}^{(i)},$$(2) $$R_{\text{local }}(\mathbf{x},{\mathbf{x}}_i)\equiv {\delta }^{(3)}(\mathbf{x}-{\mathbf{x}}_i)P_{\mathrm{P}\mathrm{O}\mathrm{S}},$$(3)

where B(x) is the magnetic field, and ${\mathbf{n}}_{\text{local }}^{(i)}$ are the observational uncertainties that contaminate our measurements. For clarity in the notation, we henceforth suppress the volume integration and write d_local − R_localB, where the volume integral is implicit. The vector x_i is the location of the ith cloud where the magnetic field is measured, and P_POS signifies a projection operator on the POS, which reflects that (mainly) the POS component of the magnetic field is measured via dust polarization, $$P_{\mathrm{P}\mathrm{O}\mathrm{S},ij}={\delta }_{ij}-{\hat{x}}_i{\hat{x}}_i^{\mathrm{T}}$$ with ${\hat{x}}_i=x_i/\parallel x_i\parallel$ (assuming the observer to be at the origin). The Dirac delta function localizes the measurements at specific known locations x_i.

The option to include the operator P_POS into the considered scenario is central to this work because it presents one of the main additions compared to Paper I. A complete projection on the POS is a pessimistic scenario because LOS information can become available by incorporating Zeeman or Faraday rotation data (Tahani et al. 2022a; Tahani et al. 2022b). A complete projection on the POS should therefore be seen as an extreme benchmarking scenario.

We note that this forward model is quite simplistic in that it assumes that accurate 3D locations are measured. Formally, this is captured by the Dirac delta function and that the locations x_i are to be assumed known up to the resolution length scale. As we show in Section 2.4, the resolution of the grid in which our reconstruction takes place is some tens of parsecs, which is comparable to the uncertainty of the cloud localization (Pelgrims et al. 2023). In this particular application, distance uncertainties are therefore not relevant because they tend to be on scales that are lower than the grid resolution. However, it should be noted that in applications of this method with grids whose resolution is below the scale of the cloud localization uncertainties, the latter should be included because they would indeed have an observable effect.

The vector ${\mathbf{n}}_{\text{local }}^{(i)}$ is assumed to be a random variable drawn from a Gaussian distribution with a known covariance N_local. When specific measurement techniques are identified, other more appropriate error distributions will be chosen. When we marginalize over the noise, the likelihood becomes $$P(\mathit{d}\mid \mathbf{B})=\mathcal{G}({\mathit{d}}_{\text{local }}-R_{\text{local }}\mathit{B},N_{\text{local }}).$$(4)

The covariance N_local was chosen to be a multiple of the identity, (N_local)_ij = σ⁻²δ_ij, where we chose $$\sigma =\frac{|\mathbf{B}{|}_{\mathrm{R}\mathrm{M}\mathrm{S}}}{2},$$(5)

where |B|_RMS = 5 μG is the RMS value of the magnitude of the ground truth. This does not imply that the noise is correlated with the GMF covariance. It was merely chosen as such in order to ensure an S/N of about 2. This S/N was chosen as a worst-case scenario for the expected GMF measurement error that arises from the measurement of the magnitude of the POS component of the GMF, which is expected to have an S/N lower than 2 (Skalidis et al. 2021). As the directional error is subdominant to that of the POS magnitude estimate (Pelgrims et al. 2024), we assumed that it was included in our choice of the worst-case scenario S/N.

In addition to local data, we explored the possible use of integrated LOS data as inferred from Faraday measurements, for instance (Hutschenreuter et al. 2024). In this case, the forward model takes the form $$d_{\mathrm{i}\mathrm{n}\mathrm{t}}^{(i)}={\left(\overline{P_{\mathrm{L}\mathrm{O}\mathrm{S}}\mathit{B}}\right)}_{L_i}+n_{\mathrm{i}\mathrm{n}\mathrm{t}}^{(i)},$$(6) $${\left(\overline{P_{\mathrm{L}\mathrm{O}\mathrm{S}}\mathit{B}}\right)}_{L_i}\equiv \frac{1}{\left|L_i\right|}{\int }_0^{\mid L_i\mid }{B}_{\parallel }(\mathit{x})d\ell ,$$(7)

where P_LOS projects a vector onto the LOS component (B_‖), and L_i is the specific LOS under consideration. Furthermore, |L_i| denotes the limit up to which we integrate. In this application, |L_i| coincides with the distance between the Earth and the intersection of L_i with the boundary of $\mathcal{V}$. Essentially, the above is equivalent to assuming that the electron density is roughly constant and known up to |L_i| and then falls to zero. While this is not a valid assumption for low Galactic latitudes, we maintained it in this proof-of-concept work. Finally, the vector $$n_{\text{int }}^{(i)}$$ corresponds to a random vector on the POS, with a covariance N_int.

The likelihood for the combined data is given by $$P(\mathit{d}\mid \mathbf{B})=\mathcal{G}({\mathit{d}}_{\text{local }}-R_{\text{local }}\mathit{B},N_{\text{local }})\mathcal{G}(d_{\text{int }}-{\left(\overline{P_{\text{LOS }}\mathit{B}}\right)}_{L_i},N_{\text{int }}).$$(8)

Similarly, we defined the covariance for the noise of the integrated measurements as $${\left(N_{\text{int }}\right)}_{ij}={\sigma }_{\text{int }}^2{\delta }_{ij}$$, where¹ $${\sigma }_{\text{int }}=\frac{1}{2}\text{μ}\mathrm{G}.$$(9)

Finally, the operator R_local, which sparsely samples the GMF, is defined as follows. After discretizing our domain to voxels (see Section 4.1), we applied a Bernoulli trial to each voxel to determine whether it was observed with a probability p and 1 − p. The probability p was chosen so that the approximate number density of measurements n_data was about $$n_{\text{data }}\sim \{\begin{array}{ll}100(\mathrm{k}\mathrm{p}\mathrm{c}{)}^{-3},& \text{ if }T\ge {10}^4\mathrm{K}\\ 1000(\mathrm{k}\mathrm{p}\mathrm{c}{)}^{-3},& \text{ if }T<{10}^4\mathrm{K},\end{array}$$(10)

where T is the corresponding gas temperature of that voxel, acquired from the same simulation that produced our ground truth. This choice of n_data reflects the decay in the number of dust clouds as a function of distance from the Galactic plane, which directly correlates with the expected number of measurements with respect to the position above the Galactic plane because the local measurements of the GMF will ultimately exist where dust clouds are located after polarized-starlight tomography has been carried out. Moreover, as most of the gas has a temperature above 10⁴ K, the number density of local measurements is 262 (kpc)⁻³. This is about a few hundred parsec per cloud in the cold medium, which is comparable to the expectation for the true ISM (Panopoulou & Lenz 2020).

2.2 Prior

As in Paper I, the only hard constraint that needed to be imposed is that regardless of the candidate magnetic field configuration B we considered, it must satisfy ∇ · B = 0 in order to be a viable candidate. To ensure that the magnetic field was divergence free, we assumed that it was related to a non-divergence-free random field φ by a divergence-cleaning operator $\mathcal{P}$. This transverse projection operator, defined in Fourier space as $${\mathcal{P}}_{ij}(\mathbf{k})={\delta }_{ij}-{\hat{k}}_i{\hat{k}}_j^{\mathrm{T}},$$(11)

projects out the degrees of freedom of the Gaussian random vector field that violate the divergence-free condition. In other words, it connects a latent field φ(x) to the true magnetic field by the harmonic space relation $${\hat{B}}_i(\mathbf{k})=\frac{3}{2}{\mathcal{P}}_{ij}(\mathbf{k}){\hat{\phi }}_j(\mathbf{k}),$$(12)

where k are wavevectors, and the Fourier modes are denoted by hatted fields. Eq. (12) ensures that ∇ · B = 0, while the factor 3/2 accounts for power loss due to reduced degrees of freedom, aligned with the original assumption of statistical isotropy for φ (Jaffe et al. 2012), and the hatted fields denote the respective Fourier transforms. Our aim is to reconstruct the local GMF B by inferring the latent field φ, which is related to the latter by Eq. (12). For φ, we assumed a Gaussian prior of the form $$\mathcal{P}(\mathit{\phi })=\frac{1}{|2\pi \mathrm{\Phi }{|}^{\frac{1}{2}}}\exp [-\frac{1}{2}\int d^{3}x{d}^3{x}^{\mathrm{\prime }}{\displaystyle \sum _{ij}}{\phi }_i(\mathbf{x}){\mathrm{\Phi }}_{ij}^{-1}(\mathbf{x},{\mathbf{x}}^{\mathrm{\prime }}){\phi }_j\left({\mathbf{x}}^{\mathrm{\prime }}\right)].$$(13)

The quantity Φ_ij is the covariance matrix, defined as $${\mathrm{\Phi }}_{ij}(\mathbf{x},{\mathbf{x}}^{\mathrm{\prime }})=⟨{\phi }_i(\mathbf{x}){\phi }_j^{\ast }\left({\mathbf{x}}^{\mathrm{\prime }}\right)⟩,$$(14)

where the symbol ⟨...⟩ signifies an average over the distribution P(φ). That is, if $\mathcal{O}(\mathbf{x})$ is some quantity of interest, then $$⟨\mathcal{O}(\mathbf{x})⟩\equiv \int d\mathit{\phi }P(\mathit{\phi })O(\mathbf{x}).$$

We note that the average was taken over field configurations.

In our analysis, we chose not to integrate any prior knowledge about the GMF geometry and statistics, and we therefore used a prior distribution exhibiting statistical isotropy, homogeneity, and mirror symmetry². This is formally encapsulated by writing the Fourier space covariance in the form $$⟨{\hat{\phi }}_i(\mathbf{k}){\hat{\phi }}_j^{\ast }\left({\mathbf{k}}^{\mathrm{\prime }}\right)⟩=(2\pi {)}^3{\delta }_{ij}{\delta }^{(3)}(\mathbf{k}-{\mathbf{k}}^{\mathrm{\prime }})P(k).$$(15)

A crucial point is that the 3D prior power spectrum P(k) is not known and is to be inferred as well. It is modeled as a sum of a power law and an integrated Wiener component (Arras et al. 2022). The defining hyperparameters and their prior PDFs (typically called hyperpriors) are summarized in Table 1, and they are also briefly discussed in Paper I.

2.3 Sampling the posterior

Equipped with the likelihood and prior, the posterior in terms of the magnetic field B is given by Eq. (1). Because the power spectrum P(k) needs to be inferred along with the configuration of the GMF, this inference problem is nonlinear and cannot be solved by a generalized Wiener filter (Pratt 1972). For this reason, a nonperturbative scheme, called geometrical variational inference (geoVI) developed by Frank et al. (2021), was used. A brief exposition on geoVI can be found in Appendix A of Paper I. For the purposes of this work, it suffices to state that we did not sample magnetic field configurations from the true posterior directly, but rather from an approximate posterior, as is usually the strategy in variational methods. For this task, we employed the numerical information field theory (NIFTy³) package in Python (Selig et al. 2013; Steininger et al. 2017; Arras et al. 2019, Edenhofer et al. 2024a). The required input is the likelihood and the prior of the original physical model, as described in Sections 2.1 and 2.2, respectively.

2.4 Procedure

We summarize the specific setting we probed below and describe how we generated the synthetic data on which the method was verified.

• Spatial domain: The modeled space was assumed to be periodic because we implemented the details of the ground truth, where periodic boundary conditions were used in order to emulate the axisymmetry of the Galaxy (see Appendix A, and Bendre & Subramanian 2022 for a broader discussion of the simulation used). Furthermore, we padded our space by a factor of two, so that the x and y directions reached an extent of ~1.6 kpc,. The resulting cube was uniformly partitioned into N_x, N_y, and N_z segments per axis, where N_x = N_y = 48, and N_z = 64, with padding. In this setting, every voxel had a linear dimension of approximately 30 pc. This accommodated the expected size of the dust clouds, as well as the uncertainty of the measurement positions (at least as an order of magnitude), as shown in Pelgrims et al. 2023.

• Data masking: We applied R_local (see Section 2.1) to the ground truth field in order to acquire the noiseless data.

• Adding noise to local data: Gaussian noise with a covariance matrix N_local (Eq. (5)) was added to each observed data vector.

• Integrated data: Optionally (see Section 3.3), the likelihood was supplemented by an additional term for the integrated local measurements, as in Eq. (8). In practice, the magnetic field was transformed from a Cartesian coordinate system to a spherical polar coordinate system with the Earth at the origin. Then, the radial component of the GMF, which is equivalent to the LOS component, was integrated along individual LOSs, resulting in a set of 2D integrated measurements that informed the model further.

• Adding noise to integrated data: Gaussian noise with a covariance N_int (Eq. (9)) was added to each pixel on the celestial sphere to contaminate the data acquired from the previous step.

• Sampling the approximated posterior: Finally, the geoVI method was applied to the true posterior distribution, resulting in samples from the approximate distribution. To all the latent fields we sampled, we again applied the projection operator (Eq. (11)) to obtain posterior samples of the divergence-free GMF.

• Application to UHECR back tracing: Through each of the GMF samples drawn from (1) in the previous step, we traced back an UHECR of known observed arrival direction θ_obs and rigidity r_* ≡ E/Z. A recording of the final velocity of the particles, in particular their original directions θ when they leave $\mathcal{V}$, essentially provides samples from the distribution P(θ|D) of the original arrival directions of the particles before they entered the GMF for the data $$D\equiv \{d,r_{\ast },{\theta }_{\text{obs }}\}.$$(16)

  To keep the discussion simple, we only considered UHECRs with a fixed rigidity r_* = 5 × 10¹⁹ eV (equivalently, protons of energy equal to the GZK limit⁴, E = 5 × 10¹⁹ eV.). As a way to benchmark the quality of our reconstructions in the context of UHECR physics, we compared the angular separation δθ between the true arrival direction θ_true and that of the back-propagated UHECR, ending up with a distribution over δθ. In this context, the true arrival direction always refers to the direction of the UHECR immediately where it entered $\mathcal{V}$. In Fig. 1, we provide a visual representation of the quantities defined in this section.

3 Results

In this section, we use NIFTy in order to sample the posterior distribution for three different scenarios: In scenario A, the observed data only consist of local measurements, and at each location, we only probe the components of the GMF that are parallel to the POS. In scenario B, all three components of the GMF (including the LOS) are probed on an equal footing for comparison. Finally, in scenario C, we use the same dataset as in scenario A, but additionally use integrated LOS information over the whole sky.

For each of these scenarios, we benchmarked the success of the reconstruction by using the GMF posterior samples in order to infer the true arrival direction of a UHECR with fiducial rigidity of r_* = 5 × 10¹⁹ eV for all possible observed arrival directions on the northern sky, as described in the previous section.

Fig. 1 Illustration and definition of quantities used in the analysis. A UHECR with a known rigidity r, enters the Galaxy with an arrival direction θtrue (red dot). Because of the GMF, it is deflected and is observed on Earth as arriving from θobs (black dot). The angular distance between θobs and θtrue is α, and it is the error that the GMF induces on the observed arrival direction. We traced back the particle through each GMF configuration sampled using NIFTy, thus ending up with a distribution of arrival directions P(θ|D), with D defined in Eq. (16). From the posterior samples drawn, we calculated the mean angular distances ⟨δθ⟩θ|D and ⟨δϕ⟩θ|D to the true and observed arrival directions, respectively, as well as the standard deviations for the former. The scales in this artificial example are larger for visual clarity, and do not correspond to an application of the method.

3.1 Scenario A: Local measurements with information on the plane of the sky alone

The local GMF information that can be acquired through starlight polarization-based tomography alone is confined to the celestial sphere (Panopoulou et al. 2019; Pelgrims et al. 2023). In this section, we therefore sample the posterior Eq. (1) for local GMF data d that are completely blind to the LOS dimension, as is the case for polarization measurements.

To do this, we worked on a spherical polar coordinate system with the Sun at the origin. The magnetic field is expressed as B(x) = (B_r, B_θ, B_φ) in this coordinate system. In Fig. 2, we reconstruct the simulated GMF described in Appendix A. In Fig. 2a, we show the ground truth. Fig. 2b depicts the synthetic local GMF data obtained from the ground truth for this scenario. The result of the reconstruction algorithm is a set of 100 posterior samples of Eq. (1), given the data of Fig. 2b. In Fig. 2c, we show the mean of the posterior samples.

In Figs. 3a and 4a, we show the mean and standard deviation of the angular distance error ⟨δθ⟩_θ|D and σ_θ|D respectively) obtained through the use of the GMF reconstructions shown in Fig. 2. ⟨δθ⟩_θ|D and σ_θ|D vary across the celestial sphere, and the specific structure of these functions depends on the specific chosen ground truth GMF. The greatest error of the reconstruction for this setting is approximately 14°. In order to judge the performance, we depict in Fig. 5a the angular error in the arrival direction assuming the observed ones were true, that is, we ignore the correction using the recovered GMF. By comparing Fig. 5a to Fig. 3a, we observe that reconstructing the local GMF for d yields a significant improvement in our ability to recover UHECR arrival directions. This result suggests that ⟨δθ⟩_θ|D is greater for UHECRs that are observed to arrive from directions in which the influence of the GMF is greater (Fig. 5a), in this case, at low longitudes. This correlation is further explored in Section 4.1.

Fig. 2 Reconstruction of the simulated 3D magnetic field with the use of local data that lack LOS field component information. The blue sphere represents the celestial sphere. Panel a: ground truth; the GMF obtained as described in Appendix A. The field is rescaled so that it has an RMS norm of 5 μG. Panel b: synthetic data based on the ground truth of Fig. 2a. The radial component of the magnetic field is not measured. Panel c: the mean of the approximating posterior distribution attained via the geoVI algorithm based on the data provided in Fig. 2b. Panel d: mean of the approximating posterior distribution attained for the local data of Fig. 2b as well as integrated measurements of the radial component (Fig. 6b). The side of the box is about 1.6 kpc, or approximately 20 times the GMF correlation length.

3.2 Scenario B: Local measurements with full 3D information at each measured location

In this section, we examine the impact of a complete lack of observation of the LOS (scenario A) on the UHECR arrival direction reconstruction. For this purpose, we performed the same inference as in Section 3.1, but the LOS component was also probed locally, just like the POS components. In Figs. 3b and 4b, we plot the mean angular error ⟨δθ⟩_θ|D and the respective standard deviation for this scenario. Compared to the results of scenario A (see Figs. 3a and 4a), the quality of the reconstruction greatly improves when local LOS information is included. While the maximum mean angular error drops by a few degrees, the improvement is dramatic in general in that the total area of the sky where the maximum bias occurs is substantially reduced. This observation also holds for the variance.

While we considered θ_obs over the whole northern hemisphere for benchmarking purposes, in real applications, only sufficiently high Galactic latitudes are relevant. The reasons for this are twofold: First, reconstructions of the GMF will mainly be available at high enough latitudes because a tomographic reconstruction is exceptionally difficult on the Galactic disk due to the high density of HI clouds. Second, as Gaia data extend to approximately two kiloparsecs, this also defines the limits of the Galactic region in which our method can be expected to yield a sufficiently good reconstruction. However, UHECRs observed as arriving at low Galactic latitudes will have traveled through a significant part of the Galaxy that is not reconstructed, and this places a limit on the applicability of our method on low enough latitudes.

We have shown that knowledge of local LOS information would yield a substantial improvement over our ability to reconstruct the GMF, at least as far as UHECR backtracing is concerned. As stellar polarization data alone cannot probe the LOS dimension, this information would have to be supplemented by additional methods (e.g., Zeeman measurements). However, the local measurement of the LOS GMF component is a notoriously difficult task. We therefore attempt to mitigate this below by including integrated LOS information in our likelihood.

Fig. 3 Mean angular error of the reconstruction (see Fig. 1) as a function of all possible arrival directions on the northern hemisphere for the case of a UHECR with a rigidity r* = 5 × 1019 eV. Panel a: magnetic field data consisting of local information with the LoS component are projected out (scenario A). Panel b: magnetic field data consisting of local information with the measured LOS component (scenario B) Panel c: as in the top left panel, but the data are supplemented by integrated LOS data (scenario C; see Fig. 6). The red and orange lines on the color bar indicate the maximum and mean values of the map, respectively.

Fig. 4 As in Fig. 3, but for the corresponding angular error standard deviations as a function of the observed arrival direction.

Fig. 5 Amount by which a UHECR of rigidity r* = 5 × 1019 eV is deflected by different GMF configurations for all possible observed arrival directions on the northern sky, θobs (the deflection map; see Fig. 1 for the definition of the relevant quantities). Panel a: true deflection map. Panel b: mean deflection over the posterior samples for scenario A. Panel c: as in 5b, but the local measurements of the GMF now contain information on the LOS component as well as the POS component (scenario B). The additional information in this case causes a greater resemblance of the posterior mean to the true field, and so the deflection map is closer to Fig. 5a. Panel d: as in Fig. 5b, but the posterior is additionally constrained by the integrated data shown in Fig. 6b (scenario C). The scale of the color bar is kept up to 30 degrees to aid visual comparison. The red line on the color bar indicates the maximum deflection for each case. The dominant central feature of Fig. 5a is recovered in Figs. 5b–5d, since it is caused by the largest-scale features of the magnetic field, which we are able to infer in every case.

3.3 Scenario C: Local measurements with information on the plane of the sky supplemented by integrated line-of-sight measurements for the whole sky

In this section, we consider the inclusion of integrated constraints on the LOS component of the GMF as shown in Fig. 6b, while the local measurements at the dust clouds, simulating those obtained through polarized starlight tomography, are still projected on the celestial sphere as in Fig. 2b. Therefore, the likelihood we used now had the full form of Eq. (8).

In Figs. 3c and 4c, we show the mean and standard deviation of the angular distance error of the inferred UHECR arrival direction using the samples that were produced through the updated posterior for both local POS data and integrated LOS data. In comparison to scenario A, shown in Figs. 3a and 4a, the improvement in the ability to reconstruct the UHECR arrival direction is substantial in that the maximum mean angular error is reduced by a factor of about 1.5, the part of the POS where the maximum mean angular error occurs is greatly reduced, and the variance of the posterior is diminished by a factor of about 1.2. For the setting we considered, we thus showed that including integrated LOS data of the GMF (which is a much more realistic expectation than full 3D local measurements of scenario B) also leads to significantly better results with regard to recovering the arrival directions of UHECRs with a rigidity r_*.

3.4 Correlation length

In addition to the arrival direction of UHECRs, we can also infer the correlation length of the field for each of the above scenarios. When we denote the fluctuating component of the GMF as δB(x), the two-point correlation function is $$⟨\delta B_i(\mathbf{x})\delta B_i^{\ast }\left({\mathbf{x}}^{\mathrm{\prime }}\right)⟩=\xi (\mathbf{x},{\mathbf{x}}^{\mathrm{\prime }}),$$(17)

where a sum over i indices is understood. When we assume homogeneity and isotropy, the function ξ only is a function of the norm of the difference vector r ≡ x − x′, that is, ξ = ξ(r), where r = |r|. In this case, the correlation length is defined as $$L_c\equiv \frac{1}{\xi (0)}{\int }_0^{\mathrm{\infty }}dr\xi (r).$$(18)

As discussed in Section 2.2, we directly sampled GMF configurations from the posterior, but also inferred a parameterized power spectrum P(k), providing us with two different ways to calculate L_c: directly from the GMF samples, and from the inferred P(k). We recall that the power spectrum was modeled as depending on the magnitude of the wavevectors alone, implying an assumption of homogeneity and isotropy. This is not because we think that the GMF obeys these conditions, but rather due to our a priori ignorance regarding their violation. In truth as well as in our simulated ground truth, the magnetic field is vertically stratified, and homogeneity therefore cannot be assumed to hold along the z direction, but only on constant z slices. The individual GMF samples are informed by the data about violations of homogeneity and isotropy, as shown in Figs. 2 and 7. Therefore, caution should be exercised when calculating the correlation length directly from the samples.

Following Bendre & Subramanian (2022), we only computed the correlation length in slices of constant z, thus ending up with 2D two-point correlation function of the field, ξ(x, y) = ξ(ρ), for fixed z, where $\rho =\sqrt{(x^2+y^2)}$. We therefore define $$L_c(z)\equiv \frac{1}{\xi (0)}{\int }_0^{\mathrm{\infty }}d\rho \xi (\rho ),$$(19)

and $$L_c\equiv {⟨L_c(z)⟩}_z,$$(20)

where ⟨...⟩_z denotes averaging over the z-axis. For the ground truth, the z-averaged correlation length is calculated to be L_c = 59 pc. In Fig. 8, we calculate L_c for all the posterior samples of the three scenarios above and compare it with the L_c calculated for the ground truth (dashed vertical line). The correlation length is systematically overestimated. This effect is expected because the reconstruction of noisily and sparsely sampled random fields necessarily loses information at small length scales, which introduces a low-pass filtering effect over the true signal and therefore leads to an increased correlation length. For this reason, the Bayesian inference of the correlation length should be seen as imposing upper bounds. In agreement with the results of the previous sections, knowledge of local LOS information on the location of individual measurements (scenario B) significantly improves our capability to infer the correlation length in this case as well.

The presence of differential shear means that even the 2D correlation function of the ground truth deviates from isotropy, and it therefore admits two distinct correlation lengths. For simplicity, however, we neglected the effect of shear and calculated the correlation length as defined above, with the understanding that the deviation from isotropy is a higher-order effect. The L_c we calculated for the ground truth lies in the range of correlation lengths computed in Bendre & Subramanian (2022), who accounted for differential shear.

Finally, for comparison, we show in Fig. 9 the results for L_c calculated directly from Eq. (18), using the inferred power spectrum P(k) of each sample. In this case, while we still obtain a correlation length of ~100pc, the variance is much greater than the results shown in Fig. 8 for all three scenarios. The tendency to overestimate L_c arises because the stratification and the subsequent large-scale correlation along the z-axis is not taken into account separately in this method, as it was in the previous method, and this manifests itself in the power spectrum by the increased measured correlation length.

Fig. 6 The result of the inference of the GMF as seen on the northern POS. Panel a: averaged LOS component of the test magnetic field, shown in Fig. 2a. Panel b: noisy integrated data that were used along with the sparse and local data shown in Fig. 2b in order to define the LOS-informed posterior distribution. The noise covariance was set to 0.5 μG2, and the density of integrated measurements was 0.1 deg−2. Panel c: averaged LOS component of the mean 3D configuration of the approximating posterior distribution for the data of Figs. 2b and 6b.

4 Discussion

4.1 Identification of a systematic bias

In Fig. 2a, we showed that the ordered component of the field primarily lies (anti)parallel to the ±ŷ direction, which corresponds to a longitude l = ±90°. In Fig. 5a, this is reflected by the fact that the observed arrival directions parallel to the ordered component, (l, b) ≃ (±90°, 0°), are minimally deflected, while the maximum deflection occurs at the arrival directions perpendicular to the ordered component of the field. We call the map of Fig. 5a the deflection map of the GMF for a UHECR with a rigidity r_*. If the deflection map of the GMF for a given of rigidity were available, we would be able to identify the regions of the celestial sphere in which observed UHECRs with this rigidity are deflected most strongly.

A comparison of Fig. 5a with Fig. 3 yields a direct correlation between the regions of the deflection map and the mean angular error of our inferred arrival directions as a function of observed arrival direction for the same rigidity. In qualitative terms, this correlation suggests that for observed arrival directions perpendicular to the GMF zero mode, where the particles must have been deflected most strongly, our inference of their true arrival direction is more prone to a systematic bias. This bias is to be understood as the angular distance of the mean of our posterior distribution with respect to the true value.

Even though we might not be able to correct for this bias using our available data, knowledge of how severely the GMF alters the UHECR trajectories can help us to characterize the regions of the POS where our reconstructions are expected to be affected by it. While the corresponding deflection of the true GMF for a value of the UHECR rigidity will not be known a priori⁵, its structure is largely dictated by the dominating mean value of the field, which is generally well captured by our algorithm, as shown in Paper I. As shown in Figs. 5b-5d, we are able to accurately recover the large-scale features of the deflection map for all three considered scenarios, and in this way, we chart the parts of the POS in which the GMF will affect the UHECR trajectories most strongly, and by extension, the regions in which our arrival direction posterior might be shifted with respect to the true value.

Fig. 7 Individual samples from two of the posterior distributions considered. Panel a: the posterior distribution conditional to the data of Figs. 2b and 6b. Panel b: the posterior distribution is conditional to the data of Fig. 2b alone.

Fig. 8 Inference of the correlation length as defined via Eqs. (19) and (20) computed directly from the posterior samples for the three scenarios. The vertical line corresponds to the correlation length computed from the ground truth.

4.2 Caveats

While tomography using starlight polarization and Gaia data can provide the location of dust clouds in the local Galaxy as well as the POS orientation of the GMF at the location of each cloud, the POS direction of the GMF is generally not known, as this inference makes use of the properties of grain alignment, which cannot infer the POS directionality of the GMF (Tassis et al. 2018).

Furthermore, the integrated measurements used here assume that the integrated Galactic LOS component has been measured or inferred. In practice, the observables that need to be measured in order to estimate these integrals is the Faraday rotation measure and the dispersion measure. This means that even if the Galactic component is separated, it will still provide an average weighted over the thermal electron density. Therefore, we practically made the simplifying assumption that the thermal electron density is constant or known. In applications to the real GMF, the electron density will be treated as an additional degree of freedom to be inferred (Hutschenreuter et al. 2024). However, it must be noted that recent research suggests the possibility that local LOS data can be available, at least in part of the dataset (Tahani et al. 2022a; Tahani et al. 2022b).

In this analysis, we only studied the case of UHECRs with a fixed rigidity of r_* = 5 × 10¹⁹ eV. This is equivalent to assuming that the UHECR particles are protons of E = 5 × 10¹⁹ eV. In general, the composition of UHECRs is unknown and is most likely mixed, especially if some of the sources have a Galactic origin (Calvez et al. 2010; Kusenko 2011; Jiang et al. 2021). The closer examination of different composition scenarios will be the subject of future work. The robustness of this method must be further confirmed in different simulations (see, e.g., Gent et al. 2024; Korpi-Lagg et al. 2024).

Fig. 9 Inference of the correlation length as defined via Eq. (18), using the inferred power spectrum of each respective field. The vertical line corresponds to the correlation length of the ground truth.

4.3 Conclusions and outlook

We extended the analysis of Paper I to the case of more realistic LOS information and local data distribution. This was motivated by the fact that in real applications, the local GMF data obtained through stellar polarization tomography will not contain LOS information, and the distribution of these measurements will follow the distribution of dust clouds, which is not homogeneous, as was assumed in Paper I.

Additionally, the ground-truth GMF that was used in order to benchmark the performance of our inference algorithm was taken from an MHD simulation with the aim of studying the effect of our Gaussian approach to magnetic field configurations whose statistical properties more closely resemble those of the real GMF. Furthermore, we supplemented the existing framework in order to also include LOS-integrated information.

Our results show that while the complete absence of LOS information in the local data diminishes the accuracy of our inferred UHECR arrival directions, we are able to significantly correct for the effect of the GMF on the observed arrival directions even in this case, at least for the rigidity considered here. The inclusion of integrated LOS data for the GMF, which can realistically be expected to be part of our available information, is enough to provide accurate enough results, however.

Even in directions in which the angular distance between the inferred arrival direction and the true direction are strongest are we still able to correct for the effect of the GMF by a factor of 3 in the setting we considered. Additionally, by our ability to reconstruct the large-scale features of the field that dominate UHECR deflection, we are able to identify the regions of the POS in which our reconstructions are most likely to be affected by the maximum error. Finally, we are able to estimate the correlation length of the GMF with two different methods, one directly from the posterior samples, and the other using the inferred power spectra.

Acknowledgements

A.T. and V.P. acknowledge support from the Foundation of Research and Technology – Hellas Synergy Grants Program through project MagMASim, jointly implemented by the Institute of Astrophysics and the Institute of Applied and Computational Mathematics. A.T. acknowledges support by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “Third Call for H.F.R.I. Scholarships for PhD Candidates” (Project 5332). V.P. acknowledges support by the Hellenic Foundation for Research and Innovation (H.F.R.I.) under the “First Call for H.F.R.I. Research Projects to support Faculty members and Researchers and the procurement of high-cost research equipment grant” (Project 1552 CIRCE). The research leading to these results has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie RISE action, Grant Agreement no. 873089 (ASTROSTAT-II). This work also benefited greatly from discussions during the program “Towards a Comprehensive Model of the Galactic Magnetic Field” at Nordita in April 2023, which is partly supported by Nord-Forsk and the Royal Astronomical Society. A.T. would like to thank Vincent Pelgrims, Raphael Skalidis, Georgia V. Panopoulou, and Konstantinos Tassis for helpful tips and stimulating discussions. G.E. acknowledges the support of the German Academic Scholarship Foundation in the form of a PhD scholarship (“Promotionsstipendium der Studienstiftung des Deutschen Volkes”). P.F. acknowledges funding through the German Federal Ministry of Education and Research for the project ErUM-IFT: Informationsfeldtheorie fuer Experimente an Großforschungsanlagen (Foerderkennzeichen: 05D23EO1). The authors wish to thank an anonymous referee for insightful comments and suggestions.

Appendix A Simulated Magnetic Field

We briefly summarize the setup and results of the Galactic dynamo simulations that have been analyzed here. A detailed description of the numerical setup is presented in Bendre et al. (2015).

These are Magnetohydrodynamic (MHD) simulations of the Galactic interstellar medium (ISM). The simulation domain is an elongated box, located roughly at the solar neighbourhood of the Milky Way. It has dimensions of approximately 1 × 1 kpc in the radial (x) and azimuth (y) direction and ranges from approximately −2 to + 2 kpc in z direction, on either side of the Galactic mid-plane. It is split in a uniform Cartesian grid with a resolution of approximately 8.3 pc, and a set of non-ideal MHD equations is solved in this domain using the NIRVANA code (Ziegler 2004) (see Eq. 1 from Bendre et al. (2015) for the set of equations we have solved). Periodic boundary conditions were used in the y direction to incorporate the axisymmetry of the Galactic disc. The flat rotation curve is incorporated by allowing the angular velocity to scale inversely with the Galactic radius as Ω *prop; 1/R, with Ω₀ = 100 km s⁻¹ kpc⁻¹ at the centre of the box. Shearing periodic boundary conditions are used in the radial x direction to accommodate the aforementioned radial dependence of angular velocity. The initial density distribution of the ISM is in hydrostatic balance with the vertical gravity pointing towards the mid-plane, such that the vertical scale-height of the initial density was approximately 300 pc, with its value in the mid-plane of approximately 10⁻²⁴ g cm⁻³. A vertical profile of gravitational acceleration is adapted from Gilmore et al. (1989). The ISM in this box is stirred by supernovae (SN) explosions, which inject the thermal energy at random locations, at a rate of approximately 7.5 kpc⁻² Myr⁻¹. The vertical distribution of the explosions scale with the mass density. A piece-wise power law, similar to Sánchez-Salcedo et al. (2002), is used to model the temperature-dependent rate of radiative heat transfer, which along with SN explosions, roughly capture the observed multi-phase morphology of the ISM. We started the simulations with negligible initial magnetic fields of strength of the order of nG, and it grew exponentially to the strengths of the order of μG, with an e-folding time of about 200 Myr, such that the final energy density of the magnetic fields reached equipartition with the kinetic energy density of the ISM turbulence (shown in the right-hand panel of Fig. A.1). The exponential amplification of the magnetic energy saturated after about a Gyr, and coherent magnetic fields of scale-height close to 500pc were sustained in the box, consistent with the typical scale-height of GMFs (shown in the left-hand panel of Fig. A.1). The initial amplification and subsequent saturation phases of the magnetic field are termed here respectively, as kinematic and dynamical phases. This refers to the fact that magnetic fields from initial kinematic strengths amplify and reach the strengths that are dynamically significant to the turbulent flow. The growth and saturation of these large-scale fields are understood in terms of a self-consistent large-scale dynamo mechanism, governed by the SN-driven stratified helical turbulence and the Galactic differential rotation (Bendre et al. 2015). Finally, regarding the correlation length, the magnetic field is anisotropic, in that it is more correlated in the direction of shear. In the direction parallel to the shear the correlation length is calculated at approximately 100 − 150 pc, while in the direction perpendicular to it, it is in the range of 20 − 70 pc (see Fig. A3 of Bendre & Subramanian 2022). Neglecting shear for simplicity, we calculate the correlation length to be 63 pc.

Fig. A.1 Left: Time evolution of the vertical (z) profile of the azimuthal component of the magnetic field averaged over x − y plane. The color code is normalized by an exponential factor to compensate for an exponential growth of magnetic fields. The mean magnetic field eventually grows to a large-sale mode symmetric with respect to the Galactic mid-plane. Right: Time evolution of various contributions to magnetic energy, normalized to the turbulent kinetic energy (which stays roughly constant in time). The black solid line corresponds to the total magnetic energy contribution, the red dashed line corresponds to the magnetic energy of mean magnetic fields (averaged over the horizontal x − y planes) and with the blue dot-dashed line to the magnetic energy in the RMS magnetic fields. The magnetic energy is amplified exponentially for about a Gyr and eventually reaches an equipartition with turbulent kinetic energy. Figure adapted from Bendre & Subramanian 2022.

References

All Tables

All Figures

Fig. 1 Illustration and definition of quantities used in the analysis. A UHECR with a known rigidity r, enters the Galaxy with an arrival direction θtrue (red dot). Because of the GMF, it is deflected and is observed on Earth as arriving from θobs (black dot). The angular distance between θobs and θtrue is α, and it is the error that the GMF induces on the observed arrival direction. We traced back the particle through each GMF configuration sampled using NIFTy, thus ending up with a distribution of arrival directions P(θ|D), with D defined in Eq. (16). From the posterior samples drawn, we calculated the mean angular distances ⟨δθ⟩θ|D and ⟨δϕ⟩θ|D to the true and observed arrival directions, respectively, as well as the standard deviations for the former. The scales in this artificial example are larger for visual clarity, and do not correspond to an application of the method.

In the text

Fig. 2 Reconstruction of the simulated 3D magnetic field with the use of local data that lack LOS field component information. The blue sphere represents the celestial sphere. Panel a: ground truth; the GMF obtained as described in Appendix A. The field is rescaled so that it has an RMS norm of 5 μG. Panel b: synthetic data based on the ground truth of Fig. 2a. The radial component of the magnetic field is not measured. Panel c: the mean of the approximating posterior distribution attained via the geoVI algorithm based on the data provided in Fig. 2b. Panel d: mean of the approximating posterior distribution attained for the local data of Fig. 2b as well as integrated measurements of the radial component (Fig. 6b). The side of the box is about 1.6 kpc, or approximately 20 times the GMF correlation length.

In the text

Fig. 3 Mean angular error of the reconstruction (see Fig. 1) as a function of all possible arrival directions on the northern hemisphere for the case of a UHECR with a rigidity r* = 5 × 1019 eV. Panel a: magnetic field data consisting of local information with the LoS component are projected out (scenario A). Panel b: magnetic field data consisting of local information with the measured LOS component (scenario B) Panel c: as in the top left panel, but the data are supplemented by integrated LOS data (scenario C; see Fig. 6). The red and orange lines on the color bar indicate the maximum and mean values of the map, respectively.

In the text

	Fig. 4 As in Fig. 3, but for the corresponding angular error standard deviations as a function of the observed arrival direction.
In the text

Fig. 5 Amount by which a UHECR of rigidity r* = 5 × 1019 eV is deflected by different GMF configurations for all possible observed arrival directions on the northern sky, θobs (the deflection map; see Fig. 1 for the definition of the relevant quantities). Panel a: true deflection map. Panel b: mean deflection over the posterior samples for scenario A. Panel c: as in 5b, but the local measurements of the GMF now contain information on the LOS component as well as the POS component (scenario B). The additional information in this case causes a greater resemblance of the posterior mean to the true field, and so the deflection map is closer to Fig. 5a. Panel d: as in Fig. 5b, but the posterior is additionally constrained by the integrated data shown in Fig. 6b (scenario C). The scale of the color bar is kept up to 30 degrees to aid visual comparison. The red line on the color bar indicates the maximum deflection for each case. The dominant central feature of Fig. 5a is recovered in Figs. 5b–5d, since it is caused by the largest-scale features of the magnetic field, which we are able to infer in every case.

In the text

Fig. 6 The result of the inference of the GMF as seen on the northern POS. Panel a: averaged LOS component of the test magnetic field, shown in Fig. 2a. Panel b: noisy integrated data that were used along with the sparse and local data shown in Fig. 2b in order to define the LOS-informed posterior distribution. The noise covariance was set to 0.5 μG2, and the density of integrated measurements was 0.1 deg−2. Panel c: averaged LOS component of the mean 3D configuration of the approximating posterior distribution for the data of Figs. 2b and 6b.

In the text

	Fig. 7 Individual samples from two of the posterior distributions considered. Panel a: the posterior distribution conditional to the data of Figs. 2b and 6b. Panel b: the posterior distribution is conditional to the data of Fig. 2b alone.
In the text

	Fig. 8 Inference of the correlation length as defined via Eqs. (19) and (20) computed directly from the posterior samples for the three scenarios. The vertical line corresponds to the correlation length computed from the ground truth.
In the text

	Fig. 9 Inference of the correlation length as defined via Eq. (18), using the inferred power spectrum of each respective field. The vertical line corresponds to the correlation length of the ground truth.
In the text

Fig. A.1 Left: Time evolution of the vertical (z) profile of the azimuthal component of the magnetic field averaged over x − y plane. The color code is normalized by an exponential factor to compensate for an exponential growth of magnetic fields. The mean magnetic field eventually grows to a large-sale mode symmetric with respect to the Galactic mid-plane. Right: Time evolution of various contributions to magnetic energy, normalized to the turbulent kinetic energy (which stays roughly constant in time). The black solid line corresponds to the total magnetic energy contribution, the red dashed line corresponds to the magnetic energy of mean magnetic fields (averaged over the horizontal x − y planes) and with the blue dot-dashed line to the magnetic energy in the RMS magnetic fields. The magnetic energy is amplified exponentially for about a Gyr and eventually reaches an equipartition with turbulent kinetic energy. Figure adapted from Bendre & Subramanian 2022.

In the text