# Resolution analysis of magnetically arrested disk simulations

L. Salas1 G. Musoke2,1 K. Chatterjee3,4 S.B. Markoff1,5 O. Porth1 M.T.P.
Liska6,7 B. Ripperda2,8-11
1Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science
Park 904, 1098 XH Amsterdam, The Netherlands
2Canadian Institute for Theoretical Astrophysics, University of Toronto, 60
St. George Street, Toronto, ON M5S 3H8, Canada
3Institute for Research in Electronics and Applied Physics, University of
Maryland, 8279 Paint Branch Drive, College Park, MD 20742, USA
4Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge, MA
02138, USA
5Gravitation and Astroparticle Physics Amsterdam Institute, University of
Amsterdam, Science Park 904, 1098 XH 195 196 Amsterdam, The Netherlands
6Center for Relativistic Astrophysics, Georgia Institute of Technology, Howey
Physics Bldg, 837 State St NW, Atlanta, GA, 30332, USA
7Institute for Theory and Computation, Harvard University, 60 Garden Street,
Cambridge, MA 02138, USA
8Dunlap Institute for Astronomy and Astrophysics, University of Toronto, 50
St. George Street, Toronto, ON M5S 3H4, Canada
9Department of Physics, University of Toronto, 60 St. George Street, Toronto,
ON M5S 1A7, Canada
10Perimeter Institute for Theoretical Physics, 31 Caroline Street North,
Waterloo, ON N2L 2Y5, Canada
11Center for Computational Astrophysics, Flatiron Institute, 162 5th avenue,
New York, NY, 10010, USA
E-mail<EMAIL_ADDRESS>
(Accepted 2024 July 25. Received 2024 July 19; in original form 2024 April 19)

###### Abstract

Polarisation measurements by the Event Horizon Telescope from M87∗ and Sgr
A∗suggest that there is a dynamically strong, ordered magnetic field, typical
of what is expected of a magnetically arrested accretion disk (MAD). In such
disks the strong poloidal magnetic field can suppress the accretion flow and
cause episodic flux eruptions. Recent work shows that General Relativistic
Magnetohydrodynamic (GRMHD) MAD simulations feature dynamics of turbulence and
mixing instabilities that are becoming resolved at higher resolutions. We
perform a convergence study of MAD states exceeding the status quo by an order
of magnitude in resolution. We use existing 3D simulations performed with the
H-AMR code, up to resolution of 5376 × 2304 × 2304 in a logarithmic spherical-
polar grid. We find consistent time-averaged disk properties across all
resolutions. However, higher resolutions reveal signs of inward angular
momentum transport attributed to turbulent convection, particularly evident
when mixing instabilities occur at the surfaces of flux tubes during flux
eruptions. Additionally, we see wave-like features in the jet sheath, which
become more prominent at higher resolutions, that may induce mixing between
jet and disk. At higher resolutions, we observe the sheath to be thinner,
resulting in increased temperature, reduced magnetisation, and greater
variability. Those differences could affect the dissipation of energy, that
would eventually result in distinct observable radiative emission from high-
resolution simulations. With higher resolutions, we can delve into crucial
questions about horizon-scale physics and its impact on the dynamics and
emission properties of larger-scale jets.

###### keywords:

accretion, accretion disks – black hole physics – MHD

††pubyear: 2024††pagerange: Resolution analysis of magnetically arrested disk
simulations–Resolution analysis of magnetically arrested disk simulations

##  1 Introduction

The accretion of matter is a fundamental astrophysical process occurring
across a variety of compact objects such as neutron stars, white dwarfs, and
black holes (BHs). Jets and wind outflows are a natural consequence of the
accretion process and are present across all scales of black hole mass– from
black hole X-ray binaries to the supermassive black holes powering the bright
compact regions at the centres of most massive galaxies (active galactic
nuclei; AGN). Accretion disks and outflows emit radiation across the
electromagnetic spectrum, from radio frequencies to high energy gamma rays
(e.g. ; , ). The relativistic outflows (jets) launched by AGN in particular
are thought to have a profound effect on their environments, playing a key
role in regulating star formation, galaxy evolution and the evolution of
galaxy clusters (e.g. Silk & Rees, 1998; Schawinski et al., 2007; Fabian,
2012). Numerical simulations have improved our understanding of the accretion
process and the mechanisms by which jets are launched and affect their
environments. Yet it remains unclear which resolutions are sufficient to
accurately capture the global evolution of accreting black holes, and jets.

Theoretical modelling and numerical simulations, combined with multi-
wavelength observations have provided critical insights into the nature of
outflows (jets) and radiation from accreting black hole systems (e.g. ; ; , ).
In particular, GRMHD simulations have become increasingly successful in
explaining the dynamics of the accretion flow, winds, and the evolution of
relativistic jets close to the black hole, therefore playing a fundamental
role in the interpretation of observational data. GRMHD simulations have been
significantly improving over the last few decades owing to the increase in
computing power and the development of efficient numerical algorithms in
multiple GRMHD codes (e.g. ; ; ; ; ; ; ; ; ; , ). Porth et al. (2019)
demonstrated the general robustness of GRMHD simulation results regardless of
the different algorithms, implementations, and grid geometries of different
numerical codes in the weakly magnetized regime (SANE). Recently, there has
been significantly more interest in a category of numerical solutions wherein
strong poloidal magnetic fields can interrupt the accretion flow, forming a
magnetically arrested disk (MAD; Bisnovatyi-Kogan & Ruzmaikin, 1974; Narayan
et al., 2003). White et al. (2019) conducted a resolution study of MAD
accretion flows, with a maximum effective 173 cells per decade in radius, 256
cells in polar angle and 512 cells in azimuthal angle. They showed that the
general large-scale structure of the flow is robust with resolution. However,
the spatial structure and Lorentz factor of the jet were not fully converged.

In recent years, the unprecedented resolution of the Event Horizon Telescope
(EHT) has enabled the first direct imaging of the plasma surrounding the
supermassive black holes M87∗ and Sgr A∗(; ). The EHT observations revealed a
bright ring, formed by the hot luminous plasma close to the black hole, along
with a dark central region, the black hole ’shadow’. Synthetic images produced
by ray-tracing of GRMHD simulations are the fiducial models used by the EHT
collaboration to compare against observations. The closest-matching GRMHD
simulations to the EHT observations suggest that both M87∗ and Sgr A∗exhibit
accretion flows in the MAD regime (e.g. ; ; ; , ). Moreover, Sgr A∗can reach
the MAD state in the inner region when poloidal magnetic fields are advected
inwards from larger scales, as demonstrated by wind-fed accretion simulations
(Ressler et al., 2023). Additionally, Liska et al. (2020) showcased that
turbulence within the accretion disk itself can generate strong poloidal
magnetic fields in-situ.

Recent work shows that MAD states feature complicated dynamics that become
resolved at higher resolutions. Consequently, simulating MAD states is
numerically challenging. Episodic flux eruptions result in the ejection of
large cavities containing flux tubes with low-density through magnetic
reconnection (Ripperda et al., 2020; Ripperda et al., 2022), which may serve
as a mechanism for flare generation (Dexter et al., 2020; Porth et al., 2021;
Ripperda et al., 2022; Hakobyan et al., 2023). The Rayleigh-Taylor instability
(RTI) mechanism acting on the boundary of a flux tube can lead to secondary
instabilities, turbulence, and magnetic reconnection (Zhdankin et al., 2023).
Ripperda et al. (2022) captured plasmoid-mediated magnetic reconnection in a
3D MAD state for the first time, by conducting a GRMHD simulation with maximum
effective resolution of 1630 cells per decade in radius, 2304 cells in polar
angle and 2304 cells in azimuthal angle. Plasmoid-mediated magnetic
reconnection can result in particle energisation through dissipation and may
explain flares from accreting black holes, in particular the TeV flares
observed from M87∗ (Aharonian et al., 2006; Acciari et al., 2010; Aliu et al.,
2012; Blanch, 2021), and the infrared and X-ray flares thought to originate
close to the event horizon of Sgr A∗(Baganoff et al., 2001; Eckart et al.,
2004; Neilsen et al., 2015; GRAVITY Collaboration et al., 2018; Gravity
Collaboration et al., 2021).

The current simulations conducted with the Graphics Processing Unit
(GPU)-accelerated GRMHD code H-AMR (Liska et al., 2022) are pushing the
boundaries of what is computationally feasible with regards to resolution.
However, conducting high resolution GRMHD simulations is computationally
extremely expensive. In this paper, we perform a convergence study on a wider
range of resolutions than has ever been previously explored (see Table 1),
employing existing 3D ideal GRMHD simulations with H-AMR of an accreting black
hole in the MAD regime (Ripperda et al., 2022). We determine how well the
physical properties of the accretion flow at different spatial resolutions
agree. We identify the processes that are only recovered at high resolutions.
We constrain the minimum resolution that ensures consistent and converged
results, which will enable the GRMHD community to minimise the computational
cost of these simulations, without sacrificing the accuracy of the global
evolution of the disk-jet-black hole system.

The structure of this paper is as follows: Section 2 provides the numerical
setup for simulations and an overview of the H-AMR code. The impact of
resolution on the disk and the jet is analysed in Section 3. The conclusions
of the study are presented in Section 4.

##  2 Simulations

###  2.1 GRMHD equations

The equations of GRMHD comprise the particle number conservation equation:

$\left(nu^{\mu}\right)_{;\mu}=0\text{ ,}$ (1)

the energy momentum conservation equations:

${T^{\mu}}_{\nu;\mu}=0\text{ ,}$ (2)

and the Maxwell’s equations:

${F^{\ast\mu\nu}}_{;\nu}=0\text{ ,}$ (3)

In ideal MHD, the dual of the electromagnetic field tensor is
$F^{\ast\mu\nu}=b^{\mu}u^{\nu}-b^{\nu}u^{\mu}$ and the stress-energy tensor
is:

${T^{\mu}}_{\nu}\equiv\left(\rho+u_{g}+p_{g}+b^{2}\right)u^{\mu}u_{\nu}-b^{\mu}b_{\nu}+\left(\frac{1}{2}b^{2}+p_{g}\right){\delta^{\mu}}_{\nu}\text{
,}$ (4)

Here $n$ is the particle number density, $\rho=mn$ is the rest mass
density111Normalised to a maximum density $\rho_{\textup{max}}=1$., $m$ is the
mean rest mass per particle, $u_{g}$ is the gas energy density, $p_{g}$ is the
gas pressure, and $b^{2}/2$ is the magnetic energy density. The 4-velocity and
the magnetic field 4-vector are $u^{\mu}$ and $b^{\mu}$ respectively.
${\delta^{\mu}}_{\nu}$ is the Kronecker delta. The gas pressure is
proportional to the fluid-frame rest-mass density and proton temperature
($p_{g}\propto\rho T_{p}$) and the magnetic pressure $p_{b}=b^{2}/2$, where
$b$ is the magnetic field strength in the frame co-moving with the fluid. The
proton temperature ($T_{p}$), in CGS units:

$T_{p}=\frac{m_{p}c^{2}p_{g}}{k_{B}\rho}\text{ ,}$ (5)

where $m_{p}$ is the proton mass and $k_{B}$ is the Boltzmann constant. We use
geometrized units with gravitational constant, black hole mass, and speed of
light $G=M=c=1$, and a factor of $1/\sqrt{4\pi}$ is absorbed in the definition
of the magnetic field. Greek indices run through $[0,1,2,3]$ and Roman indices
through $[1,2,3]$. The metric determinant is $g$.

###  2.2 H-AMR implementation

H-AMR (Liska et al., 2022), a state-of-the-art 3D GRMHD code, builds upon the
foundation of the original HARM2D code (Gammie et al., 2003). H-AMR converts
conserved quantities, such as particle number density and energy-momentum
density, into primitive variables, including rest-mass density, internal
energy density, and velocity components, utilising the inversion scheme by
Noble et al. (2006). The Adaptive mesh refinement (AMR), local adaptive time-
stepping (LAT), and GPU acceleration within a hybrid CUDA-OpenMP-MPI framework
utilised in H-AMR speed-up computations. The achieved resolutions are $\approx
10\times 10\times 5$ times higher than those in the previous MAD study (White
et al., 2019).

H-AMR utilises a staggered grid for constrained transport of magnetic fields
as described in (Gardiner & Stone, 2005), and solves the GRMHD equations of
motion in conservative form in arbitrary (fixed) spacetimes. H-AMR uses a
finite volume, shock-capturing Godunov-based HLLE scheme, with third order
accurate spatial reconstruction of cell variables (PPM, Colella & Woodward,
1984) on cell faces and second order accurate time evolution. The simulations
are performed on a logarithmic spherical-polar grid in a Kerr-Schild
foliation, $\log(r)$, $\theta$ and $\phi$. Since cells get squeezed near the
pole, the timestep in all spherical grids is reduced by an additional factor
proportional to the resolution in the $\phi$-direction. To remedy this issue,
multiple levels of static mesh de-refinement in the $\phi$-direction are used
to keep the aspect ratio of the cells close to uniform at high latitudes. This
method prevents the squeezing of cells near the pole from reducing the global
timestep, while maintaining high accuracy in all three dimensions (see section
3.4 in Liska et al., 2022). The simulations presented in this work utilise
outflow boundary conditions in the radial direction, transmissive boundary
conditions in the $\theta$-direction and periodic boundary conditions in the
$\phi$-direction.

###  2.3 Numerical setup

