A self-consistent framework to study magnetic fields with strong gravitational lensing and polarized radio sources

ABSTRACT

We introduce a unified approach that, given a strong gravitationally lensed polarized source, self-consistently infers its complex surface brightness distribution and the lens galaxy mass–density profile, magnetic field and electron density from interferometric data. The method is fully Bayesian, pixellated, and three-dimensional: the source light is reconstructed in each frequency channel on a Delaunay tessellation with a magnification-adaptive resolution. We tested this technique using simulated interferometric observations with a realistic model of the lens, for two different levels of source polarization and two different lensing configurations. For all data sets, the presence of a Faraday rotating screen in the lens is supported by the data with strong statistical significance. In the region probed by the lensed images, we can recover the rotation measure and the parallel component of the magnetic field with an average error between 0.6 and 11 rad m${}^{-2}$ and 0.3 and 3 nG, respectively. Given our choice of model, we find the electron density is the least well-constrained component due to a degeneracy with the magnetic field and disc inclination. The background source total intensity, polarization fraction, and polarization angle are inferred with an error between 4 and 10 per cent, 15 and 50 per cent, and 1–12 deg, respectively. Our analysis shows that both the lensing configuration and the intrinsic model degeneracies play a role in the quality of the constraints that can be obtained.

1 INTRODUCTION

Magnetic fields play an important role in the dynamics and star formation processes of galaxies and, therefore, in the evolution of these objects (e.g. Hennebelle & Inutsuka 2019; Martin-Alvarez et al. 2020). They are believed to grow through the dynamo effect (e.g. Parker 1992; Kulsrud 1999; Kulsrud & Zweibel 2008; Martin-Alvarez et al. 2018; Rodrigues et al. 2019; Seta & Federrath 2020) from magnetic-seed fields, whose origin is still speculative. To this day, little is understood about their origin as well as their amplification (e.g. Rees 2006).

To answer the many open questions regarding the formation and evolution of magnetic fields and how they affect various astrophysical processes, it is necessary to study their structural and statistical properties at early epochs. Owing to the sensitivity and angular resolution limitations of current observing facilities, direct measurements of the magnetic field structure in individual galaxies have so far mainly been limited to the local Universe (e.g. Beck & Wielebinski 2013; Beck 2015; Lopez-Rodriguez et al. 2023, and references therein).

Faraday rotation, the process by which the plane of polarization for linearly polarized radiation propagating through a magneto-ionic medium is rotated (Burn 1966; Sokoloff et al. 1998; Beck & Wielebinski 2013), is commonly used as an indirect tool for detecting and probing cosmic magnetic fields in distant galaxies (e.g. Kronberg & Perry 1982; Welter, Perry & Kronberg 1984; Oren & Wolfe 1995; Bernet et al. 2008; Kronberg et al. 2008; Bernet, Miniati & Lilly 2013; Farnes et al. 2017). However, this approach yields only a statistical measure of the magnetic field, which is hindered by the fact that the contributions from the Milky Way and the intergalactic medium (IGM) to the observed Faraday rotation remain unknown.

Strong gravitational lensing provides an independent channel to study magnetic fields in the individual distant lens galaxies (Greenfield, Roberts & Burke 1985; Narasimha & Chitre 2008; Mao et al. 2017) and in the background sources (Geach et al. 2023; Chen et al. 2024; de Roo et al. 2024). The basic idea is that strong gravitational lensing conserves the surface brightness, is achromatic and non-polarizing (Dyer & Shaver 1992), and, therefore, conserves the polarization properties of the background object. This allows one to study magnetic fields in distant lensed sources with a signal-to-noise ratio (SNR) and angular resolution that would not be otherwise possible at high redshifts (Geach et al. 2023; Chen et al. 2024; Roo et al. 2024). On the other hand, a non-homogeneous magnetized medium within the lens galaxy will induce a differential change in the observed polarization properties of the lensed images. This propagation effect can be used to constrain the magnetic field and electron density in the lens galaxy, free from contamination by the Milky Way (at least for galaxy-scale lenses) and is independent of the intrinsic properties of the source (Narasimha & Chitre 2008; Mao et al. 2017).

Due to the limited number of currently known strongly lensed polarized sources (less than 10; Greenfield et al. 1985; Biggs et al. 1999; Patnaik et al. 2001; Narasimha & Chitre 2008; Mao et al. 2017; Biggs 2023), this is a field that is still in its infancy and that stands to significantly gain from the large number (⁠$\sim {10}^5$⁠) of strong gravitational lens systems that are expected to be discovered with the Square Kilometre Array (SKA; McKean et al. 2015) or the next generation Very Large Array (ngVLA). For example, assuming that the MeerKAT International GigaHertz Tiered Extragalactic (MIGHTEE) source population is representative of what will be observed with the SKA and extrapolating the fraction of polarized sources recently reported in the polarization early-science data release (Taylor et al. 2024), one should expect $\sim 2\times {10}^3$ lensed polarized sources detectable with the SKA at a surface brightness limit of 20 $\mu$Jy beam${}^{-1}$⁠.

Taking full advantage of future data will require lens modelling techniques that can infer the surface brightness distribution of the polarized source, as well as the mass distribution and magneto-ionic properties of the lens for a diverse population of lenses and sources. Burns (2002) has proposed a method based on the misalignment between the polarization vector and the source morphology (or any 2-d vector on the source plane) that is induced by the presence of a gravitational lens (see also Kronberg et al. 1991). This approach only applies to background sources, such as radio jets and shocks, where these two quantities are expected to be intrinsically aligned on the source plane. Focusing on observations of lensed quasars, Greenfield et al. (1985) and Mao et al. (2017) have decoupled the lens and the Faraday rotation analyses from each other and treated the lensed images of the same background unresolved object as they were, to a certain extent, independent. This procedure is unsuitable for resolved sources, as the multiple lensed images of the same source component cannot be trivially identified without a lens model. Moreover, it assumes that the magnification does not change across the extent of the source and between the different lensed images, and it relies on cleaned maps for which the noise is correlated. In their analysis of a strongly lensed star-forming galaxy observed with the Atacama Large (sub)Millimetre Array (ALMA), Geach et al. (2023) have first derived a model for the lens mass distribution using ancillary optical and millimetre data. This model was then used to lens forward a Gaussian model for the Stokes $\mathit{Q}$ and $\mathit{U}$ source surface brightness, assuming a constant polarization angle. Due to the relatively high frequency of their observations (242 GHz), they could neglect the effect of Faraday rotation in the lens galaxy or along the line of sight to the background source.

In this paper, we introduce the first Bayesian gravitational lens modelling technique that, from the same data, self-consistently infers the lens galaxy mass distribution, magnetic field, and electron density, and the polarized complex surface brightness distribution of the background source. The latter is modelled using a pixellated surface brightness distribution, requiring no strong assumptions on the morphological and polarization properties of the source. Our technique is, therefore, suitable for modelling a large variety of lensed objects (e.g. galaxies and radio jets). The data are fitted in the visibility space where the noise is well approximated by an uncorrelated Gaussian distribution. The method is three-dimensional (one frequency and two spatial dimensions) and fully forward; therefore, it automatically takes into account beam and bandwidth depolarization.

