Synchrotron Polarization Signatures of Surface Waves in Supermassive Black Hole Jets

Supermassive black holes in active galactic nuclei are known to launch relativistic jets, which are observed across the entire electromagnetic spectrum and thought to be efficient particle accelerators. Their primary radiation mechanism for radio emission is polarized synchrotron emission produced by a population of nonthermal electrons. In this Letter, we present a global general relativistic magnetohydrodynamical (GRMHD) simulation of a magnetically arrested disk (MAD). After the simulation reaches the MAD state, we show that waves are continuously launched from the vicinity of the black hole and propagate along the interface between the jet and the wind. At this interface, a steep gradient in velocity is present between the mildly relativistic wind and the highly relativistic jet. The interface is, therefore, a shear layer, and due to the shear, the waves generate roll-ups that alter the magnetic field configuration and the shear layer geometry. We then perform polarized radiation transfer calculations of our GRMHD simulation and find signatures of the waves in both total intensity and linear polarization, effectively lowering the fully resolved polarization fraction. The telltale polarization signatures of the waves could be observable by future very long baseline interferometric observations, e.g., the next-generation Event Horizon Telescope.


INTRODUCTION
Accreting supermassive black holes can produce highly relativistic electromagnetically collimated outflows called jets, observed across the electromagnetic spectrum.These jets can be observed up to kilo-parsec scale in the case of active galactic nuclei (AGN).At radio Corresponding author: J. Davelaar jdavelaar@flatironinstitute.org and mm frequencies, the primary emission mechanism is synchrotron emission.Very Long Baseline Interferometric (VLBI) observations probe the jet substructure and reveal edge-brightened morphology, often referred to as limb brightening, see, e.g., (Walker et al. 2018;Kim et al. 2018;Giovannini et al. 2018;Janssen et al. 2021).However, the mechanism responsible for energizing the radiating electrons along the jet surface remains an active debate.The upcoming next-generation VLBI facilities will bring higher resolving power and dynamic range, allowing for better-resolved AGN jet images and polarization maps.In this work, we investigate the imprint of instabilities in the jet boundary and the effect of non-thermal tails in electron distribution functions on polarized emission features of AGN jets.
Since synchrotron emission is intrinsically linearly polarized (Rybicki & Lightman 1979), where the polarization vector is perpendicular to the magnetic field lines, the observed polarization from AGN can be used to study the magnetic field structure of jets.These systems generally show very low linear polarization fractions across the entire radio band, see, e.g., Zavala & Taylor (2003); Hada et al. (2016); Walker et al. (2018); Park et al. (2019); EHT MWL Science Working Group et al. (2021).These low fractions indicate that an external Faraday screen depolarizes the jet's emission before it reaches us or that the emission from the jet is not generated in large-scale coherent magnetic field geometries.In the case of M87 at sub-mm wavelengths, observations by Hada et al. (2016) show linear polarization fractions of up to 20%, revealing regions of coherent field geometry when observed at higher spatial resolution.
The highest resolution polarized observations of a lowluminosity AGN (LLAGN) to date are by the EHT collaboration (Event Horizon Telescope Collaboration et al. 2021).The EHT showed that in M87 * at horizon scales, the polarization vector shows a helical pattern, which is typically reproduced by simulations of accretion flows in the magnetically arrested disk (MAD) state (Igumenshchev et al. 2003;Narayan et al. 2003;Tchekhovskoy et al. 2011).MAD accretion flows can be studied with general relativistic magnetohydrodynamical (GRMHD) simulations.To study the emission generated by the GRMHD simulations, they can be post-processed with general-relativistic radiative transfer codes (Dexter et al. 2012;Mościbrodzka et al. 2017;Davelaar et al. 2018;Wong et al. 2021;Cruz-Osorio et al. 2022;Fromm et al. 2022) after postulating an electron temperature prescription.A magnetically arrested flow typically reaches a limit where the magnetic pressure balances the gas pressure due to accumulated magnetic flux on the event horizon.This limit is often identified with ϕ mad = Φ B / Ṁ ≈ 15 (Tchekhovskoy et al. 2011), here Φ B is the magnetic flux threading the horizon and Ṁ the mass accretion rate.If this threshold is reached, the antiparallel magnetic field lines in the northern and southern jet are compressed to form a thin current sheet that reconnects and expels the magnetic flux (Dexter et al. 2020;Ripperda et al. 2020;Ripperda et al. 2022).During such magnetic flux eruptions, the magnetic field undergoes no-guide-field reconnection, resulting in the expulsion of a flux tube consisting of a vertical (poloidal) magnetic field.This flux tube can push the accretion disk away, effectively arresting a part of the incoming flow.The flux tubes can orbit at sub-Keplerian velocities in the disk, where they can propagate to a few tens of gravitational radii.Within one orbit, the lowdensity fluid in the flux tube gets mixed into the higherdensity disk through magnetic Rayleigh-Taylor instabilities (RTI) (Ripperda et al. 2022;Zhdankin et al. 2023).The magnetic flux eruptions are conjectured to power high energy flares through reconnection near the horizon (Ripperda et al. 2020;Dexter et al. 2020;Ripperda et al. 2022;Hakobyan et al. 2023) and through reconnection induced by RTI at the boundary of the orbiting flux tube (Porth et al. 2021;Zhdankin et al. 2023).
In this Letter, we will use GRMHD simulations in Cartesian Kerr-Schild coordinates with adaptive mesh refinement to study the large-scale properties of the jets.Using Cartesian coordinates in combination with adaptive mesh refinement allows us to better resolve the shear layer separating the highly-relativistic bulk velocity jet and the mildly-relativistic disk.We will refer to this region as the jet-wind shear layer.Our simulation shows that the magnetic flux eruptions associated with MAD flows can drive waves along the jet-wind shear layer.At larger radii, the waves show roll-ups that mix highdensity wind material with the low-density magnetized jet.In this non-linear phase, the waves may trigger magnetic reconnection and turbulence, as was found in local-box simulations (Sironi et al. 2021).To quantify the imprint of the waves on observables, we ray-trace our simulation with our polarized radiative transfer code RAPTOR (Bronzwaer et al. 2018(Bronzwaer et al. , 2020)).We find that the waves depolarize the observed synchrotron emission.
This letter is structured as follows.Section 2 describes our numerical setup and summarizes how we compute synthetic polarized images.Section 3 explains our GRMHD and radiative transfer results, highlights linear and circular polarization properties and provides evidence that the existence of waves results in depolarization.Finally, in section 4, we discuss and summarize our main conclusions.

