Probing Accretion Turbulence in the Galactic Center with EHT Polarimetry

Magnetic fields grown by instabilities driven by differential rotation are believed to be essential to accretion onto black holes. These instabilities saturate in a turbulent state; therefore, the spatial and temporal variability in the horizon-resolving images of Sagittarius A* (Sgr A*) will be able to empirically assess this critical aspect of accretion theory. However, interstellar scattering blurs high-frequency radio images from the Galactic center and introduces spurious small-scale structures, complicating the interpretation of spatial fluctuations in the image. We explore the impact of interstellar scattering on the polarized images of Sgr A* and demonstrate that for credible physical parameters, the intervening scattering is non-birefringent. Therefore, we construct a scattering mitigation scheme that exploits horizon-resolving polarized millimeter/submillimeter VLBI observations to generate statistical measures of the intrinsic spatial fluctuations and therefore the underlying accretion flow turbulence. An optimal polarization basis is identified, corresponding to measurements of the fluctuations in magnetic field orientation in three dimensions. We validate our mitigation scheme using simulated data sets and find that current and future ground-based experiments will readily be able to accurately measure the image-fluctuation power spectrum.


INTRODUCTION
Black holes have been implicated as the engines of active galactic nuclei (AGN) and X-ray binaries. Within these objects, both their extreme luminosities and growth rate are presumably due to the interaction with the accretion of nearby matter. This occurs via accretion disks, through which material orbits, cools and falls inward toward the central object. Accretion flows are generic features in astronomical systems, from the formation of planets to the powering of AGN, and thus understanding the processes by which they operate informs astrophysics broadly.
The supermassive black hole at the Galactic Center, Sagittarius A* (Sgr A*), offers us a laboratory in which to study accretion flows in detail. Located 8 kpc from the Earth, with a mass of 4.3 × 10 6 M (Boehle et al. 2016;Gillessen et al. 2009;GRAVITY Collaboration et al. 2018), Sgr A* has now been resolved on event horizon scales with the Event Horizon Telescope (EHT Event Horizon Telescope Collaboration Corresponding author: Chunchong Ni chunchong.ni@uwaterloo.ca et al. 2022a,b,c,d,e,f). These observations present an unprecedented opportunity to probe the nature and characteristics of the hot plasma orbiting Sgr A* under the extreme conditions near the horizon.
The EHT is a global array of millimeter and submillimeter telescopes that achieves resolutions of 20 µas at a wavelength of 1.3 mm (230 GHz) via very-long baseline interferometry (VLBI). This resolution is sufficient to resolve the event horizons of Sgr A* and M87*, silhouetted against the emission from the surrounding hot plasma. In comparison, the typical angular size of the shadow for Sgr A* and M87 is around 50µas. Therefore, it has now become possible to probe accretion physics on scales comparable to those relevant for MHD turbulence.
Four observing campaigns have been completed by the EHT, in April 2017, and the first M87 and Sgr A* EHT results have been published (Event Horizon Telescope Collaboration et al. 2019a,b,c,d,e,f, 2022a. The full complement of Stokes parameters were measured at 230 GHz. Future development of the array will include the ability to observe at 345 GHz. Sgr A* presents a natural target for studies of the role played by MHD turbulence in black hole accretion be-cause of its short timescale and lack of an obvious relativistic jet. However, interpreting the small-scale brightness fluctuations, presumably associated with MHD turbulence within the accretion flow, is complicated by the interstellar scattering observed toward the Galactic center (Lo et al. 1998;Frail et al. 1994;Lazio & Cordes 1998a). This scattering is believed to be a result of variations in the electron density along the line of sight (Goldreich & Sridhar 1995;Lazio & Cordes 1998b;Cordes & Lazio 2002;Rickett 1990). Typically, the origin of the scattering is abstracted to a thin scattering screen, for which detailed models exist Johnson et al. 2018;Issaoun et al. 2021;Cho et al. 2022). For Sgr A*, two aspects of the scattering are of interest, corresponding to different regimes: diffractive and refractive scattering (Narayan 1992;Johnson & Gwinn 2015).
The diffractive scattering is the consequence of the combined effect of small-scale fluctuations in the interstellar electron density, whose impact is to blur the image with a nearly-Gaussian kernel (Johnson & Narayan 2016;Issaoun et al. 2021). This angular broadening is formally reversible, i.e., images of Sgr A* may be effectively "deblurred" by applying the appropriate multiplicative correction in the Fourier (visibility) domain (Fish et al. 2014;Lu et al. 2018;. The impact of refractive scattering is more subtle. Associated with the large-scale fluctuations of the interstellar electron density, refraction induces coherent and variable substructures in the image ). These additional variations in the image are extrinsic to the source, and indicative of the interstellar scattering screen. Unlike diffractive scattering, it is not formally invertible, and may not be simply removed during image generation. In principle, it may be modeled, leveraging the modestly different timescales between the refractively induced extrinsic substructure and the intrinsic brightness fluctuations induced by MHD turbulence.
In this paper, we demonstrate that the action of the scattering screen is expected to be independent of polarization. Based on this, we develop a new scattering mitigation scheme that exploits this non-birefringence of the scattering screen. We demonstrate that for the angular scales accessible to the EHT and ngEHT, it is possible to effectively eliminate the impact of interstellar scattering on the estimators of the intrinsic structural polarimetric fluctuations, and thus probe MHD turbulence instrinsic to the near-horizon emission region directly. While a full spatiotemporal characterization of the turbulence is highly desirable (see, e.g., Georgiev et al. 2022), we focus on mitigating the spatial distortions resulting from scattering here, leaving the construction of the temporal component of the power spectrum for future work.
In Section 2 we review scattering in the thin-screen approximation, assess the impact on polarized emission, and demonstrate that scattering can be implemented as a nonbirefringent, tensor convolution that may be inverted. In Section 3, we construct toy models that mimic the gross properties of Sgr A*, and test the feasibility of scattering mitigation. In Section 4, we apply our scheme to a representative simulation from the existing EHT general relativistic magnetohydrodynamic (GRMHD) simulation library, and confirm we are able to extract intrinsic information about the structural variability in spite of the intervening scattering. Finally, we conclude in Section 5.

SCATTERING AND OBSERVATION OF POLARIZED LIGHT
We begin with a summary of the action of an intervening scattering screen upon the emission from a compact source observed via a local interferometer. This is appropriate, e.g., for observations of Sgr A* by the EHT and ngEHT. We will follow presentation of Johnson & Gwinn (2015) where possible and refer the reader there for a detailed description.

