# Bridging Scales in Black Hole Accretion and Feedback: Magnetized Bondi
Accretion in 3D GRMHD

Hyerin Cho (조혜린) Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden
Street, Cambridge, MA 02138, USA  Black Hole Initiative at Harvard University,
20 Garden Street, Cambridge, MA 02138, USA  Ben S. Prather CCS-2, Los Alamos
National Laboratory, PO Box 1663, Los Alamos, NM 87545, USA  Ramesh Narayan
Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA  Black Hole Initiative at Harvard University, 20
Garden Street, Cambridge, MA 02138, USA  Priyamvada Natarajan Black Hole
Initiative at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA
Department of Astronomy, Yale University, Kline Tower, 266 Whitney Avenue, New
Haven, CT 06511, USA  Department of Physics, Yale University, P.O. Box 208121,
New Haven, CT 06520, USA  Kung-Yi Su Center for Astrophysics $|$ Harvard &
Smithsonian, 60 Garden Street, Cambridge, MA 02138, USA  Black Hole Initiative
at Harvard University, 20 Garden Street, Cambridge, MA 02138, USA  Angelo
Ricarte Center for Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street,
Cambridge, MA 02138, USA  Black Hole Initiative at Harvard University, 20
Garden Street, Cambridge, MA 02138, USA  Koushik Chatterjee Center for
Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA  Black Hole Initiative at Harvard University, 20 Garden Street, Cambridge,
MA 02138, USA

###### Abstract

Fueling and feedback couple supermassive black holes (SMBHs) to their host
galaxies across many orders of magnitude in spatial and temporal scales,
making this problem notoriously challenging to simulate. We use a multi-zone
computational method based on the general relativistic magneto-hydrodynamic
(GRMHD) code KHARMA that allows us to span $7$ orders of magnitude in spatial
scale, to simulate accretion onto a non-spinning SMBH from an external medium
with Bondi radius $R_{B}\approx 2\times 10^{5}\,GM_{\bullet}/c^{2}$, where
$M_{\bullet}$ is the SMBH mass. For the classic idealized Bondi problem,
spherical gas accretion without magnetic fields, our simulation results agree
very well with the general relativistic analytic solution. Meanwhile, when the
accreting gas is magnetized, the SMBH magnetosphere becomes saturated with a
strong magnetic field. The density profile varies as $\sim r^{-1}$ rather than
$r^{-3/2}$ and the accretion rate $\dot{M}$ is consequently suppressed by over
2 orders of magnitude below the Bondi rate $\dot{M}_{B}$. We find continuous
energy feedback from the accretion flow to the external medium at a level of
$\sim 10^{-2}\dot{M}c^{2}\sim 5\times 10^{-5}\dot{M}_{B}c^{2}$. Energy
transport across these widely disparate scales occurs via turbulent convection
triggered by magnetic field reconnection near the SMBH. Thus, strong magnetic
fields that accumulate on horizon scales transform the flow dynamics far from
the SMBH and naturally explain observed extremely low accretion rates compared
to the Bondi rate, as well as at least part of the energy feedback.

Accretion (14), Active galactic nuclei (16), Bondi accretion (174),
Schwarzschild black holes (1433), Supermassive black holes (1663),
Magnetohydrodynamical simulations (1966)

##  1 Introduction

Most nearby galaxies are found to harbor central supermassive black holes
(SMBHs) whose masses are correlated with properties of the stellar components
of their hosts (e.g., Magorrian et al., 1998; Ferrarese & Merritt, 2000;
Gebhardt et al., 2000; Kormendy & Ho, 2013). However, the details of how gas
flows into the galactic nucleus from large cosmic distances and how the SMBH
in turn imparts feedback into the galaxy remain coupled unresolved problems.

Due to limitations in computational power and resolution, typically there has
been a partitioning of scales in tackling the SMBH feeding and the feedback
problem. In galaxy-scale simulations, the gas flows from kpc scales and the
co-evolution of SMBH and galaxy are explicitly modeled (e.g., Sijacki et al.,
2015; Rosas-Guevara et al., 2016; Weinberger et al., 2018; Ricarte et al.,
2019; Ni et al., 2022; Wellons et al., 2023). However, even in idealized
simulations of an isolated galactic nucleus, prescriptions need to be adopted
as “sub-grid models” for including accretion and feedback (e.g., Li & Bryan,
2014; Fiacconi et al., 2018; Anglés-Alcázar et al., 2021; Talbot et al., 2021;
Weinberger et al., 2023) usually at a scale $\lesssim 10$ pc, which translates
to $\sim 3\times 10^{4}r_{g}$ for the $\sim 6.5\times 10^{9}\,M_{\odot}$ SMBH
in M87; here $r_{g}=GM_{\bullet}/c^{2}$ is the gravitational radius, where
$M_{\bullet}$ is the SMBH mass. Meanwhile, probing smaller scales, general
relativistic magnetohydrodynamic (GRMHD) simulations focus on the inner few
tens or hundreds of $r_{g}$, self-consistently tracing gas flows and feedback
on these limited spatial scales (e.g., Gammie et al., 2003; Tchekhovskoy et
al., 2011; Porth et al., 2019; Chatterjee et al., 2023). Typically, GRMHD
simulations use idealized settings and disregard the larger-scale cosmological
environment. Bridging these vastly different scales while self-consistently
following gas flows and incorporating both the GRMHD effects and cosmological
evolution remains challenging.

Recently, some attempts to connect these disparate scales have been made using
nested simulations where each smaller-scale simulation is initialized from the
next larger-scale simulation (e.g., Hopkins & Quataert, 2010; Ressler et al.,
2020; Guo et al., 2023); or by adopting Lagrangian hyper-refinement methods
(e.g., Anglés-Alcázar et al., 2021; Hopkins et al., 2023); or by pushing out
the simulation regime to larger scales (e.g., Lalakos et al., 2022; Kaaz et
al., 2023). While successful in studying the process of accretion, feedback
from the SMBH has not been followed out to galaxy scales in prior work because
either the communication between scales is only directed inwards or because of
the inability to include the entire galaxy and the black hole.

We employ here a multi-zone computational method which represents a first
attempt to not only follow accretion flows but also trace feedback across 7
orders of magnitude from the event horizon to galactic scales. With this
method, we study purely hydrodynamic Bondi accretion as well as the magnetized
Bondi problem.

The outline of this Letter is as follows: In Section 2 we give an overview of
the numerical scheme that captures a wide dynamic range simultaneously. In
Section 3, we present a purely hydrodynamic simulation and test our numerical
scheme by comparing against the general relativistic (GR) analytical solution.
In Section 4 we include magnetic fields for a more realistic representation of
the environment of SMBHs and study how Bondi accretion is modified. We discuss
our findings of feedback via convection in Section 5 and summarize conclusions
of our study in Section 6. Additional details are presented in a set of short
appendices: various GRMHD quantities are defined in Appendix A, the details of
the simulation set-up are outlined in Appendix B, and a resolution and initial
condition study is presented in Appendix C.

Throughout, we adopt units where $GM_{\bullet}=c=1$, so our unit of length
$r_{g}=1$ and time $t_{g}\equiv r_{g}/c=1$. Gas temperature $T$ is a
dimensionless quantity which is normalized by $\mu c^{2}/k_{B}$, where $\mu$
is the mean molecular weight and $k_{B}$ is the Boltzmann constant; thus gas
pressure $p_{g}=\rho T$ with $\rho$ the density. We denote the time average of
a quantity $X$ as $\overline{X}$ and the density-weighted, time- and shell-
averaged value of $X$ as $\langle X\rangle$ (Appendix A). Finally, the Bondi
radius is $R_{B}\equiv GM_{\bullet}/c_{s,\infty}^{2}$, where
$c_{s,\infty}=\sqrt{\gamma_{\rm ad}T_{\infty}}$ is the sound speed in the
external medium far outside $R_{B}$, $T_{\infty}$ is the temperature there,
and $\gamma_{\rm ad}$ is the gas adiabatic index. The free-fall time at the
Bondi radius is $t_{B}\equiv(R_{B}/r_{g})^{3/2}t_{g}$. For the models
presented here, $\gamma_{\rm ad}=5/3$, $T_{\infty}=3.4\times 10^{-6}$
($\approx 2\times 10^{7}\,{\rm K}$), $R_{B}={1.8\times 10^{5}}\,r_{g}$, and we
consider a non-spinning SMBH (Schwarzschild metric).

##  2 Numerical Methods

Our numerical scheme utilizes the GRMHD code KHARMA111https://github.com/AFD-
Illinois/kharma, a performance-portable C++ implementation based on iharm3D
(Prather et al., 2021); iharm3D is itself an extension of HARM, an efficient
second-order conservative finite-volume scheme for solving MHD equations on
Eulerian meshes in stationary curved space-times (Gammie et al., 2003). KHARMA
offers a more flexible, portable, and scalable implementation suitable for
multiple uses, by leveraging the Parthenon Adaptive Mesh Refinement Framework
and the Kokkos programming model (Grete et al., 2023; Trott et al., 2022).

The approach adopted thus far to handle a wide range of spatial and temporal
scales has been to run a series of nested simulations sequentially, starting
from the outside and proceeding down to the BH (e.g., Ressler et al., 2020;
Guo et al., 2023). This set-up is appropriate for studying the feeding of the
BH from the external galactic medium (or stellar winds), but it cannot tackle
feedback from the BH to the galaxy. To handle both processes, we run nested
simulations from radii $r\gg R_{B}$ down to $r_{g}$ as well as from $r_{g}$
out to $r\gg R_{B}$. In addition, we cycle in and out hundreds of times to
ensure self-consistent, two-way communication between the BH and the external
medium.

