# The $EB$-correlation in Resolved Polarized Images: Connections to
Astrophysics of Black Holes

Razieh Emami Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden
Street, Cambridge, MA 02138, USA  Sheperd S. Doeleman Center for Astrophysics
$|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Black
Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138,
USA  Maciek Wielgus Max-Planck-Institut für Radioastronomie, Auf dem Hügel 69,
D-53121 Bonn, Germany  Dominic Chang Center for Astrophysics $|$ Harvard &
Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Black Hole Initiative
at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA  Koushik
Chatterjee Black Hole Initiative at Harvard University, 20 Garden Street,
Cambridge, MA 02138, USA  Randall Smith Center for Astrophysics $|$ Harvard &
Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Matthew Liska Center
for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA
02138, USA  James F. Steiner Center for Astrophysics $|$ Harvard &
Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Angelo Ricarte Black
Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138,
USA  Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA  Ramesh Narayan Center for Astrophysics $|$ Harvard &
Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Black Hole Initiative
at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA  Grant
Tremblay Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA  Douglas Finkbeiner Center for Astrophysics $|$
Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Lars
Hernquist Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA  Chi-Kwan Chan Steward Observatory and Department of
Astronomy, University of Arizona, 933 N. Cherry Ave., Tucson, AZ 85721, USA
Data Science Institute, University of Arizona, 1230 N. Cherry Ave., Tucson, AZ
85721, USA  Program in Applied Mathematics, University of Arizona, 617 N.
Santa Rita, Tucson, AZ 85721  Lindy Blackburn Center for Astrophysics $|$
Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Black Hole
Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA
Ben S. Prather CCS-2, Los Alamos National Laboratory, P.O. Box 1663, Los
Alamos, NM 87545, USA  Paul Tiede Center for Astrophysics $|$ Harvard &
Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Black Hole Initiative
at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA  Avery E.
Broderick Perimeter Institute for Theoretical Physics, 31 Caroline Street
North, Waterloo, ON, N2L 2Y5, Canada  Department of Physics and Astronomy,
University of Waterloo, 200 University Avenue West, Waterloo, ON, N2L 3G1,
Canada  Mark Vogelsberger Department of Physics, Kavli Institute for
Astrophysics and Space Research, Massachusetts Institute of Technology,
Cambridge, MA 02139, USA  Charles Alcock Center for Astrophysics $|$ Harvard &
Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Freek Roelofs Black
Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138,
USA  Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA

(May 2023)

###### Abstract

We present an in-depth analysis of a newly proposed correlation function in
visibility space, between the $E$ and $B$ modes of the linear polarization,
hereafter the $EB$-correlation, for a set of time-averaged GRMHD simulations
compared with the phase map from different semi-analytic models as well as the
Event Horizon Telescope (EHT) 2017 data for M87* source. We demonstrate that
the phase map of the time-averaged $EB$-correlation contains novel information
that might be linked to the BH spin, accretion state and the electron
temperature. A detailed comparison with a semi-analytic approach with
different azimuthal expansion modes shows that to recover the morphology of
the real/imaginary part of the correlation function and its phase, we require
higher orders of these azimuthal modes. To extract the phase features, we
propose to use the Zernike polynomial reconstruction developing an empirical
metric to break degeneracies between models with different BH spins that are
qualitatively similar. We use a set of different geometrical ring models with
various magnetic and velocity field morphologies and show that both the image
space and visibility based $EB$-correlation morphologies in MAD simulations
can be explained with simple fluid and magnetic field geometries as used in
ring models. SANEs by contrast are harder to model, demonstrating that the
simple fluid and magnetic field geometries of ring models are not sufficient
to describe them owing to higher Faraday Rotation depths. A qualitative
comparison with the EHT data demonstrates that some of the features in the
phase of $EB$-correlation might be well explained by the current models for BH
spins as well as electron temperatures, while others may require a larger
theoretical surveys.

M87* – Linear Polarization – EHT – E-mode – B-mode – Spin

††journal: ApJ

##  1 Introduction

The recently published image of the giant supermassive black hole (SMBH) at
the center of the elliptical galaxy Messier 87 (hereafter M87*; Event Horizon
Telescope Collaboration et al., 2019a, b, c, d, e, f) by the Event Horizon
Telescope Collaboration (EHTC) reveals an asymmetric ring-like structure with
a diameter of 42 $\pm$ 3 $\mu$as, consistent with the predicted shadow size of
a BH from the Einstein theory of general relativity (GR) (Event Horizon
Telescope Collaboration et al., 2019e, f). The ring-like structure appears to
be a persistent image feature on timescales of years, further corroborating
this interpretation (Wielgus et al., 2020). The EHT image is generally
consistent with the common picture in which the environment of M87* is made of
a hot, relativistic, and magnetized plasma (Rees et al., 1982; Narayan & Yi,
1995; Narayan et al., 1995). However, from the intensity map, it is not
possible to constrain the morphology of the magnetic field nor identify
whether the compact emission might have originated from inflowing matter or
rather from an outflow, corresponding to thick accretion or jet/wind-based
emission, respectively.

The spatially resolved linear polarization pattern reported by the Event
Horizon Telescope Collaboration et al. (2021a, b) provides very useful
information about the orientation of the electric vector polarization angle
(EVPA). Furthermore, it also reveals an (almost) azimuthally symmetric pattern
for EVPAs. This motivated Palumbo et al. (2020a) to propose a particular
decomposition for the linear polarization in terms of the azimuthal modes, as
specified with complex coefficient $\beta_{m}$. They showed that $m=2$
provides the dominant contribution in the EVPAs. Moreover, they found that the
amplitude and phase of the $\beta_{2}$ coefficient are linked to the magnetic
field geometry and BH spin. However, they did not characterize the key driver
of $\beta_{2}$. Emami et al. (2022) made a comprehensive study of the origin
of the twisty patterns in linear polarization using a set of spatially
resolved (time-averaged) images of General Relativistic Magneto-Hydro
Dynamical (GRMHD) simulations and identified the magnetic field geometry as
the key contributor to creating the twisty patterns in linear polarization, as
approximated with $\beta_{2}$. While $\beta_{2}$ well-characterizes the
salient rotationally-symmetric features of the linearly polarized structure,
it neglects other asymmetric features of the image.

In this work, we propose a new correlation function, the $EB$-correlation in
visibility space for the EHT 2017 data as well as the time-averaged images of
GRMHD simulations. In the latter case, we use different approaches including
both a full-numerical algorithm as well as a semi-analytic method, comparing
the real and the imaginary parts of the correlation function using different
cut-offs in the azimuthal expansion of the $E$ and $B$ modes. We find that
higher order terms are indeed required in the azimuthal expansion of $E$ and
$B$ modes. We use a Zernike reconstruction algorithm to recover the
$EB$-correlation phase map and propose an empirical metric to extract features
from the phase map. Our analysis demonstrates that models that are
qualitatively similar exhibit different patterns which may be eventually
useful for recovering the BH spin. We use a family of different geometrical
ring models, each with different $B$ and $V$ field morphologies, and identify
a case that almost captures the features in both the image and visibility
spaces. We make a qualitative comparison between the $EB$-correlation function
from the EHT 2017 data as well as the time-averaged images of GRMHD
simulations and show some of the features in the phase map can be well
described by the limited library of images, as considered in this study, while
others may require additional explorations as is left for some future studies.

The paper is structured as follows. Section 2 presents the key steps in
defining the $EB$-correlation in the visibility space. In Section 3 we show an
immediate application to the EHT 2017 data. Section 4 constructs the
$EB$-correlation phase for a set of few toy models. In Section 5 we present
our numerical tools for making the polarized BH images. Section 6 introduces
the Zernike polynomials to extract the features from the correlation phase. In
Section 7 we infer the correlation phase for the time-averaged images. Section
8 presents a semi-analytic model based on m-order expansions of the $E$ and
$B$ modes and compares them with the full-numerical analysis. In Section 9 we
present a geometrical ring model that best captures both the time-averaged
images and the $EB$-correlation phase. Section 10 provides a comparison with
the EHT 2017 data. In Section 11 we present the conclusion of the paper.
Section 12 discusses future directions for this new correlation function. We
present some technical details in Appendices A-C.

##  2 $EB$-correlation function in the visibility space

In the following, we define the $E$ and $B$ modes and construct their cross-
correlation function in visibility space.

###  2.1 Definition: $E$ mode and $B$ mode

As the first step, here we construct the flat sky based $E$ and $B$ modes in
visibility space. We follow the notation of Kamionkowski & Kovetz (2016);
Event Horizon Telescope Collaboration et al. (2021b) in which $E$ and $B$
modes are defined naturally in the visibility space $(u,v)$, where $u$ and $v$
are sampled based on the interferometer vector. In this representation, the
$E$ and $B$ visibilities (note that throughout our analysis we have dropped
the tilde from $E$ and $B$ modes as they are naturally defined in the
visibility space) are related to the Stokes parameters
$\tilde{\mathcal{Q}}(u,v)$ and $\tilde{\mathcal{U}}(u,v)$ by a local rotation
of $2\theta$ in the Fourier space:

$\begin{pmatrix}E(u,\theta)\\\ \ B(u,\theta)\end{pmatrix}=\begin{pmatrix}\cos
2\theta&\sin 2\theta\\\ -\sin 2\theta&\cos
2\theta\end{pmatrix}\begin{pmatrix}\tilde{\mathcal{Q}}(u,\theta)\\\
\tilde{\mathcal{U}}(u,\theta)\end{pmatrix},$ (1)

where the rotation angle is defined as 111In this paper, we use the convention
taken in the ehtim package (Chael et al., 2018).:

$\theta\equiv\arctan{\left(\frac{u}{v}\right)}.$ (2)

In addition, in our analysis, we use the relation between the real and the
Fourier space quantities,

$\tilde{X}(u,\theta)=\int\int X(\rho,\phi)e^{-2i\pi\rho
u\cos{(\theta-\phi)}}\rho d\rho d\phi,$ (3)

where $X$ is a generic function in the image plane, e.g. the Stokes
parameter.222Notice that in the image plane, we have defined $\phi$ from $y$
axis rather than $x$. .

While for our analysis we merely focus on the $E$ and $B$ modes in the
visibility space, it is important to also comment on their definitions in real
space; where the $E$ and $B$ modes refer to the gradient and curl components
of the polarization field, respectively, (see Eqs. (A2) and (A3) of Event
Horizon Telescope Collaboration et al., 2021b, for more details). Quite
importantly, while the real space $Q$ and $U$ depend on the choice of the
coordinate system, the real space $E$ and $B$ modes are scalar quantities;
i.e. invariant under the change of the coordinate system.

Below, we make an in-depth analysis of constructing the $EB$-correlation
function using different approaches; aiming to extract information about
plasma astrophysics.

###  2.2 Construction: $EB$-correlation function

We define the normalized $EB$-correlation function in visibility space as:

$\rho_{\mathcal{EB}}(u,v)\equiv\frac{E(u,v)B^{*}(u,v)}{\sqrt{E(u,v)E^{*}(u,v)}\sqrt{B(u,v)B^{*}(u,v)}}.$
(4)

As this is a complex quantity, we infer its phase from:333To capture the sign
of the real and the imaginary part of the $EB$-correlation, in practice, we
use the $\mathrm{arctan2(imaginary,real)}$ function from the python numpy
library.

$\theta(u,v)\equiv\arctan{\left(\frac{Im(\rho_{\mathcal{EB}}(u,v))}{Re(\rho_{\mathcal{EB}}(u,v))}\right)}.$
(5)

While Eq. (5) specifies the correlation phase for a given function
$\rho_{\mathcal{EB}}(u,v)$, it does not immediately provide an algorithm for
how to compute this function for the time-averaged images of GRMHD
simulations. This is important as computing the phase average may be ill-
defined due to phase wrapping (Event Horizon Telescope Collaboration et al.,
2021a). To avoid this issue, below we introduce two different metrics for
calculating the correlation phase from the time-averaged images of GRMHD
simulations.

Figure 1:  The polarized image of M87* from EHT 2017 data on April 6 and April
11 as well as the $EB$-correlation phase in the visibility space (bottom).

####  2.2.1 Time-averaged Correlation phase: First method

In this method, we first compute the time-averaged $E$ and $B$ modes from the
GRMHD simulations. We subsequently infer the correlation function using Eq.
(5). The phase is therefore given by:

$\bar{\theta}_{\mathrm{1st}}(u,v)=\arctan{\left(\frac{Im(\rho_{\mathcal{\langle
E\rangle\langle B\rangle}}(u,v))}{Re(\rho_{\mathcal{\langle E\rangle\langle
B\rangle}}(u,v))}\right)},$ (6)

where in Eq. (6), $\mathcal{\langle E\rangle}$ and $\mathcal{\langle
B\rangle}$ refer to the time-averaged $E$ and $B$ modes, respectively.
Throughout this paper, we use this method as the main approach for computing
the $EB$-correlation phase. Hereafter we call this the $M_{1}$ method.

####  2.2.2 Time-averaged Correlation phase: Second method

In this approach, we first calculate the time-averaged correlation function
and subsequently infer the correlation phase using Eq. (5):

$\bar{\theta}_{\mathrm{2nd}}(u,v)=\arctan{\left(\frac{\langle
Im(\rho_{\mathcal{EB}}(u,v))\rangle}{\langle
Re(\rho_{\mathcal{EB}}(u,v))\rangle}\right)},$ (7)

where in Eq. (7), $\langle Im(\rho_{\mathcal{EB}}(u,v))\rangle$ refers to the
time-averaged correlation function.

As stated above, in what follows we mainly focus on the first approach.
However, in Section 7.3 we infer the phase adopting the second approach to
make a comparison and justify that differences are not substantial. Hereafter
we call this the $M_{2}$ method.

##  3 $EB$-correlation for EHT 2017 data

In this section, we take the first look at the phase of the $EB$-correlation
function using the fully polarized EHT 2017 data for M87*. We focus on the
observations performed on April 6 and April 11 in 2017. The data set was first
calibrated using the EHT-HOPS data pipeline (Blackburn et al., 2019; Event
Horizon Telescope Collaboration et al., 2019c), with additional polarization
and leakage calibration as detailed in Event Horizon Telescope Collaboration
et al. (2021a); Issaoun et al. (2022); Jorstad et al. (2023). There were minor
updates to the EHT calibration pipeline with respect to Event Horizon
Telescope Collaboration et al. (2019a); Event Horizon Telescope Collaboration
et al. (2021a), as described in Event Horizon Telescope Collaboration et al.
(2022a).

Figure 1 presents the intensity map with the EVPA ticks in the image plane
(top row) as well as the $EB$-correlation phase in the visibility space
(bottom row), respectively. Distinct features are visible in the discrete
phase maps. These features appear generally consistent between April 6 and
April 11. Motivated by this, in the remainder of this paper, we make synthetic
data, using the time-averaged images of the GRMHD simulations of H-AMR using
the EHT coverage from April 11 and we perform an in-depth analysis of the
$EB$-correlation function, comparing results with those derived from EHT 2017
data.

Figure 2: A family of different Gaussian models. For a single Gaussian model,
top panel, the phase of the $EB$-correlation is either zero or 180 deg. In a
double Gaussian model, middle panel, with the same EVPA ticks, the phase
remains the same despite the brightness asymmetry. Finally, for a double
Gaussian model with different EVPA ticks, the bottom panel, the phase gets
non-vanished owing to the polarization anisotropies.

##  4 $EB$-correlation: construction of toy models

As the first theory step, here we study the $EB$-correlation function for a
set of different geometrical toy models, including different Gaussian models
as well as a subset of geometrical ring models.

###  4.1 Gaussian models

To get an intuition on what may drive non-zero correlation phase, here we
study a set of different Gaussian models. Our study covers a single Gaussian
(top row), a double Gaussian with the same ticks of the EVPA (middle row) and
a double Gaussian model with two different sets of the EVPA ticks (down row)
in Figure 2, respectively. It is seen that the phase of the $EB$-correlation
function remains zero or 180 deg for the first two cases while it deviates
from these values for the last case. This illustrates that the asymmetries in
the EVPA ticks, as against the asymmetries in the intensity, drive non-zero
$EB$-correlation phases for Gaussian models.

Figure 3: The $EB$-correlation phase for models with zero incliation and fluid
speed (top row), models with $\mathrm{i}$=17 deg while no fluid speed (middle
row) and models with $\mathrm{i}$=17 deg and non-zero fluid speed bottom row,
respectively. I each row, from the left to right we study different magnetic
field morphologies including a radial, toroidal and a mixed magnetic field,
respectively. In each panel we overlay the ring image, on the top right corer,
to build an intuition about the geometry of the ring. The phase is zero for
cases with zero inclination, while it is non-zero for models with non-zero
inclination. It is also degenerate between models with different B and v field
geometries.

###  4.2 Geometrical Ring models

Here we extend the toy model set by studying a set of geometrical ring models
with the spin of $a=+0.5$ and various geometries for the magnetic and the
velocity fields as well as the source inclination. Figure 3 presents
$EB$-correlation phase for models with zero incliation, hereafter
$\mathrm{i}$, and zero fluid speed, referred as v, (top row), models with
$\mathrm{i}$=17 deg while no fluid speed (middle row) and models with
$\mathrm{i}$=17 deg and non-zero fluid speed bottom row, respectively. In each
row, from the left to right we study different magnetic field morphologies
including a radial, toroidal and a mixed magnetic field, respectively. In each
panel we overlay the ring image, on the top right corer, to build an intuition
about the geometry of the ring.

