Constraining the Existence of Axion Clouds in M87* with Closure Trace Analyses

Black holes can amplify incoming bosonic waves via rotational superradiance, inducing bound states of ultralight bosons around them. This phenomenon has the potential to confine the parameter spaces of new bosons. Axions and axion-like particles (ALPs) are candidate beyond-standard-model particles that can form such clouds around supermassive black holes (SMBHs) and impact the polarization signal in a similar fashion to Faraday rotation via axion–photon coupling. Prior efforts have used polarized images from the Event Horizon Telescope (EHT) M87 2017 observations to limit the dimensionless axion–photon coupling to previously unexplored regions. However, with the novel calibration-insensitive quantities, closure traces, and conjugate closure trace products, it is possible to constrain the existence of axion clouds while avoiding the dominant sources of systematic uncertainties, e.g., station gains and polarization leakages. We utilize a simple geometric model for the polarization map of M87* to fit the model parameters with both simulated and real data sets and reach a comparable level of constraint in the accuracy with which an axion cloud may be excluded in M87. Future applications of our approach include subsequent M87* and Sgr A* observations by EHT and the next-generation EHT that are expected to produce stronger constraints across a wider range of axion and ALP masses. Because it does not require imaging, closure trace analyses may be applied to target active galactic nuclei for which imaging is marginal, extending the number of SMBHs from which axion limits may be obtained significantly.