The paper is organized as follows. In Section 2, we provide a detailed description of our formalism, which we then test using simulated data that is generated according to Section 3 and analysed using the modelling strategy outlined in Section 4. We discuss the performance of the technique in Section 5 and its current limitations in Section 6. We present our findings and a future outlook for the method in Section 7.

2 DESCRIPTION OF THE METHOD

This section introduces our approach to analysing strong gravitational lensing observations with a polarized source and in the presence of a magnetized plasma in the lens galaxy.

We assume that the foreground lensing galaxy is not a source of emission at the observed frequencies (for example, 90 per cent of the foreground lenses found as part of the Cosmic Lens All-Sky Survey are radio-quiet; Browne et al. 2003) and that the contribution to the Faraday rotation from the line of sight can be safely neglected. The latter assumption is justified by the fact that magnetic fields in the cosmic web are expected to be significantly smaller than those in galaxies and by the fact that the difference in the light-ray path lengths is much smaller than the path length itself. For example, numerical simulations as well as observations constrain the rotation measure (RM) arising from the IGM of the order of $\sim {10}^{-1}$ to $\sim 10$ rad m${}^{-2}$⁠, hence several orders of magnitude smaller than what is observed in galaxies (e.g. Vernstrom et al. 2019; O’Sullivan et al. 2020; Carretti et al. 2022).

As the angular separation of the multiple lensed images is small (a few arcsec on the plane of the observer), we also assume the effect of the Milky Way to be the same for all of the lensed images and, therefore, not significant. For simplicity, we limit our focus to the contribution of a regular large-scale magnetic field and defer the inclusion of a random component to a follow-up paper. We also ignore the effect of a possible field of axion-like particles on the observed polarization angle (Basu et al. 2021). Under these conditions, we can treat any relative difference in the lensed image polarization properties as arising from a single external Faraday rotating screen at the plane of the lens. At this stage, we assume the source will be constant throughout the observations and ignore the physical processes internal to this object. As discussed in the following, our formalism can be extended to include physical assumptions on the nature of the source in a self-consistent fashion.

Our method builds upon pronto, the Bayesian grid-based software initially developed for optical data by Vegetti & Koopmans (2009). It was then further adapted for large interferometric data by Powell et al. (2021) and extended to the three-dimensional domain by Rizzo et al. (2018). Below, we describe in detail all new aspects of the method and refer the reader to the above papers as well as Rybak et al. (2015a, b) and Ritondale et al. (2019) for further details on the original implementation and its subsequent developments.

2.1 The response operator

Given the electric field $\mathit{e}={(e_{\mathrm{R}}(\mathit{y}),e_{\mathrm{L}}(\mathit{y}))}^T$ of the background polarized source, the corresponding complex surface brightness distribution as a function of the position on the source plane $\mathit{y}$ is given by the following coherency vector (Hamaker, Bregman & Sault 1996; Smirnov 2011a)

Here, we have assumed that the instrument has circularly polarized feeds, which measure the right-handed (R) and left-handed (L) components of the incoming wave. Rather than working in the space of four complex components, we hereafter represent the source surface brightness distribution using the four real-valued Stokes parameters $\mathit{I}$⁠, $\mathit{Q}$⁠, $\mathit{U}$⁠, and $\mathit{V}$⁠:

where ${\mathsf{F}}_{\mathrm{R}\mathrm{L}}$ is the ‘feed operator’ that forms the Stokes parameters from the appropriate linear combinations of raw RR, RL, LR, and LL correlations (Hamaker et al. 1996). Our method thus extends trivially to the case of orthogonal linearly polarized (X and Y) feeds.

The gravitational effect of a strong gravitational lens galaxy between the source and the observer can be expressed in terms of the so-called lensing operator $\mathsf{L}$ (see Vegetti & Koopmans 2009), which is a function of the lens projected mass density distribution (see Section 3 for a definition). It allows one to relate, via the lens equation, a position on the lens plane $\mathit{x}$ to its corresponding position on the source plane $\mathit{y}$ while taking into account the conservation of surface brightness. As described by Vegetti & Koopmans (2009), the latter is defined on a Delaunay tessellation, which is automatically adaptive with the lensing magnification. The lens plane is defined on a regular Cartesian grid of arbitrary resolution, though one should be careful to ensure Nyquist sampling. In the absence of a magnetized medium in the lens galaxy, the source light ${\mathit{s}}_j$ and the observed visibility data ${\mathit{d}}_j$ in each frequency channel j are related to each other as follows:

where the NUFFT operator ${\mathsf{D}}_j$ includes, from left to right, a gridding operation, a Fast Fourier Transform, an apodisation correction, and a zero-padding/masking operation (see Powell et al. 2021). Here, ${\mathit{n}}_j$ is a vector encoding the observational noise, which is assumed to be Gaussian and uncorrelated between visibility points (Wucknitz 2002). We assume the data ${\mathit{d}}_j$ to be already calibrated and refer the reader to Smirnov (2011a, b) for a description of how to extend the response operator further to include calibration processes. From equations (2) and (3), it can be seen that the conservation of the polarization properties (angle and fraction) of the background source by an intervening lens directly follows from the fact that lensing conserves surface brightness. Throughout this paper, we work in the basis of Stokes brightnesses, where the data ${\mathit{d}}_j$ are assumed to have been transformed into Stokes visibilities via, for example, ${\mathit{d}}_j={\mathsf{F}}_{\mathrm{R}\mathrm{L}}{\mathit{d}}_{\mathrm{RL},j}$⁠.

The main novelty of this paper is that we now include the effect of the lens galaxy magnetized interstellar medium (ISM), which acts as a single external Faraday screen and induces a differential rotation of the polarization angle between the multiple lensed images of the same source (Mao et al. 2017). Taking into account this effect leads to a modified version of the forward model described by equation (3) as follows:

where we have introduced an external Faraday screen operator ${\mathsf{S}}_j$⁠. In the Stokes basis, ${\mathsf{S}}_j$ takes the form

For a simple screen, the amount of rotation of the plane of polarization is directly related to the Faraday depth in $\mathrm{r}\mathrm{a}\mathrm{d}{\mathrm{m}}^{-2}$ (Brentjens & de Bruyn 2005; Amaral, Vernstrom & Gaensler 2021), so that

Equation (7) makes the wavelength dependence (⁠${\lambda }_j$⁠) of the polarization rotation explicit; it is by modelling multiple frequency channels simultaneously that we can recover the properties of the Faraday screen. The RM in the observer’s frame in rad m${}^{-2}$ is

Here, $n_{\mathrm{e}}(l,\mathit{x})$ is the electron density, $B_{\mathrm{L}\mathrm{O}\mathrm{S}}(l,\mathit{x})$ is the line-of-sight component of the lens magnetic field, $\mathrm{d}l$ is the path length, $z_{\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{s}}$ is the redshift of the lensing galaxy, c is the speed of light, and $m_{\mathrm{e}}$ and e are the rest mass and charge of the electron, respectively.

2.2 Bayesian inference