NUMERICAL SETUP
In this section, we will describe our GRMHD simulation setup, our radiative transfer code, and our electron thermodynamics model.

GRMHD
To model the dynamics of the accretion flow around a Kerr black hole, we use the Black Hole Accretion Code (BHAC) (Porth et al. 2017;Olivares et al. 2019), which solves the ideal GRMHD equations in curved spacetime.
The equation of state is assumed to be an ideal gas law, described via the specific enthalpy h(ρ, P gas ) = 1 + γ ad γ ad − 1 with gas pressure P gas , mass density ρ in the fluid frame, and the adiabatic index is set to γ ad = 13/9.We initialize our simulation with a Fishbone & Moncrief (1976) torus in spherical Boyer-Lindquist coordinates, (t, r, θ, ϕ), here θ is the angle with respect to the spin axis of the black hole and ϕ the azimuthal angle around the black hole spin axis.The initial conditions are then transformed to Cartesian Kerr-Schild coordinates, (t, x, y, z).The initial disk has an inner radius of r in = 20 r g , with r g the gravitational radius given by r g = GM /c 2 , with G Newtons constant, M the black hole mass, and c the speed of light.We set the pressure maximum of the disk at r max = 40 r g .The initial magnetic field is given by the vector potential A ϕ = ρ(r/r in sin θ) 3 e −r/400 − 0.2, where r is the radial coordinate, and θ the polar angle.The vector potential follows iso-contours of density, ρ.These initial conditions are chosen such that the simulation will reach the saturated state of a MAD accretion flow.
We set the dimensionless black hole spin parameter to a = 15/16, which results in a black hole horizon size of r h ≈ 1.34 r g .The domain size is (−4000 r g , 4000 r g ) in all three spatial Cartesian directions.Our base resolution is 192 3 cells.We employ nine additional levels of static mesh refinement, resulting in an effective uniform Cartesian resolution of 98, 304 3 .The highest resolution grid is centered at the horizon and has a resolution of ∆x i = 0.08 r g .The simulation is run for 10, 000r g /c.We introduce a maximum for the cold magnetization parameter σ = b 2 ρ , where b is the magnetic field strength, and we inject mass to maintain σ ≤ 100.Additionally, we ensure that β −1 p ≤ 10 3 , where β p = 8πP gas /b 2 is the plasma beta parameter.We use floor profiles for density as well as pressure, given by ρ floor = 10 −4 r −2 , and P floor = 10 −6 ρ γ ad floor .Due to our Cartesian grid, we do not have an inner radius where we can employ outflowing boundary conditions, we, therefore, introduce an artificial treatment of the fluid variables inside the black hole event horizon, known as "excision" in numerical relativity, to limit the accumulation of energy and density, which otherwise could numerically diffuse out of the event horizon when accumulated to too high values.
In our case, we introduce a ceiling on density (ρ max = 6) as well as pressure (P max = 2) when r < r crit where our critical radius is set to be r crit = 1 r g .Given the location of the critical radius, we have four cells between the event horizon and the critical radius.

