# Photon rings of spherically symmetric black holes and robust tests of non-
Kerr metrics

Maciek Wielgus<EMAIL_ADDRESS>Black Hole Initiative at Harvard
University, 20 Garden Street, Cambridge, MA 02138, USA  Center for
Astrophysics $|$ Harvard & Smithsonian, 60 Garden Street, Cambridge, MA 02138,
USA

###### Abstract

Under very general assumptions on the accretion flow geometry, images of a
black hole illuminated by electromagnetic radiation display a sequence of
photon rings (demagnified and rotated copies of the direct image) which
asymptotically approach a purely theoretical critical curve – the outline of
the black hole photon shell. To a distant observer, these images appear
dominated by the direct emission, which forms a ring whose diameter is
primarily determined by the effective radius of the emitting region. For that
reason, connecting the image diameter seen by a distant observer to the
properties of the underlying spacetime crucially relies on a calibration that
necessarily depends on the assumed astrophysical source model. On the other
hand, the diameter of the photon rings depends more on the detailed geometry
of the spacetime than on the source structure. As such, a photon ring
detection would allow for the spacetime metric to be probed in a less model-
dependent way, enabling more robust tests of General Relativity (GR) and the
Kerr hypothesis. Here we present the photon ring structure of several
spherically symmetric black hole spacetimes and perform comparisons with the
Schwarzschild/Kerr case. We offer our perspective on future tests of the
spacetime metric with photon rings, discussing challenges and opportunities
involved.

††preprint: APS/123-QED

##  I Introduction

The observational appearance of a black hole (BH) has been studied extensively
for many years from a purely theoretical standpoint [1, 2, 3, 4]. This
appearance depends on numerous astrophysical details, such as geometry and
optical depth of the emission region [5, 6, 7]. The feature depending
exclusively on the spacetime geometry is the critical curve [8], an outline of
the photon shell on a distant observer’s screen. However, since the null
geodesics comprising the photon shell are unstable orbits [9], the critical
curve remains a purely mathematical concept. Unless the accretion flow itself
is spherically symmetric [4, 6] (which is typically not an astrophysically
plausible scenario), the image of a BH contains a sequence of discrete photon
rings [8, 10]. Each photon ring is a copy of the direct image of the emission
region, corresponding to photon trajectories lensed around the photon shell,
with the $n$-th photon ring photons executing $n/2$ loops before reaching the
observer’s screen [5, 10], see Fig. 1. This sequence of photon rings quickly
converges to the critical curve, and the diameter and shape of the $n\\!=\\!2$
photon ring is already strongly dominated by the spacetime geometry – hence it
can be used for precise tests of the theory of gravity, with only very general
assumptions on the astrophysical model of the accretion flow [10, 11]. In the
case of non-Kerr spacetimes, most of the previous research has focused on the
theoretical critical curve (commonly referred to as shadow, occasionally also
as photon ring, which has a different meaning in this paper), which has been
characterized for a broad range of non-Kerr spacetimes, e.g., [12, 13, 14, 15,
16]. Relatively less work has been dedicated to the observable photon ring
structure. Additionally, the observable ring-like direct emission feature is
commonly conflated with the critical curve in the literature. One of the aims
of this paper is to clarify that issue.

Figure 1: A simple model of the emission around a spherically symmetric static
black hole. The observer is located to the right, at $r\rightarrow\infty$,
viewing the accretion disk face-on. Black hole has a horizon radius $r_{\rm
h}$, there is a photon sphere located at $r_{\gamma}$ (blue circle). The
emission from the accretion disk (represented with orange ellipses) is
dominated by the flux emitted at the effective radius $r_{\rm e}$. The photons
emitted at $r_{\rm e}$ may reach the observer along infinite number of
trajectories, corresponding to adding half loops around the photon sphere.
First 3 such trajectories are shown. They correspond to the direct image
($n=0$), and first and second photon ring ($n=1$ and $n=2$). Radii of $n$-th
photon rings in the observer’s plane converge rapidly to the critical impact
parameter, $\ell_{n}\rightarrow\ell_{\rm c}$, which defines the critical curve
on the observer’s screen.

As noticed by [10], the technique of very long baseline interferometry (VLBI)
is uniquely suited to probe the photon ring structure of black hole images. As
sharp image domain features, photon rings die-off slowly in the Fourier
domain, dominating the interferometric signal in the range of very high
spatial frequencies, corresponding to baseline lengths of tens and hundreds
billion wavelengths. While such long baselines were never observed so far and
this feat will surely require significant progress and new developments in the
field of space VLBI, the projects that would enable robust detection of the
photon ring structure are already being discussed and proposed, e.g., [11, 10,
17, 18, 7, 19], and algorithms to enable photon ring detection in VLBI data
are already being implemented and tested with synthetic data [11, 20, 21].

Only recently first resolved observational images of a BH were obtained by the
Event Horizon Telescope VLBI array (EHT; [22, 23, 24, 25, 26, 27]), presenting
the supermassive BH M 87∗ in the center of the M 87 galaxy with a resolution
corresponding to $\sim 3\,r_{\rm S}$ (Schwarzschild radii; [24, 25, 26]). The
measured diameter of a bright ring feature found in the EHT images was used to
test the null hypothesis of consistency with the Kerr metric [28, 27] under
the BH mass prior derived from stellar dynamics [29]. To formulate this test,
the scaling relationship between the ring diameter and the BH mass assuming
Kerr metric was evaluated using 103 synthetic images generated with general
relativistic magnetohydrodynamic simulations (GRMHD; [26, 27]). These
synthetic images, representing variety of viable theoretical GRMHD models of
BH accretion, allowed to perform diameter-mass calibration and quantify the
null hypothesis test. The results of this test were later on reinterpreted as
EHT constraints on the diameter of the critical curve to be 17 % consistent
with that of the Kerr spacetime [30, 31]. However, the exact value of the
uncertainty depends on the selection of simulated images and underlying
parameters of the models used to perform calibration, and since the images are
computed under Kerr spacetime assumption, the uncertainties are not directly
transferable to tests of non-Kerr spacetimes. Thus, the 17 % uncertainties
assumed by [30, 31] can only be considered a proxy to demonstrate the
principle of testing spacetime geometry with the EHT observations, but one
should be extremely careful with the quantitative interpretation of any
constraints provided in this way. While numerical GRMHD constitutes the most
sophisticated modeling framework we have at our disposal to interpret the EHT
observations, there are good reasons to question whether GRMHD simulations can
represent the appearance of an accreting BH accurately enough to be used for a
reliable diameter-mass calibration (the fact that certain GRMHD models are
consistent with the EHT observations only means fulfilling one necessary
condition). Unlike objections related to the applications to non-Kerr
spacetimes, these concerns challenge the conclusions obtained within Kerr
paradigm as well. Some of these concerns were enumerated by [32] and we will
not repeat them here. Instead, we focus on demonstrating how the measurement
of the direct emission-dominated ring image diameter translates (under a prior
on BH mass) to a measurement of the effective emission region size, largely
independent of the underlying metric, while photon rings exhibit much more
sensitivity to the underlying spacetime geometry. To that end we investigate
the photon ring structure of several spherically symmetric BH spacetimes
(spacetimes with an event horizon) and compare them with Schwarzschild/Kerr
metrics.

##  II Astrophysical constraints

VLBI observations have a potential to resolve supermassive BHs in the centers
of galaxies. Two targets of largest angular diameter are M 87∗ [22] and
Sagittarius A∗ (Sgr A∗), a supermassive black hole in the center of our own
galaxy [33]. A number of other supermassive BH sources could be resolved in a
foreseeable future with space VLBI [17]. Apart from the concerns about the
angular resolution, a potential source must emit sufficiently strong radiation
at relevant frequencies ($10^{11}-10^{12}$ Hz). Spectral observations
unambiguously indicate [e.g., 34, 35] that the accretion flow in the objects
of interest corresponds to a hot and radiatively inefficient accretion flow
(RIAF; also commonly referred to as the advection dominated accretion flow -
ADAF) model [36, 37, 38, 39, 40], where the mass accretion rate is very low,
and the radiation observed by VLBI instruments originates from a synchrotron
process. A classic thin accretion disk model [41] represents a very different
accretion regime (geometrically thinner, optically thicker, less hot, more
dense, emitting thermal radiation) and therefore _should not_ be used for
quantitative interpretations of the VLBI images of supermassive BHs. A
particularly relevant difference seems to be that in the analytic thin disk
model radiative flux terminates at the innermost stable circular orbit (ISCO),
which is not the case for the RIAF models – thus the VLBI images of BHs are
not expected to necessarily retain clear information about the ISCO location.
While we expect the RIAF synchrotron emission at all radii above the event
horizon (located at radius $r_{\rm h}$), in the horizon limit photons are
infinitely redshifted. Hence, for the axisymmetric flow, regardless of the
details of the accretion model, there must be an effective radius of emission
$r_{\rm e}>r_{\rm h}$, where majority of photons reaching the distant observer
originate from.

In the simple model discussed in this paper we technically assume that all of
the emission originates at particular $r_{\rm e}$, only computing images of an
infinitely thin ring of radius $r_{\rm e}$ and disregarding the spatially
extended character of the emission in a more physically realistic scenario. A
definition of $r_{\rm e}$ may be more ambiguous ambiguous in a realistic case,
and depending on the details of the emission profile. However, that is exactly
the reason why observing photon rings, far less dependent on $r_{\rm e}$ than
the direct image, is a tantalizing perspective.

A separate concern is related to the optical depth – the observations need to
penetrate the outer layers of the accretion flow and allow to peer into its
innermost region in the immediate vicinity of the event horizon. Synchrotron
emission has two regimes, self-absorbed optically thick regime at lower
observing frequencies, and optically thin regime for higher frequencies with a
steep decrease of flux as a function of frequency [42]. These regimes are
separated by the emission maximum at a turn-over frequency $\nu_{t}\propto
M_{\rm BH}^{-1}\dot{M}^{2/3}$ [43] for black hole mass $M_{\rm BH}$ and mass
accretion rate $\dot{M}$. Hence, the observing frequency needs to be close to
the turn-over frequency $\nu_{t}$ to reach the optically thin regime, but
still collect sufficient amount of flux to enable VLBI detections.
Fortunately, that seems to be the case for both M 87∗ and Sgr A∗ at the
current EHT observing frequency of $2.3\times 10^{11}$ GHz [33, 22].

Figure 2: Transfer functions comparison between Schwarzschild (thick lines)
and maximally charged Reissner-Nordström (thin lines) spacetimes. Assuming
known $M/D$, EHT observations of M 87∗ interpreted as a ring dominated by the
$n\\!=\\!0$ emission component allow to calculate the effective emission
radius $r_{\rm e}$ (horizontal magenta line). The EHT ring diamater can be
explained by $r_{\rm e}=4.61\,M$ in Schwarzschild or by $r_{\rm e}=4.70\,M$ in
RN spacetime. Unlike direct emission $n=0$ transfer function, higher order
transfer functions are separated along the observed radius horizontal axis,
and hence may be used as a robust discriminant between the two spacetimes.