Scattering and the Visibility Function
The primary observable quantity in interferometric radio observations, like those made by the EHT and ngEHT, is the "visibility", V (b), constructed by crosscorrelating signals at antennae separated by a projected baseline b. This quantity is directly given by the Fourier transform of the intensity map, i.e., where I(x) is the intensity map projected at the source distance, x is an angular location on the sky, and λ is the observing wavelength (see, e.g., Thompson et al. 2001). As a consequence, the V (b) encodes the degree of source structure on an angular scale of λ/|b|, oriented along the direction of b.
Scattering is frequently modeled in the thin-screen approximation (Bower et al. 2014). The physical picture is presented in Figure 1, which shows the relative position of the emitting source, an intervening thin screen, and the observer on Earth. Thick scattering screens include additional complication, and may be required toward Sgr A* (for example, see Pen & Levin (2014)). Nevertheless, in many cases, these extended scattering regions may be abstracted to a sequence of thin screens. Thus, we will focus on the latter.
The impact of scattering in the thin-screen limit, is to impart a random phase shift at the screen, φ(x). That is, the observed visibilities are where the Fresnel radius, is the characteristic radius at the observer on which the spherical nature of the approaching radio wave become important, and provides a useful scale for scattering phenomena. In Equation 2, we have introduced V obs (b) for the visibility that is observed after scattering and V int (b) for the visibility that would have been observed in the absence of scattering. It is, fundamentally, V int (b) that is of interest to studies of the compact astronomical sources.
In the ensemble average regime, obtained after averaging V obs (b) over many realizations of the scattering screen, the observed visibility is given by where . . . ea denotes ensemble averaging (Johnson & Gwinn 2015). The D φ (b) is the structure function of the phase fluctuations on the scattering screen, defined in the normal way: where for a statistically isotropic screen, like that we assume here, the absolute position x does not matter. Equation 4 is the well-known diffractive limit, in which scattering imparts only a multiplicative correction to the appropriately averaged intrinsic visibilities, and in which images may be deblurred in the normal sense (Johnson & Gwinn 2015). Analogous observable quantities can be constructed for polarized emission, and we do so for the Stokes maps, , where here Stokes V refers to the excess right-handed circular polarization, not the visibility 1 . In general, the phase shifts may depend on the particular polarization under consideration. In practice, for credible values of the magnetic field strength in the interstellar medium, the scattering screen is non-birefringent.

(Non-)Birefringence of Scattering in the ISM
The degree to which we may assume that an intervening scattering screen is non-birefringent depends on the magnitudes of two closely related quantities, the angular deflections experienced by radio waves passing through the screen (refraction), and the phase shifts imparted on those radio waves during their passage (dispersion). Here, we show that for models of interstellar scattering in which both are due to turbulent fluctuations in the magnetized interstellar plasma, both are sufficiently small that we may treat scattering as independent of the polarization of the radio waves under consideration. We will begin by analyzing the properties of the scattering in the unmagnetized limit, and thus produce estimates for the unpolarized case, followed by an analysis of a weakly Figure 1. Schematic illustration of the interstellar scattering at a thin screen that imparts random phase fluctuations. The paths of two rays, differing in polarization, are shown, indicating the slight difference in the refractive deflection angles θ±. The source-screen (R) and screen-observer (D) distances are indicated. magnetized screen to estimate the disparate impact on different radio wave polarizations. Before starting, we define the following dimensional variables, X and Y, following the conventions in plasma physics, where ω P = 4πn e e 2 /m e 1/2 is the plasma frequency, n e is the free electron density, e the electron charge, and m e the electron mass, ω B = eB/m e c is the cyclotron frequency, and z is the line-of-sight direction. X and Y are proportional to the plasma density and the magnetic field strength at given frequency, ω. Detailed computations for screens composed of turbulent magnetized plasmas reach similar conclusions and are presented in Appendix A.

Refraction in an Unmagnetized Screen
In the absence of a magnetic field, the long-wavelength dispersion relation for electromagnetic waves in an electron-ion plasma is, The associated equations of motion for radio wave are then given by Hamilton's equations obtained by setting H(x, k) = ω(x, k), where k is the wavevector and related to the photon momentum by p = k. This equation of motion describes how electromagnetic waves refract when travelling through the scattering screen. In the weakdeflection limit, the perpendicular momentum gained after propagating through the scattering screen is, where we have made use of the approximation that the line-of-site velocity of the wave is c. Thus, the deflection angle is approximately the ratio between k r and k z , where we used that k z ≈ ω/c. The typical size of scattered compact sources, broadened by the diffractive scattering, places a constraint on the magnitude of θ 0 ≈ θ diff , and therefore, the lineof-sight integrated transverse gradients of the fluctuating electron density within the screen. It is against this value that we will normalize the impact of non-zero magnetic fields, and thus the degree of birefringence for a weakly magnetized scattering screen.

Refraction in a Weakly Magnetized Screen
In the presence of a weak magnetic field, the dispersion relation of electromagnetic waves travelling through a plasma is slightly modified, becoming where there are now two propagating modes, ordinary and extraordinary, signified by the ± sign. In the quasitransverse limit (propagation along the magnetic field), the two polarization modes are nearly circular.
We again obtain the equations of motion from Hamilton's equations, though in this instance they differ for the two polarization modes, The deflection angle after integration through the screen is, therefore, where we have defined Y previously. The mean deflection of the two modes is just that associated with the unmagnetized plasmā The differential deflection, and thus the disparity in the impact of the weakly magnetized scattering screen on the two polarizations, is given by Upon averaging over a random magnetic field orientation, and thus sign of Y , this would vanish. However, the variance of δθ, and thus its typical value, does not vanish, where σ 2 Y ∼ Y 2 is the variance in Y , associated with the magnitude of the magnetic field fluctuations. Thus, typically, the difference in the deflection of the two polarization modes, is reduced by a factor of Y 2 , which at 1.3 mm is small for credible interstellar magnetic field strengths ( 1 mG)

Phase Shift Induced by a Weakly Magnetized Screen
The typical phase fluctuations imparted by the scattering screen are related to, but distinct from, the refraction. Given the dispersion relation in Equation 12, the phase shift for the transverse electromagnetic modes grow as, and thus the difference between the phases of the the polarization modes is, This expression suggests that the differentially accumulated phase and typical differential deflection angle differs by a factor of L/λ, where L is the typical correlation length within the plasma. This suspicion is born out by detailed calculations for a variety of magnetic field and electron density fluctuation spectra presented in Appendix A. Similar results are presented in the past literature (e.g. see Macquart & Melrose (2000)).
In particular, for power-law electron fluctuation spectra, we show that L may be associated with the inner scale, r in , in Psaltis et al. (2018) and Johnson et al. (2018). Therefore, the typical differential deflections may be related to the typical differential phase fluctuations induced by a thin scattering screen by, Given a typical value for θ diff ≈ 10 µas, an interstellar magnetic field of 1 µG, and an inner scale within the Galactic center scattering screen of r in ≈ 800 km Issaoun et al. 2021), we estimate that the root-mean-square phase difference between different polarization modes at 1.3 mm is of order Note that the wavelength dependence of δφ 2 is now nominally dropped. But the implicit dependence of the wavelength is within θ diff , which is the typical value at 230GHz. As a result, for the purposes of the EHT, we may safely assume that the Galactic center scattering screen is non-birefringent.

Characterizing Turbulent Substructure in Images
In this paper we are primarily interested in the measurement and characterization of the statistical properties of small-scale fluctuations in the underlying image of Sgr A*, presumably arising due to turbulent structures in the accretion and/or jet launching region. Typically, these are expressed in terms of power spectra, which measure the degree of fluctuations on each spatial scale.
As the Fourier transform of the sky brightness map, the V (b) already directly contain a measure of the degree of structure on various spatial scales. Therefore, it is natural to construct statistical measures of the variability on different spatial scales from the V (b). For reasons that will become clear in following sections, we choose to do this with the full visibilities, and characterize the spatial power spectrum of image fluctuations via, where the superscript S indicates which Stokes parameter is used to construct the visibilities, and . . . turb indicates averages over timescales long in comparison to the turbulent timescales in the source on the spatial scales of interest (typically many hours or longer). We will presume henceforth that all time averages will include both turbulent averages and ensemble averages, i.e., independent averages over realizations of the stochastic source structure and the intervening scattering screen, dropping the specifier in what follows. Note that this does not subtract the mean V (b), and therefore contains contributions from the variable and static components of the image. Nevertheless, we will find that this is the more convenient power spectrum for scattering mitigation. It is also defined consistently with the power spectral densities in Georgiev et al. (2022).