We use a grid based on spherical coordinates and split the simulation volume
between $r=r_{g}$ and $r\gg R_{B}$ into a set of $n$ overlapping annuli evenly
spaced in $\log(r)$, each of which we label as zone-$i$. For
$i\in\\{0,1,...,(n-1)\\}$, zone-$i$ extends from an inner radius $r_{i,{\rm
in}}=8^{i}\,r_{g}$ to an outer radius $r_{i,{\rm out}}=8^{i+2}\,r_{g}$. We
start by initializing and running zone-$(n-1)$, the largest annulus, for some
fraction of its free-fall time. We then simulate zone-$(n-2)$, inheriting the
final state of zone-$(n-1)$ as the initial condition over half the new domain.
This continues inward to zone-0 and back out in a sort of “V-cycle.” Our
method can be thought of as defining a global domain, but “pausing” everything
outside the particular active region, shifting the active region up and down
the range of scales a large number of times (in our fiducial MHD run, the
V-cycle was traversed 254 times) until the full system converges to a steady
state.222A few details of the numerical method are given in Section B. A
complete description is deferred to a future paper.

##  3 Hydrodynamic Bondi accretion

Figure 1:  Time-averaged radial profiles of selected quantities for our HD
(red) and MHD (black) Bondi accretion simulations with $R_{B}={1.8\times
10^{5}}\,r_{g}$ (gray vertical dashed line). In this and forthcoming plots,
negative values are plotted with dotted rather than solid lines. While most
naturally expressed in gravitational units (lower x-axis), we also present
distances scaled to M87 (upper x-axis). Where available, analytic solutions of
relativistic Bondi accretion are plotted with thick gray lines. Magnetic
fields lead to a shallower density slope, a dramatically suppressed accretion
rate, and positive energy feedback.

We first test the accuracy of our method by numerically simulating
hydrodynamic (HD) spherical Bondi (1952) accretion on a Schwarzschild BH; a
relativistic analytical solution is available for this problem (Michel, 1972;
Shapiro & Teukolsky, 1983). We use a total of $n=7$ annuli with an outer
radius of $8^{8}\sim 2\times 10^{7}\,r_{g}$ and a total runtime of $t_{\rm
run}(R_{B})\sim 9\,t_{B}$ ($\sim 0.7\,{\rm Myr}$ for M87) at the Bondi radius.
We confirmed convergence by checking that the results at $t_{\rm
run}=9\,t_{B}$ remain unchanged even when we run up to $t_{\rm
run}>56\,t_{B}$. The resolution is $128^{3}$ per annulus giving a total
simulation resolution of $512\times 128\times 128$. Details of the set-up are
presented in Appendix B.1.

Figures 1(a,b) show the radial profiles of $t,\theta,\varphi$-averaged gas
density $\rho$ and temperature $T$ (red lines). The simulation results show
excellent agreement with the GR Bondi analytic solution (thick gray line) over
7 orders of magnitude in radius. Further, the time-averaged accretion rate
$\overline{\dot{M}(r)}$ in Figure 1(c) matches the analytical Bondi accretion
rate $\dot{M}_{B}$ very well, $\overline{\dot{M}}/\dot{M}_{B}=1$, across all
radii within $R_{B}$, confirming a steady state. We note that Guo et al.
(2023) have previously simulated Bondi accretion under Newtonian gravity
spanning 5 orders of magnitude. Here we have successfully reproduced the GR
version over a larger range of scales.

Analogous to $\overline{\dot{M}(r)}$ we can define a time-averaged energy
inflow rate towards the BH, $\overline{\dot{E}(r)}$ (Equation (A5)), but this
is not useful for our purposes since it includes the rest mass energy of the
accreted gas. Instead we consider $\dot{E}_{\rm
net}(r)\equiv\overline{(\dot{M}-\dot{E})}$ (Equation (A6)), which removes the
rest mass contribution and is flipped in sign such that it is the net rate of
outflow of energy (i.e., feedback) from the BH to large radii. Normalizing by
$\overline{\dot{M}}_{10}\equiv\overline{\dot{M}}(r=10\,r_{g})$, we then obtain
a feedback efficiency $\eta^{\rm tot}$ (Equation (A7))

$\eta^{\rm tot}\equiv(\overline{\dot{M}-\dot{E}})/\overline{\dot{M}}_{10}\,.$
(1)

Following Appendix A, and noting that there is no magnetic field and the
system is spherically symmetric, $\eta^{\rm tot}\approx 1+(1+\gamma_{\rm
ad}u/\rho)u_{t}$. At large radii, $u\to\rho T_{\infty}/(\gamma_{\rm ad}-1)$
and $u_{t}\to-1$, hence

$\eta^{\rm tot}\rightarrow-\gamma_{\rm ad}T_{\infty}/(\gamma_{\rm
ad}-1)=-1.5/R_{B}\approx-8.4\times 10^{-6}.$ (2)

The agreement of the simulation results in Figure 1(d) with this theoretical
expectation is excellent.

##  4 Magnetized Bondi Accretion

Figure 2:  Top: Accretion rate (black) and magnetic flux parameter (blue) as a
function of time for our MHD simulation. The mean accretion rate
$\overline{\dot{M}}_{10}\sim 5\times 10^{-3}\dot{M}_{B}$ and the mean magnetic
flux at the horizon $\overline{\phi_{b}}\sim 30$ are plotted in gray and blue
horizontal lines respectively. We average values over the second half of the
simulations (pink background), when the simulation has reached steady state.
Bottom: Slices of $\beta$ (left) and $\rho$ (right) in our simulation,
spanning 8 orders of magnitude in spatial scale. Lines of constant magnetic
flux representing $\varphi$-averaged magnetic fields are overlaid in black.

Accretion flows in galactic nuclei involve magnetized plasma, and it has been
long known that magnetic fields strongly perturb the energetics (Shvartsman,
1971; Meszaros, 1975) and dynamics (Igumenshchev & Narayan, 2002; Pen et al.,
2003; Lalakos et al., 2022) of Bondi accretion. Apart from a toy model with
purely radial magnetic field and a non-spinning BH, where no dynamical
interaction occurs between the fluid and the magnetic field (Gammie et al.,
2003), there is no known analytical solution to the magnetized accretion
problem. It can be studied only via numerical simulations. Our computational
technique is well-suited for this.

To simulate the magnetohydrodynamic (MHD) problem, we use the same external
density and temperature as adopted for the HD case (Section 3), but now we
include an initial magnetic field in the $z$-direction with
plasma-$\beta\equiv p_{g}/p_{b}\approx 1$ where $p_{b}$ is the magnetic
pressure. Our choice of $\beta$ is motivated by the Milky Way ISM (e.g.,
Ferrière, 2020; Guerra et al., 2023), but we have verified that the results
remain unchanged if we set the initial $\beta\approx 10$ (Appendix C). We use
$n=8$ zones, which allows us to probe out to $r\sim 10^{8}r_{g}$ (and reach
steady state out to nearly $10^{7}r_{g}$). Our effective resolution is
$576\times 128\times 128$ over the whole domain. Additional set-up details are
given in Appendix B.2.

Unlike the laminar flow that we find in the HD simulation, the MHD simulation
gives a dynamic and turbulent flow which is driven by magnetic stresses and
reconnection. We need to time-average the results to obtain meaningful steady
state profiles. The black lines in Figure 1 show time-averaged radial profiles
of a number of quantities after the simulation was run for a physical duration
of $t_{\rm run}(R_{B})\sim 9t_{B}$ at the Bondi radius. The temperature
profile is fairly similar to the HD Bondi case (compare black and red lines),
but the density profile is significantly modified due to the presence of the
magnetic field. The radial scaling is now $\rho\propto r^{-1}$, which is
similar to the slope reported in several previous GRMHD simulations of
rotating accretion flows (e.g., Ressler et al., 2020; Chatterjee & Narayan,
2022) and even some hydrodynamical simulations (Guo et al., 2023). It is also
consistent with observational constraints in M87* and Sgr A* from X-ray and
EHT observations (e.g. Chatterjee & Narayan, 2022). There is a curious dip in
the density and a corresponding bump in the temperature near the Bondi radius,
which we explain in Section 5.

Despite running the simulation for an effective time of $\sim 9t_{B}$, we were
unable to obtain a reliable measurement of $\overline{\dot{M}}$ at $r\gtrsim
10^{2}\ r_{g}$ as seen in Figure 1(c), unlike the HD case. This is because the
medium exhibits violent convective turbulence (Section 5), which produces
large random velocities that overwhelm the mean accretion velocity. To
illustrate the difficulty, we show the average mass inflow rate $\dot{M}_{\rm
in}$ (solid blue) and mass outflow rate $\dot{M}_{\rm out}$ (dotted cyan),
defined in Appendix A, as a function of radius. With increasing distance from
the BH, both rates become much larger than the net $\overline{\dot{M}}$, which
is the difference of the two very similar numbers. It is therefore challenging
to calculate $\overline{\dot{M}}$. Here, we restrict ourselves to smaller
radii, where we believe our $\overline{\dot{M}}$ estimates are reliable. At
$r=10\,r_{g}$ 333It is customary in GRMHD simulations not to measure $\dot{M}$
at the horizon, where numerical floors might affect the fidelity of the
computed density, but to measure it at a larger radius. We determined that
$r=10r_{g}$ is a safe choice for our work. we find
$\overline{\dot{M}}_{10}\approx 0.005\dot{M}_{B}$, i.e., the accretion rate is
suppressed by more than a factor of 100 relative to the HD Bondi case. This
estimate is consistent with the observed accretion rate of $\sim
10^{-3}-10^{-2}\dot{M}_{B}$ in M87* (when comparing the BH accretion rate
$\dot{M}\sim 10^{-3}\,M_{\odot}$/yr in Event Horizon Telescope Collaboration
et al. (2021) with the Bondi accretion rate $\dot{M}_{B}\sim
0.2\,M_{\odot}$/yr in Russell et al. (2015)). It is also consistent with the
observed accretion rate for Sgr A* $10^{-4}-10^{-2}\,\dot{M}_{B}$ (Wang et
al., 2013).