From the plot it is seen that the phase of the $EB$-correlation almost
vanishes for the case with zero inclination and zero fluid speed, while it is
non-zero for non-zero inclination with no fluid motion. This illustrates the
importance of non-zero inclination in launching the phase. 444We have also
tried models with zero inclination and non-zero fluid speed and confirmed that
the phase remains almost zero as well.

Furthermore, it is inferred that different B-field geometries lead to distinct
phase patterns which is also distinct for the case with zero (middle row) and
nonzero (bottom row) fluid speed, respectively. Consequently, we argue that
the magnetic and velocity field geometries are both important in driving the
$EB$-correlation phase. In a future work we make more direct statement on the
role of changing the B and v field geometries in making various non-zero
phases.

##  5 Numerical methodology

In this section, we present our numerical modeling and describe how polarized
images are formed during the calculation of radiative transfer.

Figure 4: Time averaged images of MAD simulations with different BH spins.
From the top to bottom, we present the total intensity (I), fractional linear
polarization (P/I), and the Stokes Q parameter. As expected images show a
clear-handedness in their EVPAs.

###  5.1 GRMHD simulations

We employ numerical simulations in modeling the plasma flow, using an ideal
and GRMHD simulation set in the Kerr metric, with BH spin taken as a free
parameter (Koide et al., 1999; Gammie et al., 2003; Anninos et al., 2005; Del
Zanna et al., 2007), which is fixed at $a=(\pm 0.94,\pm 0.5,0)$. We integrate
the GRMHD equations in 3D using the H-AMR algorithm (Liska et al., 2022),
taking into account both a magnetically arrested disk (MAD) (Bisnovatyi-Kogan
& Ruzmaikin, 1974; Igumenshchev et al., 2003; Narayan et al., 2003) as well as
a standard and normal evolution (SANE) case (De Villiers et al., 2003; Gammie
et al., 2003; Narayan et al., 2012), with different magnetic fluxes and
different adiabatic indices: $\Gamma_{\mathrm{ad}}$ of 13/9 (MAD) and 5/3
(SANE), respectively. The initial conditions for the GRMHD simulations are
taken as a torus; i.e. with a constant angular momentum (Fishbone & Moncrief,
1976), where the orbital angular momentum is parallel or antiparallel with
respect to the BH spin. The initial torus is seeded by a weak and poloidal
magnetic field. The particular choice of the initial torus is motivated by the
emitting flow, near the event horizon, that remains in a steady state,
decoupled from the flow at larger distances. Upon starting the evolution, the
plasma experiences instabilities including the magnetorotational instability
(MRI) Balbus & Hawley (1992) as well as other instabilities such as the
magnetic Rayleigh-Taylor (RT) instability (Marshall et al., 2018), leading to
turbulence and the emergence of a low-density, highly magnetized bubbles in
MAD simulations. These play a very important role in driving angular momentum
transport. These instabilities lead to an inward accretion flow of the matter
directed toward the black hole. Subsequently, the plasma tends to a state with
(i) a mildly magnetized midplane, (ii) a coronal component in which the gas to
magnetic pressure $\beta\equiv P_{g}/P_{B}\simeq 1$, and (iii) a very strongly
magnetized funnel region near the BH poles with $\sigma\equiv B^{2}/(4\pi\rho
c^{2})\gg 1$.

###  5.2 BH Imaging

In making the BH images, we use the general relativistic radiative transfer
(GRRT) algorithm, implemented in ipole (Mościbrodzka & Gammie, 2018). Each
image is generated using a field of view (FOV) of 200 $\mu$as with a
resolution of 400 $\times$ 400 pixels. For the images, the impact of
synchrotron emission, self-absorption, Faraday rotation, and Faraday
conversion are all taken into account. Since the GRRT is not scale-invariant,
to perform the ray-tracing, we set the characteristic length scale ($r_{g}$)
(which is determined with the BH mass $M_{BH}$ as $L=GM_{BH}/c^{2}$ where $G$
and $c$ refer to the gravitational constant and the speed of light,
respectively) as well as the mass-density. Throughout our analysis, we fix the
source to be M87* with $M_{BH}=6.2\times 10^{9}M_{\odot}$, located at the
distance of $D=16.9$ Mpc from an observer on earth. Furthermore, we fix the
mass density by adopting the observed flux at 230 GHz to be $F_{\nu}$ = 0.5 Jy
(Event Horizon Telescope Collaboration et al., 2019d). Following Event Horizon
Telescope Collaboration et al. (2019e), we fix the source inclination at 17
(163) degrees for retrograde (prograde) spins, respectively. Finally, we
rotate the simulated images to match the observed position angle of the M87*
forward jet at -72 deg (e.g., Walker et al., 2018).

Since the current GRMHD simulations assume the same temperature for electrons
and ions, we incorporate the electron temperature through a post-processing
approach, in which the thermal equilibrium is replaced with a collisionless
plasma. Consequently, the electrons and ions may have different temperatures
(Shapiro et al., 1976; Narayan & Yi, 1995). Following the approach of
Mościbrodzka et al. (2016); Event Horizon Telescope Collaboration et al.
(2019e); Event Horizon Telescope Collaboration et al. (2021b), we modulate the
ratio of the ion-to-electron temperature as:

$\frac{T_{i}}{T_{e}}=R_{\mathrm{high}}\frac{\beta^{2}}{1+\beta^{2}}+R_{\mathrm{low}}\frac{1}{1+\beta^{2}}.$
(8)

where $\beta$ refers to the gas-to-magnetic pressure ratio, while
$R_{\mathrm{low}}$ and $R_{\mathrm{high}}$ are both free parameters. To reduce
the scope of this analysis, we restrict these parameters to
$R_{\mathrm{low}}=1$ and $R_{\mathrm{high}}=20$, respectively.

Figure 5: The time-averaged images of SANE simulations with different BH
spins. From top to bottom, we present the total intensity (I), fractional
linear polarization (P/I), and the Stokes Q parameter. As the system is
controlled by turbulence, there is not any clear-handedness in the EVPAs.

###  5.3 Time-averaged BH images

Having introduced the main set of GRMHD simulations as well as the approach
for generating polarized BH images, here we present the time-averaged images
of M87*. Figures 4 and 5 show the time-averaged intensity, fractional linear
polarization as well as the Stokes $Q$ parameter for the MAD and SANE
simulations, respectively. In each figure, from left to right, we increase the
BH spin as $a=(-0.94,-0.5,0.0,0.5,0.94)$. As expected, MAD simulations have
more coherent EVPA patterns and clearly exhibit organized structures. Figure 4
shows that the structure of the EVPAs (at the photon ring and outside) are
most similar among $a=(-0.94,+0.5,+0.94)$ versus among $a=(-0.5,0.0)$. As we
will see in the following, the phase of the $EB$-correlation function is also
similar for the aforementioned cases in MAD simulations. SANE simulations, on
the contrary, have mostly chaotic EVPA patterns without a clear handedness.
From Figure 5, it is inferred that their EVPAs differ substantially compared
with MAD simulations. Below, we show that it also leads to substantial
differences in their $EB$-correlation function.

###  5.4 Impact of Faraday Rotation

In the presence of a magnetized plasma, linear polarization is altered owing
to the Faraday rotation with an amount depending on the intervening plasma
density, path length, magnetic field, and the temperature. To assess the
impact of the Faraday rotation on the phase of the $EB$-correlation function,
we infer images with the Faraday rotation switched off (hereafter NFR), where
we turn off the Faraday Rotation coefficient, $\rho_{V}=0$, and compute the
time-averaged phase of the $EB$-correlation while comparing different models.
The third row in Figures 7 and 8 present the time-averaged phase map for the
NFR case in MAD and SANE simulations, respectively. It is inferred from the
figures that while the Faraday Rotation does not have any significant
contributions in MAD simulations, it plays an important role for SANEs.

Figure 6: The phase map for the $EB$-correlation function in the visibility
space from the time-averaged GRMHD simulations using the EHT 2017 coverage for
M87*. In the top(bottom) panel, we have used the MAD(SANE) simulations to make
the synthetic data.

##  6 Zernike polynomials

Here we present a new algorithm to reconstruct $EB$-correlation maps in
visibility space and show how to extract the BH features from our proposed
$EB$-correlation phase map. We employ the Zernike polynomials (Noll, 1976) in
a unit circle as a complete and orthogonal basis for the expansion. This
method is being widely used in adaptive optics (Lakshminarayanan & Fleck,
2011) to extract features from an image. Recently, it has also been used in
the galaxy cluster community (Capalbo et al., 2021). In our analysis, we
follow the convention in Noll (1976) for the positive and negative Zernike
polynomials as:

$\begin{split}Z_{n}^{m}(\rho,\theta)~{}=&~{}N_{n}^{m}R_{n}^{m}(\rho)\cos{m\theta},\\\
Z_{n}^{-m}(\rho,\theta)~{}=&~{}N_{n}^{m}R_{n}^{m}(\rho)\sin{m\theta}.\end{split}$
(9)

In Eq. (9) only the absolute value of $m$ appears. More explicitly, we have
$N_{n}^{m}=N_{n}^{|m|}$ and $R_{n}^{m}=R_{n}^{|m|}$ and the same in the
$\sin{m\theta}=\sin{|m|\theta}$. Furthermore, the normalization factor is
given by:

$N_{n}^{m}=\sqrt{\frac{2(n+1)}{1+\delta_{m0}}}.$ (10)

Also, the Zernike radial term is expressed as:

$R_{n}^{m}(\rho)=\sum_{s=0}^{(n-m)/2}\frac{(-1)^{s}(n-s)!}{s!\left(\frac{n+m}{2}-s\right)!\left(\frac{n-m}{2}-s\right)!}\rho^{n-2s},$
(11)

where $0\leq\rho\leq 1$ is the normalized radial distance, and the azimuthal
angle lies in the range $0\leq\theta\leq 2\pi$. Zernike polynomials are only
non-zero for $n-m=\mathrm{even}$. We may express any arbitrary function
$\phi(\rho,\theta)$ as the weighted sum of Zernike polynomials as:

$\phi(\rho,\theta)=\sum_{n=0}^{\infty}\sum_{m=0}^{n}c_{nm}Z_{n}^{m}(\rho,\theta),$
(12)

where $c_{nm}$ refers to the Zernike moments, given by:

$c_{nm}=\left(\frac{1}{\pi}\right)\int_{0}^{1}\int_{0}^{2\pi}\phi(\rho,\theta)Z_{n}^{m}(\rho,\theta)\rho
d\rho d\theta.$ (13)

In the above expansion, we have used the orthogonality of the Zernike
polynomials:

$\int_{0}^{1}\int_{0}^{2\pi}Z^{m}_{n}(\rho,\theta)Z^{m^{\prime}}_{n^{\prime}}(\rho,\theta)\rho
d\rho d\theta=\pi\delta_{nn^{\prime}}\delta_{mm^{\prime}}.$ (14)

Below, we use the Zernike polynomials as a novel way to reconstruct and
extract the features from the $EB$-correlation phase map. We explicitly show
that the Zernike polynomials might be very useful in distinguishing features
from maps that are qualitatively similar to one another. The generality of
these patterns must be checked against a wider library of images and is left
to a future study.

Figure 7: The time-averaged phase map for the $EB$-correlation for MAD
simulations (top-row), the reconstructed phase map using the Zernike
polynomials (middle-row) and for MAD simulations with No Faraday Rotation
(NFR) (bottom-row). To infer the phase map we adopt the $M_{1}$ method
presented in Section 2.2.1 in which we first compute the time-averaged $E$ and
$B$ modes and then compute the correlation phase using these averaged
quantities. The Zernike expansion is truncated at $n=40$ order. As the Zernike
polynomials are defined on a unit circle, we have multiplied them to 8 to be
similar to the original phase map from the GRMHD simulations.

Figure 8: The time-averaged phase map for the $EB$-correlation for SANE
simulations (top-row), the reconstructed phase map using the Zernike
polynomials (middle-row) and for SANE simulations with No Faraday Rotation
(NFR) (bottom-row). To infer the phase map we adopt the $M_{1}$ method
presented in Section 2.2.1 in which we first compute the time-averaged $E$ and
$B$ modes and then compute the correlation phase using these averaged
quantities. The Zernike expansion is truncated at $n=40$ order. As the Zernike
polynomials are defined on a unit circle, we have multiplied them to 8 to be
similar to the original phase map from the GRMHD simulations.

##  7 Phase correlation Analysis

In this section, we compute the $EB$-correlation phase using the synthetic
data generated from the time-averaged GRMHD simulations and by adopting the
EHT 2017 coverage for M87* (Event Horizon Telescope Collaboration et al.,
2021a). In our analysis, we investigate three different cases. First, we
calculate the $EB$-correlation phase using the discrete EHT coverage. We then
extend this analysis to grids in the $uv$-space. Finally, while in the first
two cases, we use the $M_{1}$ method in Section 2.2.1, as a consistency check,
we also perform the analysis using the $M_{2}$ method in Section 2.2.2.

###  7.1 $EB$-correlation phase: Discrete EHT coverage

We infer the $EB$-correlation phase map for the time-averaged GRMHD
simulations using the EHT 2017 coverage (Event Horizon Telescope Collaboration
et al., 2021a). In our analysis, we follow the $M_{1}$ method described in
Section 2.2.1.

Figure 6 presents the $EB$-correlation phase map for the synthetic
observational data made using the time-averaged images of GRMHD simulations,
adopting the EHT 2017 coverage for the M87* source. Despite the sparse $uv$
coverage, distinct patterns are visible from different columns. This motivates
us to perform an in-depth analysis making some grids in the visibility space.

###  7.2 $EB$-correlation phase: EHT gridded coverage

To get a better sense of different patterns in the $EB$-correlation phase map
as well as to circumvent the sparse $uv$-coverage in the EHT data, here we use
a grid of $uv$-space coverage to infer the time-averaged $E$ and $B$ modes. We
continue using the $M_{1}$ method in Section 2.2.1 to calculate the time-
averaged correlation phase.

Figure 9: The Zernike moments vs. the $n$ index for MAD (top) and SANE
simulations (bottom). From left to the right, we increase the BH spin
$a=(-0.94,-0.5,0.0,+0.5,+0.94)$. In each panel, we sum over all of the non-
zero $m$ indices for a given $n$ splitting the positive and negative terms
with the solid-magenta and dashed-green lines, respectively. Different
simulations show distinct patterns in $C_{n}$ (as defined in Eq. 15) vs. $n$
that can eventually be used to distinguish them from one another. In each
panel, the subplot refers to the case with an extended value of $n$s,
indicating that in most cases we do not need to consider n-orders above 20.

Figure 10: The Zernike moments vs. the absolute value of the $m$ index for MAD
(top) and SANE (bottom) simulations. From left to the right, we increase the
BH spin in $a=(-0.94,-0.5,0.0+0.5,+0.94)$. In each panel, we compute the
summed $C_{nm}$ for all $n$ indices that include the associated $m$ index,
splitting the positive and the negative terms with solid-magenta and dashed-
green lines, respectively. Different simulations show different patterns in
their $C_{|m|}$ (as defined in Eq. 16) vs $|m|$ that might be used to
distinguish them from one another. In each panel, the subplot refers to the
case with an extended value of $|m|$s that again indicate we do not need to
consider $m$ orders above 20.

Figures 7 and 8 present the $EB$-correlation phase map for time-averaged
simulations with different BH spins. In each figure, from the top to bottom,
different rows present the original simulation (MAD and SANE), the
reconstructed phase maps using the Zernike polynomials, and the phase map for
the case with No Faraday Rotation (NFR). In each row, from the left to right,
we increase the BH spin in $a=(-0.94,-0.5,0.0,+0.5,+0.94)$.

A dipolar pattern is visible in all different maps owing to the conjugate
symmetry in a baseline between the two stations $i$ and $j$. More explicitly,
it can be shown that $E^{*}_{ij}=E_{ji}$ and $B^{*}_{ij}=E_{ji}$, see Eq. 7 of
Event Horizon Telescope Collaboration et al. (2021c) for more details.
Consequently, the phase flips under $[u,v]\rightarrow[-u,-v]$. However, it
does not lead to any flips when we only flip either of $u$ or $v$ values in
the visibility space.

Focusing on MAD simulations, there are two distinct types of phase maps. In
the first case, which includes $a=(-0.94,+0.5,+0.94)$, the phase map contains
visibly distinct plus and minus sign phases, which are mainly centered around
$\theta\geq 90$ deg. Comparing the phase map with the EVPA patterns from
Figure 4, it is observed that there are some similarities between their EVPAs
as well. In the second case, which includes $a=(-0.5,0.0)$, the phase map is
mainly centered on rather small values, $\theta\leq 30$ deg. Comparing this
with the EVPA patterns from Figure 4, it is observed that their EVPAs show
very similar patterns as well.

For the SANE simulations, in the third-fourth rows, it is clearly seen that
various spins behave rather differently. Comparing this with the EVPA patterns
in Figure 5, it is observed that they have very different patterns as well. We
speculate that it might owe to the turbulence in the SANE simulations.