Impact of Scattering on the Power Spectrum
After scattering via a thin screen, as described in Section 2.1, the power spectrum is modified. The resulting expression may be found in Eq. 32 of Johnson & Gwinn (2015), which reads is the Fourier transform of the structure function. Hidden within Equation 23 is a convolution that, like for diffractive scattering, expresses the impact of refractive scattering as a linear operator, defined by the tensor convolution kernel, Expressing Equation 23 in terms of K(y), we obtain, where both indices of K(y) are summed over. Equation 26 is equivalent to Eq. 16 of Johnson & Narayan (2016) after identifyingD φ with their Q, up to appropriate scalings. Within Equation 26, the impact of diffractive and refractive scattering are clearly delineated by the first and second terms, respectively.

Characteristic scales for K(y)
The typical scales for K(y) may be inferred from its definition and the approximate limiting expressions for D φ (x), and thusD φ (y). On very small scales, D φ (x) is generically quadratic, smoothly vanishing at x = 0. Assuming isotropy of the phase screen, on very large scales, D φ (x) is a power law fixed by the nature of the turbulence within the ISM that gives rise to the scattering screen. Therefore, following Johnson & Gwinn (2015), we express D φ (x) in terms of these two regimes, separated by a spatial scale within the scattering screen, r 0 , which we will assume is much smaller than the scale at which the ISM turbulence is damped, the "inner scale", r in : Typically, the longest baselines accessible to the EHT array are well into the power-law regime, assuming the inner scale of 800km.
In this case, the above expressions simplify to where The correspondingD φ (y) is given in Eq. 34 of Johnson & Gwinn (2015), The collection of constants in front of r 2 diff evaluate to -7.09 for α sc = 1.38, appropriate for Sgr A* (Issaoun et al. 2021). Inserting this into Equation 25, the refractive scattering kernel is approximately

Approximate Inversion of Refractive Scattering
Equation 26 provides P S obs (b) in terms of a linear operation upon P S int (b). Because it is the latter that is of particular interest here, we need to invert this relation, giving P S int (b) in terms of P S obs (b). How to do this is discussed in detail in Appendix B. If the refractive term is small, this inversion can be constructed perturbatively, yielding to first order The accuracy of this approximation is dependent on the magnitude of P S which must be small. This ratio is approximately the fraction of the observed power due to refractive scattering, which is what we are explicitly expanding in.
It also indicates the origin of the scattering-mitigation strategy pursued by combining multiple polarization modes. The idea is to minimize this ratio. It is now clear how this may be done. Very red or blue P S obs (b) will distribute power from large or small scales, respectively, throughout P S int (b). Therefore, suppressing strong variations with spatial frequency in P S obs (b) is primary way in which the choice of polarization mode can impact the fidelity of the P S int (b) reconstruction. At the same time, it explicitly identifies that this approximation cannot be satisfied at all baseline lengths. At sufficiently long baselines, P S obs (b) is exponentially suppressed by the diffractive scattering. As a result, the power distributed from short baselines by the convolution term will inevitably dominate. As a result, it is exponentially difficult to push toward longer baselines in the exponentially suppressed regime. If we define b diff to be the shortest baseline for which D φ [b diff /(1+M )] ≈ 1, then the perturbative expansion will be poorly justified for baseline lengths above |b max |, defined by If we assume that at long baselines the intrinsic spatial power spectrum has a power-law fall off, i.e., P S int (b) ∼ |b| −α , then with D φ (b) andD φ (y) given by Equation 28, Equation 30, and Equation 31 where K 0 encapsulates the α sc -dependent coefficient preceding the factor of r −4 F in Equation 31. Therefore, the condition that the refractive and diffractive contributions to P S obs are similar becomes, and therefore, the maximum baseline length at which we may perturbatively invert Equation 26 to obtain P S int is where we have ignored a logarithmic term that scales as For α sc = 1.38, R = 2 kpc and D = 6 kpc, r F ≈ 10 5 km, and b min ≈ 3 Gλ (corresponding to the longaxis of the diffractive scattering kernel at 1.3 mm), this gives that b max ≈ 3.4b diff (1 + M ) ≈ 13.6 Gλ.
Note that because b diff grows as λ −2 and r F grows as λ 1/2 , b max grows more rapidly than the maximum baseline as λ shrinks. For example, if λ decreases from 1.3 mm to 0.87 mm, b diff grows by a factor of 2.2, while the remainder of the coefficient in Equation 36 decreases to 3.1, and thus b max increases to 32.8 Gλ. In comparison, the maximum Earth-bound baseline grows from 8.5 Gλ to 12.8 Gλ. Hence, at 0.87 mm, the linear approximation improves dramatically.
Henceforth, we will assume that we may utilize the first order inversion approximation in Equation 32 to recover P S int (b) from observations of P S obs (b).

Exploiting Polarization
The dominant impact of scattering on P S int (b) is the reduction of power at long baselines due to diffractive scattering. However, because the action of scattering is independent of the polarization of the observed radio wave, combinations of polarized power spectra may be constructed such that the diffractive suppression is cancelled identically. That is, with Equation 32 applied to the power spectra measured for Stokes parameters S and S , we have which is impacted only by refractive scattering. In this way, diffractive scattering can effectively be mitigated, even without an explicit model for the scattering process (i.e., a specific choice of D φ (x)).
The reconstruction of P S int (b) from P S obs (b) in Equation 32 improves dramatically when P S obs (b) is small on short baselines, and therefore there is less refractive contamination at long baselines to the recovered intrinsic power spectrum. This is reflected by the smaller contributions from the convolution term in Equation 32. When this may be neglected, Equation 37 reduces to Given the measurement of the full Stokes maps, it is generally possible to construct specific polarization modes for which the assumptions underlying Equation 38 are satisfied. Because the P S obs (b) are typically "red", it will suffice to select the polarization modes to preferentially suppress the short-baseline (large-scale) power. By doing so, the contamination of the estimated P S int (b) at long baselines (small scales) from the P S obs (b) at short baselines (large scales) can be nearly eliminated. That is, by minimizing the large scale power in the P S int (b), we are able to minimize the magnitude of the ratio of P S , rendering Equation 38 an excellent approximation at sufficiently long baselines. We begin by constructing the Stokes vector, S 0 , associated with the time-averaged values of zero-baseline Stokes visibility maps Q, U , and V (corresponding the source-integrated polarization). This vector is shown in Figure 2 after projecting it onto the Poincare sphere. We construct two additional polarization modes, S 1 and S 2 , chosen to be orthogonal to S 0 . As a direct result, S 1 = S 2 = 0 identically, and therefore P S1 obs (0) and P S2 obs (0) are generally small. We choose to construct S 1 from V Q (0) and V U (0), thereby ensuring that it corresponds to a linearly polarized mode. The resulting S 1 , after enforcing orthogonality with S 0 is shown in Figure 2. The second polarization is then unique defined up to a sign by the requirement that S 2 be orthogonal to both S 1 and S 0 , as shown in Figure 2. Generally, S 2 will be an elliptical polarization mode.
Explicitly, in terms of S 0 = (q 0 , u 0 , v 0 ), S 1 and S 2 can be expressed as where q 1 , u 1 , q 2 , u 2 and v 2 are the projected components of the Stokes vector onto Q, U and V axes: The power spectra for the two polarization modes can be written down as the linear combination of the power spectra associated with the polarization components Q, U and V , with the coefficients q 1 , u 1 , q 2 , u 2 and v 2 , i.e., The benefits of these two constructed polarization modes are two-fold. First, the Stokes vectors of S 1 and S 2 are perpendicular to S 0 , so they both have zero-mean identically, which ensures minimal impact from scattering. First, the absence of a mean value for S 1 and S 2 renders Equation 38 an excellent approximation of the relationship between the observed and intrinsic power spectra. Second, they have a clear physical meaning.
Fluctuations in S 1 correspond to variations in the observed electric vector position angle (EVPA), i.e., the orientation of the linear component of the polarization. For synchroton sources, as Sgr A* is believed to be, this maps the projected orientation of the net magnetic field as measured on different spatial scales. Fluctuations in S 2 correspond to variations in the observed ellipticity, Figure 2. Relative orientations of S0, S1 and S2 are shown on the Poincare sphere. S1 is orthogonal to S0 and is restricted to lie the Q-U plane. S2 is orthogonal to both S0 and S1. See the text for explicit definitions.
i.e., the degree of circular polarization relative to that of the linear polarization. For synchroton emission from ion-electron plasmas, again, anticipated to be the case in Sgr A*, this directly maps to the angle between the magnetic field and the line of sight. Therefore, these two polarization modes are intrinsically probing the stochastic variability in the three-dimensional magnetic field, integrated throughout the emission region.
MHD turbulence is expected to generate large variations in the plasma density, magnetic field strength, and magnetic field orientation. Thus, given global simulations of the emitting plasma, testable predictions for P S1 int (b) and P S2 int (b), and hence their ratio, may be generated. That is, not only is P S1 int (b)/P S2 int (b) technically easier to measure, but it is precisely the quantity that is expected to provide direct insight into the astrophysical processes within the source.