We plot the time-averaged energy efficiency $\eta^{\rm tot}(r)$ (Equation 1)
from the simulation in Figure 1(d). There are two striking features. First,
the radial profile of $\eta^{\rm tot}$ is nearly constant (ideally it should
be perfectly constant), which indicates that the simulation energetics have
converged well.444The simulations indicate a factor $\sim 2$ increase in
$\eta^{\rm tot}$ close to the BH. We believe this is an artifact, caused by
the need to apply numerical floors near the BH horizon. Second, $\eta^{\rm
tot}$ is positive, i.e., energy flows out from the accretion flow to the
external medium. In other words, the simulation exhibits bona fide feedback.
It is worth emphasizing that we are able to study feedback only because our
numerical technique is specifically designed to allow two-way communication
between the BH and the external medium (Section 2). The amount of feedback
energy is $\sim 10^{-2}\overline{\dot{M}}_{10}c^{2}\sim 5\times
10^{-5}\,\dot{M}_{B}c^{2}$, which is lower than the values adopted in
cosmological simulations $\gtrsim 10^{-3}\,\dot{M}_{B}c^{2}$ (e.g., Booth &
Schaye, 2009; Weinberger et al., 2018; Davé et al., 2019; Ni et al., 2022).
However, the feedback in our runs operates without BH spin and without net
rotation in the accreting gas555Because of turbulent fluctuations, the gas at
each radius tends to have a random amount of rotation, at the level of a few
percent of Kepler. But the rotation is too weak to have any dynamical effect..
We speculate therefore, that it represents the minimum feedback that any hot
quasi-spherical magnetized accretion flow will produce. The actual feedback
will likely be larger when more favorable conditions are present (e.g., a
spinning BH).

The radial profiles of two magnetic field related quantities are shown in the
bottom panels of Figure 1: the plasma $\beta$, and the normalized enclosed
magnetic flux $\phi_{b}$ (Tchekhovskoy et al., 2011),

$\phi_{b}(r)\equiv\sqrt{\pi\big{/}\overline{\dot{M}}_{10}}\iint|B^{r}|\sqrt{-g}\,d\theta\,d\varphi.$
(3)

We average the inverse of $\beta$ for stability, which is a commonly used
technique. The profile we obtain has a flat plateau with a nearly constant
value $\beta\sim 3$ between $r=10^{2}-10^{5}$. Over this extended volume, the
field is in equi-partition with the gas and the magnetic flux $\phi_{b}$
varies $\propto r$, indicating that the organized poloidal field varies as
$B_{\rm pol}\propto r^{-1}$. The pressure in the organized field is also in
equi-partition: $B_{\rm pol}^{2}\propto\rho T\propto r^{-2}$. We find that for
radii $r<10^{2}$, $\beta$ declines rapidly to well below unity and $\phi_{b}$
saturates at $\sim 30$. Our results are very reminiscent of the magnetically
arrested disk (MAD) model of MHD accretion (Tchekhovskoy et al. 2011;
Igumenshchev et al. 2003; Narayan et al. 2003, see also Bisnovatyi-Kogan &
Ruzmaikin 1974). In previous work, Ressler et al. (2020) found that the entire
stellar wind-fed accretion flow in their model of Sgr A* approached the MAD
state. In our case, MAD accretion extends from the Bondi radius
$R_{B}\approx{1.8\times 10^{5}}\,r_{g}$ down to the BH horizon at $2\,r_{g}$.
One minor difference is that $\phi_{b}\sim 30$ at the BH horizon in our
simulation, whereas typical GRMHD simulations which consider orbiting gas
(instead of spherical infall) find $\phi_{b}\sim 50$ for a non-spinning BH
(Tchekhovskoy et al., 2012; Narayan et al., 2022).

##  5 Feedback via Reconnection-driven Convection

Figure 3:  Figure 4: Top: Dissecting contributions to the total energy outflow
($\eta^{\rm tot}$), including the fluid contribution ($\eta^{\rm fl}$), the
advective outflow efficiency ($\eta^{\rm adv}$), and the convective outflow
efficiency ($\eta^{\rm conv}$). Since $\eta^{\rm fl}\approx\eta^{\rm tot}$,
the electromagnetic field does not contribute much directly to the energy
outflow. Instead, the energy is transported primarily via convection until
near the Bondi radius (gray vertical dashed line). Positive (negative) values
are shown in solid (dotted) lines. Bottom: Azimuthally-averaged energy fluxes
in the $r,\theta$ plane, normalized by the accretion rate at $10\;r_{g}$.
Outflowing energy is shown in red and inflowing energy in blue.

As noted, our simulation permits us for the first time to study the coupled BH
feeding and feedback problem. Following the time evolution of the BH mass
accretion rate $\dot{M}$ (in units of the Bondi rate $\dot{M}_{B}$) and the
magnetic flux parameter $\phi_{b}$ at the horizon (top panel of Figure 2)
reveals large amplitude fluctuations in both quantities around a fairly stable
mean value, $\overline{\dot{M}}_{\rm 10}\sim 0.005\dot{M}_{B}$ and
$\overline{\phi}_{b}\sim 30$. The time scale of the largest fluctuations is a
few Bondi times, as expected for feeding from the Bondi radius. The cause of
the fluctuations is violent turbulence, evidence for which can be seen in the
multi-scale snapshots in the middle and lower rows of Figure 2.

To investigate the mechanism driving the turbulence, we explore the physics of
the outward flow of energy. The black curves in the upper panels of Figure 4
show the radial profile of the efficiency parameter $\eta^{\rm tot}$. The
other curves show individual physical components: $\eta^{\rm fl}$ (blue, left
top panel), the energy carried in the fluid (as opposed to the magnetic
field), $\eta^{\rm adv}$ (green, middle), the advected part of the fluid
energy, and $\eta^{\rm conv}$ (red, right), the part of the fluid energy
transport that can be attributed to convection. The definitions of these
components are given in Appendix A and are identical to those used by Ressler
et al. (2021).

At all radii $r\lesssim R_{B}$, we find that $\eta^{\rm fl}$ very closely
matches $\eta^{\rm tot}$, i.e., essentially all the energy flux is carried by
the fluid, reminiscent of the convection-dominated accretion model of Narayan
et al. 2000; Quataert & Gruzinov 2000), while the electromagnetic contribution
is negligible. The electromagnetic energy transport would likely become
important if the BH had spin, since field lines near the horizon would be
frame dragged and produce a relativistic jet. In a smaller scale simulation
with $R_{B}=GM_{\bullet}/c_{s,\infty}^{2}=100$ and a spinning BH
($a_{*}=0.9375$), Ressler et al. (2021) found that the electromagnetic $\eta$
exceeded $\eta^{\rm fl}$. For the non-spinning case here, we find the
opposite.

When we divide $\eta^{\rm fl}$ into an advective part $\eta^{\rm adv}$ and a
convective part $\eta^{\rm conv}$, we find that convection dominates by a
substantial margin between $r=r_{g}$ and $\sim 0.1R_{B}$; energy transport is
principally via convection. From the work of Igumenshchev & Narayan (2002),
the convection is initiated by magnetic reconnection near the BH, which
transforms the dissipated magnetic energy into thermal energy (entropy) and
kinetic energy (Ripperda et al., 2022). Convective turbulence develops and
extends all the way to $r\sim 0.1R_{B}$, carrying the dissipated energy
outward. Note that our definition of $\eta^{\rm conv}$ closely follows Ressler
et al. (2021). However, it is possible that some fraction of what we call
convection corresponds to a violently varying wind (e.g., Yuan et al., 2015).
Interestingly, although magnetic reconnection drives the convection, the field
does not transport the energy because field lines remain largely fixed in
radius in the MAD state.

The bottom three panels in Figure 4 present detailed information on the fluid,
advective and convective energy fluxes, $F_{E}^{\rm fl}$, $F_{E}^{\rm adv}$,
$F_{E}^{\rm conv}$ (defined in Appendix A), in the logarithmic radius - polar
angle ($\log_{10}{r}$ \- $\theta$) plane. The advective flux carries energy
toward the BH near the midplane but away from the BH near the poles. When
integrated over $\theta$, the two zones nearly cancel. On the other hand,
convection transfers energy outward at all $\theta$ and overwhelms the
advective energy at all $r\lesssim 0.1R_{B}$ so that the net fluid flux
$F_{E}^{\rm fl}$ is outward everywhere. While Ressler et al. (2021) concluded
that convection was inefficient in their spinning BH set-up, we find exactly
the opposite for our non-spinning BH. Pen et al. (2003) found energy inflow
rather than outflow in their non-relativistic simulation, which they termed
“frustrated convection.” The difference from our result might be because their
magnetic field did not reach the saturated MAD limit and/or their spatial
dynamic range was only a few tens.

Finally, we see in Figure 4 that the convective energy transport ceases at
$R_{B}$, as expected since convection requires gravity and buoyancy, while
gravity becomes sub-dominant for $r>R_{B}$. Thus, in our model, the feedback
energy carried outward by convection is dumped near the Bondi radius and has
nowhere to go. This explains the bump in the temperature profile of the MHD
run near the Bondi radius in Figure 1, and the corresponding dip in density
(required by pressure balance).