In this section, we provide a summary of the Bayesian inference approach that allows us to simultaneously recover the source polarized surface brightness distribution $\mathit{s}$⁠, the lens mass distribution (see Section 3.2.1 for a definition), the magnetic field and the electron density parameters (see Section 3.2.2 for a definition). For simplicity, we refer to the lens parameters (describing its mass distribution, magnetic field, and electron density) collectively as $\mathit{\eta }$⁠. Our inference procedure is the same as the one introduced by Vegetti & Koopmans (2009) and further developed by Rizzo et al. (2018) and Powell et al. (2021). For completeness, we report the main key points here. In the following, we refer to $\mathit{s}$ and $\mathit{d}$ as the concatenation of the source and data vectors in the different frequency channels, that is, $\mathit{s}=\{{\mathit{s}}_0,..,{\mathit{s}}_{{\mathrm{n}}_{\mathrm{c}\mathrm{h}}}\}$ and $\mathit{d}=\{{\mathit{d}}_0,..,{\mathit{d}}_{{\mathrm{n}}_{\mathrm{c}\mathrm{h}}}\}$⁠.

2.2.1 Source reconstruction

For a given choice of the lens parameters $\mathit{\eta }$⁠, we can infer the maximum posterior (MAP) source ${\mathit{s}}_{\mathrm{M}\mathrm{P}}$ by maximizing the following probability:

Here, $P(\mathit{d}|\mathit{s},\mathit{\eta })$ is the likelihood and $P(\mathit{s}|{\mathit{\lambda }}_{\mathrm{s}},\mathsf{R})$ is a regularizing prior of strength ${\mathit{\lambda }}_{\mathrm{s}}$ (one per correlation) and order $\mathsf{R}$ on the source. Assuming the noise has a Gaussian distribution, we can express the likelihood function as:

where ${\mathsf{C}}_j$ is the diagonal covariance matrix for the frequency channel j, and

The response operator ${\mathsf{M}}_j$ is given by equation (5). Vegetti & Koopmans (2009) and Powell et al. (2021) have assumed a quadratic form for the source prior $P(\mathit{s}|{\mathit{\lambda }}_{\mathrm{s}},\mathsf{R})$ with a zero mean, that is,

with

and

Under this assumption, maximizing the posterior probability in equation (9) to obtain ${\mathit{s}}_{\mathrm{M}\mathrm{P}}$ reduces to solving the following system of linear equations for ${\mathit{s}}_j$

Rizzo et al. (2018) have further developed this approach to include an analytical and physically motivated model ${\mathit{s}}_{\mathrm{h}\mathrm{y}\mathrm{p}}$ as a hyper-prior to the pixellated source. In this case, equations (13) and (16) are respectively modified as follows:

and

The hyper-parameters ${\mathit{\eta }}_{\mathrm{h}\mathrm{y}\mathrm{p}}$ defining ${\mathit{s}}_{\mathrm{h}\mathrm{y}\mathrm{p}}$ also become free (hyper) parameters of the model. Including a hyper-prior allows us to retain the freedom of a pixellated source while introducing a physical prior in a forward and self-consistent way. For example, ${\mathit{s}}_{\mathrm{h}\mathrm{y}\mathrm{p}}$ could be related to the intrinsic properties of the source and its magnetic field, as well as include the effect of internal physical processes, such as time variations in the source polarization properties. In this paper, we follow Powell et al. (2021) and assume ${\mathit{s}}_{\mathrm{h}\mathrm{y}\mathrm{p}}=\mathbf{0}$⁠. For generality, we include the source hyper-prior term in the relevant equations in the following section.

2.2.2 Non-linear parameters and model comparison

We obtain the most probable parameters for the lens mass, magnetic field, and electron density, and the source regularization level by maximizing the following posterior probability,

If we assume a uniform prior for $\mathit{\eta }$⁠, $\log {\mathit{\lambda }}_{\mathrm{s}}$⁠, and ${\mathit{\eta }}_{\mathrm{h}\mathrm{y}\mathrm{p}}$⁠, maximizing equation (19) is equivalent to maximizing the following expression (see Rizzo et al. 2018, for a derivation),

with $E=E_{\mathrm{D}}+{\mathit{\lambda }}_{\mathrm{s}}{E}_{\mathrm{R}}$⁠. For each step in this non-linear optimization, we calculate the corresponding maximum a posterior source ${\mathit{s}}_{\mathrm{M}\mathrm{P}}$ by solving the linear system (16) or (18), depending on the choice of hyper-prior.

As the magneto-ionic properties of distant galaxies are still poorly understood, it is important to consider and compare different models. Within our Bayesian inference approach, we can achieve this goal by comparing different assumptions for the magnetic field and electron density in terms of their Bayesian evidence, that is, the denominator in equation (19), which we calculate using MultiNest by Feroz et al. (2019). This also allows us to explore the multidimensional posterior distribution and identify potential degeneracy. Moreover, by comparing the Bayes factor of a model that includes the Faraday screen with one that does not, we can quantify the statistical significance with which the screen is detected in the first place.

3 SIMULATED DATA

This section describes the process for creating eight simulated interferometric observations. Fig. 1 shows the steps in the data simulation pipeline, going from the input surface brightness distribution of the polarized source to the (not deconvolved) images produced by an interferometer for two different data sets. In Section 3.1, we introduce our models for the background source, which we then lens forward via the mass–density model described in Section 3.2.1. The lensed images are then Faraday rotated according to the magnetic field and electron density models provided in Section 3.2.2. Table 1 lists our choice of input lens and source parameters. Details on the simulated observations are given in Section 3.3.

3.1 Source model

The input source model has a Gaussian surface brightness profile with a total flux-density $S_{\nu }$ that varies with frequency $\nu$ via the power law $S_{\nu }\propto {\nu }^{\alpha }$⁠, where the spectral index is set to $\alpha =-0.5$⁠. This choice of spectral index is consistent with a flat-spectrum radio source but is not expected to strongly influence our results since it mainly ensures a close to uniform SNR as a function of frequency. Our goal is to constrain the RM in different scenarios. Therefore, we consider two different source polarization fractions (⁠$P_{\mathrm{f}}$⁠) and two different source polarization angles (⁠$P_{\theta }$⁠), which are defined as,

and

respectively. For the former, we consider a low (3 per cent) and a high (12 per cent) polarization fraction case, and for the latter, we use a polarization angle of 87 and 78 deg. We also consider two different lensing configurations derived from the same input mass model. These are obtained from two different source sizes and orientations of the major axis relative to the lens mass density distribution. For the rest of this paper, we refer to the source that has the smaller major-axis (10 mas) with a PA that is misaligned relative to the lens (PA = 70 deg) as $s_1$⁠, while $s_2$ has a larger major-axis (15 mas) and is aligned with the mass distribution (PA = –30 deg) to maximize the radial structure in the lensed images (see Fig. 1). Our choice of source properties and lens parameters leads to a larger fraction of $s_1$ than $s_2$ being quadruply imaged. In all cases, we assume a linearly polarized source and set Stokes $\mathit{V}$ to zero. For each Stokes parameter, we generate a total of 13 frequency channels between 2 and 4 GHz. We note that the source properties, besides its polarization fraction or, equivalently, its polarized SNR do not affect our constraints on the lens galaxy Faraday screen. This is due to the fact that the measurement is differential in nature. In this context, the source simply acts as a backlight from which differences in the observed lensed surface brightness can be used to determine the contribution of the lens magneto-ionic medium.