INTRODUCTION
Being among the most luminous astronomical objects in the universe, active galactic nuclei (AGN) are observed and studied across the entire electromagnetic spectrum.For example, in the radio regime, AGN are identified as extragalactic sources of synchrotron emission, which arises from the interaction of relativistic electrons with magnetic fields.The radio emission from AGN is often highly variable and can be used to study the properties of the relativistic outflows that are launched from the vicinity of the SMBH (Boccardi et al. 2017).Furthermore, new technologies in radio observations of AGNs are constantly pushing the limit of the resolution of images of those cosmic giants, revealing their complicated magnetic structures.
Among the multiple pathways in the radio study of AGN, polarization observation of those targets and their central engines plays a crucial role in unveiling the inner structures, especially the intrinsic magnetic fields of these cosmic giants [Johnson et al. 2015, M87 * Paper VII].Fortunately, the advent of very long baseline interferometry (VLBI) has enabled us to overcome the challenge posed by their extremely small angular sizes, allowing for detailed analyses of their vicinity (Thompson et al. 2017).Polarization signals from those massive giants encode essential information about the dynamics and magnetic structures in the surrounding environments.In 2021, the Event Horizon Telescope Collaboration presented the first horizon-resolving polarimetric images of M87* shadow (M87 * Paper VII), enabling unprecedented detailed analysis of its physics.
Besides the study of magnetic structures (M87 * Paper VIII), the acceleration of the bulk flow, and the relativistic heating of jet constituents (Hada et al. 2018), observations of SMBHs offer a unique opportunity to probe beyond-standard model (BSM) physics and place constraints on speculative particles, such as QCD (quantum chromodynamics) axions, which have been proposed as a solution to the strong CP (charge-parity) problem (Peccei & Quinn 1977), and axion like particles (ALPs), which serve as candidates for dark matter (Preskill et al. 1983).
The underlying mechanism is simple as follows.Axions would accumulate around black holes via black hole superradiance, under which traveling light bosons form bound states around SMBHs and grow exponentially as long as superradiance condition, ω < mΩ H , is satisfied.ω is the frequency of the incoming wave, Ω H is the angular frequency of the black hole and m is the azimuthal quantum number.In addition, the most significant superradiant rate occurs when m = 1, the lowest possible energy state that satisfies the condition (Dolan 2007).The superradiant rate is maximized when the axion's reduced Compton wavelength λ c is comparable to the black hole size (Arvanitaki et al. 2015).In the presence of a background axion field, the modification of a photon's equation of motion due to the axion-photon interaction will introduce a periodic oscillation to the EV-PAs (electric vector position angles) of linearly polarized photons (Carroll et al. 1990;Carroll & Field 1991;Harari & Sikivie 1992).Theoretically, modelling the variations of EVPAs around the black hole (BH) shadow of an SMBH should indicate novel constraints on light axions aside from previous results1 .
An experiment was then proposed when Chen et al. (2020) addressed that this axion-induced birefringence effect can cause the EVPAs on the photon ring of an SMBH to vary with time and position.In the condition of photon propagating in a background axion field, axion-photon interaction modifies the Lagrangian as (1) in which g aγ is the axion-photon coupling constant, F µν is the electromagnetic tensor and F µν is its dual tensor.A naive deduction would be that this axion-induced effect would function very similarly to the normal Faraday effect.Successive derivations with the assumption that the photon frequency is much larger than axion's reduced Compton wavelength λ c reached a simple expression for the EVPA rotation of a linearly polarized photon: (2) The rotation only depends on the difference between the emitting and observing axion field values.With further appropriate assumptions and simplifications that describe the EVPA variation in Boyer-Lindquist coordinates of the black hole, the final form can be parametrized like in which a polar coordinate (r, ϕ) centered at the black hole shadow is used.For a detailed derivation of the previous equations, see Chen et al. (2022a).
Recently, a practical attempt by the same group (Chen et al. 2022b) utilize the polarimetric information of M87* (M87 * Paper VII) to reach a previously unexplored parameter space of light axions.Specifically, they implement constraints on the joint distribution of the axion mass m a and the dimensionless axion-photon coupling constant c aγ ≡ 2πg aγ f a .However, one caveat is that the parametrized form of EVPA variation relies on the assumption of a saturated axion cloud.While it is impossible to know if the axion cloud is saturated given the unknown recent astrophysical history of M87, we adopt this assumption henceforth (see Chen et al. (2022a) for a full discussion of this issue).
Despite its many advantages, polarized radio interferometry, which is used to generate the polarimetric images of M87*, suffers from multiple systematic uncertainties, including station gains that modulate both the amplitude and the phase of the observed signal, and polarimetric leakages from one channel to another, e.g., from right-hand circular polariztion to left-hand circular polarization (M87 * Paper III).The former one can be circumvented by "closure" quantities (closure phase and closure amplitude) on multiple sites that are immune to station-based effect because station-based gain errors are cancelled out (Jennison 1958;Twiss et al. 1960;Readhead et al. 1980), while the latter one is calibrated with suboptimal assumptions (M87 * Paper VII).Luckily, Broderick & Pesce (2020) introduced a novel calibration-insensitive closure quantity: closure trace (CT), defined on station quadrangles, that is insensitive to station gains and polarimetric leakages (D-terms) and comprises all calibration-independent information in the visibility data; and a subsequent quantity, the conjugate closure trace product (CCTP), that can serve as a robust indicator of polarization sources.
Therefore, based on the information of axion-photon coupling and CT analysis, we present a new method to constrain light axions and demonstrate that it can reach a roughly similar level of constraint as Chen et al. (2022b).Our method, being insensitive to a much larger class of errors than ameliorated by gains and D-terms, circumvents the time-consuming direct calibration process to the final polarized images of M87* and is complete in retrieving polarimetric information.This novel method can be further applied to future EHT and nextgeneration EHT (ngEHT) observations (Johnson et al. 2023), for which more quadrangles will be available, to reach stronger constraints on axions, ALPs and other ultralight bosons.
The structure of this paper is as follows.In Section 2, we introduce our model of EVPA variations, explore and discuss its feasibility.In Section 3, we present our Monte Carlo Markov Chain (MCMC) fits for the model and its implications.Finally, in Section 4 we make conclusions.

A SIMPLE GEOMETRIC MODEL FOR THE LINEAR POLARIZATION MAPS OF M87*
The most striking feature of the EHT images of M87* (Figure 1) is the clear presence of a ring, presumably encircling the black hole shadow.Upon this ring lies a  "twisty" linear polarization pattern with the EVPA exhibiting a nearly linear dependence on azimuth.Given these gross features, we generate a simple geometric model for the polarized images of M87* with which to explore the constraints made possible using the measured CCTPs.