Signatures of Temporal Variability
In Equation 22, we intentionally did not construct the more common power spectrum of the fluctuations about the mean, (Johnson & Narayan 2016). Nevertheless, observed and intrinsic estimates for ratios of the µ S (b) may also be constructed: and their relationship to the corresponding ratios of P S (b) provide evidence for intrinsic source variability.
If the intrinsic source is stationary, and the ratio of P S (b) and |µ S | 2 (b) are identical. In contrast, when the intrinsic source is variable, where Σ 2 S (b) is the variance due to temporal fluctuations in the source structure on spatial scales b. Thus, when the source is intrinsically variable, P S int (b) is strictly larger than |µ S int | 2 (b). Unfortunately, in the absence of prior knowledge about the nature and magnitude of the variability in the two Stokes parameters, we do not know a priori if the ratios are larger, smaller or equal. Nevertheless, if at a statistically significant degree, then the source must be intrinsically temporally variable. That is, the variability cannot simply be due to the impact of scattering.

VALIDATION WITH SIMPLE SOURCE STRUCTURE
In this section, we present numerical experiments of scattering mitigation with simple source structures for which we have full control of the relevant power spectra. The purpose of this section is to demonstrate: 1. The ability to construct simple source structures with reasonable power spectra that are similar to the target source, Sgr A*,

That the approximation in Equation 38
is well justified and successfully permits reconstruction of probes of the intrinsic variability.
We begin with a description of how toy image models with different power spectra and polarization modes are constructed. This is followed by a set of simulated observations in which scattering is incorporated using eht-imaging (Chael et al. 2018;Chael et al. 2022). Power spectra are constructed and Equation 37 for various choices of S and S are compared for the intrinsic (pre-scattered) and observed (post-scattered) images. In all cases we assumed a wavelength of 1.3 mm.

Constructing Structured Intrinsic Image
The toy model is comprised of a Gaussian delta ring envelope and set of over-imposed fluctuations with a known power spectrum. In more detail: 1. A mean background image is chosen, G(x). For all experiments reported in this section, we adopt a Gaussian delta ring with radius of 25 µas and with width of 5 µas 2. A power spectrum for the fluctuations is chosen, i.e., for some normalization σ ℘ and fluctuation spectral index α, where X is some maximum spatial scale. Because both Sgr A* and GRMHD simulations exhibit fluctuations dominated by those on the largest spatial scales, we will assume that α < 0 generally.
3. A realization of fluctuations are constructed from a set of zero-mean, unit variance Gaussian random variables (GRVs). That is, on a grid in the Fourier domain, at each k we choose two GRVs, N 1 and N 2 , from which the Fourier components of the fluctuation map are given by where we have suppressed the remaining arguments of the fluctuation power spectrum. The spatial f (x) is constructed by the FFT of the F k .
4. In the image domain, the desired model for the total intensity is obtained via where we exponentiate f (x) to ensure positivity.
5. Polarized images are generated in a similar fashion as described above, with two additional Gaussian random fields, l 1 (x) and l 2 (x) constructed similarly to f (x) but with independent α L1 and α L2 , from which L 1 (x) = l 1 (x)I(x) and L 2 (x) = l 2 (x)I(x). (50) These have zero mean by construction.
The definition and the setup the toy model are controlled by five parameters, which encodes five aspects of the desired properties of the toy model: the amplitude of the fluctuations, σ ℘ , the spectral indexes of the power spectra for the total intensity, α, and two polarization models α L1 and α L2 , and a spatial scale on which power spectrum flattens, X, which we will set to 25 µas. Now, we have successfully generated the toy model, with both the total intensity and the polarized components, and the next step is to simulate observations.

Simulated Ensemble of Observations
Observations of Sgr A* are impacted by two additional effects: variability and the interstellar scattering we seek to mitigate. To apply the scattering we make use of the Stochastic Optics package within eht-imaging, which implements the scattering model described in Johnson (2016) with the parameters measured in Johnson et al. (2018) and Issaoun et al. (2021). For all simulated Stokes map (I, L 1 , and L 2 ), a single realization of the scattering screen is employed for each simulated instantaneous image.
Variability, both within the source and the scattering screen, is incorporated by producing a large collection of scattered images, each with a unique randomly constructed scattering screen and set of intrinsic fluctuations. In this way, we generate a statistical ensemble of observed images. Note that this procedure ignores potential temporal correlations within the image that will present themselves in the following section, when GRMHD models are considered.
Specifically, in this numerical experiment, we chose the following values for the parameters for this simple source structure model. The amplitude of the fluctuation of the Gaussian delta ring envelope, σ ℘ , is 100. The power index for the total intensity, α, is −4. The power indices associated with the two polarization modes, α L1 and α L2 , are −4 and −2, respectively. In conclusion, the simple source model has a red power spectrum, dominated by large-scale fluctuation. The associated two polarization modes also have independent red spectra, but with different power indices.
With this source model, we first calculate its intrinsic visibility, for the total intensity and the polarization modes, noted as V I int (u), V L1 int (u) and V L2 int (u). Second, we generate the scattered image using the scattering model implemented in eht-imaging. Each snapshot of the source is scatted with an independent realization of the scattering screen. Third, we calculate the scattered visibility, noted as V I obs (u), V L1 obs (u) and V L2 obs (u). Then, the intrinsic and scattered visiblities are averaged down over 100 realizations of the source and the scattering screen. Finally, we take the ratios among averaged scattered visibilities and averaged intrinsic visibilites, respectively.
For this first set of validation tests, we ignore measurement uncertainties, e.g., those associated with thermal fluctuations. Thus, the primary source of uncertainty is sampling error associated with the finite number of simulated images in the ensemble. In Section 4, we include more realistic assessments of array performance for EHT and ngEHT.