3.2 Lens model

3.2.1 Mass model

We consider a two-component model for the mass density distribution of the lens galaxy, made of a single elliptical power law plus an edge-on disc for all of the simulated data. The lensing convergence for the power-law term is defined as,

where ${\kappa }_0$ is the amplitude, $\gamma$ is the 3D slope, q is the projected axial ratio and $r_{\mathrm{c}}$ is the core radius. Together with the PA $\theta$ and the position coordinates x and y, they represent the geometrical parameters of the lens mass model. Following Keeton (2011), we approximate the convergence of an elliptical edge-on disc component as the sum of Kuzmin discs, such that,

with the coefficients $s_i$ and ${\kappa }_i$ re-fitted using those parameters of the lens system B0712$+$472, as determined by Hsueh et al. (2017), which is known to have an edge-on disc that crosses the lensed images. In the expression above, ${\kappa }_{0,\mathrm{d}}$ is the central convergence, $R_{\mathrm{d}}$ is the scale length, and $q_{\mathrm{d}}$ is the axis-ratio. We also include an external shear component defined by a strength $\mathrm{\Gamma }$ and a PA ${\mathrm{\Gamma }}_{\theta }$⁠. Our choice of input parameters for the lensing mass distribution is based on modelling the real gravitational lens system B0712$+$472 (Hsueh et al.2017).

For consistency, we use the same lens mass model for all of the data sets investigated here since Powell et al. (2021) have already demonstrated that our lens modelling technique can robustly recover the lens model parameters. Due to the different source surface brightness distributions, this same mass model will produce different configurations of the lensed emission (see Fig. 1).

3.2.2 Magnetic field and electron density model

We adopt a bisymmetric azimuthal magnetic-field model with a vertical extension for all simulated data sets to account for the disc component (see Fig. 2). Our expression is a slight modification from that presented by Stepanov et al. (2008), and in cylindrical coordinates, is defined as,

where p and $\chi$ are the pitch and tilt angles, respectively. $B_{\mathrm{m}}$ and $\chi$ are respectively given by,

and

The latter determines the magnetic field vertical structure and strength. Following Stepanov et al. (2008), we consider a Gaussian profile for the electron density, such that,

The disc component in both $\mathit{B}$ and ${\mathit{n}}_{\mathbf{e}}$ have a PA that matches the disc component of the lens mass density distribution. We set the inclination angle to 90 deg for half of the data sets and 45 deg for the other half. This leads to a rotating screen that is stronger (peak RM${}_{\mathrm{m}\mathrm{a}\mathrm{x}}\sim$ 300 rad m${}^{-2}$⁠) but is also more localized (i.e. it affects a smaller portion of the lensed images, and it fluctuates on smaller spatial scales) in the first case, and weaker (peak RM${}_{\mathrm{m}\mathrm{a}\mathrm{x}}\sim$ 120 rad m${}^{-2}$⁠), but more extended and slowly changing for the second one. For the disc component with a 45 deg inclination angle, we do not re-align the mass model with $\mathit{B}$ and ${\mathit{n}}_{\mathbf{e}}$ because this would change the lensing configuration, preventing a direct comparison in terms of the Faraday screen. In all cases, we choose a value of $B_0$⁠, consistent with the mean coherent magnetic field strength inferred by Mao et al. (2017) for the gravitational lens system B1152+199. All other parameter values are taken from Stepanov et al. (2008). The adopted values for the magnetic field and electron density models are listed in Table 1.

It is important to note that whether these models are a good representation of distant lens galaxies is currently unknown. Our choice for the magnetic field model was guided by the only existing published work on the magnetic field of a high-redshift lens galaxy by Mao et al. (2017). They showed that a bisymmetric profile is a possible, if not the preferred, model for the main deflector in the gravitational system B1152+199. In the future, one can explicitly test these assumptions from observations of the many strong gravitational lens systems expected to be discovered. This paper aims to test our modelling approach in a reasonably realistic scenario for which the above models are sufficient.

3.3 Observations

The final part of our simulation pipeline adds the response of the interferometer. This is done to simulate an observed data set that is representative of the angular resolution and sensitivity of next-generation instruments, which are currently being constructed or are in development. However, the actual array configuration and receiver deployment have still to be finalised for the SKA-MID and the ngVLA. From the model surface brightness distributions presented in Fig. 1, we see that the inclusion of the disc component results in large-scale changes to the Stokes $\mathit{Q}$ and $\mathit{U}$ emission that are on the order of 100 to 200 mas in size; we use this to inform the minimum angular resolution needed to identify these variations in the observed data, which can be achieved through a combination of the baseline length and the choice of observing frequency.

We generated visibility data with an effective baseline length of about 150 km; this is consistent with the maximum baseline length of the Array Assembly 4 of the SKA, which is expected to be completed by 2030 and should be an order of magnitude more sensitive than the VLA. We note that observing at a higher frequency can compensate for a shorter maximum baseline length. Our frequency choice is based on having a sufficient fractional bandwidth (⁠$>50$ per cent) from which the Faraday rotation of the lensed emission can be reliably detected while providing a high enough angular resolution (150 mas) so that beam de-polarization does not overly affect our analysis. For this, we choose observing frequencies between 2 and 4 GHz. However, once the final array configuration and receiver deployment of the SKA-MID are fixed, we will perform more realistic simulations to determine the appropriate observational set-up in terms of observing frequency, bandwidth and depth. Here, we aim to validate our method while providing some indication of the accuracy and precision with which the magneto-ionic medium within intermediate redshift galaxies can be recovered.

We create the data sets with the Common Astronomy Software Applications (casa; CASA Team et al. 2022), using a template that is based on observations taken with the VLA, except that the maximum baseline length was increased from 36 to 150 km, by re-scaling all baselines by the same value. This interferometer has 27 antennas that provide a good snapshot uv-coverage and produces 351 baselines per frequency and time sampling. As discussed above, we used 13 discrete narrow channels spread over the total bandwidth between 2 and 4 GHz to keep the wavelength coverage and to avoid de-polarization. In principle, using all of the channels (2 MHz in width) across the full bandwidth is possible but computationally more expensive. The visibility integration time is 2 s, and the total integration time is 6 min. Overall, this produced a data set with around 820 000 visibilities per correlations (RR, LL, RL, and LR), which is sufficiently small for testing our method.

4 MODELLING STRATEGY

We model each data set using the pixellated and regularized method described in Section 2. The free parameters of the model are those describing the lens mass density distribution, magnetic field, and electron density, and the source surface brightness profile. At every stage of the modelling procedure, we keep the core radius of the mass density fixed, and we force the magnetic field and electron density to be aligned and centred with each other. As the normalizations $B_0$ and $n_0$ [see equations (26) and (28)] are degenerate with each other, we refactor them into a single free parameter $A_0=B_0\times n_0$⁠.