The fundamental physics behind RIAF flow and synchrotron radiation does not
depend on the spacetime metric. However, the same can not be said about the
detailed dynamics of the accretion flow and relative importance of the
gravitational redshift, both influencing the location of the effective flux
maximum $r_{\rm e}$. The location of $r_{\rm e}$ is also strongly affected by
physics of plasma dissipation and electron heating, which remains highly
uncertain in the relevant regime. While these aspects of physics are not
directly dependent on the spacetime geometry, they remain insufficiently
understood and are prescribed rather than self-consistently calculated in the
GRMHD simulations [26]. Hence, the impact of the accretion model on the
$r_{\rm e}$ may easily be mistaken for the effect of the spacetime geometry.
This is illustrated in an example shown in Fig. 2. We consider the case of a
spherically symmetric spacetime. Hence, critical curve in the observer’s
screen corresponds to a circle with radius equal to a critical impact
parameter $\ell_{\rm c}$ (the exact scaling requires the knowledge of the mass
and distance of the observed BH), see Fig. 1. We also limit ourselves to the
case of an observer viewing the accretion disk around a static black hole
face-on (along the spin axis of the accretion disk). In that idealized case
photon rings are simply concentric circles in the distant observer’s screen.
Figure 2 presents transfer functions [5, 8] connecting the emission radius
$r_{\rm e}$ (in BH mass units) and the impact parameter $\ell_{n}$ of $n$-th
photon ring for Schwarzschild and Reissner-Nordström (RN) metric with charge
parameter $q\\!=\\!1$ (see Section III). [31] excluded RN metric with
$q\\!=\\!1$ based on the aforementioned EHT $17\,\%$ Kerr-GRMHD-tuned error
budget. Let us critically examine this conclusion. For simplicity, let us
assume that the mass and distance of M 87∗ are known perfectly well, with no
uncertainties, yielding ratio of $M/D=3.8\,\mu$as ([22]; accounting for the
additional uncertainties would only strengthen our argument). Then the EHT
observation of a 42 $\mu$as ring gives the impact parameter of $\ell_{\rm
EHT}=5.53M$, a thick vertical magenta line in Fig. 2. Since the observed
emission is strongly dominated by the $n\\!=\\!0$ direct component [5, 10], we
can use the $n\\!=\\!0$ transfer function to translate this measurement to the
effective radius of emission $r_{\rm e}$ equal to 4.61 $M$ for Schwarzschild
(thick horizontal magenta line), and 4.70 $M$ for RN spacetime (thin
horizontal magenta line, the details of the calculations are presented in
Section IV). These values are very close to one another, because transfer
functions of a direct image ($n=0$, red tilted lines) are very similar in both
spacetimes. In both cases they are well approximated by a simple linear
relation (tilted dashed black line; [8, 44])

$r_{\rm e}\approx\ell-M\ .$ (1)

Hence, the argument against RN spacetime with $q=1$ hinges on the assertion
that we understand the astrophysics of the accreting system well enough to
confidently claim that $r_{e}$ in RN spacetime must be of value lesser than
$4.70M$.

The position we take in this paper is that such claim is unsubstantiated,
given the serious uncertainties about the physics governing the radiative flux
profile (where is the effective maximum $r_{\rm e}$ is located) and lack of
dedicated accretion simulations performed in the RN spacetime (or in non-Kerr
BH spacetimes in general, with a notable exception of [45]). On the other
hand, notice that a detection of the $n\\!=\\!1$ (or, better, $n>1$) photon
ring would immediately allow to distinguish between these two spacetimes (at
least in the case of a known $M/D$), as there is no overlap between the
observable impact parameters of $n\\!=\\!1$ photon rings in the two spacetimes
(at least for $r_{\rm e}<20M$), irrespectively of the underlying astrophysics.

##  III Non-Kerr models of black holes

To demonstrate the properties of photon rings and the direct emission ring,
let us consider a general spherically symmetric spacetime with a metric in a
form

$\mathrm{d}s^{2}=g_{\mu\nu}\mathrm{d}x^{\mu}\mathrm{d}x^{\nu}=-f(r)\mathrm{d}t^{2}+\frac{g(r)}{f(r)}\mathrm{d}r^{2}+r^{2}\mathrm{d}\Omega^{2}$
(2)

The square root of $f(r)$ is commonly referred to as a lapse function and
$\mathrm{d}\Omega^{2}=\mathrm{d}\theta^{2}+\sin^{2}\theta\mathrm{d}\phi^{2}$
corresponds to a line element on a unit sphere. Given the spherical symmetry
we will restrict ourselves to $\theta=\pi/2$ plane, where
$\mathrm{d}\Omega=\mathrm{d}\phi$. All solutions are parametrized with the ADM
mass of the central object $M$. Our reference solution is that of a
Schwarzschild spacetime, a special case of Kerr spacetime in the static limit
(BH spin $a_{*}=a/M=0$), with

$f_{\rm S}(r)=1-\frac{2M}{r}\ \ ;\ \ g_{\rm S}(r)\equiv 1\ .$ (3)

Following the Birkhoff theorem, it constitutes a unique spherically symmetric
static vacuum solution to Einstein field equations (EFE) of GR, with a horizon
at the Schwarzschild radius $r_{\rm h}=r_{\rm S}=2M$. In this paper we discuss
the photon ring structure of several other spacetimes, with metrics
parametrized with a single charge parameter $q$, cast in the units of $M$,
reducing to a Schwarzshild solution for $q=0$. Our selection of models is by
no means complete and roughly corresponds to a subset of spacetimes considered
in [31, 46]. In this section we briefly introduce these spacetimes.

Figure 3: Geodesic structure of several spherically-symmetric BH spacetimes.
The radius $r_{0}$ corresponds to the emission radius $r_{\rm e}$ for which
the direct emission ring image has a radius corresponding to that of a
Schwarzschild critical curve ($3\sqrt{3}M$). Note that $r_{0}$ is only very
weakly dependent on the spacetime metric and is consistent with $4.3\pm 0.1M$
for all considered cases.

Reissner-Nordström metric (RN) is a classic solution to electrovacuum EFE,
representing a static electrically charged black hole in GR. While there exist
astrophysical arguments against accumulation of large electric charge $q$ by
BHs, the RN solution remains a useful framework to test the $\propto r^{-2}$
deviations in $f(r)$. Moreover, RN solution also arises in the Randall-Sundrum
model of the brane cosmology [47], where $q$ is reinterpreted as a tidal
charge, representing gravitational effects from the higher dimensions. RN
solution is defined with

$f_{\rm RN}(r)=1-\frac{2M}{r}+\frac{q^{2}}{r^{2}}\ \ ;\ \ g_{\rm RN}(r)\equiv
1\ .$ (4)

Diameter of the RN black hole critical curve was studied by [48].

Regular black hole solutions attempt to avoid the presence of the central
singularity, where the curvature is unbounded and GR loses its validity. First
such solution corresponded to a magnetically charged BH and was proposed by
Bardeen [49] with

$f_{\rm B}(r)=1-\frac{2Mr^{2}}{\left(r^{2}+q^{2}\right)^{3/2}}\ \ ;\ \ g_{\rm
B}(r)\equiv 1\ .$ (5)

This postulated metric inspired regular BH solutions of EFE coupled to
nonlinear electrodynamics, e.g., by Ayón-Beato and García [50],

$f_{\rm
AG}(r)=1-\frac{2Mr^{2}}{\left(r^{2}+q^{2}\right)^{3/2}}+\frac{q^{2}r^{2}}{\left(r^{2}+q^{2}\right)^{2}}\
\ ;\ \ g_{\rm AG}(r)\equiv 1,$ (6)

and by Hayward [51], with

$f_{\rm H}(r)=1-\frac{2Mr^{2}}{r^{3}+2q^{2}M}\ \ ;\ \ g_{\rm H}(r)\equiv 1\ .$
(7)

The critical curves of regular black holes were studied by [52] and [53].

Kazakov-Solodukhin black hole solution (KS) considers the deformation of the
Schwarzschild solution due to spherically symmetric quantum fluctuations,
yielding a metric with

$f_{\rm KS}(r)=-\frac{2M}{r}+\frac{1}{r}\left(r^{2}-q^{2}\right)^{1/2}\ \ ;\ \
g_{\rm KS}(r)\equiv 1.$ (8)

The critical curve diameter of the KS black hole was considered by [54].

Einstein-Maxwell-dilaton black hole (EMd) is an electromagnetically charged BH
solution in the Einstein-Maxwell-dilaton theory, where the dilaton field
couples to the electromagnetic field [55, 56]. Following [46] we consider a
following static EMd metric

$\displaystyle f_{\rm EMd}(r)=1-\frac{\sqrt{4M^{2}r^{2}+q^{4}}-q^{2}}{r^{2}}\
;$ $\displaystyle g_{\rm EMd}(r)=\frac{4M^{2}r^{2}}{4M^{2}r^{2}+q^{4}}.$ (9)

The critical curves of EMd BHs has been studied in [57].

Differences between these spacetimes are shown in Fig. 3, where we compare
their properties by plotting characteristic radii relevant for the geodesic
dynamics of particles [58] – event horizon radius $r_{\rm h}$, photon sphere
radius $r_{\gamma}$, marginally bound orbit radius $r_{\rm mb}$, and ISCO
radius $r_{\rm isco}$ as a function of the charge parameter $q$. The relevant
radii can be found as roots of the following equations

$\displaystyle r_{\rm h}\,:\,f=0,$ (10) $\displaystyle
r_{\rm\gamma}\,:\,f^{\prime}r-2f=0,$ (11) $\displaystyle r_{\rm
mb}\,:\,f^{\prime}r-2f(1-f)=0,$ (12) $\displaystyle r_{\rm
isco}\,:\,f^{\prime\prime}r-2r(f^{\prime})^{2}+3ff^{\prime}\\!=\\!0\ .$ (13)

Additionally, we plot a line $r_{0}$, representing the effective radius of
emission $r_{\rm e}$, that would result in a direct image of ring with a
radius of $\ell_{0}=3\sqrt{3}M$, which is the critical curve radius in
Schwarzshild spacetime. The charge parameter $q$ has a very limited influence
on the value of $r_{0}$, which remains in the $4.3\pm 0.1M$ range for all the
considered spacetimes. This is in contrast to some other characteristic radii,
for example between $q=0$ and $q=1$ ISCO location in RN spacetime changes by
over 30%, while $r_{0}$ changes by less than 4%. Moreover, Kerr spacetime
viewed face-on also gives consistent $r_{0}$ for all considered values of spin
($0\leq a/M<1)$, Fig. 4. This shows that in a hypothetical case of a direct
emission ring measurement perfectly consistent with the Schwarzschild critical
curve, one can not claim a metric test (constraint on $a$ or $q$) without
additional strong assumptions on the accretion model (what $r_{\rm e}$ is
allowed for which spacetimes). On the other hand, one could claim a constraint
on the object compactness, broadly independent on the underlying spacetime
metric. In case of the EHT observations of M 87∗, assuming mass and distance
priors [27], one would find the M 87∗ upper bound on compactness limited by
$r_{\rm e}/M=4.5\pm 0.7$, where uncertainties are dominated by the object mass
and distance uncertainties [29, 27].

Figure 4: Similar as Fig. 3 and Fig. 5, but for the the Kerr spacetime,
parametrized with a spin parameter $a/M$, with the accretion disk viewed from
the BH spin axis, aligned with the accretion disk spin axis. Characteristic
radii represent the equatorial slice of the Kerr spacetime [58].

Figure 5: Image domain structure of the selected spacetimes. Locations of
characteristic radii are denoted with prime. Dashed lines represent the
$\ell=r_{\rm e}+M$ formula. Blue bands indicate the range of possible radii
for the observed photon rings. For $n=1$ lighter band represents full range of
$r_{\rm h}<r_{\rm e}<\infty$, and darker band assumes $r_{\rm e}<10M$. The
narrow $n=2$ band covers whole range of $r_{\rm e}$ and contains the critical
curve.

##  IV Photon Ring Calculations

For the case of Kerr spacetime, the null geodesic structure has been studied
extensively, and geometry of the photon rings can be calculated analytically
without resorting to ray-tracing or explicit numerical integration of the
geodesic equation, see, e.g., [5, 8]. In a spherically symmetric case of a
metric given by Eq. 2 with two unspecified functions $f(r)$ and $g(r)$ we
resort to the numerical integration of the deflection angle $\delta$ to find
the location of the photon rings, noticing that the angular path between the
emission point and the distant observer is equal to

$\delta_{n}=\frac{\pi}{2}+n\pi\,$ (14)

in our simplified model of the system geometry (see Fig. 1). Specific angular
momentum $\ell=-u_{\phi}/u_{t}$, is conserved along geodesics due to the
Killing symmetries of the static axisymmetric spacetimes, and corresponds to
the impact parameter of the photon as seen on a screen of a distant observer,
Fig. 1. Hence, using the deflection angle formula, we can write down an
implicit relation between the effective emission radius $r_{\rm e}$ and the
impact parameter of the $n$-th photon ring $\ell_{n}$ for a given metric
$g_{\mu\nu}$

$\displaystyle F_{n}(\ell_{n};r_{\rm e},g_{\mu\nu})$
$\displaystyle=\Bigg{|}\int\displaylimits^{r_{\rm
e}}_{\infty}\frac{u^{\phi}}{u^{r}}\mathrm{d}r\Bigg{|}-\delta_{n}$
$\displaystyle=\Bigg{|}\int\displaylimits^{r_{\rm
e}}_{\infty}\frac{\mathrm{d}r}{\sqrt{\frac{r^{2}}{g(r)}\left[\frac{r^{2}}{\ell_{n}^{2}}-f(r)\right]}}\Bigg{|}-\delta_{n}\
.$ (15)