Figure 11: The time-averaged phase map for the $EB$-correlation for MAD and
SANE simulations using the $M_{2}$ method described in Section 2.2.2, where we
first compute the $EB$-correlation at every snapshot and then calculate the
time-averaged phase map. Comparing the results with Figures 7 and8, it is seen
that both maps are qualitatively similar. Figure 12: The time averaged phase
map for the $EB$-correlation in MAD simulations from the GRMHD simulations of
iharm. To infer the phase map we follow $M_{1}$ method described in Section
2.2.1. It is clearly seen that the phase map is very similar to that of the
MAD simulations from H-AMR simulations, which demonstrates that the observed
patterns in MAD simulations are robust and rather independent of initial
conditions.

Comparing the original phase maps, from the first row, with the reconstructed
ones from the Zernike polynomials, as shown in the second row, in Figures 7
and 8, it is demonstrated that the Zernike polynomials do a good job at
reconstructing the main features in the phase map. Consequently, in the
following, we will use this basis to make more quantitative comparisons
between the phase maps from different simulations.

####  7.2.1 Feature Extraction from the phase map using the Zernike
polynomials

Here we use the reconstructed $EB$-correlation phase map from Zernike
polynomials and extract BH features using two distinct empirical metrics.

$\blacksquare$ As the first metric, for every index $n$, we sum over the $m$
indices that appeared in the Zernike Moment $c_{nm}$, splitting the positive
(pos) and the negative (neg) $m$ indices:

$C_{n}\equiv\sum_{m}c_{nm}~{}~{}~{}~{}~{},~{}~{}~{}~{}~{}m=[\mathrm{pos},\mathrm{neg}].$
(15)

Figure 9 presents $C_{n}$ for MAD and SANE simulations. In each panel, we
present the $C_{n}$ vs. $n$ index in the range $0\leq n\leq 12$, as key
players in the expansion. To illustrate the impact of higher index $n$s, in
each panel, we add a sub-panel with $n\leq 40$. Finally, we illustrate the
positive and negative $m$ indices with solid-magenta and dashed-green lines,
respectively. From the plot, it is evident that simulations with different BH
spin as well as accretion types are distinct based on their $C_{n}$ profiles,
which also differ between the positive and negative $m$ indices. Furthermore,
in most cases, higher $n$ orders are suppressed compared with the lower values
(with an exception of a SANE simulation with $a=+0.94$). More interestingly,
it is seen that using the Zernike moments we may break some of the
degeneracies between the cases that are qualitatively very similar. For
instance, the positive $m$ index for $a=-0.5$ in MAD simulations is not the
same as for MAD with $a=0.0$. The same is true between MADs with spin
$a=-0.94,+0.5,+0.94$, which were qualitatively very similar. Consequently, we
argue that the $C_{n}$ vs. $n$ diagram might be very useful in distinguishing
different cases. The validity of this statement must be checked against
different electron temperature profiles and is left to a future study.

Figure 13: The comparison between the real part of the $EB$-correlation in MAD
simulations from the full-numeric as well as the $m$-order cut. From the top
to bottom, we present the full numeric, the case with $m_{\mathrm{cut}}$ =2,
$m_{\mathrm{cut}}$ =3 and $m_{\mathrm{cut}}$ =4, respectively.

$\blacksquare$ As the second metric, we sum over all $n$ indices that include
a given $m$ index, splitting out the positive (pos) and negative (neg) terms:

$C_{|m|}\equiv\sum_{n}c_{nm}~{}~{}~{}~{}~{},~{}~{}~{}~{}~{}m=[\mathrm{pos},\mathrm{neg}].$
(16)

Figure 10 presents $C_{|m|}$ vs. $|m|$ index from different simulations. In
all cases, the $C_{|m|}$ is suppressed at higher $|m|$ values. Furthermore,
each simulation has its own pattern which differs from other cases. This is
also true for degenerate cases. Combining this with the results from Figure 9,
it is observed that degeneracies between various models are broken if we use
the Zernike polynomials. This suggests that the combination of these two
metrics might be very informative in quantifying structural differences
between various cases.

###  7.3 $EB$-correlation: Second approach

So far we have only computed the $EB$-correlation phase map using the $M_{1}$
method from Section 2.2.1. Here we make a consistency check by computing the
phase map using the $M_{2}$ method as described in Section 2.2.2. Figure 11
presents the map using the $M_{2}$ method. Comparing the map with the one from
Figures 7 and 8, it is evident that they are quite similar. Motivated by this,
in the remainder of the paper we focus on the $M_{2}$ method.

Figure 14: Comparison between the real part of the $EB$-correlation in SANE
simulations from the full-numeric as well as the $m$-order cut. From the top
to bottom, we present the full numeric, the case with $m_{\mathrm{cut}}$ =2,
$m_{\mathrm{cut}}$ =3 and $m_{\mathrm{cut}}$ =4, respectively.

###  7.4 Comparison with iharm simulation

Here we make a comparison between the $EB$-correlation phase map inferred
using the H-AMR simulation, as done in this work, and the one from the
standard library of GRMHD simulations of iharm by Gammie et al. (2003);
Prather et al. (2021); Emami et al. (2022). Figure 12 presents the phase map
for the time-averaged GRMHD simulations of iharm. The phase map in the MAD
case is very similar to the one from H-AMR. However, it differs for SANE
simulations as the level of the turbulence is different between the iharm and
H-AMR simulations. Consequently, to avoid any confusions, we skip making the
comparison for SANE simulations.

##  8 Semi-analytic estimation of the $EB$-correlation function

Having presented the general approach in computing the
$\rho_{\mathcal{EB}}(u,v)$, here we follow a semi-analytical algorithm from
Event Horizon Telescope Collaboration et al. (2021b) and make an $m$-order
expansion for the $E$ and $B$ modes, in the visibility space. We then compare
that with the full numerical results as inferred from Section 7.2.

In this method, we expand the $E(u,v)$ and $B(u,v)$ modes as:

$\begin{split}E(u,v)=&\sum_{m=-\infty}^{\infty}i^{-m}Re\bigg{[}e^{i(m-2)\arctan{(u/v)}}W_{m}(\sqrt{(u^{2}+v^{2})})\bigg{]},\\\
B(u,v)=&\sum_{m=-\infty}^{\infty}i^{-m}Im\bigg{[}e^{i(m-2)\arctan{(u/v)}}W_{m}(\sqrt{(u^{2}+v^{2})})\bigg{]},\end{split}$
(17)

where $W_{m}(X)$ is defined as:

$\displaystyle W_{m}(X)\equiv$
$\displaystyle\int_{\rho_{\mathrm{min}}}^{\rho_{\mathrm{max}}}\int_{0}^{2\pi}\bigg{(}Q(\rho,\phi)+iU(\rho,\phi)\bigg{)}J_{m}(2\pi\rho
X)$ $\displaystyle~{}~{}~{}\times e^{-im\phi}\rho d\rho d\phi.$ (18)

Using Eq. (17) and following the $M_{1}$ method in Section 2.2.1, we first
compute the time-averaged $E$ and $B$ modes and then compute the
$EB$-correlation function. However, since the $E$ and $B$ expansions in Eq.
(17) includes infinite terms, we must truncate the expansion up to a specific
order. Former literature (see for example Event Horizon Telescope
Collaboration et al., 2021b) has only considered the linear polarization at
$m_{\mathrm{cut}}=2$ order. Here we go several steps beyond this
approximation, incorporating the new terms up to $m_{\mathrm{cut}}=\pm 4$
order, and explicitly check out whether the above expansion is sufficient at
recovering pattern morphology. As we will show, the $EB$-correlation function
is sensitive to these higher-order terms in the expansion. In Eq. (C), we
expand the $E$ and $B$ modes and compute the real and imaginary parts of the
$EB$-correlation function in Eqs. (C2) and (C3), respectively. It is evident
that while the real part of the $EB$-correlation is sensitive to even-even and
odd-odd terms, the imaginary part of the correlation function is only
sensitive to the even-odd terms and so the former approach adopted in Event
Horizon Telescope Collaboration et al. (2021b) does not work here.

Figure 15: The comparison between the imaginary part of the $EB$-correlation
in MAD simulations from the full-numeric as well as the $m$-order cut. From
the top to bottom, we present the full numeric, the case with
$m_{\mathrm{cut}}$ =2, $m_{\mathrm{cut}}$ =3 and $m_{\mathrm{cut}}$ =4,
respectively.  Figure 16: The comparison between the imaginary part of the
$EB$-correlation in SANE simulations from the full-numeric as well as with
different $m$-order cuts. From the top to bottom, we present the full numeric,
the case with $m_{\mathrm{cut}}$ =2, $m_{\mathrm{cut}}$ =3 and
$m_{\mathrm{cut}}$ =4, respectively.

