© The Authors 2023

1 Introduction

Interferometric techniques (e.g., Thompson et al. 2017) have proven essential in obtaining high-resolution 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 very-long-baseline interferometry (VLBI), particularly in polarimetric observations, since the resolutions are so high that it is very difficult to find point-like 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 (so-called “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 inter-ferometric observations since almost the beginning of phasesensitive long-baseline astronomical radio interferometry (e.g., Pearson & Readhead 1984) and most of the image reconstruction algorithms rely on them, either indirectly (via “self-calibration” 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 newly-introduced 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 parallactic-angle 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 12-m diameters each, with unpro-jected 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 Albentosa-Ruiz (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 root-mean-squared (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 bottom-left 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 bottom-right corner.

2.1 Simulations and model parameters

In order to investigate how the polarization structure of a M87-like 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 full-polarization 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⁷ rad m⁻², compared to the value of ~1.5 × 10⁵ rad m⁻² for M87*, obtained from the ALMA EHT data (see Goddi et al. 2021; Albentosa-Ruiz 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).

2.2 Prior discussion on the behaviour of the closure traces

Before exploring the χ² 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., DA50-DA52-DV05-DV12, see Fig. 3) for the ALMA observations of the M87 AGN. As the figure shows, both the closure-trace 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 magneto-optical phenomenon characterized by the variation of the polarization angle with frequency. However, another not-so-obvious cause is the structure factor of the source (i.e., the radial sampling of Fourier space by a baseline, related to the baseline’s re-projection in multi-frequency 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 closure-trace data.

Our final objective is to obtain the parameter values of our source polarization model by model-fitting 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 closure-based 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 least-squares 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 χ² is highly non-linear 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 Levenberg-Marquardt; 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; Foreman-Mackey 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 χ² (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)!] 4-antenna 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 DA50-DA52-DV05-DV12.

3.1 Modeling the core-only 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, λ₀ is the reference wavelength, and (p₀, p₁, p₂) is the set of parameters needed to characterize the polarization of the core; p₀ is the polarization angle (the EVPA), p₁ is the fractional polarization in units of the flux density of the core (i.e., p₁ = m · I_core, where m is the fractional polarization) in Jy, and p₂ is the RM of the core in units of rad m⁻², 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₂ fixed to 0 (i.e., there is no RM) as we analyzed each spectral window separately. Since this core-only model only has two fitting parameters, we have opted to perform a grid search method (i.e., the 2D exploration of χ²) 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 χ² for a simple model in a small parametric space. To retrieve the RM, we need to perform a linear fitting of the polarization angle ϕ = ϕ₀ + RMλ² as a function of the observation wavelength λ (i.e., for each spw).

We explored both parameters (p₀, p₁), giving different values of the EVPA and the fractional polarization of the core (in units of the core flux density, i.e., p₁ = m · I_core), for each spectral window, obtaining different χ² values, computed with Eq. (4.1.1), for each (p₀, p₁) value and for each observed frequency. This procedure provides 2D grids where the χ² 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 χ², it is possible to estimate the correct temperature of the χ² (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 χ² 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 Model-fitting 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₂ 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 grid-search approach. Hence, having studied the grid-search 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₀, p₁, p₂) and (p₃, p₄, p₅) are the sets of parameters needed to characterize the polarization of the core and the knot A. Also, p₀ and p₃ are the polarization angle (the EVPA) and p₁ and p₄ are the fractional polarization in units of the flux density, of the core and the knot A (i.e., p₁ = m • I, where m is the fractional polarization); then, p₂ and p₅ are the RM of the core and the knot A in units of rad m⁻².

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 grid-search method becomes computationally unfeasible for the extended set of parameters. We use the Python tool emcee (see Foreman-Mackey et al. 2013) for a affine-invariant MCMC approach to model-fitting, 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 model-fitting, we set λ₀ = 0, and therefore the p₀ value to which the MCMC method converges corresponds to the intrinsic polarization angle of the core, ϕ₀, 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)

where and , and

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₀ + nπ, where n is an integer, corresponds to the same q₀, for a given p₁). 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₂ and p₅ 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 χ² 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 closure-trace 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 (grid-search) exploration of the parametric space of the χ² 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 χ² 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 core-jet 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⁷ rad m⁻², as shown in Fig. 6, which is compatible with the expected value of (0.35 ± 0.04) × 10⁷ rad m⁻². 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 χ² 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 χ² for different values of the polarization angle, ϕ, and the fractional polarization, m. Additionally, we notice that the χ² 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 grid-search approach increases for a high-dimensional 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 model-fitting, 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⁷ rad m⁻², which is close to (albeit a bit lower than) the expected value of (0.35 ± 0.04) × 10⁷ rad m⁻² (see Fig. 6). This result proves that, even though the χ² parametric space presents several minima when model-fitting 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 closure-based 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 grid-search exploration of the χ2 parametric space, by fitting the closure-trace 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) × 107 rad m−2, compatible with the expected value of (0.35 ± 0.04) × 107 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 model-fitting 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; Albentosa-Ruiz 2021). We get a RM posterior distribution peaked around RM ~ (0.8–1.6) × 10⁵ rad m⁻², 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 core-jet 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 closure-based modeling (due to the computational limitations), which has a negative impact on the achievable signal-to-noise 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 closure-based 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 antenna-based and direction-independent 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 180-degree 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 M87-like 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) × 107 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 × 105 rad m−2, from the image analysis (see Goddi et al. 2021; Albentosa-Ruiz 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 parallactic-angle 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 calibration-independent 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 least-squares minimization approach, based on a non-degenerate set of fitting parameters in trace space; namely, focusing on differential (or contrast) polarization changes across the source structure to retrieve source-polarization 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 closure-trace 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 PRE2020-092200 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 PID2019-108995GB-C22. We acknowledge support from the Astrophysics and High Energy Physics programme by MCIN, with funding from European Union NextGenerationEU (PRTR-C17I1) 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 D-terms) 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 antenna-based phase errors and atmospherical gains, their usefulness is limited as they are not invariant to non-degenerate calibration errors (e.g., instrumental polarization effects modeled as the D-terms). 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¹.

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 non-redundant 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 Albentosa-Ruiz 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 Q2 and U2, taking the trace of the double source with Q2 = U2 = 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 wave-length λ for a source separation of 1.3 radian. The white cross in the top panel marks the values of Q2 and U2 for which we get the maximum closure-trace 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 Albentosa-Ruiz (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₁; Q₁; U₁; V₁) = (1; 0.3; 0.15; 0), and we set I₂ = I₁ = 1 and V₂ = V₁ = 0. In doing so, we can set values to the Stokes parameters (Q₂; U₂) and compute the closure traces, exploring the complete (Q₂; U₂) parametric space. Finally, we computed the closure trace ratio, for each (Q₂; U₂) value, relative to the closure trace of the double source with Q₂ = U₂ = 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₂ and U₂.

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 frequency-dependent closure traces in an observation, it might be an indication of a frequency-dependent differential polarization (i.e., a spatially-resolved 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

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 root-mean-squared (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 bottom-left 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 bottom-right 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 DA50-DA52-DV05-DV12.
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 grid-search exploration of the χ2 parametric space, by fitting the closure-trace 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) × 107 rad m−2, compatible with the expected value of (0.35 ± 0.04) × 107 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 180-degree 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 M87-like 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) × 107 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 × 105 rad m−2, from the image analysis (see Goddi et al. 2021; Albentosa-Ruiz 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 Q2 and U2, taking the trace of the double source with Q2 = U2 = 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 wave-length λ for a source separation of 1.3 radian. The white cross in the top panel marks the values of Q2 and U2 for which we get the maximum closure-trace 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