In order to study convergence in GRMHD simulations of MAD accretion flows
around a spinning black hole, we analyse a total of five 3D GRMHD simulations.
The simulation with the highest resolution, R5, is performed with an
effective222To prevent the squeezing of cells, three internal and four
external derefinement levels in $\phi$ are used to reduce the resolution from
$N_{\phi}=128-2304$ at $30^{\circ}<\theta<150^{\circ}$ to $N_{\phi}=16-18$
within $0.5^{\circ}-7.5^{\circ}$ of each pole (Liska et al., 2022). resolution
of $5376\times 2304\times 2304$ and the other four simulations we present are
a factor $\approx 2\times 2\times 2$ less in resolution (see Table 1). The
radial domain of the simulation grid is $r=\left[1.2-2000\right]r_{g}$. Each
simulation is evolved for at least $10^{4}r_{g}/c$ (see fifth column of Table
1). In all runs the disk is initialised using a torus in hydrostatic
equilibrium (Fishbone & Moncrief, 1976) around a Kerr black hole with
dimensionless spin $a=0.9375$. The inner edge of the torus is located at
$r=20r_{g}$ and the pressure maximum at $r=41r_{g}$. The torus is threaded
with a single poloidal magnetic field loop, defined by the $\phi$-component of
the vector potential
$A_{\phi}\propto\textup{max}\left[\rho/\rho_{\textup{max}}\left(r/r_{\textup{in}}\right)^{3}\sin^{3}\theta\exp\left(-r/400\right)-0.2,0\right]$,
and normalised to obtain a gas-to-magnetic-pressure ratio
$p_{g,\mathrm{max}}/p_{b,\mathrm{max}}=100$. We adopt an equation of state for
an ideal gas with an adiabatic index of $\gamma_{\mathrm{ad}}=13/9$, which
corresponds approximately to a gas mixture of relativistic electrons
$\gamma_{\mathrm{ad},e}=4/3$ and non-relativistic ions
$\gamma_{\mathrm{ad},p}=5/3$. A semi-relativistic gas with
$\gamma_{\mathrm{ad}}\sim 1.55$ was found for BH accreting at $\sim 10^{-6}$
the Eddington accretion rate, in radiative two-temperature GRMHD simulations
(Liska et al., 2024). Within the spine region the following floor and ceiling
values are employed; the rest-mass density floor is
$\rho_{\textup{fl}}=\mathrm{MAX}[b^{2}/25,10^{-7}r^{-2},10^{-20}]$, the gas
energy density floor is
$u_{g,\textup{fl}}=\mathrm{MAX}[b^{2}/750,10^{-9}r^{-26/9},10^{-20}]$ and the
magnetisation ceiling is $\sigma_{\textup{max}}=25$ where $\sigma=b^{2}/\rho$.
The disk is initialised with the same energy density perturbation for all
resolutions, $u_{g}(1+0.04(\textrm{rand([0, 1])}-0.5))$, where rand([0, 1]) is
a random deviate between 0 and 1.

Table 1: Number of cells in the radial $N_{r}$, polar $N_{\theta}$ and azimuth $N_{\phi}$ directions for each simulation, together with the final time of the simulation $t_{\mathrm{sim}}$. Resolution | $N_{r}$ | $N_{\theta}$ | $N_{\phi}$ | $t_{\mathrm{sim}}(r_{g}/c)$
---|---|---|---|---
R1 | 288 | 128 | 128 | 20420
R2 | 580 | 288 | 256 | 22790
R3 | 1280 | 576 | 512 | 12330
R4 | 2240 | 1056 | 1024 | 10780
R5 | 5376 | 2304 | 2304 | 10000

##  3 Results

There are various definitions of jets and unbound outflows in the literature
(see e.g. Narayan et al., 2012; Yuan et al., 2015). One option to define the
jet region is by using $\left(-{T^{r}}_{t}/(\rho u^{r})\right)^{2}-1>1$, where
$-{T^{r}}_{t}$ is the radial ($\mu=1$) - temporal ($\nu=0$) component of the
stress-energy tensor (Equation 4), representing the outward radial energy
flux. This is equivalent to the region where the total energy per unit rest-
mass is greater than $\sqrt{5}c^{2}$ (; ). Alternatively the jet region can
also be defined as in Mościbrodzka & Falcke (2013), where the jet is taken to
be the unbound gas outflowing with a minimum bulk velocity $v_{min}/c=0.2$ as
measured in a normal observer frame, and the jet spine is taken to be the
region at which $\sigma>0.1$. In this paper, we follow the definitions in
Davelaar et al. (2018) in Kerr-Schild coordinates, where the jet is
characterised by two components: (1) The jet spine– a relativistic and
strongly magnetized outflow ($\sigma>1$). (2) The jet sheath– a mildly
relativistic outflow defined using the Bernoulli parameter
$Be=-\bar{h}u_{t}>1.02$ where $\bar{h}=(\rho+u_{g}+p_{g})/\rho$ is the
specific gas enthalpy and $u_{t}$ is the time component of the inverse four-
velocity.

For a qualitative perspective, Fig. 1 shows 2D slices through azimuthal angle
$\phi=0$ for density $\rho$, proton temperature $T_{p}$, and magnetisation
$\sigma=b^{2}/\rho$ for the simulations conducted at different resolutions. As
the simulation resolution increases, wave-like features in the jet sheath
become more prominent and plasmoids in current sheets close to the event
horizon are captured for resolution R5. Additionally, the jet sheath gets
hotter and less magnetised.

Figure 1: Meridional 2D slices through $\phi=0$ for (log) density $\rho$ (left
column), proton temperature $T_{p}$ in Kelvin (middle column), and
magnetisation $\sigma=b^{2}/\rho$ (right column) during a flux eruption event,
when the normalised magnetic flux on the event horizon $\varphi_{BH}$
approaches a local minimum (Fig. 3). From top to bottom: resolution R1, R3 and
R5. The grey dashed and black lines correspond to surfaces where $\sigma=1$
and $Be=1.02$, respectively. The jet has two components: the spine (defined as
the region in which $\sigma>1$) and the sheath (defined as the region in which
$\sigma<1$ and $Be>1.02$). The dominance of numerical floors undermines the
reliability of the jet spine. We partially cover the spine with a zero screen
for the plot of $\rho$ and time-average screens for $T_{p}$ and $\sigma$.
These screens are activated for $\sigma>3$, allowing visibility of plasma
within the range of $1<\sigma<3$.

###  3.1 Disk properties

In order to determine quantitative differences between the simulations we
compare various fluid parameters, as defined in the following subsections. The
averaged profile of a variable $X$ is calculated by integrating over $\theta$
and $\phi$ with a combination of three conditions to select the disk; the disk
is taken to be the region that satisfies density
$\rho/\rho_{\textup{max}}=\rho>10^{-4}$, magnetisation $\sigma<1$ and
Bernoulli parameter $Be<1.02$. Additionally, we include a density weight to
give more relevance to regions of the disk with higher density. 333Note the
difference in notation, $\left\langle\right\rangle$ without any sub/super
script is time average.

$\left\langle X\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}=\frac{\iint
X\rho\left(\rho>10^{-4}\right)\left(\sigma<1\right)\left(Be<1.02\right)\sqrt{-g}d\theta
d\phi}{\iint\rho\left(\rho>10^{-4}\right)\left(\sigma<1\right)\left(Be<1.02\right)\sqrt{-g}d\theta
d\phi}\text{ .}$ (6)

####  3.1.1 Fluxes

In this section we determine the impact of resolution on the temporal
evolution of the mass, magnetic and energy fluxes in the simulations. The
fluxes are defined as follows444Subscripts indicate where in the radial domain
the fluxes are analysed, e.g. $\dot{M}_{BH}$ at the event horizon and
$\dot{M}_{5r_{g}}$ at $5r_{g}$.. The mass accretion rate $\dot{M}$ is given
by:

$\dot{M}\equiv-\iint\rho u^{r}\sqrt{-g}d\theta d\phi\text{ .}$ (7)

$\dot{M}_{BH}$ is contaminated by artificial matter injection due to density
floors above and below the black hole poles up to a radius of $\approx
5r_{g}$. Therefore, $\dot{M}_{5r_{g}}$ is a better representation of the real
accretion process. The energy flux $\dot{E}$ is given by:

$\dot{E}\equiv-\iint{T^{r}}_{t}\sqrt{-g}d\theta d\phi\text{ .}$ (8)

The energy outflow efficiency of the jet spine is the energy return rate
divided by the time averaged mass accretion rate:

$\eta_{j}\equiv\frac{\dot{M}-\dot{E}_{\mathrm{spine}}}{\left\langle\dot{M}\right\rangle}\text{
.}$ (9)

To calculate the energy flux in the jet spine only $\dot{E}_{\mathrm{spine}}$,
we add an extra condition before doing the integral, which is the
magnetisation greater than 1 ($\sigma>1$). The magnetic flux is defined as:

$\Phi\equiv\frac{1}{2}\iint\left|B^{r}\right|\sqrt{-g}d\theta d\phi\text{ .}$
(10)

To study MAD simulations is convenient to use the normalised magnetic flux:

$\varphi\equiv\frac{\Phi}{\sqrt{\left\langle\dot{M}\right\rangle}}\text{ ,}$
(11)

known as the “MAD parameter” which, for spin $a=0.9375$ and torus scale height
$H/R\approx 0.3$, has the critical value $\varphi_{max}\approx 15$ (within the
units adopted here; Tchekhovskoy et al. (2012)). $H$ and $R$ are the full
height and cylindrical radius of the disk, respectively.

Space-time diagrams of the mass accretion rate (Equation 7) and normalised
magnetic flux (Equation 11) are presented in Fig. 2 for the five resolutions
of Table 1. A quasi-stationary state of the mass flux is obtained after
$t=5\times 10^{3}r_{g}/c$ until a radius of around $50r_{g}$. For $r\approx
50r_{g}$, there is a balance between inflow and outflow so
$\left\langle\dot{M}_{50r_{g}}\right\rangle_{[5-10]\times
10^{3}r_{g}/c}\approx 0$, that marks the point where the torus begins to
spread outwards due to viscosity. The relativistic jet and outflows launched
from the inner disk have short dynamical timescales. At $\approx 80r_{g}$, the
negative $\dot{M}$ is predominantly due to unbound outflows (spine and
sheath). At $\approx 100$ and $200r_{g}$, unbound outflows contribute $50\%$
and $10\%$ of the total $\dot{M}$, respectively. For example,
$\left\langle\dot{M}_{100r_{g}}\right\rangle_{[5-10]\times
10^{3}r_{g}/c}\approx-40$ in code units, with equal contribution from unbound
outflows and viscous spreading of the disk. In the analysis of averaged
profiles, the properties of the simulations are calculated as time averaged
between $t=\left[8-10\right]\times 10^{3}r_{g}/c$ to avoid the dependency on
the initial conditions for $r>100r_{g}$ and $t<8\times 10^{3}r_{g}/c$. The
disk properties still depend on the initial conditions of the Fishbone-
Moncrief torus for radii larger than $40r_{g}$.

Figure 2: Space-time diagrams of the mass accretion rate (top panel) and
normalised magnetic flux (bottom panel). From left to right: resolution R1 to
R5. In MAD states, the magnetic flux on the event horizon can build up until
it is saturated and gets ejected. In these episodic flux eruption events, the
magnetic pressure interrupts the accretion flow, causing flux tubes of low-
density material with vertical field to be pushed into the disk until a radius
in between $15$ and $40r_{g}$ (e.g. see the white boxes on the bottom panels).
Two major eruption events occur for R1, several small eruptions for R2 and R4
and three eruptions for R3 and R5. This behaviour is seen as well in the plot
of normalised magnetic flux versus time in the middle panels of Fig. 3.

The time evolution of the mass accretion rate $\dot{M}$ at $5r_{g}$, the
normalised magnetic flux $\varphi$ and the jet outflow efficiency $\eta_{j}$
at the event horizon are presented in Fig. 3. The magnetic flux grows until
saturation level at around $t=5\times 10^{3}r_{g}/c$. The fluctuations are
caused by accumulation and escape of field line bundles in the vicinity of the
BH. The magnetic pressure of those escaping field lines stops the accretion
process for a short time, therefore, a drop in the magnetic flux is correlated
to a sudden decrease of the accretion rate. We find that the mass accretion
rate decreases with increasing resolution when $\dot{M}$ is computed at a
radius of $5r_{g}$. However, we note that there is no clear tendency for
$\dot{M}$ when computed at the horizon– this is due to the artificial matter
injection that occurs close to the poles of the BH due to flooring (see Fig. 4
and Table. 2).

The peak, average value and standard deviation of $\dot{M}$, $\varphi$ and
$\eta_{j}$ are presented in Table 2. We use the modulation index MI as a
standardised measure of time variability of the accretion rates. It is defined
as the ratio of the standard deviation to the mean. In Fig. 4, the MI of
$\dot{M}_{5r_{g}}$ initially decreases from resolution R1 to R3, and from R3
it increases. Similarly, MI of $\varphi_{BH}$ increases from R2. However, the
differences in MI are not statistically significant. The accuracy of the MI
calculation is dependent on the length of the simulation time window of
analysis, $t=[5-10]\times 10^{3}r_{g}/c$. Longer simulations can provide a
more precise measure of the variability. Additionally at higher resolutions,
we find greater variability of proton temperature, magnetisation and opening
angle of the jet sheath close to the black hole at $\approx 40r_{g}$ (see
Section 3.2). All these differences could affect the predicted multiwavelength
spectra, relevant to the EHTC analysis of Sgr A∗where the two most promising
MAD models show higher light-curve variability that historical observations at
230 GHz (EHTC Sgr A∗2017 V).

