capbtabboxtable[][]

# Numerical Modeling for Shoulder Injury Detection Using Microwave Imaging

Sahar Borzooei  Pierre-Henri Tournier  Victorita Dolean  Christian Pichot
Nadine Joachimowicz  Helene Roussel  Claire Migliaccio

###### Abstract

A portable imaging system for the on-site detection of shoulder injury is
necessary to identify its extent and avoid its development to severe
condition. Here, firstly a microwave tomography system is introduced using
state-of-the-art numerical modeling and parallel computing for imaging
different tissues in the shoulder. The results show that the proposed method
is capable of accurately detecting and localizing rotator cuff tears of
different size. In the next step, an efficient design in terms of computing
time and complexity is proposed to detect the variations in the injured model
with respect to the healthy model. The method is based on finite element
discretization and uses parallel preconditioners from the domain decomposition
method to accelerate computations. It is implemented using the open source
FreeFEM software.

Microwave imaging, Parallel computation, Inverse problem, Regularization,
Shoulder injury.

11footnotetext: Sahar Borzooei, Victorira Dolean, Christian Pichot and Claire
Migliaccio are with Université Côte d’Azur, Nice, France. (email:
sahar.borzooei@etu.univ-cotedazur.fr)22footnotetext: Pierre-Henri Tournier is
with Sorbonne Universitè, CNRS, Université Paris Cité, Inria, Laboratoire
Jacques-Louis Lions (LJLL), F-75005 Paris, France.33footnotetext: Hélène
Roussel is with Sorbonne Université, CNRS, Laboratoire de Génie Electrique et
Electronique de Paris, 75252, Paris, France.44footnotetext: Nadine
Joachimowicz is with Sorbonne Université, CNRS, Laboratoire de Génie
Electrique et Electronique de Paris, 75252, Paris, France, Université Paris
Cité, F-75006 Paris France.

##  1 Introduction

Rotator cuff tear (RCT) is one of the most common shoulder injuries. It
accounts for $70\%$ of shoulder pain and dysfunction in adults [1]. It can be
caused by a variety of factors, including age-related degeneration, overuse
and acute injury [2]. Rotator cuff is a group of muscles and their tendons
that work together to stabilize the shoulder, elevate and rotate the arm and
keep the head of the humerus securely in the shoulder socket [3]. Most of the
RCTs occur in the supraspinatus tendon [4]. A front view of the shoulder joint
and the tendon injury is shown in Figure 1 [5]. There are two types of RCTs:
partial and complete. A partial tear is when the tendon of the rotator cuff is
damaged but not completely severed. A complete tear is when the soft tissue is
detached from the bone. Complete tears can be categorized as small, medium or
large [6]. The severity of RCT does not correlate with the pain experience and
it can also be without symptoms, making it challenging to diagnose [7]. On the
other hand, early detection of a partial RCT can help in preventing its
development to full tear and may allow for nonsurgical treatment options [8].

Figure 1: Shoulder anatomy [5]

Medical imaging plays a critical role in defining the extent of RCTs, which
has important implications in clinical decision making and surgical planning
[9]. Magnetic resonance imaging (MRI), magnetic resonance angiography (MRA)
and ultrasound (US) are standard imaging modalities being used to assess the
presence and size of RCTs [10]. However, these diagnostics techniques are not
always accurate in depicting the size and number of involved tendons [11]. MRI
as the first-choice imaging modality for the detection of RCTs is costly,
bulky and is not suitable for on-site early detection. There is frequent
intra-operative reports claiming to find tears much larger than determined on
MRI, or even the lack of a tear [12]. MRA is an invasive procedure that
requires the injection of a contrast agent [13]. US is operator-dependent and
may be limited by the patients body habitus [14]. As a result, an alternative
low-cost, portable and non-invasive method is demanding, specially on-site for
competitive athletes such as swimmers [15]. Electromagnetic imaging (EMI)
systems have shown promising results in various medical applications [16, 17,
18, 19]. Fully circular tomographic-based EMI systems are designed for
different applications, such as brain [20, 21], breast [22] and knee imaging
[23]. Circular data acquisition effectively helps in increasing the cross-
range resolution of the reconstructed images [24].
To the best of our knowledge, there is no EMI system to detect shoulder
injuries. Designing an EMI system for the shoulder is challenging due to
following reasons. Firstly, the complex anatomy of the shoulder prevents
designing a fully circular phased array for spherical scan and full data
acquisition. As a consequence, the antenna array has to be defined conformal
to the shoulder geometry. Another challenge would be the electrically large
size of the shoulder along with the heterogeneous nature of the tissues,
characterized by high losses (also a characteristic of the matching medium),
making it difficult to achieve a high resolution three-dimensional (3D)
reconstruction. Shoulder being located in the near-field of the antenna array,
an efficient 3D EM-modelling is required to consider coupling effects and
near-field interactions between the imaging system and the shoulder. Besides,
the variability in the shoulder anatomy among individuals prevents us to
consider a particular shoulder structure as a priori knowledge. Finally, the
accurate knowledge of the dielectric properties of the shoulder tissues, and
more specifically the synovial fluid, remains a challenge. At the early stage
design, when assessing the potential of EMI on a new application, using
numerical modeling to help design the system presents several advantages: it
allows to accurately model the human body complexity and electrical
properties, and to simulate various anatomical scenarios in a flexible way,
while reducing cost and time. In this paper, we propose a feasibility study
based on numerical assessments to design a shoulder injury detection system
and demonstrate its performance. We make use of an anthropomorphic numerical
model of the shoulder and a built-in EM modeling based on the open source code
FreeFEM, which offers flexibility in dealing with complex geometries and
electrically large structures. The paper is organized as follows: Section 2
discusses the dielectric properties of the shoulder tissues, while Section 3
describes the choices for a first shoulder imaging system with a dense array
of antennas to achieve high resolution. The mathematical framework of the EM
modeling is outlined in Section 4, and numerical solutions are presented in
Section 5. Section 6 investigates the reduction of the number of antennas in
the EMI system.

##  2 Dielectric properties of the shoulder tissues

Once the RTCs occurs, synovial fluid (SF) accumulates at the location of the
tear [25, 26] that leads to a change in the dielectric properties in the
shoulder joint [27]. It is reported that the volume of SF on aspiration prior
to arthroscopic rotator cuff repair correlates with tear size. The mean
aspirate volume of partial tears is $0.76\pm 0.43\leavevmode\nobreak\ mL$,
small tears $1.46\pm 1.88\leavevmode\nobreak\ mL\leavevmode\nobreak\ $, medium
tears $3.04\leavevmode\nobreak\ \pm 2.21mL\leavevmode\nobreak\ $, and large
tears $6.60\leavevmode\nobreak\ \pm 3.23\leavevmode\nobreak\ mL$ [26]. When
there is a difference in the dielectric properties between different shoulder
tissues and the SF, microwave imaging can detect this contrast. The larger the
difference, the easier the detection. Thus, it is essential to measure the
dielectric properties of the SF in relation to the other tissues. In a recent
work with the purpose of detecting knee injury [23], high values of SF are
considered compared to the rest of the tissues; however, to the best of our
knowledge, no measurement of dielectric properties of SF is published in this
or any other literature. It is reported that low frequencies ranging between
$0.5-3\leavevmode\nobreak\ GHz$ are suitable to achieve deeper tissue
penetration as well as acceptable imaging quality [28]. In this section, the
dielectric properties of SF, muscle, bone cortical and tendons are measured in
the frequency range of $1-4\leavevmode\nobreak\ GHz$ that almost include and
extend the recommended frequency range. In the case of SF, four patients with
knee or shoulder pain were sampled to obtain the real human SF. We computed
the average value of the measured dielectric properties of the four SF
samples, which had different temperature and pathology. The mixtures of the
other tissues were made based on the SUPELEC RECIPES website [29]. An
optimization code based on Kraszweski’s binary law gives the concentration of
TritonX-100 and salt required to produce mixtures, whose dielectric properties
are close to those given by a Cole-Cole model, for each biological tissue of
the shoulder region [30]. Measurements were performed using a homemade coaxial
probe connected to a ZVH8 Vector Network Analyzer [31]. Fig. 2 gathers all
measurements and comparisons with values from reference websites
111https://itis.swiss/database,222http://www.ifac.cnr.it. Our approach is
validated by the good agreement with the references.

Figure 2: Dielectric properties of the tissues

The relative contrast between SF and other tissues over the frequency range of
interest is around 30%, except for the muscle for which there is a 20%
difference in the real part, making the tear detection feasible but
challenging. The anthropomorphic shoulder model is shown in Figure 3. The
complex permittivity values of the shoulder tissues at $1\leavevmode\nobreak\
GHz$ are given in Table 1 and are used in the simulations. In this Table, the
SF value corresponds to our measurements, but the values of other tissues have
been already available and are taken from reference websites${}^{1},^{2}$. The
complex permittivity of the matching medium is chosen equal to that of the
muscle.

Region | Complex Permittivity | Color in Fig. 3
---|---|---
Bone cortical | $12.4-2.79j$ | Blue
Tendon | $45.6-13.66j$ | Yellow
Muscle | $54.8-17.43j$ | Transparent
Skin | $40.9-16.17j$ | Green
SF | $68.0-29.0j$ | Red

Table 1: Complex dielectric properties at $1\leavevmode\nobreak\ GHz$ Figure
3: Shoulder Phantoms, real permittivity values

##  3 Electromagnetic imaging System

We first design an EMI system with a dense array of antennas that illuminates
the shoulder from different angles. This multi-perspective approach helps in
reconstructing a comprehensive and accurate representation of the internal
structures [32]. Results will be used as reference for the optimized system
described in section 6. The EMI system consists of 96 ceramic
($\varepsilon_{r}=59$) loaded open-ended waveguides arranged on $2$ metallic
fully-circular and $2$ metallic half-circular layers. The two sides of the
imaging chamber are open. This wearable imaging system with open sides is
adapted to the real shoulder structure and designed to surround it partially,
as shown in Figure 4. The width and height of the rectangular waveguides are
$2.1\leavevmode\nobreak\ cm$ and $0.075\leavevmode\nobreak\ cm$, respectively.
Their frequency bandwidth is $0.93-1.85\leavevmode\nobreak\ GHz$. Here, the
operating frequency of $1GHz$ is chosen as a good compromise regarding
penetration depth and low specific absorption rate (SAR) [33].

Figure 4: EMI system for the shoulder

Considering the matching medium as the reference permittivity, the wavelength
in this medium is $\lambda=4.2\leavevmode\nobreak\ cm$. The diameter of the
chamber is $7.14\lambda$, the larger length is $3\lambda$ and the distance
between antenna layers is $0.47\lambda$. The larger size of the modeled large
and partial tears is $0.69\lambda$ and $0.397\lambda$, respectively. Note that
we can expect a better resolution than $0.25\lambda$, because EMI system is
operating in the near field [34]. The space between antennas and the shoulder
is filled with a matching medium to overcome air-skin reflections [35].