Synthetic polarized images
To generate synthetic polarized images of the accretion flow, we use the general-relativistic ray tracing code RAPTOR.RAPTOR solves the polarized radiative transport equations along null geodesics.The geodesic equation is solved starting from a virtual camera outside the simulation domain (r cam = 10 4 r g ).We employ an adaptive camera as described in Davelaar & Haiman (2022).The adaptive camera allows a varying resolution over the image plane, resulting in a computational benefit since most of the resolution can be focused on the event horizon scale, which shows small-scale emission varying on short timescales, while the larger scale can be fully resolved with a lower resolution.We use a base resolution of 625 2 pixels and double the resolution four times, within 60 r g , 40 r g , 20 r g , and 10 r g , respectively, bringing the effective resolution to 10, 000 2 pixels 1 .We use an adaptive Runga-Kutta-Fehlberg method to integrate the geodesic equation and a fourth-order finite difference method to compute the metric derivatives needed for the Christoffel coefficients.The stepsize in RAPTOR is estimated based on the wavevector in the previous step; see Davelaar et al. (2019).For this work, we developed an additional stepsize criterion based on a Courant-Friedrichs-Lewy (CFL) condition, where we require a minimum of eight steps per cell to ensure convergence of all Stokes parameters.We compute synthetic images between 80 GHz and 100 GHz with a frequency spacing of 3 GHz.We also compute time-averaged spectra from the image-integrated total intensity at 25 frequencies with a logarithmic spacing between 10 10 − 10 15 Hz.
The ideal GRMHD equations are dimensionless and do not describe the evolution and thermodynamics of electrons.To this end, we need to introduce a mass scaling and a prescription for the electron temperature.To scale our simulation to a specific black hole, we use a unit of length L = r g , a unit of time T = r g /c, and a unit of mass M. The unit of mass is related to the mass accretion rate via Ṁ = Ṁsim M/T , where Ṁsim is the accretion rate in simulation units.Combinations of these units are then used to scale all relevant fluid quantities.The density is scaled by ρ 0 = M/L 3 , the internal energy by u 0 = ρ 0 c 2 , and the magnetic fields by B 0 = c √ 4πρ 0 (where B 0 is expressed in Gaussian units).Since the black hole mass is often constrained observationally, the only free parameter in our system is M, which can be used to set the energy budget of the simulation such that the total emission produced equals a user-set target.We follow Mościbrodzka et al. (2016), to parameterize the ratio between electron temperature T e and proton temperatures T p via where U the internal energy, m p the proton mass, m e the electron mass, and Θ e the dimensionless electron temperature.The variable R is a free parameter that sets the temperature ratio in regions where β p ≫ 1, allowing for different emission morphology depending on the choice of R, e.g., disk dominated if R = 1 or jet dominated if R ≫ 1.Here, we will limit ourselves to R = 20, e.g., a more jet-dominated model.Note that MADs are relatively insensitive to the exact value of R given that most of the emission originates from a region where β p ≲ 1.Finally, we must choose the electron distribution function's shape and spatial variation.We consider two models, one where the distribution function (DF) is limited to a thermal relativistic Maxwell-Jüttner DF (MJ-DF), and one where we combine the MJ-DF with a κ-DF.The κ-DF deviates from an MJ-DF by having a thermal core and a power law at high Lorentz factors.The κ-DF (Xiao 2006) is given by, where the free parameters are w, which sets the width of the distribution, and κ, which sets the slope of the power law, and N is a normalization constant such that ∞ 1 dne dγ dγ = n e .The κ parameter is related to the power-law index p of the non-thermal tail of the DF via κ = p + 1, such that for γ ≫ 1, dne dγ ∝ γ 1−κ .For the width w, we follow Davelaar et al. (2019) by enforcing that the energy in the DF equals the total available internal energy of the electrons, where Θ e is computed with equation 2. Note that this formula requires κ > 3 (p > 2).
We then compute emission coefficients c ν (emission, absorption, and rotation coefficients), using a prescription introduced in Event Horizon Telescope Collaboration et al. (2022) that combines thermal and κ coefficients via, here σ 0 sets the transition point for the magnetization, which we set to σ 0 = 0.5.The function η(β p , σ) smoothly transitions from thermal to non-thermal component if β p < 1 and σ/σ 0 > 1, representing sites with a large reservoir of magnetic energy that can dissipate to accelerate particles, e.g., in the jet's shear layer.The polarized radiation transfer coefficients are computed via fit formula.For the thermal distribution function, we use the emission and absorption coefficients presented in Dexter (2016) and the rotativities from Shcherbakov (2008), for the κ distribution function, we follow Pandya et al. (2016); Marszewski et al. (2021).
As an archetype LLAGN, we use model parameters consistent with M87*, meaning a black hole mass of M = 6.5 × 10 9 M ⊙ , and a distance of d = 16.8Mpc (Event Horizon Telescope Collaboration et al. 2019a).We set the angle between the BH spin axis and the observer to i = 160 • , following Walker et al. (2018).We set M such that the average flux at 86 GHz is F 86GHz = 1 Jy, as observed by (EHT MWL Science Working Group et al. 2021).The spectrum obtained by EHT MWL Science Working Group et al. ( 2021) also shows a spectral slope in the optically thin part of the spectrum in the near-infrared (NIR) of F ν ∝ ν −1 .Given that for optically thin synchrotron emission, the flux follows F ∝ ν −(p−1)/2 = ν −(κ−2)/2 , to match the observed spectral slope in the NIR, we need to use κ = 4.
To exclude the emission from the spine (interior of the jet, which in GRMHD simulations is typically dominated by artificial density floors), we exclude all the emission if σ > 5. We expect our results to be insensitive to this choice for larger σ values.However, lower σ values would result in a smaller emission region at the jet-wind interface.We also exclude the larger scale disk, x 2 + y 2 > 150r g , which has not settled in a steady state for the runtime of our simulation.However, given the viewing angle of i = 160 • , this choice does not strongly affect our results.

RESULTS
In this section, we summarize our results.In subsection 3.1, we find that surface waves are continuously present along the jet-wind shear layer after our simulation reaches the MAD state.In subsection 3.2, we fit our GRMHD model to the spectrum of M87 and show that the κ-jet model recovers the low-frequency radio and the NIR.In subsection 3.4, we show that the magnetic flux eruptions and waves imprint themselves in various linear polarization quantities, serving as potential tell-tale signatures that could be used to test M87's putative MAD accretion flow state.In subsection 3.5, we highlight that circular polarization plays a minor role.However, we see sign reversals in circular polarization maps that could indicate magnetic field reversals caused by the surface waves.We measure the Faraday rotation measure of our models in 3.6.Lastly, in subsection 3.7, we provide evidence that the polarization signatures are connected to the waves seen in the GRMHD simulation.