The $n$-th photon ring corresponds to $F_{n}(\ell_{n};r_{\rm e},g_{\mu\nu})=0$
if between the emission and the distant observer trajectory does not go
through a pericenter of the null geodesic. In other case, the pericenter
radius $r_{\rm p}$ can be found from the radial four-velocity $u^{r}=0$
condition, the total angular path along the fully extended photon trajectory
is then

$\Delta\phi_{\rm tot}=2\Bigg{|}\int\displaylimits^{r_{\rm
p}}_{\infty}\frac{u^{\phi}}{u^{r}}\mathrm{d}r\Bigg{|}$ (16)

and the impact parameter of the photon ring $\ell_{n}$ can be found by
comparing $F^{\prime}_{n}$ to zero

$F^{\prime}_{n}(\ell_{n};r_{\rm e},g_{\mu\nu})=\Delta\phi_{\rm
tot}-F_{n}(\ell_{n};r_{\rm e},g_{\mu\nu})-2\delta_{n}\ .$ (17)

That is the case for $n\\!=\\!1,2$ (but not for $n\\!=\\!0$) trajectories
shown in Fig. 1. In the figure we indicate the locations of the pericenters
$r_{\rm p}$. A practical criterion for using Eq. 17 rather than Eq. 15 is that
$r_{\rm p}$ exists and $\Delta\phi_{\rm tot}<2\delta_{n}$, so that the
pericenter is located between the emitter and the observer. The critical curve
location can be found without resorting to numerical integration as

$\ell_{\rm c}=\frac{r_{\gamma}}{\sqrt{f(r_{\gamma})}}\ ,$ (18)

where the photon sphere radius $r_{\gamma}$ is given as a solution to Eq. 11.
We illustrate the image domain characteristic radii for each spacetime in Fig.
5. Here, the location of direct images of characteristic radii from Fig. 3 are
shown (black continuous lines) and denoted with a prime. In particular,
$r_{\rm h}^{\prime}$ can be associated with the central brightness depression
of [40] or with the inner shadow of [59]. Dashed lines represent the
approximated location of direct images calculated with the Eq. 1. This trivial
formula is reasonably accurate for all considered spacetimes, yielding at most
$\sim 20\%$ deviation in case of $r_{\rm h}^{\prime}$ (that is $r_{\rm
e}\gtrapprox r_{\rm h}$) and a large charge, and much less in all other cases.
This once again shows that the direct image is not a robust and direct probe
of the metric, and instead it rather unambiguously probes the location of the
emission $r_{\rm e}/M$ independently of the details of the metric (note that
the accuracy of this measurement is ultimately limited by the uncertainties in
the distance $D$ and simplified geometry of the model). The figure also
indicates possible range of radii of the $n=1$ and $n=2$ photon rings as seen
by a distant observer. Similar plot for the Kerr spacetime parametrized with
spin $a/M$ is given in the bottom panel of the Fig. 4. We only consider the
case of the Kerr BH spin aligned along the the line of sight, and for the
calculations we follow procedures outlined in [8]. It is now evident that
while the direct image approximately reproduces the structure of the emission
(shifted by $1\,M$ following Eq. 1), photon rings transfer different emission
radii into small range of radii in the observer’s plane, see also transfer
functions example discussed in Section II. For all considered cases $n=1$
approximates the critical curve within 15%, while $n=2$ closer than 2.2%,
independently of the emission profile. Hence, photon rings are sensitive to
the spacetime geometry.

Figure 6: Lyapunov exponents $\gamma$ calculated for the considered families
of spacetimes with charges $q$. For Schwarzschild case we find $\gamma=\pi$.
For spacetimes with smaller $\gamma$, we expect relatively more flux in the
subsequent photon rings. For $\gamma=2.45$ the flux ratio difference with
respect to the Schwarzschild case corresponds to a factor of 2.  Figure 7:
Left: range of possible photon ring radii in the observer’s plane,
corresponding to $n=1$ photon ring, assuming the emission radius $r_{\rm
e}<10M$. Solid lines correspond to critical curves in different spacetimes.
Gray band represents $n=1$ photon ring radii for Kerr metric with spin $0\leq
a/M<1$. Middle: same as left, but assuming $3.8M<r_{\rm e}<5.2M$ following the
EHT constraints on the M 87∗. Right: same, but for the $n=2$ photon ring and
$r_{\rm h}<r_{\rm e}<\infty$.

Another interesting characteristic of a spacetime, related to its photon ring
structure, is provided by the Lyapunov exponents $\gamma$ [10], quantifying
the instability of the bound photon shell geodesics. In case of spherically
symmetric spacetimes, if photon’s radial position is perturbed from the
spherical geodesic $r_{\gamma}$ radially by an infinitesimal $\delta r_{0}$,
after $n$ half-orbits the photon will appear at the radius $\delta r_{n}$ away
from $r_{\gamma}$ and

$\delta r_{n}=e^{\gamma n}\delta r_{0}.$ (19)

We evaluate the Lyapunov exponents numerically as $\gamma\approx\ln(\delta
r_{1}/\delta r_{0})$ for the considered families of spacetimes, see Fig. 6.
For Schwarzschild case $\gamma=\pi$ [3]. Lyapunov exponents control the
relative flux in subsequent photon rings, and hence constitute an observable
quantity. While for Schwarzschild case flux in subsequent rings is
extinguished by a factor of $e^{-\pi}\approx 0.04$, for a spacetime with
$\gamma=2.45$ we obtain a twice less severe extinction – in such a spacetime
the photon rings are expected to be significantly brighter.

##  V Photon Ring Metric Tests

Is the photon ring insensitivity to the accretion model sufficient to perform
a robust test of the metric? In Fig. 7 we show the possible range of photon
ring radii for all considered spacetimes. For the first photon ring $n=1$
(left panel) all bands corresponding to charged BH spacetimes significantly
overlap with the full possible range of the first photon ring radius in Kerr
spacetimes. Hence, only some limited combinations of large charge and
particularly small (or large, for the KS spacetime) emission radius would
allow us to conclude inconsistency with the Kerr hypothesis. It makes a big
difference if we can further restrict the emission radius $r_{\rm e}$. As an
example, the middle panel of Fig. 7 shows the $n=1$ photon ring location
assuming that $r_{\rm e}=4.5\pm 0.7M$, which is consistent with the
constraints on the emission radius based on the EHT observations of M 87∗
reported in Section III. Although the critical curve is no longer within the
range of possible $n=1$ radii, it is apparent that the measurement of a first
photon ring combined with a prior on $r_{\rm e}$ from the measurement of the
direct image radius could lead to constraints on charges as well as to the
Kerr hypothesis test. The situation appears even more encouraging in the case
of the $n=2$ photon ring, right panel of Fig. 7.

So far we neglected the impact of the object mass and distance uncertainties
on the photon ring radius measurement. However, the related uncertainties may
be significant. In case of the M 87∗ assumed priors resulted in $M/D$
uncertain to about 12% [29, 27]. Clearly, this level of uncertainty could
already dominate the error budget and disable metric tests with photon rings,
as it is comparable with the whole range of possible photon ring radii for all
considered metrics shown in Fig. 7. It has been suggested by [21] that the
Kerr metric spin could be probed using a dimensionless ratio of the photon
rings radii, independently of the $M/D$ scaling factor. Similarly, we could
envision using such ratios to put constraints on the non-Kerr spacetime
models. In Fig. 8 we show the ratios of radii of second to first
($\ell_{2}/\ell_{1}$) and first to zeroth (direct image; $\ell_{1}/\ell_{0}$)
photon rings as a function of the emission radius $r_{\rm e}$. For clarity, we
only compare charged BH spacetimes with largest deviation from the
Schwarzschild spacetime, that is RN and KS with $q=1$, as well as the Kerr
spacetime with a large spin ($a/M=0.999$). The remaining problem is that the
$r_{\rm e}$ parameter is not known. However, as we discussed in Section III,
the size of the direct image $n=0$ can give us a very decent, metric
independent, estimation of the $r_{\rm e}$, within the $M/D$ uncertainties.
Hence, a combined method, using $M/D$ independent $\ell_{2}/\ell_{1}$ ratio
and $r_{\rm e}$ prior (that is not free of the $M/D$ uncertainties), could be
employed. Given the uncertainties in case of the M 87∗ observations we assume
$r_{\rm e}=4.5\pm 0.7M$, shaded region in Fig. 8. Then it is clearly possible
to constrain charges and test Kerr hypothesis, particularly with the
$\ell_{2}/\ell_{1}$ ratio. It may even be possible to constrain the spin
value, following [21].

Figure 8: Ratios of photon ring radii as a function of the emission region for
selected spacetimes as a function of the emission radius $r_{\rm e}$. The
shaded region corresponds to the $r_{\rm e}$ limits inferred from the EHT
results [22].

##  VI Conclusions

We studied structure of lensed BH images in several spherically symmetric BH
spacetimes representing deviation from the Schwarzschild spacetime
parametrized with a charge $q$, as well as in the Kerr spacetime viewed along
the spin axis. The radius of the direct image of the accretion disk is not
directly sensitive to the spacetime metric. Instead, it is sensitive to the
effective radius of the emission, which in turn may be be sensitive to the
metric, but it is also strongly dependent on the astrophysics of the accretion
flow. Based on the M 87∗ direct image radius reported by the EHT [22], we
measure the M 87∗ effective emission radius to be $4.5\pm 0.7M$, independently
of the underlying metric, indicating extreme compactness of this object.
Higher order photon rings are sensitive to the spacetime metric, and may be
used in the future to probe the metric independently (or with a strongly
bounded dependence) of the underlying model of the accretion flow. We have
demonstrated how first or second photon ring measurement could be used to
constrain BH charges in different spacetime models, particularly when such
measurement is coupled with a prior on $r_{\rm e}$ from the direct image
radius measurement.

Presented results assume a very simplified setup, ignoring the spin of charged
BHs, assuming face-on view of the accretion disk, and emission from an
equatorial thin ring. In order to properly interpret observational
measurements, a more general model needs to be employed. At the moment, theory
is a number of years ahead of the observations, so these much more detailed
and technical calculations may be performed in the future work – the relevant
conclusion is that the photon ring structure is only weakly influenced by the
assumed source model. While accounting for the inclination increases the
number of degrees of freedom in the model, it also unlocks the possibility of
performing metric tests based on non-circular photon ring geometry [11, 60].

Finally, we sensitise the readers against the straightforward interpretation
of the EHT measurements of the direct image properties (diameter, circularity)
as statements about the critical curve geometry. These are different entities,
and in order to make a connection between them, a model of the emission source
structure needs to be employed. It should be clearly stated when the inferred
properties are depending on such a model, and what are its restrictions. It is
both interesting and useful to calibrate EHT direct image observations with a
limited set of theoretical models of the accretion flow, such as GRMHD
simulations in [26], while recognizing the limitations of this framework.
Assumptions should be explicitly stated when a measurement of any physical
parameter is performed.

###### Acknowledgements.

We thank Zack Gelles, Sam Gralla, Michael Johnson, Prashant Kocherlakota,
Roman Konoplya, Alex Lupsasca, Daniel Palumbo, Bart Ripperda, and last but not
least Frederic Vincent for enlightening discussions and comments. We also
thank Alexandra Elbakyan for her contributions to the open science. This work
was supported by the Black Hole Initiative at Harvard University, which is
funded by grants from the John Templeton Foundation and the Gordon and Betty
Moore Foundation to Harvard University.