##  4 Mathematical Framework

###  4.1 Finite Element Mesh Generation

Finite element 3D mesh generation of the complete system is a challenging step
due to both the complex geometry of the real body from Computer-Aided Design
(CAD) models and the imaging system components. We have chosen to use
realistic surface CAD models for shoulder profile and bones including humerus
and scapula from a library of 3D models related to the anatomy [36]. A simple
model of rotator cuff tendons was then built surrounding the shoulder joint.
The skin is considered with a thickness of $2\leavevmode\nobreak\ mm$
surrounding the muscle geometry. The injury is modeled as an ellipsoid in the
approximate location of the tear in the supraspinatus tendon, with two
different size configurations: small tear with volume of $1mL$, and large tear
with volume of $5mL$. Note that in the healthy case, the ellipsoid is replaced
with muscle. The regions corresponding to the different tissues are visible in
Figure 3. The remaining area (excluding bone, injury, tendon, skin)
corresponds to muscle. The conformal 3D mesh of the complete system must
satisfy non-overlapping elements and nodal-match conditions. At the interfaces
between domains, special attention should be given to ensure that the mesh is
well-aligned, conformal and continuous [37]. Taking into account that the
surface meshes of some tissues intersect the imaging system boundaries, the
mesh generation procedure is done by extracting and rebuilding boundary curves
corresponding to each intersection for all intersecting surfaces, see Figure
5. In this work, we use the open-source finite element software FreeFEM, which
allows to generate 3D meshes and solve partial differential equations (PDEs)
with the finite element method (FEM) [38]. FreeFEM is interfaced with the MMG
remeshing library [39, 40], which makes it possible to generate adapted
tetrahedral meshes of complex geometric configurations, such as the one
considered in this paper.

Figure 5: The intersecting surfaces of EMI system with shoulder

###  4.2 Forward Modeling

The 3D domain ($\Omega$) includes the imaging system and the shoulder as a
heterogeneous dissipative non-magnetic medium of complex permittivity
$\varepsilon_{r}=(\varepsilon_{r}^{{}^{\prime}}-\frac{\sigma
j}{\omega\varepsilon_{0}})$, with $\sigma$ the conductivity and
$\varepsilon_{r}^{{}^{\prime}}$ the electrical permittivity of each tissue,
$\varepsilon_{0}$ the electrical permittivity of free space, and $\omega$ the
angular frequency. In the frequency domain, the electric field
$\xi(\mathbf{x},\mathbf{t})=\Re(\mathbf{E(x)}e^{\mathrm{i}\omega\mathbf{t}})$
has harmonic dependence on time of angular frequency $\omega$, where
$\mathbf{E(x)}$ is its complex amplitude depending on the space variable
$\mathbf{x}$. The boundary value problem is defined in equation 1.

$\begin{cases}\nabla\times(\nabla\times\mathbf{E})-k^{2}\mathbf{E}=0,&\text{in
}\Omega,\\\ \mathbf{E}\times\textbf{n}=0,&\text{on }\Gamma_{m}\\\
\nabla\times\mathbf{E}\times\textbf{n}+\mathrm{i}\beta\textbf{n}\times(\mathbf{E}\times\textbf{n})=g&\text{on
}\Gamma_{t}\\\
\nabla\times\mathbf{E}\times\textbf{n}+\mathrm{i}\beta\textbf{n}\times(\mathbf{E}\times\textbf{n})=0&\text{on
}\Gamma_{r}\\\
\nabla\times\mathbf{E}\times\textbf{n}+\mathrm{i}k\textbf{n}\times(\mathbf{E}\times\textbf{n})=0&\text{on
}\Gamma_{o},\end{cases}$ (1)

where $\mathbf{E(x)}$ is the solution of equation (1) for each transmitting
antenna and $k=\omega\sqrt{\varepsilon_{r}\varepsilon_{0}\mu_{0}}$ is the
complex wavenumber of the inhomogeneous medium at each point of the 3D space,
with $\mu_{0}$ the permeability of free space. $\beta$ is the propagation
wavenumber along the waveguide which corresponds to the propagation of the
dominant transverse electric mode $TE_{10}$. The excitation term is defined as
$g=2\mathrm{i}\beta\mathbf{E}^{TE_{10}}$, imposing an incident wave
corresponding to the excitation of the dominant mode of the transmitting
waveguide. $\Gamma_{t}$ is the port of the transmitting waveguide,
$\Gamma_{r}$ corresponds to ports of receiving waveguides, $\Gamma_{m}$
represents the metallic surfaces of the waveguides and the walls between them,
and $\Gamma_{o}$ represent the three open sides (right, left and bottom) of
the chamber and the boundaries of shoulder profile. Let us consider
$V=\\{\mathbf{v}\leavevmode\nobreak\ \in\leavevmode\nobreak\
H(\mbox{curl},\Omega),\leavevmode\nobreak\
\mathbf{v}\times\mathbf{n}=0\leavevmode\nobreak\ \text{on}\leavevmode\nobreak\
\Gamma_{m}\\}$, where $H(\mbox{curl},\Omega)$ is the space of square
integrable functions whose curl is also square integrable, and $v$ is the test
function. The variational formulation of equation (1) is: find
$\mathbf{E}\leavevmode\nobreak\ \in\leavevmode\nobreak\ V$ such that:

$\displaystyle\int_{\Omega}[(\nabla\times\mathbf{E})\cdot(\nabla\times\mathbf{v})-k^{2}\mathbf{E}\cdot\mathbf{v}]+\leavevmode\nobreak$
(2)
$\displaystyle\int_{\Gamma_{t}\cup\Gamma_{r}}\mathrm{i}\beta(\mathbf{E}\times\mathbf{n})\cdot(\mathbf{v}\times\mathbf{n})+\leavevmode\nobreak$
$\displaystyle\int_{\Gamma_{o}}\mathrm{i}k(\mathbf{E}\times\mathbf{n})\cdot(\mathbf{v}\times\mathbf{n})$
$\displaystyle\leavevmode\nobreak\
=\int_{\Gamma_{t}}g\cdot\mathbf{v},\leavevmode\nobreak\
\forall\leavevmode\nobreak\ \mathbf{v}\leavevmode\nobreak\
\epsilon\leavevmode\nobreak\ V.$

Through the successive solution of variational problem (2) for each
transmitting antenna, we can compute the scattering matrix, a set of complex-
valued reflection and transmission coefficients, given in equation (3):

$S_{ij}=\frac{\int_{\Gamma_{r}}\mathbf{E}\cdot\mathbf{E}^{TE_{10}}}{\int_{\Gamma_{r}}|{\mathbf{E}^{TE_{10}}}|^{2}}.$
(3)

The computed scattering parameters form a complex matrix of size $96\times
96$, which is used to produce synthetic data by adding noise to the real and
imaginary parts of the coefficients, using a multiplicative white Gaussian
noise. The Gaussian noise is applied independently to the real and imaginary
parts of each $S_{ij}$ coefficient as independent random variables drawn from
the standard normal distribution $N(0,1)$. In this work, we have corrupted the
data $S_{ij}$ with $23\leavevmode\nobreak\ dB$ noise.

###  4.3 Domain decomposition

The finite element discretization of our problem leads to a large ill-
conditioned linear system $A\mathbf{u}=\mathbf{b}$ for each transmitting
antenna. Domain decomposition methods (DDMs) are efficient tools to solve such
large systems in parallel, both in terms of convergence and computing time
[39]. These approaches rely on the partition of the computational domain into
smaller subdomains, leading to subproblems of smaller size which are
manageable by direct solvers [40]. A Krylov iterative solver (GMRES) along
with an Optimized Restricted Additive Schwarz (ORAS) preconditioner is chosen
to solve our problem. The domain decomposition preconditioner is implemented
in the HPDDM library [41], an open source high-performance unified framework
for domain decomposition methods which is interfaced with the FreeFEM
software.

###  4.4 Inverse Problem

Let $\kappa=k^{2}$ be the unknown complex parameter of the inverse problem in
each point of $\Omega$. In this step an optimization problem, including a fit-
to-data term and a regularizing term, is defined with the following cost
function:

$J{(\kappa)}=\frac{1}{2}\sum_{j=1}^{N}\sum_{i=1}^{N}\frac{|S_{ij}(\kappa)-S_{ij}^{syn}|^{2}}{|S_{ij}^{empty}|^{2}}+\alpha
R(\kappa).$ (4)

Where

  * •

$S_{ij}^{syn}$ are the scattering coefficients obtained from the forward
problem and are referred to as synthetic data in the rest of the paper.

  * •

$S_{ij}^{\kappa}$ are the scattering coefficients computed for the unknown
$\kappa$ at the current iteration.

  * •

$S_{ij}^{empty}$ are the coefficients computed from the simulation when the
domain is filled only with the homogeneous matching medium, used for
normalization.

  * •

$N$ is the number of antennas of the system.

  * •

the Tikhonov regularization method is chosen to reduce the noise effect,
defined as $R(\kappa)=\frac{1}{2}\int_{\Omega}{|\nabla{\kappa}|^{2}}$.

  * •

the regularization parameter $\alpha$ is chosen empirically equal to $10^{-6}$
to reach a good compromise between denoising as well as achieving suitable
image quality with respect to certain properties such as smoothness or
sparsity [42].

Minimizing functional (4) with respect to the parameter $\kappa$ is done by
computing its differential and using the adjoint approach in order to simplify
its expression [43]. The gradient of the resulted adjoint functional is used
in the limited-memory Broyden-Fletcher-Goldfarb-Shanno (L-BFGS) optimization
algorithm. We refer to [20] for a detailed description of the inverse modeling
that we have followed in this work. Note that we avoid inverse crime by adding
noise to the synthetic data (explained in Section 4.2) and not using a priori
knowledge of the body structure. Eliminating a priori knowledge is done by
using a mesh that includes the geometry and structure of the body for
generating synthetic data (Figure 6, left) but defining a different
homogeneous mesh that is limited to the imaging chamber (Figure 6, right) for
the inverse problem. Using Nedelec edge finite elements (FE) of first
polynomial degree ($r=1$) for the FE discretization results in 510531 unknowns
for the direct mesh and $357535$ unknowns for the inverse mesh. The spatial
resolution for both generated meshes can be defined as the number of points
per wavelength $n_{\lambda}$. Here we choose $n_{\lambda}=\frac{\lambda}{6}$,
where $\lambda$ corresponds to the wave propagation in the matching
medium/muscle tissue.

Figure 6: Left: Synthetic data, Right: Inverse problem.

##  5 Numerical Results