GRMHD
To assess when the simulation reaches the MAD state, we compute at the black hole's event horizon: the mass accretion rate, ṁ, and magnetic flux threading the horizon, Φ B , both in dimensionless units.The mass accretion rate is defined as ṁ = r=r h ρu r √ gdθdϕ, where u r is the radial component of the velocity and √ g the determinant of the metric.The magnetic flux is defined as Φ B = 0.5 r=r h | B r | √ γdθdϕ, here B r is the radial component of the magnetic field and √ γ the determinant of the spatial part of the metric2 .We also define the MAD parameter ϕ mad = Φ B / √ ṁ, which was introduced by Tchekhovskoy et al. (2011) to quantify that a simulation reaches the MAD state when ϕ mad ∼ 15.All three quantities, ṁ, Φ B and ϕ mad are shown in Figure 1.The accretion flow reaches for the first time the MAD state at ≈ 3000 r g /c, when the magnetic flux on the horizon saturates as ϕ mad ∼ 15.Our MAD simulation shows globally similar properties to standard simulations on spherical grids; see Appendix A for a comparison.
In Figure 2, we show two-dimensional maps of the logarithm of density sliced along the spin axis, top row, or along the equatorial plane, bottom row.Initially, the jet is more laminar in the left column as flux is still building up on the horizon, and no eruptions have occurred.In the middle column, after the system has reached the MAD limit, flux tubes can be seen in the x-y plane, indicated by lower densities in the disk.The exhaust from magnetic flux eruptions generates these flux tubes.During the eruptions, magnetic energy is dissipated via large equatorial current sheets generated when the disk becomes magnetically arrested, and the northern and southern jets get in direct contact (Ripperda et al. 2022).The exhaust of these sheets forms flux tubes containing vertical magnetic fields that spiral outwards in the disk In the middle panel, small amplitude waves propagate along the shear layer, interfacing the higher-density disk wind and low-density jet in the top panel.Large amplitude waves propagate outwards in the right column because the system produces strong magnetic flux eruptions; see Figure 1 at t = 8800r g /c.The panels correspond to the vertical lines in the middle panel of Fig. 1.The variability introduced by the flux eruptions at the base of the jet acts like a forced oscillator, introducing waves that propagate and grow from near the event horizon to the shear layer between the jet and the disk.The waves propagate to large scales, growing in size while shearing magnetic field lines and generating field reversals, shown in Fig. 3.
To assess the potential effect of the waves on the polarized emission, we show slices along the x-z plane in Figure 4 at t = 8850 r g /c (right panel in Figure 2) of several quantities relevant to the radiation transport.We find that distinct patches of magnetization close to σ ∼ 1 can be seen along the jet-wind shear layer coinciding with the location of the waves in the top right panel in Fig. 1.The higher magnetization is important since particle acceleration is more efficient at higher magnetization values (Sironi & Spitkovsky 2014).We use a passively advected tracer to study the acceleration of wind-based material to relativistic speeds at the jet-wind surface.The tracer is set to zero when the den-Figure 2: Density profiles in the x-z plane (top row) and x-y plane (bottom row, notice different axis scaling compared to the top row).The left column shows the simulation at t = 3200r g /c, well before the magnetic flux is saturated and a MAD state is reached.At this point, the jet-wind shear layer is more laminar.Middle column: at t = 4000r g /c, the accretion flow reaches the MAD state for the first time.In the bottom panel, small under-density regions are visible, indicative of potential flux eruptions, which can be seen as drops in the integrated Φ B time series shown in Figure 1 A first sign of waves can be seen in the jet-wind shear layer.Right column: simulation snapshot at t = 8850r g /c, large under-densities in the bottom panel is seen near the horizon, making the accretion flow non-axisymmetric.The under-densities correlate with large dissipation events of Φ B ; see again middle panel Fig. 1.Additionally, large-scale waves are present in the jet-wind shear layer.sity or pressure is set to floors and is set to one in the initial torus.We then evolve this quantity as a passive scalar, tracing the advection of disk-based material into regions that were at least once set to floors.We find that matter from the disk (un-floored material) around the location of the jet-wind shear layer waves is accelerated to high bulk Lorentz factor, and emission generated in the waves is emitted from on un-floored matter originating from the accretion disk, shown by non-zero tracer and Γ ≥ 2 in the right-top panel in Fig. 4. Looking at the bottom left panel, the waves also correlate with regions with a large fraction of non-thermal electrons since η ∼ 1, which is an evident result of the high magnetization (and low plasma-β), and our choice of electron distribution model (Eqn.6).Finally, the electrons are relativistically hot (i.e., Θ e > 1), as shown in the bottom right panel.We associate the relativistic temperatures and the large fraction of non-thermal electrons with the heating of the plasma by the waves due to the dissipation of magnetic energy, as predicted in Sironi et al. (2021).In summary, the shear layer that is at the surface of the jet interfacing with the wind has a magnetization of order unity, has a high relativistic electron temperature, is moving at relativistic bulk Lorentz factors (Γ ∼ 3), and is likely dominated by non-thermal electrons.

Spectral distribution functions
In Figure 5, we show the spectral distribution functions of our thermal-DF model (thermal-jet) and κ-DF model (κ-jet) models.The top panel shows the total intensity (Stokes I) as a function of frequency.The thermal-jet model recovers the low-frequency part of the spectrum accurately, e.g., ν < 10 12 GHz, however at higher frequencies, it drops off too fast, which is consistent with Davelaar et al. (2019).In the case of the κ-jet model, the high-frequency emission is enhanced and obtains a spectral slope of α ≈ −1, consistent with the observations.In contrast, the thermaljet model underestimates the near-infrared flux.The κ-jet model predicts that non-thermal electrons emit energetic photons pre-dominantly in the jet boundary and are, therefore, a probe of dissipation of magnetic energy due to wave dynamics.To match the observed flux at 86 GHz F 86GHz = 1 Jy we set the mass scaling to M = 1.5×10 25 g for both the thermal and κ-jet models.These units of mass correspond to a mass accretion rate of ≈ 10 −4 M ⊙ /year and a jet power of ≈ 10 42−43 ergs/s, both similar to values obtained in previous works (Chael et al. 2019;Event Horizon Telescope Collaboration et al. 2019b;Cruz-Osorio et al. 2022), and consistent with jet powers inferred from observations of M87 (Prieto et al. 2016).
The two bottom panels show linear polarization (LP) and circular polarization (CP).For LP, the thermal-jet achieves similar fluxes as the κ-jet at a lower frequency (ν ≲ 10 13 Hz), while at a higher frequency, a clear power-law is visible in the κ-jet case.A similar powerlaw is visible for CP, but at a lower frequency, the κ-jet model is comparable to the thermal case.Subsequent  The κ-jet model obtains similar LP and CP as the thermaljet at lower frequencies (ν ≲ 10 13 Hz), while at higher frequencies it produces higher fluxes.
analysis is done at 86 GHz.This frequency was chosen since it probes emission structures at the base of the jet (Hada et al. 2016;Walker et al. 2018).

Total intensity
In the top panels of Figure 6, we show total intensity maps for our κ-jet model at 86 GHz.The images correspond to t = 8500 r g /c and t = 9300 r g /c.At the core of the image, a darkening is visible, corresponding to the "black hole shadow" (Luminet 1979;Falcke et al. 2000), recently observed at 86 GHz by Lu et al. (2023).At the later stage (right panel), when the surface waves are prominently visible in the GRMHD simulations, helical wave-like substructures can be distinguished within the jet at larger scales.

Linear polarization
This subsection summarizes the LP results computed at 86 GHz.In unresolved LP fraction, we find that enhanced emission regions generate loops in a Stokes Q−U diagram during the magnetic flux eruptions.The resolved LP fraction is then inversely proportional to the magnetic flux on the horizon, resulting in an enhancement of the fraction as the flux goes down.Furthermore, we find that the waves seen in the GRMHD simulation imprint themselves in LP maps as they lower the LP fraction.The effect of non-thermal κ-DF on the LP is minor, although we see a slight decrease compared to the thermal model for the LP fraction in the jet.