Figure 3: Mass accretion rate, normalised magnetic flux and jet outflow efficiency as a function of time in solid lines, together with the time averaged values in horizontal dashed lines. The magnetic flux saturates after $t=5\times 10^{3}r_{g}/c$. The vertical dashed lines indicate the points where $\varphi_{BH}$ approaches a local minimum (Fig. 1). The colour convention in all figures is resolution R1: green, R2: blue, R3: red, R4: black, and R5: purple. Table 2: Peak, average value and standard deviation of the mass accretion rate at the event horizon and at $5r_{g}$, the normalised magnetic flux, and the energy outflow efficiency of the jet spine. Calculation interval $t=[5-10]\times 10^{3}r_{g}/c$. Res | $\dot{M}_{BH,max}$ | $\left\langle\dot{M}_{BH}\right\rangle$ | $\dot{M}_{5r_{g},max}$ | $\left\langle\dot{M}_{5rg}\right\rangle$ | $\varphi_{BH,max}$ | $\left\langle\varphi_{BH}\right\rangle$ | $\eta_{j,max}$ | $\left\langle\eta_{j}\right\rangle$
---|---|---|---|---|---|---|---|---
R1 | 75.9 | 41.6 $\pm$ 10.5 | 47.1 | 20.1 $\pm$ 7.3 | 13.6 | 10.6 $\pm$ 1.5 | 2.30 | 1.12 $\pm$ 0.33
R2 | 59.2 | 36.6 $\pm$ 7.7 | 35.2 | 17.2 $\pm$ 5.8 | 13.0 | 10.9 $\pm$ 1.2 | 2.01 | 1.09 $\pm$ 0.29
R3 | 70.0 | 35.0 $\pm$ 7.8 | 46.2 | 17.0 $\pm$ 5.6 | 13.3 | 10.5 $\pm$ 1.3 | 2.87 | 1.15 $\pm$ 0.33
R4 | 65.0 | 39.8 $\pm$ 9.4 | 38.2 | 16.2 $\pm$ 5.9 | 14.3 | 11.3 $\pm$ 1.4 | 2.01 | 1.02 $\pm$ 0.28
R5 | 72.7 | 40.1 $\pm$ 9.5 | 38.9 | 15.7 $\pm$ 6.7 | 14.2 | 11.4 $\pm$ 1.8 | 2.64 | 0.99 $\pm$ 0.29

Tchekhovskoy et al. (2011) ran simulations with resolution $288\times
128\times\left\\{32,64,128\right\\}$, seeded with large poloidal magnetic flux
to achieve the MAD state, enabling the production of efficient outflows. The
authors demonstrated the extraction of net energy from a spinning black hole’s
ergosphere, via the Penrose–Blandford–Znajek mechanism (Blandford & Znajek,
1977). The authors found an increase of the average energy outflow efficiency
from $\left\langle\eta_{j}\right\rangle\approx 0.30$ to $\approx 1.40\pm 0.15$
when the BH spin is increased from $a=0.5$ to $0.99$, respectively. McKinney
et al. (2012) extended those results by conducting a parametric study of MAD
flows with a similar resolution of $272\times 128\times 256$ around a rotating
black hole with different values of spin, poloidal and toroidal magnetic field
geometries. McKinney et al. (2012) confirmed that the Poynting flux outflow is
powered by the Blandford-Znajek mechanism when the poloidal magnetic flux
saturates near the BH, independently of the initial poloidal seed.
Additionally the authors found outflow efficiencies
$\left\langle\eta_{j}\right\rangle\gtrsim 1$ if $\left|a\right|\gtrsim 0.9$.
The simulations we present in this work are conducted with $a=0.9375$, and we
find outflow efficiencies
$\left\langle\eta_{j}\right\rangle=\left[0.99-1.15\right]\pm 0.30$. We find
that the outflow efficiency is converged within statistical uncertainty for
the five resolutions considered (Table 2), where resolution R1 is similar to
those used in Tchekhovskoy et al. (2011); McKinney et al. (2012).
Additionally, we find that $\eta_{j}$ fluctuates with time (Fig. 2) and has
maximum values up to $2.87$ (Table 2).

Figure 4: Resolution dependence of average values and standard deviations of
mass accretion rate at the event horizon (top panel) and at $r=5r_{g}$ (middle
panel) and normalised magnetic flux (bottom panel). The horizontal grey line
in each panel is the average value of all five resolutions, respectively. The
modulation index is shown on the right. The error bars represent the standard
deviations. Calculation interval $t=[5-10]\times 10^{3}r_{g}/c$.

####  3.1.2 Flux eruptions

The variation in the amount of major flux eruptions observed across different
resolutions, as depicted in Fig. 2 and 3, raises two key points: 1)
Periodicity may not be directly linked to resolution, and 2) variations in
periodicity suggest that it might not be inherently indicative of a periodic
physical phenomenon. However, the statistics of these findings are limited by
the relatively short duration of our high resolution simulations (see e.g.
Narayan et al., 2012).

Fig. 5 displays the mid-plane density structures of the disks for all
resolutions at different times. Low-density bubbles supported by magnetic
pressure initially rise buoyantly away from the black hole and tend to orbit
at a radius of $\approx 40r_{g}$. The sizes of these bubbles are relatively
constant regardless of resolution. Smaller flux eruptions, such as shorter
ones, tend to yield smaller bubbles, while larger flux eruptions typically
lead to the formation of larger bubbles. Gravity naturally pulls inward,
directing denser material towards the center. In this case, the lower density
bubbles are prone to gravitational mixing instabilities such as RTI
(Tchekhovskoy et al., 2011; Porth et al., 2021; Ripperda et al., 2022;
Zhdankin et al., 2023). At the lowest resolution R1, numerical diffusion
causes the bubbles to mix. At higher resolution, finer structures are resolved
because of mixing instabilities occurring at the surfaces of the bubbles.

There is a mixing of the low density bubbles with the higher density plasma in
the accretion disk. This mixing causes mass to be advected inside the bubbles
leading to a subsequent increase in the bubble density. At resolution R2,
single plumes of high density plasma develop, and these plumes often split a
bubble into two separate bubbles. From resolution R3 and higher, increasingly
more plumes develop simultaneously in a single bubble.

Additional work is needed to quantitatively determine the resolution
dependence on the time evolution of the volume, magnetisation, temperature and
energy contained within the bubbles (see e.g. Porth et al., 2021), the exact
role of the mixing instabilities, how mixing impacts the polarisation
signatures of magnetic flux eruptions in the radio band (e.g. Jia et al.,
2023; Davelaar et al., 2023; Najafi-Ziyazi et al., 2024), or non-thermal
particle acceleration (e.g. Zhdankin et al., 2023)

Figure 5: Equatorial 2D slices through $\theta=\pi/2$ of density $\rho$
showing low density and high magnetisation flux tubes penetrating the
accretion disk up to a radius of $\approx 40r_{g}$. From left to right:
resolution R1 to R2 (top) and R3 to R5 (bottom).

####  3.1.3 Viscosity Parameters

In this section, we study turbulence in the magnetised disk, which serves as a
means to transfer angular momentum and energy. By making a decomposition of
the velocity
$\left\langle\vec{u}\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}+\delta\vec{u}$
(disk average plus fluctuation) in the momentum conservation equation, a
turbulent stress tensor is obtained. This tensor has two components: Reynolds
and Maxwell stress tensors, which represent the statistical correlations of
velocity and magnetic field components respectively. Angular momentum
transport and extraction of available energy from the flow to fuel the
turbulent fluctuations is provided by the turbulent stress tensor in the
orbiting plasma (Balbus & Hawley, 1998).

A viscosity parameter is often related to the strength of the radial and
azimuthal turbulent speed correlations, $\delta{u^{r}}=u^{r}-\left\langle
u^{r}\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}$ and
$\delta{u^{\phi}}=u^{\phi}-\left\langle
u^{\phi}\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}$, respectively. We
recast the turbulent stresses in the classical alpha-models as follows (e.g.
Shakura & Sunyaev, 1973; Balbus & Hawley, 1998; Liska et al., 2019; Chatterjee
& Narayan, 2022). The Reynolds viscosity is defined as:

$\alpha_{r}\equiv\frac{\left\langle\left(\rho+u_{g}+p_{g}+b^{2}\right)\delta{u^{r}}\sqrt{g_{rr}}\delta{u^{\phi}\sqrt{g_{\phi\phi}}}\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}}{\left\langle
p\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}}\text{ ,}$ (12)

and the Maxwell viscosity as:

$\alpha_{M}\equiv-\frac{\left\langle
b^{r}\sqrt{g_{rr}}b^{\phi}\sqrt{g_{\phi\phi}}\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}}{\left\langle
p\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}}\text{ .}$ (13)

where, $p=b^{2}/2+p_{g}$ is the total pressure. The time averaged viscosities
in the disk are presented in bottom left panel of Fig. 6. We find that the
Maxwell component dominates over the Reynolds component. The negative values
of $\alpha_{r}$ indicate inward angular momentum transport and suggest
convection-like behaviour (e.g. Begelman et al., 2022). The relatively low
values of $\alpha_{r}$ indicate that the angular momentum transport due to
turbulent convection is subdominant (e.g. Chatterjee & Narayan, 2022). As
demonstrated in previous works (e.g. Liska et al., 2019), the combination of
both Reynolds and Maxwell viscosities is not sufficient to explain the total
angular momentum transport, provided by the effective viscosity:

$\alpha_{\textup{eff}}\equiv-\frac{\left\langle
v^{r}\sqrt{g_{rr}}v^{\phi}\sqrt{g_{\phi\phi}}\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}}{\left\langle
c_{s}^{2}\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}}\text{ ,}$ (14)

where $v^{r}=u^{r}/u^{t}$ and $v^{\phi}=u^{\phi}/u^{t}$ are the physical
radial and azimuthal velocity components, and the sound speed $c_{s}$ is:

$c_{s}\equiv\left|\frac{\gamma_{\mathrm{ad}}(\gamma_{\mathrm{ad}}-1)u_{g}}{\rho+u_{g}+p_{g}}\right|^{1/2}\text{
.}$ (15)

The magnetorotational instability (MRI) is responsible for the generation of
MHD turbulence that leads to enhanced outward angular momentum transport in
accretion discs (Shakura & Sunyaev, 1973; Balbus & Hawley, 1998; Pessah et
al., 2006). In this context, the Maxwell stress is associated to MRI. As seen
in Fig. 6, for $r=\left[3-50\right]r_{g}$, there is a significant contribution
of MRI to transport angular momentum, since the Maxwell viscosity is only a
factor of $\approx 3$ lower than the effective viscosity. For the very inner
part of the accretion disk, $r<3r_{g}$, $\alpha_{\textup{eff}}$ keeps growing
while $\alpha_{M}$ drops significantly. For MAD models, magnetic flux
eruption–driven disk winds cause a strong vertical flow of angular momentum
(Chatterjee & Narayan, 2022) and jets can remove angular momentum close to the
black hole (e.g. Tchekhovskoy et al., 2012; Narayan et al., 2022).

Figure 6: Time averaged, disk averaged profiles of variables as a function of
radius. Fist panel: plasma beta $\beta=p_{g}/p_{b}$, gas-to-magnetic pressure
ratio. Second panel: density. Third panel: viscosity parameters, Maxwell,
Reynolds and Effective viscosity. Fourth panel: full height and critical
wavelength. Fifth panel: quality factors. All the variables are averaged in
time $t=[8-10]\times 10^{3}r_{g}/c$. For plasma beta, shaded region depicts
the range of variation within one standard deviation for R1 and R5. In flux
eruption events, highly magnetised plasma from the jet spine penetrates the
disk mid-plane ($r\lesssim 40r_{g}$) and propagates above and below the mid-
plane ($40r_{g}\lesssim r\lesssim 100r_{g}$). As a consequence, the time
averaged magnetic pressure (plasma beta) in the disk increases (decreases)
with resolution.

####  3.1.4 Magnetorotational Instability

The MRI quality factor components $Q^{i}$ in the $i=[r,\theta,\phi]$-direction
are defined as the number of cells available for resolving the fastest-growing
MRI mode wavelength:

$Q^{i}\equiv\frac{2\pi rv_{A}^{i}}{\Delta x^{i}v}\text{ ,}$ (16)

where the size of the cell is $\Delta x^{i}\equiv dx^{i}\sqrt{g_{ii}}$. The
Alfvén speed, which is the characteristic speed of transverse shear waves in a
magnetized fluid (Balbus & Hawley, 1998), is defined as:

$v_{A}^{i}\equiv\left(\frac{\left|b^{i}b^{i}g_{ii}\right|}{\rho+u_{g}+p_{g}+b^{2}}\right)^{1/2}\text{
,}$ (17)

and the speed of the fluid is:

$v\equiv\frac{\left(u^{r}u^{r}g_{rr}+u^{\theta}u^{\theta}g_{\theta\theta}+u^{\phi}u^{\phi}g_{\phi\phi}\right)^{1/2}}{u^{t}}\text{
.}$ (18)

The quality factors in all three directions are presented in the bottom panel
of Fig. 6. For all simulations, the quality factors are large enough to
sustain the MRI-driven turbulence. The only exception is the polar component
$Q^{\theta}$ both near the event horizon ($r<3r_{g}$) and further away in the
disk ($40r_{g}\leq r\leq 100r_{g}$) for the lowest resolution R1.
$Q^{\theta}<10$ due to the coarse size of the cell and to the relatively lower
Alfvén speed in the polar direction. However, this does not significantly
affect the time-averaged accretion rate or normalised magnetic flux, as they
consistently converge across all five resolutions.