One also has to infer the regularization level ${\mathit{\lambda }}_{\mathrm{s}}=\{{\lambda }_{\mathrm{s},\mathit{I}},{\lambda }_{\mathrm{s},\mathit{Q}},{\lambda }_{\mathrm{s},\mathit{U}}{\lambda }_{\mathrm{s},\mathit{V}}\}$ for a given form of regularization $\mathsf{R}$⁠. Given the significant difference in SNR between the Stokes parameters, we found that leaving the regularization constants to change independently from one another is necessary to obtain a satisfactory source reconstruction for each polarization component.

From equation (6), it can be seen that the Stokes $\mathit{I}$ surface brightness, which provides the constraints to the lensing potential, is unaffected by the presence of an ionized medium in the lens. We note that interstellar scattering or free–free absorption from free electrons can affect the Stokes $\mathit{I}$ emission, but these are the focus of a companion paper (Kaung et al., in preparation) and are not considered further here. As a result, we expect the parameters describing the lens mass–density distribution and its magneto-ionic properties to be largely independent of each other. Hence, we model the data in four stages following the Bayesian inference approach described in Section 2.2. (i) First, we neglect the contribution of the lens as a Faraday rotating screen and infer the MAP lens mass density distribution, source and source regularization level ${\lambda }_{\mathrm{s},\mathit{I}}$ from the Stokes $\mathit{I}$⁠. (ii) We then keep the lens mass parameters fixed to the results obtained from step (i), include the effect of the Faraday-rotating screen and obtain the MAP of the lens magnetic field, electron density, the source and its regularization levels by modelling all Stokes parameters, as described in Section 2.2. (iii) We use MultiNest to quantify the statistical errors on all of the parameters assuming a uniform prior for $\mathit{\eta }$ and $\log {\mathit{\lambda }}_{\mathrm{s}}$⁠. (iv) By comparing the Bayesian evidence of those models with and without the contribution of the Faraday screen, we obtain a measure of the statistical significance of the magnetized plasma within the foreground lens galaxy.

5 RESULTS

In this section, we present our results for the simultaneous modelling of the mass distribution and magneto-ionic medium of the lens, as well as the properties of the recovered background source. We quantify the quality of our reconstructions of each model parameter x in terms of the root-mean-square residuals (i.e. $\mathrm{\Delta }x$⁠) and the root-mean-square residuals normalized to the peak input value (i.e. ${\sigma }_x$⁠). Both quantities are averaged over the region where the lensed images are at least 20 per cent of the peak surface brightness. Our results are summarized in Figs 3 to 9 and in Tables 2 and 3 (see also Figs A1 to A10 of the Appendix).

5.1 Lens mass–density distribution

From modelling the Stokes $\mathit{I}$ emission, we recover the lens mass density distribution parameters. Fig. 3 shows the inferred posterior distribution for data sets D${}_6$ and D${}_8$⁠, as an example. Table 2 lists the mean parameter values with the corresponding 68 per cent confidence-level uncertainties for the same data sets. An equal level of constraint is obtained for all data sets with the same lensing configuration, as expected from their matching data quality (angular resolution and SNR).

5.2 Rotation measure

As indicated by the Bayes factor $\mathrm{\Delta }\log \mathcal{E}=\log {\mathcal{E}}_{\mathrm{s}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}-\log {\mathcal{E}}_{\mathrm{n}\mathrm{o}\mathrm{s}\mathrm{c}\mathrm{r}\mathrm{e}\mathrm{e}\mathrm{n}}$⁠, the presence of a Faraday rotating screen in the lens galaxy is detected with a high statistical significance for all data sets (above 180$\sigma$⁠) and even more so for those generated from the source $s_2$ (above 260$\sigma$⁠; see Table 3). The latter result is expected because a larger source probes a larger fraction of the lens plane. At a fixed polarization fraction, we find that different lensing configurations are better at probing different RM distributions: the screen is more strongly detected when it has an inclination angle of 90 and 45 deg for $s_1$ and $s_2$⁠, respectively. This is due to the more localized effect of the 90-deg inclination model. The source $s_2$⁠, unlike $s_1$⁠, leads to a counter-image that only marginally overlaps with a positive peak in the RM distribution. On the other hand, the 45-deg inclination model is more extended and hence is more strongly detected with a larger source.

Fig. 4 displays our best (D${}_1$⁠; ${\sigma }_{\mathrm{R}\mathrm{M}}$ = 0.2 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$ = 0.6 rad m${}^{-2}$⁠) and worst (D${}_8$⁠; ${\sigma }_{\mathrm{R}\mathrm{M}}$ = 8.8 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$  = 11.1 rad m${}^{-2}$⁠) results, respectively. It can be seen that even in the best-case scenario, where the RM is reconstructed with high accuracy at the location of the lensed images, departures from the truth are possible outside this area, indicating, as expected, that we have no constraint there. These results highlight the advantage of using resolved lensed sources over unresolved ones. In general, a combination of the underlying true RM distribution and the lensing configuration sets the constraints on the recovered RM. For example, data sets D${}_2$ and D${}_4$ are generated with the same level of source polarization fraction and share the same underlying model for $\mathit{B}$ and ${\mathit{n}}_{\mathbf{e}}$⁠. However, the former leads to better results (⁠${\sigma }_{\mathrm{R}\mathrm{M}}$ = 0.7 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$ = 2.2 rad m${}^{-2}$ versus ${\sigma }_{\mathrm{R}\mathrm{M}}$ = 3.4 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$ = 11.1 rad m${}^{-2}$⁠) because of the different locations of the counter image with respect to the relevant features in the RM map (see Fig. 5). The same effect explains the relative performance between D${}_3$ (⁠${\sigma }_{\mathrm{R}\mathrm{M}}$ = 2.3 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$ = 7.5 rad m${}^{-2}$⁠) and D${}_1$ (⁠${\sigma }_{\mathrm{R}\mathrm{M}}$ = 0.2 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$ = 0.6 rad m${}^{-2}$⁠). On the other hand D${}_5$ (⁠${\sigma }_{\mathrm{R}\mathrm{M}}$ = 2.6 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$ = 3.3 rad m${}^{-2}$⁠) and D${}_7$ (⁠${\sigma }_{\mathrm{R}\mathrm{M}}$ = 2.3 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$ = 3.0 rad m${}^{-2}$⁠), and D${}_6$ (⁠${\sigma }_{\mathrm{R}\mathrm{M}}$ = 8.8 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$ = 11.1 rad m${}^{-2}$⁠) and D${}_8$ (⁠${\sigma }_{\mathrm{R}\mathrm{M}}$ = 8.8 per cent, $\mathrm{\Delta }\mathrm{R}\mathrm{M}$ = 11.1 rad m${}^{-2}$⁠) lead to a similar level of constraints because the less localized RM structure of the 45-degree inclination angle model can be probed equally well by both lensing configurations.

At a fixed inclination angle and lensing configuration, data sets with a larger polarization fraction lead to better results due to their higher SNR, as expected.

5.3 Magnetic field and electron density profiles