Figures 13 and 14 present the real part of the $EB$-correlation function,
while Figures 15 and 16 show the imaginary part of the $EB$-correlation
function in MAD and SANE simulations, respectively. In each plot, from top to
bottom, we present the full-numerical results, as inferred using the $M_{1}$
method in Section 2.2.1, as well as the results of the expansion using
$m_{\mathrm{cut}}=2$, $m_{\mathrm{cut}}=3$ and $m_{\mathrm{cut}}=4$,
respectively.

From the plots, it is evident that in all cases, there is a considerable
contribution from modes beyond a truncation at $m_{\mathrm{cut}}=2$.
Furthermore, it is seen that in MAD simulations, the real part is almost well
represented up to $m_{\mathrm{cut}}=3$. However, even in MAD simulations
$m_{\mathrm{cut}}=3$ is not completely sufficient for the imaginary part of
the correlation function and we need to truncate it up to
$m_{\mathrm{cut}}=4$. The situation becomes more complicated for SANE
simulations where neither the real nor the imaginary parts can be sufficiently
described up to $m_{\mathrm{cut}}=3$ and we need to consider up to at least
$m_{\mathrm{cut}}=4$. There are also some cases where higher order terms are
required to fully describe the $EB$-correlation function, see for example
$a=(+0.5,+0.94)$ in Figures 14 and 16.

Figure 17: The geometrical ring model vs. the time-averaged GRMHD MAD
simulations. From the top to bottom we present the polarized images, the real
part and the phases of the $EB$-correlation for the ring model and the phase
from the time-averaged GRMHD simulation, respectively. It is inferred that the
ring model recovers the full phase features in MADs.  Figure 18: The
geometrical ring model vs. the time-averaged GRMHD SANE simulations. From the
top to bottom we present the polarized images, the real part and the phases of
the $EB$-correlation for the ring model and the phase from the time-averaged
GRMHD simulation, respectively. It is inferred that the ring model does not
seem to recover the full phase features in SANEs.

##  9 Constructing the geometrical Ring model

Having presented the $EB$-correlation map for the synthetic data from the
time-averaged GRMHD simulations of H-AMR, here we use a geometrical ring
model, with free magnetic field geometry as well as the velocity field
orientation in the equatorial plane. The aim is to vary the inclination of the
magnetic field (hereafter as $i_{B}$ ), its orientation in the equatorial
plane (specified with $\eta$) together with the velocity field orientation
(identified with $\chi$) which are defined as:

$\begin{split}i_{B}\equiv&\arctan{\left(B_{z}/\sqrt{B^{2}_{x}+B^{2}_{y}}\right)},\\\
\eta\equiv&\arctan{\left(B_{y}/B_{x}\right)},\\\
\chi\equiv&\arctan{\left(V_{y}/V_{x}\right)}.\\\ \end{split}$ (19)

Subsequently, we identify a subset of models that are well-matched with the
time-averaged images and the $EB$-correlation functions in visibility space.
Our analysis shows that out of a family of roughly 2000 different ring models,
only about a few percent of them demonstrate great similarities with the GRMHD
simulations. Table 1 presents the ring parameters associated with the
consistent model with the time-averaged GRMHD simulations.

| a = -0.94 | a = -0.5 | a = 0.0 | a = +0.5 | a = +0.94
---|---|---|---|---|---
Model | $i_{B}$ | $\eta$ | $\chi$ | $i_{B}$ | $\eta$ | $\chi$ | $i_{B}$ | $\eta$ | $\chi$ | $i_{B}$ | $\eta$ | $\chi$ | $i_{B}$ | $\eta$ | $\chi$
MAD | 59 | 0 | 108 | 59 | 0 | 144 | 59 | 0 | 144 | 59 | -72 | -144 | 59 | -108 | -108
SANE | 150 | 108 | 144 | 150 | 108 | 36 | 59 | 0 | -144 | 150 | -108 | -36 | 150 | 180 | -144
Table 1: The ring model parameters are similar to the time-averaged GRMHD
simulations of MAD (top) and SANE (bottom). From left to right we increase the
BH spin in $a=(-0.94,-0.5,0.0,+0.5,+0.94)$. In each case, we present the
magnetic field inclination, $i_{B}$, its orientation in the equatorial plane,
$\eta$, together with the velocity field orientation, $\chi$, in the
equatorial plane.

Figures 17 and 18 present a ring-model approximation to the time-averaged
GRMHD simulation for MAD and SANE cases, respectively. In each plot, from top-
to-bottom, we present the ring image, the real part, and the phase of the
$EB$-correlation from the ring model and the $EB$-correlation phase from the
time-averaged H-AMR GRMHD simulations, respectively. In each row, from left to
right, we increase the BH spin as $a=(-0.94,-0.5,0.0,+0.5,+0.94)$.

From the plots, it is inferred that the phase map in MAD simulations is almost
fully recovered while it is qualitatively more distinct for SANE models, which
is generally consistent with the results of Emami et al. (2022).

##  10 Comparison with the EHT 2017 data

As already stated above, this newly proposed $EB$-correlation function
provides a complementary constraint to the previous analyses. Using image-
integrated net linear polarization $|m|_{\mathrm{net}}$, image-integrated net
circular polarization $|v|_{\mathrm{net}}$, image-averaged linear polarization
$\langle|m|\rangle$, and a coherent azimuthal structure measurement
$\beta_{2}$ (proposed by Palumbo et al., 2020b), Event Horizon Telescope
Collaboration et al. (2021c) shows strong evidence that M87* is in a MAD
state. Specifically, SANE models tend to under produce $|m|_{\mathrm{net}}$
and $|\beta_{2}|$, and are disfavored. Therefore, we focus on MAD models here.

Having presented the EHT 2017 data for M87*, we make gridding in the
visibility space using the EHT reconstructed image and we compare it with the
time-averaged GRMHD simulations from H-AMR.

In Figure 19 we present the real, imaginary, and the phase of the
$EB$-correlation function for EHT 2017 on April 6 (1st-row), April 11 (2nd-
row), and then from the time-averaged GRMHD simulations with $a=0.0$ (3rd-row)
and $a=+0.5$ (4th-row), respectively.

From the figure it is inferred that the EHT phase map contains both of the
dark and bright regions, associated with low and intermediate phases. However
the map of MAD case with $a=0$ mostly contains dark regions while the map from
MAD with $a=+0.5$ merely contains bright regions. This demonstrates that
either the BH spin sits somewhere in between and/or the electron temperature
profile must be different in a favorite model. We leave this extra exploratory
investigation to a future work with a goal to extend over the range of models
and make a more direct comparison with the EHT data.

Figure 19: The real, imaginary, and the phase of the gridded $EB$-correlation
function for the EHT 2017 data for 3598 (1st-row) and 3601 (2nd-row) and the
time-averaged GRMHD simulations in MAD with $a=0.0$ and $a=+0.5$ spins,
respectively. The phase map in the EHT data contains both of the dark and
bright regions which is similar to the MAD with $a=0.0$ in the dark regions
while similar to MADs with $a=+0.5$ for bright regions.

##  11 Conclusions

We introduced a new correlation function in the visibility space, the
$EB$-correlation, and made an in-depth analysis of its real/imaginary parts as
well as the phase using the EHT 2017 data for the M87* source together with
the synthetic data from a subset of different toy models and the time-averaged
GRMHD simulations from H-AMR. In the following, we review our key findings.

$\bullet$ The phase map for the EHT 2017 data in April 6 and April 11 from
Section 3 illustrates some patterns in the $EB$-correlation phase in Figure 1
that includes both of dark and bright regions in the phase map.

$\bullet$ It is inferred from a set of different toy models considered in
Section 4 that asymmetries in the EVPAs, as shown in Figure 2, play important
roles in launching non-zero phases. Furthermore, the phase structure from
different geometrical ring models, shown in Figure 3 depends on different
factors including the source inclinations together with the magnetic and the
velocity field structures.

$\bullet$ The phase map from different GRMHD simulations as studied in Section
7.2 present distinct features, as shown in Figures 7 and 8, with some
qualitative similarities in some cases that are predominant in MAD
simulations.

$\bullet$ Employing the Zernike polynomials done in Section 6 as a new
reconstruction algorithm for the $EB$-correlation phase map breaks
degeneracies, shown in Figures 9 and 10, between models that are otherwise
qualitatively similar.

$\bullet$ An azimuthal expansion of the $E$ and $B$ modes in terms of
$m$-orders done in Section 8 demonstrates that to reconstruct both the real
and the imaginary parts of the $EB$-correlation function we require
considering higher order terms in the azimuthal expansion. More explicitly, in
most cases considering the expansion up to $m_{\mathrm{cut}}=4$, as shown in
Figures 15 and 16, recovers the majority of the features in the
$EB$-correlation function map.