Power Spectra Estimation Results
We first examine the impact of the scattering on the power spectra of the total intensity, which is shown in the left panel of Figure 3. The central peak is associated with the net source structure, i.e., a unit Jy Gaussian delta ring with radius of 25 µas. The ringing is associated with the Fourier transform of the the delta ring as the envelope. The extended plateau at baselines longer than 15 Gλ is associated with two effects: first the small- Figure 3. Power spectra for the total intensity and two linear polarization modes. From left to right, the three panels are for total intensity, first linear polarization and second linear polarization. The black and red lines are the intrinsic and scattered power spectra, respectively. The blue dotted lines are the diffractively scattered power spectra. The green dashed line in the leftmost panel is the power spectrum of the Gaussian delta ring with radius of 25 µas, which is the same envelope for the simple source model. scale variable structures that we seek to recover, and second the refractive scattering.
It is, in fact, the latter of these two that overwhelmingly dominates, as evidenced by the impact of diffractive scattering, shown by the dashed blue line, which strongly suppresses the contributions from the intrinsic source structure at u 10 Gλ.
The power spectra associated with the polarization, shown in the center and right panels of Figure 3, do not have a prominent central peak because the total polarization flux vanishes on average (though while small, is non-zero for any given realization). Thus, while again the diffractive scattering suppresses the power at long baselines due to intrinsic structure, there is much less contamination from refractive scattering, which dominates only for u 20 Gλ.
In all cases, the power spectra after scattering deviate substantially from those intrinsic to the source. That is, as anticipated, due to both diffractive and refractive effects, the observed power spectra are themselves a poor proxy for the intrinsic power spectra, becoming worse as the spatial scale decreases.
The ratios of observed power spectra, as defined in Equation 38 and shown in Figure 4, produce remarkable agreement with those from the intrinsic images (i.e., prior to scattering). When the ratio is made with P I , the intrinsic source structure and positive definite nature of the intensity responsible for the large bump in the left panel of Figure 3 introduces a large depression, overwhelming any structure that may be attributed to the small-scale fluctuations. In addition, the dominance of refractive scattering at long baselines in P I results in a significant departure of the observed from the intrinsic power spectra ratios by u = 12 Gλ as anticipated at the end of Section 2.3.2.
However, in stark contrast, the power spectra ratio of the two polarized modes, show in the right panel of Figure 4, match well out to u ≈ 17 Gλ. While beyond u ≈ 20 Gλ, refractive scattering drives large deviations from intrinsic power spectra ratio, this is well beyond the baselines accessible to EHT and ngEHT at 1.3 mm.
Similar experiments were performed for a variety of choices of the α, α L1 , and α L2 , with similar results.
While the above is schematic, involving only a very simple source structure, nevertheless a number of immediate conclusions can be drawn that we will see reflected in the more physically applicable demonstrations that follow: • Observed power spectra are poor estimators in the presence of interstellar scattering for the intrinsic variable structures.
• For all polarization modes, diffractive scattering suppresses long-baseline observed power spectra.
• For polarization modes with large net flux (e.g., Stokes I), refractive scattering substantially contaminates the power spectra.
• Ratios of power spectra generally produced better estimates of the corresponding intrinsic quantities.
• Ratios of power spectra associated with polarization modes with zero net flux are substantially more accurate estimates, extending well beyond the spatial frequency range accessible from the ground at 1.3 mm.
Based on the above, we conclude that when the polarization modes are well chosen, the approximation in Equation 38 is well motivated.

VALIDATION WITH GRMHD SIMULATIONS
In contrast to the simple, phenomenological models discussed in the previous section, GRMHD simulations provide a natural astrophysically-motivated set of complex source stuctures. While GRMHD simulations do . Power spectra ratio between the total intensity and the first polarization mode (left), the total intensity and the second polarization mode (middle), and the two polarization modes (right) for the simple source structure defined in Section 3. The black line shows the intrinsic power spectra ratio, while the red line is the power spectra ratio observed after scattering. The green band represents the size of Earth-bound baselines at 230 GHz.
not afford the freedom to arbitrarily modify the input fluctuation spectra, they do incorporate credible realizations of the anticipated turbulence and magnitude of the polarized flux. As a result, it is possible to reasonably assess the practical limitations imposed by thermal noise noise and limited number of observations to be averaged. Here we repeat the kinds of tests performed in Section 3 for one such GRMHD simulation, taken from the set presented in Event Horizon Telescope Collaboration et al. (2022e).

GRMHD Simulated Intrinsic Image and Simulated Observation
We employ a SANE, a = 0, i = 30 • , R high = 40 simulation from the set presented in Event Horizon Telescope Collaboration et al. (2022e), to which we direct the reader for information about the simulation particulars. For our purposes, it is important only that the simulation presents a physically self-consistent realization of the kind of turbulence, degree of net polarization (∼3%), and typical polarization fractions (∼20%) appropriate for Sgr A*.
An arbitrary snapshot drawn from the simulation is shown in Figure 5 with its four Stokes components, and in Figure 6 with the constructed polarization modes S 0 , S 1 and S 2 , as stated in Section 2.4.
Note that because our goal here is not to predict the statistics of GRMHD simulation images, but rather to demonstrate the ability to faithfully retrieve statistical elements of the underlying intrinsic images in the presence of scattering, this single GRMHD simulation is sufficient for our purposes.
The simulated emission is assumed to arise from synchrotron emission due to a population of hot electrons. The simulation data contains the total intensity and polarization maps at 1.3 mm for accretion flow parameters relevant for Sgr A*, i.e., the four Stokes parameters, I, Q, U and V , for 3000 individual snapshots (Event Horizon Telescope Collaboration et al. 2022e, and references therein). From these intrinsic images, scattered images are produced using the Stochastic Optics package within eht-imaging in manner identical to that used in Section 3.
We follow a nearly identical procedure to generate simulated observations from the GRMHD simulations as that described in Section 3.2. This procedure differs in that the polarization maps are now identified with the three polarized Stokes maps (Q, U , V ). Where these are scattered, we produce a new scattering screen realization for each image. For each frame, we construct the single-snapshot estimate of the spatial power spectra associated with I, Q, U , and V are constructed, i.e., V I 2 , V Q 2 , V U 2 , and V V 2 , respectively. This procedure is repeated for every image in the GRMHD simulation to generate estimates for the ensemble-and turbulenceaveraged estimates of the spatial power spectra.

Spatial Power Uncertainty Estimate
To assess if the observed and intrinsic spatial power spectra are distinguishable, we require an estimate of the anticipated uncertainty on the spatial power spectra. This arises from multiple potential origins. However, there are two irreducible contributors: thermal noise associated with the individual stations within the EHT and ngEHT, and the sampling uncertainty due to an insufficiently complete ensemble. Here we describe how we estimate each of these. Note that we make aggressive assumptions to reduce both sources of uncertainty, thereby enforcing a stricter limit on the required fidelity of the spatial power spectra ratios.