##  6 Summary and Conclusions

We have simulated general relativistic accretion on a non-spinning SMBH using
a numerical technique which allows us to maintain two-way communication across
the entire range of spatial scales from the BH horizon to beyond the Bondi
radius, $R_{B}\approx{1.8\times 10^{5}}r_{g}$. For spherical accretion without
magnetic fields, we find perfect agreement with the GR analytical Bondi
solution, i.e., for $r<R_{B}$, the density and temperature of the simulated
gas scale as $\rho\propto r^{-3/2}$, $T\propto r^{-1}$ respectively, as
expected, and the measured accretion rate also matches the analytical Bondi
accretion rate $\dot{M}_{B}$. The thermal energy of the external gas, which is
$1.5/R_{B}\approx 10^{-5}$ (Equation 2) of rest-mass energy, is advected into
the BH and there is no feedback in the opposite direction.

Adding a magnetic field of realistic strength to the external gas changes the
results drastically. The density slope inside $R_{B}$ changes to $\rho\propto
r^{-1}$, and the mass accretion rate $\dot{M}$ into the BH drops dramatically
to $\approx 5\times 10^{-3}\dot{M}_{B}$, consistent with observations of M87*
and Sgr A*. Equi-partition of magnetic and gas pressure is achieved
($\beta\sim 3$) over the radius range $10^{2}r_{g}\lesssim r\lesssim R_{B}$,
and the magnetic pressure becomes much more dominant close to the BH. Except
for the lack of net rotation, the solution closely resembles the magnetically
arrested disk (MAD) model.

Significantly, the magnetized simulation shows net energy feedback from the
accretion flow to the external medium at a rate $\approx 10^{-2}\dot{M}c^{2}$,
i.e., with an efficiency $\eta^{\rm tot}\sim 0.01$. Expressed in terms of the
Bondi accretion rate, the energy outflow rate is $\sim 5\times
10^{-5}\,\dot{M}_{B}c^{2}$. This positive feedback occurs even in the absence
of BH spin or net gas rotation. The energy outflow is not strongly collimated;
it is quasi-isotropic. Advection carries energy inwards in the mid-plane and
transports energy outwards near the poles, while convection transfers energy
outwards at all $\theta$. We find that the level of energy feedback is
significantly lower than that typically invoked in cosmological simulations,
so we view the feedback estimate in the present work as representing a lower
limit to energy outflow. Additional ingredients, especially BH spin and gas
angular momentum, are needed to study the problem more fully.

The dominant outward energy transport mechanism is convection driven by
magnetic reconnection in the magnetosphere near the BH. Even though a strong
magnetic field is essential for this mechanism to work, the energy is carried
almost entirely in the form of thermal and kinetic energy of the gas, while
the electromagnetic energy flux is negligible. Thus the feedback mechanism
seen here is different from the electromagnetically-dominated outflows found
in relativistic jets and rotation-driven disk winds.

These new insights are possible because our numerical technique permits two-
way communication between the BH and the external medium over many decades of
spatial and temporal scales. In future work, we plan to extend our
investigation to the case of a spinning BH that also includes the angular
momentum of the gas.

## Acknowledgements

The authors HC, RN, PN, KS, AR, and KC acknowledge support from the Gordon and
Betty Moore Foundation and the John Templeton Foundation via grants to the
Black Hole Initiative at Harvard University. This work used Delta at the
University of Illinois at Urbana-Champaign through allocation PHY230079 from
the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support
(ACCESS) program, which is supported by National Science Foundation grants
#2138259, #2138286, #2138307, #2137603, and #2138296. This research used
resources provided to BP by the Los Alamos National Laboratory Institutional
Computing Program, which is supported by the U.S. Department of Energy
National Nuclear Security Administration under Contract No. 89233218CNA000001.
We thank the referee for a thorough reading of the manuscript and constructive
comments that have helped improve it.

##  Appendix A GRMHD Primer and Definitions

Grid-based GRMHD codes such as KHARMA used in this work evolve a magnetized
fluid in a stationary spacetime (often a BH) described by a spacetime metric
$g_{\mu\nu}$. We focus on a non-spinning BH, and hence use the Schwarzschild
metric. The magnetized fluid is treated within the ideal MHD approximation
with the stress-energy tensor $T^{\mu}_{\nu}$ written as (see Anile 1989;
Komissarov 1999; Gammie et al. 2003 for details)

$T^{\mu}_{\nu}=\left(\rho+\gamma_{\rm
ad}u+b^{2}\right)u^{\mu}u_{\nu}+\left[p_{g}+(b^{2}/2)\right]\delta^{\mu}_{\nu}-b^{\mu}b_{\nu},$
(A1)

where $u^{\mu}$ is the four-velocity of the gas, $\rho$ is its rest mass
density, $u$ is its internal energy density, and $b^{\mu}$ is a four-vector
which describes the magnetic field in the fluid frame (Anile, 1989;
Komissarov, 1999) with $b^{\mu}u_{\mu}=0$ and $b^{\mu}b_{\mu}=b^{2}$. The gas
pressure $p_{g}=(\gamma_{\rm ad}-1)u$, where $\gamma_{\rm ad}=5/3$ is the
adiabatic index. The temperature $T$ is defined in relativistic units such
that $p_{g}=\rho T$, i.e., $T=(\gamma_{\rm ad}-1)u/\rho$. Given the initial
conditions of a system, the GRMHD code evolves the equation of mass
conservation, $(\rho u^{\mu})_{;\mu}=0$, and the equations of energy-momentum
conservation, $T_{\ \nu;\mu}^{\mu}=0$, along with the ideal MHD induction
equation. The output data describe the time evolution of the primitive
quantities, $\rho$, $u$, $u^{\mu}$, $b^{\mu}$, on each grid cell. At a given
instant of time and at a given location in the grid, the quantity $\rho
u^{r}\sqrt{-g}\,d\theta d\varphi$ measures the rate at which rest mass flows
out radially within a solid angle $d\theta d\varphi$, where $g={\rm
det}(g_{\mu\nu})$. The instantaneous net mass accretion rate through a sphere
of radius $r$ is then:

$\dot{M}(r)\equiv-\iint\rho u^{r}\sqrt{-g}\,d\theta\,d\varphi,$ (A2)

where the negative sign is because we associate accretion with inflow, which
corresponds to a negative value of $u^{r}$. Because the accreting gas has
large fluctuations in MHD (see the top panel in Figure 2), we are rarely
interested in the instantaneous $\dot{M}$. Rather, we average over an extended
period of time (e.g., the latter half of the entire simulation), and consider
the time-averaged mass accretion rate, which we write as
$\overline{\dot{M}}(r)$ with a bar. Furthermore, when we quote a unique
accretion rate, we use the value at $r=10r_{g}$, which we write as
$\overline{\dot{M}}_{10}$. Sometimes it is useful to distinguish between mass
that is flowing in at a given radius versus that flowing out (e.g., Fig. 1).
In this case, we separately define:

$\dot{M}_{\rm in}(r)\equiv\overline{-\iint\rho\ {\rm
Min}(u^{r},0)\sqrt{-g}\,d\theta\,d\varphi},\qquad\dot{M}_{\rm
out}(r)\equiv\overline{-\iint\rho\ {\rm
Max}(u^{r},0)\sqrt{-g}\,d\theta\,d\varphi}.$ (A3)

Per this definition, $\dot{M}_{\rm in}\geq 0$ (inflow) and $\dot{M}_{\rm
out}\leq 0$ (outflow). Following the convention used in this paper,
$\dot{M}_{\rm in}$ and $\dot{M}_{\rm out}$ are shown in Figure 1 as solid
lines (positive values) and dotted lines (negative values), respectively.

For the analyses in the paper, we define time-averaged and $\rho$-weighted
shell average of the quantity $X$ as:

$\langle X\rangle=\overline{\left(\frac{\iint
X\rho\sqrt{-g}~{}d\theta\,d\varphi}{\iint\rho\sqrt{-g}~{}d\theta\,d\varphi}\right)}.$
(A4)

The one exception is $\langle\rho\rangle$, which we compute as a straight
shell-average of the density without a further $\rho$-weighting. In the shell
averages, we leave out one closest cell from each of the two poles because
these are affected by the polar boundary conditions.

For computing the flow of energy, we note that $T_{t}^{r}\sqrt{-g}\,d\theta
d\varphi$ measures the rate at which energy flows in radially within a solid
angle $d\theta d\varphi$. Thus the instantaneous energy inflow rate through a
sphere at radius $r$ is given by

$\dot{E}(r)\equiv\iint T^{r}_{t}\sqrt{-g}\,d\theta\,d\varphi\,.$ (A5)

When energy flows toward the BH, $\dot{E}$ is positive. For studying feedback
we are interested in the flow of energy away from the BH. We are also not
interested in the flow of rest mass energy, which has no effect on the
thermodynamics. We thus define the time-averaged net rate of outflow of
“useful” energy as:

$\dot{E}_{\rm
net}(r)\equiv\overline{(\dot{M}-\dot{E})}=-\iint\overline{(T^{r}_{t}+\rho
u^{r})}\sqrt{-g}\,d\theta\,d\varphi\,,$ (A6)

and the corresponding efficiency $\eta^{\rm tot}$ as:

$\eta^{\rm tot}(r)=\dot{E}_{\rm
net}(r)/\overline{\dot{M}}_{10}=\overline{(\dot{M}-\dot{E})}/\overline{\dot{M}}_{10}\,.$
(A7)

To study the contribution of the fluid alone (ignoring the magnetic field),
the relevant terms in the stress-energy tensor are:

$\left(T^{\rm fl}\right)^{r}_{t}\equiv(\rho+\gamma_{\rm ad}u)u^{r}u_{t}.$ (A8)

Then the time-averaged net fluid energy outflow rate is:

$\dot{E}^{\rm fl}_{\rm net}(r)\equiv-\iint\overline{\left(T^{\rm fl}+\rho
u^{r}\right)^{r}_{t}}\sqrt{-g}\,d\theta\,d\varphi=\iint\overline{Be~{}\rho
u^{r}}\sqrt{-g}\,d\theta\,d\varphi\,,$ (A9)

where the relativistic Bernoulli parameter is (e.g., Penna et al., 2013, but
not including the magnetic field)

$Be=-\left(1+\gamma_{\rm ad}u/\rho\right)u_{t}-1.$ (A10)

Following Ressler et al. (2021), we further sub-divide the fluid energy
outflow rate into an advective part,

$\dot{E}^{\rm adv}_{\rm net}(r)\equiv\iint\overline{Be}~{}\overline{\rho
u^{r}}\sqrt{-g}\,d\theta\,d\varphi\,,$ (A11)

and a convective part,

$\dot{E}^{\rm conv}_{\rm net}(r)\equiv\dot{E}^{\rm fl}_{\rm
net}(r)-\dot{E}^{\rm adv}_{\rm net}(r)\,.$ (A12)

In the same manner as in (A7), we define the corresponding efficiencies,
$\eta^{\rm fl}(r)$, $\eta^{\rm adv}(r)$, $\eta^{\rm conv}(r)$.

In addition to the time-averaged and shell-integrated radial profiles
$\dot{E}(r)$ defined above, we also define local time- and azimuth-averaged
energy fluxes as a function of $r$ and $\theta$:

$\displaystyle F_{E}^{\rm fl}(r,\theta)$ $\displaystyle\equiv$
$\displaystyle(1/2\pi)\int\overline{Be~{}\rho u^{r}}\sqrt{-g}\,d\varphi,$
(A13) $\displaystyle F_{E}^{\rm adv}(r,\theta)$ $\displaystyle\equiv$
$\displaystyle(1/2\pi)\int\overline{Be}~{}\overline{\rho
u^{r}}\sqrt{-g}\,d\varphi,$ (A14) $\displaystyle F_{E}^{\rm conv}(r,\theta)$
$\displaystyle\equiv$ $\displaystyle F_{E}^{\rm fl}(r,\theta)-F_{E}^{\rm
adv}(r,\theta).$ (A15)

##  Appendix B Numerical Set-up

###  B.1 Hydrodynamic Simulation Setup

In the hydrodynamic simulation described in Section 3, we use a domain with 7
zones ($i=0-6$) covering the range $[1,1.7\times 10^{7}]\;r_{g}$. Zone-6 is
initialized with a static fluid of constant density and temperature, and no
perturbations, so as to recover the analytic Bondi solution without rotation.
We use spherical Kerr-Schild coordinates with coordinate four-vector
$(x^{t},x^{r},x^{\theta},x^{\varphi})$, where the radial coordinate is
$x^{r}\equiv\log{r}$, sometimes referred to as “eKS” coordinates. Rather than
allow vacuum to arise in the inner regions before the initial material from
zone-6 accretes, new areas are filled with material before they are first
simulated, with the same density and temperature as the innermost nonzero
density zone present when they are initialized. Initialization in this way
prevents disruption of the simulation, but bears no resemblance to the Bondi
solution, validating that the scheme can converge to the correct result
regardless of its initialization. After the first pass, no further
initialization is needed. During the rest of the simulation, each active
annulus receives half its initial data from the previous active annulus and
half from the last output of the current active annulus. In the process, the
inner and outer radial boundaries of the active annulus end up with new values
of $\rho$, $u$ and $u^{\mu}$, which are treated as fixed boundary conditions
while the annulus is evolved. Once the system has reached its final steady
state, conservation of mass and energy between annuli are satisfied in a time-
averaged sense.

###  B.2 Magnetohydrodynamic Simulation Setup

The magnetized simulation in Section 4 is run over a larger domain with 8
zones $[1,1.3\times 10^{8}]\;r_{g}$. Zone-7 is initialized with a similar
density and temperature as in the hydrodynamical model, but the two differ in
several ways. The density is now initialized in the same way in all zones, to
follow $\rho_{\rm init}(r)=\rho_{0}(r+R_{B})/r$ such that it is $\propto
r^{-1}$ interior to $R_{B}$ and constant outside of $R_{B}$. This is motivated
by simulations (e.g. Ressler et al., 2020; Chatterjee & Narayan, 2022; Guo et
al., 2023) and observations (Chatterjee & Narayan, 2022) that find a
$\rho\propto r^{-1}$ relation in the presence of strong magnetic fields.
However, the results are insensitive to the choice of $\rho_{\rm init}$, as
discussed in Appendix C. The temperature is initialized to the Bondi
analytical solution, and the velocities are initialized to zero. The magnetic
fields are initialized with a vector potential
$A_{\varphi}(r,\theta)=b_{z}r(r+R_{B})\sin^{2}{\theta}/2$ such that the
plasma-$\beta$ parameter is approximately constant over all radii. The
resulting magnetic field is vertical outside $R_{B}$, and becomes more radial
inside $R_{B}$. The field is normalized with the parameter $b_{z}$ such that
the initial $\beta_{\rm init}(r)\sim 1$, motivated by observations of the
Milky Way ISM (Ferrière, 2020; Guerra et al., 2023). A per-zone random
perturbation of maximum 5% is applied to the internal energy when initializing
fluid at all radii, to help break axisymmetry. Since the magnetic field picks
out a rotation axis for the resulting flow parallel to the coordinate axis, we
choose a coordinate system which compresses the $\theta$ coordinate away from
coordinate axes and toward the midplane. We adopt a simplified version of
“funky modified Kerr-Schild” (FMKS) coordinates (see e.g. $\theta_{j}$ in the
Appendix in Wong et al., 2021), eliminating the cylindrification term to
obtain:

$\theta=Ny\left(1+\frac{(y/\chi_{t})^{\alpha}}{\alpha+1}\right)+\frac{\pi}{2},$
(B1)

where $N\equiv\pi/2(1+\chi_{t}^{-\alpha}/(1+\alpha))^{-1}$ is a normalization
factor, $y=2x^{\theta}-1$, $\chi_{t}=0.8$, $\alpha=16$, and
$x^{\theta}\in[0,1]$. The simplified system retains the primary benefit of
FMKS coordinates: zones near the coordinate pole are widened, relaxing the
Courant condition and allowing a greater simulation timestep, by a factor of 2
or more, dramatically reducing simulation cost. As we are interested primarily
in steady state behavior rather than resolving particular structures, e.g., a
stable jet, the cylindrification is unnecessary for our case.

In addition to changing the coordinates, we lower the Courant factor for
safety. We also introduce dynamical floors: in order to maintain
$\rho>10^{-6}\,r^{-3/2}$, internal energy $u>10^{-8}\,r^{-5/2}$, temperature
$u/\rho<100$, and magnetization $b^{2}/\rho<100$, we introduce material and
energy in the coordinate frame. We also reduce the fluid velocity when the
Lorentz factor measured by the Eulerian observer is larger than $\gamma_{\rm
max}=10$. While running each active annulus, we keep $\rho$, $u$ and $u^{\mu}$
fixed at the inner and outer boundaries, as in the hydro case. In addition, we
also keep $b^{\mu}$ fixed and make sure that we preserve $\nabla\cdot B=0$.

We have validated the ‘multi-zone’ method described in Section 2 in the
presence of strong magnetic fields, by comparing a smaller 4-zone simulation
with $R_{B}\approx 460$ against a single-zone simulation of the same problem.
We obtained substantially similar averaged radial profiles. These comparisons,
along with detailed descriptions of the method, boundary conditions, etc.,
will be presented in a forthcoming longer paper.

##  Appendix C Resolution and initial condition study

Figure 5: Comparison between simulations with varying resolutions, coordinate
systems, and initial magnetic strength, described in Table 1. Each profile is
averaged between $t_{\rm run}(R_{B})=4.5\,t_{B}-9\,t_{B}$ measured at the
Bondi radius (gray vertical dashed line). Positive (negative) values are shown
in solid (dotted) lines. Overall, we find good agreement between these
simulations.

We have presented our fiducial MHD simulation results in Section 4 for a
resolution of $128^{3}$ per annulus with the FMKS coordinate system and
initialized with $\rho_{\rm init}\propto r^{-1}$ at $r<R_{B}$. Here we show
that the results remain the same even when we choose different resolutions,
coordinate systems, or initial conditions. First, we introduce a new “Modified
Kerr-Schild” (MKS) coordinate system (see e.g. $\theta_{g}$ in the Appendix of
Wong et al., 2021) where the resolution is more focused in the midplane but
the zones near the pole are not widened:

$\theta=\pi x^{\theta}+\frac{1}{2}(1-h)\sin{(2\pi x^{\theta})}$ (C1)

Here we use $h=0.3$. In the case of $h=1$, it is equivalent to the uniform
$\theta$ grid of the eKS coordinates in Section B.1.