## References

  * Darwin [1959] C. Darwin, The Gravity Field of a Particle, Proceedings of the Royal Society of London Series A 249, 180 (1959).
  * Bardeen [1973] J. M. Bardeen, Timelike and null geodesics in the Kerr metric., in  _Black Holes (Les Astres Occlus)_ , edited by C. Dewitt and B. S. Dewitt (1973) pp. 215–239.
  * Luminet [1979] J.-P. Luminet, Image of a spherical black hole with thin accretion disk, Astron. Astrophys. 75, 228 (1979).
  * Falcke  _et al._ [2000] H. Falcke, F. Melia, and E. Agol, Viewing the Shadow of the Black Hole at the Galactic Center, Astrophys. J. Lett. 528, L13 (2000), arXiv:astro-ph/9912263 [astro-ph] .
  * Gralla  _et al._ [2019] S. E. Gralla, D. E. Holz, and R. M. Wald, Black hole shadows, photon rings, and lensing rings, Phys. Rev. D 100, 024018 (2019), arXiv:1906.00873 [astro-ph.HE] .
  * Narayan  _et al._ [2019] R. Narayan, M. D. Johnson, and C. F. Gammie, The Shadow of a Spherically Accreting Black Hole, Astrophys. J. Lett. 885, L33 (2019), arXiv:1910.02957 [astro-ph.HE] .
  * Vincent  _et al._ [2021] F. H. Vincent, M. Wielgus, M. A. Abramowicz, E. Gourgoulhon, J. P. Lasota, T. Paumard, and G. Perrin, Geometric modeling of M87* as a Kerr black hole or a non-Kerr compact object, Astron. Astrophys. 646, A37 (2021), arXiv:2002.09226 [gr-qc] .
  * Gralla and Lupsasca [2020] S. E. Gralla and A. Lupsasca, Lensing by Kerr black holes, Phys. Rev. D 101, 044031 (2020), arXiv:1910.12873 [gr-qc] .
  * Teo [2003] E. Teo, Spherical Photon Orbits Around a Kerr Black Hole, Gen. Relativ. Gravit. 35, 1909 (2003).
  * Johnson  _et al._ [2020] M. D. Johnson, A. Lupsasca, A. Strominger, G. N. Wong, S. Hadar, D. Kapec, R. Narayan, A. Chael, C. F. Gammie, P. Galison, D. C. M. Palumbo, S. S. Doeleman, L. Blackburn, M. Wielgus, D. W. Pesce, _et al._ , Universal interferometric signatures of a black hole’s photon ring, Science Advances 6, eaaz1310 (2020), arXiv:1907.04329 [astro-ph.IM] .
  * Gralla  _et al._ [2020] S. E. Gralla, A. Lupsasca, and D. P. Marrone, The shape of the black hole photon ring: A precise test of strong-field general relativity, Phys. Rev. D 102, 124004 (2020), arXiv:2008.03879 [gr-qc] .
  * Zakharov [1994] A. F. Zakharov, Particle capture cross sections for a Reissner–Nordström black hole, Class. Quantum Gravity 11, 1027 (1994).
  * Bambi [2013] C. Bambi, Can the supermassive objects at the centers of galaxies be traversable wormholes? The first test of strong gravity for mm/sub-mm very long baseline interferometry facilities, Phys. Rev. D 87, 107501 (2013), arXiv:1304.5691 [gr-qc] .
  * Medeiros  _et al._ [2020] L. Medeiros, D. Psaltis, and F. Özel, A Parametric Model for the Shapes of Black Hole Shadows in Non-Kerr Spacetimes, Astrophys. J. 896, 7 (2020), arXiv:1907.12575 [astro-ph.HE] .
  * Wielgus  _et al._ [2020] M. Wielgus, J. Horák, F. Vincent, and M. Abramowicz, Reflection-asymmetric wormholes and their double shadows, Phys. Rev. D 102, 084044 (2020), arXiv:2008.10130 [gr-qc] .
  * Bacchini  _et al._ [2021] F. Bacchini, D. R. Mayerson, B. Ripperda, J. Davelaar, H. Olivares, T. Hertog, and B. Vercnocke, Fuzzball Shadows: Emergent Horizons from Microstructure, arXiv e-prints , arXiv:2103.12075 (2021), arXiv:2103.12075 [hep-th] .
  * Johnson  _et al._ [2019] M. Johnson, K. Haworth, D. W. Pesce, D. C. M. Palumbo, L. Blackburn, K. Akiyama, D. Boroson, K. L. Bouman, J. R. Farah, V. L. Fish, M. Honma, T. Kawashima, M. Kino, A. Raymond, M. Silver, _et al._ , Studying black holes on horizon scales with space-VLBI, in  _Bulletin of the American Astronomical Society_ , Vol. 51 (2019) p. 235, arXiv:1909.01405 [astro-ph.IM] .
  * Pesce  _et al._ [2019] D. Pesce, K. Haworth, G. J. Melnick, L. Blackburn, M. Wielgus, M. D. Johnson, A. Raymond, J. Weintroub, D. C. M. Palumbo, S. S. Doeleman, and D. J. James, Extremely long baseline interferometry with Origins Space Telescope, in  _Bulletin of the American Astronomical Society_ , Vol. 51 (2019) p. 176, arXiv:1909.01408 [astro-ph.IM] .
  * Gurvits  _et al._ [2021] L. I. Gurvits, Z. Paragi, V. Casasola, J. Conway, J. Davelaar, H. Falcke, R. Fender, S. Frey, C. M. Fromm, C. G. Miró, M. A. Garrett, M. Giroletti, C. Goddi, J.-L. Gómez, J. van der Gucht, J. C. Guirado, Z. Haiman, F. Helmich, E. Humphreys, V. Impellizzeri, M. Kramer, M. Lindqvist, H. Linz, E. Liuzzo, A. P. Lobanov, Y. Mizuno, L. Rezzolla, F. Roelofs, E. Ros, K. L. J. Rygl, T. Savolainen, K. Schuster, T. Venturi, M. C. Wiedner, and J. A. Zensus, THEZA: TeraHertz Exploration and Zooming-in for Astrophysics, Experimental Astronomy 10.1007/s10686-021-09714-y (2021), arXiv:1908.10767 [astro-ph.IM] .
  * Broderick  _et al._ [2020] A. E. Broderick, D. W. Pesce, P. Tiede, H.-Y. Pu, and R. Gold, Hybrid Very Long Baseline Interferometry Imaging and Modeling with THEMIS, Astrophys. J. 898, 9 (2020).
  * Broderick  _et al._ [2021] A. E. Broderick, P. Tiede, D. W. Pesce, and R. Gold, Measuring Spin from Relative Photon Ring Sizes, arXiv e-prints , arXiv:2105.09962 (2021), arXiv:2105.09962 [astro-ph.HE] .
  * Event Horizon Telescope Collaboration  _et al._ [2019a] Event Horizon Telescope Collaboration _et al._ , First M87 Event Horizon Telescope Results. I. The Shadow of the Supermassive Black Hole, Astrophys. J. Lett. 875, L1 (2019a), arXiv:1906.11238 [astro-ph.GA] .
  * EHTC  _et al._ [2019a] EHTC, K. Akiyama, A. Alberdi, W. Alef, K. Asada, R. Azulay, A.-K. Baczko, D. Ball, M. Baloković, J. Barrett, and et al., First M87 Event Horizon Telescope Results. II. Array and Instrumentation, Astrophys. J. Lett. 875, L2 (2019a).
  * EHTC  _et al._ [2019b] EHTC, K. Akiyama, A. Alberdi, W. Alef, K. Asada, R. Azulay, A.-K. Baczko, D. Ball, M. Baloković, J. Barrett, and et al., First M87 Event Horizon Telescope Results. III. Data Processing and Calibration, Astrophys. J. Lett. 875, L3 (2019b).
  * Event Horizon Telescope Collaboration  _et al._ [2019b] Event Horizon Telescope Collaboration _et al._ , First M87 Event Horizon Telescope Results. IV. Imaging the Central Supermassive Black Hole, Astrophys. J. Lett. 875, L4 (2019b).
  * Event Horizon Telescope Collaboration  _et al._ [2019c] Event Horizon Telescope Collaboration _et al._ , First M87 Event Horizon Telescope Results. V. Physical Origin of the Asymmetric Ring, Astrophys. J. Lett. 875, L5 (2019c), arXiv:1906.11242 [astro-ph.GA] .
  * Event Horizon Telescope Collaboration  _et al._ [2019d] Event Horizon Telescope Collaboration _et al._ , First M87 Event Horizon Telescope Results. VI. The Shadow and Mass of the Central Black Hole, Astrophys. J. Lett. 875, L6 (2019d), arXiv:1906.11243 [astro-ph.GA] .
  * Psaltis [2019] D. Psaltis, Testing general relativity with the Event Horizon Telescope, Gen. Relativ. Gravit. 51, 137 (2019), arXiv:1806.09740 [astro-ph.HE] .
  * Gebhardt  _et al._ [2011] K. Gebhardt, J. Adams, D. Richstone, T. R. Lauer, S. M. Faber, K. Gültekin, J. Murphy, and S. Tremaine, The Black Hole Mass in M87 from Gemini/NIFS Adaptive Optics Observations, Astrophys. J. 729, 119 (2011), arXiv:1101.1954 [astro-ph.CO] .
  * Psaltis  _et al._ [2020] D. Psaltis, L. Medeiros, P. Christian, F. Özel, and EHT Collaboration, Gravitational Test beyond the First Post-Newtonian Order with the Shadow of the M87 Black Hole, Phys. Rev. Lett. 125, 141104 (2020), arXiv:2010.01055 [gr-qc] .
  * Kocherlakota  _et al._ [2021] P. Kocherlakota, L. Rezzolla, H. Falcke, C. M. Fromm, M. Kramer, Y. Mizuno, A. Nathanail, H. Olivares, Z. Younsi, and EHT Collaboration, Constraints on black-hole charges with the 2017 EHT observations of M87*, Phys. Rev. D 103, 104047 (2021), arXiv:2105.09343 [gr-qc] .
  * Gralla [2021] S. E. Gralla, Can the EHT M87 results be used to test general relativity?, Phys. Rev. D 103, 024023 (2021), arXiv:2010.08557 [astro-ph.HE] .
  * Doeleman  _et al._ [2008] S. S. Doeleman, J. Weintroub, A. E. E. Rogers, R. Plambeck, R. Freund, R. P. J. Tilanus, P. Friberg, L. M. Ziurys, J. M. Moran, B. Corey, K. H. Young, D. L. Smythe, M. Titus, D. P. Marrone, R. J. Cappallo, D. C. J. Bock, G. C. Bower, _et al._ , Event-horizon-scale structure in the supermassive black hole candidate at the Galactic Centre, Nature 455, 78 (2008), arXiv:0809.2442 [astro-ph] .
  * Narayan  _et al._ [1995] R. Narayan, I. Yi, and R. Mahadevan, Explaining the spectrum of Sagittarius A∗ with a model of an accreting black hole, Nature 374, 623 (1995).
  * Di Matteo  _et al._ [2003] T. Di Matteo, S. W. Allen, A. C. Fabian, A. S. Wilson, and A. J. Young, Accretion onto the Supermassive Black Hole in M87, Astrophys. J. 582, 133 (2003), arXiv:astro-ph/0202238 [astro-ph] .
  * Rees  _et al._ [1982] M. J. Rees, M. C. Begelman, R. D. Blandford, and E. S. Phinney, Ion-supported tori and the origin of radio jets, Nature 295, 17 (1982).
  * Narayan and Yi [1994] R. Narayan and I. Yi, Advection-dominated Accretion: A Self-similar Solution, Astrophys. J. Lett. 428, L13 (1994), arXiv:astro-ph/9403052 [astro-ph] .
  * Abramowicz  _et al._ [1995] M. A. Abramowicz, X. Chen, S. Kato, J.-P. Lasota, and O. Regev, Thermal Equilibria of Accretion Disks, Astrophys. J. Lett. 438, L37 (1995), arXiv:astro-ph/9409018 [astro-ph] .
  * Yuan and Narayan [2014] F. Yuan and R. Narayan, Hot Accretion Flows Around Black Holes, Annu. Rev. Astron. Astrophys. 52, 529 (2014), arXiv:1401.0586 [astro-ph.HE] .
  * Bronzwaer  _et al._ [2021] T. Bronzwaer, J. Davelaar, Z. Younsi, M. Mościbrodzka, H. Olivares, Y. Mizuno, J. Vos, and H. Falcke, Visibility of black hole shadows in low-luminosity AGN, Mon. Notices Royal Astron. Soc. 501, 4722 (2021), arXiv:2011.00069 [astro-ph.HE] .
  * Novikov and Thorne [1973] I. D. Novikov and K. S. Thorne, Astrophysics of black holes., in  _Black Holes (Les Astres Occlus)_ (1973) pp. 343–450.
  * Rybicki and Lightman [1979] G. B. Rybicki and A. P. Lightman,  _Radiative processes in astrophysics_ (1979).
  * Falcke  _et al._ [2004] H. Falcke, E. Körding, and S. Markoff, A scheme to unify low-power accreting black holes. Jet-dominated accretion flows and the radio/X-ray correlation, Astron. Astrophys. 414, 895 (2004), arXiv:astro-ph/0305335 [astro-ph] .
  * Gates  _et al._ [2020] D. E. A. Gates, S. Hadar, and A. Lupsasca, Maximum observable blueshift from circular equatorial Kerr orbiters, Phys. Rev. D 102, 104041 (2020), arXiv:2009.03310 [gr-qc] .
  * Mizuno  _et al._ [2018] Y. Mizuno, Z. Younsi, C. M. Fromm, O. Porth, M. De Laurentis, H. Olivares, H. Falcke, M. Kramer, and L. Rezzolla, The current ability to test theories of gravity with black hole shadows, Nature Astronomy 2, 585 (2018), arXiv:1804.05812 [astro-ph.GA] .
  * Kocherlakota and Rezzolla [2020] P. Kocherlakota and L. Rezzolla, Accurate mapping of spherically symmetric black holes in a parametrized framework, Phys. Rev. D 102, 064058 (2020), arXiv:2007.15593 [gr-qc] .
  * Dadhich  _et al._ [2000] N. Dadhich, R. Maartens, P. Papadopoulos, and V. Rezania, Black holes on the brane, Phys. Lett. B 487, 1 (2000), arXiv:hep-th/0003061 [hep-th] .
  * Zakharov [2014] A. F. Zakharov, Constraints on a charge in the Reissner-Nordström metric for the black hole at the Galactic Center, Phys. Rev. D 90, 062007 (2014), arXiv:1407.7457 [gr-qc] .
  * bar [1968] _Presented at GR5, Tbilisi, U.S.S.R._ (ACM, 1968).
  * Ayón-Beato and García [1998] E. Ayón-Beato and A. García, Regular Black Hole in General Relativity Coupled to Nonlinear Electrodynamics, Phys. Rev. Lett. 80, 5056 (1998), arXiv:gr-qc/9911046 [gr-qc] .
  * Hayward [2006] S. A. Hayward, Formation and Evaporation of Nonsingular Black Holes, Phys. Rev. Lett. 96, 031103 (2006), arXiv:gr-qc/0506126 [gr-qc] .
  * Abdujabbarov  _et al._ [2016] A. Abdujabbarov, M. Amir, B. Ahmedov, and S. G. Ghosh, Shadow of rotating regular black holes, Phys. Rev. D 93, 104004 (2016), arXiv:1604.03809 [gr-qc] .
  * Stuchlík and Schee [2019] Z. Stuchlík and J. Schee, Shadow of the regular Bardeen black holes and comparison of the motion of photons and neutrinos, European Physical Journal C 79, 44 (2019).
  * Konoplya [2020] R. A. Konoplya, Quantum corrected black holes: Quasinormal modes, scattering, shadows, Phys. Lett. B 804, 135363 (2020), arXiv:1912.10582 [gr-qc] .
  * Gibbons and Maeda [1988] G. W. Gibbons and K.-I. Maeda, Black holes and membranes in higher-dimensional theories with dilaton fields, Nuclear Physics B 298, 741 (1988).
  * Garfinkle  _et al._ [1991] D. Garfinkle, G. T. Horowitz, and A. Strominger, Charged black holes in string theory, Phys. Rev. D 43, 3140 (1991).
  * Wei and Liu [2013] S.-W. Wei and Y.-X. Liu, Observing the shadow of Einstein-Maxwell-Dilaton-Axion black hole, J. Cosmol. Astropart. Phys. 2013 (11), 063, arXiv:1311.4251 [gr-qc] .
  * Bardeen  _et al._ [1972] J. M. Bardeen, W. H. Press, and S. A. Teukolsky, Rotating Black Holes: Locally Nonrotating Frames, Energy Extraction, and Scalar Synchrotron Radiation, Astrophys. J. 178, 347 (1972).
  * Chael  _et al._ [2021] A. Chael, M. D. Johnson, and A. Lupsasca, Observing the Inner Shadow of a Black Hole: A Direct View of the Event Horizon, arXiv e-prints , arXiv:2106.00683 (2021), arXiv:2106.00683 [astro-ph.HE] .
  * Paugnat  _et al._ [2021] H. Paugnat, F. H. Vincent, A. Lupsasca, S. E. Graala, and M. Wielgus, Variability of the photon ring shape in Kerr black hole test, in prep. (2021).