Thermal Noise Estimates
The thermal error at EHT stations arises from a number of potential sources, including the atmosphere, side lobes picking up the local environment, and the electronics within receiver. The combination is typically characterized by a system equivalent flux density, SEFD. In Figure 5. Arbitrary snapshot drawn from the GRMHD simulation used. The figure includes the Stokes I, Q, U and V components (left) and their corresponding scattered versions (right). Bottom row: the two constructed polarization modes, S1 and S2, are shown, both intrinsic and scattered.
terms of these, the uncertainty on a visibility measured by stations A and B with a bandwidth ∆B and coherently averaged over a time τ , SEFDs for EHT stations are listed in Table 2 of (Event Horizon Telescope Collaboration et al. 2019c). These range from 74 Jy for ALMA to 19300 Jy for SPT. We adopt, for illustration, PV and APEX, for which the SEFDs are 1900 Jy and 4700 Jy, respectively, the intermediate SEFDs in the EHT. Were we to adopt the two stations with the highest SEFDs in the EHT, SMT and SPT, the estimated thermal noise could be a factor of 30 higher. The median thermal noise across the EHT is Figure 6. The same arbitrary snapshot drawn from the GRMHD simulation used in Figure 5. The figure includes the two constructed polarization modes, S1 and S2, both intrinsic and scattered.
roughly an order of magnitude larger than the minimal value. Finally, we adopt a bandwidth of 4 GHz, corresponding to the combination of high-and low-band data from the EHT, and an integration time of 10 min, corresponding to a typical scan time, yielding σ th ≈ 0.4 mJy. Upon averaging N independent observations, the effect thermal noise is reduced by a further factor of N −1/2 . Note that this is independent of the particular polarization mode under consideration.
The thermal noise on the spatial power spectra after averaging N independent observations is obtained via standard error propagation, The thermal noise on the ratio of spatial power spectra after averaging N observations is, The sampling noise describes the uncertainty associated with having a finite number of samples in the estimate of the ensemble and turbulence averages. When the number of samples, N , is large, the central limit theorem implies that this is related to the variance of the visibility amplitude, i.e.,

Sampling Noise Estimates
Note that unlike the thermal noise, this differs between the various Stokes parameters, which may exhibit different degrees of variability. From Σ V A , the uncertainties on the spatial power spectrum and spatial power spectra ratios can be immediately constructed via the standard error propagation, and

Power Spectra Estimation Results
The mean power spectra associated with Stokes I, Q, and S1 are shown in Figure 7. As with the toy model presented in Section 3, the non-zero total flux results in a peak at u = 0 Gλ, and a deficit associated with the diffractive component of the scattering at long baselines. As with Figure 3, refractive scattering lessens the reduction, seen most prominently in P I because of the comparatively large net value. In all cases, scattering significantly suppresses the power spectra relative to their intrinsic values, as anticipated.
The mean P Q /P I and P S2 /P S1 power spectra ratios are shown in Figure 8, and are directly comparable to those in Figure 4. In addition, the sampling and thermal error scales for N = 25 independent observations of the source and scattering screen are indicated, providing a natural assessment of the accuracy of the approximation in Equation 38. Apart from differences in the underlying source structure, e.g., the clearly evident oscillations associated with the lensed emission ring, the gross properties noted in Section 3 remain: the suppression at short baselines by the non-zero mean total flux in the P Q /P I ratio (similar to P U /P I ), and lack of such a suppression in the Stokes basis defined by S 1 and S 2 . Even at the longest ground-based baseline lengths -10 Gλ at 230 GHz -the impact of the scattering screen is effectively mitigated.
While Figure 8 shows only the ray through the uvplane along the u-axis, this improvement is generic. Figure 9 presents P S2 /P S1 for radial rays at different orientations. While the magnitude of the discrepancies between the observed and intrinsic power spectra ratios and the location where they begin to differ varies, in all cases, at all baselines relevant for EHT and ngEHT, these discrepancies are small in comparison to relevant uncertainties.
In addition to vastly increasing the number and density of baselines available, the ngEHT envisions receivers with increased bandwidth and improved detector efficiency that will improve sensitivity across the array, reducing the thermal noise contributions to the visibility uncertainties (Raymond et al. 2021). For nominal SEFDs associated with the reference ngEHT array in Raymond et al. (2021) and a bandwidth of 16 GHz, we show P S2 /P S1 at the observing frequencies under discussion for the ngEHT in Figure 10. As shown in the middle panel of Figure 10, for the ngEHT, thermal noise ceases to be the chief impediment to measuring the intrinsic power spectra ratios below ∼ 15 Gλ, covering all Earth-sized baselines. The right panel of Figure 10 shows P S2 /P S1 at 345 GHz with ngEHT telescope sensitivities and bandwidths; again all ground-base baselines are dominated by the sampling uncertainty associ- Figure 8. Ratios of the power spectra of different polarization modes. From left to right, the three panels show the power spectra ratio for P Q /P I , P U /P I and P Q /P U . The black and red line represents before and after the imposition of scattering. The blue error bars are the thermal noise associated with telescopes. We used PV and APEX, which have the intermediate sensitivities among all EHT telescopes, to generate the error bars. The grey error bands are the sampling noise associated with the intrinsic spatial variability of the intrinsic source, which we averaged down assuming 25 independent observations. The green band represents the size of Earth-bound baselines at 230 GHz. Figure 9. Intrinsic (black) and observed (red) P S2 /p S1 along rays at different orientations relative to the u axis. From left to right: 0 • (u axis), 45 • , and 90 • (v axis). The sampling noise and thermal errors associated with N = 25 observations are indicated by the grey band and blue error bars, respectively. Figure 10. Mean intrinsic (black) and observed (red) ratio of the power spectra associated with polarization modes defined by the Stokes vectors S1 and S2 at 86 GHz (left), 230 GHz (middel) and 345 GHz (right). The sampling noise and thermal errors associated with N = 25 observations and based on ngEHT telescope sensitivities are indicated by the grey band and blue error bars, respectively. The green band represents the size of Earth-bound baselines at 86 GHz, 230 GHz and 345 GHz, from left to right. ated with source variability. Out to ∼ 50 Gλ, nearly twice the 30 Gλ region shown in Figure 10, the observed mean power spectra ratio provides an excellent estimate of the intrinsic ratio. Within ∼ 100 Gλ observed and intrinsic power spectra ratio differ by less than the sampling uncertainty. These imply that future high-frequency polarimetric observations by EHT, ngEHT, and space-based mm-VLBI experiments will all be able to accurately statistically probe the turbulent structures in horizon-scale targets on scales from 7 µas (ground) to 2 µas (space). Whereas at lower observation frequencies, the impact of scattering is enhanced. In the left panel of Figure 10, which is at 86 GHz, although scattering deviates the power spectrum ratio from the intrinsic at as early as 3 Gλ, which is the Earth-based VLBI baseline limit at this frequency, they are within the allowed sampling uncertainty until 20 Gλ. This implies the possibility to apply this scattering mitigation scheme to data from GMVA and other telescopes with lower observation frequencies.
The ability explicitly distinguish between intrinsic source variability and extrinsic evolution in the scattering screen using power spectra and mean ratios is demonstrated in Figure 11. Like P S2 /P S1 , the ratio of the means, |µ S2 | 2 /|µ S1 | 2 , are insensitive to the scattering for baselines relevant for current and future Earthbound VLBI arrays. When the source is static, here taken to be a single GRMHD snapshot, the two sets of ratios are identical out to baselines well in excess to those accessible on the ground, ultimately limited by diffractive scattering. In contrast, for variable sources, the two ratios differ at high significance. Therefore, the comparison of these two ratios, i.e., P S2 /P S1 and |µ S2 | 2 /|µ S1 | 2 , presents a way in which to find direct evidence for intrinsic source evolution.
Which component of the uncertainty dominates depends on baseline length. At sufficiently large u the strong suppression due to diffractive scattering reduces the signal precipitously. The lower intrinsic S/N of the polarized data result in characteristically larger thermal uncertainties on the power spectra ratio. At short baselines, the sampling error dominates. These two regimes conceptually differ in the manner that measurements can be improved. The thermal noise can be reduced by improvements to station sensitivity (e.g., increased bandwidth, dish size, phase referencing, etc.). In contrast, the sampling noise is a consequence of the intrinsic variability alone, and can only be improved by repeated observation. Because additional observation epochs also reduce the thermal noise at the same rate, where the sampling noise dominates, it will do so regardless of the number of observations.
The location of the transition from sampling dominated to thermally dominated noise depends on the nature of the variability and the sensitivity of the individual stations (see Figure 12). For the simulation considered here and the median thermal noise, this transition Figure 11. Ratios of the power spectra (solid lines) and means (dashed lines) for the two constructed polarization modes, S1 and S2. Top: a single, static intrinsic image, represented by a single GRMHD snapshot, viewed through an evolving scattering screen. Bottom: an evolving intrinsic image and scattering screen. In both panels, the sampling noise and thermal errors associated with N = 25 observations and based on ngEHT telescope sensitivities are indicated by the grey band and blue error bars, respectively. The green band represents the size of Earth-bound baselines at 230 GHz. occurs just beyond 10 Gλ. However, for the baseline with the maximum thermal noise, the uncertainty at all baselines is dominated by the thermal component, implying that sensitivity of existing EHT stations is likely to be the limiting factor in the accuracy with which the power spectra ratios can be measured.