Results are obtained on the Université Côte d’Azur’s High-Performance
Computing (HPC) center. In this HPC center, cluster is composed of 48 CPU
computing nodes, including 32 nodes with Dual Intel Xeon Gold processor,
providing 40 cores per node and 192 GB of memory and 16 nodes with 2 AMD Epyc
processors, providing 32 cores per node and 256 GB of memory. The simulations
presented in this paper were carried out using $480$ cores. Each
reconstruction starts from an initial guess of homogeneous matching medium.
The reconstruction results shown here are obtained after 60 iterations; the
residual is decreased by a factor of $10^{-2}$. Subsequent iterations do not
provide any further noteworthy decrease. Figures 7 and 8 show the 3D view of
the exact (top) and reconstructed (bottom) results for healthy, partially and
fully injured shoulder on the regions of interest. Note that in the healthy
shoulder case, the ellipsoid is filled in with muscle as shown in Figs 7, 8
for comparison purpose. The rest of the muscle is not shown for the sake of
visibility. The tear is visible in reconstructed real and imaginary parts. For
a better view, we also plot the difference between reconstructed permittivity
of healthy case and those for the partially and largely injured cases, known
as differential images [44]. They are shown for the real part in Figure 9,
proving that the inverse algorithm succeeds in detecting the injury, with
respect to size and location, even for the partial tear. Results for the
imaginary part are similar, but are not presented here for the sake of
brevity.

Figure 7: Cross-section of the real part of the permittivity. Top: exact
permittivity, bottom: reconstructed permittivity. From left to right: Healthy,
partial tear and large tear.

Figure 8: Cross-section of the imaginary part of the permittivity. Top: exact
permittivity, bottom: reconstructed permittivity. From left to right: Healthy,
partial tear and large tear. Figure 9: The differential image for partial and
large tear.

$err_{absolute}=|{\varepsilon_{r}^{reconstructed}-\varepsilon_{r}^{exact}}|$
(5)

Figure 10: The absolute error on real part, from left to right: Healthy
shoulder, partial and large tear Figure 11: The absolute error on imaginary
part, from left to right: Healthy shoulder, partial and large tear

The reconstruction is further assessed by computing the absolute error on the
real and imaginary parts according to equation (5) applied on each pixel of
the reconstructed domain. Results are shown in Figure 10 and 11. Error
distribution is not uniform across tissues. Highest values are on the edges of
the bone area, due to the higher contrast between the complex permittivity of
bone and muscle. However, inside the bone in the areas that are not in the
vicinity of the edges, the error decreases as expected, due to better
performance of the inversion algorithm in a homogeneous region compared to the
interfaces. The minimum values of the absolute error in the injury area are
$0.78$ and $0.435$, respectively for the real and imaginary parts of the
partial tear. These values for the large tear are $0.0088$ and $0.065$,
respectively.

##  6 An improved design of EMI system

In this section, an efficient design of the imaging system in terms of reduced
number of antennas is proposed. With less transmitting and receiving
waveguides, faster computing time as well as cost-effective design can be
achieved. The system is optimized according to the following guidelines:

  * •

The number of antennas is a power of two to ease practically feeding signals
into a switching matrix.

  * •

We introduce a spatial shift between antennas of adjacents rings as shown in
Fig. 12, in order to preserve the spatial diversity.

Different design configurations are assessed based on the value of the $L^{2}$
norm error, referred to as $err$ according to equation (6), as well as the
mean value of the subtraction between the reconstructed injured model and
reconstructed healthy model in the injury area, defined in equation (7). Three
different configurations are considered, with 16, 32, and 64 antennas ; the
reference configuration corresponds to the complete system with $96$ antennas.
The values of $err$ for each configuration are computed in the whole domain as
well as in the injury area for the partially injured model and are gathered in
Table 2. Looking at the $err$ values on both real and imaginary parts, there
is no significant increase when reducing the number of antennas by 64
($N=32$). However, further removing $16$ antennas and going from $N=32$ to
$N=16$ antennas leads to a higher increase in $err$. The trend is confirmed by
looking at the contrast values $ctr$ which are shown in Table 3 for both
partial and large tear and with different levels of noise. In the partial tear
case, the value of $ctr$ remains positive notable for $N=32$ for all noise
levels, and we expect to be able to detect the injury. In contrast, $ctr$ is
either negligible or negative for $N=16$, resulting in poor differential
images. A first conclusion is that the best design compromise appears to be
the $N=32$ configuration, which still allows detection in the most difficult
partial tear case. This was further confirmed by extracting and comparing the
differential images for each configuration, but due to space limitation the
images are not shown in the paper. The optimal $N=32$ antennas configuration
is depicted in Figure 12 and has $11$ and $10$ antennas in both full layers
and $6$ and $5$ antennas on both half layers, from left to right. Note that
the computing time decreases to about half of the required time for the full
$N=96$ configuration.

Figure 12: Optimal shoulder EMI system

$err=\frac{||\Re(\varepsilon_{r}^{reconstructed})-\Re(\varepsilon_{r}^{exact})||_{2}}{||\Re(\varepsilon_{r}^{exact})||_{2}}$
(6)

$ctr=\varepsilon_{r}^{injured}-\varepsilon_{r}^{healthy}$ (7)

Table 2: $err$ in the injury area and total domain for the partial tear case,
and computing time.

N | Total | Injury | time
---|---|---|---
real | imaginary | real | imaginary | (min:s)
96 | $8.8\%$ | $21.4\%$ | $14\%$ | $24.4\%$ | $20:09$
64 | $9.0\%$ | $21.2\%$ | $10.4\%$ | $19.2\%$ | $16:24$
32 | $9.4\%$ | $20.6\%$ | $16.8\%$ | $33.3\%$ | $12:02$
16 | $10.8\%$ | $17.2\%$ | $19.3\%$ | $45.9\%$ | $11:27$

Table 3: $ctr$ in the injury area for partial and large tears with different
noise levels.

| N | Noise of 23 dB | Noise of 15 dB | Noise of 10 dB
---|---|---|---|---
| real | imaginary | real | imaginary | real | imaginary
Partial | 96 | $1.803$ | $1.910$ | $1.934$ | $1.767$ | $1.973$ | $1.696$
64 | $1.605$ | $1.673$ | $1.664$ | $1.654$ | $2.233$ | $1.503$
32 | $0.890$ | $0.690$ | $0.843$ | $0.669$ | $0.599$ | $0.681$
16 | $0.133$ | $0.1918$ | $-0.136$ | $0.0193$ | $0.3617$ | $0.256$
Large | $96$ | $5.254$ | $5.101$ | $5.182$ | $5.230$ | $5.412$ | $5.033$
$64$ | $5.053$ | $4.312$ | $5.145$ | $4.310$ | $5.396$ | $4.244$
$32$ | $2.322$ | $2.241$ | $2.467$ | $2.191$ | $2.524$ | $2.204$
$16$ | $0.750$ | $0.711$ | $0.692$ | $0.681$ | $0.787$ | $0.780$

##  7 Conclusion

We have conducted the first numerical study on the feasibility to image the
shoulder joint with an EMI system. The first part of the study consisted in
the investigation of the dielectric properties of the shoulder tissues, in
particular the synovial fluid which was measured for the first time. Secondly,
we have shown that it is possible to reconstruct 3D images of the shoulder
joint based on an anthropomorphic numerical model of the shoulder with an
imaging system composed of a dense array of $96$ antennas. The reconstruction
takes $20$ minutes and $9$ seconds on $480$ computing cores and shows great
promise for rapid diagnosis or medical monitoring. Finally, we take advantage
of numerical modeling to optimize the number of antennas in the EMI system. We
achieve a drastic reduction from $96$ to $32$ antennas. After having proved
the relevance of microwave imaging for the detection of shoulder injury, the
next step will be to manufacture the imaging system and conduct measurements
on built-in phantoms of the different tissues in order to further validate the
imaging system proposed in this work.

## Acknowledgment

This project has received funding from the European Union’s Horizon 2020
research and innovation programme under the Marie Skłodowska-Curie grant
agreement No 847581 and is co-funded by the Région Provence-Alpes-Côte d’Azur
and IDEX $UCA^{JEDI}$and supported by National Research Agency (ANR) under
reference number ANR-15-IDEX-01. The authors are grateful to the OPAL
infrastructure from Université Côte d’Azur and the Université Côte d’Azur’s
Center for High-Performance Computing for providing resources and support. The
authors like to thank Dr Eric GIBERT, rheumatologist at the Pitié-Salpêtrière
Hospital in Paris, for providing the SF samples for dielectric
characterization.