# Three-dimensional shape and connectivity of physical networks

Luka Blagojević<EMAIL_ADDRESS>Department of Network and Data Science,
Central European University, Vienna, Austria  Márton Pósfai<EMAIL_ADDRESS>Department of Network and Data Science, Central European University, Vienna,
Austria

###### Abstract

Data describing the three-dimensional structure of physical networks is
increasingly available, leading to a surge of interest in network science to
explore the relationship between the shape and connectivity of physical
networks. We contribute to this effort by standardizing and analyzing 15 data
sets from different domains. Treating junction points as nodes and connections
between them as links, we divide the networks into three categories: lattice-
like networks, trees, and linked trees. We find that the degree distribution
of physical networks is bounded, with most nodes having degree one or three.
Characterizing the physical properties of links, we show that links have an
elongated shape and tend to follow a nearly straight trajectory, while a small
fraction of links follow a winding path. These typical node and link
properties must be reflected by physical network models. We also measure how
confined a link is in space by comparing its trajectory to a randomized null-
model, showing that links that are central in the abstract network tend to be
physically confined by their neighbors. The fact that the shape and
connectivity of the physical networks are intertwined highlights that their
three-dimensional layout must be taken into account to understand the
evolution and function of physical networks.

##  I Introduction

Recently available detailed maps of physical networks provide an opportunity
to systematically investigate the relationship between their shape and their
network properties. For example, open-data large-scale experimental projects
provide a neuron-level mapping of biological neural networks [34, 26], high-
throughput MRI measurements map out detailed vascular networks [10], or
mycelia mapping projects explore local morphological and mechanical properties
of fungi by representing them as networks of filaments [11]. The nodes and
links of these physical networks have two distinguishing features: (i) they
are characterized by a complex three-dimensional shape and (ii) they
physically interact with each other, for example, nodes and links obey volume
exclusion, i.e., they cannot overlap. To understand how physicality affects
network evolution and function, we must extend the toolset of network science
to take into account these features.

Recent work investigated artificial spatial embeddings of complex networks
that obey volume exclusion [6], the entanglement of physical links [15],
models of physical network growth [20, 19], and the effect of physical shape
on the dynamics on networks [19]. However, systematic exploration of the
three-dimensional shape, the network properties and the relationship between
them in real networks is still lacking. Such exploration is hampered by the
lack of standardized representation. First, there is the technical difficulty
that experimental maps of physical networks, like neural or molecular
networks, are collected, processed, and analyzed with domain-specific
methodology. Therefore, any investigation of physical networks must be
preceded by the time-consuming and computationally burdensome task of data
pre-processing. Second, even seemingly simple questions like what is a node
and a link in a physical network carry a level of ambiguity: a physical
network is a continuous object in space; to represent it as a network, we must
discretize it into nodes and links. The definition of nodes and links in turn
affects, for example, what properties of the network we can study or the right
choice of null models.

Here, we compile and standardize 15 data sets from various domains. Each of
these physical networks is composed of tube-like objects bound together at
junction points; motivating us to treat the junction points as nodes and the
tubes connecting them as physical links. We characterize both physical shape
and the abstract network structure and the correlations between them. For
this, we calculate standard descriptors such as the degree distribution of the
abstract network or the fractal dimension of the layout. We also introduce a
measure of link confinement to understand the role of volume exclusion, which
compares the physical links to a null model that randomizes link routing. The
remainder of the paper is organized as follows: In the next section, we
describe the data sets we collected and their standardization. In. Secs.
III-V, we analyze the data sets’ abstract network properties, their physical
shape, and the emergent correlations between network and shape. Finally, Sec.
VI provides a brief discussion.

Figure 1: Physical networks. (a) A physical network composed of three tube-
like links bound together at a single junction point. (b) The skeleton
representation of (a) approximates the original structure as a collection of
vertices (black points) and edges (colored segments). (c) The combinatorial or
abstract network of (a) captures the connectivity of the system without the
physical structure: nodes represent junction points and terminal points, with
a link between two the nodes that they are directly connected by a physical
link. (d) The skeleton representation allows us to approximate the original
volume of real physical networks – e.g., vascular network (top) and
mitochondrial network (bottom). (e) We compiled a set of 15 physical networks
from various domains. The size of the networks varies greatly: the number of
skeleton segments $N_{\textnormal{seg}}$ capturing the shape of each networks
spans 4 orders of magnitude.

##  II Data

Our goal is to systematically study the three-dimensional layout of physical
networks and to understand the relationship between their physical properties
and their network structure. For this, we collected 15 data sets from various
domains, including individual neurons [23] [3] [12] [30], biological neural
networks [25], plant roots [17], vascular networks [10], a mitochondrial
network [31], and the imprint of an anthill (See SI Sec. I). Before any
analysis, however, we must uniformly represent and standardize these data
sets. In the following sections, we propose the use of a labeled skeleton
representation, which efficiently captures both the physical shape and
abstract network of physical networks.

###  II.1 Skeleton representation

Experimental imaging techniques that are used to capture the shape of physical
networks, such as scanning electron microscopy or magnetic resonance imaging,
typically output a three-dimensional image composed of voxels. Hence, this
voxel representation is the most accurate description available of network
layouts. A three-dimensional voxel image, however, is difficult to handle both
computationally and analytically; therefore a more compact representation of
the data is needed.

Physical networks are typically composed of tube-like objects bound together
at junction points (Fig. 1a), making them suitable to be approximated by a
series of straight segments and radii of the network at the endpoints of the
segments (Fig. 1b). The process of creating these segments from a raw data
format is called skeletonization [22], often employed in the fields of
computer graphics [29] and neuroscience [27, 2]. Formally, a skeleton
representation ${\mathcal{S}}$ of a physical network is a graph whose set of
vertices ${\mathcal{V}}$ correspond to points in space and set of edges
${\mathcal{E}}$ correspond to segments connecting point pairs. Therefore each
vertex $i\in{\mathcal{V}}$ has a set of coordinates ${\bf
r}_{i}=(x_{i},y_{i},z_{i})$ and a radius $\rho_{i}$ associated to it. Figure
1b shows the skeleton of Fig. 1a, the black circles represent the skeleton
vertices.

The radius $\rho_{i}$ captures the thickness of the physical network at each
vertex $i$. Therefore, we can approximate the original occupied volume of the
network as a union of truncated cones: for each edge $(i,j)$ in
${\mathcal{S}}$ we add a truncated cone with axis $({\bf r}_{i},{\bf r}_{j})$
and parallel faces with radii $\rho(i)$ and $\rho(j)$ which has volume:

$V_{\textnormal{seg}}(i,j)=\frac{1}{3}\cdot\pi\cdot(\rho_{i}^{2}+\rho_{j}^{2}+\rho_{i}\cdot\rho_{j})\cdot|{\bf
r}_{i}-{\bf r}_{j}|.$ (1)

Three-dimensional physical network data obtained from experiments is routinely
skeletonized, and the skeleton of the network is published together with the
raw data. In fact, all but one of the 15 data sets that we study here was
skeletonized by the original authors, the only exception is the anthill
imprint provided to us as a surface mesh which we skeletonized using the
Tangent-ball [33] algorithm from the Skeletor Python module [2]. The
experimental setup and the choices made during the skeletonization may affect
our analysis which is performed on the skeleton. For example, increasing the
number of skeleton segments $N_{\textnormal{seg}}$ increases how well the
skeleton approximates the original shape of the network. However, increasing
$N_{\textnormal{seg}}$ also increases the size of the data set, hence
increasing the computational burden of the analysis. To improve the uniformity
of the data sets, we perform two pre-processing steps:

  1. 1.