$\bullet$ A set of geometrical ring models considered in Section 9, with
distinct magnetic and velocity field geometries, shows that only a small set
of models survive in recovering the patters in the $EB$-correlation function.
Consequently, we argue that the $EB$-correlation function is capable of
breaking the degeneracy between ring models that are otherwise very similar.

$\bullet$ A qualitative comparison between the $EB$-correlation from the EHT
2017 data with the time-averaged GRMHD simulations, done in Section 10,
demonstrates that the phase map from the EHT data contains, as shown in Figure
19 a mixture of both of the dark and bright regions which partially appear in
MAD simulations with $a=0.0$ (dark spots) while in MAD with $a=+0.5$ (bright
spots). This suggests that to find out the favorite model we might need to
change the BH spin as well as trying alternative electron temperature
profiles.

##  12 Future direction

While the current work focuses on a set of GRMHD simulations with 5 different
BH spins with a limited set of electron temperature profile, in a future study
we extend this analysis to also include more spins as well as different
electron temperature profiles. Furthermore, we aim to study cases with
electron temperature being computed using 2 temperature fluid simulations
(e.g., Sądowski et al., 2017). The impact of the a tilted accretion disk would
be also intriguing to be figured out.

## Data Availability

Data directly corresponding to this manuscript and the figures are available
to be shared on reasonable request from the corresponding author. The ray
tracing of the simulation done in this work was performed using the ipole
method (Mościbrodzka & Gammie, 2018). We have used the library of H-AMR
simulations by Liska et al. (2022) from the standard library of 3D time-
dependent GRMHD simulations in Event Horizon Telescope Collaboration et al.
(2022b).

## acknowledgement

It is a great pleasure to thank Geoffrey Bower, Andrew A. Chael, Michael
Johnson, Daniel Palumbo, George Wong and the entire of the EHT collaboration
for very fruitful conversations and comments. We also acknowledge the EHT
internal review panel for constructive comments which improved the quality of
this paper. Razieh Emami acknowledges the support from grant numbers
21-atp21-0077, NSF AST-1816420, and HST-GO-16173.001-A as well as the
Institute for Theory and Computation at the Center for Astrophysics. Dominic
Chang acknowledges the support of the Black Hole Initiative at Harvard
University. We thank the supercomputer facility at Harvard University where
most of the simulation work was done. This research was made possible through
the support of grants from the Gordon and Betty Moore Foundation and the John
Templeton Foundation. It was also supported by the National Science Foundation
grants AST 1935980, AST 2034306, and OISE 1743747. The opinions expressed in
this publication are those of the author(s) and do not necessarily reflect the
views of the Moore or Templeton Foundations.

Software: matplotlib (Hunter, 2007), numpy (van der Walt et al., 2011), scipy
(Oliphant, 2007), seaborn (Waskom et al., 2020), pandas (Reback et al., 2021),
h5py (de Buyl et al., 2016).

##  Appendix A Impact of the resolution

While the main text solely focuses on the unblurred time-averaged GRMHD
simulations, here we study the impact of adding the EHT finite resolution to
the $EB$-correlation phase. Figure 20 presents the time-averaged GRMHD
simulation on top of the EHT coverage. The EHT finite resolution of 20 $\mu$as
is taken into account in computing the time-averaging. Comparing this with the
results from Figure 6, it is seen that the phase map is very similar between
these simulations as all what blurring do is to only shift the phase in the
long baselines with no impacts for the short baselines. Consequently, the EHT
resolution does not affect the $EB$-correlation phase.

Figure 20: The phase map for the $EB$-correlation for the time-averaged GRMHD
simulations using the EHT coverage for M87* and taking into account the EHT
finite resolution of 20 $\mu$as. Comparing this with Figure 6 demonstrates no
noticeable differences.

Figure 21: The impact of time-variability as a variance in the phase map of
the time-averaged $EB$-correlation map for MAD and SANE simulations with
different BH spins.

##  Appendix B Impact of the time-variability in the mean phase

In the first order, the time variability in GRMHD simulations leads to an
averaged change in the correlation phase map. To model this effect, we infer
the phase variance and compute the phase map with this extra change being
incorporated. Here we estimate the phase variance as:

$\sigma(u,v)=f(a,b)\sqrt{\left(\frac{\partial f(a,b)}{f\partial a}\Delta
a\right)^{2}+\left(\frac{\partial f(a,b)}{f\partial b}\Delta b\right)^{2}},$
(B1)

where we have defined:

$\begin{split}f(a,b)\equiv&\arctan{\left(\frac{\langle
Im(\rho_{\mathcal{EB}}(u,v))\rangle}{\langle
Re(\rho_{\mathcal{EB}}(u,v))\rangle}\right)},\\\
a\equiv&Im(\rho_{\mathcal{EB}}(u,v)),\\\
b\equiv&Re(\rho_{\mathcal{EB}}(u,v)).\end{split}$ (B2)

More explicitly, we compute the phase map as:

$\Theta=\bar{\theta}(u,v)\pm\sigma(u,v).$ (B3)

Figure 21 presents the $EB$-correlation phase map including the above
variance. Descending rows refer to MADs and SANEs with phase variance being
subtracted and added to the mean phase, respectively.

From the plot, it is inferred that different simulations with various BH spins
lead to substantially different phase maps with distinct features between the
MAD and SANEs as well as different spins. Consequently, we argue that the
phase of the $EB$-correlation could be a very good parameter to probe the BH
spin as well as the accretion types near the BHs. Furthermore, the inclusion
of the phase variance does not change the shape of the phase significantly.
This means that at the leading order, we are not dominated by noise.

##  Appendix C Phase estimation from $\beta$ expansion

In this appendix, we walk through some of the steps in estimating the phase of
the $EB$-correlation. Since the $E$ and $B$ modes are suppressed for higher
values of $m$s, we truncate the expansion at $m_{\mathrm{cut}}=\pm 4$. The $E$
and $B$ modes can then be expressed as:

$\displaystyle E=$ $\displaystyle
E^{r}_{-4}+E^{i}_{-3}+E^{r}_{-2}+E^{i}_{-1}+E^{r}_{0}+E^{i}_{1}+E^{r}_{2}+E^{i}_{3}+E^{r}_{4},$
$\displaystyle B=$ $\displaystyle
B^{r}_{-4}+B^{i}_{-3}+B^{r}_{-2}+B^{i}_{-1}+B^{r}_{0}+B^{i}_{1}+B^{r}_{2}+B^{i}_{3}+B^{r}_{4}.$
(C1)

where the super-index $r$ and $i$ refer to the real and imaginary parts of
individual modes, respectively. While the sub-index numbers describe the
$m$-order in our expansion.

From Eq. (C) it is evident that the phase of the $EB$-correlation function
will be non-zero only if we have non-equal $m$-modes in $E$ and $B$.
Consequently, we split the $EB$-correlation into its real and imaginary
components. Where the real part is computed as:

$\displaystyle E(u,v)B^{*}(u,v)\bigg{|}_{\mathrm{real}}$
$\displaystyle=\bigg{(}E^{r}_{-4}+E^{r}_{-2}+E^{r}_{0}+E^{r}_{2}+E^{r}_{4}\bigg{)}\bigg{(}B^{r}_{-4}+B^{r}_{-2}+B^{r}_{0}+B^{r}_{2}+B^{r}_{4}\bigg{)}-\bigg{(}E^{i}_{-3}+E^{i}_{-1}+E^{i}_{1}+E^{i}_{3}\bigg{)}\bigg{(}B^{i}_{-3}+B^{i}_{-1}+$
(C2) $\displaystyle B^{i}_{1}+B^{i}_{3}\bigg{)},$

while the imaginary component is computed as:

$\displaystyle E(u,v)B^{*}(u,v)\bigg{|}_{\mathrm{imag}}$
$\displaystyle=-\bigg{(}E^{r}_{-4}+E^{r}_{-2}+E^{r}_{0}+E^{r}_{2}+E^{r}_{4}\bigg{)}\bigg{(}B^{i}_{-3}+B^{i}_{-1}+B^{i}_{1}+B^{i}_{3}\bigg{)}+\bigg{(}E^{i}_{-3}+E^{i}_{-1}+E^{i}_{1}+E^{i}_{3}\bigg{)}\bigg{(}B^{r}_{-4}+B^{r}_{-2}+$
(C3) $\displaystyle B^{r}_{0}+B^{r}_{2}+B^{r}_{4}\bigg{)}.$

As already stated above, the $EB$-correlation function mixes the even and odd
components of the $E$ and $B$ modes. As is explained in the main text, it
requires us to include higher order terms in the $E$ and $B$ mode expansion.

