Polarimetric properties of Event Horizon Telescope targets from ALMA

We present the results from a full polarization study carried out with ALMA during the first VLBI campaign, which was conducted in Apr 2017 in the $\lambda$3mm and $\lambda$1.3mm bands, in concert with the Global mm-VLBI Array (GMVA) and the Event Horizon Telescope (EHT), respectively. We determine the polarization and Faraday properties of all VLBI targets, including Sgr A*, M87, and a dozen radio-loud AGN. We detect high linear polarization fractions (2-15%) and large rotation measures (RM $>10^{3.3}-10^{5.5}$ rad m$^{-2}$). For Sgr A* we report a mean RM of $(-4.2\pm0.3) \times10^5$ rad m$^{-2}$ at 1.3 mm, consistent with measurements over the past decade, and, for the first time, an RM of $(-2.1\pm0.1) \times10^5$ rad m$^{-2}$ at 3 mm, suggesting that about half of the Faraday rotation at 1.3 mm may occur between the 3 mm photosphere and the 1.3 mm source. We also report the first unambiguous measurement of RM toward the M87 nucleus at mm wavelengths, which undergoes significant changes in magnitude and sign reversals on a one year time-scale, spanning the range from -1.2 to 0.3 $\times\,10^5$ rad m$^{-2}$ at 3 mm and -4.1 to 1.5 $\times\,10^5$ rad m$^{-2}$ at 1.3 mm. Given this time variability, we argue that, unlike the case of Sgr A*, the RM in M87 does not provide an accurate estimate of the mass accretion rate onto the black hole. We put forward a two-component model, comprised of a variable compact region and a static extended region, that can simultaneously explain the polarimetric properties observed by both the EHT and ALMA. These measurements provide critical constraints for the calibration, analysis, and interpretation of simultaneously obtained VLBI data with the EHT and GMVA.