Merging segments: If two consecutive segments $({\bf r}_{i},{\bf r}_{j})$ and
$({\bf r}_{j},{\bf r}_{k})$ are approximately parallel to each other, we merge
them into a single segment $({\bf r}_{i},{\bf r}_{k})$ (See SI Sec. II B).

  2. 2.

Skeleton healing: Due to noisy data, a skeleton may be disconnected even when
it represents a single continuous object in reality. For example, skeleton of
a neuron may appear to have multiple components. To remedy this, we add a
segment to connect the two nearest skeleton vertices from two disconnected
components. We repeat this step until the skeleton becomes connected.

Following these pre-processing steps, the number of segments
$N_{\textnormal{seg}}$ in a skeleton, as shown in Fig. 1e, spans three orders
of magnitude from mitochondrial networks, which have approximately $10^{3}$
segments, to fruit fly neural networks which have up to $10^{6}$ segments.

###  II.2 Network structure

A physical network is a continuous object embedded in Euclidean space, to
characterize this object as a network we must separate it into discrete nodes
and links. For this, first note that all 15 data sets that we collected can be
seen as a collection of tube-like objects bound together at junction points.
Motivated by this observation, we define the junction and terminal points of
the tubes as physical nodes and the non-branching tubes pairwise connecting
these terminal and junction points as physical links. A motivation for this
definition is that cutting a physical link (i.e., a tube) at any point along
its length causes the same disruption to the connectivity of the network.

More formally, for a skeleton representation ${\mathcal{S}}$ of a physical
system, we define each physical node to correspond to a skeleton vertex $i$
with degree $k(i)\neq 2$, and each physical link to corresponds to a path in
the skeleton given by the ordered set
${\mathcal{T}}(i_{0},i_{l})=[(i_{0},i_{1}),(i_{1},i_{2}),\ldots,(i_{l-1},i_{l})]$,
such that $k(i_{0}),k(i_{l})\neq 2$ and $k(i_{j})=2$ for $j=1,2,\ldots,l-1$.

With the above definition of nodes and links, we can talk about the abstract
or combinatorial network ${\mathcal{G}}$ of the system which captures its
connectivity without the physical structure. The skeleton ${\mathcal{S}}$ is
one possible physical realization of the abstract network ${\mathcal{G}}$;
however, there are many possible physical realizations of the same
${\mathcal{G}}$. In general, we are interested in understanding the
relationship between the physical layout captured by ${\mathcal{S}}$ and the
network structure captured by ${\mathcal{G}}$.

As an example consider the physical network shown in Fig. 1a which consists of
three tubes bound together at a single junction point. Its skeleton
representation (Fig. 1b), therefore, has three vertices with degree 1
corresponding to the terminal points, one vertex with degree $k=3$
corresponding to the junction point, and several vertices with degree $k=2$
tracing the trajectory of the tubes. This means that the network consists of
four physical nodes and three physical links, and its abstract network is a
star (Fig. 1c).

Finally, note that for a given skeleton ${\mathcal{S}}$, the above definition
of physical nodes and links is not the only viable definition. For example, in
a neural network, it is natural to treat a neuron as a physical node and
synapses between them as links, as individual neurons can have complex three-
dimensional shapes that can be represented by a skeleton itself. More
generally, subgraphs of a large ${\mathcal{S}}$ may represent functional units
in a physical network and it can be useful to treat these functional units as
physical nodes [19]. Note, however, that our definition of the abstract
network ${\mathcal{G}}$ provides the most detailed picture of the system’s
connectivity, and other definitions can be thought of as coarse-grained
versions of ${\mathcal{G}}$.

With the skeleton representation and the definition of the abstract network at
hand, we are in the position to start our analysis. In the following sections,
we first explore the structure of the abstract networks of the 15 data sets,
then we continue by characterizing their physical properties, and finally, we
investigate the relation between the two.

##  III Abstract network properties

The above definition of physical nodes and links allows us to explore the
properties of the abstract networks capturing the connectivity of physical
networks without their three-dimensional structure; we focus on the degree
distribution and motif frequencies.

The nodes in our physical networks are terminal and junction points, meaning
that by construction nodes cannot have degree $k=2$, only degree $k=1$ or
$k>2$. Figure 2a shows the degree distribution $p(k)$ of all 15 networks, our
main observation is that most nodes have degree $k=1$ or $k=3$, and nodes with
a larger degree are exceedingly rare, i.e., $p(1)+p(3)\approx 1$ This means
that junctions tend to be bifurcation points along the tube-like physical
links making up the network. This is in line with previous empirical
observations and theoretical predictions for neurons [18, 7] and transport
networks [13, 8]. The observed narrow degree distribution is in contrast with
degree-heterogeneous networks typically in the focus of network science and
should be accounted for by mathematical models of physical networks [20].

In graph topological terms, 7 out of 15 collected physical networks are trees:
the individual neurons, the anthill imprint, and the plant roots. Each
bifurcation point with degree $k=3$ in a tree creates one new leaf node with
degree $k=1$; therefore the fact that these networks are trees together with
the observation that most junctions are bifurcation points, completely
determines their degree distribution as $p(1)\approx p(3)\approx 1/2$. The
remainder of the networks contain cycles: The vascular networks have no
terminal points, apart from a few appearing due to finite sample size;
therefore are almost completely composed of $k=3$ nodes. The mitochondrial
network representing a network of molecular strands has the highest fraction
of $k>3$ nodes. Finally, the fruit fly neural networks represent a collection
of individual neurons which are trees, bound together by synapses, and their
degree distribution closely resembles that of trees.

To explore the local loop structure of the networks, we calculate the
abundance of observed 4-node motifs. Namely, we focus on two motifs: the star
motif and the 4-cycle. To quantify their abundance, we calculate their Z-score
compared to their degree preserving randomized counterparts:

$z_{\textrm{s/c}}=\frac{n_{\textrm{s/c}}-\left\langle
n_{\textrm{s/c}}\right\rangle}{\sigma_{\textrm{s/c}}},$ (2)

where $n_{\textrm{s/c}}$ is the number of occurrences of the star and 4-cycle
motifs in the original networks, and the expected value $\left\langle
n_{\textrm{s/c}}\right\rangle$ and standard deviation $\sigma_{\textrm{s/c}}$
is estimated by creating 200 independent randomizations[16], while keeping the
degree sequence fixed. Figure 2b shows the scatter plot of $z_{\textrm{s}}$
and $z_{\textrm{c}}$ for the 15 networks. As expected, in tree networks
containing no cycles the star motif is slightly over-represented
$z_{\textrm{s}}>0$, and the cycle motif is slightly underrepresented
$z_{\textrm{s}}<0$, while for networks containing cycles, we find the
opposite. The highest abundance of 4-cycles is observed for the networks
representing different brain regions of the fruit fly brain.

Based on the degree distributions and motif profile of the abstract networks,
the physical networks fall into three broad categories: (i) topological trees
(the individual neurons, the root systems, and the anthill tunnel imprint),
(ii) lattice-like networks that are characterized by a loopy structure and few
terminal points (the vascular networks and the mitochondrial network), and
(iii) linked trees which are a collection of trees bound together by
additional links (the fruit fly brain regions).

Figure 2: Abstract network properties (a) The abstract networks are composed
of terminal and bifurcation points; therefore, their degree distribution is
mostly concentrated on $k=1$ and $k=3$. More specifically, lattice-like
networks, such as vascular and mitochondrial networks, are mostly made up of
branching nodes (degree $k=3$). For the rest of the network, nodes with $k=1$
and $k=3$ are approximately evenly split, as expected for tree networks. (b)
We calculate the $z$-scores of four-node star and cycle motifs of the original
networks compared to random networks with the same degree distribution.
Markers are the average of $n_{\textnormal{t}}=200$ independent
randomizations, error bars representing the standard error of the mean are
smaller than the marker size. Both axes are symlog axes, i.e., linear from -1
to 1, log otherwise for both positive and negative values. The star motif is
slightly over and the cycle motif is slightly underrepresented in tree
networks (green markers), as expected. In contrast. In contrast, stars are
under- and cycles are over-represented in both lattices (blue markers) and
linked trees (red markers). The high $z$-scores observed for the 4-cycle
indicate a lattice-like structure.

##  IV Physical properties

In the previous section, we measured properties of the abstract network
capturing connectivity in the 15 physical networks. We continue our
investigation focusing on the physical properties of the system: we
characterize the three-dimensional shape of the system without considering the
abstract network.

###  IV.1 Space filling and fractal dimension

A fundamental property of a physical network is its density or space-filling,
i.e., the amount of volume it occupies from the available space. We expect
that networks that are tightly packed in space are strongly affected by
physicality [6, 21, 19]. Here, we measure a local space-filling of the 15
physical networks by dividing the bounding box of each physical network into
cubes with side length $l_{b}$ and calculating the space-filling within each
cube $i$ as

$\phi(i)=\frac{V_{\textrm{occ}}(i)}{l_{b}^{3}},$ (3)

where $V_{\textrm{occ}}(i)$ is the volume of the intersection of the networks
and cube $i$. The distribution of $\phi(i)$ depends on the choice of $l_{b}$;
therefore, to ensure consistency, we set $l_{b}$ separately for each dataset
such that the cubes are most likely to contain one or two links (For details,
check SI. Sec LABEL:physical_density_SI). Figure 3 shows that the distribution
of $\phi$ for all 15 physical networks is not peaked around a single value,
but the physical networks fill out the space unevenly with both densely packed
and sparse regions. Therefore we expect that physicality will also affect the
network structure unevenly: volume exclusion can limit the number and shape of
links in dense regions.

To further characterize the shape of the networks, we calculate their box-
counting fractal dimension $D_{\textnormal{f}}$, which compactly describes the
space-filling of a physical object on multiple scales [28, 9, 32], and is
widely used to characterize the shape of complex biological systems [32, 4].
Possible values of $D_{\textnormal{f}}$ for connected networks range between 1
and the embedding dimension $D=3$, Fig. 3b shows $D_{\textnormal{f}}$ for the
15 networks, we find that both lattices and linked trees have
$D_{\textnormal{f}}\gtrapprox 2$, while trees are typically characterized by
$D_{\textnormal{f}}\lessapprox 2$, except for the anthill imprint. The fractal
scaling spans at least two orders of magnitude of length scales (see SI
Sec.LABEL:fractal_dimension_SI), again pointing towards regions of high and
low physical density at different resolution levels.

Figure 3: Space-filling and fractal dimension. (a) The distribution of local
space-filling $\phi$ shows the coexistence of mostly sparse and fewer dense
regions. (b) The fractal dimension of physical networks $D_{\textrm{f}}$
ranges between $1$ and the embedding dimension $D=3$. Linked trees (red) and
lattices (blue) have fractal dimension values between $D_{\textrm{f}}\approx
2.0$ and $D_{\textrm{f}}\approx 2.3$, while trees (green) have more variation,
ranging from $D_{\textrm{f}}\approx 1.5$ to $D_{\textrm{f}}\approx 2.0$
(except for from the anthill imprint). We estimate $D_{\textrm{f}}$ using the
box-counting method, the error bars indicate the standard deviation of the
local scaling estimates (see SI Sec.LABEL:fractal_dimension_SI).

###  IV.2 Link volume and shape

The distribution of space-filling and the fractal dimension characterize the
shape of physical networks as a whole. In this section, we continue by
quantifying the shape of individual links; we focus on their volume and their
straightness.

The skeleton describing the three-dimensional shape of a physical network is
composed of straight segments connecting pairs of skeleton vertices and a
radius associated to each skeleton vertex, allowing us to approximate the
volume belonging to a segment in Eq. (1) as a truncated cone. The total volume
of a link $(i,j)$ is then:

$V_{\textrm{link}}(i,j)=\frac{1}{V_{\textrm{total}}}\sum_{(v,w)\in(i,j)}{V_{\textnormal{seg}}(v,w)},$
(4)

where $V_{\textnormal{seg}}(v,w)$ is the volume of each segment tracing the
link $(i,j)$. We normalize the link volume by the total volume of the network
$V_{\textrm{total}}$, setting the unit of measurement. We find that linked
trees, or the fruit fly neural networks, have consistently high link volume
heterogeneity, as their distributions span 6 to 10 orders of magnitude, which
is higher compared to most lattices and trees (Fig.4a).