## References

  * Anninos et al. (2005) Anninos, P., Fragile, P. C., & Salmonson, J. D. 2005, ApJ, 635, 723, doi: 10.1086/497294
  * Balbus & Hawley (1992) Balbus, S. A., & Hawley, J. F. 1992, ApJ, 400, 610, doi: 10.1086/172022
  * Bisnovatyi-Kogan & Ruzmaikin (1974) Bisnovatyi-Kogan, G. S., & Ruzmaikin, A. A. 1974, Ap&SS, 28, 45, doi: 10.1007/BF00642237
  * Blackburn et al. (2019) Blackburn, L., Chan, C.-k., Crew, G. B., et al. 2019, ApJ, 882, 23, doi: 10.3847/1538-4357/ab328d
  * Capalbo et al. (2021) Capalbo, V., De Petris, M., De Luca, F., et al. 2021, MNRAS, 503, 6155, doi: 10.1093/mnras/staa3900
  * Chael et al. (2018) Chael, A. A., Johnson, M. D., Bouman, K. L., et al. 2018, ApJ, 857, 23, doi: 10.3847/1538-4357/aab6a8
  * de Buyl et al. (2016) de Buyl, P., Huang, M.-J., & Deprez, L. 2016, arXiv e-prints, arXiv:1608.04904.  https://arxiv.org/abs/1608.04904
  * De Villiers et al. (2003) De Villiers, J.-P., Hawley, J. F., & Krolik, J. H. 2003, ApJ, 599, 1238, doi: 10.1086/379509
  * Del Zanna et al. (2007) Del Zanna, L., Zanotti, O., Bucciantini, N., & Londrillo, P. 2007, A&A, 473, 11, doi: 10.1051/0004-6361:20077093
  * Emami et al. (2022) Emami, R., Ricarte, A., Wong, G. N., et al. 2022, arXiv e-prints, arXiv:2210.01218, doi: 10.48550/arXiv.2210.01218
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2019a, ApJ, 875, L1, doi: 10.3847/2041-8213/ab0ec7
  * Event Horizon Telescope Collaboration et al. (2019b) —. 2019b, ApJ, 875, L2, doi: 10.3847/2041-8213/ab0c96
  * Event Horizon Telescope Collaboration et al. (2019c) —. 2019c, ApJ, 875, L3, doi: 10.3847/2041-8213/ab0c57
  * Event Horizon Telescope Collaboration et al. (2019d) —. 2019d, ApJ, 875, L4, doi: 10.3847/2041-8213/ab0e85
  * Event Horizon Telescope Collaboration et al. (2019e) —. 2019e, ApJ, 875, L5, doi: 10.3847/2041-8213/ab0f43
  * Event Horizon Telescope Collaboration et al. (2019f) —. 2019f, ApJ, 875, L6, doi: 10.3847/2041-8213/ab1141
  * Event Horizon Telescope Collaboration et al. (2021a) Event Horizon Telescope Collaboration, Akiyama, K., Algaba, J. C., et al. 2021a, ApJ, 910, L12, doi: 10.3847/2041-8213/abe71d
  * Event Horizon Telescope Collaboration et al. (2021b) —. 2021b, ApJ, 910, L13, doi: 10.3847/2041-8213/abe4de
  * Event Horizon Telescope Collaboration et al. (2021c) —. 2021c, ApJ, 910, L12, doi: 10.3847/2041-8213/abe71d
  * Event Horizon Telescope Collaboration et al. (2022a) Event Horizon Telescope Collaboration, Akiyama, K., Alberdi, A., et al. 2022a, ApJ, 930, L13, doi: 10.3847/2041-8213/ac6675
  * Event Horizon Telescope Collaboration et al. (2022b) —. 2022b, ApJ, 930, L16, doi: 10.3847/2041-8213/ac6672
  * Fishbone & Moncrief (1976) Fishbone, L. G., & Moncrief, V. 1976, ApJ, 207, 962, doi: 10.1086/154565
  * Gammie et al. (2003) Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444, doi: 10.1086/374594
  * Hunter (2007) Hunter, J. D. 2007, Computing in Science and Engineering, 9, 90, doi: 10.1109/MCSE.2007.55
  * Igumenshchev et al. (2003) Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 1042, doi: 10.1086/375769
  * Issaoun et al. (2022) Issaoun, S., Wielgus, M., Jorstad, S., et al. 2022, ApJ, 934, 145, doi: 10.3847/1538-4357/ac7a40
  * Jorstad et al. (2023) Jorstad, S., Wielgus, M., Lico, R., et al. 2023, ApJ, 943, 170, doi: 10.3847/1538-4357/acaea8
  * Kamionkowski & Kovetz (2016) Kamionkowski, M., & Kovetz, E. D. 2016, ARA&A, 54, 227, doi: 10.1146/annurev-astro-081915-023433
  * Koide et al. (1999) Koide, S., Shibata, K., & Kudoh, T. 1999, ApJ, 522, 727, doi: 10.1086/307667
  * Lakshminarayanan & Fleck (2011) Lakshminarayanan, V., & Fleck, A. 2011, Journal of Modern Optics, 58, 1678, doi: 10.1080/09500340.2011.633763
  * Liska et al. (2022) Liska, M. T. P., Chatterjee, K., Issa, D., et al. 2022, ApJS, 263, 26, doi: 10.3847/1538-4365/ac9966
  * Marshall et al. (2018) Marshall, M. D., Avara, M. J., & McKinney, J. C. 2018, MNRAS, 478, 1837, doi: 10.1093/mnras/sty1184
  * Mościbrodzka et al. (2016) Mościbrodzka, M., Falcke, H., & Shiokawa, H. 2016, A&A, 586, A38, doi: 10.1051/0004-6361/201526630
  * Mościbrodzka & Gammie (2018) Mościbrodzka, M., & Gammie, C. F. 2018, MNRAS, 475, 43, doi: 10.1093/mnras/stx3162
  * Narayan et al. (2003) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2003, PASJ, 55, L69, doi: 10.1093/pasj/55.6.L69
  * Narayan et al. (2012) Narayan, R., Sądowski, A., Penna, R. F., & Kulkarni, A. K. 2012, MNRAS, 426, 3241, doi: 10.1111/j.1365-2966.2012.22002.x
  * Narayan & Yi (1995) Narayan, R., & Yi, I. 1995, ApJ, 444, 231, doi: 10.1086/175599
  * Narayan et al. (1995) Narayan, R., Yi, I., & Mahadevan, R. 1995, Nature, 374, 623, doi: 10.1038/374623a0
  * Noll (1976) Noll, R. J. 1976, Journal of the Optical Society of America (1917-1983), 66, 207
  * Oliphant (2007) Oliphant, T. E. 2007, Computing in Science and Engineering, 9, 10, doi: 10.1109/MCSE.2007.58
  * Palumbo et al. (2020a) Palumbo, D. C. M., Wong, G. N., & Prather, B. S. 2020a, ApJ, 894, 156, doi: 10.3847/1538-4357/ab86ac
  * Palumbo et al. (2020b) —. 2020b, ApJ, 894, 156, doi: 10.3847/1538-4357/ab86ac
  * Prather et al. (2021) Prather, B., Wong, G., Dhruv, V., et al. 2021, The Journal of Open Source Software, 6, 3336, doi: 10.21105/joss.03336
  * Reback et al. (2021) Reback, J., jbrockmendel, McKinney, W., et al. 2021, pandas-dev/pandas: Pandas 1.3.2, v1.3.2, Zenodo, doi: 10.5281/zenodo.5203279
  * Rees et al. (1982) Rees, M. J., Begelman, M. C., Blandford, R. D., & Phinney, E. S. 1982, Nature, 295, 17, doi: 10.1038/295017a0
  * Shapiro et al. (1976) Shapiro, S. L., Lightman, A. P., & Eardley, D. M. 1976, ApJ, 204, 187, doi: 10.1086/154162
  * Sądowski et al. (2017) Sądowski, A., Wielgus, M., Narayan, R., et al. 2017, MNRAS, 466, 705, doi: 10.1093/mnras/stw3116
  * van der Walt et al. (2011) van der Walt, S., Colbert, S. C., & Varoquaux, G. 2011, Computing in Science and Engineering, 13, 22, doi: 10.1109/MCSE.2011.37
  * Walker et al. (2018) Walker, R. C., Hardee, P. E., Davies, F. B., Ly, C., & Junor, W. 2018, ApJ, 855, 128, doi: 10.3847/1538-4357/aaafcc
  * Waskom et al. (2020) Waskom, M., Botvinnik, O., Ostblom, J., et al. 2020, mwaskom/seaborn: v0.10.0 (January 2020), v0.10.0, Zenodo, doi: 10.5281/zenodo.3629446
  * Wielgus et al. (2020) Wielgus, M., Akiyama, K., Blackburn, L., et al. 2020, ApJ, 901, 67, doi: 10.3847/1538-4357/abac0d