Image Domain Model
Besides the outstanding feature of an annulus in the EHT images of M87*, we clearly notice the position angle of the maximum Stokes I map lies in the southern direction, the position angle of the maximum linear polarization lies in the SW direction and the EVPA position angle around the shadow appears to deviate from an overall radial EVPA pattern (Figure 1).Based on this information and putting time-dependent EVPA variations aside, we begin by creating movie objects to model the flux and the polarization map, that can be processed by eht-imaging (Chael et al. 2018).The parameters include the diameter d of the annular emission region in µas, annulus width w, linear polarization fraction m, the number of periods n cyc per 24 h, the position angle Φ I of the maximum brightness in Stokes I map and the position angle Φ pol of the maximum polarized flux, where the latter two angles start from 0 in the south direction on the image and increase clockwise, i.e., and (5) Given an EVPA, θ(t, ϕ), we can then construct the remaining Stokes parameters via We convolve the polarimetric images by a Gaussian filter and then generate movies that reflect the time evolution of EVPAs that results from axion-photon coupling.Combining the most prominent features of the polarized images of M87* and the parametrization of EVPA variation by Chen et al. (2022b), we construct a simple geometric model of the EVPA oscillation by introducing an additional time-dependent position angle θ on top of a fundamental time-independent radial EVPA pattern (Figure 2) around the shadow.
In this model, θ 1 is the oscillation amplitude and is proportional to c aγ and ω is the axion angular frequency proportional to m a , later substituted by n cyc , as in Table 1.The rest of the parameters are defined as follows: r ring = √ 27r g is the radius of the shadow (assuming a=0) of M87*, θ 0 is the initial value of θ, ϕ is the azimuthal angle on the annulus and δ describes an arbitrary initial phase.The inclination angle of 17 • between the spin direction of M87* and the sky plane is introduced by Chen et al. (2022b) to account for the time delay of emissions.θ 1 is independent of the azimuthal angle because the azimuthal variation is rather marginal, as seen in Chen et al. (2022b).Figure 2 shown below is the fundamental radial EVPA plot we started from, while in Figure 3 we added θ and clipped a total of eight frames from the EVPA movie of our model with morphologically estimated parameters of 2017 April 11 to fully cover its evolution.This simple geometric model captures the characteristics of time-dependent EVPA variation induced by axion-photon coupling appearing as propagating waves around the shadow.
Chen et al. (2022b) extended the "forbidden" region of axions on the c aγ -m a parameter space by directly fitting the differential EVPA on the polarimetric images across 2017 April 5, 6, 10, 11, highlighting the possibility of exploring new physics with the unprecedented EHT images on SMBHs (although assuming the existence of a saturated, superradiantly generated axion cloud around M87*).However, it is important to note that this approach neglect the systematic uncertainties within the EVPA, like station gains and D-terms, which will be discussed below, associated with calibrations to the final polarimetric images.Our new method, on the other hand, takes advantage of novel calibration-insensitive quantities for radio astronomy, to "observe" the timedependent movies generated by this model, and reach an approximately comparable level of constraints on light axions, without the need for full polarimetric image reconstruction, with the attendant systematic uncertainties.

Constructing EHT Observables
The primary EHT data products are visibilities, corresponding to the Fourier transform of the Stokes I, Q, U, and V brightness maps.We obtain these using the eht-imaging package via the observe function, generating visibilities in the more common circular basis, RR, LL, RL, LR (see M87 * Paper VII), commonly collected into a coherence matrix, where the indicies ij refer to the two stations defining a particular baseline.All baselines relevant for EHT are generated at the same time.However, the observed V ij are degraded by aforementioned station gains and  polarimetric leakages (D-terms), which means the true coherence matrix are modified as where G and D contain the gain terms and leakage terms for each station respectively, and † is the complex conjugate transpose.Polarimetric imaging requires full reconstructions of the time-dependent station characteristics (i.e., the G and D).However, the immunity/insensitivity of the CT and CCTP to G and D present an opportunity to avoid this process.
The CTs are generated from combinations of V ij on quadrangles, sets of four stations whose baselines form a closed circuit (Broderick & Pesce 2020).
These are invariant to all linear, station-based corruptions of the coherence matrices, including atmospheric phase delays, receiver gains, and polarization leakage.
The full complement of T ijkl contain all of the residual information in the coherence matrices apart from the standard calibration quantities, and are a superset of the more familiar closure phases and closure amplitudes (Thompson et al. 2017).
While the T ijkl are generally sensitive to polarized and unpolarized image substructure, a specific combination of them, the CCTP defined by are sensitive only to polarization substructure; in the absence of polarization substructure, the CCTP is identically unity.Therefore, the CCTPs are valuable probes of polarization in much the same way that closure phases are sensitive to asymmetry.For M87*, the only quadrangle exhibiting CCTPs with the most significant departures from unity, and therefore the most significant evidence for polarized structure, was that combining APEX, ALMA, LMT, and SMT (see Fig 13 of M87 * Paper VII).One potential explanation is that ALMA and LMT have very larger apertures, and are thus the most sensitive stations.Another possible reason is that the ordered polarization pattern indicates that the most power is on scales of a few Gλ, which corresponds to the Chile(ALMA, APEX)-LMT and Chile-SMT baselines (usually below 4 Gλ).Therefore, we construct the CCTPs from the complex visibilities generated by eht-imaging on this quadrangle for exploration and comparison with the 2017 EHT data on M87*.