We also measure the aspect ratio of physical links
$a(i,j)=\rho_{\textrm{link}}(i,j)/l_{\textrm{link}}(i,j)$, where
$\rho_{\textrm{link}}$ is the average radius and $l_{\textrm{link}}(i,j)$ is
the length of link $(i,j)$. We find the largest average aspect ratios for one
of the fruit fly networks and the anthill ($med(a)=\frac{1}{3}$), while
$med(a)$ is substantially lower for other data sets (See SI Sec. I A).
Overall, this confirms that physical links are elongated tube-like objects.

Since physical links are tube-like objects, we can capture most of their shape
by characterizing their one-dimensional trajectory. Here, we calculate the
deviation of the link trajectories from a straight line, quantifying how
curved a link is. For this we rely on a measure of straightness introduced
originally in the context of geographical networks [5], namely we calculate
the complimentary straightness for each link

$\bar{S}(i,j)=1-\frac{\lvert{\bf r}_{i}-{\bf
r}_{j}\rvert}{l_{\textnormal{link}}(i,j)},$ (5)

where $\lvert{\bf r}_{i}-{\bf r}_{j}\rvert$ is the Euclidean distance between
nodes $i$ and $j$ and
$l_{\textnormal{link}}(i,j)=\sum_{(v,w)\in(i,j)}\lvert{\bf r}_{v}-{\bf
r}_{w}\rvert$ is the length of the physical link $(i,j)$. The complimentary
straightness $\bar{S}(i,j)$ is 0 if the physical link is straight and close to
1 if it follows a winding trajectory much longer than the straight path
between the two points.

Calculating the median of complimentary link straightness distribution,
$\text{med}(\bar{S})$, reveals that links in all of the 15 physical networks
tend to follow a trajectory close to a straight line: most networks have
$\text{med}(\bar{S})\approx 0.1$, meaning that the length of links is most
often less than 10% longer than the optimal straight trajectory. Similarly to
link volume heterogeneity, linked trees tend to cluster together and are among
the networks with the straightest links with $\text{med}(\bar{S})=0.05$.
Reference [21] introduced random linear physical networks, a minimal model
that constructs physical network from straight cylinders. The fact that we
observed an abundance of straight or close-to-straight links lends support for
using such linear physical network models to understand the role of
physicality in real networks. Note, however, that although most links are
close to straight, the distribution of $\bar{S}$ is right-skewed as seen in
Fig. 4, which points to a smaller fraction of links that significantly deviate
from a straight trajectory. (See SI Sec. III and SI Sec. IV D).

Finally, we computed the correlations between link straightness $S(i,j)$ and
the total link length $l_{\textrm{link}}(i,j)$ and volume
$V_{\textrm{link}}(i,j)$ for each dataset using Kendall’s rank correlation
coefficient $\tau$ [1] Figure 4b shows that for all networks, there is a
positive rank correlation $\tau>0$ between $\bar{S}(i,j)$ and
$l_{\textrm{link}}(i,j)$, indicating that longer links tend to follow a more
winding path. We also observe a positive correlation $\tau>0$ between
$\bar{S}(i,j)$ and $V_{\textrm{link}}(i,j)$, since longer links tend to have
larger volume. The only exception to this are the fruit fly neural networks,
for which $\bar{S}(i,j)$ and $V_{\textrm{link}}(i,j)$ are negatively
correlated $\tau<0$. To explain this, note that fruit fly neural networks are
composed of neurons that have large somas, which are represented in the
dataset as short, yet high-volume segments, hence the largest volume of
physical links are short and straight (See SI Sec. I G).

Figure 4: Link shape, length and volume. (a) Box plots of $\bar{S}$ and
$V_{\textrm{link}}$, where the vertical lines indicate the median, the boxes
span the $25$th-$75$th percentile range, and the whiskers extend from the
$0.1$th to the $99.9$th percentile. We find that the bulk of the distributions
fall in the range between $\textrm{med}(\bar{S})=0$ and
$\textrm{med}(\bar{S})=0.1$, indicating that the networks are mostly composed
of close-to-straight links. Link volume $V_{\textrm{link}}$ distributions span
a wider range for linked trees, and are the most narrow for lattices. (b) By
computing Kendall rank correlation $\tau$ between link complementary
straightness $\bar{S}$, link volume $V_{\textrm{link}}$ and link trajectory
length $l_{\textrm{link}}$, we observe a consistent trend of $\tau>0$, meaning
that longer and more voluminous links tend to have more winding paths. This
trend is only reversed for the fruit fly networks, which have $\tau<0$ between
$\bar{S}$ and $V_{\textrm{link}}$: neuron somas are represented as high-
volume, straight links composed of a small number of skeleton segments.

###  IV.3 Link confinement

In the previous sections, we found that space-filling and link properties are
heterogeneously distributed: most regions of space are sparse and most links
are straight, yet there exists dense regions of the network and a small
fraction of links follow paths that deviate from a straight line
significantly. This suggests that volume exclusion or other repulsive physical
interactions may also play an uneven role in shaping the network. To further
investigate this hypothesis we devise a quantity that captures the confinement
of a link by other components of a network.

A link $(i,j)$ following a trajectory ${\mathcal{T}}(i,j)$ in a real physical
network obeys volume exclusion: it does not overlap with other links. Our
strategy to quantify the role of repulsive forces that may shape
${\mathcal{T}}(i,j)$ is to calculate the number of overlaps with other links
for a random ensemble of synthetic links that follow similar trajectories to
${\mathcal{T}}(i,j)$. If the synthetic links typically overlap with many other
links, the trajectory ${\mathcal{T}}(i,j)$ is an outlier and must be shaped by
forces not captured by the random ensemble.

The trajectory of a physical link $(i,j)$ in our datasets is given by the
ordered set ${\mathcal{T}}(i,j)$ composed of a series of oriented three-
dimensional segments. To generate the random trajectory
${\mathcal{T}}_{\textrm{r}}(i,j)$, we shuffle the order of the segments while
maintaining their orientation and length, creating a uniform random
permutation of ${\mathcal{T}}(i,j)$. The randomization preserves the endpoints
and the total length of the link, but otherwise removes any correlation
between subsequent segments; therefore, the possible link trajectories
${\mathcal{T}}_{\textrm{r}}(i,j)$ have the same complimentary straightness
$\bar{S}$ as the original link. Next, we estimate $I(i,j;l,k)$, the expected
number of intersections between the randomized link
${\mathcal{T}}_{\textrm{r}}(i,j)$ and each other non-randomized link
${\mathcal{T}}(k,l)$. To quantify the confinement of the link $(i,j)$, we sum
up the expected number of intersections with other links:

$C(i,j)=\sum_{k,l\neq i,j}I(i,j;l,k)+I(l,k;i,j),$ (6)

where the first term corresponds to intersections when link $(i,j)$ is
randomized, and the second term corresponds to intersections when link $(l,k)$
is randomized. Note that the summation in Eq. (6) excludes links that share an
endpoint with $(i,j)$. We do this to exclude trivial intersections form the
count, since adjacent links $(i,j)$ and $(j,k)$ necessarily overlap at the
junction point $j$ even for non-randomized link trajectories.

The procedure of calculating $C(i,j)$ is illustrated by Fig. 5a: we start with
a link (green) surrounded by two neighboring links (red and blue). The figure
shows $n_{\textnormal{t}}=2$ randomization trials of the green link: In trial
one, the randomized link intersects the red link, but not the blue. In trial
two, the randomization creates an intersection with the blue link; the two
links, however, are adjacent (they share a junction point), hence the
intersection is not counted. In this particular example, the contribution to
the link confinement measure from the randomization of the green link will be:

$\displaystyle I(\textrm{green},\textrm{red})$
$\displaystyle=\frac{0+1}{2}=0.5$ (7) $\displaystyle
I(\textrm{green},\textrm{blue})$ $\displaystyle=\frac{0+0}{2}=0$ (8)

To complete the calculation of $C(\textrm{green})$, we also need to randomize
the red link to estimate $I(\textrm{red},\textrm{green})$ in the same manner,
while $I(\textrm{blue},\textrm{green})=0$ by definition. Finally, the
confinement of the green link is obtained by summing up the contributions,
i.e.,
$C(\textrm{green})=I(\textrm{green},\textrm{red})+I(\textrm{red},\textrm{green})$.

Collision detection between link trajectories is a computationally expensive
task, in practice we randomize each link $n_{\textrm{t}}=20$ times (and
$n_{\textrm{t}}=6$ for the fruit_fly_1 network) and we rely on an efficient
collision detection algorithm leveraging kd-trees [24] (See SI Sec. IV F for
details).

Figure 5b shows a large variation in link confinement $C$. Across all
networks, physical links are characterized by $C\approx 0$, indicating that
these links are not affected by the physical proximity of other components of
the network. However, we also find highly confined links with $C>10$ and even
$C>100$ expected intersections, suggesting again that physicality tends to
play a heterogeneous role in forming networks. In particular, the linked tree
networks typically have heavy-tailed link confinement distributions (See SI
Sec. V). In terms of absolute counts, the linked trees or the fruit fly neural
networks show the highest values of link confinement. This can be explained by
the fact that these networks are composed of multiple neurons, hence in these
networks we have more complete information about the physical environment of
the links, compared to networks that describe single neurons.

Figure 5: Link confinement. (a) To quantify volume exclusion at the link
level, we introduce link confinement $C$, the expected number of intersections
after randomizing the link trajectory. To estimate $C$ of the green link, we
randomize its trajectory twice. In trial #1, the green link intersects the red
link, while in in trial #2, it intersects the blue link. The blue link,
however, shares and endpoint with the green link; therefore their intersection
is not counted. (b) The distribution of $C$ for each physical network. There
are many links with link confinement values close or equal to $C=0$, and
typical values of link confinement are around $C\sim 10$, and only for the
fruit fly networks, there are links with link confinement values $C>100$,
indicating highly confined links for the linked tree networks.

##  V Link confinement correlation profiles

In the previous section, we defined the link confinement $C$ as the expected
number of intersections if a link would follow a random trajectory, allowing
us to identify links whose trajectory is most affected by repulsive forces in
the network. Here, we characterize the properties of such confined links by
calculating the Kendall rank correlation $\tau$ between the link confinement
$C$ and other link properties for each data set. Specifically, we focus on the
(i) physical properties, complementary straightness $\bar{S}$ and link volume
$V_{\textrm{link}}$ (ii) abstract network properties link betweenness
$B_{\textnormal{link}}$ and link degree $k_{\textnormal{link}}(i,j)$, where
the latter is defined as the sum of the degrees of the endpoints link $(i,j)$.

Fig. 6 shows the correlation profiles of all 15 networks. A persistent pattern
we observe across the data sets is the positive correlation between link
confinement $C$ and link volume $V_{\textrm{link}}$. This is expected, since
larger links have more opportunities to intersect or be intersected by
neighboring links. On the other hand, correlations between link confinement
and straightness show a more curious pattern: we observe that lattice-like
networks tend to have positive and significant correlations $\tau$ between the
link confinement $C$ and link complementary straightness $\bar{S}$, indicating
that more winding links are also more confined. While for linked trees (fruit
fly neural networks) we find a negative $\tau$ between $C$ and $\bar{S}$. We
explain this by the fact that somata of neurons are represented as high-
volume, but short and straight links in fruit fly data sets (see discussion in
Sec. IV.2). This means they will be confined more due to their large volume
(positive correlations between $C$ and $V$), while also following a straighter
path (negative correlations between $C$ and $\bar{S}$).

For correlations between link confinement $C$ and the abstract network
properties, such as link betweenness $B_{\textnormal{link}}$ and degree
$k_{\textnormal{link}}$, we find consistent and significant positive
correlations for fruit fly neural networks. This indicates emergent
correlations between the three-dimensional layout and abstract network
properties of physical networks: more central links in the abstract network
tend to be more confined in physical space. For lattice-like networks and
trees, we find less consistent and weaker positive correlations. Overall, we
are able to show that the abstract network structure and physical layout are
intertwined for networks where we have sufficient information about the
surrounding environment of the physical links, such as the fruit fly neural
networks, which contain multiple neurons in close proximity with each other.
In contrast, in cases with limited information about the we identify positive
correlations for networks where