Six extra simulations were run with parameters outlined in Table 1. Using MKS
coordinates and other set-up details the same as in our fiducial run, we ran
simulations at (i) $48^{3}$ resolution per annulus and (ii) $64^{3}$
resolution per annulus. Run (iii) with $96^{3}$ resolution per annulus used
FMKS but with a different $\alpha=14$ which is less extreme than $\alpha=16$
in the fiducial run. Run (iv) used a weaker initial magnetic field with
$\beta\sim 10$ and used MKS coordinates with $64^{3}$ resolution per annulus.
Runs (v) and (vi) were initialized with a different initial density profile
compared to the other runs and used MKS with $64^{3}$ resolution. Run (v) was
initialized with a piecewise constant density profile as in the HD run. Run
(vi) was initialized with the analytical HD Bondi solution where the density
scales as $\rho_{\rm init}\propto r^{-3/2}$ for $r<R_{B}$.

Radial profiles of all the quantities of interest for the six extra
simulations are compared with the fiducial run in Figure 5. Overall, the
profiles are very similar between the various runs implying that our
conclusions from the fiducial run hold true regardless of the choice of
coordinate system, resolution, initial magnetic field strength, or initial
density profile. The accretion rate at $10\,r_{g}$ of all five runs lies in
the range of $\overline{\dot{M}}\approx(1-5)\times 10^{-3}\,\dot{M}_{B}$,
which is consistent with observations of M87* (Event Horizon Telescope
Collaboration et al., 2021; Russell et al., 2015) and Sgr A* (Wang et al.,
2013). The outflow efficiency is around $\eta^{\rm tot}\approx 1-3\,\%$, where
the highest value is from the lowest resolution ($48^{3}$) run (i). Even
though initialized with a weaker magnetic field, run (iv) still shows magnetic
field saturation at $\beta\sim 3$ at radii $r=10^{2}-10^{4}$. Near the
horizon, run (iv) accreted the least amount of magnetic field as can be
inferred from its larger $\beta$ and smaller $\phi_{b}$. However, the
difference is not large. The $\beta$-parameter decreases more steeply towards
the BH for run (iii) and the fiducial run, both of which use FMKS coordinates.
We suspect this is because the FMKS coordinates are characterized by larger
cells near the poles where $\beta^{-1}$ tends to be larger, and this might
have affected the shell average. Runs (v) and (vi), despite their very
different initial density profiles, converge to a similar final density
profile with $\rho\propto r^{-1}$. This demonstrates that our results are not
dependent on the initial conditions.

Label | Resolution | Coordinates | Initial $\beta$ | Initial $\rho$ $(r<R_{B})$
---|---|---|---|---
$128^{3}$, fmks (fiducial) | $128^{3}$ | FMKS ($\chi_{t}=0.8$, $\alpha=16$) | $\sim 1$ | $r^{-1}$
(i) $48^{3}$ | $48^{3}$ | MKS ($h=0.3$) | $\sim 1$ | $r^{-1}$
(ii) $64^{3}$ | $64^{3}$ | MKS ($h=0.3$) | $\sim 1$ | $r^{-1}$
(iii) $96^{3}$, fmks | $96^{3}$ | FMKS ($\chi_{t}=0.8$, $\alpha=14$) | $\sim 1$ | $r^{-1}$
(iv) $64^{3}$, $\beta 10$ | $64^{3}$ | MKS ($h=0.3$) | $\sim 10$ | $r^{-1}$
(v) $64^{3}$, const | $64^{3}$ | MKS ($h=0.3$) | $\sim 1$ | piecewise constant
(vi) $64^{3}$, bondi | $64^{3}$ | MKS ($h=0.3$) | $\sim 1$ | $r^{-3/2}$
Table 1: Simulation set-up for different runs.

## References

  * Anglés-Alcázar et al. (2021) Anglés-Alcázar, D., Quataert, E., Hopkins, P. F., et al. 2021, ApJ, 917, 53, doi: 10.3847/1538-4357/ac09e8
  * Anile (1989) Anile, A. M. 1989, Relativistic fluids and magneto-fluids: With applications in astrophysics and plasma physics
  * Bisnovatyi-Kogan & Ruzmaikin (1974) Bisnovatyi-Kogan, G. S., & Ruzmaikin, A. A. 1974, Ap&SS, 28, 45, doi: 10.1007/BF00642237
  * Bondi (1952) Bondi, H. 1952, MNRAS, 112, 195, doi: 10.1093/mnras/112.2.195
  * Booth & Schaye (2009) Booth, C. M., & Schaye, J. 2009, MNRAS, 398, 53, doi: 10.1111/j.1365-2966.2009.15043.x
  * Chatterjee & Narayan (2022) Chatterjee, K., & Narayan, R. 2022, ApJ, 941, 30, doi: 10.3847/1538-4357/ac9d97
  * Chatterjee et al. (2023) Chatterjee, K., Chael, A., Tiede, P., et al. 2023, Galaxies, 11, 38, doi: 10.3390/galaxies11020038
  * Davé et al. (2019) Davé, R., Anglés-Alcázar, D., Narayanan, D., et al. 2019, MNRAS, 486, 2827, doi: 10.1093/mnras/stz937
  * Event Horizon Telescope Collaboration et al. (2021) Event Horizon Telescope Collaboration, Akiyama, K., Algaba, J. C., et al. 2021, ApJ, 910, L13, doi: 10.3847/2041-8213/abe4de
  * Ferrarese & Merritt (2000) Ferrarese, L., & Merritt, D. 2000, ApJ, 539, L9, doi: 10.1086/312838
  * Ferrière (2020) Ferrière, K. 2020, Plasma Physics and Controlled Fusion, 62, 014014, doi: 10.1088/1361-6587/ab49eb
  * Fiacconi et al. (2018) Fiacconi, D., Sijacki, D., & Pringle, J. E. 2018, MNRAS, 477, 3807, doi: 10.1093/mnras/sty893
  * Gammie et al. (2003) Gammie, C. F., McKinney, J. C., & Tóth, G. 2003, ApJ, 589, 444, doi: 10.1086/374594
  * Gebhardt et al. (2000) Gebhardt, K., Bender, R., Bower, G., et al. 2000, ApJ, 539, L13, doi: 10.1086/312840
  * Grete et al. (2023) Grete, P., Dolence, J. C., Miller, J. M., et al. 2023, 37, 465, doi: 10.1177/10943420221143775
  * Guerra et al. (2023) Guerra, J. A., Lopez-Rodriguez, E., Chuss, D. T., Butterfield, N. O., & Schmelz, J. T. 2023, AJ, 166, 37, doi: 10.3847/1538-3881/acdacd
  * Guo et al. (2023) Guo, M., Stone, J. M., Kim, C.-G., & Quataert, E. 2023, ApJ, 946, 26, doi: 10.3847/1538-4357/acb81e
  * Hopkins & Quataert (2010) Hopkins, P. F., & Quataert, E. 2010, MNRAS, 407, 1529, doi: 10.1111/j.1365-2966.2010.17064.x
  * Hopkins et al. (2023) Hopkins, P. F., Grudic, M. Y., Su, K.-Y., et al. 2023, arXiv e-prints, arXiv:2309.13115.  https://arxiv.org/abs/2309.13115
  * Igumenshchev & Narayan (2002) Igumenshchev, I. V., & Narayan, R. 2002, ApJ, 566, 137, doi: 10.1086/338077
  * Igumenshchev et al. (2003) Igumenshchev, I. V., Narayan, R., & Abramowicz, M. A. 2003, ApJ, 592, 1042, doi: 10.1086/375769
  * Kaaz et al. (2023) Kaaz, N., Murguia-Berthier, A., Chatterjee, K., Liska, M. T. P., & Tchekhovskoy, A. 2023, ApJ, 950, 31, doi: 10.3847/1538-4357/acc7a1
  * Komissarov (1999) Komissarov, S. S. 1999, MNRAS, 303, 343, doi: 10.1046/j.1365-8711.1999.02244.x
  * Kormendy & Ho (2013) Kormendy, J., & Ho, L. C. 2013, ARA&A, 51, 511, doi: 10.1146/annurev-astro-082708-101811
  * Lalakos et al. (2022) Lalakos, A., Gottlieb, O., Kaaz, N., et al. 2022, ApJ, 936, L5, doi: 10.3847/2041-8213/ac7bed
  * Li & Bryan (2014) Li, Y., & Bryan, G. L. 2014, ApJ, 789, 153, doi: 10.1088/0004-637X/789/2/153
  * Magorrian et al. (1998) Magorrian, J., Tremaine, S., Richstone, D., et al. 1998, AJ, 115, 2285, doi: 10.1086/300353
  * Meszaros (1975) Meszaros, P. 1975, A&A, 44, 59
  * Michel (1972) Michel, F. C. 1972, Ap&SS, 15, 153, doi: 10.1007/BF00649949
  * Narayan et al. (2022) Narayan, R., Chael, A., Chatterjee, K., Ricarte, A., & Curd, B. 2022, MNRAS, 511, 3795, doi: 10.1093/mnras/stac285
  * Narayan et al. (2000) Narayan, R., Igumenshchev, I. V., & Abramowicz, M. A. 2000, ApJ, 539, 798, doi: 10.1086/309268
  * Narayan et al. (2003) —. 2003, PASJ, 55, L69, doi: 10.1093/pasj/55.6.L69
  * Ni et al. (2022) Ni, Y., Di Matteo, T., Bird, S., et al. 2022, MNRAS, 513, 670, doi: 10.1093/mnras/stac351
  * Pen et al. (2003) Pen, U.-L., Matzner, C. D., & Wong, S. 2003, ApJ, 596, L207, doi: 10.1086/379339
  * Penna et al. (2013) Penna, R. F., Kulkarni, A., & Narayan, R. 2013, A&A, 559, A116, doi: 10.1051/0004-6361/201219666
  * Porth et al. (2019) Porth, O., Chatterjee, K., Narayan, R., et al. 2019, ApJS, 243, 26, doi: 10.3847/1538-4365/ab29fd
  * Prather et al. (2021) Prather, B. S., Wong, G. N., Dhruv, V., et al. 2021, 6, 3336, doi: 10.21105/joss.03336
  * Quataert & Gruzinov (2000) Quataert, E., & Gruzinov, A. 2000, ApJ, 539, 809, doi: 10.1086/309267
  * Ressler et al. (2021) Ressler, S. M., Quataert, E., White, C. J., & Blaes, O. 2021, MNRAS, 504, 6076, doi: 10.1093/mnras/stab311
  * Ressler et al. (2020) Ressler, S. M., White, C. J., Quataert, E., & Stone, J. M. 2020, ApJ, 896, L6, doi: 10.3847/2041-8213/ab9532
  * Ricarte et al. (2019) Ricarte, A., Tremmel, M., Natarajan, P., & Quinn, T. 2019, MNRAS, 489, 802, doi: 10.1093/mnras/stz2161
  * Ripperda et al. (2022) Ripperda, B., Liska, M., Chatterjee, K., et al. 2022, ApJ, 924, L32, doi: 10.3847/2041-8213/ac46a1
  * Rosas-Guevara et al. (2016) Rosas-Guevara, Y., Bower, R. G., Schaye, J., et al. 2016, MNRAS, 462, 190, doi: 10.1093/mnras/stw1679
  * Russell et al. (2015) Russell, H. R., Fabian, A. C., McNamara, B. R., & Broderick, A. E. 2015, MNRAS, 451, 588, doi: 10.1093/mnras/stv954
  * Shapiro & Teukolsky (1983) Shapiro, S. L., & Teukolsky, S. A. 1983, Black holes, white dwarfs, and neutron stars : the physics of compact objects
  * Shvartsman (1971) Shvartsman, V. F. 1971, Soviet Ast., 15, 377
  * Sijacki et al. (2015) Sijacki, D., Vogelsberger, M., Genel, S., et al. 2015, MNRAS, 452, 575, doi: 10.1093/mnras/stv1340
  * Talbot et al. (2021) Talbot, R. Y., Bourne, M. A., & Sijacki, D. 2021, MNRAS, 504, 3619, doi: 10.1093/mnras/stab804
  * Tchekhovskoy et al. (2012) Tchekhovskoy, A., McKinney, J. C., & Narayan, R. 2012, in Journal of Physics Conference Series, Vol. 372, Journal of Physics Conference Series, 012040, doi: 10.1088/1742-6596/372/1/012040
  * Tchekhovskoy et al. (2011) Tchekhovskoy, A., Narayan, R., & McKinney, J. C. 2011, MNRAS, 418, L79, doi: 10.1111/j.1745-3933.2011.01147.x
  * Trott et al. (2022) Trott, C. R., Lebrun-Grandié, D., Arndt, D., et al. 2022, 33, 805, doi: 10.1109/TPDS.2021.3097283
  * Wang et al. (2013) Wang, Q. D., Nowak, M. A., Markoff, S. B., et al. 2013, Science, 341, 981, doi: 10.1126/science.1240755
  * Weinberger et al. (2018) Weinberger, R., Springel, V., Pakmor, R., et al. 2018, MNRAS, 479, 4056, doi: 10.1093/mnras/sty1733
  * Weinberger et al. (2023) Weinberger, R., Su, K.-Y., Ehlert, K., et al. 2023, MNRAS, 523, 1104, doi: 10.1093/mnras/stad1396
  * Wellons et al. (2023) Wellons, S., Faucher-Giguère, C.-A., Hopkins, P. F., et al. 2023, MNRAS, 520, 5394, doi: 10.1093/mnras/stad511
  * Wong et al. (2021) Wong, G. N., Du, Y., Prather, B. S., & Gammie, C. F. 2021, ApJ, 914, 55, doi: 10.3847/1538-4357/abf8b8
  * Yuan et al. (2015) Yuan, F., Gan, Z., Narayan, R., et al. 2015, ApJ, 804, 101, doi: 10.1088/0004-637X/804/2/101