Sensitivity of CCTPs to Axion Signal
A fiducial set of parameters, listed in Table 1, is chosen to qualitatively match the general morphology of the polarized images on 2017 April 11. Figure 4 shows the CCTP calculated from the fiducial parameters plotted against read CCTP of 2017 Apr 11.Specifically, we choose the position angle of the maximum brightness in Stokes I map to lie in the south, the position of the maximum polarized flux to lie in the southwest, and fixed polarization fraction of m = 25%.The fiducial period is selected so that α ≡ r g /λ c = 0.4, consistent with that shown in Figure 1 of Chen et al. (2020), which corresponds to a period of roughly 131 h, or n cyc = 0.183 cycles in a 24 h period.The CCTPs produced from this model match those observed on 2017 Apr 11 quantitatively.
Figure 5 shows the impact of increasing the variable component of the EVPA pitch angle, θ 1 , associated with the axion cloud density.For reference the 2017 Apr 11 high band CCTP measurements are shown in light gray, which indicate that CCTP-based constraints on θ 1 of the order of 10 • should be possible.In contrast, variations in θ 0 (not shown) have no impact on the CCTPs; changing θ 0 induces a total rotation in the EVPA across the image to which the closure traces are invariant.
Due to the long fiducial period, modifications of the axion phase offset δ can have a large impact on the sensitivity of CCTPs to a putative axion signal.Shown in Figure 6 are two models at very different locations in their respective oscillations.The very different sensitivity to θ 1 is easy to understand given the form of θ given in Equation 7. When δ = 0 • , for small times θ evolves only slowly.However, when δ = ±90 • , the evolution is maximized.For small times θ depends linearly on time, and thus the variation in the CCTPs is a combination of both θ 1 and the temporal evolution throughout a single night.Unsurprisingly, this behavior is strongly dependent on n cyc as shown in Figure 7. Non-zero θ 1 would induce peaks as substructures on the overall CCTP curves, with the number of peaks proportional to n cyc .In addition, for shorter periods, higher n cyc , the magnitude of the variation in the CCTPs driven by θ 1 is generally increased, maximizing when n cyc ≳ 1.
The dependence on δ and n cyc motivates a number of practical concerns for constraining the existence of axion clouds in EHT targets.For M87*, the large SMBH mass requires a large period, and thus small n cyc -constraints will benefit from coherently combining CCTP measurements across many nights.This is further justified by the long intrinsic timescale for astrophysical changes in the emission region.
However, for Sgr A*, the much lower mass requires a smaller axion mass to drive the superradiant instability, and thus shorter orbital period and larger n cyc .Simply rescaling by the SMBH mass and keeping other physical parameters the same, our fiducial M87* model produce n cyc = 270.29 for Sgr A*, comparing to n cyc = 0.183 for M87*.This is shown in Figure 7, from which we conclude that even a single observation night may strongly constrain the existence of an axion cloud in the Galactic center.
3. IMPLICATIONS OF M87* CCTP DATA Armed with an expectation that CCTPs may be a sensitive discriminator among axion models in M87* and Sgr A*, we performed a number of fits to simulated and real data sets.We begin with a summary of the fitting procedure, followed by demonstration on simulated data sets, and ending with an application to the 2017 EHT M87* data.