Core
In Figure 7 top panel, we show a time series of the unresolved LP fraction m net for both the thermal-jet (solid lines) as well as the κ-jet models (dashed lines), defined as The image integrated net LP fraction does not depend on telescope beam size since it is incoherent addition of the Stokes parameters.We obtain an average value of m net = 0.026, consistent with the low values found by observations (Hada et al. 2016;Walker et al. 2018).
The thermal and the κ-jet models show almost identical values.
In Figure 7 bottom panel we show the resolved LP fraction, m, defined as Here we follow the definition of resolved LP fraction from Eqn. 8 of (Event Horizon Telescope Collaboration et al. 2021), which is an image-averaged linear polarization fraction, taking into account some telescope beam size.As a first step we assumed that the telescope fully resolves the image, taking the beam size to be much smaller than the intrinsic features we are interested in.Due to the coherent addition the resulting resolved LP fraction m is substantially higher than the unresolved LP fraction m net .The reason for this is that we preserve the sign of Q and U in the summation, in the case of m net , but the sign is dropped for m.In reality, Q and U will be convolved with the telescope's beam, resulting in incoherent addition.This effect can be seen in the bottom panel of Fig. 7, where we compute the convolved LP fraction m conv for varying telescope beam sizes 20, 40, 80, 160 µas.This computation is done by blurring the original images with a Gaussian filter where the full-width half maximum represents the beam size.As the telescope beam size increases, the underlying substructure in both Q and U is being averaged out, resulting in an overall drop in LP fraction asymptotically approaching m net as the beam becomes comparable to the core size.
In the fully resolved LP fraction m as well as the 20 − 80 µas m conv cases in Fig. 7, at t = 9000 r g /c, an increase in LP fraction is visible.This increase coincides with a flux eruption at the horizon, visible from the over-plotted Φ B curve (solid green line, identical to middle panel of Figure 1).The m shows a fractional increase by 20%, while Φ B shows a fractional decrease of 20%, indicating an inverse relationship between the two quantities.A similar correlation can also be seen at the smaller eruption at t = 8400 r g /c.The κ-jet model reaches identical LP fractions as the thermal-jet.
Close to the flux eruption starting around t ∼ 8700 r g /c, we show the unresolved Stokes parameters in a Q/I − U/I diagram.After the onset of flux eruption, we see a clear clockwise loop with an LP excess of m net ∼ 0.06.The loop we find is similar to the loops found by Marrone et al. (2008); Wielgus et al. (2022) at 230 GHz during X-ray flares for Sagittarius A*.Najafi-Ziyazi et al. ( 2023) finds evidence that these loops could be generated by enhanced polarized emission in the accretion disk by orbiting flux bundles ejected into the disk after a flux eruption.The enhanced emission leads to a local polarization excess: as the emission increases, the polarized emission increases.This enhanced emission region, often called a hot spot, orbits through the local magnetic field.Since the spot only lights up a small region of the accretion disk, this hot spot acts as a probe for the underlying magnetic field geometry.The

Jet
To study the large-scale jet only, we exclude the image plane's inner r g ; in other words, we exclude the near-horizon emission.The resulting m net and m are shown in Figure 9.The unresolved LP fraction m net ranges from 5-20%.The resolved fraction, m, reaches values of around 50%.The κ-jet model shows slightly lower values than the thermal-jet model in the case of m.In the second row of Figure 6, we show a map of m.In these maps, alternating regions of high and low linear polarization are visible, coinciding with enhanced emission in the total intensity panels in the top row.This indicates a potential correlation between LP fraction and the presence of the waves seen in the GRMHD simulation.This potential correlation will be investigated further in Section 3.8.The LP substructure seen in our simulation in Cartesian coordinates is absent in a low-resolution MAD simulation in spherical coordinates at an effective resolution of 192 × 96 × 96 cells in r, θ, ϕ respectively, see Appendix B.
The jet stands out as a high-intensity emission region, while the accretion disk does not contribute to the emis-sion at larger scales (r > 30 r g ).Comparing the total intensity map with the LP fraction map, the foreground disk surrounding the jet shows high values of m (close to unity).To understand this behavior we evaluate the asymptotic limit of our fitting formula for the emission coefficients J S (with S indicating one of the Stokes parameters), see Eqn. 31 in Pandya et al. (2016).Given that the disk has weak magnetic fields and low temperatures; we have ν/ν c ≪ 1, with ν c = eB/(2πm e c 2 ), and Θ e ≪ 1, we find that J I /|J Q | → 1.0.This makes physical sense since due to the low temperature the thermal distribution function is narrow, which means that there is a quasi-mono-energetic population of electrons that is causing the emission so the polarization fraction should go to unity.This, however, does not alter the computation of our image integrated m, m conv and m net , since this outer region has low intensity both polarized as well as Stokes I, so they don't contribute to the numerator and denominator of Eqn 7 and 8. Additionally, we also exclude regions of very low intensity from our map, where we set the Stokes parameters of a pixel to zero if S I /max(S I ) < 10 −6 .
Lastly, we compute polarization angle maps as a function of beam width.These can be seen in Fig. 10, where we overplotted ticks of the polarization vector on top of the total intensity map from the top left in Fig 6 .The length of the ticks is set by the linear polarization fraction, while the angle is set by the EVPA, as defined by χ = 1 2 tan −1 (S Q /S U ), here we also exclude linear polarization fractions when S I /max(S I ) < 10 −6 , so tick lengths are set to zero.For the case of zero beam size the polarization vector clearly follows the ridges of enhanced intensity.As the beam size increases, the correlation becomes weaker, however re-orientation of the polarization vector is also here visible, e.g. at x ∼ 0 r g y ∼ 50 r g the polarization pattern switches from diagonal, to horizontal, and back to diagonal.We only show beam sizes up to 60 µas which is the expected resolution of the ng-EHT at 86 GHz (Issaoun et al. 2023), assuming identical baselines to the current EHT array (using, θ ∼ 20µas(λ/1mm) with λ = 3 mm).