Apart from the saturation of magnetic flux, another characteristic of the MAD
state is the suppression of the MRI due to the strong magnetic fields close to
the the event horizon. Igumenshchev et al. (2003); Igumenshchev (2008) in
pseudo-Newtonian MHD, and Tchekhovskoy et al. (2011); McKinney et al. (2012)
via GRMHD simulations, demonstrated the suppression of MRI when poloidal
magnetic flux saturates near the BH, forming a MAD. White et al. (2019) showed
that the suppression of MRI is robust with resolution. From the general-
relativistic dispersion relation connecting the Alfvén frequency (as seen by
an observer at infinity) to oscillation frequency, Gammie (2004) proposed an
instability criterion for the critical wavelength:

$\lambda_{crit}\equiv 2\pi
v_{A}^{z}\frac{r^{3/2}+a}{3^{1/2}\left(1-2r^{-1}+a^{2}r^{-2}\right)^{1/2}}\text{
.}$ (19)

The MRI is suppressed if the critical wavelength is bigger than the full
height of the disk,
$\left\langle\lambda_{crit}\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}>H$
(White et al., 2019). The full height is defined geometrically as $H\equiv
r\left\langle\left|\theta-\pi/2\right|\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}$.

As a consequence of the disk magnetic field strength increasing with
resolution, the MRI is more suppressed according to the criterion
$\left\langle\lambda_{crit}\right\rangle_{\rho[\theta,\phi]}^{\textrm{disk}}>H$.
The degree of the suppression of MRI in GRMHD simulation depends on whether
the simulation is performed in 2D or 3D, alongside the numerical resolution
used (see e.g. White et al., 2019). In axisymmetric 2D simulations, the MRI
turbulence decays after time intervals of $\sim 1000r_{g}/c$ (e.g. Hide &
Palmer, 1982; Balbus & Hawley, 1991; Guan & Gammie, 2008) and the most
unstable axisymmetric linear mode of MRI is stabilised (e.g. McKinney et al.,
2012; Begelman et al., 2022). Begelman et al. (2022) found that MRI is not
completely suppressed in non-axisymmetric MAD flows, where the saturation of
magnetic flux is driven by radial convective/interchange instabilities
triggered by a dominant toroidal field. The MRI is a MHD instability
originally introduced for axisymmetric, magnetized accretion disks (Balbus &
Hawley, 1998). Goedbloed & Keppens (2022) argue that axisymmetric MRIs
represent only a finite amount of unstable modes and that turbulence,
accretion and dynamo activity in 3D disks is governed by the excitation of
nonaxisymmetric super-Alfvénic rotational instabilities (SARIs).

####  3.1.5 Plasmoids

During an eruption in MAD states, the expulsion of magnetic flux goes through
the magnetic reconnection of field lines in a current sheet in the equatorial
plane. If this current sheet is sufficiently thin (see criteria for current
sheet thickness $\delta$ below) such that plasmoid–mediated reconnection is
resolved, then the current sheet can show signs of the tearing instability
(Ripperda et al., 2020; Ripperda et al., 2022). Evidence of the plasmoid
instability acting in the current sheet signifies that the rate at which
magnetic energy is converted to heat, the reconnection rate, has converged to
a value of $v_{\textup{rec}}=0.01v_{A}$ (Bhattacharjee et al., 2009; Uzdensky
et al., 2010), where $v_{A}=\sqrt{\sigma/(\sigma+1)}c$ is the Alfvén speed.
The plasmoid instability is typically triggered when the ratio between the
alfvénic timescales and diffusion timescales, the Lundquist number $S\equiv
v_{A}\omega/\eta\gtrsim 10^{4}$, where $\omega$ is the length of the current
sheet (Bhattacharjee et al., 2009; Uzdensky et al., 2010). In this picture the
thickness of the reconnection layer, $\delta$, or rather its aspect ratio
$\delta/\omega$, must match the asymptotic reconnection rate to minimally
resolve the tearing instability,
$v_{\textup{rec}}=v_{in}/v_{A}=0.01=\delta/\omega$, where $v_{in}$ is the
inflow velocity into the layer, equal to the
$\mathbf{E}\times\mathbf{B}/B^{2}$-velocity in the direction perpendicular to
the layer.

In ideal (GR)MHD, a current sheet, i.e., a magnetic null in the case of zero
guide field as applicable to the equatorial current sheet in a MAD state
(Ripperda et al., 2022), is only captured by a single cell (since it is
governed by a numerical resistivity, see below). Therefore, for a current
sheet that tears at $r=2r_{g}$ (approximately at the ergosphere), one requires
1 cell per $\delta=0.01r_{g}$ for a current sheet of length $\omega=1r_{g}$.
Comparing to the domain size (poloidal angle $\theta\in[0,\pi]$), this results
in $2\pi r_{g}/0.01r_{g}\gtrsim 600$ cells at $2r_{g}$, implying that only
simulations R4 and R5 have enough resolution in $\theta$ to minimally resolve
the tearing instability near the horizon for current sheets of length
$\omega\geq 1r_{g}$, as indicated by our results below.

In ideal GRMHD simulations, such as those conducted here, reconnection is
mediated via numerical resistivity $\eta_{\textup{num}}$ instead of an
explicit resistivity $\eta$ (Ripperda et al., 2019). We can determine whether
our simulations are computed with sufficient numerical accuracy in order to
resolve plasmoid-mediated reconnection, and thus form plasmoids by calculating
the numerical Lundquist number $S\equiv v_{A}\omega/\eta_{\textup{num}}$,
where $\eta_{\textup{num}}=\hat{C}\Delta x$ and the coefficient $\hat{C}$
depends on the accuracy and order of the reconstruction scheme (Zhang et al.,
2003). In the plasmoid instability regime, the rate of magnetic flux decay on
the horizon converges to $\varphi_{BH}\propto e^{-t/500}$ (Bransgrove et al.,
2021; Ripperda et al., 2022). In simulations R4 and R5 we observe plasmoids in
current sheets of length $\omega\sim 10r_{g}$ close to the event horizon,
indicating that resolutions of a minimum of $2240\times 1056\times 1024$ are
required to resolve the plasmoid instability, see Fig. 4 and Appendix C of
Ripperda et al. (2022). Based on these results, the constant $\hat{C}\approx
0.03$ is obtained by setting $10^{4}=10/\eta_{\textup{num,R4}}$, where
$\eta_{\textup{num,R4}}$ is the numerical resistivity at $10r_{g}$ for R4. The
Alfvén speed $v_{A}=\sqrt{25/26}c\approx c=1$ for the reconnecting plasma
because the reconnection is fed by the plasma in the jet at
$\sigma_{\textup{max}}=25$. Due to the use of a spherical grid, it is harder
to resolve plasmoid-mediated reconnection at larger radii and in smaller
current sheets, due to the reduction in resolution at larger radii.

####  3.1.6 Disk profiles

The time averaged profiles of the disk plasma beta (ratio of gas-to-magnetic
pressure $\beta=p_{g}/p_{b}$), density, viscosity parameters (Maxwell,
Reynolds and effective viscosity), full height, critical wavelength and
quality factors are presented in Fig. 6. All variables presented in Fig. 6 are
averaged between $8$ to $10\times 10^{3}r_{g}/c$.

2D slices through azimuthal angle $\phi=0$ of plasma beta $\beta$ are shown in
Fig. 7. We find that the jet spine gets more magnetised with increasing
resolution as demonstrated in Section 3.2 (Fig. 8). During flux eruption
events, highly magnetized flux tubes originating from the jet spine penetrate
the disk and propagate until $r\lesssim 40r_{g}$. The pressure maximum of the
initialised Fishbone-Moncrief torus is located at $r=41r_{g}$– at around this
radius, the flux tubes split into two and are ejected as outflows above and
below the disk mid-plane. These ejections increase the time-averaged disk
magnetic pressure, and consequently reduce the plasma beta in the disk, as
observed in Fig. 6. Plasma beta in the disk for $r\gtrsim 3r_{g}$ is
convergent for resolutions R3 and higher. For resolution R5, plasma beta for
$r\lesssim 3r_{g}$ is slightly lower than the rest of simulations. Very close
to the BH $r\lesssim 3r_{g}$, the Maxwell viscosity drops significantly for
all resolutions. In this inner part of the accretion disk, the vertical
magnetic pressure is responsible for launching winds that transport angular
momentum outwards (e.g. Chatterjee & Narayan, 2022).

Figure 7: Meridian 2D slices through $\phi=0$ for plasma beta
$\beta=p_{g}/p_{b}$. From left to right: resolution R1, R3 and R5. The grey
dashed and black lines correspond to surfaces where $\sigma=1$ and $Be=1.02$,
respectively. We partially cover the spine with a zero screen activated for
$\sigma>3$, allowing visibility of plasma within the range of $1<\sigma<3$. As
resolution increases, the polar axis features a smaller volume and therefore
there is less dissipation of magnetic energy along the pole. In flux eruption
events, highly magnetised plasma from the jet spine penetrates the disk,
causing the time averaged plasma beta in the disk to decrease with resolution.

###  3.2 Jet Properties

The averaged profile of a jet variable does not involve any density or
magnetisation weight. For the jet spine:

$\left\langle X\right\rangle_{[\theta,\phi]}^{\textrm{spine}}=\frac{\iint
X\left(\sigma>1\right)\sqrt{-g}d\theta
d\phi}{\iint\left(\sigma>1\right)\sqrt{-g}d\theta d\phi}\text{ ,}$ (20)

and for the jet sheath:

$\left\langle
X\right\rangle_{[\theta,\phi]}^{\textrm{sheath}}=\frac{\iint_{\left|\pi-\theta\right|\leqslant\pi/6}X\left(\sigma<1\right)\left(Be>1.02\right)\sqrt{-g}d\theta
d\phi}{\iint_{\left|\pi-\theta\right|\leqslant\pi/6}\left(\sigma<1\right)\left(Be>1.02\right)\sqrt{-g}d\theta
d\phi}\text{ .}$ (21)

We follow the methodology used in Chatterjee et al. (2019) to investigate the
instabilities occurring in the jet sheath-spine interface. Decomposing the
stress-energy tensor
${T^{r}}_{t}\equiv{\left(T_{H}\right)^{r}}_{t}+{\left(T_{M}\right)^{r}}_{t}$
(Equation 4), with the hydro component:

${\left(T_{H}\right)^{r}}_{t}\equiv\left(\rho+u_{g}+p_{g}\right)u^{r}u_{t}\text{
,}$ (22)

and magnetic component:

${\left(T_{M}\right)^{r}}_{t}\equiv b^{2}u^{r}u_{t}-b^{r}b_{t}\text{ ,}$ (23)

We define the total specific energy, measuring the maximum Lorentz factor if
all forms of energy are converted into kinetic energy,

$\mu\equiv\frac{\left\langle{T^{r}}_{t}\right\rangle_{[\theta,\phi]}^{\textrm{spine}}}{\left\langle-\rho
u^{r}\right\rangle_{[\theta,\phi]}^{\textrm{spine}}}\text{ ,}$ (24)

the magnetisation, as a measure of conversion efficiency of kinetic to
magnetic energy,

$\hat{\sigma}\equiv\frac{\left\langle{\left(T_{M}\right)^{r}}_{t}\right\rangle_{[\theta,\phi]}^{\textrm{spine}}}{\left\langle{\left(T_{H}\right)^{r}}_{t}\right\rangle_{[\theta,\phi]}^{\textrm{spine}}}\text{
,}$ (25)

the specific enthalpy,

$h\equiv\frac{\left\langle{T^{r}}_{t}\left(u_{g}+p_{g}\right)/\rho\right\rangle_{[\theta,\phi]}^{\textrm{spine}}}{\left\langle{T^{r}}_{t}\right\rangle_{[\theta,\phi]}^{\textrm{spine}}}\text{
,}$ (26)

and the Lorentz factor,

$\gamma\equiv\frac{\left\langle{T^{r}}_{t}u^{t}\sqrt{-1/g^{tt}}\right\rangle_{[\theta,\phi]}^{\textrm{spine}}}{\left\langle{T^{r}}_{t}\right\rangle_{[\theta,\phi]}^{\textrm{spine}}}\text{
.}$ (27)

The energy equation can be obtained by combining the four previously defined
variables:

$\mu=\gamma\left(\hat{\sigma}+1\right)\left(h+1\right)\text{ .}$ (28)

####  3.2.1 Jet spine

The top panel of Fig. 8 shows the total specific energy, magnetisation,
enthalpy, Lorentz factor and half opening angle of the upper jet spine
(Equation 20), each of which are averaged over the time interval
$t=[8-10]\times 10^{3}r_{g}/c$. We find that the Lorentz factor is nearly the
same for all resolutions presented, and increases from approximately $2$ to
$7$ with radius $20-10^{3}r_{g}$. The total specific energy is convergent and
nearly constant with radius. However, the magnetisation and enthalpy are
convergent only for resolutions R4 and R5 and are not physically trustworthy.