# The sub-millimetre non-uniformity measurement of residual and coil-generated
field in the magnetic shield using atomic vapor cell

Chen Liu  Haifeng Dong∗<EMAIL_ADDRESS>Junjun Sang

###### Abstract

Magnetic field source localization and imaging happen at different scales. The
sensing baseline ranges from meter scale such as magnetic anomaly detection,
centimeter scale such as brain field imaging to nanometer scale such as the
imaging of magnetic skyrmion and single cell. Here we show how atomic vapor
cell can be used to realize a baseline of 109.6 $\mu$m with a magnetic
sensitivity of<EMAIL_ADDRESS>and a dynamic range of 2062-4124nT. We use
free induction decay (FID) scheme to suppress low-frequency noise and avoid
scale factor variation for different domains due to light non-uniformity. The
measurement domains are scanned by digital micro-mirror device (DMD). The
currents of 22mA, 30mA, 38mA and 44mA are applied in the coils to generate
different fields along the pumping axis which are measured respectively by
fitting the FID signals of the probe light. The residual fields of every
domain are obtained from the intercept of linearly-fitting of the measurement
data corresponding to these four currents. The coil-generated fields are
calculated by deducting the residual fields from the total fields. The results
demonstrate that the hole of shield affects both the residual and the coil-
generated field distribution. The potential impact of field distribution
measurement with an outstanding comprehensive properties of spatial
resolution, sensitivity and dynamic range is far-reaching. It could lead to
capability of 3D magnetography for small stuffs and/or organs in millimeter or
even smaller scale.

††preprint: AIP/123-QED

Measuring the distribution of spin polarization and magnetic field in atomic
vapor cell has been pursued for long time. Optical magnetic resonance imaging
(OMRI) is used to detect the spin polarization distribution since
1990s1997Three ; 1997Optical ; 1997OpticalRb ; 1998Alkali . As a magnetic
gradient field for space coding is used in OMRI, the measured spin
polarization distribution is usually generated by pumping light profile
instead of field non-uniformity1999Optical ; 2000Visualization ; 2000Diffusion
; article ; Savukov2015Gradient . In 2003, photodetector array is firstly used
to measure the field distribution in atomic vapor cell with a baseline of 3
mmkominis2003subfemtotesla . Since then, this technology has been verified and
applied in magnetoencephalography by many groups with a baseline of 3 mm, 4
mm, 5 mm and 7.5 mm, respectively2009Three ; 2018High ;
Cort2010Magnetoencephalography ; 2017A . Charge coupled device (CCD) is also
used to record fluorescence emitted from Cs atoms and detect the field
distribution in the vapor cell caused by magnetic nanoparticle sample2016A ;
2017Characterizing . In 2017, Shuji Taue et.al. use DMD to measure the radio
frequency (RF) field distribution in atomic vapor cell2020Signal . While the
baseline is sub-millimetre ($\sim$260$\mu$m), the field sensitivity and drift
is unknown. Our group measured the spin spatial frequency response of atomic
vapor cell using DMD in 20192019Dong . Furthermore, we demonstrate that spin
image with a linewidth of 13.7 $\mu$m can be achieved in atomic vapor cell,
which is much smaller than the corresponding diffusion crosstalk free
distance2019Spin . DMD is also used in measuring the coil-generated field
distribution under SERF regime, with a sensitivity of $\sim$20fT@30-40Hz and a
dynamic range of $\sim$8nTfang2020high .

Here we tried to obtain residual and coil-generated field distribution in the
magnetic shield simultaneously with high spatial resolution. We measure non-
uniformity of these fields using thermal atomic vapor cell heated to 85 ∘C.
Although the effect of shield on the coil coefficient has been realized and
measured recentlyMa2021 , we first demonstrate that the distribution of coil-
generated field in the shield has strong correlation with the distribution of
the residual field due to the hole of the shield.

Figure 1: Experimental setup for magnetic field distribution measurement. (P:
polarizer; AOM: Acoustic Optical Modulator; $\lambda/4$: quarter-wave plate;
$\lambda/2$: half-wave plate; L: lens; PBS: polarization beam splitter; DAQ:
Data Acquisition )

The experiment setup is depicted in Fig. 1. We use orthogonal pump-probe
arrangement. Cs atoms are contained in a $\rm 2.5cm\times 2.5cm\times 2.5cm$
cubic cell. The cell is filled with 600 Torr of $\rm{}^{4}He$ buffer gas and
150 Torr quenching gas $\rm N_{2}$. It is placed in an oven heated by non-
magnetic twisted wires. We use a set of three nested cylindrical magnetic
shields with a shielding factor larger than $10^{4}$. A set of saddle coil is
used to generate radio frequency magnetic field along $y$ direction. Another
set of circular coil is used to generate the measured field along $x$
direction. Cs atoms are optically

(a)

pumped by a 180mW distributed feedback (DFB) diode laser of 894nm
(PH895DBR240TS, Photodigm). Another DFB laser of 894nm is used to probe the
spin polarization in $z$ direction. The measurement scheme is shown in Fig.
LABEL:ModulatedPumpingSignal(a). DMD (Vialux V-7001) is used to scan the spin
polarization of different domains. It is comprised of 1024$\times$768 digital
micromirrors, each of which is a square reflectors with a side length of 13.7
$\mu$m. We set 8$\times$8 micromirrors as one scanning domain, which
corresponds to a baseline of 109.6$\mu$m. Balanced polarimetry is used to
measure the rotation of the plane of polarizaiton of probe light corresponding
to every selected domain. The differential signal is recorded by data
acquisition card NI4461.