Figure 6: Link confinement correlations.We compute Kendall rank correlations
$\tau$ between link confinement $C$ and physical and abstract network
properties of links, obtaining a correlation profile for each physical
network. distinct correlation profiles for each abstract network category,
which are separated into three sections separated by red lines - top for
lattices, middle for trees and bottom for linked trees. Across many datasets
(excluding some trees), link volume $V_{\textrm{link}}$ and link confinement
$C$ has a statistically significant ($p<0.01$) small to moderate positive
relationship. In addition to this general trend, there are weaker trends in
terms of category-specific statistically significant correlations between link
confinement $C$ and other measures - $\bar{S}$ for lattices and $B_{link}$ for
trees, with a stronger trend for linked trees for all other measures -
$\bar{S}$, $k_{link}$ and $B_{link}$.

##  VI Discussion

Experimental data describing the three-dimensional shape of physical networks
is increasingly becoming available, and the growth in the number and size of
these data sets is expected to continue: connectome of the human brain cosists
of $\approx 10^{9}$ neurons and fungal mycorrhizal networks are estimated to
span $\approx 10^{17}$ km in Earth’s soil [14]. The new data calls for
extending the tool set of network science to analyze, model and understand how
the three-dimensional layout and physical interactions shape the structure and
function of physical networks. Here, we contributed to this effort in three
distinct ways: (i) We collected and standardized 15 data sets describing the
three-dimensional layout of physical network from diverse domains. (ii) We
characterized the structure of both the abstract network and physical layout
of the 15 systems using descriptors such as the degree distribution and
fractal dimension. (iii) We introduced link confinement as a method to
quantify how physical interactions shape link routing in physical networks,
allowing us to investigate emergent correlations between physical and abstract
network properties.

Our work may support future research on physical networks in several ways.
First, we promote the use of labeled skeleton graphs to represent both the
layout and the connectivity of physical networks. The skeleton captures the
shape of the network, while the labelling identifies the physical objects
corresponding the the nodes and links of the abstract network. Here, we
focused on treating junction points in the skeleton as nodes and sequences of
segments connecting them as links; however, the labeled skeletons are not
limited to such interpretation. For example, sub-graphs representing larger
functional units, such as neurons in the brain, can be identified as physical
nodes.

Second, our results also inform theoretical models of physical network growth.
Recent work that modeled physical nodes as spheres and links between them as
tubes [6, 20]. We found that most physical links follow close to straight
trajectories, suggesting that linear physical network models where links are
straight cylinders are indeed a useful class of models to understand
physicality in networks. On the other hand, these physical network models
generalize the classic Erdős-Rényi and Barabási-Albert models to physical
space and thus do not restrict the node degree. We, however, found that
junction points in real physical networks almost exclusively have degree
three, a fact that must be accounted for by future models. Note that to obtain
real physical networks with non-trivial degree distributions one must abandon
identifying junction points as nodes, instead we must identify larger sub-
graphs of the skeleton as physical nodes. In more formal terms, these networks
can be modeled as a network-of-networks: we represent each physical node as a
skeleton that has junction points with degree 3, and these physical sub-
networks are bound together to form a network-of-networks with no restriction
on the number of connections a sub-network can make with other sub-networks
[19]. Finally, we quantified the physical confinement of individual links by
comparing the path that links follow to randomized trajectories, allowing us
to identify correlations between physical and abstract network properties.
Future work may rely on other spatially randomized null-models to further
probe the relation of the two.

Our work is limited by the scope of the available data sets and computational
constraints. First, our data sets do not contain information about the
environment the networks are embedded in; therefore, we can only investigate
interactions between the components of the network and not interactions
between the networks and their surroundings. For example, we found the
strongest relationship between link confinement and abstract network structure
for the fruit fly neural network data sets and we found a weaker or no
relationship for individual neurons. This is likely due to the fact that the
fruit fly data sets contain multiple neurons, thus capture more of the
environment of individual physical links. Future work may consider more
complete data sets as they become available or theoretical models of network
growth could incorporate non-trivial environments. Also, investigating hard-
core exclusion in physical networks requires collision detection between
objects with complex shapes, which is a computationally difficult problem.
Future work may improve on these collision detection algorithms.

## References

  * [1] Hervé Abdi.  The kendall rank correlation coefficient.  Encyclopedia of Measurement and Statistics. Sage, Thousand Oaks, CA, pages 508–510, 2007.
  * [2] Oscar Kin-Chung Au, Chiew-Lan Tai, Hung-Kuo Chu, Daniel Cohen-Or, and Tong-Yee Lee.  Skeleton extraction by mesh contraction.  ACM transactions on graphics (TOG), 27(3):1–10, 2008.
  * [3] Péter Barthó, Andrea Slézia, Viktor Varga, Hajnalka Bokor, Didier Pinault, György Buzsáki, and László Acsády.  Cortical control of zona incerta.  Journal of Neuroscience, 27(7):1670–1681, 2007.
  * [4] Armin Bunde and Shlomo Havlin.  Fractals in science.  Springer, 2013.
  * [5] Paolo Crucitti, Vito Latora, and Sergio Porta.  Centrality measures in spatial networks of urban streets.  Physical Review E, 73(3):036125, 2006.
  * [6] Nima Dehmamy, Soodabeh Milanlouei, and Albert-László Barabási.  A structural transition in physical networks.  Nature, 563(7733):676–680, 2018.
  * [7] Paheli Desai-Chowdhry, Alexander B Brummer, and Van M Savage.  How axon and dendrite branching are guided by time, energy, and spatial constraints.  Scientific Reports, 12(1):20810, 2022.
  * [8] Marc Durand.  Architecture of optimal transport networks.  Physical Review E, 73(1):016116, 2006.
  * [9] Kayhan Foroutan-pour, Pierre Dutilleul, and Donald L Smith.  Advances in the implementation of the box-counting method of fractal dimension estimation.  Applied mathematics and computation, 105(2-3):195–210, 1999.
  * [10] Louis Gagnon, Sava Sakadžić, Frédéric Lesage, Joseph J Musacchia, Joël Lefebvre, Qianqian Fang, Meryem A Yücel, Karleyton C Evans, Emiri T Mandeville, Jülien Cohen-Adad, et al.  Quantifying the microvascular origin of bold-fmri from first principles with two-photon microscopy and an oxygen-sensitive nanoprobe.  Journal of Neuroscience, 35(8):3663–3675, 2015.
  * [11] Mohammad Retaful Islam, Gregory Tudryn, R Bucinell, L Schadler, and RC Picu.  Morphology and mechanics of fungal mycelium.  Scientific reports, 7(1):13070, 2017.
  * [12] Christof Koch and Allan Jones.  Big science, team science, and open science for neuroscience.  Neuron, 92(3):612–616, 2016.
  * [13] Michael LaBarbera.  Principles of design of fluid transport systems in zoology.  Science, 249(4972):992–1000, 1990.
  * [14] JR Leake and DJ Read.  Mycorrhizal symbioses and pedogenesis throughout earth’s history.  In Mycorrhizal mediation of soil, pages 9–33. Elsevier, 2017.
  * [15] Yanchen Liu, Nima Dehmamy, and Albert-László Barabási.  Isotopy and energy of physical networks.  Nature Physics, pages 1–7, 2020.
  * [16] Ron Milo, Shai Shen-Orr, Shalev Itzkovitz, Nadav Kashtan, Dmitri Chklovskii, and Uri Alon.  Network motifs: simple building blocks of complex networks.  Science, 298(5594):824–827, 2002.
  * [17] Mizue Ohashi, Hidetoshi Ikeno, Kotaro Sekihara, Toko Tanikawa, Masako Dannoura, Keitaro Yamase, Chikage Todo, Takahiro Tomita, and Yasuhiro Hirano.  Reconstruction of root systems in cryptomeria japonica using root point coordinates and diameters.  Planta, 249(2):445–455, 2019.
  * [18] Gérard Percheron.  Quantitative analysis of dendritic branching. i. simple formulae for the quantitative analysis of dendritic branching.  Neuroscience Letters, 14(2-3):287–293, 1979.
  * [19] Gábor Pete, Ádám Timár, Sigurdur Örn Stefánsson, Ivan Bonamassa, and Márton Pósfai.  A network-of-networks model for physical networks.  arXiv preprint arXiv:2306.01583, 2023.
  * [20] Márton Pósfai, Balázs Szegedy, Iva Bačić, Luka Blagojević, Miklós Abért, János Kertész, László Lovász, and Albert-László Barabási.  Understanding the impact of physicality on network structure.  arXiv preprint arXiv:2211.13265, 2022.
  * [21] Márton Pósfai, Balázs Szegedy, Iva Bačić, Luka Blagojević, Miklós Abért, János Kertész, László Lovász, and Albert-László Barabási.  Impact of physicality on network structure.  Nature Physics, pages 1–8, 2023.
  * [22] Punam K Saha, Gunilla Borgefors, and Gabriella Sanniti di Baja.  A survey on skeletonization algorithms and their applications.  Pattern recognition letters, 76:3–12, 2016.
  * [23] Carlos Santos-Filho, Camila M de Lima, César AR Fôro, Marcus A de Oliveira, Nara GM Magalhães, Cristovam Guerreiro-Diniz, Daniel G Diniz, Pedro F da C Vasconcelos, and Cristovam WP Diniz.  Visuospatial learning and memory in the cebus apella and microglial morphology in the molecular layer of the dentate gyrus and ca1 lacunosum molecular layer.  Journal of Chemical Neuroanatomy, 61:176–188, 2014.
  * [24] Johannes Schauer and Andreas Nüchter.  Collision detection between point clouds using an efficient kd tree implementation.  Advanced Engineering Informatics, 29(3):440–458, 2015.
  * [25] Louis K Scheffer, C Shan Xu, Michal Januszewski, Zhiyuan Lu, Shin-ya Takemura, Kenneth J Hayworth, Gary B Huang, Kazunori Shinomiya, Jeremy Maitlin-Shepard, Stuart Berg, et al.  A connectome and analysis of the adult drosophila central brain.  Elife, 9:e57443, 2020.
  * [26] Alexander Shapson-Coe, Michał Januszewski, Daniel R Berger, Art Pope, Yuelong Wu, Tim Blakely, Richard L Schalek, Peter H Li, Shuohong Wang, Jeremy Maitin-Shepard, et al.  A connectomic study of a petascale fragment of human cerebral cortex.  BioRxiv, pages 2021–05, 2021.
  * [27] Kaleem Siddiqi and Stephen Pizer.  Medial representations: mathematics, algorithms and applications, volume 37.  Springer Science & Business Media, 2008.
  * [28] Jacques Soddell and Robert Seviour.  Using box counting techniques for measuring shape of colonies of filamentous micro-organisms.  Complex systems: mechanism of adaptation. Edited by RJ Stonier and XH Yu. IOS Press, Amsterdam, pages 313–319, 1994.
  * [29] Andrea Tagliasacchi, Thomas Delame, Michela Spagnuolo, Nina Amenta, and Alexandru Telea.  3d skeletons: A state-of-the-art report.  In Computer Graphics Forum, volume 35, pages 573–597. Wiley Online Library, 2016.
  * [30] MM Torvund, TS Ma, VP Connaughton, F Ono, and RF Nelson.  Cone signals in monostratified and bistratified amacrine cells of adult zebrafish retina.  Journal of Comparative Neurology, 525(7):1532–1557, 2017.
  * [31] Matheus P Viana, Aidan I Brown, Irina A Mueller, Claire Goul, Elena F Koslover, and Susanne M Rafelski.  Mitochondrial fission and fusion dynamics generate efficient, robust, and evenly distributed network topologies in budding yeast cells.  Cell systems, 10(3):287–297, 2020.
  * [32] Tamás Vicsek.  Fractal growth phenomena.  World scientific, 1992.
  * [33] Brian Scott Whited.  Tangent-ball techniques for shape processing.  Georgia Institute of Technology, 2009.
  * [34] C Shan Xu, Michal Januszewski, Zhiyuan Lu, Shin-ya Takemura, Kenneth J Hayworth, Gary Huang, Kazunori Shinomiya, Jeremy Maitin-Shepard, David Ackerman, Stuart Berg, et al.  A connectome of the adult drosophila central brain.  BioRxiv, pages 2020–01, 2020.