Since numerical dissipation is proportional to the cell size, as the
resolution increases the cells get smaller. Consequently, cells near the polar
axis have a smaller volume which results in less dissipation of magnetic
energy and heating along the pole. As a result the jet spine gets more
magnetised, $\hat{\sigma}$ drops notably slower with radius, and the internal
gas energy and pressure decrease with resolution as shown in Fig. 1.
Properties such as magnetisation and enthalpy in the jet become unreliable if
number of cells per half opening angle drops below 10 (see bottom panel of
Fig. 8). Additionally, the polar axis introduces considerable uncertainty into
thermodynamics– while the transmissive boundary conditions used at the poles
effectively capture the flow of energy, they artificially change the
conversion of one energy source into another.

Figure 8: Jet spine ($\sigma>1$) properties as a function of $r$, averaged in
time $t=[8-10]\times 10^{3}r_{g}/c$. Total specific energy, Lorentz factor,
magnetisation and enthalpy (top panel). Half opening angle (middle panel) and
cells per half opening angle (bottom panel). Shaded regions depict the range
of variation within one standard deviation for R1 and R5. $\hat{\sigma}$ and
$h$ are not physically trustworthy. Additionally, the jet spine is wider for
R4 and R5.

The half opening angle of the jet spine is slightly higher for R4 and R5,
especially close to the base of the jet. The standard deviation of the half
opening angle decreases with radius. Close to the black hole, $r<50r_{g}$, we
see high fluctuations in the width of the jet spine and sheath due to the
magnetic flux eruptions in the MAD states.

####  3.2.2 Jet sheath

The top panel of Fig. 9 shows the total specific energy, magnetisation,
enthalpy, Lorentz factor, and half opening angle of the jet sheath, all
averaged over the time interval $t=[8-10]\times 10^{3}r_{g}/c$. In general,
all these variables in the sheath are lower than the corresponding values in
the spine, except the enthalpy for $r\lesssim 60r_{g}$ where
$h_{\mathrm{sheath}}>h_{\mathrm{spine}}$ for R4 and R5. The enthalpy for R1
has an inflection point around $100r_{g}$, where $h$ starts to increase due to
enhanced numerical dissipation caused by a lack of resolution in the polar
direction.

When a flux eruption occurs, $\hat{\sigma}$ drops locally at the base of the
jet generating a wave that propagates outwards. This wave is accompanied by a
simultaneous increase of enthalpy, creating a thermal pressure gradient that
slows down the jet outflow (see Fig. 10). However, the Lorentz factor of the
jet $\gamma$ only decreases a few percent for all resolutions. In contrast,
Chatterjee et al. (2019) found a significant reduction of $\gamma$ at
$r\approx 10^{3}r_{g}$ when pinch instabilities cause magnetic dissipation,
that raise the jet specific enthalpy to order unity. The fact that we do not
see such a significant reduction of $\gamma$ in our simulations is a
consequence of the much smaller disc that collimates the jet out to
approximately $500r_{g}$. For $r>500r_{g}$, the diminishing confining pressure
exerted by the disk wind triggers a lateral expansion of the jet, effectively
suppressing pinch instabilities (see, e.g., Porth & Komissarov, 2015). This
expansion of the jet leads to acceleration and increase of $\gamma$.

Mixing instabilities due to shearing motions across the interface between the
jet and the accretion disk can potentially cause the dissipation of magnetic
energy into heat. Resolving these mixing instabilities is crucial for
comprehending the energetic dynamics of jets at large radii, as mixing will
lead to mass loading of the jet. Close to BH, the jet is primarily loaded by
material from density floors. We see wave-like features in the jet sheath,
which become more prominent (less diffusive) at higher resolutions. These
wave-like features may induce mixing between jet, sheath and disk (see Fig. 1
and 7). In addition to affecting the degree of mass loading in the jet, the
wave-like features in the jet sheath may also affect the magnetic field in the
proximity of the jet sheath: Davelaar et al. (2023) found that waves
propagating along the jet-wind shear layer alter the orientation of the
magnetic field lines.

For $r\lesssim 100r_{g}$, the half opening angle of the sheath features strong
fluctuations and on average the sheath is thinner at resolutions R3, R4 and R5
(see bottom panel of Fig. 9). Between $100$ and $\approx 500r_{g}$, the sheath
is thinner at resolution R5. This slimming of the jet sheath can result in an
increase of enthalpy and reduction of magnetisation in the sheath, while
keeping an approximately constant total specific energy for all resolutions:
These results have a direct impact on the temperature, as described in Section
3.2.3.

Figure 9: Jet sheath properties as a function of $r$, averaged in time
$t=[8-10]\times 10^{3}r_{g}/c$. Total specific energy, Lorentz factor,
magnetisation and enthalpy (top panel). Half opening angle of the jet sheath
measured from the jet spine (bottom panel). Shaded regions depict the range of
variation within one standard deviation for R1 and R5. Mixing instabilities
occurring in the jet sheath-spine interface cause a reduction of magnetisation
and increase of enthalpy with resolution. Figure 10: Space-time diagrams of
magnetisation (top panels) and enthalpy (bottom panels) of the jet sheath.
From left to right: resolution R1 to R5. A localised decrease in
$\hat{\sigma}$ is accompanied by a simultaneous rise in enthalpy, generating a
wave that propagates outwards.

####  3.2.3 Proton temperature

The proton temperature in the jet sheath and disk is presented in Fig. 11. The
protons in the disk are near virial temperature $\propto r^{-1.3}$, and
$T_{p}$ in the disk is converged for all resolutions we consider. On the other
hand, $T_{p}$ in the sheath is found to increase with increasing resolution.
Processes causing a mixing of the colder jet sheath with the hotter jet spine
can include: (1) mixing instabilities due to shearing motion across the
interface between the jet sheath and funnel, (2) waves generated at the jet
base in MAD states that propagate along the jets and, (3) wobbling of the jet
itself (see e.g. Wong et al., 2021; Davelaar et al., 2023). In addition, the
increase we see in the magnetisation of the jet spine as the resolution is
increased also impacts $T_{p}$ in the sheath. Highly magnetised plasma from
the jet spine supplies matter to the current sheet, with
$T_{p}\propto\sigma_{spine}$, and the reconnection exhaust deposits this hot
plasma in the jet sheath up to at least $20r_{g}$ for resolutions R4 and R5
(Ripperda et al., 2022). For $r<100r_{g}$, the sheath is thinner at higher
resolution (see Fig. 9), the amount of dissipation/heating gets distributed
across a smaller volume of plasma, leading to higher maximum temperature. The
jet sheath (proton) temperatures are non fully convergent with resolution.
Despite this, $T_{p}$ varies by less than a factor of $\sim 2$ between the
lowest and highest resolutions.

Figure 11: Proton temperature in the jet sheath (top panel) and in the disk
(bottom panel). All averaged in $t=[8-10]\times 10^{3}r_{g}/c$. Shaded regions
depict the range of variation within one standard deviation for R1 and R5.

The proton temperature in the jet sheath follows a shallower profile, $\propto
r^{-0.3}$, than in the disk $\propto r^{-1.3}$. The jet is well collimated by
the large disk when the numerical resolution is adequate for capturing the jet
at large radii away from the black hole. The proton temperature profile in the
jet sheath is different for the lowest resolution level R1. In this case,
$T_{p}$ remains nearly constant for radii between 100 and 500 $r_{g}$, and
then follows a power law of $\propto r^{-0.3}$ at larger radii. This
discrepancy arises from a lack of resolution in the polar direction at
$r>100r_{g}$ along the jet axis for R1 compared to the other resolutions we
present.

Our simulations show that the adiabatic decompression of the jet sheath begins
at around $500r_{g}$, where the temperature gradient is much steeper than
$r^{-0.3}$ (Fig. 11). It is worth mentioning that the change in steepness of
the temperature gradient depends on the size of the initial torus, and it is
unknown whether this adiabatic decompression happens in real systems. For
instance, observations utilising the Very Large Array (VLA) have revealed a
parabolic collimation profile in the M87∗ jet extending up to $10^{5}r_{g}$
(Asada & Nakamura, 2012; Hada et al., 2013; Mertens et al., 2016; Kim et al.,
2018; Nakamura et al., 2018).

##  4 Conclusions

GRMHD codes face significant numerical challenges when used to perform
simulations in the MAD regime. The higher magnetic flux characteristic of the
MAD regime leads to new dynamics, including interchange-type accretion modes,
suppression of MRI and flux eruptions governed by plasmoid dominated magnetic
reconnection. The occurrence and dynamics of these phenomena may be affected
by the numerical resolution of GRMHD MAD simulations. Recently, White et al.
(2019) conducted a resolution study of MAD accretion flows, with simulations
reaching a maximum effective resolution of $512\times 256\times 512$. They
showed that the general large-scale structure of the accretion flow is robust
with resolution, for the resolutions they considered. However, they found that
the spatial structure and Lorentz factor of the jet, small-scale features of
the turbulence and variability of modelled synchrotron emission were not fully
converged with resolution.

We have conducted a study of MAD states and jets across a wider range of
resolutions than has been previously explored, up to resolution of $5375\times
2304\times 2304$ in a logarithmic spherical-polar grid, using existing
simulations in Ripperda et al. (2022). Such extreme resolution is needed to
achieve convergence in the reconnection rate in the plasmoid-dominated regime
in GRMHD, which is important for the timescale of variability. We divided the
system in three components: the jet spine (relativistic and strongly
magnetised outflow), the jet sheath (mildly relativistic outflow) and the
disk. Below, we elucidate how our findings might influence ongoing studies
within the black hole accretion community.

We find that the time-averaged disk properties are consistent for all
resolutions with only small differences in flux variability, plasma beta and
inferred viscosity parameters close to the BH. For resolution $2240\times
1056\times 1024$ and higher, we find that the time-averaged Reynolds
viscosity, attributed to turbulent convection, indicates inward angular
momentum transport around $10r_{g}$. Higher resolutions resolve finer
structures when mixing instabilities occur at the surfaces of flux tubes
during flux eruptions (see Fig. 5). This mixing may function as mechanism for
energy dissipation, non-thermal particle acceleration and flare generation
across various wavelengths.

For the jet spine, the total specific energy and Lorentz factor is nearly the
same for all resolutions we consider. We find that the magnetisation and
enthalpy are convergent for resolution $2240\times 1056\times 1024$ and
higher. However, these four jet spine properties do depend on, and are
governed by, the chosen magnetisation ceiling in the code setup. For this
reason, we recommend caution in drawing conclusions based on jet dynamics in
general, as the jet spine is dominated by numerical floors and the
magnetisation and enthalpy are not physically reliable.

We find that the jet-disk interface, the "sheath" is the most sensitive region
to resolution. There are large fluctuations in the width of the jet sheath due
to magnetic flux eruptions for $r\lesssim 50r_{g}$, more prominent at higher
resolutions. These fluctuations directly influence the variability of
properties within the jet sheath, notably the temperature and magnetisation.
At higher resolutions, we see more resolved wave-like features in the jet
sheath, that may induce mixing between the sheath and plasma within
$1\lesssim\sigma\lesssim 3$ (see Fig. 1 and 7). Mixing processes may include
(1) shearing motion across the interface between the jet sheath and spine, (2)
propagation of waves generated at the jet base, and (3) wobbling of the jet
itself. Additionally for resolution $2240\times 1056\times 1024$ and
$5376\times 2304\times 2304$, highly magnetised plasma from the jet spine
supplies matter to the current sheets, and the reconnection exhaust deposits
this hot plasma into the jet sheath (Ripperda et al., 2022). At higher
resolutions, we find that the jet sheath gets thinner, resulting in increased
temperature, reduced magnetisation, and greater variability of $T_{p}$,
$\hat{\sigma}$ and opening angle, while keeping an approximately constant
total specific energy for all resolutions. These differences could affect the
predicted multiwavelength spectra coming from the jet sheath, relevant to e.g.
the EHT collaboration and VLBI imaging of large-scale jets. In the next paper
on jet images, we will analyse whether a hotter, thinner jet sheath could
result in increased limb-brightening in the observed images (e.g. Kim et al.,
2018; Janssen et al., 2021) or potentially lead to filamentary structures
observed in jets (Fuentes et al., 2023).

Capturing plasma mixing in flux tubes and along the jet-disk interface is
essential for improved comparisons between observations and simulations. The
presence of mixing instabilities could significantly influence the dissipation
of energy, that could result in distinct multiwavelength emission from high-
resolution simulations (e.g. Sironi et al., 2021; Zhdankin et al., 2023). In
forthcoming research we will investigate the influence of mixing instabilities
on radiative emission in both total intensity and polarisation. We will focus
on analysing non-thermal electron distribution functions on the dissipative
regions. We will aim to simulate three-dimensional systems capable of
accurately capturing jets extending up to $10^{5}r_{g}$ or more. This research
is relevant to studies of AGN, such as M87∗. By doing so, we aspire to deepen
our comprehension of the fundamental horizon-scale physical processes and
their implications for the dynamics and emission properties of larger-scale
jets.

## Acknowledgements

We gratefully recognise Sasha Tchekhovskoy and Doosoo Yoon for insightful
discussions. L.S and S.M. were supported by a Dutch Research Council (NWO)
VICI award, grant No. 639.043.513 and by a European Research Council (ERC)
Synergy Grant "BlackHolistic" grant No. 101071643. In addition, L.S was
supported by the Colfuturo Scholarship. G.M. was supported by a Canadian
Institute of Theoretical Astrophysics (CITA) postdoctoral fellowship and by a
Netherlands Research School for Astronomy (NOVA), Virtual Institute of
Accretion (VIA) postdoctoral fellowship. K.C. was supported in part by grants
from the Gordon and Betty Moore Foundation and the John Templeton Foundation
to the Black Hole Initiative at Harvard University, and by NSF award
OISE-1743747. O.P. acknowledges funding from VIA within NOVA. ML was supported
by the John Harvard, ITC and NASA Hubble Fellowship Program fellowships. B.R.
is supported by the Natural Sciences & Engineering Research Council of Canada
(NSERC) and by a grant from the Simons Foundation (MP-SCMPS-00001470).
Research at the Flatiron Institute is supported by the Simons Foundation.