## References

  * [1] Barbara Juliette Mera  “Current perspectives on rotator cuff disease”  In _Osteology_ 2.2 MDPI, 2022, pp. 62–69
  * [2] Richard W Nyffeler, Nicholas Schenk and Philipp Bissig  “Can a simple fall cause a rotator cuff tear? Literature review and biomechanical considerations”  In _International orthopaedics_ 45 Springer, 2021, pp. 1573–1582
  * [3] Jeffrey R Dugas et al.  “Anatomy and dimensions of rotator cuff insertions”  In _Journal of shoulder and elbow surgery_ 11.5 Elsevier, 2002, pp. 498–503
  * [4] Reza Omid and Brian Lee  “Tendon transfers for irreparable rotator cuff tears”  In _JAAOS-Journal of the American Academy of Orthopaedic Surgeons_ 21.8 LWW, 2013, pp. 492–501
  * [5] Dr. Fernando Cinagava  “https://consultoriodoombro.com.br/por-que-meus-ombros-doem-com-video-explicativo-sobre-tendinite-do-ombro/”
  * [6] Umile Giuseppe Longo et al.  “Histopathology of rotator cuff tears”  In _Sports Medicine and Arthroscopy Review_ 19.3 LWW, 2011, pp. 227–236
  * [7] Atsushi Yamamoto et al.  “The impact of faulty posture on rotator cuff tears with and without symptoms”  In _Journal of shoulder and elbow surgery_ 24.3 Elsevier, 2015, pp. 446–452
  * [8] Ali Abdelwahab et al.  “Traumatic rotator cuff tears-Current concepts in diagnosis and management”  In _Journal of Clinical Orthopaedics and Trauma_ 18 Elsevier, 2021, pp. 51–55
  * [9] Gabrielle P Konin  “Imaging of the Rotator Cuff”  In _Massive Rotator Cuff Tears: Diagnosis and Management_ Springer, 2014, pp. 31–51
  * [10] Joseph P Iannotti et al.  “Accuracy of office-based ultrasonography of the shoulder for the diagnosis of rotator cuff tears”  In _JBJS_ 87.6 LWW, 2005, pp. 1305–1311
  * [11] Sharlene A Teefey et al.  “Detection and quantification of rotator cuff tears: comparison of ultrasonographic, magnetic resonance imaging, and arthroscopic findings in seventy-one consecutive cases”  In _JBJS_ 86.4 LWW, 2004, pp. 708–716
  * [12] Xusheng Li et al.  “The accuracy of MRI in diagnosing and classifying acute traumatic multiple ligament knee injuries”  In _BMC Musculoskeletal Disorders_ 23 Springer, 2022, pp. 1–7
  * [13] Ashley Yip et al.  “Magnetic resonance imaging compared to ultrasonography in giant cell arteritis: a cross-sectional study”  In _Arthritis research & therapy_ 22.1 BioMed Central, 2020, pp. 1–8
  * [14] Mark Adrian  “MRI: understanding its limitations”  In _Br Columb Med J_ 47, 2005, pp. 359–361
  * [15] JT Stahnke, M Schenk and KP Sharpe  “Magnetic Resonance Imaging of Symptomatic Swimmer’s Shoulder”  In _Int J Radiol Imaging Technol_ 6, 2020, pp. 070
  * [16] Kazuki Kanazawa, Kazuki Noritake, Yuriko Takaishi and Shouhei Kidera  “Microwave imaging algorithm based on waveform reconstruction for microwave ablation treatment”  In _IEEE Transactions on Antennas and Propagation_ 68.7 IEEE, 2020, pp. 5613–5625
  * [17] Daniel Oloumi et al.  “Microwave imaging of breast tumor using time-domain UWB circular-SAR technique”  In _IEEE Transactions on Medical Imaging_ 39.4 IEEE, 2019, pp. 934–943
  * [18] Amir Mirbeik-Sabzevari et al.  “Synthetic ultra-high-resolution millimeter-wave imaging for skin cancer detection”  In _IEEE Transactions on Biomedical Engineering_ 66.1 IEEE, 2018, pp. 61–71
  * [19] Kamel S Sultan, Beadaa Mohammed, Mohamed Manoufali and Amin M Abbosh  “Portable electromagnetic knee imaging system”  In _IEEE Transactions on Antennas and Propagation_ 69.10 IEEE, 2021, pp. 6824–6837
  * [20] Pierre-Henri Tournier et al.  “Numerical modeling and high-speed parallel computing: New perspectives on tomographic microwave imaging for brain stroke detection and monitoring”  In _IEEE Antennas and Propagation Magazine_ 59.5 IEEE, 2017, pp. 98–110
  * [21] Markus Hopfer et al.  “Electromagnetic Tomography for Detection, Differentiation, and Monitoring of Brain Stroke: A Virtual Data and Human Head Phantom Study.”  In _IEEE Antennas and Propagation Magazine_ 59.5 IEEE, 2017, pp. 86–97
  * [22] Caroline Zellweger et al.  “Breast computed tomography: diagnostic performance of the maximum intensity projection reformations as a stand-alone method for the detection and characterization of breast findings”  In _Investigative Radiology_ 57.4 Wolters Kluwer, 2022, pp. 205–211
  * [23] Kamel Sultan  “Design and implementation of electromagnetic knee imaging systems”, 2022
  * [24] Lulu Wang  “Microwave Imaging and Sensing Techniques for Breast Cancer Detection”  In _Micromachines_ 14.7 MDPI, 2023, pp. 1462
  * [25] Athanasios Papatheodorou et al.  “US of the shoulder: rotator cuff and non–rotator cuff disorders”  In _Radiographics_ 26.1 Radiological Society of North America, 2006, pp. e23–e23
  * [26] Michael Stone et al.  “Synovial fluid volume at the time of arthroscopic rotator cuff repair correlates with tear size”  In _Cureus_ 12.7 Cureus, 2020
  * [27] John Kiel and Eamon Olwell  “It might be a tumor: a unique presentation of a chronic rotator cuff tear”  In _African Journal of Emergency Medicine_ 10.4 Elsevier, 2020, pp. 288–290
  * [28] Natalia K Nikolova  “Microwave imaging for breast cancer”  In _IEEE microwave magazine_ 12.7 IEEE, 2011, pp. 78–94
  * [29] “SUPELEC RECIPES, Standard Uwb Phantom ELECtromagnetic Recipes [online]. Available http://applis.iut.univ-paris-diderot.fr/CNRS/.”
  * [30] Sami Gabriel, RW Lau and Camelia Gabriel  “The dielectric properties of biological tissues: III. Parametric models for the dielectric spectrum of tissues”  In _Physics in medicine & biology_ 41.11 IOP Publishing, 1996, pp. 2271
  * [31] Soroush Abedi, Hélène Roussel and Nadine Joachimowicz  “Standardized Phantoms”  In _Electromagnetic Imaging for a Novel Generation of Medical Devices: Fundamental Issues, Methodological Challenges and Practical Implementation_ Springer, 2023, pp. 1–32
  * [32] Chang Liu et al.  “Multi-perspective ultrasound imaging technology of the breast with cylindrical motion of linear arrays”  In _Applied Sciences_ 9.3 MDPI, 2019, pp. 419
  * [33] Kamel S Sultan, Ahmed Mahmoud and Amin M Abbosh  “Textile electromagnetic brace for knee imaging”  In _IEEE Transactions on Biomedical Circuits and Systems_ 15.3 IEEE, 2021, pp. 522–536
  * [34] Tie Jun Cui, Weng Cho Chew, Xiao Xing Yin and Wei Hong  “Study of resolution and super resolution in electromagnetic imaging for half-space problems”  In _IEEE Transactions on Antennas and Propagation_ 52.6 IEEE, 2004, pp. 1398–1411
  * [35] Hoi-Shun Lui, Andreas Fhager and Mikael Persson  “Matching medium for biomedical microwave imaging”  In _2015 International Symposium on Antennas and Propagation (ISAP)_ , 2015, pp. 1–3  IEEE
  * [36] “Plasticboy, Blender Rigged V07 Anatomy 3D models store [online]. Available https://www.plasticboy.co.uk/store/index.html”
  * [37] Richard Michael Jack Kramer and David R Noble  “Guaranteed Quality Conformal Mesh Decomposition.”, 2015
  * [38] Frédéric Hecht  “New development in FreeFem++”  In _Journal of numerical mathematics_ 20.3-4 De Gruyter, 2012, pp. 251–266
  * [39] Niall Bootland, Sahar Borzooei, Victorita Dolean and Pierre-Henri Tournier  “Numerical assessment of PML transmission conditions in a domain decomposition method for the Helmholtz equation”  In _arXiv preprint arXiv:2211.06859_ , 2022
  * [40] Sahar Borzooei, Victorita Dolean, Pierre-Henri Tournier and Claire Migliaccio  “Solution of time-harmonic Maxwell’s equations by a domain decomposition method based on PML transmission conditions”  In _arXiv preprint arXiv:2211.04912_ , 2022
  * [41] Pierre Jolivet, Frédéric Hecht, Frédéric Nataf and Christophe Prud’Homme  “Scalable domain decomposition preconditioners for heterogeneous elliptic problems”  In _Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis_ , 2013, pp. 1–11
  * [42] Yuhui Zheng, Min Li, Jianwei Zhang and Jin Wang  “Selection of regularization parameter in GMM based image denoising method”  In _Multimedia Tools and Applications_ 77 Springer, 2018, pp. 30121–30134
  * [43] P.-H. Tournier et al.  “Microwave tomographic imaging of cerebrovascular accidents by using high-performance computing”  In _Parallel Computing_ 85, 2019, pp. 88–97  DOI: https://doi.org/10.1016/j.parco.2019.02.004
  * [44] R Scapaticci, OM Bucci, I Catapano and L Crocco  “Differential microwave imaging for brain stroke followup”  In _International Journal of Antennas and Propagation_ 2014 Hindawi, 2014


# The optical appearance of a nonsingular de Sitter core black hole geometry
under several thin disk emissions

I. De Martino<EMAIL_ADDRESS>Universidad de Salamanca, Departamento
de Fisica Fundamental, P. de la Merced, 37008 Salamanca, Spain  R. Della
Monica<EMAIL_ADDRESS>Universidad de Salamanca, Departamento de Fisica
Fundamental, P. de la Merced, 37008 Salamanca, Spain  D. Rubiera-Garcia
<EMAIL_ADDRESS>Departamento de Física Teórica and IPARCOS, Universidad
Complutense de Madrid, E-28040 Madrid, Spain

(August 30, 2024)

###### Abstract

We consider the optical appearance under a thin accretion disk of a regular
black hole with a central de Sitter core implementing
$\mathcal{O}(l^{2}/r^{2})$ far-corrections to the Schwarzschild black hole. We
use the choice $l=0.25M$, which satisfies recently found constraints from the
motion of the S2 star around Sgr A∗ in this model, and which leads to
thermodynamically stable black holes. As the emission model, we suitably adapt
ten samples of the Standard Unbound emission profile for a monochromatic
intensity in the disk’s frame, which have been previously employed in the
literature within the context of reproducing General Relativistic Magneto-
Hydrodynamic simulations of the accretion flow. We find the usual central
brightness depression surrounded by the bright ring cast by the disk’s direct
emission as well as two non-negligible photon ring contributions. As compared
to the usual Schwarzschild solution, the relative luminosities of the latter
are significantly boosted, while the size of the former is strongly decreased.
We discuss the entanglement of the background geometry and the choice of
emission model in generating these black hole images, as well as the
capability of these modifications of Schwarzschild solution to pass present
and future tests based on their optical appearance when illuminated by an
accretion disk.

##  I Introduction

In the last decade direct proofs for the existence of astrophysical black
holes have come from both the observations of dozens of gravitational waves
out of binary black hole [1] and neutron star [2] mergers, and from the
imaging brought by the superheated plasma in orbit around the supermassive
objects at the heart of the M87 [3] and Milky Way [4] galaxies, M87* and Sgr
A*, respectively. Both types of observations are compatible with the rough
expectations brought by the Kerr solution of General Relativity (GR), which is
solely characterized by mass and angular momentum [5]. At the same time, these
observations have provided a great opportunity for testing both modified black
holes and black hole mimickers. The former include black holes with additional
fields or hair (i.e. hairy black holes [6]) as well as black holes beyond GR
[7], while the latter typically refer to horizonless compact objects of
different types, such as traversable wormholes [8] or boson stars [9], see
e.g. [10] for a review of their current observational status.

A major argument in favour of these alternatives to the Kerr black hole comes
from the unavoidable existence of space-time singularities inside it, as
provided by the singularity theorems [11]. Such singularities are associated
to the existence of at least one incomplete geodesic, and their presence
threaten the classical determinism and predictability of GR itself. Since
geodesic incompleteness typically correlates with the blow up of a certain set
of curvature scalars, common wisdom in the field appeals to quantum-gravity
regularization mechanisms when curvature nears the Planck scale and restoring,
in turn, the geodesic completeness of the geometry. This way, the finding of
regular black holes has become a popular trend in the community, and many
different strategies to achieve this end have been concocted [12]. For the
sake of this paper, we consider the one of de Sitter (dS) cores, in which the
central point-like singularity of a spherically symmetric (Schwarzschild)
black hole is replaced by one such core, a strategy which manages to keep all
curvature scalars finite [13]. The regular de Sitter core black hole can be
later set to rotate [14, 15], in order to find the counterpart of the Kerr
solution.