We recover the electron density profile with an accuracy between ${\sigma }_{{\mathrm{n}}_{\mathrm{e}}}$ = 5.4 per cent and $\mathrm{\Delta }{\mathrm{n}}_{\mathrm{e}}=2.0\times {10}^{-6}$ cm${}^{-3}$ in the best case (for data set D${}_5$⁠), and ${\sigma }_{{\mathrm{n}}_{\mathrm{e}}}$ = 319 per cent and $\mathrm{\Delta }{\mathrm{n}}_{\mathrm{e}}=8.48\times {10}^{-4}$ cm${}^{-3}$ in the worst case (for data set D${}_2$⁠). Our choice of a Gaussian profile with a scale height that is almost twice as large as the Einstein radius essentially leads to a constant value of the electron density at the location of the lensed images. As a consequence, we find this quantity to be largely insensitive to the lensing configuration and, hence, will likely not be well constrained in general by galaxy-scale gravitational lensing with typical Einstein radii of 0.15 to 1.5 arcsec (or 0.5 to 5 kpc in linear size, for a lens redshift of 0.5). On the other hand, we see a clear trend with the inclination angle; data sets with $i=45$ deg perform consistently better when compared to those with $i=90$ deg. This is due to the magnetic field profile sharply dropping to zero along the vertical direction in the first model and from the degeneracy between the two scale heights $h_{\mathrm{n}\mathrm{e}}$ and h. For an inclination angle of $i=45$ deg and at a fixed lensing configuration, data sets with a higher SNR (i.e. a larger polarization fraction) provide better constraints.

For the line-of-sight magnetic field profile, ${\mathit{B}}_{\mathrm{L}\mathrm{O}\mathrm{S}}$⁠, we recover the input model within the region probed by the lensed images with an accuracy that ranges from ${\sigma }_{{\mathrm{B}}_{\mathrm{L}\mathrm{O}\mathrm{S}}}$ = 2.9 per cent and $\mathrm{\Delta }{\mathrm{B}}_{\mathrm{L}\mathrm{O}\mathrm{S}}$ = 0.3 nG to ${\sigma }_{{\mathrm{B}}_{\mathrm{L}\mathrm{O}\mathrm{S}}}$  = 151 per cent and $\mathrm{\Delta }{\mathrm{B}}_{\mathrm{L}\mathrm{O}\mathrm{S}}$ = 3.0 nG. Similar to the RM profile, the lensing configuration plays a substantial role in the 90-deg inclination angle model, but less so for the 45-deg case. This aligns with the fact that structures in the RM profile come from the magnetic field model and are just rescaled by the almost constant electron density model. For both inclination angles, data sets with a larger polarization fraction lead to a higher accuracy for a given lensing configuration. As above, we find that reliable constraints can only be obtained within the region probed by the lensed images.

We note that a good accuracy in the RM model does not necessarily translate to good quality constraints for the two separate quantities of ${\mathit{n}}_{\mathbf{e}}$ and ${\mathit{B}}_{\mathrm{L}\mathrm{O}\mathrm{S}}$⁠, which is due to the inherent degeneracy arising from equation (8). In this respect, the magnetic field and electron density models are latent parameters that are not directly measurable. Hence, one should be careful when interpreting the inferred RM in the context of real data. This issue is likely to be model dependent, and we will investigate it more broadly in a follow-up study.

5.4 Magnetic field and electron density model degeneracies

In Table 3, we give the mean values and the 68 per cent confidence-level uncertainties derived with MultiNest for all magnetic field and electron density model parameters. The corresponding posterior distributions for data sets D${}_2$ and D${}_7$⁠, our best-recovered models for a disc inclination of respectively 90 and 45 deg, are presented in Figs 6 and 7, with all remaining data sets presented in Figs A3 and A4 of the Appendix.

From the posterior distributions, we see that all the expected degeneracies are recovered. In the case of the RM, which is an integrated quantity, we find that the parameter $A_0$ is underestimated by, on average, 30 per cent for all data sets and by a factor of 2 in the worst case (for data set D${}_6$⁠). This is due to $r_{{}_0}$ and h of the magnetic field model being overestimated by, on average, a factor of 1.6 and 2, respectively, and given the degeneracy we see between $A_0$ and $r_{{}_0}$⁠, as opposed to an issue with the recovered RM (e.g. see Fig. 4).

The relative spacing between the positive and negative RM features is affected by both the pitch p and phase ${\varphi }_0$ angles. While the former is relatively well constrained within about 2 deg from the truth, the latter can be off by up to 70 deg. For the 90-deg inclination angle, we also recover the degeneracy between the magnetic field and electron density scale heights. The latter is underestimated by a factor of 1.4. As discussed already above, the input scale radius of the electron density distribution is about twice as large as the Einstein radius of the lens model, which in combination with our adopted model for ${\mathit{n}}_{\mathbf{e}}$⁠, makes our simulated data sets insensitive to values of $r_{\mathrm{n}\mathrm{e}}>R_{\mathrm{e}\mathrm{i}\mathrm{n}}$⁠; any value above this limit results in the same Faraday rotating screen within the region where it is measured. As a result, we significantly overestimate this parameter by up to a factor of 15. Both the inclination and PAs are well recovered, with at most a difference of 5 to 15 deg from the input model, respectively. As expected, they are both more tightly constrained for inclination angles of 90 deg due to the stronger and much more localized effect of the screen when compared to the 45-deg case. The limit on the tilt angle ${\chi }_0$ is at most 8 deg off from the input value.

We can conclude that the accuracy with which each parameter is inferred depends on the many factors at play, namely the lensing configuration and how it probes different inclination angle models, the SNR of the data and the intrinsic degeneracies, not only between the electron density and the magnetic field models, but also between the parameters within each model. For example, we find that data set D${}_7$ is the one that performs best in terms of recovering the largest number of parameters with best accuracy (between 2 and 13 per cent). This is likely related to the low-inclination angle and the high SNR of that data set. On the other hand, D${}_6$ provides the worst constraints for most parameters, for example, by underestimating $A_0$ by almost a factor of 2 and overestimating h by nearly a factor of 4. This particular data set has the smallest source size and polarization fraction and, therefore, provides the worst constraints in terms of the structure of the lensed images and their SNR. Thus, although the RM can be well reproduced with this method, the individual model parameters of the electron density and magnetic field models may not be without additional information that can break some of the degeneracies at play.

5.5 Background source

We now consider the ability of our method to recover the intrinsic properties of the background source, namely the recovered surface brightness distribution, polarization fraction and polarization angle. We re-iterate here that although we have used a smooth parametric model to describe the background source surface brightness distribution in the Stokes $\mathit{I}$⁠, $\mathit{Q}$⁠, and $\mathit{U}$⁠, our recovered model has no prior information on this, but is instead based on a regularized pixellated model. The recovered source can be used to infer the robustness of the method since the presence of an intervening magneto-ionic medium will result in a different Faraday rotation of each of the lensed images, which will produce a change in the surface brightness in the Stokes $\mathit{Q}$ and $\mathit{U}$ that is inconsistent with lensing. This will result in a mismatch between the observed and modelled surface brightness in the lensed emission, which will be compensated for with an inaccurate model for the background source. Also, measuring the polarization properties of the gravitationally lensed background source will be in itself interesting for resolving the magnetic field distribution within high-redshift sources (e.g. Roo et al. 2024), which is applicable whether there is an effective intervening magneto-ionic medium (low frequencies; LOFAR, SKA and ngVLA) or not (high frequencies; ALMA).