Fit Procedure
In all cases, we make use of the Monte Carlo Markov Chain sampler emcee to explore the posterior space of the model parameters (Foreman-Mackey et al. 2013).
To fit the geometric model with the CCTPs of the entire 4 days of observations, both real and simulated data are grouped into two sets, each contains two consecutive days of observations, i.e Apr 5, 6 and Apr 10, 11, allowing the variability of M87* in a timescale of a week.In addition, a fit directly to Apr 5, 6, 10, 11 is also performed.Beginning with the fiducial model parameters in Table 1, we calculate the corresponding CCTP values and then the log-probability by comparing the predicted and observed CCTP values, along with the corresponding errors, and summing over all datapoints, and update them within 200000 samples for a total of 40 walkers in parallel, i.e., a total of 8 million samples.Triangle plots are produced from chains (Foreman-Mackey 2016), and plots for the joint posterior distributions with axes converted to log(c aγ ) and α in alignment with Chen et al. (2022a) are generated as the axion limits derived from CCTP analyses.In particular, we transform the θ 1 -n cyc plot into the log(c aγ )-m a plot with the conversion factors discussed in the appendix.b Measured clockwise relative to radial direction.
c Number of full periods in 24 h, proportional to ω, the axion angular frequency.
d Defined such that θ = 0 • at 0 UT on 2017 Apr 5 or on 2017 Apr 10, depending on which dataset we are fitting.
e The ranges of parameters used for Figure 5, Figure 6 and Figure 7.

Simulated Data Sets
In order to explore the impact of large θ 1 and the significant errorbars that came with the observation data, a total of 4 simulated data sets were created based on the actual mass of M87* and two choices of predetermined θ 1 .Posterior distribution plots for the simulated data sets are shown in Figure 8.
The null case (θ 1 = 0) shows a clear excluded region of axion parameters in the upper right of the posterior plot, resembling the feature in Figure 4 of Chen et al. (2022a), although our boundary is weaker.However, the case with non-trivial θ 1 = 60 • indicates a significant clump of high probabilities located roughly between log(c aγ ) = 1 and log(c aγ ) = 2, where the actual parameters fall in.The "truth" is indeed narrowly constrained.Therefore, we believe the fitting works as expected since we clearly see anticipated behaviors in Figure 8, which recovers reasonable values of M87* mass and θ 1 , and uncertainties.Overall, considering the fact that simulated data sets were constructing based on the method of control variable, patterns deviating from the null case posterior plot can indicate a possible detection of light axions that can trigger EVPA variations with magnitude above the level which can be explained by the errorbars.

2017 EHT M87* Campaign
As stated above, fits to EHT data were performed on the 2017 Apr 5, 6 and Apr 10, 11 data sets separately.In addition, a fit to the entire data set (2017 Apr 5, 6, 10, 11) was also executed.The resulting 2σ band of fits is overplotted on the EHT data in Figure 9. Posterior distribution plots, including 2σ and 3σ regions, are shown in Figure 10.The posterior plot for Apr 10, 11 clearly matches a non-detection of light axions at even 2σ and is consistent with the null case posterior plot in Figure 8.Although it may seem that on Apr 5, 6, there is a phantom of detection at 2σ level, there is no evidence for any detection at 3σ.This interpretation, is confirmed by the fit to the entire data set.That is, we see no evidence for an axion cloud from the 2017 M87* EHT campaign.