The main aim of this work is to characterize one of these dS-core black holes
via its optical appearance when illuminated by a (optically and geometrically)
thin accretion disk. The chosen dS geometry is characterized by the fact that
in the far-away range it introduces corrections to the Schwarzschild black
hole of order $\mathcal{O}(l^{2}/r^{2})$, where $l$ is a typical scale of the
new geometry. This makes this dS core one of the strongest departures from the
Schwarzschild solution within this class. Indeed, the strong deviations in the
orbits of test particles, that characterize these dS-core black holes,
motivated a Monte Carlo Markov chain analysis aimed to constrain the scale
length $l$ using astrometric and spectroscopic data of positions and velocity
of the star S2 that orbits around the supermassive compact object Sgr A* at
the centre of the Milky Way [16, 17]. These data provided an upper bound on
the scale length $l$ of $\lesssim 0.47M$ (in units of $G=c=1$) at the 95% of
the confidence level. For the sake of this work we shall take the choice
$l=0.25M$, which lies within such a constraint and in the lower end of the
branch of thermodynamically stable black holes within this class.

To carry out our analysis we shall suitably adapt ten samples of the Standard
Unbound (SU) Johnson’s distribution for a monochromatic (in the frame of the
disk) emission profile [18], some of which have been used in several recent
surveys of General Relativistic Magneto HydroDynamic (GRMHD) simulations of
the accretion flow within simpler settings. Such simulations consistently
report the optical appearance of black holes to be largely dominated by a
bright ring of radiation enclosing a central brightness depression, the latter
closely tracking, in the observer’s plane, a critical curve resulting from the
projection of the unstable bound geodesics (the photon sphere), and typically
known as the shadow [19]. Furthermore, within this thin accretion disk
framework, such a ring is broken into an infinite sequence of self-similar
rings labelled by an integer $n$ corresponding to the number of half-loops of
light trajectories around the black hole, creating, on the asymptotic
observer’s screen, a similar number of photon rings [20]. These rings produce
universal signatures in very-long base interferometry (VLBI) detectors [21],
potentially allowing to distinguish between different black hole metrics [22].
However, due to the fact that such rings are exponentially decreased in their
corresponding luminosities, typically those beyond $n\geq 2$ are dismissed
from the images since only the $n=1$ and (perhaps) the $n=2$ ones might be
observed in future observational devices, such as the next generation Event
Horizon Telescope (ngEHT) [23]. In addition, in this setting the outer edge of
the central brightness depression is associated with the (gravitational
lensed) effective region of emission (typically the horizon) and thus it can
be strongly reduced as compared to the usual shadow of the geometrically thick
case [24].

The features above have been confirmed in many other studies in the field, see
e.g. [25, 26, 27, 28, 29, 30, 32, 31, 33, 34, 35, 36] for a non-extensive
list. The dS-core black holes images studied in this work with the SU emission
profiles reproduce this behaviour, but significantly decrease the luminosity
extinction rate between the $n=1$ and $n=2$ photon rings, making the latter to
appear comparatively much brighter in the images. At the same time, a strong
reduction in the shadow’s size is also observed. We discuss the interplay
between the ten SU models and the dS-core background geometry in generating
the photon ring and central brightness depression features of the images. We
furthermore comment on the resemblances and differences of these images with
respect to its Schwarzschild counterpart as well as the chances of such black
holes to pass present and future observational constraints of this kind.

This work is organized as follows. In Sec II we introduce the theoretical
background and the thin accretion disk model and emission profiles used
throughout the paper. Images are found and discussed in Sec. III for both dS
black holes and naked cores, and in Sec. IV we further comment on our results.
Appendix A contains the main ingredients of null geodesics for the generation
of images. Throughout the article we shall work in units $G=c=1$.

##  II Framework

###  II.1 Background geometry

We consider an asymptotically flat, spherically symmetric geometry written as

$ds^{2}=-A(r)dt^{2}+\frac{dr^{2}}{A(r)}+r^{2}d\Omega^{2}\,,$ (1)

where $d\Omega^{2}=d\theta^{2}+r^{2}\sin^{2}\theta d\phi^{2}$ is the volume
element on the unit two-spheres. A large family of dS cores parameterized by
the metric function (here $l$ is a constant)

$M(r)=M\left(\frac{r^{q}+l^{q}}{r^{q}}\right)^{\frac{p}{q}}\,,$ (2)

was introduced in [14], since it contains as particular members the well known
Bardeen [37] ($p=3,q=2$) and Hayward solutions [38] ($p=q=3$), and are
typically found as solutions of non-linear electrodynamics (possibly coupled
to scalar fields as well); the condition $p\geq 3$ must be imposed to
guarantee regularity at the center $r=0$. For the sake of this work we take
the choice $p=3,q=1$, that is

$M(r)=M\frac{r^{3}}{(r+l)^{3}}\,.$ (3)

This choice represents an asymptotically flat space-time, behaving in this
limit as $A(r)\approx 1-2M/r+\mathcal{O}(l^{2}/r^{2})$, while at the center
$r\to 0$ it implements a dS core behaviour of the form $A(r)\approx
1-2Mr^{2}/l^{3}+\mathcal{O}(r^{3}/l^{3})$. The motivation for this model lies
therefore in the fact that it provides the strongest polynomial correction to
the Schwarzschild black hole at far distances within this class, while
preserving the dS core at the center and the regularity of the curvature
scalars it brings with it.

The metric (3) represents different configurations depending on the choice of
the parameter $l$. In this sense, for $0\leq l<l_{c}$, where
$l_{c}=\tfrac{8M}{27}\approx 0.2969M$, it describes two-horizons black holes
with an external event horizon $r_{+}$ and an internal one $r_{-}$. These two
horizons merge into a single (degenerate) one (i.e. an extreme black hole) at
$l=l_{c}$. For $l>l_{c}$ the horizon disappears, leaving a regular naked dS
core. Furthermore, by analyzing the thermodynamic properties of the black
holes above, one finds that the stability of the system, as given by the
positivity of the specific heat (at constant $l$) can only be achieved within
the range $0.245\leq l\leq l_{c}$, thus leaving a narrow gap of stable black
holes within this model. Furthermore, the metric (3) has been constrained
using the astrometric and spectroscopic data of the orbital motion of the S2
star around Sgr A∗ [16]. Their analysis constrains the value of $l$ to be
$\lesssim 0.47M$ at the 95% of the confidence level, therefore also allowing
for the existence of naked dS cores compatible with such an observation.

For our aims, the main object of interest will be the generation of images of
a sample of two-horizons black hole with a dS core (dSBH), taking the choice
of $l=0.25M$, which yields the event horizon radius $r_{+}\approx 1.059M$ and
an inner horizon radius $r_{-}=0.25M$. For completeness, we shall also discuss
briefly the images of a naked dS core (NdS), for which we take the choice
$l=0.30M$; obviously in such a case no horizon radius is present.

A brief description of light trajectories in spherically symmetric space-times
is provided in Appendix A. For our analysis below the important concept is the
one of photon sphere, namely, the locus of unstable bound geodesics, and the
critical curve, namely, the projection on the observer’s plane of such a
photon sphere.

###  II.2 Accretion disk emission

We consider our source of emission to be an optically thin, equatorial thin
disk near the black hole. Such a disk is typically modelled as composed of
individual emitters that follow approximately circular (Keplerian) orbits
until reaching the innermost stable circular radius (ISCO) radius; after that
point the particles of the disk plunge into the black hole following
Cunningham’s prescription [39]. We assume the disk to define a specific
monochromatic intensity $I_{\nu_{e}}$, with $\nu_{e}$ the frequency in the
frame of the disk, and the total intensity emitted is thus only dependent on
the radial distance from the central mass, $I_{e}(r)\equiv I_{\nu_{e}}$. The
reconstruction of the image requires, in general, solving the radiative
transfer equation [40]

$\frac{d\mathcal{I}_{\nu}}{ds}=\frac{j(\nu)}{\nu^{2}}-\nu\mathcal{I}_{\nu}\chi(\nu)\,,$
(4)

where $\mathcal{I}_{\nu}\equiv I_{\nu}/\nu^{3}$ is the invariant intensity,
while $j_{\nu}$ is the emissivity and $\chi(\nu)$ the absorption, the latter
effectively taken to zero by the optically thin assumption. The invariant
character of $\mathcal{I}_{\nu}$ means that in the frame of the observer with
frequency $\nu_{o}$ the observed intensity will behave as

$I_{\nu_{o}}=\frac{\nu_{o}^{3}}{\nu_{e}^{3}}I_{\nu_{e}}=\frac{g(r_{o})}{g(r_{\infty})}I_{\nu_{e}}\,,$
(5)

where the last equality follows from Eq.(1) with $g(r)=A(r)^{1/2}$. Since we
are dealing with asymptotically flat space-times with an observer located far
away from the black hole, then $g(r_{\infty})=1$. The total observed intensity
will be then given by the integration

$I_{o}(r)=\int d\nu_{o}I_{\nu_{0}}(r)=\int
d\nu_{e}g^{3}(r)I_{\nu_{e}}=g^{4}I(r)\,.$ (6)

Every time a light ray intersects the (optically thin) disk it will increase
its luminosity as dictated by the disk’s local intensity, and thus the final
image brightness will be given by the addition of all the contributions for
every possible intersection with the disk, in other words, we write

$I_{o}(b)=\sum_{n}\xi_{n}(A^{2}I)|_{r=r_{n}(b)}\,,$ (7)

where the function $r_{n}(b)$ with $n=0,1,2,\ldots$ stores the information
about the radial location of the $n$-th intersection with the disk outside the
event horizon, and it is dubbed as the transfer function. In this sense, $n=0$
corresponds to the direct emission of the disk, namely, those photons emitted
from the disk and directly reaching the observer without undergoing any
additional turns around the black hole; those with $n\geq 1$ correspond to the
photon rings, i.e., photons that have completed $n$ half-turns around the
black hole (since the disk is located both in the front and the back of it).
Additionally, we have included a fudge factor $\xi_{n}$, which has been shown
to increase the compatibility of the simple modelling of the emission employed
here with the results of time-averaged GRMHD simulations in non-zero thickness
disks [18, 24]; depending on the sought properties of the disk this numerical
factor can be modelled according to different prescriptions. For our purposes
we shall conform to the thin-disk scenario and thus take $\xi_{n}=1$ for
$n=0,1,2$ and $\xi_{n}=0$ otherwise.

On theoretical grounds, the contribution of the photon rings to the total
luminosity can be proven to be exponentially-suppressed. To see this, one
considers a photon that starts slightly above the photon sphere $r_{m}$ as
$r_{0}=r_{m}+\delta r_{0}$, where $\delta r_{0}\ll r_{m}$. After $n$ half-
turns, the photon’s location will be given by [41]