To suppress low-frequency noise and avoid scale factor variation for different
domains due to light non-uniformity, we use free induction decay to measure
the field. Single beam FID magnetometry has been researched in Ref. 2018Free,
. There are also double-beam FID schemes where amplitude-modulated pumping
light is usually used to polarize the spins and the probe light is used to
detect the FID signal2011An ; 2021All . To enhance the signal amplitude, we
use RF stimulated FID scheme. Fig. LABEL:ModulatedPumpingSignal(b)-(e) shows
the time sequence diagram of the measurement. Cs atoms are pumped from 0 to
$t_{0}$ and the pumping light is shut-off then by an acoustic optic modulator
(TEM-80-10-894, BRIMROSE). After that RF signal with a frequency of
$\omega_{0}=\gamma B_{0}$ is turned on from $t_{0}$ to $t_{1}$ so that the
spin polarization turns from $x$ to $y$, as show in Fig.
LABEL:ModulatedPumpingSignal(a). After that FID spin polarization signal is
probed by Faraday rotation of the linearly-polarized probe light. Both the
duration $\Delta t=t_{1}-t_{0}$ and amplitude of RF are optimized to obtain
the maximum FID amplitude.

The evolution of spin polarization is described by the Bloch
equation1946Nuclear

$\frac{\mathrm{d}}{\mathrm{d}t}\mathbf{P}=\gamma\mathbf{P}\times\mathbf{B}+R_{p}(s\hat{z}-\mathbf{P})-R_{rel}\mathbf{P}$
(1)

where $\gamma$ is the gyromagnetic ratio, $R_{p}$ is the pumping rate,
$R_{rel}$ is the relaxation rate.

From 0 to $t_{0}$, $B_{x}=B_{0}$, $B_{y}=0$, $B_{z}=0$ and $R_{p}=R_{p0}$. The
spins are polarized along $x$ direction and the amplitude at $t_{0}$ is

$P_{x}(t_{0})=\dfrac{R_{p0}}{R_{p0}+R_{rel}}$ (2)

From $t_{0}$ to $t_{1}$, $R_{p}=0$, $B_{x}=B_{0}$, $B_{y}=B_{\rm
rf}\cos(\omega t)$ and $B_{z}=0$. The RF magnetic field duration $\Delta
t=t_{1}-t_{0}=\dfrac{\pi}{2\omega_{\rm rf}}$, where $\omega_{\rm rf}$ is lamor
precession frequency of Cs under $B_{\rm rf}$. The spin polarization is
rotated from $x$ direction to $z$ direction by the RF field. We set $B_{\rm
rf}=1250{\rm nT}$, which corresponds to a spin polarization loss of less than
10$\%$. The polarization along $z$ direction at the end of $t_{1}$ is

$P_{z}(t_{1})=\dfrac{R_{p0}}{R_{p0}+R_{rel}}{\rm e}^{-\dfrac{\pi}{2\gamma
B_{\rm rf}}R_{rel}}$ (3)

After that, RF magnetic field is shut off at $t_{1}$, then the spin
polarization precess and decay in the $yOz$ plane according to Eq. (4).

$\begin{cases}P_{y}(t)=P_{z}(t_{1}){\rm e}^{-R_{rel}t}\sin(\gamma B_{0}t)\\\
P_{z}(t)=P_{z}(t_{1}){\rm e}^{-R_{rel}t}\cos(\gamma B_{0}t)\end{cases}$ (4)

Figure 3: (a) 0.025s FID signal. (b) Magnetic field data obtained by fitting
FID signals. (c) PSD of the data in (b).

As the signal from balance polarimetry is proportional to $P_{z}$, we use
model Eq. (4) to fit the recorded free induction decay signal.

FID signal of 30s for one domain at current of 44mA is recorded and the data
from 0 to 0.025s is shown in Fig. 3(a). After fitting, we obtain about 10710
field measurement values, as shown in Fig. 3(b). Fig. 3(c) is the power
spectral density (PSD) which shows a magnetic sensitivity of
<EMAIL_ADDRESS>In the scanning measurement, we record five seconds data
for every domain under currents of 22, 30, 38 and 44 mA, respectively. The
corresponding fields are from 2062 to 4124nT. The total fields of every domain
corresponding to these four currents are obtained by FID signal fitting and
averaging the five second data.

Figure 4: The total fields and the coil-generated fields in different domains
corresponding to currents of 22mA and 44mA. The coil-generated fields are
obtained by deducting the residual fields from the total fields.

By linearly fitting the measured field corresponding to 22, 30, 38 and 44mA at
every domain, we obtain the residual field distribution from the intercepts.
After deducting the residual fields from total fields for every domains, the
coil-generated field distribution is also obtained. To verify whether the
residual field is deducted effectively from the total field, we list all the
domain’s total fields and coil-generated fields corresponding to 22mA and 44mA
from small to large In Fig. 4. The black triangle represents total fields and
the red square represents coil-generated fields which are equal to total
fields minus residual fields. The blue dash line in Fig. 4 shows the fitting
result of $B_{44}=2.00026B_{22}$ with $\rm R^{2}$=0.9732. The reason why the
scale factor is not equal to 2 is that the current accuracy is 0.01mA
(Thorlabs LDC 205C). It is clearly seen from Fig. 4 that the residual field is
deducted effectively.

Figure 5: The residual magnetic field distribution in the middle of the
shield.

The image of the residual field distribution is shown in Fig. 5 and the image
of coil-generated field distribution corresponding to 44mA is shown in Fig. 6.
The solid black points in Fig. 5 and Fig. 6 indicate the middle positions of
the measurement domains.

Comparing Fig. 5 with Fig. 6, we can see that the residual field distribution
has strong correlation with coil-generated field distribution. The correlation
coefficient is -0.98416. The minus sign means that the hole of the shield has
opposite effect on the residual field and coil-generated field distribution.
This phenomenon has never been observed. However, it is reasonable considering
that the shield can affect the coefficient of the magnetic coil Ma2021 . So
for any experiment and device where the distribution of coil-generated field
is important, one should consider the effect of shield structure beyond the
Biot-Savart law simulation.

Figure 6: Coil-generated field distribution at a current of 44mA in the middle
of the shield.

In summary, we have demonstrated a magnetic field distribution measurement
scheme with a baseline of 109.6 $\mu$m, sensitivity of<EMAIL_ADDRESS>and a dynamic range of 2062-4124nT. The residual field distribution and coil-
generated field distribution in middle of the magnetic shield are obtained
simultaneously. We found that there are strong correlation between the
residual and coil-generated field distribution at sub-millimeter scale. One of
the possible reasons is that the hole of the shield affects both of these two
field distributions. With further improvement of the sensitive domain
deployment, this scheme can be used in magnetic field source localization and
imaging of millimeter scale objects.

This work is supported by the National Natural Science Foundation of China
(61973021).

Data Availability The data that support the findings of this study are
available from the corresponding author upon reasonable request.

## References

  * (1) A. R. Young, S. Appelt, B. A. Baranga, C. Erickson, and W. Happer,  Applied Physics Letters 70, 3081 (1997).
  * (2) J. Skalla, G. Wackerle, and M. Mehring,  Optics Communications 143, 209 (1997).
  * (3) J. Skalla, G. Wackerle, M. Mehring, and A. Pines,  Physics Letters A 226, 69 (1997).
  * (4) B. A. Baranga, S. Appelt, C. J. Erickson, A. R. Young, and W. Happer,  Physical Review A 58, 2282 (1998).
  * (5) K. Ishikawa, Y. Anraku, Y. Takahashi, and T. Yabuzaki,  Journal of the Optical Society of America B-optical Physics 16, 31 (1999).
  * (6) K. Ishikawa, S. Morii, K. Enomoto, and T. Yabuzaki,  Josa B 17, 182 (2000).
  * (7) D. Giel, G. Hinz, D. Nettels, and A. Weis,  Opt. Express 6, 251 (2000).
  * (8) A. Weis, D. Giel, and D. Nettels,  Laser Physics 11, 470 (2001).
  * (9) Savukov and I.,  Journal of Magnetic Resonance 256, 9 (2015).
  * (10) I. Kominis, T. Kornack, J. Allred, and M. V. Romalis,  Nature 422, 596 (2003).
  * (11) A. Gusarov, D. Levron, E. Paperno, R. Shuker, and B. A. Baranga,  IEEE Transactions on Magnetics 45, 4478 (2009).
  * (12) K. Nishi, Y. Ito, and T. Kobayashi,  Optics Express 26, 1988 (2018).
  * (13) P. D. D. S. Cort Johnson and M. Weisend,  Applied Physics Letters 97, 413 (2010).
  * (14) Y. Mamishin, Y. Ito, and T. Kobayashi,  IEEE Transactions on Magnetics 53, 1 (2017).
  * (15) V. Dolgovskiy et al.,  Applied Physics Letters 109 (2016).
  * (16) A. Weis et al.,  Journal of Physics Conference 793, 012032 (2017).
  * (17) S. Taue and Y. Toyota,  Optical Review 27 (2020).
  * (18) H.-F. Dong et al.,  Journal of Applied Physics 125, 023908 (2019).
  * (19) H. F. Dong et al.,  Journal of Applied Physics 125, 243904.1 (2019).
  * (20) X. Fang et al.,  Optics Express 28, 26447 (2020).
  * (21) D. Ma et al.,  Journal of Physics D: Applied Physics 54, 275001 (2021).
  * (22) D. Hunter et al.,  Physical Review Applied 10, 014002 (2018).
  * (23) A. Gusarov, D. Levron, B. A. Baranga, E. Paperno, and R. Shuker,  Journal of Applied Physics 109, 130801 (2011).
  * (24) Applied Physics B 127, 1 (2021).
  * (25) F. Bloch, W. W. Hansen, and M. Packard,  Physical Review 69, 127 (1946).