Circular polarization
In this subsection, we study the CP fractions computed at 86 GHz.We find that the resolved CP fraction decreases during a flux eruption as the inner accretion disk is ejected, resulting in a more dilute plasma to perform Faraday conversion.Overall, we find low unresolved and resolved CP fractions for our thermal and κjet models.The CP maps show sign reversals, indicative of alternating magnetic field orientation that coincides with the features seen in the linear polarization maps.In Figure 11 we show the unresolved and resolved CP fractions, v net and v, defined as Both v net and v are small in value, as expected from synchrotron radiation (Rybicki & Lightman 1979).In Figure 11   version converts linear to circular polarization.Due to the ejection of the inner part of the accretion disk during a flux eruption, there is more dilute plasma in this region, resulting in a drop in the CP fraction.

Jet
In Figure 12, we compute v net and v, but now also exclude the inner 30r g to exclude the near horizon emission and focus on the larger scale jet only.This exclusion results in smaller fractions than the entire imageintegrated values since Faraday conversion happens in high-density regions.The third row of Fig. 6 shows maps of v where the CP fraction is smaller at larger radii.In the map, we preserved the sign of Stokes V, which is set by the direction of the magnetic field along the line of sight.The image shows reversals in the sign of v and additional ridges of low CP fraction.Since Stokes V carries information on the direction of the magnetic field along the line of sight, this indicates a potential orientation switch in the underlying magnetic field geometry.This reversal could be caused by the waves, as shown in Figure 3; we will further investigate the sign reversal in Section 3.8.

Faraday rotation
Figure 6 bottom row shows maps of rotation measures computed between 80 and 100 GHz.The rotation measure is defined as where χ ν is the electric vector position angle (EVPA) computed at a specific frequency ν defined as χ = 1 2 tan −1 (S Q /S U ), and λ the wavelength.Large RMs are visible along the jet, and in the core, we find values as large as 10 4 − 10 5 rad/m 2 .However, the core dominates the total RM, which is expected since the Faraday depth is larger due to higher density and lower temperatures.The relatively large value for the RM in the jet is somewhat surprising, given that the jet does not exhibit large Faraday depths.Given the small Faraday depth, the change in EVPA is not caused by Faraday rotation but is caused by the transverse gradients in n e , Θ e and B in the emitting shear layer.These gradients result in emission at different frequencies to peak at different depths in the shear layer, which have different plasma properties, e.g., different magnetic field orientations will result in different orientations of the EVPA, giving rise to non-zero RM values.

Properties along a ray
To test if the waves cause linear depolarization, we identified representative light rays showing high or low polarization fractions.The selected geodesics are indicated with the red, blue, and green dots in the top panel of Figure 13.We then compute the linear polarization fraction as a function of the Cartesian Kerr-Schild x ′ coordinate, meaning smaller values of x ′ are closer to the spin axis.The result of this is shown in the bottom of Figure 13.We show the geodesics only as they approach the jet-wind surface and only show segments when the local magnetization σ < 5, meaning no radiation transfer is applied when the geodesic is inside the jet.The red-colored ray is terminated early, meaning that for larger x ′ its net polarization fraction is higher, indicating that the shear layer at that point is thinner.The blue and green-colored rays have a larger travel path, meaning the polarization starts to drop.In the case of the green-colored ray, the situation is even more interesting.The line of this ray is interrupted twice, indicating it crossed into two regions of high magnetization but then left.We interpret this as the ray crossing through a rolling wave, similar to the wave seden at y = −100 r g in Figure 2. We test this by also computing the magnetization along the ray, this is shown also in the bottom panel of Fig. 13 by the black line.This line crosses our magnetization threshold (σ > 5) twice.In general, the waves alter the thickness of the jet-wind shear layer, and their presence results in varying path lengths inside the shear layer for different rays, which affects the linear polarization fraction.

Correlation with the jet-wind shear waves
To connect the waves in GRMHD with the structures seen in the synthetic 86 GHz images, we compute emissivity weighted averages of the magnetization σ, the pitch angle between the wave vector and the magnetic where q represents the weighted quantity, the result of this computation is shown in Figure 14.The top left panel, σ, shows that the images' higher intensity features also have a larger magnetization.This agrees with the waves having a larger magnetization in Figure 4. Additionally, the same patterns are visible in the higher electron temperatures (top right panel in Figure 14), and the over-densities (bottom left panel in Figure 14), also in agreement with the properties of the waves, as shown in Section 3.1.
For Stokes V, we finally compare the emissivity weighted average of cos(θ B ) (bottom right panel in Figure 14).The overall map of cos(θ B ), where θ B is the angle between the wave vector and the magnetic field vector, shows the same sign as Stokes V, implying that reversals in Stokes V, are caused by the shearing of the magnetic fields leading to reversed field orientations.