$\delta r_{n}=e^{\gamma_{L}n}\delta r_{0}\,,$ (8)

where $\gamma_{L}$ is dubbed as the Lyapunov exponent, a universal quantifier
of a given background geometry in the limit $n\to\infty$. This exponential
drift of the photon’s location induces a similar exponential decay for the
luminosities of the corresponding photon rings, hence its relevance in
connecting properties of the background geometry with actual observables (up
to the non-trivial modelling of the accretion disk, as we shall see below).
From an observational point of view, however, a VLBI observatory can only
measure luminosity contrasts between the direct emission $n=0$ and the first
photon ring $n=1$, while future experiments such as the ngEHT could allow to
measure up the second photon ring $n=2$ [23]. Therefore, we shall keep up to
the $n=2$ ring in the generation of images and the computation of the Lyapunov
exponent above.

Model | $\mu$ | $\gamma$ | $\sigma/M$ | Comments
---|---|---|---|---
SU1 | $2r_{+}$ | $-2$ | $3/2$ | Softest decay of [42]
SU2 | $17r_{+}/6$ | $-2$ | $1/4$ | This is GLM3 model of [43] for a Schwarzschild black hole.
SU3 | $\frac{3}{2}r_{+}$ | $2$ | $1/4$ | Narrow rings of [44]
SU4 | $r_{-}$ | $-2$ | $3/2$ | Overestimated of [44]
SU5 | $r_{+}$ | $0$ | $1$ | Typical 1 of [44]
SU6 | $r_{+}$ | $2$ | $1/2$ | Underestimated of [44]
SU7 | $r_{-}$ | $-3/2$ | $1/2$ | Considered in [18], GLM1 model of [43]
SU8 | $r_{-}$ | $0$ | $1/2$ | This is GLM2 model of [43]
SU9 | $r_{-}$ | $2$ | $1$ | Typical 2 of [44]
SU10 | $r_{-}$ | $2$ | $1/4$ | Hardest decay of [42]
Table 1: The SU class of emission profiles defined in Eq.(9) for the ten
choices of the parameters $\mu$, $\gamma$, and $\sigma/M$ characterizing them.
Here $r_{\pm}$ denote the horizon and inner radius of the black hole; for
Schwarzschild this is $r_{+}=2M$ and $r_{-}=0$, while for dSBH this is
$r_{+}\approx 1.059M$ and $r_{-}=0.25M$. In the rightmost column we display
the references each model is inspired from.

Figure 1: The ten (normalized) emission profiles of Table 1 for the
Schwarzschild (left) and dSBH (right) configurations in double-logarithmic
scale. The three dashed vertical lines denote the event horizon radius $r_{+}$
($r_{h}$ in the Schwarzschild case representing the single horizon in this
geometry), the photon sphere radius $r_{m}$ and the ISCO radius $r_{\rm
ISCO}$, respectively, for each configuration; for the dSBH a fourth vertical
dashed line denotes the location of the inner horizon radius $r_{-}$.
Obviously, only those parts of the emission profile at radius $r>r_{+}$ are
relevant for the generation of images.

As our class of emission profiles we select Johnson’s SU models, given by the
expression [18]

$J_{SU}=\frac{e^{-\frac{1}{2}\left[\gamma+\mathop{\rm
arcsinh}\nolimits\left(\frac{r-\mu}{\sigma}\right)\right]^{2}}}{\sqrt{(r-\mu)^{2}+\sigma^{2}}}\,,$
(9)

which is characterized by three parameters; $\mu$ is related to the location
of the peak of the emission: small (big) values tend to move the peak away
from (behind) the event horizon; $\gamma$ governs the profile asymmetry:
negative (positive) values tend to locate the steep part of the profile inside
(outside) the horizon; $\sigma$ (typically given in units of $M$) is related
to the profile width; small (big) values tend to have a steeper (broader)
decay with narrower (broader) rings. For the sake of this work we shall take
ten choices for these values covering a range of qualitatively different
profiles, and which are taken from a number of different sources where the
parameters range from $\mu\in[r_{-},2r_{+}]$, $\gamma\in[-2,2]$ and
$\sigma/M\in[0.25,1.5]$, where $r_{\pm}$ are the outer and inner horizons of a
rotating black hole. Despite the lack of rotation in our model, we shall use
the outer/inner horizons of the dSBH configuration as a substitute for them,
while in the Schwarzschild black hole where only the event horizon is present,
we obviously have $r_{+}=2M$ and $r_{-}=0$.

Figure 2: The transfer function $r_{n}(b)$ for the Schwarzschild black hole
(left), dSBH (middle) and NdS (right) solutions. The contributions of the
direct emission correspond to blue straight lines, while the $n=1$ and $n=2$
photon rings are shown in orange and green, respectively.

This pool of parameters is detailed in Table 1, where the different models are
presented in decreasing order in terms of the radius $r_{e}$ of the region
where they possess the effective source of emission i.e. its maximum. The
corresponding radial dependence of these profiles is depicted in Fig. 1 for
the Schwarzschild black hole (left) and the dSBH (right), including their
respective values of the event horizon radius $r_{+}$, the inner one $r_{-}$
(in the dSBH case), the photon sphere radius $r_{m}$ (associated to the locus
of unstable geodesics, see Appendix A) and the ISCO. Three of these choices
are adapted from those originally introduced by Gralla, Marrone and Lupsasca
(GLM) in Ref. [18], and employed for the generation of images in several
papers by some of us, such as in [43]; here they are dubbed as SU2, SU7, and
SU8, respectively. Five of the remaining models were introduced by Paugnat et
al. in [44] as representative cases of a survey performed over black hole
spins and inclinations for combination of values of the emission profile
parameters $\mu=\\{r_{-},r_{+}/2,r_{+},3r_{+}/2,2r_{+}\\}$,
$\gamma=\\{-2,-1,0,1,2\\}$ and $\sigma/M=\\{0.25,0.5,1,1.5\\}$, which cover
the range of results found in the corresponding simulations performed by such
a collaboration. Such models represent a “Typical 1” profile with a $n=2$ ring
diameters slightly above average values (SU5), a “Typical 2” profile with
diameters slightly below average values (SU9), an overestimated profile (SU4)
and an underestimated one (SU6), and finally a narrow profile with the
smallest diameters (SU3). The remaining two choices come from the analysis of
[42] and, in particular, from their Table 2, in which they present over a
hundred of emission profiles, highlighting three of them. The first one (SU1)
corresponds to the largest radius of emission and a slow decay, while the
second one (SU10) lies in the opposite end, with a very quick decay. The third
model is the GLM1 one already included in our list as the SU7 one.

This pool of models is flexible enough to cover different scenarios of
behaviour of the radial dependence around the effective region of emission of
the disk. In this sense, four models peak outside the event horizon of the
Schwarzschild black hole: SU1 and SU2 near the ISCO radius with soft and hard
decays, respectively, and SU3 and SU4 slightly above the Schwarzschild radius,
with hard and soft decays, respectively. Another one (SU5) peaks almost
exactly at $r=r_{+}$. The remaining SU6/SU7/SU8/SU9/SU10 models peak inside
the event horizon, so the effective radius of emission is marked by their
values at $r=r_{+}$. These last five models have increasingly harder decays,
which means that for those in the end of the sequence (SU8, SU9, and SU10),
there will be little emission outside the event horizon, something that should
have its reflection in the corresponding optical appearances. Given the fact
that $r_{\pm}$ are different in the Schwarzschild/dSBH solutions, there are
some differences in the shape of the corresponding profiles; this seems to be
an unavoidable trouble of comparing different background geometries: the
behaviour of the disk is not independent of the shape of the metric, see [15]
for a recent discussion on this. Accepting this, we normalize all our
intensity profiles to their maximum values for as an optimal comparison of the
Schwarzschild/dSBH images as possible.

##  III Generation of images

Model | $\gamma_{L}$ | $e^{\gamma_{L}}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU1}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU2}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU3}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU4}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU5}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU6}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU7}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU8}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU9}$ | $\left(\frac{I_{1}}{I_{2}}\right)^{SU10}$
---|---|---|---|---|---|---|---|---|---|---|---|---
Sch | $3.1507$ | $23.35$ | $28.94$ | $27.83$ | $21.57$ | $27.21$ | $23.34$ | $21.14$ | $24.74$ | $23.45$ | $22.81$ | $22.20$
dSBH | $2.266$ | $9.64$ | $14.06$ | $12.51$ | $8.00$ | $13.14$ | $10.07$ | $8.20$ | $10.82$ | $9.61$ | $9.14$ | $8.47$
Table 2: The Lyapunov index $\gamma_{L}$, the theoretical extinction rate
$e^{\gamma_{L}}$, and the observational extinction rate $I_{1}/I_{2}$ between
the $n=1$ and $n=2$ photon rings for the ten SU models listed in Table 1, as
obtained for the Schwarzschild and the dSBH models.

The observational appearances of the Schwarzschild, dSBH and NdS
configurations are built using a ray-tracing procedure: all trajectories
arriving to the observer’s plane image are backtracked into the emitting
region via the geodesic equation (recall that a short derivation of the
equations of null geodesics is provided in the Appendix A)

$\frac{d\phi}{dr}=-\frac{b}{r^{2}}\frac{1}{\sqrt{1-b^{2}\frac{A(r)}{r^{2}}}}\,.$
(10)

The relevant quantities to carry out this analysis and to characterize the
resulting features are the critical impact parameter and photon sphere radius.
For Schwarzschild these are $b_{c}=3\sqrt{3}M\approx 5.196M$ and $r_{m}=3M$,
while for dSBH $b_{c}\approx 3.614M$ and $r_{m}\approx 1.767M$, and, finally,
for NdS one has $b_{c}\approx 3.108M$ and $r_{m}\approx 1.329M$.

We place the observer’s at a radius of $r_{\infty}=1000M$ to recover the
asymptotically Schwarzschild character of the space-time there, and run the
integrations to collect the info on the transfer function, which is depicted
in Fig. 2 for the Schwarzschild (left), dSBH (middle) and NdS (right)
configurations; to this end we assume the disk to be seen in a face-on
orientation. This figure displays the contributions of the direct emission
(blue straight line) and the $n=1$ (orange) and $n=2$ (green) photon rings:
the effective source of emission $r_{e}/M$ determines the parts of these
curves above which their luminosity can actually reach the asymptotic
observer. Since the slopes of the curves are associated with the
demagnification of their corresponding emissions in the images, this figure
states the dominant role of the direct emission in the image via its nearly
constant slope; the $n=1$ and $n=2$ photon rings will be sub-dominant, while
$n>2$ rings will be ignored. In the NdS case we observe that both photon rings
extend much deeper into the inner region of the configuration; this is due to
the absence of an event horizon which allows photons to reach the innermost
regions of the effective potential, where an infinite divergence is found,
thus repelling light rays that get to it. We shall see in a while the
repercussions of this fact for the corresponding optical appearances.