Fig. 8 presents the input and recovered source surface brightness distribution for Stokes $\mathit{I}$⁠, and the polarization fraction and polarization angle for the data sets D${}_2$ and D${}_7$ at 2 and 3 GHz, as an example (results for the other data sets can be found in Figs A5 to A10 of the Appendix). The corresponding lensed surface brightness distributions for Stokes $\mathit{I}$⁠, $\mathit{Q}$⁠, and $\mathit{U}$ can be seen in Fig. 9. We do not show the results for Stokes $\mathit{V}$⁠, but we find it consistent with zero across all frequency channels within the noise. We see from inspecting the residuals presented in Fig. 8 that the method recovers the input source surface brightness very well for both data sets due to the high SNR of the total intensity emission. We find that the absolute polarization fraction and angle are recovered to within 40 and 14 per cent and to within 4 and 3 deg across the recovered source for $s_1$ and $s_2$⁠, respectively. The most significant deviations in the polarization fraction are seen for data set D${}_2$⁠, which has a factor of 4 lower SNR in polarized flux. This again demonstrates how the robustness of the results will be affected by the quality of the data.

We have also calculated the source residual RMS within the region of the highest SNR, where the input Stokes $\mathit{I}$ surface brightness is at least 20 per cent of the peak. As the total intensity is not affected by the Faraday screen, we find that the quality of the reconstruction solely depends on the lensing configuration, with a relative RMS of about 7 and 10 per cent for $s_1$ and $s_2$⁠, respectively. This is likely related to the fact that a larger fraction of $s_1$ is within the quadruply imaged region, providing more constraints for the source reconstruction. The same effect is observed in Fig. 9 for Stokes $\mathit{Q}$ and $\mathit{U}$ at fixed polarization fraction and for the inclination angle of the underlying RM model. A higher level of polarization fraction leads to a higher SNR and, hence, lower RMS residuals for Stokes $\mathit{Q}$ and, in particular, $\mathit{U}$⁠. We note that for a fixed lensing configuration and polarization fraction, the Stokes $\mathit{Q}$ surface brightness is better recovered when the underlying RM model has a lower inclination angle.

6 DISCUSSION

We have developed the first gravitational lens modelling code that can simultaneously model the lensing mass distribution and the effect of an intervening magneto-ionic medium and recover the polarized properties of the background source, all in the native visibility plane of the data. In this section, we discuss the prospects of our method for constraining the magneto-ionic medium of intermediate-redshift galaxies and the future developments that we plan to carry out.

6.1 Probing the magnetic field structure of intermediate redshift galaxies

The magneto-ionic medium of galaxies is extremely challenging to constrain, except for the very nearby Universe, where the resolved polarized emission from different types of galaxies can be mapped (Fletcher et al. 2011; Heald et al. 2016; Mulcahy, Beck & Heald 2017). However, it has been recently shown that gravitational lens systems with a simple point-like lensed source can provide a new channel for such studies at intermediate redshifts (Mao et al. 2017).

The key observables when it comes to polarized emission are the Stokes $\mathit{Q}$ and $\mathit{U}$ surface brightnesses, from which the electric-vector polarization angle can be determined, which if known as a function of frequency allows the RM to be measured; it is this measurement that links the observed emission to the astrophysics via the dependence of the magnetic field strength and the electron density (Reissl et al. 2023). Indeed, from our simulations, we have shown that the observed surface brightness of the lensed images probes different parts of the foreground lens galaxy within a few to a few tens of kpc in projection, allowing a relative change in the magnetic field structure to be detected. Since gravitational lensing conserves the surface brightness of the lensed source, any change in the polarization angle that is seen in the lensed emission must be due to the presence of a magneto-ionic medium along the line-of-sight (in the absence of variability or an axion-like dark matter particle; see below). This is found at high confidence from our simulations when we compared the results from modelling Faraday rotated data with and without the inclusion of a Faraday screen.

Therefore, this methodology, when applied to a large sample of polarized sources that are gravitationally lensed, can potentially detect evidence of such an intervening medium in the first instance and provide a method for selecting and modelling those objects of interest for detailed follow-up at higher SNR and angular resolution. This will be particularly relevant for the large-area surveys to be carried out with, for example, LOFAR2.0 and SKA-MID in the near-term (McKean et al. 2015). These instruments are expected to find a large population of lensed radio sources that are also polarized. However, in addition to detecting evidence for a magneto-ionic medium, we also find that the RM at the location of the lensed emission is well-recovered, with an average error of around 0.6 to 11 rad m${}^{-2}$⁠; this is an interesting result because, not only is it precise, it has been determined through modelling the magneto-ionic medium and the extended lensed emission simultaneously, as opposed to making independent measurements from the different lensed images to infer the former (e.g. Mao et al. 2017).

Our ability to recover the underlying structure of the magneto-ionic medium and separate the contributions of the magnetic field and the electron density profiles was more challenging than measuring the RM. This is a general problem for polarization studies that also use non-lensed objects (O’Sullivan et al. 2023) as there is an inherent degeneracy between these two components. However, we have found that the product of the magnetic field strength and the electron density normalization could be recovered to within a factor of two, which is extremely encouraging given that there is no other method that can probe the inner parts of such galaxies at intermediate redshifts.

We find that the RM profile can be measured at the locations of the lensed images. However, due to the degeneracies that various combinations of magnetic field structure can have, different models can fit the data equally well. This is because the lensed emission still only probes a limited area of the galaxy structure, and since the lens magnification is mainly in a tangential direction (for a close to isothermal mass distribution), there is only limited information in the radial direction for a single-component polarized source. Therefore, the angular scale of any detectable variation in the magnetic field profile will be limited by the resolution of the data and the configuration of the lensed emission. That said, for the models tested here, we do find that the properties of the magnetic field are relatively well-recovered. This motivates further investigation into what types of RM profiles can be constrained for a wide range of magnetic field structures and lensing configurations. This will be the subject of future work. However, we note that combining with other data, for example, deep optical and infrared imaging, which can be used as a prior on the ellipticity and PA of a magneto-ionic disc or its structure via the observed morphology of the spiral arms, we may be able to break some of the degeneracies.

Finally, we are also able to recover the underlying polarization properties of the lensed source, corrected for any propagation effects as the emission passes through the magneto-ionic medium of the foreground lens. Overall, the surface brightness distribution of the lensed source (its size, ellipticity, and PA) are well-recovered since we take advantage of the Stokes $\mathit{I}$ emission to simultaneously constrain the lens model parameters. This means that the polarized emission is also well recovered because any uncertainties due to the lens model parameters are minimized. We find that for a good SNR, the polarization fraction and the electric-vector polarization angle of the background source are recovered to about 15 per cent and 1 deg, respectively, from the ground truth. The similarity of the polarized source with our input model also demonstrates how well the foreground magento-ionic medium in the lens has been modelled with our method. We note that our method can also be applied to study high-redshift galaxies at mm-wavelengths where propagation effects within the lens are thought to be negligible (e.g. Geach et al. 2023).