DISCUSSION AND CONCLUSION
In this Letter, we present a global GRMHD simulation in Cartesian Kerr-Schild coordinates in the MAD regime that shows the formation of waves along the jetwind shear layer.We post-process this simulation with our polarized radiative transfer code and compute polarized spectral energy distributions, times series of polar-ization quantities, and synthetic images at 86 GHz.We find observational signatures of the surface waves seen in the GRMHD simulation in the polarization information.As the waves propagate outwards, they show up as bright features along the jet that alter the polarization signature of the jet at larger scales.As the waves alter the orientation of the magnetic field lines, the linear polarization fraction drops due to the cancellation of subsequently rotated Stokes vectors.
During magnetic flux eruptions, we find an inversed relation between Φ B and the LP fraction, meaning that as Φ B drops, the LP fraction increases.We see the opposite for the CP, where the fraction decreases as Φ B decreases.Both effects can be explained by the flux eruption ejecting the disk near the horizon; as the strong poloidal field arrests parts of the disk, the density drops, and the disk becomes more optically and Faraday thin, leading to lower Faraday rotation and conversion.This results in a decrease in CP and an increased LP fraction.
At the largest magnetic flux eruption in our simulation, we find that the unresolved Stokes Q and U at 3 mm show a clockwise loop in a Q − U diagram with a polarization excess of m net ∼ 0.06.Loops like this were previously identified in the case of our galactic supermassive black hole SgrA*, either observationally (Marrone et al. 2008;Wielgus et al. 2022;The GRAVITY Collaboration et al. 2023), or in theoretical works, e.g.(Vos et al. 2022;Najafi-Ziyazi et al. 2023).A key difference is that the time scales in the case of M87 are longer, which puts the period of our loop at ∼ 3 months, compared to ∼ 1 hour in the case of SgrA*.
Our model recovers resolved polarization fractions that are too high compared to the ones measured by Hada et al. (2016).However, when convolved with a more realistic telescope beam, we show that the LP fraction substantially drops since Q and U are averaged out due to patches within the beam having opposing signs.We find consistent rotation measures without invoking an external Faraday screen, and we recover the observed spectral shape from radio to optical frequencies (EHT MWL Science Working Group et al. 2021).Although we limit ourselves to M87, our results generally apply to other LLAGNs reaching the MAD state since the waves result from the underlying flow geometry and the flux eruptions typical for such systems.We expect these polarization signatures to be independent of black hole masses and accretion rates.
In the literature, studies of wave instability at jetwind shear layers are typically limited to analytical studies, e.g., linear analysis (Ferrari et al. 1978;Sobacchi & Lyubarsky 2018;Chow et al. 2022) or with numerical MHD/Particle-in-Cell studies of local idealized se-tups (Hardee et al. 2007;Perucho et al. 2010;Sironi et al. 2021).The overall conclusions of these works are that jets, if in the right conditions, can be prone to the excitation of waves due to linear instability, e.g., KH waves.These waves are asymmetric, meaning they have different plasma properties on either side of the shear layer.Previous work by Sironi et al. (2021) showed that particles can be accelerated to high energies in mildly relativistic, magnetized asymmetric shear flows.However, evidence of wave instabilities in global simulations is sparse and often underresolved in 3D simulations due to the restrictions on the large-scale resolution in spherical coordinates; see Chatterjee et al. (2019); Wong et al. (2021).Observationally, some evidence for wave-like perturbations at large distances from the central engine is found by, e.g., Perucho & Lobanov (2007);Pasetto et al. (2021); Issaoun et al. (2022).
In this work, we did not perform a rigorous linear analysis to determine if the waves could be grown from linear scales and what instability is driving them.Visually, the jet is initially stable and shows no waves, while when the system reaches the MAD state and the first flux eruptions occur, waves travel outwards along the jet-wind boundary.Therefore, the waves we see are more likely to grow by forced oscillations of the jet base due to accretion variability and ejecta from magnetic flux eruptions.The waves are already non-linear within a few gravitational radii, which would require short linear growth times.A more likely scenario is that the jet base's variability efficiently drives the waves' growth and becomes non-linear at larger scales.A study of the conditions under which these waves are growing by either applying linear analysis (Chow et al. 2022) to local conditions extracted from our simulation or by performing local idealized simulations will be done in future works.
Compared to previous global simulations of LLAGN jets, our simulations stand out due to the Cartesian nature of the grid, allowing us to resolve the jet to larger distances compared to more standard simulations with spherical grids as used in, e.g., Event Horizon Telescope Collaboration et al. (2019c).This higher resolution enables us to follow the perturbations at the jet base to larger scales.However, the waves we see in our simulations are not a result of our choice of coordinates and can be found in spherical simulations if run at sufficiently high resolution (Ripperda et al. 2022), as shown in Appendix A.
The evidence we find for the shear layer waves in our simulation may have implications for particle energization.The waves could introduce a source of turbulence or reconnection in the shear layer (Sironi et al. 2021).These processes could lead to electron acceleration, re-sulting in non-thermal emission that could explain the edge brightening of AGN jets.Additionally, reconnection induced by the waves could drive injection of high energy ions originating from the disk into shear-driven acceleration, potentially producing ultra high energy Cosmic Rays, see, e.g., Caprioli (2015); Rieger (2019); Mbarek & Caprioli (2021).
In this work, we discovered a correlation between jetwind surface waves and polarized emission properties.We find evidence that the substructure in the jet, in the form of waves, imprints itself on the Stokes I and LP maps.We identify ridges and alternating low and high linear polarization fractions as tell-tale signatures of these waves.Although currently below the achievable resolution of VLBI arrays, this effect might be resolvable by future next-generate arrays such as the nextgeneration EHT (ng-EHT) (Ricarte et al. 2023;Issaoun et al. 2023).If the ng-EHT operates at 86 GHz, it will achieve a resolution of ∼ 60 µas, or 20 r g scaled to M87, which would be sufficient for resolving the features we find in this study.We cross-compared our CKS MAD simulation with a low-resolution spherical MKS simulation to further strengthen our conclusions.The MKS simulation has a base resolution of [96,48,48] in r, θ, ϕ, and one additional level of AMR.The simulation was run up to t = 10, 000 r g /c with BHAC.Due to the low resolution in this simulation's jet region, no waves are present along the jet-wind surface.Note that due to the low resolution, the physical solution of this simulation is far from resolved and, therefore, in an unrealistically low regime of Reynolds number.We only use it here to compare a laminar flow to a flow where the jet-wind surface shows wave instabilities.We ray trace the spherical simulation over the final 2000 r g /c, with the same model and camera parameters as the κ-jet model presented in the main manuscript.Comparing the resolved linear polarization fraction with the higher resolution Cartesian case shows a substantially higher fraction, namely at m ∼ 0.7, compared to ∼ 0.5, see Figure 17

Figure 1 :
Figure 1: Time series in dimensionless units of horizon integrated mass accretion rate ṁ (top panel), magnetic flux Φ B (middle panel, gray lines correspond to the slices shown in Figure 2), and the MAD parameter ϕ mad = Φ B / √ ṁ (bottom panel).The MAD parameter saturates at ϕ mad ≈ 15, corresponding to the horizontal black line.

Figure 3 :
Figure 3: The density profile in the x-z plane overplotted with the magnetic field lines in the same plane.The magnetic field lines are sheared in the jet-wind shear layer, and roll-ups and magnetic islands are visible.