This research was enabled by support provided by grant no. NSF PHY-1125915
along with a INCITE program award PHY129, using resources from the Oak Ridge
Leadership Computing Facility, Summit, which is a US Department of Energy
office of Science User Facility supported under contract DE-AC05- 00OR22725,
as well as Calcul Quebec (http://www.calculquebec.ca) and Compute Canada
(http://www.computecanada.ca). The computational resources and services used
in this work were partially provided by facilities supported by the Scientific
Computing Core at the Flatiron Institute, a division of the Simons Foundation;
and by the VSC (Flemish Supercomputer Center), funded by the Research
Foundation Flanders (FWO) and the Flemish Government – department EWI. This
research is part of the Frontera (Stanzione et al., 2020) computing project at
the Texas Advanced Computing Center (LRAC-AST20008). Frontera is made possible
by National Science Foundation award OAC-1818253. This work used the Dutch
national e-infrastructure with the support of the SURF Cooperative using grant
no. EINF-3036 and EINF-5383, which is (partly) financed by the Dutch Research
Council (NWO), for post-processing of simulation data.

## Data Availability

The simulation post-processed data used to plot the images in this work are
available in Zenodo at http://doi.org/10.5281/zenodo.10996780

## References

  * Acciari et al. (2010) Acciari V. A., et al., 2010, ApJ, 716, 819
  * Aharonian et al. (2006) Aharonian F., et al., 2006, Science, 314, 1424
  * Akiyama et al. (2023) Akiyama K., et al., 2023, ApJ, 957, L20
  * Aliu et al. (2012) Aliu E., et al., 2012, ApJ, 746, 141
  * Anninos et al. (2005) Anninos P., Fragile P. C., Salmonson J. D., 2005, ApJ, 635, 723
  * Asada & Nakamura (2012) Asada K., Nakamura M., 2012, ApJ, 745, L28
  * Baganoff et al. (2001) Baganoff F. K., et al., 2001, Nature, 413, 45
  * Balbus & Hawley (1991) Balbus S. A., Hawley J. F., 1991, Astrophysical Journal, 376, 214
  * Balbus & Hawley (1998) Balbus S. A., Hawley J. F., 1998, Reviews of Modern Physics, 70, 1
  * Begelman et al. (2022) Begelman M. C., Scepi N., Dexter J., 2022, MNRAS, 511, 2040
  * Bhattacharjee et al. (2009) Bhattacharjee A., Huang Y.-M., Yang H., Rogers B., 2009, Physics of Plasmas, 16, 112102
  * Bisnovatyi-Kogan & Ruzmaikin (1974) Bisnovatyi-Kogan G. S., Ruzmaikin A. A., 1974, Ap&SS, 28, 45
  * Blanch (2021) Blanch O., 2021, The Astronomer’s Telegram, 14483, 1
  * Blandford & Znajek (1977) Blandford R. D., Znajek R. L., 1977, Monthly Notices of the Royal Astronomical Society, 179, 433
  * Bransgrove et al. (2021) Bransgrove A., Ripperda B., Philippov A., 2021, Phys. Rev. Lett., 127, 055101
  * Chatterjee & Narayan (2022) Chatterjee K., Narayan R., 2022, ApJ, 941, 30
  * Chatterjee et al. (2019) Chatterjee K., Liska M., Tchekhovskoy A., Markoff S. B., 2019, Monthly Notices of the Royal Astronomical Society, 490, 2200
  * Colella & Woodward (1984) Colella P., Woodward P. R., 1984, JCoPh, 54, 174
  * Davelaar et al. (2018) Davelaar J., Mościbrodzka M., Bronzwaer T., Falcke H., 2018, Astronomy & Astrophysics, 612, A34
  * Davelaar et al. (2023) Davelaar J., et al., 2023, ApJ, 959, L3
  * Del Zanna et al. (2007) Del Zanna L., Zanotti O., Bucciantini N., Londrillo P., 2007, A&A, 473, 11
  * Dexter et al. (2020) Dexter J., et al., 2020, MNRAS, 497, 4999
  * EHT MWL Science Working Group et al. (2021) EHT MWL Science Working Group Algaba J. C., Anczarski J., Asada K., Baloković M., Chandra S., et al. 2021, ApJ, 911, L11
  * Eckart et al. (2004) Eckart A., et al., 2004, A&A, 427, 1
  * Etienne et al. (2015) Etienne Z. B., Paschalidis V., Haas R., Mösta P., Shapiro S. L., 2015, Classical and Quantum Gravity, 32, 175009
  * Event Horizon Telescope Collaboration et al. (2019a) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Asada K., et al. 2019a, ApJ, 875, L1
  * Event Horizon Telescope Collaboration et al. (2019b) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Asada K., et al. 2019b, ApJ, 875, L5
  * Event Horizon Telescope Collaboration et al. (2021) Event Horizon Telescope Collaboration Akiyama K., Algaba J. C., Alberdi A., Alef W., et al. 2021, ApJ, 910, L13
  * Event Horizon Telescope Collaboration et al. (2022a) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Algaba J. C., et al. 2022a, ApJ, 930, L12
  * Event Horizon Telescope Collaboration et al. (2022b) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Algaba J. C., et al. 2022b, ApJ, 930, L13
  * Event Horizon Telescope Collaboration et al. (2022c) Event Horizon Telescope Collaboration Akiyama K., Alberdi A., Alef W., Algaba J. C., et al. 2022c, ApJ, 930, L16
  * Fabian (2012) Fabian A., 2012, Annual Review of Astronomy and Astrophysics, 50, 455
  * Fishbone & Moncrief (1976) Fishbone L. G., Moncrief V., 1976, The Astrophysical Journal, 207, 962
  * Fuentes et al. (2023) Fuentes A., et al., 2023, Nature Astronomy, 7, 1359
  * GRAVITY Collaboration et al. (2018) GRAVITY Collaboration et al., 2018, A&A, 618, L10
  * Gammie (2004) Gammie C. F., 2004, The Astrophysical Journal, 614
  * Gammie et al. (2003) Gammie C. F., McKinney J. C., Tóth G., 2003, ApJ, 589, 444
  * Gardiner & Stone (2005) Gardiner T. A., Stone J. M., 2005, Journal of Computational Physics, 205, 509
  * Goedbloed & Keppens (2022) Goedbloed H., Keppens R., 2022, ApJS, 259, 65
  * Gravity Collaboration et al. (2021) Gravity Collaboration Abuter R., et al., 2021, Astronomy and Astrophysics-A&A, 654, A22
  * Grete et al. (2021) Grete P., Glines F. W., O’Shea B. W., 2021, IEEE Transactions on Parallel and Distributed Systems, 32, 85
  * Guan & Gammie (2008) Guan X., Gammie C. F., 2008, ApJS, 174, 145
  * Hada et al. (2013) Hada K., et al., 2013, ApJ, 775, 70
  * Hakobyan et al. (2023) Hakobyan H., Ripperda B., Philippov A. A., 2023, ApJ, 943, L29
  * Hide & Palmer (1982) Hide R., Palmer T., 1982, Geophysical & Astrophysical Fluid Dynamics, 19, 301
  * Igumenshchev (2008) Igumenshchev I. V., 2008, ApJ, 677, 317
  * Igumenshchev et al. (2003) Igumenshchev I. V., Narayan R., Abramowicz M. A., 2003, ApJ, 592, 1042
  * Janssen et al. (2021) Janssen M., et al., 2021, Nature Astronomy, 5, 1017
  * Jia et al. (2023) Jia H., Ripperda B., Quataert E., White C. J., Chatterjee K., Philippov A., Liska M., 2023, MNRAS, 526, 2924
  * Kim et al. (2018) Kim J. Y., et al., 2018, A&A, 616, A188
  * Liska et al. (2019) Liska M., Tchekhovskoy A., Ingram A., van der Klis M., 2019, MNRAS, 487, 550
  * Liska et al. (2020) Liska M., Tchekhovskoy A., Quataert E., 2020, MNRAS, 494, 3656
  * Liska et al. (2022) Liska M. T. P., et al., 2022, ApJS, 263, 26
  * Liska et al. (2024) Liska M. T. P., Kaaz N., Chatterjee K., Emami R., Musoke G., 2024, ApJ, 966, 47
  * McKinney et al. (2012) McKinney J. C., Tchekhovskoy A., Blandford R. D., 2012, MNRAS, 423, 3083
  * Mertens et al. (2016) Mertens F., Lobanov A. P., Walker R. C., Hardee P. E., 2016, A&A, 595, A54
  * Mościbrodzka & Falcke (2013) Mościbrodzka M., Falcke H., 2013, Astronomy & Astrophysics, 559, L3
  * Najafi-Ziyazi et al. (2024) Najafi-Ziyazi M., Davelaar J., Mizuno Y., Porth O., 2024, MNRAS, 531, 3961
  * Nakamura et al. (2018) Nakamura M., et al., 2018, ApJ, 868, 146
  * Narayan et al. (2003) Narayan R., Igumenshchev I. V., Abramowicz M. A., 2003, PASJ, 55, L69
  * Narayan et al. (2012) Narayan R., Sadowski A., Penna R. F., Kulkarni A. K., 2012, Monthly Notices of the Royal Astronomical Society, 426, 3241
  * Narayan et al. (2022) Narayan R., Chael A., Chatterjee K., Ricarte A., Curd B., 2022, MNRAS, 511, 3795
  * Neilsen et al. (2015) Neilsen J., et al., 2015, The Astrophysical Journal, 799, 199
  * Noble et al. (2006) Noble S. C., Gammie C. F., McKinney J. C., Zanna L. D., 2006, The Astrophysical Journal, 641, 626
  * Pessah et al. (2006) Pessah M. E., Chan C.-K., Psaltis D., 2006, MNRAS, 372, 183
  * Porth & Komissarov (2015) Porth O., Komissarov S. S., 2015, MNRAS, 452, 1089
  * Porth et al. (2017) Porth O., Olivares H., Mizuno Y., Younsi Z., Rezzolla L., Moscibrodzka M., Falcke H., Kramer M., 2017, Computational Astrophysics and Cosmology, 4, 1
  * Porth et al. (2019) Porth O., et al., 2019, The Astrophysical Journal Supplement Series, 243, 26
  * Porth et al. (2021) Porth O., Mizuno Y., Younsi Z., Fromm C. M., 2021, MNRAS, 502, 2023
  * Prather et al. (2021) Prather B., Wong G., Dhruv V., Ryan B., Dolence J., Ressler S., Gammie C., 2021, The Journal of Open Source Software, 6, 3336
  * Ressler et al. (2023) Ressler S. M., White C. J., Quataert E., 2023, MNRAS, 521, 4277
  * Ripperda et al. (2019) Ripperda B., et al., 2019, ApJS, 244, 10
  * Ripperda et al. (2020) Ripperda B., Bacchini F., Philippov A. A., 2020, ApJ, 900, 100
  * Ripperda et al. (2022) Ripperda B., Liska M., Chatterjee K., Musoke G., Philippov A. A., Markoff S. B., Tchekhovskoy A., Younsi Z., 2022, ApJ, 924, L32
  * Schawinski et al. (2007) Schawinski K., Thomas D., Sarzi M., Maraston C., Kaviraj S., Joo S.-J., Yi S. K., Silk J., 2007, Monthly Notices of the Royal Astronomical Society, 382, 1415
  * Shakura & Sunyaev (1973) Shakura N. I., Sunyaev R. A., 1973, Astronomy and Astrophysics, 24, 337
  * Silk & Rees (1998) Silk J., Rees M. J., 1998, Astronomy and Astrophysics, 331, L1
  * Sironi et al. (2021) Sironi L., Rowan M. E., Narayan R., 2021, ApJ, 907, L44
  * Sądowski et al. (2014) Sądowski A., Narayan R., McKinney J. C., Tchekhovskoy A., 2014, MNRAS, 439, 503
  * Stanzione et al. (2020) Stanzione D., West J., Evans R. T., Minyard T., Ghattas O., Panda D. K., 2020, in Practice and Experience in Advanced Research Computing. ACM, New York, NY, USA, pp 106–111, doi:10.1145/3311790.3396656
  * Takahashi et al. (2016) Takahashi H. R., Ohsuga K., Kawashima T., Sekiguchi Y., 2016, ApJ, 826, 23
  * Tchekhovskoy et al. (2011) Tchekhovskoy A., Narayan R., Mckinney J. C., 2011, Monthly Notices of the Royal Astronomical Society: Letters, 418
  * Tchekhovskoy et al. (2012) Tchekhovskoy A., McKinney J. C., Narayan R., 2012, in Journal of Physics: Conference Series. p. 012040
  * Uzdensky et al. (2010) Uzdensky D. A., Loureiro N. F., Schekochihin A. A., 2010, Phys. Rev. Lett., 105, 235002
  * White et al. (2019) White C. J., Stone J. M., Quataert E., 2019, The Astrophysical Journal, 874, 168
  * White et al. (2023) White C. J., Mullen P. D., Jiang Y.-F., Davis S. W., Stone J. M., Morozova V., Zhang L., 2023, ApJ, 949, 103
  * Wong et al. (2021) Wong G. N., Du Y., Prather B. S., Gammie C. F., 2021, ApJ, 914, 55
  * Yuan & Narayan (2014) Yuan F., Narayan R., 2014, ARA&A, 52, 529
  * Yuan et al. (2015) Yuan F., Gan Z., Narayan R., Sadowski A., Bu D., Bai X.-N., 2015, The Astrophysical Journal, 804, 101
  * Zhang et al. (2003) Zhang Y.-T., Shi J., Shu C.-W., Zhou Y., 2003, Phys. Rev. E, 68, 046709
  * Zhdankin et al. (2023) Zhdankin V., Ripperda B., Philippov A. A., 2023, Physical Review Research, 5, 043023


w1+∞ and the Celestial Sphere

Andrew Strominger

Center for the Fundamental Laws of Nature, Harvard University, Cambridge, MA
02138, USA

It is shown that the infinite tower of tree-level plus-helicity soft graviton
symmetries in asymptotically flat 4D quantum gravity can be organized into a
single chiral 2D Kac-Moody symmetry based on the wedge algebra of w1+∞ which
naturally acts on the celestial sphere at null infinity. The infinite towers
of soft photon or gluon symmetries also transform irreducibly under w1+∞.

###### Contents

  1. 1 Introduction
  2. 2 Gravity
  3. 3 Gauge theory
  4. 4 Speculations

##  1 Introduction

A central problem in physics is to find all the fundamental nontrivial
symmetries of nature implied by all the experimentally verified physical laws.
By this we mean symmetries associated via Noether’s therorem to conservation
laws with observable consequences. A hundred years ago the answer to this
question would have been a short list including Poincare symmetries, (which
lead to energy-angular momentum conservation) plus a few more.111We do not
consider here the rich variety of emergent symmetries also observable in
nature.

In the 1960s BMS [1, 2] showed that the answer cannot be so simple because
there is no sense in which the diffeomorphism group of general relativity (GR)
in asymptotically flat spacetimes can be reduced to the Poincare group. They
did not however (due to then-uncertainties about the structure of asymptotic
infinity) either identify an alternate larger asymptotic symmetry group of the
full past and future spacetime or associate any observable conservation laws.
Recently this problem has been translated into the language of quantum field
theory and Feynman diagrams where it becomes equivalent to identifying soft
theorems. A soft theorem is a linear relation between scattering amplitudes
(in asymptotically past and future flat spacetimes) in which one particle
becomes ‘soft” in that its energy is taken to zero. Such theorems can always
be recast as a conservation-law-implying symmetry[3, 4, 5]. Completing the
program of BMS, an infinite number of exact conservation laws were thereby
discovered which relate arbitrary moments of ingoing and outgoing energy-
momentum fluxes to measurable gravitational memory effects [6].

These developments were satisfying but far from the end of the story. Soft
theorems abound in both gauge theory and gravity, with more being discovered
only recently, and each associated to an infinite number of symmetries and
measurable conservation laws. Moreover, the known symmetries do not close
under commutation, implying an infinite tower of soft theorems [7, 8, 9, 10].

Hence despite all the progress, finding all the symmetries and conservation
laws implied by the standard model plus GR remains an outstanding open
problem.

In this letter we solve the problem in the limited context of the tree-level
approximation with vanishing cosmological constant. Moreover we make the
significant restriction to symmetries associated with plus-helicity soft
particles only. Under these conditions we show that the soft symmetries can be
succinctly described by a certain well-known infinite-dimensional $w$-symmetry
group.

We rely heavily on recent progress made in [11] for tree-level gravity and
gauge theory. In this work plus-helicity soft symmetries were compactly
represented by 2D (higher-spin) currents in the celestial conformal field
theory (CCFT)222This paper employs a bottom-up approach in which the CCFT is
simply defined by a Mellin transform of gauge and gravity scattering
amplitudes whose conformal properties are deduced from soft theorems. This
differs from the top-down approach usually employed in the AdS/CFT
correrspondence, where the boundary CFT can be constructed from first
principles using string theory. A top down approach is not available here as
we strive to describe the real world for which the complete microscopic
fundamental laws are unknown. living on the celestial sphere at null infinity.
An infinite tower of such currents and their algebra were derived at tree
level using positive-helicity soft theorems and the celestial operator product
expansion (OPE).

The somewhat lengthy results of [11] were displayed in a basis with manifest
covariance under the global Lorentz-conformal group of the celestial sphere.
The currents include an ‘$SL(2,{\mathbb{R}})_{\rm w}$’ Kac-Moody algebra
arising from the subleading soft graviton theorem. In this paper we find
dramatic simplifications by reorganizing the generators according to
representations of this $SL(2,{\mathbb{R}})_{\rm w}$, which is accomplished
with a version of the light ray transform [12]. The entire tower of currents
is assembled into a chiral Kac-Moody symmetry of the wedge algebra of w1+∞!
This same Kac-Moody algebra has appeared previously including in the Penrose
twistor construction [13], in discrete states of the $c=1$ string [14] and
w1+∞-gravity[15]. Moreover we find, unlike in [11], in the
$SL(2,{\mathbb{R}})_{\rm w}$ covariant presentation only positive half-
integral weights appear.

In the next section we present explicity the algebraic transformation from the
conventional to the $SL(2,{\mathbb{R}})_{\rm w}$ covariant basis, and show
that the entire algebra is the Kac-Moody symmetry of the wedge algebra of
w1+∞. In section 3 the infinite tower of soft symmetries for gauge theory are
written in $SL(2,{\mathbb{R}})_{\rm w}$ covariant form and again found to
dramatically simplify. Again only positive $SL(2,{\mathbb{R}})_{\rm w}$
weights appear. These results suggest that w1+∞ \- or perhaps its quantization
W1+∞ [15]\- provides an organizing principle for CCFT. We conclude with
speculations on this role in section 4.

##  2 Gravity

Let us recap the basic results and notation of [11] for gravity. Let
$G_{\Delta}^{+}(z,\bar{z})$ denote the positive-helicity conformal-primary
graviton operator with 2D conformal weight $\Delta$ which crosses the
celestial sphere at a point $(z,{\bar{z}})$. Define a discrete family of
conformally soft333The conformally soft gravitons defined here differ from the
usual energetically soft gravitons in that the conformal weight, rather than
the enegy, is taken to a limiting value. Nevertheless they obey soft theorems
which mirror their energetically-defined counterparts [16].positive-helicity
gravitons

$H^{k}=\lim_{\varepsilon\to 0}\varepsilon G_{k+\varepsilon}^{+},\ \ \ \
k=2,1,0,-1,\ldots,$ (2.1)

with weights

$(h,\bar{h})=\left(\frac{k+2}{2},\frac{k-2}{2}\right),$ (2.2)

and a consistently-truncated antiholomorphic mode expansion444Outside the
specified range of $n$, the $SL(2,{\mathbb{R}})_{R}$-invariant norm vanishes.
Such operators may still have contact interactions but in this paper operators
are always at distinct points.

$\displaystyle H^{k}(z,\bar{z})$
$\displaystyle=\sum_{n=\frac{k-2}{2}}^{\frac{2-k}{2}}\frac{H^{k}_{n}(z)}{\bar{z}^{n+\frac{k-2}{2}}}.$
(2.3)

Each $H^{k}_{n}(z)$ is a 2D symmetry-generating conserved current whose Ward
identity is a soft theorem. The factor of $\varepsilon$ in (2.1) cancels the
soft poles, allowing a finite OPE. Throughout this paper $z$ and $\bar{z}$ are
treated as independent, which amounts to continuing (3,1) Minkowski space to
(2,2) Klein space, the celestial sphere to the celestial torus [17] and
Lorentz-$SL(2,\mathbb{C})$ to $SL(2,\mathbb{R})_{L}\otimes
SL(2,\mathbb{R})_{R}$. The $n$-index in (2.3) then transforms in the $3-k$
dimensional representation of the $SL(2,{\mathbb{R}})_{R}$. The simplest
example of the $k=1$ term generates supertranslations. Expanding
$H^{1}_{\pm{1\over 2}}(z)=\sum_{m}H^{1}_{\pm{1\over 2},m}z^{-m-3/2}$, the four
modes $H^{1}_{\pm{1\over 2},\pm{1\over 2}}$ generate the four global
translations. $k=0$ is related to superrotations and includes Lorentz
transformations. Defining the commutator for holomorphic objects555Note that
this is a 2D celestial commutator on a 1D circle, not to be mistaken for a 4D
commutator on a 3D slice [18].

$\left[A,B\right](z)=\oint_{z}\frac{dw}{2\pi i}A(w)B(z),$ (2.4)

the soft current algebra found in [11] for gravity is666We note that the
$(2,0)$ current $H^{2}_{0}(z)$ commutes with all other generators and is a
central term in the supertranslation algebra. It has been taken to vanish in
most applications but is natural to include here.

$\left[H^{k}_{m},H^{l}_{n}\right]=-\frac{\kappa}{2}\left[n(2-k)-m(2-l)\right]\frac{(\frac{2-k}{2}-m+\frac{2-l}{2}-n-1)!}{(\frac{2-k}{2}-m)!(\frac{2-l}{2}-n)!}\frac{(\frac{2-k}{2}+m+\frac{2-l}{2}+n-1)!}{(\frac{2-k}{2}+m)!(\frac{2-l}{2}+n)!}H^{k+l}_{m+n},$
(2.5)

To write the algebra (2.5) in a simpler form we define

$w^{p}_{n}={1\over\kappa}(p-n-1)!(p+n-1)!H^{-2p+4}_{n}.$ (2.6)

This is essentially the light-transform or right-shadow [12] in one dimension
adapted to finite $SL(2,{\mathbb{R}})$ representations, and $p$ is essentially
the righ-shadowed weight. (2.5) then becomes

$[w^{p}_{m},w^{q}_{n}]=\left[m(q-1)-n(p-1)\right]w^{p+q-2}_{m+n}.$ (2.7)

Since the index $k$ on $H^{k}_{m}$ runs over $k=2,1,0,....$, $p$ runs over the
positive half integral values

$p=1,{3\over 2},2,{5\over 2},....$ (2.8)

The restriction ${k-2\over 2}\leq m\leq{2-k\over 2}$ becomes777For example for
the $p=2$ Virasoro case this restricts to the $SL(2,\mathbb{R})$ current
subalgebra.

$1-p\leq m\leq p-1.$ (2.9)

This is one of our main results.

The commutators (2.7) were first written down by Ioannis Bakas in 1989 [19].
In this work $m$ is an arbitrary integer and the resulting algebra is now
referred to as w1+∞. The closed $p=2$ subalgebra of w1+∞ is the $c=0$ Virasoro
algebra. This algebra in (2.7) is a subgroup of w1+∞ with $m$ restricted
according to (2.9). The resulting restricted algebra is known as the wedge
subalgebra of w1+∞, and can also be rewritten as
$GL(\infty,{\mathbb{R}})$[15].

Moreover, each $w_{m}^{p}(z)$ acts as a current on the celestial sphere. Hence
we have a $GL(\infty,{\mathbb{R}})$ Kac-Moody algebra. The modes of the Kac-
Moody currents act on the states of the 2D CCFT associated by the state-
operator correspondence to operator insertions in the celestial sphere [18].
In fact precisely this Kac-Moody algebra has been studied previously both as
the symmetry group of $c=1$ string theory [14] and 2D w1+∞-gravity
[20].888Here we take the classical limit W1+∞ $\to$w1+∞ of the quantum
symmetry W1+∞. A review of this and other fascinating aspects of W algebras
can be found in [21]. We return to this connection in the last section.

The closed subgroup generated by

$L_{m}\equiv w_{m}^{2},~{}~{}~{}m=-1,0,1$ (2.10)

is an ‘$SL(2,{\mathbb{R}})_{\rm w}$’ current algebra implied by the subleading
soft graviton theorem. It is related to the original generators in (2.5) by

$L_{0}={1\over\kappa}H^{0}_{0},~{}~{}L_{\pm 1}={2\over\kappa}H^{0}_{\pm 1}.$
(2.11)

$w^{q}_{n}$ transforms under this $SL(2,{\mathbb{R}})_{\rm w}$ as

$[L_{m},w^{q}_{n}]=\left[m(q-1)-n\right]w^{q}_{m+n}.$ (2.12)

like the modes of an $SL(2,{\mathbb{R}})_{\rm w}$ primary operator of
$SL(2,{\mathbb{R}})_{\rm w}$ weight $q$. Since $q\geq 1$ here, these are all
positive. However, unlike the $H^{k}_{n}$, they do not transform canonically
under the original $SL(2,{\mathbb{R}})_{R}$ because of the mode-dependent
relation between the two. Indeed the normalization of the $H^{k}_{n}$ modes
was chosen so that they transform canonically under the original
$SL(2,{\mathbb{R}})_{R}$ conformal generators with weights
$\bar{h}=\frac{k-2}{2}$.

Evidently $SL(2,{\mathbb{R}})_{R}$ and $SL(2,{\mathbb{R}})_{\rm w}$ are not
quite the the same thing. $w_{n}^{p}$ on the LHS of (2.6) lies in a positive
weight $h_{w}=p$ representation of $SL(2,{\mathbb{R}})_{\rm w}$, while
$H_{n}^{-2p+4}$ on the RHS lies in a $2p+1$ dimensional, negative weight
$\bar{h}=1-p$ representation of $SL(2,{\mathbb{R}})_{R}$. (2.6) is the
relation between them. Both representations are $2p+1$ dimensional, because
finite weight $h$ $SL(2,{\mathbb{R}})$ representations are $2h-1$ dimensional
for positive half integral $h$ and $-2h-1$ dimensional for negative half
integral $h$.

Significantly, the $SL(2,{\mathbb{R}})_{\rm w}$ representations which appear
here are all positive weight.

##  3 Gauge theory

Now we turn to nonabelian gauge theory. Let ${O}_{k}^{a,+}(z,\bar{z})$ denote
a positive helicity, conformal weight $k$ gluon operator with adjoint group
index $a$ at the point $(z,\bar{z})$ on the celestial sphere. Mode-expanding
in $\bar{z}$ on the right

${O}_{k}^{a,+}(z,\bar{z})=\sum_{n}\frac{{O}^{a,+}_{k,n}(z)}{\bar{z}^{n+{k-1\over
2}}},$ (3.1)

conformally soft currents are defined by

$R^{k,a}_{n}(z):=\lim_{\varepsilon\to 0}\varepsilon
O_{k+\varepsilon,n}^{a,+}(z),~{}~{}~{}~{}~{}\ \ \ \
k=1,0,-1,-2,\ldots,~{}~{}~{}~{}~{}~{}\quad{k-1\over 2}\leq n\leq{1-k\over 2}.$
(3.2)

This has $SL(2,{\mathbb{R}})$${}_{L}\otimes$$SL(2,{\mathbb{R}})_{R}$ weights

$(h,\bar{h})=\left(\frac{k+1}{2},\frac{k-1}{2}\right).$ (3.3)

These values of conformal weights $\Delta=k$ include all the conformally soft
poles encountered in the OPE [22, 23]. The factor of $\varepsilon$
incorporated in (3.2) is needed to cancel these poles, leading to finite OPEs
for the rescaled $R^{k,a}$.

The soft current algebra for gauge theory is [11]

$\left[R^{k,a}_{n},R^{l,b}_{n^{\prime}}\right]=-if^{ab}{}_{c}\frac{(\frac{1-k}{2}-n+\frac{1-l}{2}-n^{\prime})!}{(\frac{1-k}{2}-n)!(\frac{1-l}{2}-n^{\prime})!}\frac{(\frac{1-k}{2}+n+\frac{1-l}{2}+n^{\prime})!}{(\frac{1-k}{2}+n)!(\frac{1-l}{2}+n^{\prime})!}R^{k+l-1,c}_{n+n^{\prime}}.$
(3.4)

Let us define

$S^{q,a}_{m}=(q-m-1)!(q+m-1)!R^{3-2q,a}_{m},$ (3.5)

where $q=1,{3\over 2},2,...$. One then finds the simple algebra

$\left[S^{q,a}_{n},S^{p,b}_{n^{\prime}}\right]=-if^{ab}{}_{c}S^{q+p-1,c}_{n+n^{\prime}}.$
(3.6)

Moreover using [11]

$\left[H^{k}_{m},R^{l,a}_{n}\right]=-\frac{\kappa}{2}\left[n(2-k)-m(1-l)\right]\frac{(\frac{2-k}{2}-m+\frac{1-l}{2}-n-1)!}{(\frac{2-k}{2}-m)!(\frac{1-l}{2}-n)!}\frac{(\frac{2-k}{2}+m+\frac{1-l}{2}+n-1)!}{(\frac{2-k}{2}+m)!(\frac{1-l}{2}+n)!}R^{k+l,a}_{m+n},$
(3.7)

one finds the irreducible representation

$\left[w^{p}_{m},S^{q,a}_{n}\right]=\left[m(q-1)-n(p-1)\right]S^{p+q-2,a}_{m+n}.$
(3.8)

The $SL(2,{\mathbb{R}})_{\rm w}$ transformation of $S^{q}_{m}$

$\left[L_{m},S^{q,a}_{n}\right]=\left[m(q-1)-n\right]S^{q,a}_{m+n}.$ (3.9)

is that of a primary of weight $q$. Again we find only positive weights.

##  4 Speculations

The appearance of W-algebras in the 2D celestial symmetry group connects
celestial holography to several other research areas. We find it irresistible
to speculate on what might lie ahead as these connections unfold.

w1+∞ has a natural deformation to W1+∞, a significantly more complicated
algebra. This deformation can be understood as arising from quantization. In
the context of specific classical 2D realizations of w1+∞, anomalies
encountered in quantization deform the classical w1+∞ algebra to the quantum
W1+∞[24]. While currently not well-understood, the tree-level soft algebra of
4D quantum gravity also gets deformed on quantization, as implied among other
things by one-loop corrections to the soft theorems. So it is natural to
speculate that the action of w1+∞ on this soft algebra is deformed to W1+∞ in
the 4D quantum theory of gravity.

One might have expected a Virasoro algebra rather than an
$SL(2,{\mathbb{R}})_{\rm w}$ current algebra among the symmetries. However the
two are closely related in many similar contexts. For example the Chern-Simon
formulation of AdS3 gravity, at first sight gives an
$SL(2,{\mathbb{R}})_{R}$$\otimes SL(2,{\mathbb{R}})_{L}$ current algebra on
the boundary. However, the AdS3 boundary conditions implement a Hamiltonian
reduction to Vir${}_{R}\otimes$VirL [25, 26].999 A similar reduction could
perhaps be operative here from constraints related to IR divergences. The
Virasoro generators are field-dependent $SL(2,{\mathbb{R}})$ transformations
which lie in the enveloping algebra of the $SL(2,{\mathbb{R}})_{R}$$\otimes
SL(2,{\mathbb{R}})_{L}$ current algebra. In another example, 2D gravity in
light cone gauge exhibits an $SL(2,{\mathbb{R}})_{R}$ current algebra [27],
but the same theory in conformal gauge exhibits a Virasoro symmetry [28],
again indicating a relation between the two symmetry actions.

This type of relation has an infinite-dimensional uplift to the present
context. $SL(2,{\mathbb{R}})$ is the wedge algebra of Virasoro, just as
$GL(\infty,{\mathbb{R}})$ is the wedge algebra of w1+∞[15]. Just as Virasoro
lies in the enveloping algebra of $SL(2,{\mathbb{R}})$ Kac-Moody, w1+∞ (or
W1+∞ at the quantum level) lies in the enveloping algebra of
$GL(\infty,{\mathbb{R}})$ Kac-Moody [20]. So it is natural to speculate that
W1+∞ is a celestial symmetry of 4D quantum gravity.

w1+∞ also appears as the symmetry group of classical self-dual gravity in
(2,2) signature Klein space [29], whose sole degree of freedom is given by the
Kahler potential. The self-duality suggests that the dual CCFT should be
chiral, and acted on (at bulk tree level) by w1+∞. This is a natural context
in which to study celestial holography, because the complicated interactions
of opposite chirality sectors and double soft-limit ambiguities are absent.
The $N=2$ string describes the quantization of this theory [30], where w1+∞ is
potentially deformed to W1+∞. Perhaps the holographic CCFT dual of $N=2$
string theory in 4D Klein space is given by a 2D quantum W1+∞-gravity theory
on the celestial torus.

We leave these thoughts to future explorations.

## Acknowledgements

I am grateful to Ofer Aharony, Adam Ball, Alfredo Guevara, Mina Himwich,
Sruthi Narayanan and Monica Pate for useful conversations. This work was
supported by DOE grant de-sc/0007870.

## References

  * [1] H. Bondi, M. van der Burg, and A. Metzner, “Gravitational Waves in General Relativity. 7. Waves from Axisymmetric Isolated Systems,” Proc. Roy. Soc. Lond. A 269 (1962) 21–52.
  * [2] R. Sachs, “Gravitational Waves in General Relativity. 8. Waves in Asymptotically Flat Space-times,” Proc. Roy. Soc. Lond. A 270 (1962) 103–126.
  * [3] A. Strominger, “Asymptotic Symmetries of Yang-Mills Theory,” JHEP 07 (2014) 151, arXiv:1308.0589 [hep-th].
  * [4] A. Strominger, “On BMS Invariance of Gravitational Scattering,” JHEP 07 (2014) 152, arXiv:1312.2229 [hep-th].
  * [5] T. He, V. Lysov, P. Mitra, and A. Strominger, “BMS Supertranslations and Weinberg’s Soft Graviton Theorem,” JHEP 05 (2015) 151, arXiv:1401.7026 [hep-th].
  * [6] A. Strominger and A. Zhiboedov, “Gravitational Memory, BMS Supertranslations and Soft Theorems,” JHEP 01 (2016) 086, arXiv:1411.5745 [hep-th].
  * [7] Y. Hamada and G. Shiu, “Infinite set of soft theorems in gauge-gravity theories as ward-takahashi identities,” Physical Review Letters 120 no. 20, (May, 2018) . http://dx.doi.org/10.1103/PhysRevLett.120.201601.
  * [8] Z.-Z. Li, H.-H. Lin, and S.-Q. Zhang, “Infinite soft theorems from gauge symmetry,” Physical Review D 98 no. 4, (Aug, 2018) . http://dx.doi.org/10.1103/PhysRevD.98.045004.
  * [9] S. Banerjee, S. Ghosh, and P. Paul, “Mhv graviton scattering amplitudes and current algebra on the celestial sphere,” Journal of High Energy Physics 2021 no. 2, (Feb, 2021) . http://dx.doi.org/10.1007/JHEP02(2021)176.
  * [10] S. Pasterski, A. Puhm, and E. Trevisani, “Celestial Diamonds: Conformal Multiplets in Celestial CFT,” arXiv:2105.03516 [hep-th].
  * [11] A. Guevara, E. Himwich, M. Pate, and A. Strominger, “Holographic Symmetry Algebras for Gauge Theory and Gravity,” arXiv:2103.03961 [hep-th].
  * [12] P. Kravchuk and D. Simmons-Duffin, “Light-ray operators in conformal field theory,” JHEP 11 (2018) 102, arXiv:1805.00098 [hep-th].
  * [13] R. Penrose, “Nonlinear Gravitons and Curved Twistor Theory,” Gen. Rel. Grav. 7 (1976) 31–52.
  * [14] I. R. Klebanov and A. M. Polyakov, “Interaction of discrete states in two-dimensional string theory,” Mod. Phys. Lett. A 6 (1991) 3273–3281, arXiv:hep-th/9109032.
  * [15] C. N. Pope, L. J. Romans, and X. Shen, “$W$(infinity) and the Racah-wigner Algebra,” Nucl. Phys. B 339 (1990) 191–221.
  * [16] M. Pate, A.-M. Raclariu, and A. Strominger, “Conformally Soft Theorem in Gauge Theory,” Phys. Rev. D 100 no. 8, (2019) 085017, arXiv:1904.10831 [hep-th].
  * [17] A. Atanasov, A. Ball, W. Melton, A.-M. Raclariu, and A. Strominger, “$(2,2)$ Scattering and the Celestial Torus,” arXiv:2101.09591 [hep-th].
  * [18] E. Crawley, N. Miller, S. A. Narayanan, and A. Strominger, “State-Operator Correspondence in Celestial Conformal Field Theory,” arXiv:2105.00331 [hep-th].
  * [19] I. Bakas, “The Large n Limit of Extended Conformal Symmetries,” Phys. Lett. B 228 (1989) 57.
  * [20] C. N. Pope, X. Shen, K. W. Xu, and K.-j. Yuan, “SL(infinity,R) symmetry of quantum W(infinity) gravity,” Nucl. Phys. B 376 (1992) 52–74.
  * [21] C. N. Pope, “Lectures on W algebras and W gravity,” in Summer School in High-energy Physics and Cosmology.  12, 1991.  arXiv:hep-th/9112076.
  * [22] W. Fan, A. Fotopoulos, and T. R. Taylor, “Soft Limits of Yang-Mills Amplitudes and Conformal Correlators,” JHEP 05 (2019) 121, arXiv:1903.01676 [hep-th].
  * [23] M. Pate, A.-M. Raclariu, A. Strominger, and E. Y. Yuan, “Celestial Operator Products of Gluons and Gravitons,” arXiv:1910.07424 [hep-th].
  * [24] E. Bergshoeff, P. S. Howe, C. N. Pope, E. Sezgin, X. Shen, and K. S. Stelle, “Quantization deforms w(infinity) to W(infinity) gravity,” Nucl. Phys. B 363 (1991) 163–184.
  * [25] M. Banados, “Global charges in Chern-Simons field theory and the (2+1) black hole,” Phys. Rev. D 52 (1996) 5816–5825, arXiv:hep-th/9405171.
  * [26] O. Coussaert, M. Henneaux, and P. van Driel, “The Asymptotic dynamics of three-dimensional Einstein gravity with a negative cosmological constant,” Class. Quant. Grav. 12 (1995) 2961–2966, arXiv:gr-qc/9506019.
  * [27] A. M. Polyakov, “Quantum Gravity in Two-Dimensions,” Mod. Phys. Lett. A 2 (1987) 893.
  * [28] F. David, “Conformal Field Theories Coupled to 2D Gravity in the Conformal Gauge,” Mod. Phys. Lett. A 3 (1988) 1651.
  * [29] C. P. Boyer and J. F. Plebanski, “An infinite hierarchy of conservation laws and nonlinear superposition principles for selfdual Einstein spaces,” J. Math. Phys. 26 (1985) 229–234.
  * [30] H. Ooguri and C. Vafa, “Geometry of N=2 strings,” Nucl. Phys. B 361 (1991) 469–518.