6.2 Future developments and prospects

There are still a number of improvements that can be made to our current implementation. For example, the method does not include the contribution of small-scale random magnetic fields, which are expected to dominate the ISM of elliptical galaxies (Seta et al. 2021). These currently make up the bulk of the foreground lens galaxy population (however, see Narasimha & Chitre 2008). In this regard, we plan to develop our technique further in a follow-up publication where the effect of a random magnetic field component is included in the model. We can also take into account the physical processes internal to the background source. In principle, one could fit any model to the reconstructed source to infer its intrinsic polarized and magnetic-field properties, similarly to what is commonly done in the analysis of non-lensed objects. In practice, however, this approach is problematic for strong gravitational lensing observations because the pixels are correlated on the source plane, the noise is not well characterized, and the effective resolution changes across the source due to a differential lensing magnification. These limitations could be overcome by introducing physically motivated models into the prior of the pixellated source following the formalism presented by Rizzo et al. (2018).

In this paper, we have assumed that the polarized flux density, polarization fraction, and polarization angle of the background source are constant in time. This is certainly the case for lensed star-forming galaxies (Geach et al. 2023; Chen et al. 2024; Roo et al. 2024), but for lensed AGN, any variability of the background source within the observing period and on scales smaller than the gravitational lensing time-delay (Biggs 2023) can introduce a bias in the inferred properties of the magnetized plasma in the lens as well as those of the source. If the source intrinsic polarization angle changes on time-scales comparable to the time delay, one cannot assume any more that the observed difference in the polarized properties of the multiple images is solely due to the Faraday screen in the lens, leading to an overestimation of the effect of the latter. This effect can, however, be taken into account in the lens modelling phase as the time delay is related to the mass density distribution of the lensing galaxy. We plan to investigate whether source variability can be included in the forward model. At the same time, from the expected large numbers of lensed sources to be found with next-generation instruments, it will be possible to select samples of objects that are not strongly variable in their polarized emission.

We note that a dark-matter field made of axion-like particles also induces a differential rotation of the lensed-image polarization angles (Basu et al. 2021). However, this effect is frequency independent and, therefore, separable from that caused by the lens with observations taken with an appropriate frequency coverage.

For simplicity, our simulated data sets have a single Gaussian source that is polarized in its entirety for two different source sizes. The extent of the polarized emission has no effect on our ability to recover the lens model parameters since these are constrained using the Stokes $\mathit{I}$ channel; this will always have a larger or equal extent and a larger SNR when compared to the Stokes $\mathit{Q}$ and $\mathit{U}$ emission. We do find that our ability to recover the magneto-ionic parameters is affected by the extent of the background source (and its SNR) since a larger source probes a larger part of the foreground lens. We have also limited our focus to one specific, though realistic, model for the lens galaxy (both in terms of the mass distribution and magneto-ionic properties) and data configuration (in terms of frequency coverage and angular resolution). We will provide an investigation of more generalised scenarios in the future. We will also investigate whether the data contains enough information to rank different choices for the magnetic field and electron density models. Finally, we have neglected the contribution of the lens environment. As many lens galaxies (though not all) reside within a galaxy cluster or group, one should expect a non-negligible contribution to the Faraday rotation by the IGM (Greenfield et al. 1985). In the future, the large samples of gravitational lens systems provided by LOFAR2.0 and SKA-MID could allow one to study the Faraday rotation effect as a function of the lens properties and environment. We intend to test this using simulations of SKA-like observations in the future.

7 CONCLUSIONS

Existing and upcoming observing facilities such as LOFAR2.0 and the SKA are characterized by a large sensitivity, field of view and bandwidth, and a large number density of polarized sources. They are expected, therefore, to push the boundaries of our knowledge on magnetic fields in a large variety of astrophysical systems, from the Milky Way and high-redshift radio galaxies to fast radio bursts and the cosmic web (Beck 2007, 2011; Heald et al. 2020). However, the polarized synchrotron emission from galaxies at high redshift (⁠$z\sim$2 and above) will be too faint to be easily detected even with these telescopes, and measuring their magnetic field structure directly will be challenging. Observations of polarized strongly lensed galaxies will allow us to study magnetic fields at cosmological distances that cannot otherwise be probed. At the same time, the effect of Faraday rotation within the lens will also provide a measure of the magneto-ionic properties of the deflector. Similarly to the traditional method of Faraday RM towards background polarized sources, strong gravitational lensing provides a partial view of the magnetic field in the intervening galaxy. That is, one can only probe the Faraday screen of the lens at the location of the lensed images. This is unlike direct observations of galaxies in synchrotron radiation or heated dust, where one gets a view of the full galaxy or at least where there is emission of this type. However, compared to the latter approach, it has the advantage of being a differential measurement and, hence, potentially more robust on an object-by-object basis. Moreover, as most of the lens galaxies will be massive ellipticals, strong gravitational lensing provides us with an avenue to study a population of objects for which magnetic fields are still poorly understood and can potentially provide observational insights into the fluctuation dynamo action (Seta et al. 2021). Uniquely, our method can probe the magneto-ionic properties of two objects at two different redshifts (the lens and the source) from a single observation.

In this paper, we introduced the first fully self-consistent framework to simultaneously constrain the mass density distribution, the magneto-ionic properties of strong gravitational lens galaxies, and the polarized surface brightness distribution of the background sources in the native visibility plane. We tested this modelling technique using simulated data sets with a realistic lens and source at GHz frequencies, different levels of source polarization fraction and different lensing configurations. We found that we can recover the RM from the lens galaxy with an accuracy that is dependent on the lensing configuration, the underlying magnetic field model and the data SNR. We also obtain a reliable reconstruction of the polarized source properties with better performance for higher SNRs of the data and a larger fraction of the source within the diamond caustic. In the future, we will further develop our method to address some of its current limitations and use existing and upcoming data to investigate the nature of magnetic fields at cosmological distances.

ACKNOWLEDGEMENTS

SWN, SV, and DMP have received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement no. 758853). SV thanks the Max Planck Society for support through a Max Planck Lise Meitner Group. JPM acknowledges support from the Netherlands Organization for Scientific Research (NWO) (project no. 629.001.023) and the Chinese Academy of Sciences (CAS) (project no. 114A11KYSB20170054). This work is based on the research supported in part by the National Research Foundation of South Africa (grant number: 128943).

DATA AVAILABILITY

The data used in this paper are available from the corresponding author on request.

REFERENCES

APPENDIX A: APPENDIX

A1 Reconstructions for all data sets

This section displays the modelling results for all the data sets not shown in the main text. In particular, Figs A1 and A2 display the inferred RM, electron density and line-of-sight magnetic field compared to the input models. The Figs A3 and A4 show the posterior distribution for the magnetic field and electron density parameters for the data sets with an inclination angle of 90 and 45 deg, respectively. The reconstructed and input cleaned lensed images and sources are plotted in Figs A5 to A10.