CONCLUSIONS
Interstellar scattering associated with turbulence in the ISM is nonbirefringent for physically reasonable Figure 12. Comparison of the sampling and thermal noise estimates for P S2 /P S1 at 230 GHz. The thermal noise for the minimum, median, and maximum noise estimates are shown by the dotted, solid, and dashed blue lines. The sampling noise is indicated by the grey line. magnetic field strengths. As a consequence, the effect of scattering on horizon-resolving polarization maps of Sgr A* is expected to be independent of the polarization mode being observed. This presents an opportunity to statistically separate the small-scale structures induced by refractive scattering and those intrinsic to the source, presumably due to turbulence within the near-horizon emission region.
We characterize the statistical properties of the structural variability by their power spectra, P S , defined to be the mean squared visibility associated with Stokes parameter S. This definition is convenient because the impact of scattering is a linear operator that may be expressed as a tensor convolution acting on the P S . As far as S/N permits, this convolution may be inverted via a perturbative expansion; for existing and proposed Earth-sized mm-VLBI arrays, like the EHT and ngEHT, only the first term in the expansion is required. This effectively reduces to applying constraints to the ratio of P S , and is otherwise insensitive to the details of the scattering screen.
It is possible to select the polarization modes to minimize the impact of refractive scattering and simplify the interpretation of the P S . We do this by constructing a particular basis of Stokes vectors, S 1 and S 2 , that are orthogonal to the source-integrated mean Stokes vector. These correspond to the fluctuations in the projected orientation of the magnetic fields on the sky and along the line of sight, respectively. Thus, the ratio P S2 /P S1 is a direct measure of the degree of isotropy in the MHD turbulence across spatial scale.
Using both a toy model, in which there is substantial control over the properties of the turbulent structures, and a GRMHD simulation, which contains a realistic representation of MHD turbulence, we have demonstrated the ability to reconstruct various intrinsic power spectra ratios, including P S2 /P S1 . At all baselines accessible to ground-based mm-VLBI experiments, the difference between the reconstructed and intrinsic P S2 /P S1 is small in comparison to the uncertainties due to finite sample size and measurement errors. This remains true for other power spectra ratios between polarized components (e.g., P Q /P U ) on baseline lengths of interest, and at higher observation frequencies.
The improved sensitivity expected in future mm-VLBI experiments can significantly increase the accuracy with which P S2 /P S1 may be measured. More importantly, arrays like that envisioned by the ngEHT provide a much more dense sampling of the uv-plane, and therefore the ability to measure more completely the twodimensional power spectra ratios. The thermal uncertainties may be further reduced by aggregating nearby measurements in the uv-plane and/or exploiting assumptions regarding azimuthal symmetry; a complete discussion of these will appear elsewhere.
At higher observation frequencies, e.g., 345 GHz, the reduced impact of scattering results in nearly exact mitigation out to baseline lengths of ∼ 50 Gλ and within the 25-epoch sampling variance for baseline lengths up to ∼ 100 Gλ. Thus it is possible to effectively mitigate interstellar scattering on baselines relevant for spacebased mm-VLBI concepts that place stations in low and medium Earth orbits (2,000 km and <35,000 km, respectively), and accurately probing the magnetic field power spectrum on scales as small as 2 µas and thus a quarter of the Schwarzschild radius in Sgr A*.

ACKNOWLEDGMENTS
We would like to thank Vedant Dhruv and many others from University of Illinois at Urbana-Champaign, who provided us with the GRMHD simulation. We would also like to thank Ramesh Narayan for his helpful comments. The scattering screen can be envisioned as a screen of width L comprised of many electron bubbles, which have typical size of L ⊥ . When a photon travels travels through the scattering screen, each bubble causes a slight deflection to the photon's trajectory. After averaging the photon deflection angle over the width of the scattering screen, we should be able to derive the root-mean-squared deflection angle as a function of the width.
The argument is also valid for the phase change caused by the scattering screen, as the deflection angle and the phase change are linearly correlated. Next, we are going to consider two different scenarios, where the autocorrelation function of the electron density and the magnetic field takes different forms.

A.1. Difference in the Deflection Angles for Different Polarization Modes
We have shown in the Section 2.2 that the difference in the deflection angles and the phase changes for different polarization modes can be written as an integral over functions of the magnetic field and the electron density field, as In this appendix, we will explore the properties of δθ and δφ, given different fluctuations of the electron density fields and the magnetic fields.
We assume that X and Y , the fluctuation electron density field and the magnetic field, are some independent Gaussian random fields. For the electron density field, it has mean value X 0 , while for magnetic field the mean value is zero. Decomposing X and Y into Fourier modes with cylindrical coordinate system, we have where q and r are the Fourier conjugate in the cylindrical plane, and m and z are conjugate in the axial direction.
The power spectra for the two Gaussian random fields are defined as: (A6) The independence of the two fields require the Fourier coefficients a q,m and b q,m to satisfy a q 1 ,m1 b * q 2 ,m2 = 0. For simplicity, from here on in this appendix, we will use prime to denote the derivative in the radial direction, e.g. Y = ∂Y /∂r.
Inserting the fields X and Y expanded in the Fourier domain back to Equation A1, and taking the ensemble average, we have A quick check can be done that because a q,m and b q,m to satisfy a q 1 ,m1 b * q 2 ,m2 = 0, the average of δθ is zero. Similarly, we can insert Equation A3 and Equation A4 into δθ 2 , which takes the form where the subscripts 1 and 2 of X and Y denote z 1 and z 2 dependence of X and Y , and the subscripts 1, 2, 3 and 4 of q and m denote different realizations of the Gaussian random fields.
The integrand above can be divided into nine different terms, as (A9) Each nine components of the integral can be done independently using the identities Equation A5 and Equation A6.
The first term is: The second term is: The third term is: The fourth term is: The fifth term is: (A14) The other four terms are zero, because a q 1 ,m1 b * q 2 ,m2 = 0.
Grouping nine terms together, we have (A15) We can define the auto-correlation function of fields X and Y to help simplifying Equation A15: At the line of sight, Equation A15 can be evaluated with the auto-correlation functions C (r = 0, z) and D (r = 0, z): Because the derivative is with respect of the perpendicular direction, and the integral is along the line-ofsight direction, derivative and integral can be switched:

A.2. Difference in Phase Changes for Different Polarization Modes
Similar to that Equation A1 can be expanded in the Fourier space and expressed in the compact form of the auto-correlation functions, same can be done to Equation A2.
(A22) The equation above can be evaluated in the same way as we evaluate the variance of δθ, with the help of the auto-correlation functions: Compared to Equation A19, the variance of δθ is proportional to the second order derivative of the variance of δφ: where L is the typical correlation length within the plasma.

A.3. Explicit Examples and Quantitative Estimates
Quantitatively assessing the magnitude of the phase differences between the two polarization modes requires an explicit model for the density and magnetic field fluctuations within the scattering screen. Because Equation A19 and Equation A23 depend only on the correlation functions of these underlying physical quantities, specifying the statistical properties of the density and magnetic field through their correlation functions is sufficient. Because these are not known a priori, we explore a handful of examples, beginning with a Gaussian correlation functions with a natural intrinsic scale, and culminating in the Kolmogorov models that are traditionally employed.

A.3.1. Gaussian Auto-correlation Functions
We first consider the simple case that the autocorrelation functions of the electron density, X, and the magnetic field, Y , are Gaussian. These have a clear scale, beyond which the correlation is exponentially suppressed. Both the Gaussian distributions have the standard deviation of L ⊥ , which is the typical size of fluctuations in the electron density and magnetic field strength. Within each bubble, both the electron density and magnetic field fluctuations are highly correlated, while for different bubbles they are not. The auto-correlation functions are: where the corresponding variances σ X and σ Y are of order of one.
To estimate Equation A23, we first make a change of variable:z in terms of which the variance in the phase fluctuations may be written with new upper and lower bounds: Inserting the Gaussian correlation functions above, this becomes, The scattering screen is comprised of many electron bubbles along the line-of-sight, L /L ⊥ 1. In this limit both erf L /L ⊥ is approximately unity and exp −L 2 /L 2 ⊥ is approximately zero. With these simplifications, the variance of δφ becomes, Similarly, the variance of the deflection angle, Equation A19, can be calculated given Equation A25 and Equation A26 in the limit that L /L ⊥ 1. This is facilitated by the fact that application of the second-order derivative within the screen is straightforward within the integrals. The result is, (A32) These two characterizations of the degree of birefringence are related by where the terms in parentheses are generally of order unity. This relationship between δφ 2 and δθ 2 matches the general expectation from Equation A24.

A.3.2. Broken Power-law Auto-correlation Functions
The second case we considered is when the autocorrelation functions take the form of a broken power law, which incorporates fluctuations on multiple scales. Above a minimum scale, L ⊥ , the distribution of density and magnetic field fluctuations is self similar, characterized by a power-law index α, C (r, z) = σ 2 X 1 + (z 2 + r 2 ) 1/2 /L ⊥ α (A34) D (r, z) = σ 2 Y 1 + (z 2 + r 2 ) 1/2 /L ⊥ α .
As a concrete example, we present the computation for α = 2 prior to moving on to a Kolmogorov description; the result is qualitatively similar for any α > 1. The variance in the phase perturbations are, Again we will assume L /L ⊥ 1, and thus arctan L /L ⊥ ≈ π/2 and the logarithmic term is small relative to the linear terms. In this limit, the phase fluctuation variance simplifies to which, up to factors of order unity, matches the expression found for Gaussian auto-correlation functions. The same argument holds true also for Equation A19. After taking the second derivative with respect to r and evaluating at r = 0, integrating of Equation A19 gives, This differs only by factors of order unity from the Gaussian-correlation-function case, with δφ 2 = 4π 2 L 2 ⊥ λ 2 δθ 2 σ 2 X + 2X 2 0 3σ 2 X + 2X 2 0 .
Again, this expression matches the expectation from Equation A24 up to factors that are generally of order unity. This is a little different compared with that in the Gaussian case (Equation A33), but as long as σ 2 X , σ 2 Y and X 2 0 are of order one, Equation A39 also holds true for the broken power law case.

A.3.3. Kolmogorov Turbulence-implied Auto-correlation Functions
The autocorrelation function and the structure function describing the same field are closely correlated. A more physically inspired autocorrelation function needs to be derived from a physical model of the turbulence. Past literature has suggested the power spectrum of the electron density fluctuation (e.g. see Armstrong et al. (1995); Rickett (1990)) as when the scale is between the inner and outer scale, where the factor C 2 N is the structural coefficient, and β is the spectral power index, which in the Kolmogorov case, β = 11/3. The corresponding structure function, D φ , which is the integral form of the power spectrum, is related to Equation 27 where L ⊥ = 2 1/α x 0 can be interpreted as the typical transverse size of the electron bubble.

B. GENERAL INVERSION OF REFRACTIVE SCATTERING
Equation 26 describes a linear relationship between the power spectra associated with the brightness fluctu-ations before and after scattering. However, it is written in such a fashion that the observed power spectrum is a function of the properties of the scattering screen and the intrinsic power spectrum. What is desired is an expression for P S int (b) in terms of the scattering screen and P S obs (b), i.e., to invert Equation 26. This is generally possible since |K(u)| > 0 and D φ (u) < ∞ for all u. In practice, it requires some care since despite being a linear relationship, it is nonlocal due to the presence of the convolution.
Here we develop a general inversion scheme, first using a local representation of the relationship between modified correlation functions, and second making use of a perturbative expansion that exploits the limit in which the refractive scattering kernel, K(u), is small. Equation 38 makes use of the first order approximation in this peturbative solution.

B.1. Alternate Convolution Representation
To further understand the role which K(u) plays in the process of scattering, it is straightforward to express Equation 26 in the convolution representation. In the convolution picture, we are able to acquire the solutions to the observed power spectrum P S obs (b) as required. And we can further assess the impact of the scattering kernel function K(u) on different baselines.
In the convolution representation, We define the spatial correlation functions C S obs (x) = d 2 b e −ib·x/r 2 F P S obs (b), and D S int (x) = d 2 b e −ib·x/r 2 F e −D φ [b/(1+M )] P S int (b). (B54) These are associated with the observed and diffractively suppressed spatial power spectra. In terms of these, Equation 26 may be written as defines a set of general coefficients across all Stokes parameters.

B.2. General Solution
Because H is manifestly Hermitian, Equation B55 can be inverted by solving for the eigenfunctions of (1 − H :