comprised of a variable compact region and a static extended region, that can simultaneously explain the polarimetric properties observed by both the EHT (on horizon scales) and ALMA (which observes the combined emission from both components). These measurements provide critical constraints for the calibration, analysis, and interpretation of simultaneously obtained VLBI data with the EHT and GMVA.
1. INTRODUCTION Active galactic nuclei (AGN) are known to host supermassive black holes (SMBHs), which accrete gas through a disk and drive powerful relativistic jets that are observed on scales of parsecs to megaparsecs (Blandford et al. 2019). Magnetic fields are believed to play a major role in the formation of such relativistic jets, by either extracting energy from a spinning SMBH via the Blandford-Znajek mechanism (Blandford & Znajek 1977) or by tapping into the rotational energy of a magnetized accretion flow via the Blandford-Payne mechanism (Blandford & Payne 1982).
Polarization observations are a powerful tool to probe magnetic fields and to understand their role in blackhole mass-accretion and launching and acceleration of relativistic AGN jets. In fact, the radio emission from AGN and their associated jets is thought to be produced by synchrotron processes, and thus they display high intrinsic linear polarization (LP; e.g., Pacholczyk 1970;Trippe et al. 2010;Agudo et al. 2018). LP fractions and polarization vector orientations can provide details on the magnetic field strength and topology. Besides LP, circular polarization (CP) may also be present as a consequence of Faraday conversion of the linearly-polarized synchrotron emission (Beckert & Falcke 2002), and can also help constraining the magnetic field configuration (e.g., Muñoz et al. 2012).
As the linearly polarized radiation travels through magnetized plasma, it experiences Faraday rotation of the LP vectors. The externally magnetized plasma is also known as the "Faraday screen" and the amount of Faraday rotation is known as the "rotation measure" (RM). If the background source of polarized emission is entirely behind (and not intermixed with) the Faraday screen, the RM can be written as an integral of the product of the electron number density (n e ) and the magnetic field component along the line-of-sight (B || ) via: Thus, by measuring the RM one can also constrain the electron density, n e , and the magnetic field, B || , in the plasma surrounding SMBHs. Under the assumption that the polarized emission is produced close to the SMBH and then Faraday-rotated in the surrounding accretion flow, the RM has been used in some cases to infer the accretion rate onto SMBHs (e.g, Marrone et al. 2006Marrone et al. , 2007Plambeck et al. 2014;Kuo et al. 2014;Bower et al. 2018). Alternatively, the polarized emission may be Faraday-rotated along the jet boundary layers (e.g., Zavala & Taylor 2004;Martí-Vidal et al. 2015). Therefore, Faraday rotation measurements can provide crucial constraints on magnetized accretion models and jet formation models.
RM studies are typically conducted at cm wavelengths using the Very Large Array (VLA) or the Very Long Baseline Array (VLBA; e.g., Zavala & Taylor 2004). However, cm wavelengths are strongly affected by synchrotron self-absorption (SSA) close to the central engines and can therefore only probe magnetized plasma in the optically thin regions at relatively larger distances (parsec scales) from the SMBH (Gabuzda et al. 2017;Kravchenko et al. 2017). On the other hand, emission at mm wavelengths is optically thin from the innermost regions of the jet base (and accretion disc), enabling us to study the plasma and magnetic fields much closer to the SMBH. In addition, LP can be more easily detected at mm wavelengths because the mm emission region is smaller (e.g., Lobanov 1998), and so depolarization induced by RM variations across the source (e.g., owing to a tangled magnetic field) is less significant. Finally, since Faraday rotation is smaller at shorter wavelengths (with a typical λ 2 dependence), mm-wavelength measurements more clearly reflect the intrinsic LP properties, and therefore the magnetic field of the system. Unfortunately, polarimetric measurements at mm wavelengths have so far been limited by sensitivity and instrumental systematics. The first interferometric measurements of RM at (sub)mm wavelengths were conducted towards Sgr A* with the Berkeley-Illinois-Maryland Association (BIMA) array (Bower et al. 2003(Bower et al. , 2005 and the Submillimeter Array (SMA Marrone et al. 2006Marrone et al. , 2007, which yielded a RM ∼ −5 × 10 5 rad m −2 . SMA measurements towards M87 provided an upper limit |RM| < 7.5×10 5 rad m −2 (Kuo et al. 2014). Other AGN with RM detections with mm interferometers include 3C 84 with RM = 8×10 5 rad m −2 (Plambeck et al. 2014; see also Nagai et al. 2017 for a similarly high RM measured with the VLBA at 43 GHz), PKS 1830-211 (at a redshift z = 2.5) with RM ∼ 10 7 rad m −2 (Martí-Vidal et al. 2015), and 3C 273 with RM = 5 × 10 5 rad m −2 (Hovatta et al. 2019). Additional examples of AGN RM studies with mm single-dish telescopes can be found in Trippe et al. (2012) and Agudo et al. (2018).
In order to progress in this field, polarization interferometric studies at mm wavelengths should be extended to a larger sample of AGN and it will be important to investigate both time and frequency dependence effects, by carrying out observations at multiple frequency bands and epochs. Ultimately, observational studies should be conducted at the highest possible angular resolutions in order to resolve the innermost regions of the accretion flow and/or the base of relativistic jets.
The advent of the Atacama Large Millimeter/ submillimeter Array (ALMA) as a phased array (hereafter phased-ALMA; Matthews et al. 2018;Goddi et al. 2019a) as a new element to Very Long Baseline Interferometry (VLBI) at millimeter (mm) wavelengths (hereafter mm-VLBI) has been a game charger in terms of sensitivity and polarimetric studies. In this paper, we present a complete polarimetric analysis of ALMA observations carried out during the first VLBI campaign.

mm-VLBI with ALMA
The first science observations with phased-ALMA were conducted in April 2017 (Goddi et al. 2019a), in concert with two different VLBI networks: the Global mm-VLBI Array (GMVA) operating at 3 mm wavelength (e.g., Marti-Vidal et al. 2012) and the Event Horizon Telescope (EHT) operating at 1.3 mm wavelength (Event Horizon Telescope Collaboration et al. 2019a). These observations had two "key science" targets, the SMBH candidate at the Galactic center, Sgr A*, and the nucleus of the giant elliptical galaxy M87 in the Virgo cluster, M87*, both enabling studies at horizon-scale resolution (Doeleman et al. 2008(Doeleman et al. , 2012Goddi et al. 2017;Event Horizon Telescope Collaboration et al. 2019b). In addition to those targets, VLBI observations with phased-ALMA also targeted a sample of a dozen radio-loud AGN, including the closest and most luminous quasar 3C 273, the bright γ-ray-emitting blazar 3C 279, the closest radio-loud galaxy Centaurus A (Cen A), and the best supermassive binary black hole candidate OJ287.
In 2019, the first EHT observations with phased-ALMA yielded groundbreaking results, most notably the first ever event-horizon-scale image of the M87* SMBH (Event Horizon Telescope Collaboration et al. 2019b,a,c,d,e,f). Beyond this breakthrough, EHT observations have now imaged polarized emission in the ring surrounding M87*, resolving for the first time the magnetic field structures within a few Schwarzschild radii (R Sch ) of a SMBH (Event Horizon Telescope Collaboration et al. 2021a). In addition, these new polarization images enable us to place tight constraints on physical models of the magnetized accretion flow around the M87* SMBH and, in general, on relativistic jet launching theories (Event Horizon Telescope Collaboration et al. 2021b).
Both the VLBI imaging and the theoretical modelling use constraints from ALMA observations (Event Horizon Telescope Collaboration et al. 2021a,b). In fact, besides providing a huge boost in sensitivity and uvcoverage (Event Horizon Telescope Collaboration et al. 2019c;Goddi et al. 2019b), the inclusion of ALMA in a VLBI array provides another important advantage: standard interferometric visibilities among the ALMA antennas are computed by the ALMA correlator and simultaneously stored in the ALMA archive together with the VLBI recording of the phased signal (Matthews et al. 2018;Goddi et al. 2019a). Furthermore, VLBI observations are always performed in full-polarization mode in order to supply the inputs to the polarization conversion process (from linear to circular) at the VLBI correlators, which is carried out using the Pol-Convert software (Martí-Vidal et al. 2016) after the "Level 2 Quality Assurance" (QA2) process (Goddi et al. 2019a). Therefore, VLBI observations with ALMA yield a full-polarization interferometric dataset, which provides both source-integrated information for refinement and validation of VLBI data calibration (Event Horizon Telescope Collaboration et al. 2021a) as well as observational constraints to theoretical models (Event Horizon Telescope Collaboration et al. 2021b). Besides these applications, this dataset carries valuable scientific value on its own and can be used to derive mm emission, polarization, and Faraday properties of a selected sample of AGN on arcsecond scales.

This paper
In this paper, we present a full polarization study carried out with ALMA in the λ3 mm and λ1 mm bands towards Sgr A*, M87, and a dozen radio-loud AGN, with particular emphasis on their polarization and Faraday properties. The current paper is structured as follows.
Section 3 describes the procedures of data analysis. After presenting some representative total-intensity images of Sgr A* and M87 ( §3.1), two independent algorithms to estimate the Stokes parameters of the compact cores are described ( §3.2). The Stokes parameters for each source and spectral window are then converted into fractional LP and EVPA ( §3.3.1), and used to estimate Faraday rotation ( §3.3.2) and (de)polarization effects ( §3.3.3). Finally, the CP analysis is summarised in §3.3.4.
Section 4 reports the polarimetric and Faraday properties of all the GMVA and EHT target sources, with dedicated subsections on AGN, M87, and Sgr A*.
In Section 5, the polarization properties presented in the previous sections are used to explore potential physical origins of the polarized emission and location of Faraday screens in the context of SMBH accretion and jet formation models. §5.1.1 presents a comparison between the λ3 mm and λ1.3 mm bands, including a discussion on the effects of synchrotron opacity and Faraday rotation; §5.1.2 presents a comparison between the case of blazars and other AGN; §5.1.3 discusses about depolarization in radio galaxies and its possible connection to instrumental effects. §5.2 is devoted to the special case of M87, including a discussion about the origin of the Faraday screen (internal vs. external; §5.2.1) as well as a simple two-component Faraday model ( §5.2.2). Finally, §5.3 is dedicated to the special case of Sgr A*.
Conclusions are drawn in Section 6. The paper is supplemented with a number of appendices including: the list of ALMA projects observed during the VLBI campaign in April 2017 (Appendix A), a full suite of polarimetric images (B) for all the observed targets, comparisons between multiple flux-extraction methods (C) and between the polarimetry results obtained during the VLBI campaign and the monitoring programme with the Atacama Compact Array (D), tables with polarimetric quantities per ALMA spectralwindow (E), Faraday RM plots (F), quality assessment of the circular polarization estimates (G), and mm spectral indices of all the observed targets (H). Finally, a two-component polarization model for M87, which combines constraints from ALMA and EHT observations, is presented in Appendix I. The observations with phased-ALMA were conducted as part of Cycle 4 during the 2017 VLBI campaign in ALMA Band 3 (April 1-3) and Band 6 (April 5-11), respectively. The ALMA data were acquired simultaneously with the VLBI observations (in this sense they are a "byproduct" of the VLBI operations). The ALMA array was in the compact configurations C40-1 (with 0.15 km longest baseline) and, after Apr 6, C40-3 (with 0.46 km longest baseline). Only antennas within a radius of 180 m (from the array center) were used for phasing on all days. About 37 antennas were normally phased together, which is equivalent to a telescope of 73 m diameter 1 . In both Band 3 and 6, the spectral setup includes four spectral windows (SPWs) of 1875 MHz, two in the lower and two in the upper sideband, correlated with 240 channels per SPW (corresponding to a spectral resolution of 7.8125 MHz 2 ). In Band 3 the four SPWs are centered at 86.268, 88.268, 98.328, and 100.268 GHz 3 while in Band 6 they are centered at 213. 100,215.100,227.100,and 229.100 GHz. Three projects were observed in Band 3 with the GMVA (science targets: OJ 287, Sgr A*, 3C 273) and six projects were observed in Band 6 with the EHT (science targets: OJ 287, M87, 3C 279, Sgr A*, NGC 1052, Cen A). The projects were arranged and calibrated in "tracks" (where one track consists of the observations 1 A few more antennas participated in the observations without being phased, the so-called "comparison" antennas, which are mostly used to provide feedback on the efficiency of the phasing process (see Matthews et al. 2018;Goddi et al. 2019a, for details).
2 The recommended continuum setup for standard ALMA observations in full polarization mode is somewhat different and consists of 64 channels, 31.25 MHz wide, per SPW. 3 The "uneven" frequency separation with SPW=2 is due to constraints on the first and second Local Oscillators in the ALMA's tuning system. taken during the same day/session). In Appendix A we provide a list of the observed projects and targets on each day, with the underlying identifications of (calibration and science target) sources within each project (see Tables 4 and 5 in Appendix A). More details of the observation structure and calibration sources can be found in Goddi et al. (2019a).

Data calibration and processing
During phased-array operations, the data path from the antennas to the ALMA correlator is different with respect to standard interferometric operations (Matthews et al. 2018;Goddi et al. 2019a). This makes the calibration of VLBI observations within the Common Astronomy Software Applications (casa) package intrinsically different and some essential modification in the procedures is required with respect to ALMA standard observations. The special steps added to the standard ALMA polarization calibration procedures (e.g., Nagai et al. 2016) are described in detail in Goddi et al. (2019a). The latter focus mostly on the LP calibration and the polarization conversion at the VLBI correlators (Martí-Vidal et al. 2016). In this paper we extend the data analysis also to the CP.
Only sources observed in VLBI mode were calibrated in polarization (see Section 5 in Goddi et al. 2019a). Therefore the sources exclusively observed for ordinary ALMA calibration during the VLBI schedule gaps (i.e., Flux and Gain calibrators) are excluded from this analysis (compare the source list in Tables 4 and 5 in Appendix A with Tables 4 and 6 in Goddi et al. 2019a). Two additional sources observed on Apr 7, 3C 84 and J0006-0623, are also excluded from the following analysis. These sources are in fact flagged in a final flagging step (run on the fully calibrated uv-data before imaging and data analysis), which removes visibility data points having amplitudes outside a certain range (set by three times the RMS from the median of the data) and a source elevation below 25 • . Finally, the two weakest targets observed at 1.3 mm, NCG 1052 and J0132-1654, were found to fall below the flux threshold (correlated flux density of >0.5 Jy on intra-ALMA baselines) required to enable on-source phasing of the array as commissioned (Matthews et al. 2018). Despite these two sources are detected with high signal-to-noise ratio (SNR) in total intensity (SNR> 1000) and polarized flux (SNR> 50 for J0132-1654), we recommend extra care in interpreting these source measurements owing to lower data quality.

Full-stokes imaging
All targets observed in Band 3 and Band 6 are imaged using the casa task tclean in all Stokes parameters: I, Q, U , V . A Briggs weighting scheme (Briggs 1995) is adopted with a robust parameter of 0.5, and a cleaning gain of 0.1. A first quick cleaning (100 iterations over all four Stokes parameters) is done in the inner 10 and 4 in bands 3 and 6, respectively. Providing there is still significant emission (> 7σ) in the residual maps (e.g., in M87 and Sgr A*), an automatic script changes the cleaning mask accordingly, and a second, deeper cleaning is done down to 2σ (these two clean steps are run with parameter interactive=False). A final interactive clean step (with interactive=True) is run to adjust the mask to include real emission which was missed by the automatic masking and to clean deeper sources with complex structure and high-signal residuals (this step was essential for proper cleaning of Sgr A*). No self-calibration was attempted during the imaging stage (the default calibration scheme for ALMA-VLBI data already relies on self-calibration; see Goddi et al. 2019a for details).
We produced maps of size 256×256 pixels, with a pixel size of 0. 5 and 0. 2 in Band 3 and Band 6, respectively, resulting in maps with a field-of-view (FOV) of 128 × 128 and 51 × 51 , respectively, thereby comfortably covering the primary beams of ALMA Band 3 (60 ) and Band 6 (27 ) antennas. We produced maps for individual SPWs and by combining SPWs in each sideband (SPW=0,1 and SPW=2,3), setting the tclean parameters deconvolver='hogbom' and nterms=1, as well as by combining all four SPWs, using deconvolver='mtmfs' and nterms=2. The latter achieved better sensitivity and yielded higher quality images 4 , so we used the combined SPW images for the imaging analysis presented in this paper (except for the per-SPW analysis).
Representative total-intensity images in Band 3 and Band 6 are shown in Figures 1 (Stokes I) and 2 (Stokes I + polarized intensity), whereas the full suite of images including each source observed in Band 3 and Band 6 on each day of the 2017 VLBI campaign is reported in Appendix B (Figures 8-13).
The array configurations employed during phasedarray observations yielded synthesized beams in the range [4. 7-6. 1] × [2. 4-3. 4] in Band 3 and [1. 2-3. 0] × [0. 7-1. 5] in Band 6 (depending on the day and the target). Images on different days achieve different sensitivities and angular resolutions, depending on the time on-source and baseline lengths of the phased-array. In particular, the relatively large range of beamsizes in Band 6 is due to the fact that, during the EHT campaign, progressively more antennas were moved out from the "central cluster" (with a diameter < 150 m). As a consequence, in the last day of the campaign (Apr 11) the observations were carried out with a more extended 4 The deconvolver='mtmfs' performed best when combining all four SPW, yielding on average 30-40% better sensitivity than deconvolver='hogbom' combining two SPW at a time, as expected for RMS ∝ 1/ √ ∆ν. However, deconvolver='hogbom' performed poorly when combining all four SPW, especially for steep spectral index sources, yielding up to 50% worse RMS than deconvolver='mtmfs'. array, yielding a beamsize in the range [1. 2-1. 5] × [0. 7-0. 9] (i.e., an angular resolution roughly two times better than that of other tracks). Tables 6 and 7 in Appendix B report the synthesized beamsize and the RMS achieved in the images of each Stokes parameter for each source observed in Band 3 and Band 6 on each day.

Additional ALMA polarization datasets on M87
In addition to the April 2017 data, we have also analysed ALMA data acquired during the 2018 VLBI campaign as well as ALMA archival polarimetric experiments targeting M87.
The 2018 VLBI campaign was conducted as part of Cycle 5 in Band 3 (April 12-17) and Band 6 (April 18-29), respectively. The observational setup was the same as in Cycle 4, as outlined in §2.1 (a full description of the 2018 VLBI campaign will be reported elsewhere). Three observations of M87 at λ1.3 mm were conducted on Apr 21, 22, and 25 under the project 2017.1.00841.V. For the data processing and calibration, we followed the same procedure used for the 2017 observations, as outlined in §2.2.
The archival experiments include three observations at λ3mm carried out on Sep and Nov 2015(project codes: 2013.1.01022.S and 2015 and Oct 2016 (project code: 2016.1.00415.S), and one observation at λ1.3mm from Sep 2018 (project code: 2017.1.00608.S). For projects 2013.1.01022.S and 2015.1.01170.S, we used directly the imaging products released with the standard QA2 process and publicly available for download from the ALMA archive. For projects 2016.1.00415.S and 2017.1.00608.S, we downloaded the raw visibility data and the QA2 calibration products from the ALMA archive, and we revised the polarization calibration after additional data flagging, following the procedures outlined in Nagai et al. (2016).
The data imaging was performed following the same procedures outlined in §2.3. After imaging, we found that in 2017.1.00608.S, Stokes I, Q, and U are not colocated: U is shifted ∼0.07 to the East, while Q is shifted ∼0.13 West and ∼0.07 north, with respect to I, respectively. This shift (whose origin is unknown) prevents us to assess reliably the polarimetric properties of M87. Therefore, we will not use 2017.1.00608.S in the analysis presented in this paper. The analysis and results of the other datasets will be presented in §4.2.
3. DATA ANALYSIS 3.1. Representative total intensity images The sources targeted by the GMVA and EHT are generally unresolved at arcsecond scales and their images are mostly consistent with point sources (see images displayed in Appendix B). The EHT key science targets, Sgr A* and M87, are clear exceptions, and show complex/extended structures across tens of arcseconds. We show representative images of Sgr A* (3 mm, Apr 3; 1.3 mm, Apr 6) and M87 (1.3 mm, Apr 11) in Figure 1.  The image showcases the well-known "mini-spiral" structure surrounding the central compact core, including the eastern and northern arms, the western arc, and the bar at the center. The contour levels at 1.3 mm are 5σ × 2 n where σ =0.44 mJy beam −1 and n= 0 , 1 , 2 , 3 . . . up to the peak flux-density; the contour level at 3 mm corresponds to 20σ (σ =0.8 mJy beam −1 ). The peak flux-density is 2.5 ± 0.1 (2.6 ± 0.3) Jy beam −1 and the integrated flux density across the entire source is 9.9 ± 0.5 (4.9 ± 0.5) Jy at a representative frequency of 93 (221) GHz. The FOV is given by the primary beam in Band 3 (∼60 ) and 1 pc corresponds to 24 . The beamsizes are 5. 0 × 2. 7 (P.A.= -81.1 • ) in Band 3 and 2. 2 × 1. 3 (P.A.= -77.5 • ) in Band 6, shown as a blue and yellow ovals, respectively, in the lower left corner. Right panel: Image of M87 at 1.3 mm on Apr 11 2017. The image showcases the structure of the kpc-scale relativistic jet comprised of a bright core at the nucleus and the knots along the jet labeled as D, F, A, B, C; HST-1 is not resolved from the nucleus in these images. The RMS noise level is 0.16 mJy beam −1 , and the contour levels are a factor of 10 and 40 of the RMS. The peak flux-density is 1.34 Jy beam −1 and the integrated flux-density is 1.57 Jy at a representative frequency of 221 GHz. The FOV is given by the primary beam in Band 6 (∼ 27 at 1.3 mm). 1 kpc corresponds to 12 . This observation was conducted with the most extended array during the VLBI campaign, yielding the highest angular resolution (beam size = 1. 2 × 0. 8, P.A.= 79.3 • , shown in the lower left corner). In both panels the four observing SPWs (see § 2.1) were used together for imaging. The intensity brightness is plotted using a logarithmic weighting function (starting from the 5σ-level), in order to highlight the full extent of both the mini-spiral (in Sgr A*) and the jet (in M87).
The images displayed cover an area corresponding to the primary beam of the ALMA antennas (27 in Band 6 and 60 in Band 3; the correction for the attenuation of the primary beam is not applied to these maps). The Sgr A* images clearly depict the well-known "mini-spiral" structure which traces ionized gas streams surrounding the central compact source; the mini-spiral has been studied in a wide range of wavelengths (e.g. Zhao et al. 2009;Irons et al. 2012;Roche et al. 2018). The "eastern arm", the "northern arm", and the "bar" are clearly seen in both Band 3 and Band 6, while the "western arc" is clearly traced only in the Band 3 image (it falls mostly outside of the antenna primary beam for Band 6). Similar images were obtained in the 100, 250, and 340 GHz bands in ALMA Cycle 0 by Tsuboi et al. (2016, see their Fig. 1). Since Sgr A* shows considerable variability in its core at mm-wavelengths (e.g. Bower et al. 2018), the displayed maps and quoted flux values throughout this paper should be considered as time-averaged images/values at the given epoch.
The M87 jet has been observed across the entire electromagnetic spectrum (e.g., Prieto et al. 2016), and imaged in detail at radio wavelengths from λ1 metre (with LOFAR: de Gasperin et al. 2012) through λ[15-0.7] cm (with the VLA and the VLBA: e.g., Hada et al. 2013;Walker et al. 2018) up to λ 3 mm (with the GMVA: e.g., Kim et al. 2018). VLA images at lower radio frequencies (e.g., Biretta et al. 1995) showcase a bright component at the nucleus and a kiloparsec-scale (kpcscale) relativistic jet, extending across approximately 25 (∼ 2 kpc) from the central core. Images of the kpc-scale relativistic jet were also produced with ALMA Cycle-0 observations at λ 3 mm ) and with the SMA at λ 1 mm (Tan et al. 2008;Kuo et al. 2014), but could only recover the bright central core and the strongest knots along the jet.
Our λ1.3 mm ALMA image showcases a similar structure, but the higher dynamic range (when compared with these earlier studies) allows us to recover the continuous structure of the straight and narrow kpc-scale jet across approximately 25 from the nucleus, including knots D, F, A, B, C, at increasing distance from the central core (HST-1 is not resolved from the nucleus in these images). The jet structure at larger radii ( 2 kpc) as well as the jet-inflated radio lobes, imaged in great detail with observations at lower frequencies, are not recovered in our images (see for example the NRAO 20 cm VLA image).

Extracting Stokes parameters in the compact cores
We extract flux values for Stokes I, Q, U , and V in the compact cores of each target observed in Band 3 and Band 6. We employ three different methods which use both the visibility data and the full-Stokes images. In the uv-plane analysis, we use the external casa library uvmultifit (Martí-Vidal et al. 2014). To reduce its processing time, we first average all (240) frequency channels to obtain one-channel four-SPW visibility uvfiles. We assume that the emission is dominated by a central point source at the phase centre and we fit a delta function to the visibilities to obtain Stokes I, Q, U , and V parameters in each individual SPW. Uncertainties are assessed with Monte Carlo (MC) simulations, as the standard deviation of 1000 MC simulations for each Stokes parameter. For the image-based values, we take the sum of the central nine pixels of the CLEAN model component map (an area of 3 × 3 pixels, where the pixel size is 0. 2 in Band 6 and 0. 5 in Band 3). Summing only the central pixels in the model maps allows one to isolate the core emission from the surroundings in sources with extended structure. A third independent method provides the integrated flux by fitting a Gaussian model to the compact source at the phase-center in each image with the casa task IMFIT. In the remaining of this paper, we will indicate these three methods as uvmf, 3x3, and intf.
From a statistical perspective, any fitting method in the visibility domain should be statistically more reliable than a χ 2 -based fitting analysis in the image-plane (whose pixels have correlated noise), and should therefore be preferred to image-based methods. However, we have two reasons for considering both approaches in this study: (i) some of our targets exhibit prominent emission structure at arc-second scales (see Figure 1, and the maps in Appendix B), (ii) the observations are carried out with various array configurations, resulting in a different degree of filtering of the source extended emission. Both elements can potentially bias the flux values of the compact cores extracted in the visibility domain vs. the image domain.
In Appendix C we present a comparative analysis of three flux-extraction methods to assess the magnitude of such systematic biases (reported in Table 8). The statistical analysis shows that the Stokes I values estimated with uvmf are consistent with those estimated from the images, with a median absolute deviation (MAD) ≤ 0.07% and individual offsets <1% (for both point sources and extended sources) in the case of the 3x3 method (the agreement is slightly worse for the intf method). These deviations are negligible when compared to the absolute uncertainty of ALMA's flux calibration (10% in Band 6). This consistency generally holds also for Stokes Q and U (with MAD < 1%) and other derived parameters within their uncertainties (see Table 8). We therefore conclude that, for the purpose of the polarimetric analysis conducted in this paper, the uv-fitting method uvmf provides sufficiently precise flux values for the Stokes parameters (but see Appendix C for details on M87 and Sgr A*). Goddi et al. (2019a) report the Stokes I flux values per source estimated in the uv-plane from amplitude gains using the casa task fluxscale. We assess that the Stokes I estimated from the visibilities with uvmf are consistent with those estimated with fluxscale generally within 1%. In addition, Goddi et al. (2019a) compared the fluxscale flux values (after opacity correction) with the predicted values from the regular flux monitoring programme with the ALMA Compact Array (ACA), showing that these values are generally within 10% (see their Appendix B and their Fig. 16). In Appendix D we perform a similar comparative analysis for the sources commonly observed in the ALMA-VLBI campaign and the AMAPOLA polarimetric Grid Survey, concluding that our polarimetric measurements are generally consistent with historic trends of grid sources (see Appendix D for more details and comparison plots).

Polarimetric data analysis
In this section we use the measured values of the Stokes parameters to determine the polarization properties for all targets, including the fractional LP ( § 3.3.1), the electric vector position angle (EVPA) and its variation as a function of frequency or Faraday Rotation ( § 3.3.2), the degree of depolarization ( § 3.3.3), and the fractional CP ( § 3.3.4). These polarization quantities, averaged across the four SPWs, are reported in Tables 1  and 2 for each target observed with the GMVA and  the EHT, respectively (while Table 3 summarizes all the ALMA polarimetric observations towards M87 analysed in this paper). For selected EHT targets, the polarization properties (per SPW and per day) are displayed in Figure 3.

Linear polarization and EVPA
The values estimated for Stokes Q and U can be combined to directly provide the fractional LP in the form Q 2 + U 2 /I, as well as the EVPA, χ, via the equation 2χ = arctan(U/Q). Tables 9 and 10 in Appendix E report Stokes parameters, LP, and EVPA, for each SPW. The LP has been debiased in order to correct for the LP bias in the low SNR regime (this correction is especially relevant for low polarization sources; see Appendix E for the debiased LP derivation).  The raster image and blue contour show the total intensity emission, the orange contours show the linearly polarized emission, and the black vectors showcase the orientation of the EVPAs (their length is linearly proportional to the polarized flux). The total intensity brightness is plotted using a logarithmic weighting function (starting from the 1σ-level), the blue contour corresponds to 5σ (where σ is the Stokes I map RMS), while the orange contour levels are 5σ × 2 n (where σ is the LP map RMS and n=0,1,2,3... up to the peak in the image). The LP fraction at the peak of the compact core is reported in the upper left corner in each panel. The EVPAs are plotted every 8 pixels (1. 6 or about 1 per beam) for Sgr A* and every 4 pixels (0. 8 or about 2 per beam) for M87 (in order to sample more uniformly the jet). According to the measured RM, the EVPAs towards the compact core should be rotated by −23 • (east of north) in Sgr A* and by −16 • in M87. The beamsizes (shown as an oval in the lower left corner) are 2. 2 × 1. 3 (P.A. -77 • ) and 2. 2×1. 5 (P.A. -69 • ) in the left and right panels, respectively. Note that there are several tiny EVPAs plotted across the mini-spiral, apparently locating regions with polarized flux above the image RMS noise cutoff (5σ). The LP and EVPA errors are however dominated by the systematic leakage (0.03% of I onto QU), which is not added to the images. Once these systematic errors are added, the LP flux in those points falls below the 3σ measurement threshold. Therefore we do not claim detection of polarized emission outside of the central core in Sgr A*. Besides, only the polarisation within the inner 1/3 of the primary beam is guaranteed by ALMA. The full set of 1.3 mm observations of Sgr A* and M87 are reported in Figures 8  and 9, respectively.     The estimated LP fractions range from 0.1 − 0.2% for the most weakly polarized targets (Cen A and NGC 1052) to 15% for the most strongly polarized target (3C279), consistent with previous measurements (see Appendix D). The uncertainties in LP include the fitting (thermal) error of Stokes Q and U and the (systematic) Stokes I leakage onto Stokes Q and U (0.03% of Stokes I) added in quadrature. This analysis yields LP uncertainties < 0.1%, similar to those quoted in previous studies (Nagai et al. 2016;Bower et al. 2018). Figure 2 showcases representative polarization images of Sgr A* (left panel) and M87 (right panel) as observed at 1.3 mm on Apr 6. The individual images display the measured EVPAs overlaid on the polarized flux contour images and the total intensity images. Note that the EVPAs are not Faraday-corrected and that the measured 5 magnetic field orientations should be rotated by 90 • . In Sgr A*, polarized emission is present only towards the compact core, while none is observed from the mini-spiral. In M87, the EVPA distribution appears quite smooth along the jet, with no evident large fluctuations of the EVPAs in nearby regions, except between Knots A and B. For a negligible RM along the jet, one can infer that the magnetic field orientation is first parallel to the jet axis, then in Knot A it changes direction (tending to be perpendicular to the jet), and then turns back to be parallel in Knot B, and finally becomes perpendicular to the jet axis further downstream (Knot C). This behaviour can be explained if Knot A is a standing or recollimation shock: if multiple standing shocks with different magnetic field configurations form along the jet and the latter is threaded with a helical magnetic field, its helicity (or magnetic pitch) would be different before and after the shock owing to a different radial dependence of the poloidal and toroidal components of the magnetic field (e.g., Mizuno et al. 2015). The EVPA distribution is also in good agreement with the polarization characteristics derived from observations at centimeter wavelengths with the VLA (e.g., Algaba et al. 2016). We nevertheless explicitly note that only the polarisation within the inner third of the primary beam is guaranteed by the ALMA observatory. Since we focus on the polarization properties in the core, the analysis presented in this paper is not affected by this systematics.

Rotation measure
Measuring the EVPA for each SPW (i.e., at four different frequencies) enables us to estimate the RM in the 3 mm band (spanning a 16 GHz frequency range of 85-101 GHz) and in the 1.3 mm band (spanning a 18 GHz frequency range of 212-230 GHz), respectively.
In the simplest assumption that the Faraday rotation is caused by a single external Faraday screen (i.e., it occurs outside of the plasma responsible for the polarized emission), a linear dependence is expected between the EVPA and the wavelength squared. In particular, we fit the RM and the mean-wavelength (λ) EVPA (χ) following the relation: where χ is the observed EVPA at wavelength λ andχ is the EVPA at wavelengthλ. The EVPA extrapolated to zero wavelength (assuming that the λ 2 relation holds) is: The RM fitting is done using a weighted least-squares method of χ against λ 2 . Theχ, χ 0 , and the fitted RM values are reported in the sixth, seventh, and eighth columns of Tables 1 and 2, respectively. The EVPA uncertainties quoted in Tables 1, 2, 3, 9, 10, are typically dominated by the systematic leakage of 0.03% of Stokes I into Stokes Q and U . At 1.3 mm, this results in estimated errors between 0.06 • for the most strongly polarized source (3C279) and 0.8 • for the weakest source (J0132-1654), with most sources in the range 0.1 • -0.4 • . These EVPA uncertainties imply RM propagated errors between 0.4 × 10 4 rad m −2 and 6×10 4 rad m −2 , with most sources in the range (1−3)× 10 4 rad m −2 . Similarly, at 3 mm we find EVPA uncertainties of 0.07 • -1.4 • , with a typical value of 0.2 • , and RM uncertainties in the range (0.06−1.0)×10 4 rad m −2 , with a typical value of ∼ 0.13 × 10 4 rad m −2 .
In Appendix F, we present plots of the measured EV-PAs at the four ALMA SPWs and their RM fitted models as a function of λ 2 (displayed in Figures 17,18,19,and 20) and we discuss the magnitude of thermal and systematic errors in the RM analysis.

Bandwidth depolarization
In the presence of high RM, the large EVPA rotation within the observing frequency bandwidth will decrease the measured fractional polarization owing to Faraday frequency or "bandwidth" depolarization, which depends on the observing frequency band. The RM values inferred in this study (e.g., Table 2) introduce an EVPA rotation of less than one degree within each 2 GHz spectral window, indicating that the bandwidth depolarization in these data should be very low (<0.005%). However, if there is an internal component of Faraday rotation (i.e., the emitting plasma is itself causing the RM), there will be much higher frequency-dependent (de)polarization effects (the "differential" Faraday rotation), which will be related to the structure of the Faraday depth across the source (e.g. Cioffi & Jones 1980;Sokoloff et al. 1998).
We have modeled the frequency dependence of LP using a simple linear model: where m is the observed LP at frequency ν,m is the LP at the mean frequencyν, and m is the change of LP per unit frequency (in GHz −1 ). Given the relatively narrow fractional bandwidth ( 2 GHz), the linear approximation given in Eq. 4 should suffice to model the frequency depolarization (multifrequency broadband single-dish studies fit more complex models; see for example Pasetto et al. 2016Pasetto et al. , 2018. We have fitted the values of m from a least-squares fit of Eq. 4 to all sources and epochs, using LP estimates for each spectral window from Table 10. We show the fitting results for selected sources in Fig.  3 (lower panels). There are clear detections of m for 3C279, Sgr A*, and M87; these detections also differ between epochs. Such complex time-dependent frequency effects in the polarization intensity may be indicative of an internal contribution to the Faraday effects observed at mm wavelengths.

Circular polarization
Measuring Stokes V provides, in principle, a direct estimate of the fractional CP as V /I. In practice, the polarization calibration for ALMA data in CASA is done by solving the polarization equations in the linear approximation, where parallel-hands and cross-hands visibilities are expressed as a linear function of I, Q, U , while it is assumed V = 0 (e.g., Nagai et al. 2016;Goddi et al. 2019a). A non-negligible Stokes V in the polarization calibrator will introduce a spurious instrumental Stokes V into the visibilities of all the other sources. Moreover, such a Stokes V introduces a bias in the estimate of the cross-polarization phase, β, at the reference antenna (see Appendix G), which translates into a leakage-like effect in the polconverted VLBI visibilities (see Eq. 13 in Goddi et al. 2019a). The magnitude of such a bias may depend on the fractional CP of the polarization calibrator, the parallactic-angle coverage of the calibrator, and the specifics of the calibration algorithm. In Appendix G we attempt to estimate such a spurious contribution to Stokes V by computing the cross-hands visibilities of the polarization calibrator as a function of parallactic angle (see Figures 21 and 22). This information can then be used to assess Stokes V and CP for all sources in all days (reported in Tables 11 and 12 for GMVA and EHT sources, respectively).
We stress two main points here. First, the reconstructed Stokes V values of the polarization calibrators are non-negligible and are therefore expected to introduce a residual instrumental X-Y phase difference in all other sources, after QA2 calibration. This can be seen in the dependence of the reconstructed Stokes V with feed angle in almost all the observed sources (displayed in Figure 22). The estimated X-Y residual phase offsets are of the order of 0.5 • , but they can be as high as 2 • (e.g., on Apr 5). These values would translate to a (purely imaginary) leakage term of the order of a few % in the polconverted VLBI visibilities.
The second point is that there is a significant variation in the estimated values of reconstructed stokes V across the observing week. In particular, on Apr 5, 3C279 shows a much higher value, indicating either an intrinsic change in the source, or systematic errors induced by either the instrument or the calibration. In either case, this anomalously large Stokes V in the polarization calibrator introduces a large X-Y phase difference in all other sources. This can be seen in the strong dependence of reconstructed Stokes V on feed angle for sources OJ287 and 4C 01.28 (displayed in Figure 22, upper left panel) and in their relatively high Stokes V when compared to the following days (see Table 12). Besides the anomalous value in Apr 5, it is interesting to note that the data depart from the sinusoidal model described by Eq. G2, for observations far from transit, especially on Apr 11. These deviations may be related to other instrumental effects which however we are not able to precisely quantify. For these reasons, we cannot precisely estimate the magnitude of the true CP fractions for the observed sources (see Appendix G for details). Nevertheless, our analysis still enables us to obtain order-of-magnitude values of CP. In particular, excluding the anomalous Apr 5, we report CP =[-1.0,-1.5] % in Sgr A*, CP∼ 0.3% in M87, and possibly a lower CP level (∼ 0.1 − 0.2%) in a few other AGNs (3C273, OJ287, 4C 01.28, J0132-1654, J0006-0623; see Table 12). In the 3 mm band, we do not detect appreciable CP above 0.1%, except for 4C 09.57 (-0.34%), J0510+1800 (-0.14%) and 3C273 (0.14%). We however note that the official accuracy of CP guaranteed by the ALMA observatory is < 0.6% (1σ) or 1.8% (3σ), and therefore all of these CP measurements should be regarded as tentative detections.
Their polarimetric quantities at 3 mm and 1.3 mm are reported in Tables 1 and 2, respectively, and displayed in Figure 3. Overall, we find LP fractions in the range 0.1-15% (with SN R ∼ 3 − 500σ) and RM in the range 10 3.3 − 10 5.5 rad/m 2 (with SN R ∼ 3 − 50σ), in line with previous studies at mm-wavelengths with singledish telescopes (e.g., Trippe et al. 2010;Agudo et al. 2018) and interferometers (e.g., Plambeck et al. 2014;Martí-Vidal et al. 2015;Hovatta et al. 2019). We also constrain CP to <0.3% in all the observed AGN, consistent with previous single-dish (e.g., Thum et al. 2018) and VLBI (e.g., Homan & Lister 2006) studies, suggesting that at mm-wavelengths AGN are not strongly circularly polarized and/or that Faraday conversion of the linearly polarized synchrotron emission is not an efficient process (but see Vitrishchak et al. 2008).
In Appendix B, we also report maps of all the AGN targets observed at 1.3 mm (Figures 10,11,12), and at 3 mm ( Figure 13), showcasing their arcsecond-structure at mm wavelengths.
In the rest of this section, we briefly comment on the properties of selected AGN.
3C 279 -3C 279 is a bright and highly magnetized gamma-ray emitting blazar, whose jet is inclined at a very small viewing angle ( 3 • ). At its distance (z=0.5362), 1 arcsecond subtends 6.5 kpc. 3C 279 was observed on four days at 1.3 mm and one day at 3 mm. It is remarkably highly polarized both at 1.3 mm and 3 mm. At 1.3 mm, LP varies from 13.2% on Apr 5 to 14.9% in Apr 11, while the EVPA goes from 45 • down to 40 • . At 3 mm, LP is slightly lower (∼ 12.9%) and the EVPA is 44 • .
While at 1.3 mm we can only place a 1σ upper limit of <5000 rad/m 2 , at 3 mm we measure a RM= 1790 ± 460 rad/m 2 (with a ∼ 4σ significance). Lee et al. (2015) used the Korean VLBI Network to measure the LP at 13, 7, and 3.5 mm, finding RM values in the range -650 to -2700 rad m −2 , which appear to scale as a function of wavelength as λ −2.2 . These VLBI measurements are not inconsistent with our 3 mm measurement and our upper limits at 1.3 mm, but more accurate measurements at higher frequencies are needed to confirm an increase of the RM with frequency.
The total intensity images at 1.3 mm reveal, besides the bright core, a jet-like feature extending approximately 5 towards south-west (SW) (Fig. 10); such a feature is not discernible in the lower-resolution 3 mm image (Fig. 13). The jet-like feature is oriented at approximately 40 • , i.e. is roughly aligned with the EVPA in the core. Ultra-high resolution images with the EHT reveal a jet component approximately along the same PA but on angular scales 10 5 times smaller (Kim et al. 2020).
3C 273 was observed both at 1.3 mm and 3 mm (2 days apart). Total intensity and LP are higher in the lower frequency band: F=9.9 Jy and LP=4.0% (at 3 mm) vs. F=7.6 Jy and LP=2.4% (1.3 mm). We estimate a RM = (2.52 ± 0.27) × 10 5 rad/m 2 at 1.3 mm, confirming the high RM revealed in previous ALMA observations (conducted in Dec 2016 with 0. 8 angular resolution) by Hovatta et al. (2019) who report LP=1.8% and a (twice as large) RM= (5.0 ± 0.3) × 10 5 rad/m 2 . We also report for the first time a RM measurement at 3 mm, RM = (−0.60 ± 0.14) × 10 4 rad/m 2 , about 40 times lower and with opposite sign with respect to the higher frequency band. The χ 0 changes from −82 ± 3 • at 1.3 mm to −41.9 ± 0.8 • at 3 mm. These large differences may be explained with opacity effects ( § 5.1.1; see also Hovatta et al. 2019). The EVPAs measured at 3 mm and 1.3 mm are in excellent agreement with predictions from the AMAPOLA survey (which however over-predicts LP∼3.5% at 1.3 mm; see Fig. 16).
The total intensity images both at 1.3 mm and 3 mm display, besides the bright core, a bright, one-sided jet extending approximately 20 (54 kpc) to the SW. In the higher resolution 1.3 mm image ( Fig. 12), the bright component of the jet is narrow and nearly straight, starts at a separation of ∼10 from the core and has a length of ∼10 . We also detect (at the 3σ level) two weak components of the inner jet (within ∼10 from the core) joining the bright nucleus to the outer jet. The jet structure is qualitatively similar to previous λ cm images made with the VLA at several frequencies between 1.3 and 43 GHz (e.g., Perley & Meisenheimer 2017), where the outer jet appears highly linearly polarized 6 . We do not detect LP in the jet feature.
OJ 287 -The bright blazar OJ 287 (z=0.306, 1 arcsec = 4.7 kpc) is among the best candidates for hosting a compact supermassive binary black hole (e.g. Valtonen et al. 2008). OJ287 was observed on three days at 1.3 mm and one day at 3 mm 7 . OJ287 is one of the most highly polarized targets both at 1.3 mm (LP∼ 7 − 9%) and 3 mm (LP = 8.8%). LP drops from 9% on Apr 5 down to 7% on Apr 10, while the EVPA is stable around [-59.6 • ,-61.8 • ] at 1.3 mm and -70 • at 3 mm. The LP variation and stable EVPA are consistent with the historical trends derived from the AMAPOLA survey (see Fig. 15). Its flux density is also stable. At 1.3mm, the EVPA either does not follow a λ 2 -law (Apr 5 and 11) or the formal fit is consistent with RM = 0 (Apr 10). Although we do not have a RM detection at 1.3 mm, we measure a RM = 3050 ± 620 rad/m 2 at 3 mm. A 30 years monitoring of the radio jet in OJ287 has revealed that its (sky-projected) PA varies both at cm and mm wavelengths and follows the modulations of the EVPA at optical wavelengths (Valtonen & Wiik 2012). The observed EVPA/jet-PA trend can be explained with a jet precessing model from the binary black hole which successfully predict an optical EVPA = -66.5 • in 2017 (Dey et al., submitted), consistent with actual measurements from optical polarimetric observations during 2016/17 (Valtonen et al. 2017) and close to the EVPA measured at 3 mm and 1.3 mm with ALMA.
NRAO 530 -J1733-1304 (alias NRAO 530) is a highly variable QSO (at z = 0.902; 1 arcsec = 8 kpc) that exhibits strong gamma-ray flares. It was observed on two consecutive days at 1.3 mm and one day at 3 mm. It is linearly polarized at a ∼2.4% level at 1.3 mm but only 0.9% at 3 mm. The EVPA goes from ∼ 51 • at 1.3 mm to 39 • at 3 mm, while χ 0 is stable around 51-52 • . At 3 mm, we estimate RM = −0.21 ± 0.06 × 10 5 rad/m 2 at a 3.6σ significance, which is comparable to the interband RM between 1 and 3 mm (−0.33 × 10 5 rad/m 2 ). These RM values are in agreement with those reported by Bower et al. (2018) The arcsecond-scale structure at 1.3 mm is dominated by a compact core with a second weaker component at a separation of approximately 10 from the core towards west (Fig. 12). At 3 mm, there is another feature in opposite direction (to the east), which could be a counter jet component (Fig. 13). This geometry is apparently inconsistent with the north-south elongation of the jet revealed on scales < 100 pc by recent VLBI multifrequency (22, 43 and 86 GHz) imaging (e.g., Lu et al. 2011), although the Boston University Blazar monitoring program 8 conducted with the VLBA at 7 mm has revealed significant changes in the jet position angle over the years, and possibly jet bending.
J1924-2914 is completely unresolved on arcsecond scales both at 1.3 mm and 3 mm (see Figures 11 and  13), a result consistent with images at cm-wavelengths made with the VLA (e.g., Perley 1982).
4C 01.28 -J1058+0133 (alias 4C 01.28) is a blazar at z= 0.888 (1 arcsec = 8 kpc). It was observed on three days at 1.3 mm and one day at 3 mm. The source is strongly polarized with a mean LP of 5.5% at 1.3 mm and 4.4% at 3 mm. At 1.3 mm, the LP varies by <15% while the EVPA changes from ∼ −23 • (Apr 5) to ∼ −15 • (Apr 11); the EVPA at 3 mm, measured on Apr 2, is -32 • , apparently consistent with the trend at 1.3 mm. Both the measured EVPA and LP values at 1.3 mm and 3 mm follow very closely the time evolution predicted in the AMAPOLA survey (see Fig. 15), where the LP and EVPA follow a trend parallel to the Stokes I evolution.
On Apr 11, we tentatively detect RM ∼ (0.33 ± 0.13) × 10 5 rad m −2 at the ∼ 3σ level; we however caution that on Apr 5 and 10 the EVPAs do not follow the λ 2 trend (Fig. 17), and we do not have a RM detection at 3 mm (with a 3σ upper limit of 3600 rad m −2 ).
Cen A -Centaurus A (Cen A) is the closest radio-loud AGN (at a distance of 3.8 Mpc, 1 arcsec = 18 pc). Although it is a bright mm source (with F=5.7 Jy), it is unpolarized at 1.3 mm (with a 3σ LP upper limit of 0.09%). We find a spectral index of -0.2 in the central core, consistent with a flat spectrum, as also measured between 350 and 698 GHz with (non-simultaneous) ALMA observations (Espada et al. 2017).
The total intensity images reveal a diffuse emission component around the central bright core, extending across 12 and mostly elongated north-south, and two additional compact components towards north-east (NE) separated by roughly 14 and 18 from the central core and aligned at P.A.∼50 • (see Fig. 12, bottom-right panel). The first component could be associated with the inner circumnuclear disk, mapped in CO with the SMA (Espada et al. 2009) and ALMA (Espada et al. 2017), and may indicate the presence of a dusty torus. The two additional components correspond to two knots of the northern lobe of the relativistic jet, labelled as A1 (inner) and A2 (outer) in a VLA study by Clarke et al. (1992); no portion of the southern jet is seen, consistent with previous observations (McCoy et al. 2017).
NGC 1052 -NGC 1052 is a nearby (19.7 Mpc; 1 arcsec = 95 pc) radio-galaxy that showcases an exceptionally bright twin-jet system with a large viewing angle close to 90 degrees (e.g., Baczko et al. 2016). With F∼0.4 Jy and LP<0.15%, it is the weakest mm source (along with J0132-1654) and the second least polarized AGN in our sample. The apparent discrepancy in flux-density and spectral index between Apr 6 and 7 is most likely a consequence of the low flux-density (below the threshold required by the commissioned on-source phasing mode; see § 2.2) and the much poorer data quality on Apr 7, rather than time-variability of the source.

M87
We report the first unambiguous measurement of RM toward the M87 nucleus at mm wavelengths (Table 2; Figure 17, middle panels). We measure (1.51 ± 0.30) × 10 5 rad m −2 (with a 5σ significance) on Apr 6 and tentatively (0.64 ± 0.27) × 10 5 rad m −2 (with a 2.4σ significance) on Apr 5. On the last two days we can only report best-fit values of (−0.24 ± 0.23) × 10 5 rad m −2 (with a 3σ confidence level range [−0.93, 0.45]) on Apr 10 and (−0.39 ± 0.24) × 10 5 rad m −2 (with a 3σ confidence level range [−1.11, 0.33]) on Apr 11. Although we cannot determine precisely the RM value on all days, we can conclude that the RM appears to vary substantially across days and there is marginal evidence of sign reversal.
Before this study, the only RM measurement was done with the SMA at 230 GHz by Kuo et al. (2014), who reported a best-fit RM = (−2.1 ± 1.8) × 10 5 rad m −2 (1σ uncertainty) and could therefore only provide an upper limit. In order to better constrain the RM amplitude and its time variability, in addition to the 2017 VLBI observations (which are the focus of this paper), we have also analysed the ALMA data acquired during the Apr 2018 VLBI campaign as well as additional ALMA archival polarimetric experiments (these are introduced in §2.4). For two projects (2016.1.00415.S and 2017.1.00608.S) we produced fully-calibrated uvfiles and then used the uvmf flux extraction method with uvmultifit to determine the M87 Stokes parameters. For the remaining two projects (2013.1.01022.S and 2015.1.01170.S), we used the full-Stokes images released with QA2. Since these images do not include clean component models, we used the intf method to extract the Stokes parameters in the compact core directly in the images 10 . Table 3 reports the full list of ALMA observations, project codes, and derived polarimetric parameters. In total, we have collected data from three and eight different observations at 3 mm and 1.3 mm, respectively, spanning three years (from Sep 2015 to Sep 2018). The main findings revealed by the analysis of the full dataset are the following:

The magnitude of the RM varies both at 3 mm (range
5. The RM can either be positive or negative in both bands (with a preference for a negative sign), indicating that sign flips are present both at 3 mm and 1.3 mm.
6. In Apr 2017, the RM magnitude appears to vary significantly (from non-detection up to 1.5 × 10 5 rad m −2 ) in just 4-5 days.
We will interpret these findings in Section 5.2.

Sgr A*
In this section, we analyse the polarimetric properties of Sgr A* and its variability on a week timescale based on the ALMA observations at 1.3 mm and 3 mm.
LP. -We measured LP between 6.9% and 7.5% across one week at 1.3 mm (Table 2). These values are broadly consistent with historic measurements using BIMA on several epochs in the time span 2002-2004 at 227 GHz (7.8-9.4%; Bower et al. 2003Bower et al. , 2005, SMA on several days in Jun-Jul 2005 (4.5-6.9% at 225 GHz; Marrone et al. 2007), and more recently with ALMA in Mar-Aug 2016 at 225 GHz (3.7-6.3%, 5.9% mean; Bower et al. 2018, who also report intra-day variability). Besides observations at 1.3 mm, LP variability has been reported also at 3.5 mm with BIMA (on a timescale of days- Macquart et al. 2006) and at 0.85 mm with the SMA (on a timescale from hours to days- Marrone et al. 2006). All together, these measurements imply significant time-variability of LP across timescales of hours/days to months/years. While at 1.3 mm LP ∼ 7%, at 3 mm we detect LP 1% (Table 1). It is interesting to note that the LP fraction increases from ∼ 0.5% at ∼ 86 GHz (our SPW=0,1) up to ∼ 1% at ∼ 100 GHz (our SPW=2,3; see Table 9). This trend is consistent with earlier measurements at 22 GHz and 43 GHz with the VLA, and at 86 GHz and 112-115 GHz with BIMA, yielding upper limits of LP∼ 0.2, 0.4%, 1% (Bower et al. 1999b), and 1.8% (Bower et al. 2001), respectively (but see Macquart et al. 2006, who report LP ∼ 2% at 85 GHz with BIMA observations in Mar 2004).
Variations in RM appear to be coupled with LP fraction: the lower the polarization flux density, the higher the absolute value of the RM. In particular, we find ∆LP ∼ +5% (∆RM ∼ −9%) and ∆LP ∼ +9% (∆RM ∼ −32%) in Apr 7 and 11, respectively, with respect to Apr 6. This can be understood if a larger RM scrambles more effectively the polarization vector fields resulting in lower net polarization. Although with only three data points we cannot draw a statistically significant conclusion, we note that the same trend was also seen by Bower et al. (2018) on shorter (intra-day) timescales.
We report for the first time a measurement of RM at 3 mm, with a magnitude of (−2.1 ± 0.1) × 10 5 rad m −2 (Table 1; Figure 17, upper-left panel). The RM magnitude at 3 mm (measured on Apr 3) is a factor of 2.3 (2.1) smaller than the RM value measured at 1.3 mm on Apr 6 (Apr 7). Furthermore, we note a large offset in χ 0 between the 3 mm (+135 • or -45 • for a full Figure 5. Sgr A* EVPA as a function of λ 2 observed in 2017 at 3 mm (Apr 3) and 1.3 mm (Apr 6, 7, 11). Each grey line is a linear fit to the EVPAs measured at the four ALMA Band 3 and Band 6 SPWs, yielding the RM in each epoch, and the extrapolated intercept at the Y-axis is χ0. Note the large offset in χ0 between the 3 mm and 1.3 mm bands, despite the consistency of χ0 at 1.3 mm across 6 days. 180 • wrap) and the 1.3 mm bands (∼ [−14.7, −18.8] • ), which is unlikely a consequence of time variability (given the χ 0 consistency on Apr 6-11). The comparison of RM and χ 0 in the two frequency bands (showcased in Fig. 5) indicates the presence of both Faraday and intrinsic changes of the source. We will provide an interpretation of the differences observed between the two frequency bands in § 5.3.
CP. -We report a tentative detection of CP at 1.3 mm in the range (−1.0±0.6)% to (−1.5±0.6)%. This is consistent with the first detection with the SMA from observations carried out in -2007(Muñoz et al. 2012 and with a more recent ALMA study based on 2016 observations (Bower et al. 2018). This result suggests that the handedness of the mm-wavelength CP is stable over timescales larger than 12 yr. Interestingly, historical VLA data (from 1981 to 1999) between 1.4 and 15 GHz show that the emission is circularly polarized at the 0.3% level and is consistently left-handed (Bower et al. 1999a(Bower et al. , 2002, possibly extending the stability of the CP sign to 40 years. Such a remarkable consistency of the sign of CP over (potentially four) decades suggests a stable magnetic field configuration (in the emission and conversion region).
Similarly to the RM, we also note a weak anticorrelation between LP and CP (although more observations are needed to confirm it).
We do not detect CP at 3 mm (< 0.06%, 3σ upper limit). Muñoz et al. (2012) and Bower et al. (2018) find that, once one combines the cm and mm measurements, the CP fraction as a function of frequency should be characterized by a power law with ν 0.35 . Using the measurements at 1.3 mm, this shallow power-law would imply a CP fraction at the level of ∼ 0.7−1.1% at 3 mm, which would have been readily observable. The non detection of CP at 3 mm suggests that the CP spectrum may not be monotonic.
Although the origin of the CP is not well understood, since a relativistic synchrotron plasma is expected to produce little CP, Muñoz et al. (2012) suggest that the observed CP is likely generated close to the event horizon by the Faraday conversion which transforms LP into CP via thermal electrons that are mixed with the relativistic electrons responsible for the linearly polarized synchrotron emission (Beckert & Falcke 2002). In this scenario, while the high degree of order in the magnetic field necessary to produce LP ∼ 7% at 1.3 mm naturally leads to a high CP in a synchrotron source, the absence of CP at 3 mm is consistent with the low LP measured. See Muñoz et al. (2012) for a detailed discussion of potential origins for the CP emission.
A final caveat is that based on the analysis presented in §3.3.4 and Appendix G, the physical interpretations above should be considered as tentative.
Flux-density variability -We do not report significant variability in total intensity and polarized intensity, which is about 10% in six days (comparable to the absolute flux-scale uncertainty in ALMA Band 6). Marrone et al. (2006) and Bower et al. (2018) report more significant variability in all polarization parameters based on intra-day light curves in all four Stokes parameters. This type of analysis is beyond the scope of this paper, and will be investigated elsewhere.

DISCUSSION
In this section, we review general polarization properties of AGN comparing the two (1.3 mm and 3 mm) frequency bands, different AGN classes, and depolarization mechanisms ( § 5.1); then we interpret the Faraday properties derived for M87 in the context of existing accretion and jet models as well as a new two-component polarization model ( § 5.2); and finally we discuss additional constraints on the Sgr A* polarization model from a comparison of 1.3 mm and 3 mm observations ( § 5.3).

1.3 mm vs. 3 mm
Synchrotron emission opacity -The total intensity spectral indexes for the AGN sources in the sample vary in the range α=[-0.7,-0.3] at 3 mm and α=[-1.3,-0.6] at 1.3 mm, Cen A being the only exception, with α=-0.2 (see Tables 1 and 2 and Appendix H). This contrasts with the flat spectra (α=0) typically found at longer cm wavelengths in AGN cores (e.g., Hovatta et al. 2014), corresponding to optically-thick emission. In addition, we observe a spectral steepening (with ∆α = [−0.2, −0.4]) between 3 mm and 1.3 mm; although one should keep in mind the caveat of time variability, since the observations in the two frequency bands were close in time (within ten days) but not simultaneous. Such spectral steepening can naturally be explained by decreased opacity of the synchrotron emission at higher frequencies in a standard jet model (e.g., Blandford & Königl 1979;Lobanov 1998).
LP degree -We detect LP in the range 0.9-13% at 3 mm and 2-15% at 1.3 mm (excluding the unpolarized sources NGC 1052 and Cen A). At 1.3 mm, the median fractional polarization is 5.1%, just slightly higher than the median LP at 3 mm, 4.2%, yielding a ratio of 1.2. If we consider only the sources observed in both bands, then the ratio goes slightly up to 1.3 (or 1.6 including also Sgr A*). Despite the low statistics, these trends are marginally consistent with results from previous singledish surveys with the IRAM 30-m telescope (Agudo et al. 2014  The polarimetric data analysis is based on Earth rotation polarimetry and is antenna-based, i.e. executed for each antenna of polarization at 1 mm as compared to 3 mm. This finding can be related either to a smaller size of the emitting region and/or to a higher ordering of the magnetic-field configuration (e.g., see discussion in Hughes 1991). In fact, according to the standard jet model, the size of the core region decreases as a power-law of the observing frequency, which could help explain the higher LP observed at 1 mm. Alternatively, the more ordered magnetic-field configuration could be related to a large-scale (helical) magnetic-field structure along the jets.
Faraday RM -Among the six sources observed both at 3 mm and 1.3 mm, we have RM detections at the two bands only in 3C273, where the estimated value at 3 mm is significantly lower than at 1.3 mm. For the remaining sources with RM detections at 3 mm (NRAO 530, OJ 287, and 3C279) and at 1.3 mm (J1924-2914 and 4C 01.28), their 3σ upper limits, respectively at 1.3 mm and 3 mm, still allow a larger RM at the higher frequency band.
A different 'in-band' RM in the 3 mm and 1.3 mm bands can be explained either with (i) the presence of an internal Faraday screen or multiple external screens in the beam; or with (ii) a different opacity of the synchrotron emission between the two bands. Case (i) will cause a non-λ 2 behaviour of the EVPA and a non-trivial relation between the 'in-band' RM determined at only two narrow radio bands. Evidence for non-λ 2 behaviour of the EVPA can be possibly seen in OJ 287, 4C 01.28, and J0006-0623 at 1.3 mm (Fig. 18) and J0510+1800 separately. Therefore, no interferometric polarization images are available from this study. at 3 mm (Fig. 19). In order to estimate B or n e from the RM (see Eq. 1), one would need to sample densely the EVPA over a broader frequency range and perform a more sophisticated analysis, using techniques like the Faraday RM synthesis or Faraday tomography (e.g. Brentjens & de Bruyn 2005). This type of analysis is beyond the scope of this paper and can be investigated in a future study (we refer to § 5.2.1 for evidence of internal Faraday rotation and § 5.2.2 for an example of a multiple component Faraday model for the case of M87).
Since the spectral index analysis shows that the AGN in the sample become more optically thin at 1 mm, the observed differences in the 'in-band' RM at 3 mm and 1.3 mm can be likely explained with synchrotron opacity effects alone (with the caveat of time variability since the observations are near-in-time but not simultaneous).
It is also interesting to note that we also see a sign reversal between the RM measured at 3 mm and 1.3 mm for 3C 273. RM sign reversals require reversals in B || either over time (the observations in the two bands were not simultaneous) or across the emitting region (the orientation of the magnetic field is different in the 3 mm and 1 mm regions). With the data in hand we cannot distinguish between time variability or spatial incoherence of the magnetic field (we refer to § 5.2 for a discussion on possible origins of RM sign reversals in AGN).

Blazars vs. other AGN
We find that blazars are more strongly polarized than other AGN in our sample, with a median LP ∼7.1% vs. 2.4% at 1.3 mm, respectively. Furthermore, blazars have approximately an order-of-magnitude lower RM values (on average) than other AGN, with a median value of ∼ 0.07×10 5 rad m −2 at 1.3 mm (with the highest values of ∼ 0.4 × 10 5 rad m −2 exhibited by J1924-2914), whereas for other AGN we find a median value of ∼ 0.4 × 10 5 at 1.3 mm 14 (with the highest values > 10 5 rad m −2 exhibited by M87 and 3C 273). Bower et al. (2017) used the Combined Array for Millimeter Astronomy (CARMA) and the SMA to observe at 1.3 mm two low-luminosity AGN (LLAGN), M81 and M84, finding upper limits to LP of 1%-2%. Similarly, Plambeck et al. (2014) used CARMA to observe the LLAGN 3C 84 at 1.3 and 0.9 mm, measuring an LP in the 1%-2% range, and a very high RM of ∼ (9±2)×10 5 . These low values of LP (and high values of RM) are comparable to what we find in M87, which is also classified as a LLAGN (e.g. Di Matteo et al. 2003).
When put together, these results suggest that blazars have different polarization properties at mm wavelengths from all other AGN, including LLAGN, radio galaxies, or regular QSOs 15 . These mm polarization differences can be understood in the context of the viewing angle unification scheme of AGN. A smaller viewing angle implies a stronger Doppler-boosting of the synchrotron emitting plasma in the jet, which in turn implies a higher polarization fraction for blazars. Furthermore, their face-on geometry allows the observer to reach the innermost radii of the nucleus/jet and reduces the impact of the 'scrambling' of linearly polarized radiation by averaging different polarization components within the source (e.g. Faraday and beam depolarization -see next section), also resulting in higher LP (and lower RM).

Depolarization in radio galaxies
In the previous section we point out that radio galaxies and LLAGN exhibit lower polarization degree than blazars. In particular, the radio galaxies Cen A and NGC 1052 do not show appreciable polarized intensity (LP<0.2%) at 1.3 mm. We suggest several depolarization mechanisms that may be at play in these radio galaxies (and potentially other LLAGN): a. Faraday depolarization due to a thick torus or a dense accretion flow.
b. Bandwidth depolarization due to a strong magnetic field.
c. Beam depolarization due to a disordered magnetic field.
In the following, we elaborate on these mechanisms.
a. Faraday depolarization due to a thick torus or a dense accretion flow -Radio galaxies are often surrounded by an obscuring torus. The cold gas in the torus can be photoionized by UV photons from the inner accretion disk and the mixture of thermal and non-thermal material could be responsible for the strong depolarization of the inner regions via Faraday rotation -and one speaks of Faraday depolarization. In the case of NGC 1052, Fromm et al. (2018Fromm et al. ( , 2019 created synthetic radio maps of the jets using special-relativistic hydrodynamic (SRHD) simulations and suggested that an obscuring torus can explain some of the observed properties of these jets. In fact, the presence of a massive (∼ 10 7 M ) and dense (> 10 7 cm −3 ) molecular torus has been recently demonstrated with ALMA observations (Kameno et al. 2020). A clumpy torus is also known to surround the nucleus of Cen A (e.g. Espada et al. 2017), as also suggested by our 1.3 mm map (see Fig. 12, bottom-right panel). Therefore the presence of a thick torus of cold gas could 15 Similar conclusions were reached from VLBI imaging studies of large AGN samples at cm wavelengths (e.g., Hodge et al. 2018). naturally explain the low polarization degree in both radio galaxies. A similar mechanism can be at play in LLAGN whose radio emission is thought to be powered by synchrotron radiation from a geometrically thick, hot accretion flow (e.g., Narayan & Yi 1994;Blandford & Begelman 1999;Quataert & Gruzinov 2000).
b. Bandwidth depolarization due to a strong magnetic field -A large homogeneous magnetic field implies an intrinsically large homogeneous RM (see Eq. 1), resulting in the source to appear unpolarized in broadband observations -and one speaks of bandwidth depolarization. In NGC 1052, GMVA imaging at 86 GHz helped constrain the magnetic field at Schwarzschild radius (R Sch ) scales in the range 360-70000 G (Baczko et al. 2016), providing evidence of an extremely high magnetic field near the SMBH. Coupled with its high inclination (e.g. Kadler et al. 2004), such a strong magnetic field would result in the source to appear unpolarized in ALMA broadband observations. c. Beam depolarization due to a tangled magnetic field -If the magnetic field in the emitting regions or in a foreground Faraday screen is tangled or generally disordered on scales much smaller than the observing beam, magnetic field regions with similar polarization degrees but opposite signs will cancel out and the net observed polarization degree would be significantly decreased -and one speaks of beam depolarization. We do not have evidence of such tangled magnetic field for any of the low LP sources in our sample, which will require high-resolution polarization imaging with the GMVA or the EHT. d. Thermal (non-synchrotron) emission -Multi-wavelength (MWL) studies in Cen A show that its SED is moderately inverted up to the infrared, possibly indicating a dust contribution at mm wavelengths (e.g. Espada et al. 2017). Using VLBI imaging at 229 GHz with the EHT, Janssen et al. (submitted) measured a flux of ∼ 2 Jy in the VLBI core, indicating that the EHT filters out ∼ 65% of the emission seen by ALMA. While we cannot exclude contribution from thermal emission to the total flux measured by ALMA, the flux measured with the EHT must necessarily be associated with non-thermal emission. We therefore conclude that dust emission is an unlikely explanation for the lack of LP at 1.3 mm.
An improved data analysis including spectro-polarimetry could be helpful to measure the actual RM in both NGC 1052 and Cen A and thus assess which is the dominant depolarization mechanism among the ones discussed above (this is however beyond the scope of this paper). As a final note, an interesting insight may come from a comparison between mm and IR wavelengths, where both NGC 1052 (Fernández-Ontiveros et al. 2019) and Cen A (Jones 2000; Lopez Rodriguez 2021) are highly polarized. These characteristic are similar to Cygnus A, where the low polarized core at mm-wavelengths and the high polarized core at IR wavelengths may be the signature of an ordered magnetic field in the plane of the accretion disk supporting the accretion flow and/or jet formation (Lopez-Rodriguez et al. 2018).

Physical origin of the rotation measure in M87
We can now use the polarimetric and Faraday properties of the mm emission from M87 reported in § 4.2, to constrain properties of accretion models onto the M87 SMBH. Models aiming to explain the origin of the RM in M87 should match the following key observed features (see § 4.2 for a full list of findings): i. Low LP and high RM. M87 has a rather low LP (∼ 2.3% mean at 1.3 mm) when compared to Sgr A* and other blazars in the sample (see §5.1.2), while the measured RM can be as high as a few times 10 5 rad m −2 .
ii. RM sign reversals. Observations on different dates yield large differences in the measured RM values, which can be either positive or negative (in both the 3 mm and 1.3 mm bands). This requires sign flips in B || over time and/or across the emitting region.
iii. Rapid RM time variability. In Apr 2017, the RM magnitude appears to vary significantly (from nondetection up to 1.5 × 10 5 rad m −2 ) in just 4-5 days. This suggests the presence of small-scale fluctuations in the emitting source and/or the Faraday screen.
iv. λ 2 scaling. Plots in Figures 17 and 20 clearly display a λ 2 dependence of the EVPA at 1.3 mm and 3 mm on specific days (although this is not always the case).
The MWL spectral energy distribution (SED) of the M87 core is best explained by emission from advection dominated/radiatively inefficient accretion flow (ADAF/RIAF - Reynolds et al. 1996;Di Matteo et al. 2003) or from a jet (e.g., Dexter et al. 2012;Prieto et al. 2016;Mościbrodzka 2019), or emission from a combination of these two components (e.g., Broderick & Loeb 2009;Nemmen et al. 2014;Mościbrodzka et al. 2016;Feng et al. 2016;Mościbrodzka et al. 2017;Davelaar et al. 2019;Event Horizon Telescope Collaboration et al. 2019e). In the hybrid models, the low-frequency radio emission is produced by the jet while the opticallythin mm/submm (and X-ray) emission can either come from the jet base or the inner accretion flow.
The traditional approach adopted in previous studies was to assume that the large scale (r ∼ 100 R Sch ) accretion flow itself may act as a Faraday screen and that the core emission region lies entirely behind the same portion of the Faraday screen (e.g., the core emission size is small compared to the scale of any fluctuations in large scale flow). In the framework of semi-analytic RIAF/ADAF models, the RM magnitude has then been used to estimate mass accretion rates onto black holes in Sgr A* (Marrone et al. 2006(Marrone et al. , 2007 and in 3C 84 (Plambeck et al. 2014). A similar approach has been used in M87, yielding estimates of the accretion rate in the range fromṀ < 9 × 10 −4 (Kuo et al. 2014) tȯ M ∼ [0.2, 1] × 10 −3 M yr −1 (Feng et al. 2016), where the quoted values are either upper limits or depend on specific model assumptions (e.g. black hole spin or the exact location of the Faraday screen). From the largest RM that we measured, 4 × 10 5 rad m −2 (on Apr 2018), using Equation (9) in Marrone et al. (2006), we would infer a mass accretion rate ofṀ = 7.7 × 10 −8 |RM | 2/3 ∼ 4 × 10 −4 M yr −1 assuming an inner boundary to the Faraday screen of R RM,in = 21 R Sch (as suggested by Kuo et al. 2014).
While our estimates of mass accretion rates are consistent with previous similar estimates and upper limits from the ADAF/RIAF models, the observed properties listed above, especially the time variable RM and its sign reversals, provide new constraints. In particular, the timescale of the RM variability can be set by the rotating medium dynamical time (∝ R 3 /GM ), thus constraining the radius at which the RM originates. The rapid variability observed in Apr 2017 implies that the RM should occur much closer to the SMBH (within a few R Sch ) than assumed in previous mass accretion models, which in turn suggests the possibility of a co-location of the emitting and rotating medium. In the alternative, the Faraday screen could be at further distance and the observed variability could be ascribed to rapid fluctuations in the emitting source. Therefore, both a turbulent accretion flow acting as a Faraday screen or a varying compact source with an external screen can explain the observed time-variability. Finally, the accretion flow is not the only possible source of Faraday rotation. Since simulations show that relativistic jets can have a "spinesheath" structure (e.g., McKinney 2006), the jet sheath can provide a magnetized screen surrounding the jet, and indeed it has been also suggested as a plausible source of Faraday rotation in AGN (e.g., Zavala & Taylor 2004). Therefore, either (or both) the accretion flow and/or the jet can in principle be the sources of the mm emission and/or the Faraday rotation.
All the scenarios described above imply a more complicated physical origin of the Faraday rotation than is usually assumed in traditional semi-analytic models that use the RM to infer a mass accretion rate. We conclude that, unlike the case of Sgr A*, the RM in M87 may not provide an accurate estimate of the mass accretion rate onto the black hole.
In what follows, we review clues on the location of the Faraday screen using observational constrains from ALMA ( § 5.2.1) as well as information on horizon scales from the EHT ( § 5.2.2).

Location of the Faraday screen: internal vs. external
We distinguish between two general cases: internal and external Faraday rotation.
I. Internal Faraday rotation: the accretion flow or jet can simultaneously be the source of synchrotron radiation and the Faraday screen.
a. RM rapid time variability and sign reversals.
Recent time-dependent GRMHD simulations of the M87 core (Ricarte et al. 2020) show that turbulence within the accretion flow is able to change B || in both amplitude and orientation, resulting in significant RM fluctuations and sign reversals on the dynamical time at R 2.5 − 5 R Sch , corresponding to short timescales of a few days for M87 (consistent with properties #ii and #iii).
An internal Faraday screen could cause beam depolarization of the synchrotron emission. This has been theoretically predicted by GRMHD simulations of the M87 core emission, which yield low values of LP (typically in the range 1%-3%) and large Faraday RM ( 10 5 rad m −2 ) (Mościbrodzka et al. 2017;Event Horizon Telescope Collaboration et al. 2021b), broadly consistent with the observed feature #i. We however note that the beam depolarization could be also caused by external Faraday rotation in an inhomogeneous screen.
II. External Faraday rotation: the emission region lies entirely behind (and it is not inter-mixed with) the Faraday screen.
a. λ 2 scaling of the EVPA. A λ 2 dependence is typically adopted as observational evidence of the fundamental assumption on the location of the Faraday screen as external relative to the background source. Although it can be argued that EVPA variations in a narrow frequency range could be approximated to be linear (at 1.3/3 mm the fractional λ 2 bandwidth is only 16/32% of the central wavelength), good linear fits of EVPA vs. λ 2 are also obtained from lowerfrequency observations, including the range 2, 5 and 8 GHz ( Asada et al. 2002;Gómez et al. 2008;Gabuzda et al. 2014). Systematic changes in the signs of these gradients, leading to RM sign reversals in unresolved measurements, can be explained with a number of models, including the magnetic "tower" model (Lynden-Bell 1996;Contopoulos & Kazanas 1998;Lico et al. 2017), or the "striped" jet model (Parfrey et al. 2015;Mahlmann et al. 2020;Nathanail et al. 2020).
Nevertheless, it remains difficult to explain the rapid fluctuations observed in Apr 2017 with these models.
A long-term monitoring with beam-matched simultaneous observations at multiple frequency bands would be required to assess the frequency-dependence of the RM and to conclusively discriminate between internal and external Faraday rotation. Clear evidence of λ 2 scaling in a wider frequency range would be evidence of the external scenario, while deviations from λ 2 would support the internal scenario. A time cadence from a few days to a few months would allow us to assess whether the RM sign flips are stochastic (favoring the internal scenario), or systematic (favoring the external scenario).

Two-component polarization model for M87
EHT estimates of the flux in the compact core, i.e., that arising on horizon scales, can at most account for 50% of that measured by ALMA, both in total intensity (Event Horizon Telescope Collaboration et al. 2019d) and polarized intensity (Event Horizon Telescope Collaboration et al. 2021a) emission. Natural origins for the additional emission are the larger-scale structures in the jet. The additional components may encompass many scales, be discrete features (e.g., HST-1), or be a combination thereof. In order to interpret these differences revealed by the EHT and ALMA, we adopt the simplest version of a multi-scale model permissible -a two-component model comprised of a variable compact region and static extended region (see Figure 6). We find that this is sufficient to harmonize the polarimetric properties observed by both the EHT and ALMA in Apr 2017, including the interday variability in the ALMA RMs and the EVPA variation of the compact core as observed by the EHT.
Both the compact and extended components of the two-component model consist of total intensity, spectral index, linearly polarized flux, and polarization an- gle. We consider both internal and external Faraday screen models for the compact component. In both cases, the Faraday screen for the extended component is assumed to be external. A model likelihood is constructed using the integrated EHT Stokes I, Q, and U ranges presented in Table 7 in Appendix H2 of Event Horizon Telescope Collaboration et al. (2021a), and the ALMA core Stokes I, Q, and U values for the individual SPWs in Table 10, assuming Gaussian errors. This likelihood is then sampled with the EMCEE python package (Foreman-Mackey et al. 2013) to obtain posterior probability distributions for the model components. For more details regarding the model, priors, and fit results, see Appendix I.
Across days, only the LP and EVPA of the compact component is permitted to vary. This is consistent with the extended component being associated with much larger physical structures and required by the polarimetric variability observed by the EHT (Event Horizon Telescope Collaboration et al. 2021a). There is no evidence that variability in any other component of the two-component model is required: despite static Faraday screens, permitting the polarization of the compact component to vary is capable of reproducing the rapid changes in the ALMA RMs. In this picture, the observed RMs arise in part from the wavelength-dependent competition between the two components, and thus are not directly indicative of the properties of either Faraday screen. Figure 7. Posteriors implied by the April 5, 6, 10, and 11, 2017 ALMA and EHT observations for the RMs measured on ALMA (extended) and EHT (compact) scales in the two component model when it is assumed that the compact Faraday screen is external (lower left, blue) and internal (upper right, red) to the emission region. Contours are shown for the 50%, 90%, and 99% quantiles. For comparison, the 2σ range of April 5, 6, 10, and 11, 2017 ALMA RMs reported in Table 3 are shown by the black crosses and gray bands.
Nevertheless, via this model we are able to separately constrain the RMs that are observed on ALMA and EHT scales; these are shown in Figure 7. Specifically, while the RM of the large-scale Faraday screen is comparable to the Apr 2017 values reported in Table 3, that associated with the compact component is not directly constrained by the ALMA measurements and can be factors of many larger: at 95% confidence, the compact RM is between −5.4 × 10 5 rad m −2 and −2.9 × 10 5 rad m −2 . Interestingly, the estimated range is consistent with the RM inferred from low-inclination GRMHD models of the M87 core (Ricarte et al. 2020;Event Horizon Telescope Collaboration et al. 2021b) and comparable to the Apr 2018 values reported in Table 3. This consistency suggests the possibility that in Apr 2018 ALMA may be seeing the core RM (e.g., as a consequence of a different beating of the two components). This hypothesis can be directly tested with the 2018 EHT observations which, unlike the 2017 ones, covered the full frequency spacing of ALMA (212-230 GHZ; see § 2.1), and are therefore expected to directly measure the resolved RM of the core. This will in turn allow to quantify the interplay between compact and extended components, and potentially explain the time variability observed with ALMA.

Constraints on Sgr A* model from polarization and Faraday properties at 3 mm and 1.3 mm
Measurements of Faraday rotation at radio/mm wavelengths, either towards Sgr A* itself (e.g., Marrone et al. 2007;Bower et al. 2018) or the nearby pulsar PSR J1745-2900 (e.g., Eatough et al. 2013;Kravchenko et al. 2016), have been used to probe the magnetized accretion flow in Sgr A* on scales from tens of R Sch out to the Bondi radius (∼ 10 5 R Sch ). Using the same semi-analytic RIAF/ADAF models introduced in § 5.2 (Marrone et al. 2006), from the measured RM values at 1.3 mm we obtain an accretion rate of order 10 −8 M yr −1 (assuming R RM,in = 10 R Sch ), with a maximum variation of approximately 20% across the observing week.
We have derived for the first time the polarization and Faraday properties of Sgr A* both at 3 mm and 1.3 mm in a time-window of ten days. Since the synchrotron photosphere in the accretion flow moves outwards with decreasing frequency (because of increased opacity), the polarized emission at 3 mm and 1.3 mm is expected to arise from different locations with potentially different magnetic field structures. Any variation of the intrinsic EVPA or RM with frequency can therefore provide interesting insights on the polarized source and magnetic field structure. In § 4.3 we infer that the intrinsic polarization vector is rotated between the 3 mm (χ o ∼ +135 • or -45 • assuming a full 180 • wrap) and the 1.3 mm (χ o ∼ [-15 • ,-19 • ]) bands and that the RM magnitude at 3 mm is about half of the RM value measured at 1.3 mm over a three-days separation. From earlier VLBI measurements, we know that the emission at mm wavelengths must come from very closely situated regions of the black hole, with an intrinsic (i.e. unscattered) size of ∼120 µas (or 12 R Sch ) at 3.5 mm (Issaoun et al. 2019) and ∼50 µas (or 5 R Sch ) at 1.3 mm . Therefore the radius of the 3.5 mm source is 2.4 times larger than the 1.3 mm source, i.e. very close to the ratio of the RM values measured with ALMA at the two wavelengths (see RM paragraph in § 4.3). Taken at face value, this result suggests that about half of the Faraday rotation at 1.3 mm may occur between the 3 mm photosphere and the 1.3 mm source.
Although this result would be extremely constraining for the model of Sgr A*, we should point out two caveats: i) possible presence of multiple components; ii) potential RM time variability. We explain these caveats below.
available VLBI polarimetry study at 1.3 mm has shown that the polarization structure of Sgr A* is complex (Johnson et al. 2015) and can be in principle different at the two wavelengths. Therefore, we cannot completely exclude the presence of an additional jet component to the accretion flow or a more complex morphology for the intrinsic polarization. Nevertheless, we argue that the stability of LP, RM, and CP (including their sign) observed in Sgr A* over more than a decade, unlike the case of M87, favors the presence of one single dominating polarized component.
RM time variability -At 1.3 mm, the RM changes by +1.5 × 10 5 rad/m 2 (∼ 30%) between the Apr 6/7 and Apr 11. Assuming the same rate of time variability at 3 mm and 1.3 mm, such variability is likely not responsible for a factor of two difference over three days. Likewise, the large offset in χ o observed at 3 mm and 1.3 mm is unlikely a consequence of time variability, given the χ 0 consistency on Apr 6-11. Larger variations in RM were however recorded by Marrone et al. (2007) and Bower et al. (2018) on timescales from hours to months (see for example Figure 12 in Bower et al. 2018). Since the observations in the two frequency bands were close in time but not simultaneous, we cannot definitely exclude time variability as the origin of the observed difference in RM magnitude at 3 mm and 1.3 mm.
Future simultaneous measurements over a wider wavelength range (including 3 mm and 1.3 mm) will allow us to separate time variability and source structure effects.

CONCLUSIONS
We have determined and analysed the polarization and Faraday properties of Sgr A*, the nucleus of M87, and a dozen radio-loud AGN, observed with ALMA during the 2017 VLBI campaign in the 3 mm and 1 mm bands in concert with the GMVA and the EHT, respectively.
Our main findings can be summarised as follows: 1. The AGN sources in our sample are highly polarized, with linear polarization degrees in the range 2-15% at 1.3 mm and 0.9-13% at 3 mm. The radio galaxies NGC 1052 and Cen A are the only exceptions with LP<0.2%.
2. The AGN sources have negative spectral indexes varying in the range α = [−1.3, −0.2], in contrast with the flat spectra (α=0) typically found at longer cm wavelengths in AGN cores. We also observe a spectral steepening between the 3 mm and the 1.3 mm bands, which can naturally be explained by decreased opacity of the synchrotron emission at higher frequencies in a standard jet model (e.g., Lobanov 1998).
3. We find marginal evidence for a general higher degree of polarization and RM magnitude in the 1 mm band as compared to the 3 mm band (a trend which is consistent with single-dish surveys). The increase of polarized intensity at higher frequency may be the result of an increased magnetic-field order in the inner portions of jets and/or to the smaller size of the high-frequency emitting regions. The increase of RM with frequency can be explained by opacity effects: emission at higher frequencies is generated in, and propagates along, regions with higher magnetic fields and plasma densities (e.g. Hovatta et al. 2014). Given the smallnumber statistics (eight AGN observed at 3 mm, eleven at 1.3 mm, and six in both bands) and the caveat of time variability (in a time window of ten days), simultaneous observations of a larger AGN sample at multiple frequency bands are needed to confirm these results.
4. The blazars (seven in our sample) have on average the highest level of polarization (LP ∼7.1% at 1.3 mm) and an order of magnitude lower RM (∼ 0.07 × 10 5 rad m −2 at 1.3 mm) when compared with other AGN in our sample (with LP ∼2.4% and RM ∼ 0.4 × 10 5 rad m −2 , respectively). These mm polarization differences can be understood in the context of the viewing angle unification scheme of AGN: blazars' face-on geometry implies a stronger Doppler-boosting of the synchrotron emitting plasma in the jet and reduces the effect of Faraday and beam depolarization in the accretion flow, resulting in higher LP (and lower RM). Future observations of a broader sample of sources are necessary for assessing the statistical significance of these trends. 5. We constrain the circular polarization fraction in the observed AGN to <0.3%. For Sgr A* we report CP = [-1.0,-1.5]%, consistent with previous SMA and ALMA studies. However, we explicitly note that the ALMA observatory does not guarantee a CP level <0.6% (1σ), therefore these measurements should be regarded as tentative detections.
6. We derive for the first time the polarization and Faraday properties of Sgr A* both at 3 mm and 1.3 mm in a time-window of ten days. The RM magnitude at 3 mm, (−2.1 ± 0.1) × 10 5 rad/m 2 , is about half of the RM value measured at 1.3 mm over a three-days separation, suggesting that about half of the Faraday rotation at 1.3 mm may occur between the 3 mm photosphere and the 1.3 mm source (although we cannot exclude effects related to time variability).
7. We report the first unambiguous measurement of Faraday rotation toward the M87 nucleus at mm wavelengths. At variance with Sgr A*, the M87 RM exhibits significant changes in magnitude and sign reversals. At 1.3 mm, it spans from positive values (+1.5 × 10 5 rad m −2 at a 5σ level), to < 3σ non-detections in Apr 2017, to negative values (−3 to −4 × 10 5 rad m −2 at a 10σ level) in Apr 2018. At 3 mm the RM measured values span the range from −1.2 to 0.3 × 10 5 rad m −2 from Sep 2015 to Oct 2016. The large scatter and time variability revealed by the ALMA measurements suggest a more complicated physical origin of the Faraday rotation than is usually assumed in models using the RM to infer a mass accretion rate. We conclude that, unlike the case of Sgr A*, the RM in M87 may not provide an accurate estimate of the mass accretion rate onto the black hole.
8. The observed RM in M87 may result from Faraday rotation internal to the emission region, as commonly found in GRMHD models of turbulent accretion flows or expected in a structured jet, or from a time-varying helical magnetic field threading the jet boundary layer acting as an external Faraday screen. As an alternative, we put forward a two-component model comprised of a variable compact region and static extended region. We find that this simple model is able to simultaneously explain the polarimetric properties observed in Apr 2017 by both the EHT (on horizon scales) and ALMA (which observes the combined emission from both components).
The ALMA measurements presented in this work provide critical constraints for the calibration and analysis of simultaneously obtained VLBI data. This is an essential resource for two instruments like the EHT and the GMVA which have the resolving power to reveal polarization structures and measure magnetic field strengths and particle densities on horizon scales (in the case of M87 and Sgr A*) and/or in the inner few parsecs for the AGN.
ALMA is a partnership of the European Southern Observatory (ESO; Europe, representing its member states), NSF, and National Institutes of Natural Sciences of Japan, together with National Research Council This work used the Extreme Science and Engineering Discovery Environment (XSEDE), supported by NSF grant ACI-1548562, and CyVerse, supported by NSF grants DBI-0735191, DBI-1265383, and DBI-1743442. XSEDE Stampede2 resource at TACC was allocated through TG-AST170024 and TG-AST080026N. XSEDE JetStream resource at PTI and TACC was allocated through AST170028. The simulations were performed in part on the SuperMUC cluster at the LRZ in Garching, on the LOEWE cluster in CSC in Frankfurt, and on the HazelHen cluster at the HLRS in Stuttgart. This research was enabled in part by support provided by Compute Ontario (http://computeontario.ca), Calcul Quebec (http://www.calculquebec.ca) and Compute Canada (http://www.computecanada.ca). We thank the staff at the participating observatories in the GMVA and EHT, correlation centers, and institutions for their enthusiastic support. The GMVA is coordinated by the VLBI group at the Max-Planck-Institut für Radioastronomie (MPIfR) and consists of telescopes operated by MPIfR, IRAM, Onsala, Metsahovi, Yebes, the Korean VLBI Network, the Green Bank Observatory and the VLBA. APEX is a collaboration between the MPIfR (Germany), ESO, and the Onsala Space Observatory (Sweden). The SMA is a joint project between the SAO and ASIAA and is funded by the Smithsonian Institution and the Academia Sinica. The JCMT is operated by the East Asian Observatory on behalf of the NAOJ, ASIAA, and KASI, as well as the Ministry of Finance of China, Chinese Academy of Sciences, and the National Key R&D Program (No. 2017YFA0402700) of China. Additional funding support for the JCMT is provided by the Science and Technologies Facility Council (UK) and participating universities in the UK and Canada. The LMT is a project operated by the Instituto Nacional de Astrófisica,Óptica, y Electrónica (Mexico) and the University of Massachusetts at Amherst (USA  Tables 4 and 5 list the projects observed in the 3 mm and 1.3 mm bands, ordered by date of execution. Each row reports the observing date, the ALMA project code, the science target, the source used as polarization calibrator, other sources observed in the project, and the duration of each observation. In Table 5, each row-group refers to an individual VLBI run or "Track", which includes observations of different projects carried out during the same night. The calibration of EHT projects was done per track (and not per project; see Goddi et al. 2019a). Two sources listed in Table 5, 3C 84 and J0006-0623, both observed on Apr 07, were excluded from the analysis presented in this paper: 3C 84 was observed with an elevation below 25 • , while J0006-0623 was observed for just ∼ 2 min close to an elevation of 25 • ; the resulting calibrated data display critical phase and amplitude scatter and hence were flagged before data analysis (see § 2.2).
B. POLARIMETRIC IMAGES Tables 6 and 7 report the main imaging parameters for each source observed on each day of the 2017 VLBI campaign in Band 3 and Band 6, respectively. These parameters include the on-source time, the RMS achieved in each Stokes parameter, and the synthesized beamsize. The on-source time is computed after full calibration and flagging of bad data (see §2.2). The RMS does not simply scale as T on−source but depends on several parameters such as source structure, number of observing antennas, array configuration, weather, and details of the VLBI scheduling blocks (e.g. low-elevation scans). The synthesized beamsize changes by a factor of two at 1.3 mm due to the changing array configuration during the observing week (see §2.1). The resulting images are dynamic range limited and showcase different structures on different days (depending on the array beamsize).
The full suite of polarization images for all the sources observed in the VLBI campaign are shown in  Figures 11 and 12 report 1.3 mm maps for all the other AGN sources observed with the EHT for three days and one/two days, respectively. Finally, Figure 13 shows 3 mm maps of the sources observed with the GMVA. In each plot, the black vectors showcase the orientation of the EVPAs, while their length is linearly proportional to the polarized flux. The EVPAs are plotted every 8 pixels (i.e. are spaced by 1. 6 at 1.3 mm and 4 at 3 mm) for all sources, except for M87, where the EVPAs are plotted every 4 pixels (i.e. are spaced by 0. 8), in order to sample more uniformly the jet. Note that the EVPAs are not Faraday-corrected and that the magnetic field vectors should be rotated by 90 • .

C. COMPARATIVE ANALYSIS ACROSS MULTIPLE FLUX-EXTRACTION METHODS
Some of our targets (chiefly Sgr A* and M87, but also a few other AGNs; see Appendix B) reveal extended emission at arc-seconds scales, which is differently resolved out by the different observing array configurations. Here we assess the impact of such extended emis-sion in the flux values extracted in the visibility domain for the compact cores. For this purpose, we compare the parameters derived with uvmultifit with those derived with two imaged-based methods: the sum of the nine central pixels in the clean model image (3x3) and the integrated flux (intf) from Gaussian fitting with the casa task IMFIT 16 (see §3.2). Table 8 reports the results from this comparative analysis between the image-based and the uv-based methods, showing the ratio of the four   Table 7). Note that there are several tiny EVPAs plotted across the mini-spiral, apparently locating regions with polarized flux above the image RMS noise cutoff (5σ). The LP and EVPA errors are however dominated by the systematic leakage (0.03% of I onto QU), which is not added to the images. Once these systematic errors are added, the LP flux in those points falls below the 3σ detection threshold. Therefore we do not claim detection of polarized emission outside of the central core in Sgr A*.
Stokes parameters, the LP, and the RM (3x3/uvmf and intf/uvmf, respectively), and the difference of the EVPA in degrees (3x3-uvmf and intf-uvmf, respectively). Three sources (Cen A, NGC 1052, and J0132-1654) are excluded from this statistics due to their weak polarized signal.
Overall, the analysis clearly shows that Flux(3x3) < Flux(uvmf) < Flux(intf) when considering the median of the results. On the one hand, one can interpret this as the 3×3 method being the least affected and intf being the most affected by extended emission (with a weak dependency on the array-configuration resolution). On the other hand, the 3×3 method can also be more affected by phase and amplitude calibration errors (which will remove flux from the phase-center into the sidelobes), resulting in less recovered flux. Despite this systematic deviation, the median offsets are compatible with no significant difference within the MAD values. In Stokes I, both the offset and MAD are ≤ 0.4% across the methods, which is negligible when compared to the absolute uncertainty of ALMA's flux calibration (10% in Band 6). In the case of Stokes Q and U , and resulting LP, the systematic offset and MAD between methods can be up to 1%. With the LP measured in our sample, these differences are generally comparable (in absolute terms) to the 0.03% of Stokes I leakage onto Q and U . The EVPA shows a MAD value of ∼ 0.1 deg, which results in a MAD value of up to 10% in the RM (note these are comparable or better than the observed statistical uncertainties -see Tables 9 and 10). The poorer results for Stokes V are induced by the low level of CP in the sample (although the MAD is still comfortably around 4%). Table 8 reports a fraction of "outliers" in the distributions (those cases that are 5σ away from the sample median). These are generally associated with the sources with the most prominent extended emission (other outliers are associated to parameters estimated at low-significance values). In order to better assess the magnitude of these outliers, in Figure 14 we show a detailed comparison between the uv-based and the 3x3 image-based methods for M87 (upper panels) and Sgr A* (lower panels), respectively. Analogously to Table 8, individual panels show the ratio of the four Stokes parameters, the LP, and the RM (3x3/uvmf), and the difference of the EVPA in degrees (3x3-uvmf); the results are reported per day and SPW.  Table 7).
In the case of M87, the two methods exhibit an excellent agreement, except for Stokes U . This is due to a combination of faint U emission at the core and appreciable emission from knots A and B. The discrepancy in Stokes U results generally in ∆EVPA 0.4 • comparable with the uncertainties quoted in Table 10, with one single case of ∆EVPA ∼ 0.6 • (on Apr 11 SPW=1). Given the consistency in RM (within 1σ), this worstcase scenario discrepancies in EVPA do not affect the results of the analysis of this paper.
In the case of Sgr A*, the discrepancies are much more pronounced in most parameters. Despite the prominence of the mini-spiral (see §3.1 and Fig. 8), Stokes I values actually agree within less than 1% between the two methods. It is however noteworthy that Flux(3x3)>Flux(uvmf) (i.e., the opposite trend with respect to Table 8). This inverted trend is due to a flux-decrement in the visibility amplitudes at around 25-30 kλ affecting most prominently uvmfit on Apr 11 (this can be assessed by inspecting amplitude vs. baseline-length uv-plots -not shown here -for the parallel hands on all days). Unlike Stokes I, Stokes QU are heavily affected, up to 20-30% in the worst cases, resulting in LP differences up to 10% and EVPA offsets up to 2 • . These large deviations in Stokes QU are systematic since all four SPW appear to deviate in each day by an approximate amount. This in turn results in consistent RM values within 1σ. We note that we do not observe the same systematic offset either in M87 or J1924-2914, which were also calibrated using 3C279 as the polarization-calibrator (on Apr 6 and 11), hence we discard calibration issues as being responsible for these systematics. One possible explanation is differential Stokes I leakage onto Q and U from extended unpolarized emission (see the spurious EVPA vectors matching the spiral-arms in Fig. 8). Note also that such unpolarized emission extends farther than the inner third of the image field of view, where beam-squash (Hull et al. 2020) combined with leakage may induce significant deviations. Given the results from this comparative analysis, we conclude that, for the purpose of the polarimetric analysis conducted in this paper, the uv-fitting method provides sufficiently accurate flux values of Stokes IQU in all cases, although for Sgr A* we observe deviations in Stokes QU when comparing the two flux extraction methods.

D. COMPARISON OF STOKES PARAMETERS WITH THE AMAPOLA POLARIMETRIC GRID SURVEY
For the purposes of absolute flux calibration, ALMA monitors the flux of bright sources (mainly blazars or QSOs) spread over the full range in right ascension (the Grid) by observing them together with solar system objects, the so-called Grid-Survey (GS), with a period of approximately 10 days. These observations are executed with the ACA, in Bands 3, 6, and 7. Since fullpolarization mode is adopted, it is possible to retrieve polarimetry information from the GS sources. This is done with AMAPOLA 17 , a set of CASA-friendly Python scripts used to reduce the full-Stokes polarimetry of GS observations with ACA. Some of our targets are part of the GS, including: 3C273 (J1229+0203); 3C279 (J1256-0547); 4C 01.28 (J1058+0133); 4C 09.57 (J1751+0939); J0510+1800; OJ287 (J0854+2006); J0006-0623; J1733-1304; and J1924-2914. Although the reported values by AMAPOLA are used for observation planning and the two arrays cover different uv-ranges, it is still useful to make a comparison with this database to assess any systematics or clear inconsistent variability within a week's time-frame. Figures 15 and 16 show the observed polarimetry parameters for the 2017 VLBI sessions (data points and errorbars), specifically Q, U , and V Stokes parameters, LP, and EVPA. The shaded ±1σ regions highlight the time variance of the same parameters as measured by AMAPOLA (an inflection in the trend means a time of GS observation). The colour-coding is such that blue refers to Band 3 measurements, green to Band 6, and red to Band 7 ones. These figures show that most Band 3 ALMA-VLBI measurements fall within the blue regions, while most Band 6 measurements fall inbetween the blue and red regions, indicating that our measurements are broadly consistent with the AMAP-OLA trends. Some potential conflicts can still be consistent with the inter-GS-cadence variability or differential time-variability between frequency bands (as also observed in some cases in the AMAPOLA monitoring). We conclude that, despite the different array specifications and data-reduction schemes, the ALMA-VLBI and AMAPOLA results are consistent between each other.  Figure 2 for a description of the plotted quantities). Note that for NGC 1052 and Cen A no EVPA could be reliably derived owing to their low level of LP. We also note that although we detect polarized flux in Cen A above the image RMS noise cutoff (5σ), once the systematic leakage (0.03% of I onto QU) is added to thermal noise, the LP flux would fall below the 3σ detection threshold. Therefore we do not claim detection of polarized emission in Cen A.

E. STOKES PARAMETERS PER ALMA FREQUENCY BAND (SPW)
We report the polarimetric quantities (Stokes IQU , LP, EVPA) per SPW in Tables 9 (GMVA sources) and 10 (EHT sources). The Stokes parameters were fitted directly on the visibilities using uvmultifit (see §3.2). Uncertainties are assessed with MC simulations, as the standard deviation of 1000 MC simulations for each Stokes parameter. The quoted uncertainties include in quadrature the fitting error and the Stokes I leakage onto Stokes QU (0.03% of I), as recommended by the ALMA observatory. The LP uncertainty is dominated by such systematic error, except for the weakest sources J0132-1654 and NGC 1052, for which the thermal noise starts to dominate. In order to correct for the LP bias in the low SNR regime, we estimate a debiased LP as LP 2 − σ LP 2 , where LP = Q 2 + U 2 /I. σ LP can be estimated using propagation of errors: where σ I , σ Q , σ U are the uncertainties in I, Q and U , respectively. Assuming σ QU ∼ σ QI ∼ σ U I ∼ 0, then The LP values quoted in Tables 1, 2, 9, and 10 have this debiased correction applied. The latter does not affect most of the sources studied here, but it is especially important for the low-polarization sources such as NGC 1052 and Cen A.
The absolute flux-scale calibration systematic error (5% and 10% of Stokes IQU fluxes in Band 3 and Band 6, respectively), is not included in Tables 9 and  10. The analysis on Stokes V is performed separately (see Appendix G and Tables 11, 12). Table 9. Polarization parameters of GMVA sources per frequency band (spw) and per day.  Table 8. The labels in the Y-axis indicate the observing day in April 2017 (5, 6, 7, 10, 11) and the observing SPW (s1, s2, s3, s4); the RM panel shows only the days. Non-detections in the images are indicated with a cross. Note that the plotted uncertainties are not those of individual measurements but ratios or differences between the two methods. Therefore the errors on LP and EVPA for M87 on Apr 10/11 appear comparable to other days despite the large error bars and/or non detections in Stokes U (i.e. the errors in U displayed do not propagate in the LP and EVPA plots).  a This table compares uvmultifit (uvmf) with: the sum of the nine central pixels (3x3) of the model image; and the integrated flux (intf) from imfit. b The values reported in columns are ratios between the method referred in the first column and uvmf, with the latter being in the denominator. c The EVPA column shows the angle difference (in deg) between the methods.
d In each instance, three values are displayed: the median value of the distribution; the Median Absolute Deviation or MAD (the value in brackets); the percentage of cases farther from the median than five times the MAD in quadrature with the measurement error (the value in the second row).
[     Figure 16. Same as Fig. 15, but for different sources (from left to right: 3C273, J0006-0623, NRAO 530, 4C 09.57). We note that source J0006-0623 was observed on Apr 7 but the calibrated data were of poor quality and were flagged before analysis (see § 2.2). Given its flat spectral index, Fspw=2 = Fspw=3 should be assumed for Sgr A*.
F. FARADAY RM PLOTS We display the EVPA data as a function of λ 2 and their RM fitted models in Figures 17, 18, 19, and 20; EVPAs are measured in each of the four ALMA SPWs. Let us first focus on Figure 17, showing the Faraday RM of Sgr A* (upper panels), M87 (middle panels), and 3C 279 (lower panels). The EVPAs measured in Sgr A* reveal a remarkably precise λ 2 dependence both at λ3 mm (left panel) and λ1.3 mm (second to fourth panels). It is also remarkable that Sgr A* showcases a very consistent slope across all days, while in M87 the slope appears to change sign from Apr 5/6 to Apr 10/11. 3C 279 was used as polarization calibrator on Apr 5, 6, 10, and 11, and its measured RM (consistent with zero across all days) demonstrates the stability of the polarization measurements on M87 and on Sgr A* on the same days. On Apr 7, the polarization calibrator was J1924-2914, which also shows consistent RM across days (see Figure 18, upper panels). Collectively, they demonstrate the stability of the polarization measurements in Apr 2017 on Sgr A*, M87, and, by extension, on the remaining AGN observed at λ1.3 mm (displayed in Figure 18) and at λ3 mm (displayed in Figure 19).
We notice that, for some of the targets, the RM fits are almost perfectly consistent with the measured EVPAs; in particular, the plots for Sgr A*, 3C 279 (Fig. 17), 3C 273, J1924-2914 (Fig. 18), and M87 (Fig. 20), look very unlikely given the error bars. In fact, a standard χ 2 analysis for these sources would give values close to 0 while for 4 points/2-parameter fit, we would expect χ 2 ∼2. This apparently unexpected behaviour can be explained by the fact that the EVPA uncertainties displayed in the RM plots are not just the thermal errors (which would naturally introduce scatter in the measure-ments) but also include a systematic error (0.03% of I into QU errors), which in fact dominates the total error budget (especially for the strongest sources with high Stokes I). We assessed that once such systematic error is removed, the EVPA data points are no longer consistent with the line to within their (thermal-noise-only) uncertainties. The fact that the thermal-only error bars are too small, but the error bars in the plots of this paper that include the systematic 0.03% Stokes I leakage into the QU error budget are too large, suggests that there is a real systematic error in the EVPA measurements but that it is smaller than the ALMA standard value. The fact that we are being conservative in our error estimates should ensure that we are not over-interpreting our measurements.
Finally, we notice that, for some of the targets and/or on specific days, there are > 1σ deviations between the observed EVPA and the RM fit-predicted values (e.g., OJ287, 4C 01.28, J0006-0623, J0510+1800; see selected panels in Figs. 18 and 19). This may suggest either the presence of an additional systematic error not accounted for in our error budget or that the assumption of the EVPA λ 2 -dependence is not valid in some cases. In fact, in § 3.3.3, § 5.1.1, and § 5.2 we discuss the possibility that some of the observed Faraday rotation may be partly internal, which would imply a more complex RM model than assumed in our analysis. Without a wider wavelength coverage, we cannot distinguish between the two cases. Therefore, we will not consider any additional systematic error in the RM analysis presented in this paper.

G. CIRCULAR POLARIZATION AND ASSESSMENT OF THE ALMA POLARIMETRY
A reliable detection of Stokes V in interferometric observations with linear-polarization feeds (i.e., the case of ALMA) strongly depends on the correct estimate of two instrumental quantities. On the one hand, the relative phase, ∆, between the X and Y polarizers (i.e., the cross-polarization phase) at the reference antenna (i.e., the quantity stored in the XY0.APP table, see Goddi et al. 2019a). On the other hand, the imaginary parts of the Dterms that describe the polarization leakage of the individual ALMA elements (i.e., the quantities stored either in the Df0.APP or in the Df0.ALMA table, depending on the calibration strategy, as described in Goddi et al. 2019a). Following the standard ALMA calibration procedure, ∆ is estimated by assuming a negligible Stokes V in the polarization calibrator. With this assumption, and neglecting also the effects from polarization leakage, the cross-polarization correlations (XY * and Y X * ) between two ALMA antennas observing a polarized point source are where the data are assumed to be already corrected for phase and amplitude gains, ψ and φ are the feed angle of the antennas 18 and the EVPA of the observed source (respectively), and p is the calibrator's linearly-polarized flux density (p = LP × I).
The only complex quantity that appears in Eqs. G1 is given by the factor e j∆ . Hence, the XY * (and Y X * ) visibilities may change their amplitudes as a function of time (via the changes in ψ), but their phases will remain constant and equal to ∆. The CASA-based calibration algorithm for ALMA (provided in the polcal task) takes advantage of this fact, and estimates the value of ∆ from the fit of XY * and/or Y X * (over different values of ψ) to a model with a constant phase. We call the phase estimated in this way as ∆ QA2 . If the true value of the cross-polarization phase, ∆, is offset from ∆ QA2 by an unknown quantity, β (i.e., ∆ = β + ∆ QA2 ), this offset will introduce a leakage-like effect into the polconverted VLBI visibilities, which will be described by the Dterms given in Eq. 13 of Goddi et al. (2019a). If β is very small, that equation predicts an ALMA leakage for VLBI which has two remarkable properties: 1) it is pure imaginary, and 2) it is the same for the R and L polarization hands. The value of this Dterm is directly related to β via the equation: Therefore, if the QA2 calibration procedure has introduced an offset β into the cross-polarization phase, ∆ QA2 , we predict a pure imaginary Dterm in the ALMA-VLBI visibilities. The actual Dterms estimated from the ALMA-VLBI observations during the EHT 2017 campaign are reported in Event Horizon Telescope Collaboration et al. (2021a).
The phase offset β may also introduce a spurious V signal into the data, which can be described by the equation (e.g., Hovatta et al. 2019) where V tot is the total V signal recovered from the (QA2-calibrated) data and V const is a constant V signal, independent of ψ. If the QA2 Dterm estimates for ALMA are correct, V const is equal to the true Stokes V of the source, (i.e., V const = V true ). However, if the ALMA Dterms estimated in the QA2 are offset from their true values, there is another instrumental contri-   Figure 17. Faraday RM of Sgr A* (upper panels), M87 (middle panels), and 3C 279 (lower panels). The EVPA as a function of wavelength-squared is presented with 1σ error bars for each SPW and for each day. The EVPAs are measured in the 1.3 mm band, except for the upper left panel, presenting the EVPAs measured at 3 mm in Sgr A*. The error bars are derived adding in quadrature to the thermal uncertainty of the Q and U maps (1σ image RMS) a systematic error of 0.03% of Stokes I (error bars for Sgr A* are smaller than the displayed points). The line is a linear fit to the data giving the RM reported inside the box. Plots in each row span the same vertical axis range in degrees in order to highlight differences in slope (the upper left panel is an exception). It is remarkable that Sgr A* showcases a very consistent slope across days, while in M87 the slope appears to change sign. Note also the different EVPA values between Apr 6/7 and Apr 11 in Sgr A* and between Apr 5/6 and Apr 10/11 in M87 and 3C279. The measured RM in 3C 279 (used as polarization calibrator on Apr 5, 6, 10, and 11) is consistent with zero across all days, demonstrating the stability of the polarization measurements on M87 and on Sgr A* on the same days.
bution to XY * and Y X * , which couples to V const as following In this equation, D k X is the part of the Dterm of the X polarizer of antenna k that remains uncalibrated after the QA2 (and similarly for the Y polarizer). Therefore, if the effects of the antenna Dterms are not fully removed from the data, the calibrator source will show a nonnegligible V true , while the CASA model assumes it to be null.
In Fig. 21, we show the V tot signal (computed as the average real part of (XY * − Y X * )/j among all ALMA antennas) of the polarization calibrator 3C 279, as a function of parallactic angle, for the epochs where this source was observed. We have averaged the visibilities in time bins of 120 seconds and the data taken with antenna elevations below 30 degrees have been discarded. In the figure, we also show the simple model given by Eq. G2, where V const and ∆ are the only two free parameters used in the fit.
The V signals from the polarization calibrator show a clear dependence with parallactic angle, which indicates that there are offsets, β in the estimated crosspolarization phases, ∆ QA2 . Using the values of p and φ estimated for 3C 279 on different days, we fit an X-Y phase offset of β = 0.2 − 0.5 degrees in all tracks except April 5, where β ∼ 1 − 2 degrees. We notice that these ranges assume no bias in the QA2 estimates of the ALMA Dterms.
It is interesting to note that the data depart from the sinusoidal model of Eq. G2, especially for the epochs on April 6 and 11, and for observations far from transit. These deviations may be related to other instrumental effects (e.g., the second-order leakage contributions, like O(p D) in Eq. G3).
The estimated values of V const obtained from Eq. G2 are all of the order of 1 mJy, at most. These val-ues are small compared to the V true coming from the β estimates. This is especially true for the epoch of April 5. As discussed at the beginning of this section, a low V const can be explained by the compensating ef- . RM fits for the 2017 GMVA targets with EVPA measurements at 3 mm (see Table 1). fect of (small) biases in the QA2 estimates of the ALMA Dterms (Eq. G3). Such systematics would force V const to be close to zero, regardless of the value of V true . In summary, two clear conclusions can be drawn from Fig. 21. On the one hand, there is an offset, β, in the X-Y cross-polarization phase of the reference antenna. On the other hand, there may be other systematics in the estimates of the antenna Dterms that may compensate the imprint of a true circular polarization of the calibrator (since CASA always assumes a null Stokes V in the calibrator).

G.1.2. Other sources
The Dterms and cross-polarization phases discussed in the previous subsection are applied to all sources in the data. Hence, the V Stokes from the polarization calibrator, which may have introduced a cross-polarization phase offset, β, and other biases to the Dterm estimates, will be systematically put back (after applying the polarization calibration) into the Stokes V signals of the rest of sources.
Thus, if we find different values of the fractional V Stokes among different sources, there has to be a contribution to their V const values (Eqs. G2 and G3) that is independent of that introduced by the polarization calibration. In the frame of our modelling of the instrumental polarization, such a contribution would likely be related to V true , i.e., a true circular polarization associated to the sources. In Fig. 22, we show the fractional Stokes V values for all the sources, with the exception of 3C 279, as a function of feed angle. As with Fig. 21, a time binning of 120 seconds has been applied and data with elevations lower than 30 degrees have not been used. We also show the model fitted with Eq. G2, where β is fixed to the value derived from the 3C 279 data (so the only free parameter is V const ). The only exception to this modelling is for the epoch on April 7, where J1924-291 was used instead as the polarization calibrator.
Some sources show a clear dependence of V tot with feed angle, being 3C 273 (on April 6) and OJ287 (on April 5) remarkable examples. However, the model prediction for such a variability, based on the crosspolarization phase offsets, β, estimated from 3C 279, cannot reproduce all the data. A possible explanation for this discrepancy could be, for instance, a small (within 1 degree) variation of ∆ with pointing direction (i.e., antenna elevation and azimuth) and/or an effect related to residual Dterms (see Eq. G3). A deep analysis of these possibilities is out of the scope of this paper and should indeed be carried out at the Observatory level.
In any case, the fractional circular polarizations shown in Fig. 22 are very different among sources, which is a good indicative that these are not dominated by the Dterm systematics (at least, to a first-order approximation). Some sources do show the sinusoidal dependence with feed angle, whereas others (like M87) are dominated by V const . Since we do not know the exact effects related to residual Dterms, it is not possible to derive V true from the estimated V const . A robust conclusion, however, is that there is circular polarization detected in most of the sources (with amplitudes of the order of a few 0.1%) and that, to our understanding, the instrumental effects, alone, cannot explain the results for all sources in a self-consistent way. Apr 11 -Stokes V ; β = 0.35 • Figure 22. Reconstructed fractional Stokes V of all sources (but 3C 279), as a function of feed angle, during the EHT campaign. Stokes V is computed as the real part of (XY * − Y X * )/j. Continuum lines show a simplified model based on cross-polarization phase offsets at ALMA, which have been fixed to the β values shown on top of each figure. Note that Sgr A* is not plotted because it displays significant intrinsic variability on the timescales (of hours) plotted here (e.g., Bower et al. 2018 for polarization calibration purposes (see Goddi et al. 2019a).

H. SPECTRAL INDICES OF TOTAL INTENSITY
We compute the total intensity spectral index for all the sources observed in the VLBI campaign at 3 mm and 1.3 mm. For each source, the spectral index α, de-fined as I(ν) ∝ ν α , is derived "in-band", performing a weighted least-squares fit across the four flux-density values estimated with uvmultifit in each SPW, i.e. at frequencies of 213, 215, 227, 229 GHz in the 1.3 mm band, and 86, 88, 98, and 100 GHz in the 3 mm band, respectively. For Sgr A*, the spectral index of the compact core is around 0 both at 3 mm (α = 0.01 ± 0.1) and 1.3 mm (α = [−0.03, −0.15] ± 0.06). For M87, the spectral index of the compact core at 1.3 mm is negative (α=[-1.2,-1.1]). Cycle 0 ALMA observations at 3 mm at a comparable angular resolution (2.6 × 1.4 ) yields a much flatter α=[-0.2,-0.3] , consistent with VLA measurements at radio frequency bands (8.4-43 GHz). VLBA observations revealed a flat-spectrum (α=0) compact core (e.g., Kravchenko et al. 2020). The steeper spectrum measured in the 1.3 mm band suggests that a spectral break must occur between 3 mm and 1.3 mm. Spectral steepening is also observed in other AGN sources in the sample, which vary in the range α=[-0.7,-0.3] at 3 mm and α=[-1.3,-0.6] at 1.3 mm (Cen A being the only exception, with α=-0.2). When put together, these results indicates that the spectrum of AGN cores becomes progressively more optically thin at mm wavelengths (see § 5.1.1).
H.1. Foreground absorption at 226.91 GHz toward Sgr A* Figure 23 shows the spectrum of Sgr A* in SPW=2 on 2017 Apr 07. The spectrum is virtually the same in the remainder days of observation in Apr 2017, suggesting that the absorption is probably associated with material that is in front of the galactic center core. Cyanide Radical (CN) and its hyperfine structure at 226.3600 GHz (N=2-1, J=3/2-3/2), 226.6595 GHz (N=2-1, J=3/2-1/2), and 226.8748 GHz (N=2-1, J=5/2-3/2) 19 are identified as the carriers of the absorption features. The lines predict a loss of integrated emission over the 1.8 GHz bandwidth of about 2%, in reasonable agreement with the decrements seen in SPW=2 in Fig. 23.

I. TWO-COMPONENT POLARIZATION MODEL FOR M87
Here we present the two-component polarization model for M87 summarized in Section 5.2.2 and shown schematically in Figure 6. Details of the model itself, analysis, and parameter constraints can be found below. Note that unlike the detailed image models presented in Event Horizon Telescope Collaboration et al. (2021b), these models seek to reconstruct only the Stokes I, Q, and U integrated over the EHT map and within the ALMA core.
Bandpass for Sgr A* core (after removing instrumental bandpass with the calibrator NRAO 530) for 2017 Apr 07 observations. Only baselines longer than 60m were used to remove virtually all of the extended arcsecond emission. For each channel, every antenna bandpass amplitude is plotted (represented by different colors). All three April observations show virtually the same 'bandpass'.
The motivation of this modeling is to assess whether or not the significant interday variations seen in the RMs of M87 can be accommodated by a model in which the only variable element is the intrinsic polarization of the horizon-scale emission, holding all other properties of the Faraday rotation and a putative large-scale polarized component fixed. In summary, we find that it is possible to do so, though it does require the RM of the horizon-scale component to be significantly larger than the RMs associated with the ALMA core in M87 reported in Table 2.

I.1. Model Definitions
We consider two two-component models for the polarimetric properties of M87, differing in the location of the small-scale Faraday screen. For both components, we construct the set of Stokes I, Q, U , which may depend on observation day and wavelength. The directly compared quantities that comprise the model are the integrated Stokes parameters of the compact component, I com,day,λ , Q com,day,λ , and U com,day,λ , and of the combination of the compact and extended components, I tot,day,λ = I com,day,λ + I ext,λ , Q tot,day,λ = Q com,day,λ + Q ext,λ , and U tot,day,λ = U com,day,λ + U ext,λ . A summary list of the model parameters for the models described in Sections I.1.1-I.1.3 is contained in Table 13.

I.1.2. Compact component: external Faraday screen
For each day on which M87 was observed by ALMA and the EHT (i.e., April 5, 6, 10 and 11) we specify a similar model for the compact component. When the screen is assumed to be external, we adopt a similar model to the large-scale component with the exception that the rotation is due to both screens: I Ex com,day,λ = I 0,com (λ/λ 0 ) αcom Q Ex com,day,λ = m com,day I Ex com,day,λ cos 2ψ com,day + 2(RM com + RM ext )(λ 2 − λ 2 0 ) U Ex com,day,λ = m com,day I Ex com,day,λ sin 2ψ com,day + 2(RM com + RM ext (λ 2 − λ 2 0 ) . (I5) This introduces an additional eleven parameters: I 0,com , α com , RM com , and four each m com,day and ψ com,day . The polarization fraction and EVPA of the compact component, m com,day and ψ com,day , are distinct from all of the other parameters in the model in that they vary among observation days. Therefore, where useful, we will distinguish these as "dynamic" parameters, with the remainder of the parameters being "static" in the limited sense that they do not vary across the observation campaign.

I.1.3. Compact component: internal Faraday screen model
Many simulations of M87 indicate the presence of large Faraday depths in the emission region (Broderick & McKinney 2010;Mościbrodzka et al. 2017;Ricarte et al. 2020). Therefore, we also consider a simple model for the compact component in which the emission and rotation are co-located, i.e., an internal Faraday screen. In principle, this is inextricably linked to the detailed properties of the emission region. Here we employ the gross simplification of a single-zone, or slab, model: the emission and Faraday rotation within the compact component occurs within a homogeneous region. We begin with a summary of the polarimetric properties of such a slab.
For a plane-parallel source with physical depth L and at some reference wavelength λ 0 a uniform emissivity j, polarization fraction m em at emission, EVPA at emission ψ em , Faraday rotativity R, we have total intensity, and Stokes Q and U , jm em e 2i[Rxλ 2 +ψem] dx = m em I 2RLλ 2 sin(2ψ em + 2RLλ 2 ) − sin(2ψ em ) + i m em I 2RLλ 2 cos(2ψ em ) − cos(2ψ em + 2RLλ 2 ) = 1 2 m em Isinc(RLλ 2 ) cos(2ψ em + RLλ 2 ) + i sin(2ψ em + RLλ 2 ) (I7) from which we can immediately read off Q and U . The EVPA at the top of the slab is tan(2ψ) = Q U = tan 2ψ em + RL(λ 2 − λ 2 0 ) , from which it is apparent that the effective contribution to the compact RM is RM com = RL/2. The internal Faraday screen model is then defined by I In com,day,λ = I 0,com (λ/λ 0 ) αcom Q In com,day,λ = 1 2 m com,day I In com,day,λ sinc(2RM com λ 2 ) cos 2ψ com,day + 2(RM com + RM ext )λ 2 U In com,day,λ = 1 2 m com,day I In com,day,λ sinc(2RM com λ 2 ) sin 2ψ com,day + 2(RM com + RM ext )λ 2 , (I9) where, as with the external model, we have added the Faraday rotation from the large-scale Faraday screen. As with the external Faraday screen model, this introduces eleven parameters, though the interpretations of the polarization fraction, and EVPA subtly differ, here referring to those of the emission process instead of derotated values.
As with the external Faraday screen model, where useful, we will refer to m com,day and ψ com,day , which differ among observation days, as "dynamic" parameters to distinguish them from the remaining "static" parameters.

I.2. Markov Chain Monte Carlo Analysis
From the models described above and the integrated EHT Stokes parameter ranges presented in Table 7 in Appendix H2 of Event Horizon Telescope Collaboration et al. (2021a) (I EHT , Q EHT , U EHT ) and the ALMA core Stokes parameter values for the individual SPWs in Table 10 (I spw , Q spw , U spw ), we construct a log-likelihood for each set of model parameters. These comprise sixty data values in total: three (I EHT , Q EHT , U EHT ) on each of four days from the EHT observations (3 × 4 data points), three (I spw , Q spw , U spw ) in four SPWs on each of four days presented here (3 × 4 × 4 data points).
We assume that the integrated EHT Stokes parameters are distributed normally with means and standard deviations set by the centers and half-widths of the ranges; this likely over-estimates the true uncertainty on the I EHT , Q EHT , and U EHT . The resulting log-likelihood is We use 256 independent walkers, and run for 10 5 steps, discarding the first half of the chains. Explorations with fewer walkers and steps indicate that by this time the MCMC chains are well converged.
In addition to the models described in Appendix I.1, we also considered versions of the two-component model applied to each day independently, i.e., keeping only one day in the sum in Equation I10. In these, on each day the five parameters of the external screen and five parameters of the internal screen (a single m com and ψ com ) are independently fit on each observation day. Effectively, this corresponds to a forty-parameter model across the four observation days, permitting Faraday screens and emission from both the compact and extended Faraday screens to vary independently across the four observation days.
On any given observation day, the parameters are less well constrained in this case and, with the notable exceptions of the compact component polarization fraction and EVPA, are consistent with a single set of values across all days. The variable polarization properties of the compact component matches the expectation from the EHT measurements in Event Horizon Telescope Collaboration et al. (2021a). The consistency with a single set of values for the remaining model parameters serves as a partial motivation for the more restricted variability in the models presented in Appendix I.1.

I.3. Two-component model results
Excellent fits are found for both the external and internal Faraday screen models. For 44 degrees of freedom, the best-fit external and internal screen models have reduced χ 2 = 3.96 and χ 2 = 2.54, respectively, both modestly small and possibly indicating that the uncertainties on the integrated EHT Stokes parameters are indeed over estimated by their half-range values. Thus, it is possible to reproduce the variable polarimetric properties observed by ALMA with a model in which only the compact emission evolves. Figures 24 and 25 show the joint posteriors for the external (lower left) and internal (upper right) Faraday screen models for the static and dynamic model components, respectively. To facilitate a direct comparison with the external Faraday screen model, in Figure 25, the polarization fractions and EVPAs of the internal Faraday screen model have been depolarized and rotated to show the corresponding posteriors on their observed analogs. Model parameter estimates, marginalized over all other parameters, are contained in Table 13.
After adjusting the polarization fraction and EVPA of the compact component, the properties of the twocomponent models are consistent among the external and internal Faraday screen models. Strong correlations exists between many of the compact component features. These are very strong for the polarization fractions and EVPAs on neighboring observation days, i.e., April 5 and 6, and April 10 and 11. This is anticipated by the similarities in the integrated polarimetric properties reported in Event Horizon Telescope Collaboration et al. (2021a) between neighboring observation days, which naturally constrain the two-component polarization models accordingly.
The flux normalizations, polarization fractions, and EVPAs are well constrained. RM ext is restricted to small magnitudes in both cases, typically less than 1.5 × 10 5 rad m −2 , and remains consistent with vanishing altogether. In contrast, RM com is significantly nonzero, and typically of order −5 × 10 5 rad m −2 , factors of 3-10 larger than those reported in Table 2, implying a significant degree of competition between the spectral variations in the two components and the residual Faraday rotation. In neither model are the spectral indexes of the components well constrained.
By construction, these are necessarily consistent with the results from ALMA-only RM measurements reported in Table 2. Figure 26 shows a number of realizations drawn from the posteriors for the wavelengthdependence of the EVPAs on each of the observation days in comparison to the measured values listed in Table 10. The interday evolution in the ALMA RMs are well reproduced, despite restricting the variable elements of the model to the compact component. Within the ALMA SPWs the EVPAs only weakly depart from the linear dependence on λ 2 indicative of Faraday rotation; outside of the ALMA SPWs this divergence can become considerably larger, implying that additional coincident polarimetric measurements at longer For the internal Faraday screen model, the polarization fractions have been depolarized by a factor of sinc(2RMcomλ 2 0 ) and the EVPAs have been rotated by (RMcom + RMext)λ 2 0 , corresponding to observed values, concordant with the definition for the external Faraday screen model. Contours indicate 50%, 90%, and 95% quantiles.