Discussion
In summary, fits for simulated data sets properly recover the undelying predetermined parameters, indicating the feasibility of using CCTP measurements to probe for light axions.For EHT 2017 data, at 3σ level, we conclude that there is no detection of light axions and therefore we arrive at an upper limit on the putative axion-photon coupling.
It is essential to carefully define what a "detection" and "non-detection" of axions means in this context given the substantial astrophysical uncertainties.For a non-detection it suffices to conclude that within the current accuracy that the null hypothesis, no axion cloud, is consistent with the observation.Consequently, we exclude the axion parameter spaces above an observationally-dependent bound, above which any light axion is expected to result in CCTPs that are inconsistent with the EHT data.Conversely, the EVPA variations due to lower density axion clouds are consistent within the formal uncertainties.
We note that to claim a "detection" of an axion cloud, the elimination of the null hypothesis is necessary but insufficient.In addition, all potential astrophysical origins for the EVPA oscillations (e.g., evolving magnetic field structures) must be excluded, something we have not explored here.Therefore, at present we can only produce upper-limits on the light axion cloud density associated with non-detections, and therefore the coupling (subject to assumptions about the history of the axion cloud).
Figure 11 summarizes our constraint to light axion parameters derived via the CCTP analysis on the M87* 2017 linear polarization data.The bound on the dimensionless coupling constant becomes weaker for smaller axion masses, primarily due to a smaller value of the radial wave function of the axion cloud and a smaller single-day axion field variation resulting from the longer oscillation period (Chen et al. 2022a).Despite presenting a modestly weaker limit than the preceding one, our result highlights the feasibility of CCTP analyses for constraining the existence of axion clouds.More importantly, the image-based analysis of Chen et al. (2022a) is inevitably subject to significant underlying systematic uncertainties originating from the estimation of station gains and D-terms, to which the CCTP-based analysis is insensitive.Additionally, the use of CCTP analyses on the polarization data of other SMBHs, for which full polarimetric imaging may be infeasible or impractical, holds potential for cross-validation of the existence of light axions; more is discussed in Section 4.
4. CONCLUSION Photons propagating in a light axion background around black holes would experience a variable, spatially-dependent birefringence effect, caused by axion-photon coupling, which may be detectable via polarimetric studies of the horizon-scale emission.Imagebased studies on constraining the parameters of light axion have been done by modeling the temporal EVPA variations around the supermassive black hole M87*.Unlike previous efforts, CCTPs provide a novel nonimaging method for constraining the existence of light axions.The direct analysis of closure quantities, i.e., closure traces and the CCTPs generated from them, provide two significant advantages.First, they are intrinsically immune to station-based polarimetric systematics, including atmospheric phase delays and absorption, and polarimetric leakage.Second, they are a non-imaging quantity, and therefore do not require image reconstructions and their attendant assumptions.
Application of CCTP-fitting to EHT observation of M87* finds no detection of axions, resulting in an upper limit on the coupling constant of log 10 (c aγ ) ≲ 1 for axions with mass 2.2 × 10 −21 eV < m a < 9.9 × 10 −21 eV.This conclusion is, however, based on the assumption that a saturated axion cloud exists around M87*, generated by the superradiant instability and limited only by  axion-axion annihilation (bosenova).That is, we ignore the complicated astrophysical history of this target.
A significant current uncertainty in the M87* analysis is the unknown phase of the axion cloud (δ in Equation 7).Future observations of M87* extending over multiple weeks may retire this uncertainty, though potentially at the expense of introducing additional astrophysical variability.
Applying this technique to Sgr A* would probe a very different range of axion masses, and therefore temporal scales.The expected CCTP signal is qualitatively different for Sgr A* than M87*, exhibiting rapid, periodic fluctuations throughout a single night that may be excluded by existing and future EHT observations; tests based on Sgr A* are forthcoming.Because the genera-tion of CCTPs does not require full imaging, our method may be applied to marginally resolved and possibly unresolved targets, including many AGN targets beyond M87* and Sgr A*, and thereby expanding the range of axion masses probed.
Already, following the 2017 EHT campaign, additional sites have been added to the array, increasing the number of available quadrangles significantly.The improved baseline coverage, spatial and temporal resolution of the next-generation EHT (Johnson et al. 2023) will enable more detailed geometric modeling and more precise EVPA measurements, resulting in stronger constraints on light axions.The magnitude of the improvement can be roughly estimated from the increased volume of observations (4 d to 16 d), increased bandwidth (4 GHz to 64 GHz), increased number of quadrangles (20 stations will produce 170 independent quadrangles, in comparison to the single quadrangle considered here), and by analyzing both the amplitude and phase of the CCTPs on these quadrangles.Assuming that each quadrangle produces CCTPs with similar statistical power, the limit on the axion-photon coupling constant, log 10 (c aγ ) could be improved by two orders of magnitude.Verifying this with simulated data sets is a natural future project.Software: eht-imaging, emcee APPENDIX A. APPENDIX While natural to characterize the observational impact of the axion cloud on the evolving EVPA in terms of an amplitude (θ 1 ) and number of cycles (n cyc ), these quantites are directly related to the underlying axion properties.Here we collect the various conversion relations between the two sets of quantities.
The oscillation frequency of the axion field, ω a , is proportional to the axion mass, m a (Plascencia & Urbano 2018).Incorporating fundamental constants, the num-ber of cycles in a 24 h period is then related to m a via,

