Issue 
A&A
Volume 672, April 2023



Article Number  A67  
Number of page(s)  11  
Section  Astronomical instrumentation  
DOI  https://doi.org/10.1051/00046361/202244573  
Published online  31 March 2023 
Robust analysis of differential Faraday rotation based on interferometric closure observables
^{1}
Dpt. Astronomia i Astrofísica, Universitat de València,
C/ Dr. Moliner 50,
46120
Burjassot,
Spain
email: ealruiz4@uv.es
^{2}
Observatori Astronòmic, Universitat de València,
C/ Cat. José Beltrán 2,
46980
Paterna,
Spain
Received:
22
July
2022
Accepted:
5
March
2023
Polarization calibration of interferometric observations is a costly procedure and in some cases (e.g., due to limited coverage of parallactic angle for the calibrator), it may not be possible to carry out this process. To avoid this worstcase scenario, while extending the possibilities for the exploitation of polarization interferometric observations, the use of a new set of calibrationindependent quantities (closure traces) has been proposed. However, these quantities suffer from some degeneracies, so their use in practical situations may be rather limited. In this paper, we explore the use of closure traces on simulated and real observations, showing that (with the proper selection of fitting parameters) it is possible to retrieve information on the source polarization using only closure traces and to constrain spatially resolved polarization. We carried out the first application of closure traces to the brightness modeling of real data, using the ALMA observations of M87 conducted during the April 2017 EHT campaign and quantifying a gradient in the Faraday rotation along the source structure (the M87 jet). This work opens the possibility to apply similar strategies to observations from any kind of interferometer (with a special focus on VLBI), whereby quantities such as the differential rotation measure or the spatially resolved polarization can be retrieved.
Key words: techniques: interferometric / polarization / radio continuum: general / techniques: polarimetric / methods: data analysis
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Interferometric techniques (e.g., Thompson et al. 2017) have proven essential in obtaining highresolution images of celestial radio sources, since the long baselines involved can enable the synthesis of large apertures. However, to properly analyze the interferometric data, a thoughtful calibration and reduction is required, as both atmospheric effects and instrumental limitations modulate the phase and amplitude of the observed signals, affecting the interferometric observations and the images that are generated from them.
These atmospheric and instrumental effects are especially relevant at high observing frequencies in verylongbaseline interferometry (VLBI), particularly in polarimetric observations, since the resolutions are so high that it is very difficult to find pointlike calibrators. The effect of the structure further complicates the estimation of the instrumental polarization, which is an important limiting factor that increases the difficulty of the calibration at millimeter (mm) and submillimeter (submm) wavelengths (see e.g., Event Horizon Telescope Collaboration 2021). Even though the calibration of these instrumental effects is often difficult, there are observable quantities (socalled “closures”) that encode information about the properties of the observed sources while remaining relatively insensitive to most of the atmospheric and instrumental effects. Closure quantities have been used in the analysis of interferometric observations since almost the beginning of phasesensitive longbaseline astronomical radio interferometry (e.g., Pearson & Readhead 1984) and most of the image reconstruction algorithms rely on them, either indirectly (via “selfcalibration” used as part of the hybrid imaging techniques Readhead & Wilkinson 1978) or directly (as part of the error function used in direct deconvolution algorithms, e.g., Akiyama et al. 2017a,b). A new set of complex closure quantities, dubbed “closure traces,” was recently reported (see Broderick & Pesce 2020). The closure traces are complex closure quantities constructed from the visibility matrices V measured for baselines connecting four antennas {A, B, C, D} (see more details and an analytical study of the closure traces in Appendix A), as follows: (1)
Closure traces are a generalization of the classical closure amplitude and phase (i.e., the latter can be understood as special cases of the former), but with the important addition of being sensitive to the source polarization structure in a way that is independent of any instrumental polarization. Closure traces encode information of the differential polarization throughout the structure of the source. These unique properties of the closure traces serve to testify to the latent potential of these observables to complement and improve the analysis of the polarization structure of distant radio sources. However, even though Broderick & Pesce (2020) defined the closure traces and explored their symmetries and degeneracies, given these aspects are particularly relevant to the demonstration on the insensitivity of the closure traces to arbitrary coherent (i.e., global throughout the source) rotations in the Poincare sphere, a thorough analysis of the potentially codified information within the closure traces of interferometric observations has yet to be carried out.
In this paper, we apply the newlyintroduced closure traces to extract robust information on the polarization structure of the sources and their possible dependence on the observation frequency. The importance of this study lies within the fact that with a method for obtaining the information of the source polarimetry (independent of the instrumental polarization) hidden in the closure traces, we could analyze this information and scientifically explore polarization observations where the parallacticangle coverage of the correlator is insufficient and, as a result, a proper calibration of the data is not possible.
We organize this paper as follows. In Sect. 2, we offer a detailed description of both simulated and observational data used. In Sect. 3, we present a method for extracting polarimetric information from observed radio sources by fitting parametric models to closure traces. We describe how we tested this method on realistic simulations from interferometric observations of the M87 active galactic nucleus (AGN) with the ALMA telescope and real data taken with the ALMA telescope during the EHT observations in the April 2017 campaign. Finally, we show in Sect. 4, the results of our test and discuss the possibility of using closure traces to determine Faraday rotation (FR) from broadband interferometric observations.
2 Observations and simulated data
The ALMA radio telescope (see e.g., Wootten & Thompson 2009; ALMA Partnership 2015) is the most powerful and sensitive telescope in the mm/submm wavelength regime and, consequently, it is a key component of the EHT. This paper concerns, as a test dataset, the ALMA observations conducted during the April 2017 EHT campaign (see Event Horizon Telescope Collaboration 2019). In this campaign, besides recording the signal from the observations in a VLBI backend for its subsequent processing with the rest of the EHT stations, the ALMA correlator computed all the visibilities among the phased antennas.
We used the ALMA data obtained from the observations of the M87 AGN conducted in April 2017, and properly calibrated by the ALMA observatory. To observe this source, ALMA used an array of 33 antennas, of 12m diameters each, with unprojected baseline lengths ranging from 15.1m to 160.7 m. The observations were performed in the ALMA Band 6 (ranging from 212 GHz to 230 GHz). The ALMA heterodyne receivers can sample the radio signal in various parts of the spectrum simultaneously; specifically, the spectrum accessible to ALMA in one observation is separated into two sidebands. Each of these sidebands has two spectral windows (spws) of 2 GHz each. The measurements were taken in four different spws (0, 1, 2, and 3, corresponding to frequencies centered at 213.1 GHz, 215.1 GHz, 227.1 GHz, and 229.1 GHz, respectively) for all epochs of the EHT campaign (see Goddi et al. 2019, for details on the setup of these observations).
In this paper, we focus on the M87* track observed on 2017 April 6, between 00:52 and 08:02 UT (see Goddi et al. 2021, for details). Figure 1 shows the polarization structure of the AGN and the jet of M87, obtained from the ALMA data (see Goddi et al. 2021).
From the brightness distribution of polarized intensity shown in Fig. 1, both the active core of M87 and the two peaks of polarization intensity (knots A and B) can be distinguished. Moreover, the rotation measure (RM) analysis of the M87 ALMA observations described in Goddi et al. (2021) shows a clear detection of FR from the region of the M87 core. On the other hand, the RM synthesis performed across the whole source brightness distribution shows a clear differential FR between the core and the two jet knots, A and B, with the RM at the core shown to be notably larger than those at the knots. In Fig. 2, we show the RM map of M87, obtained from the results presented in AlbentosaRuiz (2020).
Fig. 1 Image of the polarization structure of the AGN in galaxy M87. The blue scale corresponds to the mapping of the linear polarization intensity. The red lines show the EVPA, with lengths proportional to the polarization intensity. The contours quantify the total intensity (Stokes I). The contours are spaced logarithmically between the image peak (1.21 Jy beam^{−1}) and 0.325 mJy beam^{−1}, which corresponds to the rootmeansquared (rms) of the image residuals. The dashed contour marks the first negative contour in total intensity, at 0.325 mJy beam^{−1}. The full width at half maximum (FWHM) of the convolving Gaussian beam is shown at the bottomleft corner. An alternative version of this figure can be found in Goddi et al. (2021). 
Fig. 2 Rotation measure image of the M87 galaxy detected from the ALMA observations analyzed in this paper. The full width at half maximum (FWHM) of the convolving Gaussian beam is shown at the bottomright corner. 
2.1 Simulations and model parameters
In order to investigate how the polarization structure of a M87like source (and its dependence on frequency) may affect the closure traces, we generate synthetic data of a source with the more relevant components of M87 (i.e., a core and two peaks of polarization intensity) based on the exact baseline coordinates of the observations described in the previous section.
For the simulation of closure traces related to a simplified model of M87 at ALMA scales, we parameterized the polarization intensity as a set of compact (delta) components, associated to the most prominent polarization features of the true image: the core of M87 and the two polarization knots, A and B. The Stokes parameters of each of these components were set to different values for each spectral window (i.e., for each frequency) to simulate the effects of FR in the M87 visibilities. The component positions and Stokes values are given in Table 1. The details of the simulation procedure are described below.
The simulated brightness distribution of Stokes I for each component is taken from the values given by the deconvolved model obtained from the real ALMA data at spw 0. Regarding the rest of Stokes parameters (i.e., Q, U and V), the model consists of three point sources, located at the positions marked in Fig. 2 (which correspond to the jet core and the two knots, A and B). The stokes parameters Q, U and V of the two jet components are taken from the local peak intensity values of the real fullpolarization image (Fig. 1) at each spectral window.
On the other hand, values Q, U, and V of the core component are set manually (see Table 1) to simulate a source with greater EVPA rotation compared to M87, thus increasing the RM of our simulated source to ~0.35 × 10^{7} rad m^{−2}, compared to the value of ~1.5 × 10^{5} rad m^{−2} for M87*, obtained from the ALMA EHT data (see Goddi et al. 2021; AlbentosaRuiz 2021). Doing so allows us to test the performance and reliability of our method, since the high simulated RM values produce clear imprints in the closure traces.
The simulated data are obtained by computing the discrete Fourier transform of the Stokes parameters ((Ĩ, Q̃, Ũ, Ṽ)), for each (u, υ) coordinate and observed frequency v, that is, (2)
where u_{λ} and υ_{λ} are the (u, υ) coordinates in units of wavelength, (I_{i}, Q_{i}, U_{i}, V_{i}) are the Stokes parameters and l_{i}, m_{i} are the direction cosines (see e.g., Thompson et al. 2017) measured with respect to the u and υ axes, that is, the angular offsets with respect the phase center (i.e., the core of M87), in the tangent plane along the right ascension and declination axis, of each of the components of the source.
To conclude our simulation, we introduced, at each integration time, an uncorrelated Gaussian noise contribution to the real and imaginary parts of the visibilities, with zero average and a standard deviation of 0.5 mJy (similar to the average noise levels per integration time seen in the real data).
Simulated model component positions and Stokes values.
2.2 Prior discussion on the behaviour of the closure traces
Before exploring the χ^{2} parametric space with different methods using the closure traces (as proposed in Sect. 3), we analyze the behaviour of the traces for 4 ALMA stations (e.g., DA50DA52DV05DV12, see Fig. 3) for the ALMA observations of the M87 AGN. As the figure shows, both the closuretrace amplitude and phase behaviour differ for the different spws. This is evident when comparing the results for the two sidebands of the spectrum explored by ALMA (i.e., spws {0,1} versus {2,3}).
In short, the closure traces behaviour shows a variation of the amplitude and phase with the observed frequency. This variation can be partially attributed to the FR, since it is, by definition, the magnetooptical phenomenon characterized by the variation of the polarization angle with frequency. However, another notsoobvious cause is the structure factor of the source (i.e., the radial sampling of Fourier space by a baseline, related to the baseline’s reprojection in multifrequency observations). At higher observing frequencies, visibilities are observed at greater distances within the Fourier plane, resulting in a slightly different UV coordinate where the Fourier transform is being measured and, namely, in a slightly different closure trace. Both factors play a part, so this visible appreciation of the variation in the behaviour of the closure traces with frequency might be justified even if the source FR is low. To check the importance of source structure in the spectral dependence of the closure traces, we compare the model visibilities, computed from the Stokes I deconvolved model obtained from the real ALMA data, using Eq. (2) at 213.1 GHz and 229.1 GHz (i.e., spws 0 and 3). We retrieved an average deviation of the model visibilities at spws 0 and 3 of (1.2 ± 0.7)%, much lower than the average deviation computed directly from the visibilities, of (10 ± 20)%. Therefore, the frequency dependence seen in the observed closure traces of M87 may be related to FR.
3 Methodology
Our goal is to study how closure traces can be used to determine the polarimetric structure of sources from interferometric data, that is, we want to quantify the spatial variation of the Stokes parameters across the source brightness distribution. Thus, we need to define a model of the Stokes parameters and compute an error function that can be used to fit the model to the closuretrace data.
Our final objective is to obtain the parameter values of our source polarization model by modelfitting the closure traces only, that is, by searching the set of parameters that minimizes the differences between the closure traces of the simulated or observational data and the parametrized model. We initiate this analysis by fitting a source polarization model to the closure traces obtained from the simulated observations, to assess the correct convergence of (and/or the presence of degeneracies in) the closurebased fitting. Then, we can proceed to the analysis of the observational data.
First, we developed a Python module to compute the closure traces of both simulated and model data. This Python module retrieves the visibilities measured on all the baselines connecting four stations {A, B, C, D} and computes the corresponding closure traces 𝒯_{ABCD}, defined in Eq. (1). The error function used in the leastsquares minimization is simply given by the squared differences between the data closure traces and those from the parameterized model: (3)
We notice, however, that the behaviour of this χ^{2} is highly nonlinear in the parameter space of the chosen model, with many local minima. Therefore, methods based on the gradient (e.g., the Newton method or the more elaborated LevenbergMarquardt; see e.g., Moré 1978) may not converge to the global minimum, depending on the initial parameter values, and the use of methods that explore the complete parametric space, such as Markov chain Monte Carlo (see Goodman & Weare 2010; ForemanMackey et al. 2013; Bussmann et al. 2015) or with a grid search method, are necessary for the purpose of analyzing a greater number of complex sources and obtaining better results.
In VLBI observations, typically taken with arrays of around ten antennas, these methods would be useful for exploring the entire parametric space of the χ^{2} (see e.g., Broderick et al. 2020; Pesce 2021). However, for a large number of antennas, which is the case of the ALMA observations analyzed here, the computational cost may be prohibitive (for a number of n antennas, there are n!/ [4!(n – 4)!] 4antenna baselines). Therefore, in this work we have decided to use a selection of 13 antennas, distributed in such a way that there is not a significant loss of information in the Fourier space (see Fig. 4).
In the next subsections, we discuss the use of two different polarization models: one that only accounts for the core differential polarization and a more complete one that also solves for the polarization in the jet. Each model was fitted using different approaches.
Fig. 3 Closure trace phases and (logarithmic) amplitudes for the observational dataset from the ALMA observations of the M87 AGN, for each observed frequency (spws {0; 1; 2; 3}), for the quadruple of antennas DA50DA52DV05DV12. 
3.1 Modeling the coreonly polarization
The fitting model consists of a simple parameterization of the source structure, where all the Stokes parameters of the jet components are fixed to the measured values in the ALMA image, collected in Table 1, and the only free parameters are related to the core. Since the closure traces are only sensitive to differential (i.e., contrast) changes within the source (in both total intensity and polarization), we have to parameterize the core’s Stokes parameters in a way that encodes differences with respect to the fixed jet components. Otherwise, degeneracies may appear in the fitted parameters. The Stokes parameters of the core are thus formulated using the quantities of the jet components in the following way: (4a) (4b)
where λ_{i} is the wavelength of the spectral window i, λ_{0} is the reference wavelength, and (p_{0}, p_{1}, p_{2}) is the set of parameters needed to characterize the polarization of the core; p_{0} is the polarization angle (the EVPA), p_{1} is the fractional polarization in units of the flux density of the core (i.e., p_{1} = m · I_{core}, where m is the fractional polarization) in Jy, and p_{2} is the RM of the core in units of rad m^{−2}, which accounts for the rotation of the polarization angle with the wavelength (i.e., across the observed bandwidth, divided into four spws).
We fit each spectral window independently by setting the RM to 0, to avoid a direct fit of the RM (which would increase the dimension of the parameter space), reducing the computational cost. Thus, the new model (fitted to each spectral window) is given by Eq. (4) with the parameter p_{2} fixed to 0 (i.e., there is no RM) as we analyzed each spectral window separately. Since this coreonly model only has two fitting parameters, we have opted to perform a grid search method (i.e., the 2D exploration of χ^{2}) focusing on two parameters, the EVPA and the fractional polarization of the core. This method allows us to visualize and understand the behaviour of the χ^{2} for a simple model in a small parametric space. To retrieve the RM, we need to perform a linear fitting of the polarization angle ϕ = ϕ_{0} + RMλ^{2} as a function of the observation wavelength λ (i.e., for each spw).
We explored both parameters (p_{0}, p_{1}), giving different values of the EVPA and the fractional polarization of the core (in units of the core flux density, i.e., p_{1} = m · I_{core}), for each spectral window, obtaining different χ^{2} values, computed with Eq. (4.1.1), for each (p_{0}, p_{1}) value and for each observed frequency. This procedure provides 2D grids where the χ^{2} is shown, adjusted to the traces, as a function of the linear polarization intensity and the EVPA.
From the distribution of the residual closure traces at the minimum of the χ^{2}, it is possible to estimate the correct temperature of the χ^{2} (e.g., Craiu & Rosenthal 2014) that reflects the natural scatter of the data. The probability density is then recovered as
where T is the standard deviation of the closure trace residuals (evaluated for each dataset and spectral window) at the χ^{2} minimum, and assuming the same weights for all visibilities.
Fig. 4 Array of ALMA antennas used in the observations conducted on the April 2017 EHT campaign. The antennas selected for computing the closure traces are highlighted. 
3.2 Modelfitting two M87 components: the MCMC approach
To retrieve the RM directly from the closure traces, we need to include it as a new parameter in the fitting model, by using Eq. (4), without fixing p_{2} to 0. Furthermore, as discussed in Sect. 2, the polarized structure of the M87 jet seen by ALMA presents three relevant components, and therefore our simulation was modeled with three compact components.
Thus, a more accurate modeling procedure would be to allow all polarization components to be characterized by their own RM. However, the closure traces are only sensitive to differential polarization (or polarization contrast) across the source structure, which means that we must leave one of the components fixed (and force its Stokes V to zero, in order to break the latitudinal degeneracy in the Poincaré Sphere) and allow the EVPAs and intensity of the other two components to vary.
A greater number of parameters for our fitting model increases the computational cost of our gridsearch approach. Hence, having studied the gridsearch algorithm, we can go one step further and apply a Markov chain Monte Carlo (MCMC) approach to an extended set of parameters, by including in the fitting model one of the two most polarized jet components. We chose to fix the knot B because it is the jet component with lower RM (see Fig. 2). Therefore, our new fitting model consists of the following parameterization of the source Stokes parameters: (5a) (5b) (5c) (5d)
where (p_{0}, p_{1}, p_{2}) and (p_{3}, p_{4}, p_{5}) are the sets of parameters needed to characterize the polarization of the core and the knot A. Also, p_{0} and p_{3} are the polarization angle (the EVPA) and p_{1} and p_{4} are the fractional polarization in units of the flux density, of the core and the knot A (i.e., p_{1} = m • I, where m is the fractional polarization); then, p_{2} and p_{5} are the RM of the core and the knot A in units of rad m^{−2}.
The MCMC method allows us to test if the closure traces can retrieve the polarization of different components of an observed source, simultaneously, since the gridsearch method becomes computationally unfeasible for the extended set of parameters. We use the Python tool emcee (see ForemanMackey et al. 2013) for a affineinvariant MCMC approach to modelfitting, estimating the differential FR among the source components. Furthermore, as this information is obtained by using closure traces, the results will be independent of the instrumental polarization.
In our modelfitting, we set λ_{0} = 0, and therefore the p_{0} value to which the MCMC method converges corresponds to the intrinsic polarization angle of the core, ϕ_{0}, namely, the EVPA at infinite frequency. If we want to use this result to obtain the EVPA of the core at another reference frequency (e.g., at the average wavelength of the 4 spws observed by ALMA, ), we need to add the neglected FR effect,
We applied this correction to compare the polarization angle to which the MCMC method converges and the EVPA obtained from the image analysis.
An important property of the Stokes components (Q, U) parameterization in Eq. (5) is that the polarization angle is degenerated by a nπ phase, which could hinder the convergence of the MCMC method and give results affected by said degeneracy. We can break this ambiguity using circular statistics (see e.g., Fisher et al. 1987; Jupp & Mardia 2009), namely, a simple but effective change in the parameters of the fitting model, which uses the angle sum identity on Eq. (5),
to redefine the parameterization model as: (6a) (6b) (6c) (6d)
is the knot A auxiliary phase; the new parameters (to determine with the convergence of the MCMC method) are defined, respect the old parameterization, as
These new parameters are unaffected by the EVPA ambiguity, as the ambiguity is absorbed within the parameters of these trigonometric functions (e.g., p_{0} + nπ, where n is an integer, corresponds to the same q_{0}, for a given p_{1}). Consequently, we can retrieve the original parameters without the ambiguities after finding the q_{i} values to which the MCMC method converges as: (7a) (7b)
Finally, p_{2} and p_{5} are preserved in our parameterization as the RM of the core and the knot A.
4 Results and discussion
4.1 Simulated data
4.1.1 Grid search
The exploration of the parametric space of the χ^{2} using the grid search method (described in Sect. 3) provides the value of the polarization angle ϕ and the linear polarization intensity, m, that best fit the closuretrace data (hence, independently of calibration) at each spw. We have, thus, one independent set of polarization parameters per spw, which we can compare to the true Stokes parameters used in the simulation.
In Table 2, we collect the result of our (gridsearch) exploration of the parametric space of the χ^{2} for the four spws of the simulated data. We can see a clear dependence between the polarization angle of the core and the spectral window. We also notice that the fitting model, based only on closure traces, retrieves (with minor discrepancies) the expected EVPA and fractional polarization for spws 0 and 1. On the other hand, even though the results from spws 2 and 3, which correspond to the other sideband of the separated ALMA spectrum, present greater discrepancies compared to the first sideband (especially for the polarization angle), they are still close to the expected values of our simulation. Figure 5 shows the distribution of residual traces at the χ^{2} minimum, that is, at the two fitting parameters of the core (i.e., the EVPA, ϕ, and the fractional polarization, m), which becomes minimized for our fitting model, for spws 0 and 3.
From these results, we retrieve the differential FR along the corejet structure of the simulation by fitting the polarization angle as a function of the observation wavelength ϕ(λ). By doing so, we obtain a value of the RM of the core component, RM = (0.32 ± 0.04) × 10^{7} rad m^{−2}, as shown in Fig. 6, which is compatible with the expected value of (0.35 ± 0.04) × 10^{7} rad m^{−2}. The expected value is computed by fitting the polarization angle of the simulation, determined from the Stokes values of the core given in Table 1, as a function of the observation wavelength ϕ(λ), as seen in Fig. 6.
Furthermore, if we explore the complete χ^{2} parametric space, as shown in Fig. 7, we notice the closure traces degeneracies when exploring each spw independently, as there are several relative minima of the χ^{2} for different values of the polarization angle, ϕ, and the fractional polarization, m. Additionally, we notice that the χ^{2} presents a periodic behaviour in the [−π/2, π/2] interval for the polarization angle, as expected.
4.1.2 Markov chains
At this point, we have analyzed the results for a model where only the polarization of the core is being fitted for each spw independently, retrieving the polarization angle and the fractional polarization and then determining the RM with a linear fitting (as shown in Fig. 6). For this purpose, we used a fitting model with a limited number of parameters, as the computational cost of the gridsearch approach increases for a highdimensional parametric space.
However, by using the MCMC approach (described in Sect. 3.2), we can include not only a parameter to retrieve the RM of the core directly combining the visibilities observed for all the spws, but also allow all polarization components to vary in the MCMC exploration – leaving only one of the components fixed, as the closure traces are only sensitive to differential polarization across the source structure. Therefore, our new fitting model consist of the parameterization of the source core and knot B Stokes Q and U, using the polarization angle, the fractional polarization, and the RM of each component as fitting parameters.
In addition to this, for a correct MCMC exploration of the source polarization and to avoid the effect of EVPA ambiguities, we made use of the circular statistics parameterization (the fitting model described by Eq. (6)). In Fig. 8, we present the results of the MCMC modelfitting, once the polarization angle and the fractional polarization have been retrieved from the circular statistics parameters with Eq. (7).
We obtained a RM posterior distribution spread around RM ~ (0.2–0.4) × 10^{7} rad m^{−2}, which is close to (albeit a bit lower than) the expected value of (0.35 ± 0.04) × 10^{7} rad m^{−2} (see Fig. 6). This result proves that, even though the χ^{2} parametric space presents several minima when modelfitting for each spw individually, as reported in Fig. 7, due to the degeneracies of the closure traces with the EVPA and the fractional polarization, combining all the spw seems to help in breaking the degeneracies, as we retrieved an RM posterior with values close to the true one, within a factor of 2.
Finally, the discrepancies reported for both methods could be explained by the combined effect of introducing Gaussian noise to the simulated visibilities, as well as the use of circular statistics and/or allowing one knot component to vary, and the loss of information from using a subset of 13 antennas in our closurebased modeling (due to the computational limitations), which could also explain the wide width of the RM posterior distribution obtained using the MCMC approach.
Fig. 5 Histograms of the distribution of residual traces at the χ^{2} minimum obtained from the exploration of the χ^{2} parametric space carried out with the grid search method, obtained with the simulated data for spws 0 (left panel) and 3 (right panel). 
Fig. 6 Linear fitting of the polarization angle of the core, ϕ = ϕ_{0} + RMλ^{2}, as a function of the observation wavelength squared, λ^{2}. The ϕ values are determined from the gridsearch exploration of the χ^{2} parametric space, by fitting the closuretrace model to the simulated data, with the polarization of knots A and B fixed to their values at spw0 (see Sect. 2.1). We obtain a fitted value of RM = (0.32 ± 0.04) × 10^{7} rad m^{−2}, compatible with the expected value of (0.35 ± 0.04) × 10^{7} rad m^{−2} obtained from fitting the simulated polarization angles. 
4.2 Real ALMA observations
We have proven that the degeneracy of the closure traces with the EVPA and the fractional polarization observed for each spw is broken when using the MCMC method to quantify the FR (i.e., directly retrieving the RM as a parameter of the fitting model). This motivates the last part of our analysis on the use of closure traces on the real ALMA observations of M87 to retrieve the differential FR. The RM posterior distribution obtained from the MCMC modelfitting approach is shown in Fig. 9.
We go on to compare our results with the results obtained from the image analysis (see Goddi et al. 2021; AlbentosaRuiz 2021). We get a RM posterior distribution peaked around RM ~ (0.8–1.6) × 10^{5} rad m^{−2}, consistent with the values obtained from the image analysis. Therefore, from these results, we can infer that there is indeed a differential FR along the corejet structure of M87 AGN. However, the width of the RM posterior distribution is very wide (i.e., the precision of the RM is poor), likely due to the small EVPA rotation (of the order of 2 deg.) of the core across the ALMA band. In addition to the small EVPA rotation at the core, we used only a subset of 13 antennas in our closurebased modeling (due to the computational limitations), which has a negative impact on the achievable signaltonoise ratio (S/N), as compared to that of the image analysis, where all the phased ALMA antennas have been used. However, for VLBI arrays (consisting typically of 10–20 antennas), it is possible to use the whole interferometer in the closurebased fitting.
If the fractional bandwidth covered by the observations was much larger than that of these ALMA observations (e.g., in the case of the VLBI Global Observing System (VGOS), with practically continuous coverage between 2 and 14 GHz, see Petrachenko et al. 2012; Alef et al. 2019), the precision of the RM could be higher. This would make the use of closure traces for imaging VGOS data a promising task.
The results presented in this work, obtained exclusively and directly using closure traces, are independent of the instrumental effects which contaminate the signal, since the closure traces are insensitive to these effects. We have shown that the use of closure traces on simulated and real observations, with the proper selection of fitting parameters, allows to retrieve information of the source polarization. These results pave the way towards new image reconstruction algorithms that could rely on them, possibly improving their performance, as these quantities are insensitive to all sorts of antennabased and directionindependent calibration effects.
Fig. 7 Results of the exploration of the χ^{2} parametric space carried out with the grid search method, for each of the observed spectral windows. The probability density distribution of the two core fitting parameters (i.e., the EVPA, ϕ, and the fractional polarization, m), corresponding to the χ^{2} scaled by the proper temperature, is shown with six contours located (in peak units) at 0.95, 0.317 (i.e., 1σ), 0.045 (2σ), 0.0027 (3σ), 4σ and 5σ. We note the 180degree ambiguity of the χ^{2} distribution in the direction of the EVPA. 
Fig. 8 Results of the MCMC posterior sampling of the model consisting of three polarized components, fitting two of them (the core and knot B of the simulated M87like source), for the simulated data. The histogram shows the RM of the core, and the dotted line indicates the (higher) expected value, (0.35 ± 0.04) × 10^{7} rad m^{−2}, increased in our simulation, from the image analysis. 
Fig. 9 Results of the MCMC posterior sampling of the model consisting of three polarized components, fitting two of them (the core and knot A of the M87 AGN+jet structure), for the real ALMA observations. The histogram shows the RM of the core, and the dotted line indicates the expected value, of ~1.5 × 10^{5} rad m^{−2}, from the image analysis (see Goddi et al. 2021; AlbentosaRuiz 2021). 
5 Conclusions
Polarization calibration of interferometric observations is a complex procedure prone to different kinds of errors. For the calibration to be successful, strong calibrators, either strongly polarized (in the case of observations with linear feeds) or ideally weakly polarized (in the case of observations with circular feeds) need to be done, covering a wide range of parallactic angles for the participating antennas.
For observations where the parallacticangle coverage of the correlator is insufficient, a proper calibration of the data is not possible. Consequently, rather strong quality assurance (QA) conditions are imposed to polarization observations, resulting in several observing epochs ending up discarded. In this sense, the discovery of calibrationindependent quantities that are sensitive to the source intrinsic polarization (the closure traces, Broderick & Pesce 2020) opens very promising possibilities to develop algorithms able to retrieve robust polarimetric information that will be free from instrumental effects.
In this study, we use a simple leastsquares minimization approach, based on a nondegenerate set of fitting parameters in trace space; namely, focusing on differential (or contrast) polarization changes across the source structure to retrieve sourcepolarization information from a set of ALMA observations. We successfully tested our method with simulated ALMA data and, once it was applied to real observations, we were able to recover the RM of the core of M87 (as observed in ALMA Band 6 during the EHT 2017 campaign) using only closuretrace observables.
Our method can be applied to any kind of interferometer, observing in any kind of polarization basis. Possible uses could include, for instance, studies of FR in AGN jets from VLBI observations, from which longitudinal and transversal RM gradients provide important information about the structure and dynamics of the magnetic fields associated with the jet propagation (e.g., see O’Sullivan & Gabuzda 2009; O’Sullivan et al. 2012).
Acknowledgements
This work has been supported by the grant PRE2020092200 funded by MCIN/AEI/ 10.13039/501100011033 and by ESF invest in your future. This work has been partially supported by the Generalitat Valenciana GenT Project CIDEGENT/2018/021 and by the MICINN Research Project PID2019108995GBC22. We acknowledge support from the Astrophysics and High Energy Physics programme by MCIN, with funding from European Union NextGenerationEU (PRTRC17I1) and the Generalitat Valenciana through grant ASFAE/2022/018. This paper makes use of the following ALMA data: ADS/JAO.ALMA#2016.1.01154.V ALMA is a partnership of ESO (representing its member states), NSF (USA) and NINS (Japan), together with NRC (Canada), MOST and ASIAA (Taiwan), and KASI (Republic of Korea), in cooperation with the Republic of Chile. The Joint ALMA Observatory is operated by ESO, AUI/NRAO and NAOJ.
Appendix A Analytical study of the closure traces
The primary data products obtained by an interferometer are the observed visibilities, as the radio interferometer measurement equation (RIME) (see Hamaker et al. 1996; Hamaker 2000) sets a connection, by means of the Fourier transform, between the observed visibilities and the Stokes maps of astronomical sources. In modern astronomical radio interferometry, one of the greatest difficulties that needs to be solved to analyze the data obtained in the observations is to properly characterize the instrumental gain and polarization terms (also known as Dterms) that contaminate the visibility matrices, modeled by RIME as Jones matrices.
Traditionally, the closure amplitude and closure phase (see e.g., Jennison 1958; Twiss et al. 1960, respectively) have been used to help with the calibration thanks to their invariance to the antenna dependent atmospheric effects and electronic gains. Even though they are completely invariant to antennabased phase errors and atmospherical gains, their usefulness is limited as they are not invariant to nondegenerate calibration errors (e.g., instrumental polarization effects modeled as the Dterms). In this context, Broderick & Pesce (2020) reported a new closure quantity, namely, the closure trace, an interferometric observable characterized for being invariant to the effects of any set of Jones matrices contaminating the visibilities^{1}.
The closure traces are complex closure quantities constructed from the visibility matrices V measured for baselines connecting four antennas {A, B, C, D}, as follows: (A.1)
With this definition, a total of 24 complex traces can be built for the different combinations of possible antennas, with only 6 nonredundant traces (see Broderick & Pesce 2020): 𝒯_{ABCD}, 𝒯_{ABDC}, 𝒯_{ACBD}, 𝒯_{ACDB}, 𝒯_{ADBC} y 𝒯_{ADCB}.
It is on interest to understand, at least qualitatively, how differential variations of the Stokes parameters Q and U throughout the source influences the behaviour of the closure traces. To simplify our analysis, we study analytically the closure traces of the double source (see AlbentosaRuiz 2021).
A.1 Case of the double source
Limiting the problem to one spatial dimension, the polarized structure of the double source, given by the Stokes parameters distribution (I, Q, U, V), can be expressed as a sum of Dirac deltas placing the origin of coordinates at the intermediate point of the separation d of the two double source components. The Fourier transform of the Stokes maps (I, Q, U, V) is, evaluated for the spatial frequency u, in units of the observed wavelength λ, (A.2)
where (I_{i}, Q_{i}, U_{i}, V_{i}) are the Stokes parameters of the double source components 1 and 2.
Fig. A.1 Mapping of the closure trace phase difference and amplitude ratio of a double source, for different values of the Stokes parameters Q_{2} and U_{2}, taking the trace of the double source with Q_{2} = U_{2} = 0 (i.e., unpolarized component 2) as the reference trace. The antennas are distributed in the positions {A; B; C; D} → {0; 1; 4; 6}, in units of wavelength λ for a source separation of 1.3 radian. The white cross in the top panel marks the values of Q_{2} and U_{2} for which we get the maximum closuretrace amplitude ratio. On the other hand, the dashed line in the bottom panel marks the locus of points for which the electric vector position angle (EVPA) of component 1 and component 2 is aligned. 
We studied the behaviour of the closure trace calculated from the visibilities obtained in an interferometric observation with four stations, {A, B, C, D}. For circular feeds, the visibilities are related to the Fourier transform of the Stokes parameters, for a baseline A – B, via (see e.g., Smirnov 2011, and references therein), (A.3)
where u_{AB} is the coordinate of the Fourier plane where the baseline A – B is observing. Plugging these equations in the definition of the trace (Eq. 1), we obtain the analytical expression of the closure trace for the interferometric observation of a double source presented in AlbentosaRuiz (2021).
With a double source, we can analyze how the behaviour of the closure trace varies by modifying the Stokes parameters of one of the components, while keeping the polarization of the other constant. Hence, for a configuration of antennas at positions {A; B; C; D} → {0;1;4;6}, in units of wavelength λ, we fixed the Stokes parameters of component 1, (I_{1}; Q_{1}; U_{1}; V_{1}) = (1; 0.3; 0.15; 0), and we set I_{2} = I_{1} = 1 and V_{2} = V_{1} = 0. In doing so, we can set values to the Stokes parameters (Q_{2}; U_{2}) and compute the closure traces, exploring the complete (Q_{2}; U_{2}) parametric space. Finally, we computed the closure trace ratio, for each (Q_{2}; U_{2}) value, relative to the closure trace of the double source with Q_{2} = U_{2} = 0 (unpolarized component 2). Thus, after generating a synthetic double source and computing the closure traces of a simulated observation, in Fig. A.1 we show, for a certain quadruple of antennas, how both the closure traces phase and amplitude change as a function of Stokes Q_{2} and U_{2}.
With this analytical exploration on the closure traces, we can highlight two conclusions. First, a closure trace amplitude peak (maximum or minimum) is obtained when the EVPAs of the two components are equal. Second, the closure trace phases are more sensitive to polarization intensity variations when the EVPAs of both components point in the same direction. In any case, another interesting conclusion is that the closure traces are actually sensitive to the differential polarization between the two components. Extending this conclusion to the more general case, if we observe frequencydependent closure traces in an observation, it might be an indication of a frequencydependent differential polarization (i.e., a spatiallyresolved depolarization and/or RM). This paper explores this idea, studying how closure traces could allow us to extract information on the polarization structure of an observed source and its frequency dependence directly from the observational data.
References
 Akiyama, K., Ikeda, S., Pleau, M., et al. 2017a, AJ, 153, 159 [NASA ADS] [CrossRef] [Google Scholar]
 Akiyama, K., Kuramochi, K., Ikeda, S., et al. 2017b, ApJ, 838, 1 [NASA ADS] [CrossRef] [Google Scholar]
 AlbentosaRuiz, E. 2020, https://zenodo.org/record/7196017 [Google Scholar]
 AlbentosaRuiz, E. 2021, https://zenodo.org/record/7196001 [Google Scholar]
 Alef, W., Anderson, J. M., Bernhart, S., et al. 2019, in Proceedings of the 24th European VLBI Group for Geodesy and Astrometry Working Meeting, 24, eds. R. Haas, S. GarciaEspada, & J. A. López Fernández, 107 [NASA ADS] [Google Scholar]
 ALMA Partnership (Fomalont, E. B., et al.) 2015, ApJ, 808, L1 [Google Scholar]
 Broderick, A. E., & Pesce, D. W. 2020, ApJ, 904, 126 [NASA ADS] [CrossRef] [Google Scholar]
 Broderick, A. E., Gold, R., Karami, M., et al. 2020, ApJ, 897, 139 [Google Scholar]
 Bussmann, R. S., Riechers, D., Fialkov, A., et al. 2015, ApJ, 812, 43 [Google Scholar]
 Craiu, R. V., & Rosenthal, J. S. 2014, Annu. Rev. Stat. Appl., 1, 179 [NASA ADS] [CrossRef] [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2019, ApJ, 875, L1 [Google Scholar]
 Event Horizon Telescope Collaboration (Akiyama, K., et al.) 2021, ApJ, 910, L12 [Google Scholar]
 Fisher, N. I., Lewis, T., & Embleton, B. J. J. 1987, Statistical Analysis of Spherical Data (Cambridge University Press) [CrossRef] [Google Scholar]
 ForemanMackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013, PASP, 125, 306 [Google Scholar]
 Goddi, C., MartíVidal, I., Messias, H., et al. 2019, PASP, 131, 075003 [Google Scholar]
 Goddi, C., MartíVidal, I., Messias, H., et al. 2021, ApJ, 910, L14 [NASA ADS] [CrossRef] [Google Scholar]
 Goodman, J., & Weare, J. 2010, Commun. Appl. Math. Comput. Sci., 5, 65 [Google Scholar]
 Hamaker, J. P. 2000, Astrom. Astrophys. Suppl., 143, 515 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Hamaker, J. P., Bregman, J. D., & Sault, R. J. 1996, Astrom. Astrophys. Suppl., 117, 137 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Jennison, R. C. 1958, MNRAS, 118, 276 [Google Scholar]
 Jupp, P., & Mardia, K. 2009, Directional Statistics, Wiley Series in Probability and Statistics (Wiley) [Google Scholar]
 Moré, J. J. 1978, in Lecture Notes in Mathematics, 630 (Berlin: Springer Verlag), 105 [CrossRef] [Google Scholar]
 O’Sullivan, S. P., & Gabuzda, D. C. 2009, MNRAS, 400, 26 [Google Scholar]
 O’Sullivan, S. P., Brown, S., Robishaw, T., et al. 2012, MNRAS, 421, 3300 [CrossRef] [Google Scholar]
 Pearson, T.J., & Readhead, A. C. S. 1984, ARA&A, 22, 97 [NASA ADS] [CrossRef] [Google Scholar]
 Pesce, D. W. 2021, AJ, 161, 178 [NASA ADS] [CrossRef] [Google Scholar]
 Petrachenko, W. T., Niell, A. E., Corey, B. E., et al. 2012, Int. Assoc. Geodesy Symp., 136, 999 [CrossRef] [Google Scholar]
 Readhead, A. C. S., & Wilkinson, P. N. 1978, ApJ, 223, 25 [NASA ADS] [CrossRef] [Google Scholar]
 Smirnov, O. M. 2011, A&A, 527, A106 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Thompson, A., Moran, J., & Swenson, G. 2017, Interferometry and Synthesis in Radio Astronomy, Astronomy and Astrophysics Library (Springer International Publishing) [CrossRef] [Google Scholar]
 Twiss, R. Q., Carter, A. W. L., & Little, A. G. 1960, The Observatory, 80, 153 [NASA ADS] [Google Scholar]
 Wootten, A., & Thompson, A. R. 2009, IEEE Proc., 97, 1463 [NASA ADS] [CrossRef] [Google Scholar]
For detailed description on the Radio Interferometer Measurement Equation (RIME) matrix formalism, which describes the response of a radio interferometer to a signal, (see Hamaker et al. 1996; Hamaker 2000)
All Tables
All Figures
Fig. 1 Image of the polarization structure of the AGN in galaxy M87. The blue scale corresponds to the mapping of the linear polarization intensity. The red lines show the EVPA, with lengths proportional to the polarization intensity. The contours quantify the total intensity (Stokes I). The contours are spaced logarithmically between the image peak (1.21 Jy beam^{−1}) and 0.325 mJy beam^{−1}, which corresponds to the rootmeansquared (rms) of the image residuals. The dashed contour marks the first negative contour in total intensity, at 0.325 mJy beam^{−1}. The full width at half maximum (FWHM) of the convolving Gaussian beam is shown at the bottomleft corner. An alternative version of this figure can be found in Goddi et al. (2021). 

In the text 
Fig. 2 Rotation measure image of the M87 galaxy detected from the ALMA observations analyzed in this paper. The full width at half maximum (FWHM) of the convolving Gaussian beam is shown at the bottomright corner. 

In the text 
Fig. 3 Closure trace phases and (logarithmic) amplitudes for the observational dataset from the ALMA observations of the M87 AGN, for each observed frequency (spws {0; 1; 2; 3}), for the quadruple of antennas DA50DA52DV05DV12. 

In the text 
Fig. 4 Array of ALMA antennas used in the observations conducted on the April 2017 EHT campaign. The antennas selected for computing the closure traces are highlighted. 

In the text 
Fig. 5 Histograms of the distribution of residual traces at the χ^{2} minimum obtained from the exploration of the χ^{2} parametric space carried out with the grid search method, obtained with the simulated data for spws 0 (left panel) and 3 (right panel). 

In the text 
Fig. 6 Linear fitting of the polarization angle of the core, ϕ = ϕ_{0} + RMλ^{2}, as a function of the observation wavelength squared, λ^{2}. The ϕ values are determined from the gridsearch exploration of the χ^{2} parametric space, by fitting the closuretrace model to the simulated data, with the polarization of knots A and B fixed to their values at spw0 (see Sect. 2.1). We obtain a fitted value of RM = (0.32 ± 0.04) × 10^{7} rad m^{−2}, compatible with the expected value of (0.35 ± 0.04) × 10^{7} rad m^{−2} obtained from fitting the simulated polarization angles. 

In the text 
Fig. 7 Results of the exploration of the χ^{2} parametric space carried out with the grid search method, for each of the observed spectral windows. The probability density distribution of the two core fitting parameters (i.e., the EVPA, ϕ, and the fractional polarization, m), corresponding to the χ^{2} scaled by the proper temperature, is shown with six contours located (in peak units) at 0.95, 0.317 (i.e., 1σ), 0.045 (2σ), 0.0027 (3σ), 4σ and 5σ. We note the 180degree ambiguity of the χ^{2} distribution in the direction of the EVPA. 

In the text 
Fig. 8 Results of the MCMC posterior sampling of the model consisting of three polarized components, fitting two of them (the core and knot B of the simulated M87like source), for the simulated data. The histogram shows the RM of the core, and the dotted line indicates the (higher) expected value, (0.35 ± 0.04) × 10^{7} rad m^{−2}, increased in our simulation, from the image analysis. 

In the text 
Fig. 9 Results of the MCMC posterior sampling of the model consisting of three polarized components, fitting two of them (the core and knot A of the M87 AGN+jet structure), for the real ALMA observations. The histogram shows the RM of the core, and the dotted line indicates the expected value, of ~1.5 × 10^{5} rad m^{−2}, from the image analysis (see Goddi et al. 2021; AlbentosaRuiz 2021). 

In the text 
Fig. A.1 Mapping of the closure trace phase difference and amplitude ratio of a double source, for different values of the Stokes parameters Q_{2} and U_{2}, taking the trace of the double source with Q_{2} = U_{2} = 0 (i.e., unpolarized component 2) as the reference trace. The antennas are distributed in the positions {A; B; C; D} → {0; 1; 4; 6}, in units of wavelength λ for a source separation of 1.3 radian. The white cross in the top panel marks the values of Q_{2} and U_{2} for which we get the maximum closuretrace amplitude ratio. On the other hand, the dashed line in the bottom panel marks the locus of points for which the electric vector position angle (EVPA) of component 1 and component 2 is aligned. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.