###  III.1 Regular dS black hole vs Schwarzschild

Figure 3: The observed intensity $I_{o}(b)$ and the optical appearance (using
$b\in[-8,8]$) of a Schwarzschild black hole (left two figures) and a dSBH
solution (right two figures) for the (from top to bottom) SU1, SU2, SU3, SU4,
and SU5 emission models of Table 1, using a fudge factor
$\xi_{0}=\xi_{1}=\xi_{2}=1$ and $\xi_{n}=0$ for $n>2$.

Figure 4: The observed intensity $I_{o}(b)$ and the optical appearance (using
$b\in[-8,8]$) of a Schwarzschild black hole (left two figures) and a dSBH
solution (right two figures) for the (from top to bottom) SU6, SU7, SU8, SU9,
and SU10 emission models of Table 1, using a fudge factor
$\xi_{0}=\xi_{1}=\xi_{2}=1$ and $\xi_{n}=0$ for $n>2$.

In Table 2, we report the numerical results of our simulations for the
Schwarzschild and dSBH configurations. Our main interest here is to check the
reliability of the Lyapunov index $\gamma_{L}$ of Eq.(8) as a marker of the
actual exponential decay of the luminosity of successive photon rings. To this
end, we compute the extinction rate of the $n=1$ ring as compared to the
$n=2$, that is, the luminosity ratio $I_{1}/I_{2}$ for each SU model, and
compare it to the theoretical rate $e^{\gamma_{L}}$111For the latter we do not
use the $n\to\infty$ limit where (for the Schwarzschild black hole)
$\gamma_{L}=\pi$, but approximate it by the $n=2$ value (for Schwarzschild
$\gamma_{L}\approx 3.150$). This difference ($\sim 0.3\%$) is negligible
enough for our purposes here .. From this Table, we observe that in both the
Schwarzschild and dSBH geometries five SU models produce higher extinction
rates than its Lyapunov-based rate, one is left unchanged, while the other
four predict a lower rate. Such deviations from the theoretical prediction can
be quite important: in the SU1 model it can increase up to a $\sim 24\%$ for
Schwarzschild, while in the SU3/SU6 it can decrease by a factor $\sim 10\%$
for Schwarzschild and $\sim 15\%$ for dSBH. This behaviour is strongly
governed by the disk’s parameter $\gamma$: negative (positive) values of
$\gamma$ yield larger (smaller) extinction rates; in this sense, the SU5 and
SU8 models, having $\gamma=0$, produce extinction rates that track very
closely the Lyapunov prediction. This is in agreement with the results of
GRMHD simulations, which associate the size of $\gamma$ with the thickness of
the rings and, in turn, with their relative luminosity, as we have just seen
here.

In Figs. 3 and 4 we depict the observed intensity $I_{o}(r)$ and the optical
appearance for a range of values for the impact parameter $b\in[-8,8]$ for the
Schwarzschild (left two plots) and dSBH (right two plots) solutions, and for
the SU1/SU2/SU3/SU4/SU5 and SU6/SU7/SU8/SU9/SU10 models, respectively, and we
recall that such models are organized in decreasing values of their effective
region of emission. There are several features of interest in these images. We
find universally (i.e. for every SU model, and in both the Schwarzschild and
dSBH configurations) the expected presence of both a central brightness
depression and a bright photon ring surrounding it, which can be further
decomposed into up to two additional rings (the photon rings) depending mostly
on the choice of SU model but also with some differences between the
Schwarzschild/dSBH solutions.

Regarding the bright rings, their distribution is governed by both the
location of the effective source of emission and the strength of the decay
with distance of the emission profile. We can see this effect in the SU1/SU2
models, where the location of the maximum of emission is very similar but the
decay is weak (strong) in the SU1 (SU2) model: as a consequence in the former
(latter) the glow of the direct emission extends to much farther (closer)
regions of the impact parameter space. The photon rings ($n=1,n=2$) will be
clearly visible in the Schwarzschild-SU2 model, but they appear overlapped
with the direct emission in the Schwarzschild-SU1 and dSBH-SU1/SU2 models,
which means we are able to distinguish between both background geometries at
fixed SU emission model. Both SU3 and SU6 produce an unusual signature for
both Schwarzschild/dSBH models with the $n=1$ photon ring located above the
direct emission but separated from it, so in the optical appearance one neatly
sees the photon ring encircling the direct emission; this effect is more
exaggerated in the SU3 model given the fact that the decay of the profile with
distance is stronger than in SU6. In every other SU model the photon rings are
overlapped with the direct emission. In particular, for the SU4 it nears the
outer edge of the brightness depression while the direct emission extends to
farther distances given the weak decay of the intensity profile with the
distance (as in SU1). Finally, SU5, SU7, SU8, SU9, and SU10 produce quite
similar images: the direct emission dominates the image from the central
brightness depression outwards, and superimposed on it lies the boost of
luminosity of the two photon rings, whose exact locations within the direct
emission depend on the strength of the decay of the intensity profile which
defines the effective domain of the direct emission: weak (strong) decays
produce a photon ring placed far (near) the outer edge of the image.

Figure 5: The observed intensity $I_{o}(r)$ (top) and the optical appearance
(bottom) of a naked regular dS core for the SU2 (left), SU7 (middle), and SU8
(right) emission models of Table 1, using a fudge factor
$\xi_{0}=\xi_{1}=\xi_{2}=1$ and $\xi_{n}=0$ for $n>2$.

Differences between the dSBH and Schwarzschild solutions (at fixed SU model)
are both obvious at naked eye but also subtle. The first significant
difference is a strong reduction of the central brightness depression size in
the dSBH case as compared to the Schwarzschild one. This was expected on the
grounds of a much smaller horizon radius, critical impact parameter and photon
sphere radius of the dSBH, since these aspects strongly influence the actual
location of the outer edge of the brightness depression, and which in the case
where the emission goes all the way down to the event horizon is bounded by
above by the so-called inner shadow, and which, as already discussed, depends
only on the background geometry [24]. The second difference is with the photon
rings themselves: we already computed in Table 2 a significant reduction in
the extinction rate $I_{1}/I_{2}$ between both geometries (up to a factor 2),
which has a clear reflection of brighter rings as compared to the surrounding
image. Besides these aspects, the arrangement of the locations and overlapping
of the photon rings with the direct emission in the dSBH has some tiny
differences with respect to the Schwarzschild case but mostly in those SU
models with emission outside the event horizon; those disks whose emission
decreases from the horizon outwards tend to wash out these differences between
background geometries. Let us recall that these images are obtained in a face-
on orientation, $i=0$: inclined images would make the lower part of the photon
ring to tilt upon the central brightness depression introducing an asymmetry
in the images (though only noticeable at large enough inclinations). However,
this effect would be similar for both Schwarzschild and dSBH and would help
little to distinguish them beyond what we discussed here; an illustration of
this effect can be found in e.g. [43].

###  III.2 Regular naked dS core

For the sake of completeness of our analysis, let us briefly discuss an
example of a regular naked dS core, for which we choose the value $l=0.30$,
which brings it very close to becoming a (extreme) black hole. In this case,
the lack of a horizon allows light trajectories to interact with the inner
part of the effective potential, which has a divergent behaviour at the
center. This has a dramatic effect in the images of the corresponding objects.
To illustrate this, we depict in Fig. 5 the observed intensity and the optical
appearances for three choices of the SU models (SU2, SU7, and SU8),
corresponding to the GLM3/GLM1/GLM2 used in previous papers by some of us.
There one neatly sees the main difference brought by these NdS objects: the
appearance of new peaks of intensity, which manifest themselves as new photon
rings in the optical appearances. The reason for this can be seen in the
transfer function of Fig. 2 (right plot), where photons interacting with the
internal part of the potential will be repelled by it and have the chance to
circulate one more time around the photon sphere, producing the new sets of
photon rings in the observer’s screen. Furthermore, such new contributions
break the expected exponentially decay of successive photon rings, which makes
rings with $n>2$ to produce non-negligible contributions to the luminosity of
the images, and yield a multi-ring structure previously studied in the
literature [45, 46], though for simplicity we just kept up to the $n=2$ ring
to illustrate the main aspects of the optical appearances of the NdS
solutions. Indeed, the latter resemble the canonical appearance of a black
hole for the SU2 model, with the new photon rings living inside the region
previously occupied by the central brightness depression only, yielding a
clear signature to distinguish both objects. But in the SU7 and SU8 models the
resulting objects yield optical appearances that have nothing to do with the
canonical black hole ones and, as such, these kinds of objects can hardly pass
as black hole mimickers from what we know about their cast images.

##  IV Conclusion

In this work we have considered a black hole geometry replacing the central
region of the Schwarzschild black hole by a de Sitter core characterized by a
single parameter $l$, this way removing the curvature singularity lurking
there, and studied its optical appearance when surrounded by an optically and
geometrically thin accretion disk. This de Sitter core is actually the
strongest modification of the Schwarzschild solution of a certain class of
such cores which include the well known Bardeen and Hayward solutions as
particular members.

We set a value of $l=0.25M$ that is allowed by the recent analysis of the
orbital motion of the S2 star around Sgr A*, that estimated $l\lesssim 0.47M$
at the 95% of the confidence level, and corresponds to the lower end of the
sub-family of these black holes which are thermodynamically stable. Regarding
the accretion disk, we use a simplified setting for a monochromatic intensity
profile belonging to the class of models known as Standard Unbound Johnson’s
distribution, a sub-class of which has been shown in the literature to be
capable to reproduce the results of some scenarios for the accretion flow in
GRMHD simulations. In particular, we employed ten choices for the parameters
of such profiles, representatives of different emission profiles and producing
different signatures of the corresponding optical appearances.

We run our simulations in this setting, which for the background geometry
means significantly lower horizon radius, critical impact parameter, and
photon sphere radius. This is translated, within the optical appearances of
these objects, into a strong suppression of the luminosity extinction rate (up
to a factor two) between the $n=1$ and $n=2$ photon rings, which therefore
appear in the images as significantly more luminous within the direct emission
of the disk as compared to their Schwarzschild counterparts. This happens for
all the emission models, yet visible differences between dSBH/Schwarschild
configurations are more acute in those SU model whose peak of emission is
located away from the event horizon; those in which the profile decreases from
the event horizon outwards tend to wash out such differences and produce more
similar images.