Figure 1 .Figure 2 .
Figure 1.Time evolution of the fiducial average polarimetric images of M87* from EHT 2017 Campaign from Apr 5 to Apr 11.Total intensity is shown in grayscale.EVPA is shown in colorful ticks, while the length indicates linear polarization intensity magnitude, and color indicates fractional linear polarization.(Reproduced from M87 * Paper VII.)

Figure 3 .
Figure 3. Eight snapshots with phases between 0 and 2π captured from the time-dependent model, approximating the stationary model in Figure 2 with a time-dependent axion cloud.The phase of each snapshot is labeled at the bottom left.The position angles of maximal Stokes I and polarization flux are chosen to match the actual image of 2017 Apr 11.The intensity color range and the length of the polarization ticks are the same as in Figure 2

Figure 4 .
Figure 4. HI+LO band real CCTP measurements for 2017 Apr 11, with CCTP curves from the fiducial set of parameters.HI and LO refer to the high frequency band and low frequency band centered at 229.1 GHz and 227.1 GHz, respectively.LO band data are intentionally shifted to the right by 0.02 hrs for better readability.The errorbars in the CCTP measurements are estimated using Monte Carlo sampling, are purely due to the underlying thermal errors in the visibilities, and are well-approximated as Gaussian errors.

Figure 5 .
Figure 5. Change in θ1, the variable EVPA pitch angle proportional to the axion-photon coupling constant.The first line in the legend indicates the target of the simulated CCTP.θ0 = 30 • is arbitrary as changes of θ0, the mean EVPA pitch angle, have no impact on the shape of the curve.From red to purple, the spacing of θ1 between two consecutive colors is 2 • .For reference, the HI band real CCTP measurements on Apr 11 are plotted in light grey.

Figure 6 .
Figure 6.Change in δ, the axion cloud phase offset.The first line in the legend indicates the target of the simulated CCTP and the parameter to be varied.From red to purple, the spacing of θ1 between two consecutive colors is 2 • .For reference, the HI band real CCTP measurements on Apr 11 are plotted in light grey.

Figure 7 .
Figure 7. Change in ncyc, or the number of cycles which is proportional to axion mass.The first line in the legend indicates the target of the simulated CCTP and the parameter to be varied.From red to purple, the spacing of θ1 between two consecutive colors is 2 • .For reference, the HI band real CCTP measurements on Apr 11 are plotted in light grey.

Figure 8 .
Figure 8. Posterior distribution plots for simulated data.The difference in date and predetermined θ1 is labeled at the bottom left.Constraints from Chen et al. (2022a) are overplotted as red lines, while the blue lines are our estimates of excluded regions of axion parameters.Dashed lines are 3σ constraints and solid lines are 2σ constraints.The negative numbers on the colorbar imply a logscale of the probability in the corresponding posterior plot.The actual parameters used to generate those simulated cases are labeled as black dots on the posterior distribution plots or indicated by black arrows if they fall way below the lower limit.

Figure 9 .
Figure 9. 2σ band plots over the real data, depicted in light blue, estimated from the MCMC fitting procedures on all 4 days.The central blue lines indicate the average fitting result.

Figure 10 .
Figure 10.Posterior distribution plots for real data.Constraints from Chen et al. (2022a) are overplotted as red lines, while the blue lines are our estimates of excluded regions of axion parameters.Dashed lines are 3σ constraints and solid lines are 2σ constraints.The negative numbers on the colorbar imply a logscale of the probability in the corresponding posterior plot.
We thank Asimina Arvanitaki and Yifan Chen for their valuable discussions.This work was supported in part by Perimeter Institute for Theoretical Physics.Research at Perimeter Institute is supported by the Government of Canada through the Department of Innovation, Science and Economic Development Canada and by the Province of Ontario through the Ministry of Economic Development, Job Creation and Trade.A.E.B. thanks the Delaney Family for their generous financial support via the Delaney Family John A. Wheeler Chair at Perimeter Institute.A.E.B. receives additional financial support from the Natural Sciences and Engineering Research Council of Canada through a Discovery Grant.