Figure 4 :
Figure 4: Simulation snapshot at t = 8850 r g /c slice along the x-z plane of various quantities used for the radiative transfer calculation.Top left: cold magnetization σ = b 2 /ρ.The inner core of the jet, which is at σ ≥ σ cut = 5, is dominated by the simulation floors, also visible in Figure 2 by the low-density regions in the jet.Top right: the specific kinetic energy Γ − 1 multiplied by our passive scalar.The scalar is set to zero for the floors and unity for disk matter.At the jet-wind shear layer, we see an increase in Γ − 1 around the locations of the waves, indicating that matter originating from the disk is mixed in the shear layer and accelerated to high bulk Lorentz factor.Bottom left: acceleration efficiency η, where η = 1 means κ-DF only, whereas η = 0 is thermal-DF only.The jet-wind shear layer shows η = 1, coinciding with the high temperatures in the right bottom panel.Grey region indicates σ ≥ 5, which is excluded from our GRRT computations.Bottom right: the electron temperature Θ e as prescribed by Eqn.1b.The largest temperatures are found in the jet-wind shear layer.Grey region indicates σ > 5.

Figure 5 :
Figure 5: Spectra of time-averaged Stokes I (top panel), linearly polarized (LP) flux (middle panel), and circular polarized (CP) flux (bottom panel) for both the κ-jet (cyan) as well as thermal-jet (purple) models.The shaded area shows a standard deviation on the timeaveraged fluxes.The thermal-jet under-produces the NIR emission compared to observations for Stokes I, while the κ-jet recovers the observed spectral slope.The κ-jet model obtains similar LP and CP as the thermaljet at lower frequencies (ν ≲ 10 13 Hz), while at higher frequencies it produces higher fluxes.

Figure 6 :
Figure 6: Synthetic synchrotron images at 86 GHz of two GRMHD simulation snapshots at t = 8500r g /c and t = 9000r g /c.Top row: total intensity.Second row: LP fraction m, Eqn 8. Third row: CP fraction v, Eqn. 9. Bottom row: Rotation Measure, Eqn 11, normalized by 10 5 rad/m 2 .Large values of linear polarization at the edges of the jet are artificial since they are caused by regions of very low total intensity.

Figure 7 :
Figure 7: LP fraction as a function of time.Top panel, net LP fraction m net , for the thermal and κ-jet models, showing an average m net ∼ 0.026.Middle panel:resolved LP fraction m, overplotted with the magnetic flux on the horizon Φ B (green curve).Both models show a substantially higher resolved LP fraction.Additionally, at the moment of a magnetic flux eruption, t = 9000 r g /c, an increase in m is visible.Bottom panel: LP fraction as function as telescope beam size.The LP fraction decreases with increasing LP fraction due to incoherent addition.

Figure 8 :
Figure 8: Q − U diagram at 86 GHz.At the strongest magnetic flux eruption during the duration of our simulation at t ∼ 9000r g /c, we observe a pattern that manifests as a loop moving in a clockwise motion in a Q − U diagram with a linear polarization excess of m net ∼ 0.06.

Figure 9 :
Figure 9: Identical to Figure 7 but now with the inner 30 r g excised from the image plane of m net , m and m conv to compute the jet contribution only.Also for the excised LP fraction, both models obtain similar fractions.

Figure 10 :
Figure 10: Polarization angle maps as a function of beam size, overplotted on the total intensity map of the top left panel in Fig.6.Top left, zero beam size, the polarization vector shows a clear correlation with the ridges seen in total intensity.For increasing beam size, this correlation becomes weaker, however, some reorientation of the vector is still visible, see e.g.x ∼ 50 r g , y ∼ 50 r g .

Figure 11 :
Figure 11: CP fraction as a function of time.Top panel, net CP fraction v net , for the thermal and κ-jet models, showing an average v net ∼ 0.026.Middle panel: resolved CP fraction v, overplotted with the magnetic flux on the horizon Φ B (green curve).The thermal-jet model and κ-jet model obtain similar CP fractions.Additionally, at the moment of a magnetic flux eruption, t = 9000 r g /c, a decrease in v is visible.Bottom panel: cP fraction as function as telescope beam size.The CP fraction decreases with increasing beam size due to incoherent addition.
bottom panel, during the strongest flux eruption at t = 9000 r g /c, both models show a slight decrease in CP fraction.The accretion disk enhances the amount of circularly polarized emission as Faraday con-

Figure 12 :
Figure 12: Identical to Figure 11 but now with the inner 30 r g excised from the computation of v net , v and v conv.While for vnet, both models are similar, for v, the κ-jet model similar CP fractions as the thermal-jet model.

Figure 13 :
Figure 13: In the top panel we show a linear polarization map, m.Red, blue, and green dots indicate the pixels along which we compute the linear polarization as a function of the geodesics x ′ .The dependence on x ′ is shown in corresponding colors in the bottom panel.The black line indicate magnetization along the pixel corresponding to the green dot (axis on the right).

Figure 14 :
Figure 14: Emissivity weighted averages of the magnetization σ (top left), electron temperature Θ e (top right), electron number density n e (bottom left), and the angle between the magnetic field and wave vector cos(θ B ) (bottom right).

Figure 16 :
Figure16: Simulation snapshots at t = 10, 000 r g /c for the four HAMR MKS models at varying resolution.Shown are slices along the x − z plane of the logarithm of density.As the resolution increases from left to right and top to bottom, the simulations show a more extended substructure in the form of waves along the jet-wind interface, where the lowest two resolutions are more diffuse than the highest resolutions.
left panel.Looking at synthetic images of linear polarization, Figure17right panel, also no wave-like substructure, as seen in the Cartesian case, is visible along the jet-wind shear layer.This comparison, therefore, further confirms our hypothesis that the jet-wind shear layer waves, which are only captured with sufficiently high resolution, lead to the drop in LP fraction.

Figure 17 :
Figure 17: Left: LP fraction as a function of time for both the Cartesian and the Spherical simulations.Both panels exclude the core emission within 30 r g of the image origin.Top: net LP fraction m net , showing similar values.Bottom: resolved LP fraction m, showing a clear factor 1.5 difference.Middle: Stokes I map of the Spherical simulation, which shows limited substructure.Right: LP fraction map of Spherical simulation.A limited amount of structure is visible compared to the Cartesian case and shows a substantially higher LP fraction, m ∼ 0.75 versus m ∼ 0.5.