Furthermore, the dS black hole also yields a strongly reduced size of the
outer edge of the central brightness depression, neatly distinguishing it from
the Schwarzschild black hole. Under the recent results from the calibrated
measurements of the “shadow’s size” of Sgr A∗ by the EHT Collaboration
reported in [47], where it is claimed that the size of the bright ring
correlates, under certain assumptions and after proper calibration (and some
debate) with such an edge, one expects the latter to be constrained (at
$2\sigma$) within the range $4.21\lesssim r_{sh}/M\lesssim 5.66$, while in our
case this is $r_{sh}/M\approx 3.614$. This strategy was actually used in [48]
to constrain a large number of alternative spherically symmetric geometries
(further analysis of the implications of such a measurement has been recently
discussed in [49]). Given the strong reduction of the size of the central
brightness depression in our case, this makes these objects hardly compatible
with the images of Sgr A∗, pinpointing that strong deviations from the GR
solutions may have a hard time in reconciling themselves with present and
future observations. Obviously, such a criticism applies even more strongly to
the naked dS configurations since they also show a multi-ring structure that
does not resemble current observations, except for very particular emission
model scenarios.

To conclude, our analysis with a pool of well motivated emission models shows
that regular de Sitter-core black holes can be actually tested with current
and feature “shadow” observations since they introduce large enough
differences in their optical appearances as compared to canonical black hole
solutions. Even though the background geometry employed here seems to be too
strong a modification of the Schwarzschild one to be compatible with current
images, it paves the way for other dS core black holes to be tested with this
tool. Of course, in order to produce competitive images, the simplified
treatment employed here needs to be upgraded towards more realistic settings
with rotation, thick disks, absorption, inclination, and so on.

## Acknowledgments

IDM and RDM acknowledge support from the grant PID2021-122938NB-I00, and DRG
acknowledges support from the grant PID2022-138607NB-I00, both funded by
MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”.

##  Appendix A Null geodesic behaviour

In this Appendix we provide a short summary of the main ingredients (since
they have been derived many times in the literature of the subject, see e.g.
[51] for a quite general analysis) of the equations of null geodesic motion
together with their main features for the generation of images. We start from
the Lagrangian density of a particle, which is given by

$\mathcal{L}=\frac{1}{2}g_{\mu\nu}\dot{x}^{\mu}\dot{x}^{\nu}\ ,$ (11)

where a dot represents a derivative with respect to the affine parameter. For
a photon with wave function $k^{\mu}=\dot{x}^{u}$ this Lagrangian equals zero,
a reflection of the fact that photons travel along null geodesics of the
background geometry, i.e., $g_{\mu\nu}k^{\mu}k^{\nu}=0$. The spherically
symmetric character of the geometry means that we can fix $\theta=\pi/2$
without any loss of generality, while the existence of two Killing vectors
associated to translational and rotational invariances of the system implies
two conserved quantities, $E=-A\dot{t}$ and $L=r^{2}\dot{\phi}$, naturally
interpreted as the energy and angular momentum per unit mass, respectively.
This allows to write the geodesic equation under the form

$\dot{r}^{2}=\frac{1}{b^{2}}-V(r)\ ,$ (12)

where $b=L/E$ is the impact parameter, while the effective potential reads as

$V(r)=\frac{A(r)}{r^{2}}\ .$ (13)

Eq.(12) has the typical shape of a particle moving in a one-dimensional
effective potential. From here one can define the photon sphere as given by
the simultaneous fulfilment of the equations

$\displaystyle b_{c}^{2}=\frac{1}{V(r_{m})}$ (14) $\displaystyle
V^{\prime}(r)|_{r=r_{m}}=0$ (15) $\displaystyle
V^{\prime\prime}(r)|_{r=r_{m}}<0\ .$ (16)

The first condition states that a turning point has been reached, the second
its critical (extremum) character, and the third that it corresponds to a
maximum (otherwise, it is dubbed as an anti-photon sphere, which are typically
unstable [50]). Provided that all these conditions are met then the photon
sphere radius $r_{m}$ (which is projected in the observer’s screen as a
circular critical curve) yields the locus of bound unstable geodesics, with
$b_{c}$ corresponding to the critical impact parameter a light ray issued from
the observer’s screen would need in order to asymptote to $r=r_{m}$. This way,
$b_{c}$ separates those trajectories with $b>b_{c}$ that find a turning point
at some $r_{t}>r_{m}$, from those with $b<b_{c}$ that overcome the effective
potential and end up swallowed by the event horizon of the black hole. In the
backwards ray-tracing procedure, it splits the observer’s plane image into the
bright ($b>b_{c}$) and dark ($b<b_{c}$) regions, with those nearing $b\gtrsim
b_{c}$ being associated to the presence of photon rings, namely, light
trajectories that have circled the black hole $n$-half times. In order to
classify the latter and to find the optical appearance of the object one is
interested in obtaining the deflection angle; this is done using again the
conserved quantities of the system to rewrite (12) as a variation of the
azimuthal angle $\phi$ with respect to $r$, which results in Eq.(10).
Nonetheless, in thin disk scenarios where emission from inside the photon
sphere (up to the event horizon) can occur, this simple picture must be
modified in order to allow light rays to get out from the accretion disk
towards reaching the asymptotic observer; this is described in the main text.

## References

  * [1] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 116 (2016) 061102.
  * [2] B. P. Abbott et al. [LIGO Scientific and Virgo], Phys. Rev. Lett. 119 (2017) 161101.
  * [3] K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 875 (2019) L1.
  * [4] K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 (2022) L12.
  * [5] R. P. Kerr, Phys. Rev. Lett. 11 (1963) 237.
  * [6] C. A. R. Herdeiro and E. Radu, Int. J. Mod. Phys. D 24 (2015) 1542014.
  * [7] R. Carballo-Rubio, F. Di Filippo, S. Liberati and M. Visser, Phys. Rev. D 98 (2018) 124009.
  * [8] A. Simpson and M. Visser, JCAP 02 (2019) 042.
  * [9] F. H. Vincent, Z. Meliani, P. Grandclement, E. Gourgoulhon and O. Straub, Class. Quant. Grav. 33 (2016) 105015.
  * [10] V. Cardoso and P. Pani, Living Rev. Rel. 22 (2019) 4.
  * [11] J. M. M. Senovilla and D. Garfinkle, Class. Quant. Grav. 32 (2015) 124008.
  * [12] C. Lan, H. Yang, Y. Guo and Y. G. Miao, Int. J. Theor. Phys. 62 (2023) 202.
  * [13] S. Ansoldi, [arXiv:0802.0330 [gr-qc]].
  * [14] J. C. S. Neves and A. Saa, Phys. Lett. B 734 (2014) 44.
  * [15] K. Boshkayev, et al. [arXiv:2307.15003 [gr-qc]].
  * [16] M. Cadoni, et al. Phys. Rev. D 107 (2023) 044038.
  * [17] M. de Laurentis, I. De Martino and R. Della Monica, Rept. Prog. Phys. 86 (2023) 104901.
  * [18] S. E. Gralla, A. Lupsasca and D. P. Marrone, Phys. Rev. D 102 (2020) 124004.
  * [19] H. Falcke, F. Melia and E. Agol, Astrophys. J. Lett. 528 (2000) L13.
  * [20] S. E. Gralla, D. E. Holz and R. M. Wald, Phys. Rev. D 100 (2019) 024018.
  * [21] M. D. Johnson, et al. Sci. Adv. 6 (2020) eaaz1310.
  * [22] M. Wielgus, Phys. Rev. D 104 (2021) 124058.
  * [23] P. Tiede, et al., Galaxies 10 (2022) 111.
  * [24] A. Chael, M. D. Johnson and A. Lupsasca, Astrophys. J. 918 (2021) 6.
  * [25] R. Narayan, M. D. Johnson and C. F. Gammie, Astrophys. J. Lett. 885 (2019) L33.
  * [26] R. Shaikh and P. S. Joshi, JCAP 10 (2019) 064.
  * [27] F. H. Vincent, et al., Astron. Astrophys. 646 (2021) A37.
  * [28] J. Peng, M. Guo and X. H. Feng, Chin. Phys. C 45 (2021) 085103.
  * [29] Q. Gan, P. Wang, H. Wu and H. Yang, Phys. Rev. D 104 (2021) 044049.
  * [30] G. P. Li and K. J. He, JCAP 06 (2021) 037.
  * [31] S. Li, T. Mirzaev, A. A. Abdujabbarov, D. Malafarina, B. Ahmedov and W. B. Han, Phys. Rev. D 106 (2022) 084041.
  * [32] F. H. Vincent, S. E. Gralla, A. Lupsasca and M. Wielgus, Astron. Astrophys. 667 (2022) A170.
  * [33] T. Mirzaev, et al., Eur. Phys. J. C 83 (2023) 800.
  * [34] S. Staelens, D. R. Mayerson, F. Bacchini, B. Ripperda and L. Küchler, Phys. Rev. D 107 (2023) 124026.
  * [35] A. Uniyal, R. C. Pantig and A. Övgün, Phys. Dark Univ. 40 (2023) 101178.
  * [36] X. J. Wang, X. M. Kuang, Y. Meng, B. Wang and J. P. Wu, Phys. Rev. D 107 (2023) 124052.
  * [37] J.M. Bardeen, Non-singular general-relativistic gravitational collapse, in Proceedings of of International Conference GR5,Tbilisi, USSR (1968) 174.
  * [38] S. A. Hayward, Phys. Rev. Lett. 96 (2006) 031103.
  * [39] C. T. Cunningham, Astrophys. J. 202 (1975) 788.
  * [40] R. Gold, et al., Astrophys. J. 897 (2020) 148.
  * [41] V. Cardoso, A. S. Miranda, E. Berti, H. Witek and V. T. Zanchin, Phys. Rev. D 79 (2009) 064016.
  * [42] A. Cárdenas-Avendaño and A. Lupsasca, Phys. Rev. D 108 (2023) 064043.
  * [43] L. F. D. da Silva, F. S. N. Lobo, G. J. Olmo and D. Rubiera-Garcia, [arXiv:2307.06778 [gr-qc]].
  * [44] H. Paugnat, A. Lupsasca, F. Vincent and M. Wielgus, Astron. Astrophys. 668 (2022) A11.
  * [45] J. Peng, M. Guo and X. H. Feng, Phys. Rev. D 104 (2021) 124010.
  * [46] M. Guerrero, G. J. Olmo, D. Rubiera-Garcia and D. Sáez-Chillón Gómez, Phys. Rev. D 106 (2022) 044070.
  * [47] K. Akiyama et al. [Event Horizon Telescope], Astrophys. J. Lett. 930 (2022) L17.
  * [48] S. Vagnozzi, et al., Class. Quant. Grav. 40 (2023) 165007.
  * [49] A. E. Broderick, K. Salehi and B. Georgiev, [arXiv:2307.15120 [astro-ph.HE]].
  * [50] V. Cardoso, L. C. B. Crispino, C. F. B. Macedo, H. Okawa and P. Pani, Phys. Rev. D 90 (2014) 044069.
  * [51] V. Perlick and O. Y. Tsupko, Phys. Rept. 947 (2022) 1.