#  Discovering Love numbers through resonance excitation during extreme mass
ratio inspirals

Shani Avitan, Ram Brustein and Yotam Sherf
${}^{\textrm{\normalsize Department of Physics, Ben-Gurion University, Beer-
Sheva 84105, Israel}}$<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

General Relativity predicts that black holes do not possess an internal
structure and consequently cannot be excited. This leads to a specific
prediction about the waveform of gravitational waves which they emit during a
binary black hole inspiral and to the vanishing of their Love numbers.
However, if astrophysical black holes do possess an internal structure, their
Love numbers would no longer vanish, and they could be excited during an
inspiral by the transfer of orbital energy. This would affect the orbital
period and lead to an observable imprint on the emitted gravitational waves
waveform. The effect is enhanced if one of the binary companions is resonantly
excited. We discuss the conditions for resonant excitation of a hypothetical
internal structure of black holes and calculate the phase change of the
gravitational waves waveform that is induced due to such resonant excitation
during intermediate- and extreme-mass-ratio inspirals. We then relate the
phase change to the electric quadrupolar Love number of the larger companion,
which is resonantly excited by its smaller companion. We discuss the
statistical error on measuring the Love number by LISA and show that, because
of this phase change, the statistical error is small even for small values of
the Love number. Our results provide a strong indication that the Love number
could be detected by LISA with remarkable accuracy, much higher than what can
be achieved via tidal deformation effects. Our results further indicate that
resonant excitation of the central black hole during an extreme- or
intermediate-mass-ratio inspirals is the most promising effect for putting
bounds on, or detecting, non-vanishing tidal Love numbers of black holes.

##  1 Introduction

General Relativity predicts that black holes (BHs) do not possess an internal
structure. They are “bald” and can be characterize solely by their mass and
angular momentum [1]. Coalescing BHs radiate gravitational waves (GWs) which
are being detected by the LIGO and VIRGO observatories since September 2015
[2]. The calculational efforts for improving the accuracy of the general
relativistic (GR) predictions for the emitted GW waveform, could hopefully
provide an opportunity for testing the baldness of BHs. Particularly, the
inclusion of tidal interactions may allow us to probe the hypothetical
interior structure of the binary companions and quantitatively test the
predictions of GR [3, 4, 5, 6, 7, 8].

In spite of the increasing precision of ground-based detectors, their limited
frequency band enables observations of only a few cycles in the inspiral phase
of a binary-BH (BBH) coallesence event for a limited range of masses. The LISA
space detector [9, 10], whose design sensitivity is maximal in the mHZ region,
is expected to be able to detect and track many BBH coalescence events from
the early stages of the inspiral through the merger to the late post-merger
ringdown.

In GR, the interior of BHs is vacuum, except for a possibly singular core. But
is this their true description? A common expectation is that quantum effects
at the Planck length scale, or at a somewhat larger length scale, such as the
string length scale, will be sufficient to resolve all singularities. However,
there are strong indications to the contrary when it comes to the resolution
of BH singularities. First, a seemingly necessary condition for evading the
singularity theorems [11, 12] and the closely related “Buchdahl-like” bounds
[13, 14, 15, 16] is that the geometry is sourced by matter that has the
maximal negative radial pressure permitted by causality, $p_{r}=-\rho$, all
the way to the surface of the star [17]. Furthermore, if one also considers
the emitted Hawking radiation from such a quantum-regularized BH, one finds an
untenable violation of energy conservation: When the scale of resolution is
parametrically smaller than that of the Schwarzschild radius, the emitted
energy of Hawking particles will greatly exceed the original mass of the
collapsing matter [18, 19]. Thus, the tentative conclusion that we will adopt
in our following discussion, is that deviations from GR must extend throughout
the object’s interior, that is, horizon-scale deviations from GR.

The Love numbers of an object encode its response to an external tidal field.
These numbers could provide some limited information about the mass
distribution and the compactness of the object. The Love numbers determine the
possible corrections to the GW signal due to tidal interactions of the
companion in a binary system. The quadrupolar Love number $k_{2}$ identically
vanishes for GR BHs in four spacetime dimensions [20, 21, 22, 23, 24, 25, 26],
making it a key observable. Measuring non-zero values will indicate a
deviation from the GR predictions [27, 28, 29, 30, 31, 32]. If indeed horizon
scale deviations from the GR predictions occur, then the expectation is that
the Love numbers will be small, but not extremely small, suppressed only by
some additional pertrubative parameter that quantifies the strength of the
deviations. The reason for such expectation is that the Love numbers are
normalized such that they are order unity if all the dimensional scales are of
order of their radius [29, 30].

Previous studies have primarily focused on measuring the Love numbers using
tidal deformability, which constitutes a subleading correction to the emitted
GW waveform and enters at 5PN order compared to the dominant point-particle
term. Tidal-deformability effects are more pronounced at the late inspiral
phase. This makes the measurement of the Love number more challenging, since
other finite-size effects are also of similar magnitude, requiring the
construction of more accurate GW waveforms and detectors with better
sensitivity. [3, 33, 34, 35, 36].

For GR BHs the inspiral evolution is dominated by the point-particle GW
emission. For BHs which posses an internal structure, an interesting different
effect can dominate the evolution if the orbital frequency becomes comparable
to a characteristic frequency of some internal mode. In this case, this mode
is resonantly excited, resulting in a rapid energy transfer from the orbit to
the internal mode. The loss of orbital energy effectively advances the
inspiral evolution, bringing the companions to a closer and faster orbit. The
abrupt energy transfer changes significantly the emitted GW waveform compared
to the point particle waveform since it leads to an instantaneous phase jump
and a secondary subleading accumulated dephasing due to the differences in
orbital velocities. Such resonant energy transfer can only be realized when
the internal modes are non-relativistic. The reason is that the Keplerian
orbital frequency is much smaller than the relativistic frequency $c/R$ ($c$
is the speed of light and $R$ is the radius of the compact object) when the
two objects are far from each other.

Tidal resonant excitations were first discussed in the context of ordinary
polytropic stars [37], then, much later, for binaries with at least one of the
companions being a neutron star [38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48,
49, 50]. In these cases, the effect was related to the tidal Love numbers [51,
52, 53, 54, 55, 56]. However, as already mentioned, since the corrections
enter the GWs waveform at 5PN order, the effect becomes significant during the
late inspiral phase, where additional effects are also significant, making it
difficult to detect the Love number with high confidence.

More recent studies related to BH mimickers [57, 58, 59, 60], treat the BBH as
if they were horizon-less ultra-compact objects (UCOs). In [57, 59, 60], the
tidal field was exciting some additional spacetime modes of a hypothetical
spacetime structure outside the UCO. In [59], the resulting phase shift due to
the resonant excitation of these additional spacetime modes was related to the
tidal Love numbers, and the detectability of the quadrupolar Love number
$k_{2}$ using observations of ground-based GW detectors and the proposed
Einstein telescope was discussed. In [58], the detectability prospects of the
resonance effects were discussed, but without connecting the effect to the
tidal Love numbers. In this study, no evidence for resonance was found in the
observations of the first two runs of Advanced LIGO and Advanced Virgo.

Here, in contrast to previous studies, we discuss the tidal excitation of
hypothetical non-relativistic internal modes of the BH, relate the resulting
phase shift to the Love numbers and discuss the possible detectability of
$k_{2}$ in LISA observations of IMRIs and EMRIs. We find that this is, by far,
the most promising way to put bounds on, or detect Love numbers of
astrophysical BHs.

In the following, contrary to GR, we assume that astrophysical BHs do have an
internal structure that supports non-relativistic fluid modes. We keep the
calculations as model–independent as possible by expressing the model-
dependent quantities through the dimensionless tidal Love number. We follow
the discussion in [59, 29], to relate the resonance phase shift of the excited
modes to the quadrupolar tidal Love number $k_{2}$, and their relation to
internal modes frequencies of quantum black holes [61, 29, 30, 31] and recent
frozen star model results [62, 63, 64]. We estimate the statistical error in
the measurement of $k_{2}$ through resonance excitations during the inspiral
of slowly rotating EMRIs and IMRIs, using the design noise spectrum of LISA
[9, 65]. We find that the statistical error is small even for small values of
the Love number, providing a strong indication that the Love number could be
detected with impressive accuracy. We end with an explicit comparison between
the detection prospects of the Love numbers with tidal deformability and tidal
resonance, and conclude that resonance excitations are the most promising
effect for detecting the Love numbers.

##  2 Tidal-Resonance interaction

Here, we examine the tidal interaction in a binary system, focusing on the
central object that is subjected to the weak periodic tidal force exerted by
the smaller companion. Following the ideas presented in [38, 42, 54, 56] and
more recently in [30], we describe the response of the object to the tidal
force from the perspective of an asymptotic observer. The idea is that the
object possess a set of non-relativistic fluid modes which are driven by the
tidal force and can be therefore described as a collection of driven harmonic
oscillators.

The spectrum of the interior fluid modes depends on the radial number $n$ and
the angular numbers $l,m$, so their frequencies depend of the three numbers
$\omega_{nlm}$. Here, we are particularly interested in the dominant effect
which is due to the excitation of the $n=1$ mode by the quadrupolar tidal
field, so we focus on the case $l=m=2$ [42]. As for the other modes; the
spherically symmetric static $m=0$ mode cannot generate pressure gradients
that are needed for resonance excitaion and therefore is not relevant to our
discussion. The $m=1$ mode can be resonantly excited in the case that the
spin-orbit configuration is misaligned [42, 44]. Here, we restrict our
attention to spin vectors that are aligned with the orbital angular momentum.

The mode corresponding to $n,~{}l,~{}m=1,~{}2,~{}2$ is non-relativistic,
meaning that, as for neutron stars, $\omega_{122}$ is parametrically smaller
than $c/R$. The orbital frequency which determines the frequency of the
driving tidal force is determined by the Kepler law. It follows, as explained
in the Introduction, that as the smaller object gets closer to the central
object, the orbital and internal frequencies can match.

When the frequency of one of the interior modes of the central, larger,
object, matches the orbital frequency of the companion, it is resonantly
excited and efficiently absorbs energy from the orbital motion. The
instantaneous energy absorption increases the orbital velocity and shortens
the inspiral duration, thus leading to a phase difference in the emitted GW
waveform, when compared to the emitted waveform in the absence of a resonance.
To calculate the dephasing of the GW waveform, we adopt the derivation in [39,
43], resulting in the following phase evolution,

$\displaystyle\begin{cases}\Phi(t)=\Phi_{PP}(t)&t<t_{R}\\\
\Phi(t)=\Phi_{PP}(t+\Delta t)-\Delta\Phi&t>t_{R}+\Delta t,\end{cases}$ (1)

where $\Phi_{PP}(t)$ is the point particle phase, $t_{R}$ is the time at which
the resonance starts, $\Delta t$ is the resonance duration and $\Delta\Phi$ is
the instantaneous resonance phase difference, which in general depends on the
object’s properties as demonstrated below. The point particle phase
$\Phi_{PP}$, is independent, by definition, on the object’s composition. In
particular, it has the same value for a GR black hole and one endowed with an
internal excitation spectrum, such as the objects we are discussing.

Then, assuming that the resonance duration is short compared to the inspiral
duration and under adiabatic evolution, we arrive at the frequency domain
resonance phase [39, 43],

$\displaystyle\Phi(f)~{}=~{}\Phi(f)_{PP}+\Theta(f-f_{R})\left(\dfrac{f}{f_{R}}-1\right)\Delta\Phi_{Res}~{},$
(2)

where $f_{R}$ is the internal mode frequency which satisfies the resonance
condition $2\pi f_{R}=m\Omega$, $\Omega$ being the orbital angular velocity.
Resonance corrections to the phase $\Delta\Phi_{Res}$, are composed of two
terms; a dominant term that enters at 2.5PN order higher than the leading
order point-particle term and a subleading 4PN-higher contribution. The
dominant contribution, which is frequency independent and proportional to
$\Delta\Phi$, originates from the instantaneous energy absorption during
resonance. The subleading term, which is proportional to the frequency, is a
secular effect that increases towards the late stages of the inspiral.

###  2.1 The phase shift

Fluid perturbations of compact objects are described by the displacement
vector $\xi^{i}$, of a fluid element from its unperturbed position, which is
given by the orthonormal base decomposition,

$\displaystyle\xi^{i}~{}=~{}\sum_{n}a_{n}\xi^{i}_{n},$ (3)

$\xi_{n}$ being the normal displacement vectors, and $a_{n}$ are the
dimensionless displacement amplitudes.111We use relativistic units $G,c=1$. In
the presence of tidal forces, the fluid modes satisfy the damped-driven
harmonic oscillator equation [38, 44],

$\displaystyle\ddot{a}_{nlm}+2\gamma_{n}\dot{a}_{nlm}+\omega_{n}^{2}a_{nlm}~{}=~{}\mathcal{F}(t)_{nlm},$
(4)

where $\gamma_{n}=-\text{Im}~{}\omega_{n}$ is the damping rate of the mode.
The source of the damping and its precise magnitude are irrelevant for the
resulting resonant excitation and dephasing. So, $\gamma$ can be neglected
altogether (see below).

The external periodic force $\mathcal{F}(t)_{nlm}$ excites the $n$th mode
interior fluid mode is given by

$\displaystyle\mathcal{F}(t)_{nlm}~{}=~{}N_{lm}{\dfrac{\mathcal{E}_{l}Q_{nl}}{MR^{2}}e^{-im\phi(t)}}~{},$
(5)

where $M$ and $R$ are the mass and radius of the central object. The order
unity factor $N_{lm}$ is proportional to the Wigner function and is specified
below. The tidal field of the $l$ mode is denoted by $\mathcal{E}_{l}$, which
for the $l=2,m=\pm 2$ satisfies
$\mathcal{E}_{ij}x^{i}x^{j}=\mathcal{E}r^{2}Y_{2\pm 2}$. The mass moment of
the quadrupolar $n$th mode $Q_{n}$ is given by the overlap integral [39],

$\displaystyle Q_{n}=-\int d^{3}r\delta\rho_{n}r^{2}~{},$ (6)

where $\delta\rho_{n}$ is the corresponding energy density perturbation.

Next, we aim to find the instantaneous phase shift $\Delta\Phi$ and the
corresponding phase evolution in Eq. (1). We start by solving Eq. (4) for the
amplitudes $a_{n}$, which at resonance is given by [44],

$\displaystyle
a_{n}(t)~{}=~{}\left({\dfrac{\pi}{m\ddot{\phi}}}\right)^{1/2}\frac{\mathcal{F}(t)_{nlm}}{\gamma_{nl}-i\omega_{nl}}e^{-i\omega_{nl}t},$
(7)

where $\ddot{\phi}$ denotes the rate of change of the orbital frequency at
resonance. The transferred energy to the mode $nlm$ during the resonance is a
sum of kinetic and potential terms [38, 44],

$\displaystyle
E_{nlm}(t)~{}=~{}\left(\dfrac{1}{2}\dot{a}_{nlm}(t)^{2}+\dfrac{1}{2}\omega_{nl}^{2}a_{nlm}^{2}(t)\right)MR^{2}~{}.$
(8)

The total energy absorbed by the mode, neglecting $\gamma_{nl}$, is given by

$\displaystyle\Delta
E_{nlm}~{}=~{}~{}N_{lm}^{2}\dfrac{\pi}{4m\ddot{\phi}}\dfrac{(\mathcal{E}_{l}Q_{nl})^{2}}{MR^{2}}~{}.$
(9)

The resonance excitations lead to a phase shift, since the orbital energy
decreases as it excites the interior modes. Accordingly, the orbital velocity
increases and the inspiral duration decreases by a time $\Delta t$. To
estimate $\Delta t$, we follow [43]. The energy absorbed by the central
objects decreases the energy of the orbit by the same amount. In the absence
of resonance, such a decrease in energy can only occur by the emission of GW
and the time that it would take the orbit to emit GW with such energy $\Delta
t$ would be determined by the equality $\dot{E}_{GW}\Delta t=\Delta E_{nlm}$.
The rate of GW emission $\dot{E}_{GW}$ is, to a very good approximation, the
same rate as in the absence of resonance, which to leading order is given by
$\dot{E}_{GW}=\frac{32}{5}(\mathcal{M}_{c}~{}\Omega)^{10/3}$, with
$\mathcal{M}_{c}$ being the chirp mass. The resulting phase shift
$\Delta\Phi=m\Omega\Delta t$ is the following,

$\displaystyle\Delta\Phi_{nlm}~{}=~{}m\Omega\dfrac{\Delta
E_{nlm}}{\dot{E}_{GW}}=\frac{5}{32}~{}m\Omega~{}\dfrac{\Delta
E_{nlm}}{(\mathcal{M}_{c}~{}\Omega)^{10/3}}.$ (10)

For IMRIs or EMRIs $\mathcal{M}_{c}\approx M$ and $\dot{E}_{GW}\sim v^{10}$.

Using Eq. (9), we may calculate the phase shift induced by the leading order
quadrupolar mode $l=m=2$ [39, 59],

$\displaystyle\Delta\Phi_{n22}~{}=~{}\dfrac{25\pi}{1024q(1+q)}\dfrac{1}{R_{1}^{5}}\dfrac{|Q_{n2}|^{2}}{M_{1}\omega_{22}^{2}R_{1}^{2}}=\dfrac{25\pi}{2048q(1+q)}\dfrac{1}{R_{1}^{5}}\dfrac{|Q_{n2}|^{2}}{\Delta
E^{int}},$ (11)

where we used that $N_{22}=\sqrt{3/32}$. Here $q=M_{2}/M_{1}$ is the mass
ratio and $\Delta E^{int}=\frac{1}{2}M_{1}\omega_{22}^{2}R_{1}^{2}$ is the
internal energy of oscillations which is related to the energy stored in the
$n$th mode by $\Delta E^{int}=\sum\limits_{n}\Delta E_{n22}$, [54].

We wish to justify our estimate of $\Delta t$ using only $\dot{E}_{GW}$ and
neglecting other dissipation effects. In general, the time difference $\Delta
t$ should include all types of dissipation channels, mainly the dominant
dissipation due to tidal friction and the subleading tidal deformation.
However, the rate of work of tidal friction is given by [66, 67]
$\dot{E}_{TF}=\frac{1}{2}Q_{ij}\dot{\mathcal{E}}^{ij}\sim k_{2}v^{15}\nu/M$,
where $\nu$ is the kinematic viscosity giving rise to viscous dissipation. In
[67], it is demonstrated that, under reasonable assumptions, the contribution
of viscous dissipation is negligibly small compared to the leading order GW
emission and, therefore, can be ignored. For example, for cold Neutron stars,
considered to be highly viscous $\nu/M\approx 10^{-7}$, whereas for BHs
$\nu/M=1$ [68]. During the inspiral, when the orbital velocity is non-
relativistic the ratio of the different emission rates scales as
$\dot{E}_{TF}/\dot{E}_{GW}\sim v^{5}\ll 1$, which shows that the internal
dissipation effects can indeed be neglected.

##  3 Fluid-origin Love numbers

Here we follow [29, 30] to determine the relationship between the Love number
and the spectrum of internal fluid modes. We focus on the static tidal Love
number, ignoring dissipative effects.

Following [30] (see also [56, 54]), we wish to find the static response of the
object to an external tidal field. At low frequencies, away from resonance,
the amplitude in Eq. (7) reduces to

$\displaystyle a_{n}=\dfrac{\mathcal{E}Q_{n}}{M\omega_{n}^{2}R^{2}}~{}.$ (12)

Then, using the definition of the Love number, $k_{2}R^{5}=3Q/(2\mathcal{E})$,
we apply the normal mode decomposition identities $Q=\sum_{n}a_{n}Q_{n}$, and
$k_{2}=\sum_{n}a_{n}k_{2n}$, where the n$th$ mode Love number, which is
associated with the n$th$ mode quadrupolar moment, is given by

$\displaystyle k_{2n}R^{5}=\frac{3Q_{n}}{2\mathcal{E}}~{}.$ (13)

when substituting the explicit form of $a_{n}$ from Eq. (12), the Love number
becomes

$\displaystyle
k_{2}~{}=~{}\sum_{n}\dfrac{3}{2R^{5}}\dfrac{Q_{n}^{2}}{M\omega_{n}^{2}R^{2}}~{}.$
(14)

We now approximate $k_{2}$ by the first term in the sum in Eq. (14) relying on
a physically motivated assumption. The sum in Eq. (14) is dominated by the
fundamental $n=1$ mode. The justification is that the number of nodes in the
overlap integral in Eq. (6) increases as $n$ increases. It follows that the
contribution of $Q_{n}$ decreases as $n$ increases. Using the $l=2$-mode
excitation energy $\Delta E^{int}_{n}=\frac{1}{2}M\omega_{n2}^{2}R^{2}$, the
sum in Eq. (14) can be approximated as

$\displaystyle k_{2}~{}\simeq~{}\dfrac{3}{4R^{5}}\dfrac{Q_{1}^{2}}{\Delta
E^{int}_{1}}~{}.$ (15)

We now observe that a similar expression to the one in Eq. (15), appears in
Eq. (11) which determines the phase shift $\Delta\Phi_{122}$. This allows to
express $\Delta\Phi_{122}$ in terms of $k_{2}$,

$\displaystyle\Delta\Phi_{Res}~{}=~{}\dfrac{25\pi}{1536}\dfrac{k_{2}}{q(1+q)}~{}.$
(16)

We are interested in the case of small mass ratios, $q\lesssim 1/1000$ and a
small but not extremely small $k_{2}$, $k_{2}\lesssim 1/10$. Then we can
parameterize the resonance dephasing by

$\displaystyle\Delta\Phi_{Res}~{}\simeq~{}5\times\left(\dfrac{k_{2}}{10^{-1}}\right)\left(\dfrac{q}{10^{-3}}\right)^{-1}~{}.$
(17)

The resonance-induced dephasing is governed by the dimensionless tidal Love
number and the companion’s mass ratio. Generally, the detection threshold for
the instantaneous phase jump requires $\Delta\Phi_{Res}\gtrsim 1$ [69]. Thus,
for typical values of Love numbers $k_{2}\lesssim 10^{-1}$, it is more likely
to observe resonances for moderate to extreme mass-ratio binaries $10^{-3}\leq
q\leq 10^{-5}$.

We can also express $k_{2}$ in terms of the frequency
$\omega_{12}\equiv\omega_{2}$ of the $n=1$, $l=2$ mode. At resonance, from Eq.
(6), $Q\sim\Delta E^{int}$, where $\Delta
E^{int}=\frac{1}{2}M\omega_{2}^{2}R^{2}$ is the energy of the oscillating star
at resonance. Thus, on dimensional grounds, we get $Q\sim\Delta E^{int}R^{2}$.
For example, for a constant energy density perturbation $Q=\frac{3}{5}\Delta
E^{int}R^{2}$, while typical non-constant energy density profiles result in a
numerical prefactor $\lesssim 1$ [56] (see also [29]). Substituting the
expressions for $Q$ and $\Delta E^{int}$, we arrive at our final result for
the Love number

$\displaystyle k_{2}~{}\simeq~{}\mathcal{N}\omega_{2}^{2}R^{2}~{},$ (18)

where $\mathcal{N}$ is an order unity dimensionless number that depends on the
object’s energy density profile and contains the numerical factors in the
definition of the Love number [29]. We will use Eq. (18) to determine the
detectability of $k_{2}$ in the next section.

Remarkably, in [29], it is shown that the gravitational polarizability of
objects which possess a discrete spectrum of quantum mechanical energy levels
is similar to that of classical stars. This follows from the fact that the
wavelength of the oscillation is of order of the star radius. We shall refer
these objects as “quantum black holes” (QBHs) to mean the quantum state that
corresponds to a classical BH. The idea is justified on the grounds of the
Bohr correspondence principle, where at macroscopic excitations, expectation
values are replaced by classical observables. Therefore, an excited quantum
macroscopic object can be treated as a semi-classical oscillating fluid-like
object that satisfy Eq. (4). Using standard time-independent quantum
perturbation theory, the Love number of QBHs is given by [29, 30]

$\displaystyle
k_{2}\simeq\dfrac{3}{4R^{5}}\dfrac{|\langle\Psi_{0}|\hat{Q}|n=1,l=2\rangle|^{2}}{\Delta
E_{1}^{int}}~{}.$ (19)

where $\Psi_{0}$ is the QBH ground state, $\hat{Q}$ is the mass moment
operator that obeys the no-hair theorem;
$\langle\Psi_{0}|\hat{Q}|\Psi_{0}\rangle=0$. The definition of Eq. (15) is
restored by applying the Bohr correspondence principle and replacing
expectations values with classical observables,
$\langle\Psi_{0}|\hat{Q}|n,l=2\rangle\leftrightarrow Q_{n}$. In this form, Eq.
(19) can be treated in a similar way to the classical treatment of Eqs.
(15),(18), which eventually recovers the result
$k_{2}\simeq\mathcal{N}\omega_{2}^{2}R^{2}$. The result is valid for any
object of radius $R$, quantum or classical, which has a quadrupole internal
mode whose non-relativistic frequency is $\omega_{2}$.

##  4 Detectability

In this section, using the Fisher method, we give a quantitative estimation of
the statistical error in measuring the Love number. We discuss the prospects
for detection of a non-vanishing Love number with the future space LISA
detector and demonstrate that during the inspiral, it is more likely to detect
the Love number with resonances rather than tidal deformability.

We evaluate the detectability of the Love numbers through resonant excitations
with the planned space telescope LISA, which according to [9], could track and
observe moderate to extreme mass-ratio binaries from the early stages of the
inspiral and up to the merger with high SNR. Before addressing the precise
statistical analysis, we wish to emphasize that for most of the range of the
binary masses and spins and for Love numbers $k_{2}\lesssim 10^{-1}$, the
leading order 2.5PN resonance phase term is comparable to the other effects
entering at 2.5PN, such as the PP 2.5PN term and the leading order tidal
heating term. For smaller values of $q$, the resonance phase term becomes
significant. Since it is established that LISA can detect the other 2.5PN
effects, we expect that LISA could be able to detect the Love numbers with
high confidence.

To evaluate the statistical error, we employ the Fisher information method.
Assuming a signal $s(t)=h(t,\theta^{i})+n$, with the uncorrelated noise $n$, a
model signal $h(t,\theta^{i})$ with model parameters $\theta^{i}$. For high
SNR events, the posterior distribution takes the form

$\displaystyle p(\theta^{i}|s)\propto
e^{-\frac{1}{2}\Delta\theta^{i}\Delta\theta^{j}\Gamma_{ij}}~{}.$ (20)

where $\Gamma_{ij}$ is the Fisher matrix defined as

$\displaystyle\Gamma_{ij}~{}=~{}\left(\dfrac{\partial
h}{\partial\theta^{i}}\Big{|}\dfrac{\partial
h}{\partial\theta^{j}}\right)~{}.$ (21)

with the inner product defined by
$(h_{1}|h_{2})=4\text{Re}\int_{f_{min}}^{f_{max}}\frac{\tilde{h}_{1}(f)\tilde{h}^{*}_{2}(f)}{S_{n}(f)}df$,
and $S_{n}(f)$ is LISA’s design noise spectral density. We choose
$f_{max}=f_{ISCO}(\chi)$, where $f_{ISCO}$ is the orbital frequency at the
innermost stable circular orbit(ISCO) and $f_{min}=10^{-5}$Hz as the lowest
frequency in the LISA frequency band. The model parameters are
$\theta^{i}=(\ln\mathcal{A},\ln\mathcal{M}_{c},\eta,\Phi_{c},t_{c},\chi_{1},\chi_{2},k_{2})$,
where $\mathcal{A}$ is the amplitude, $\mathcal{M}_{c}$ is the chirp mass,
$\eta$ is the symmetric mass-ratio, $\Phi_{c}$ and $t_{c}$ are the phase and
time at coalescence, $\chi_{i}$ are the companions spin parameter and $k_{2}$
is the Love number given in Eq. (18). The statistical error in measuring
$k_{2}$ is related to the Fisher matrix,

$\displaystyle\sigma_{k_{2}}~{}=~{}\sqrt{\langle\Delta
k_{2}\rangle^{2}}~{}=~{}\sqrt{(\Gamma^{-1})_{k_{2}k_{2}}}$ (22)

We consider quasi-circular orbits and employ the analytical frequency domain
post-Newtonian approximation TaylorF2, which accurately describes the binary
evolution of the inspiral up to the ISCO [70, 71, 72]. The frequency domain GW
waveform describing the binary inspiral is of the form
$\tilde{h}(f,\theta_{i})~{}=~{}\mathcal{A}e^{i\Phi}$, where $\Phi$ is the
phase evolution in Eq. (2). From Eq. (18), for $q\ll 1$, the instantaneous
phase shift at resonance becomes

$\displaystyle\Delta\Phi_{Res}~{}\approx~{}\mathcal{N}~{}\dfrac{\omega_{2}^{2}R^{2}}{20q}.$
(23)

In our analysis we included correction terms up to 3PN order and neglected the
higher order tidal deformability terms that depend on the Love number and
enter at 5PN and 6PN order (See Sec. 4.1).

Figure 1: The solid blue lines correspond to a potential measurement of
$k_{2}$ for a given mass ratio $q$, with relative error
$\sigma_{k_{2}}/k_{2}=1/3$. The region above each solid line corresponds
potential measurement of $k_{2}$ with a relative error smaller than $1/3$. As
anticipated by Eq. (16), for a smaller mass ratio, the error on measuring a
specific $k_{2}$ is smaller and it is possible to measure smaller values of
$k_{2}$. The purple region describes the parameter space accessible to our
model for values of the spin parameter between 0 and 0.5, taking into account
Love-resonance-spin relation: $k_{2}\propto\omega_{ISCO}(\chi)R(\chi)$, such
that a given Love number corresponds to a specific resonance frequency and
spin parameter. The gray region describes the parameter space which is not
accessible to our model for these values of the spin parameters.

Additionally, since our model is valid only until the ISCO, the frequency
range $\omega_{2}>\omega_{ISCO}$ is not included in our analysis.
Consequently, it is beneficial to parameterize the oscillation frequencies in
terms of the ISCO frequency $\omega_{2}=\alpha\omega_{ISCO}$, where
$0<\alpha\leq 1$, and $\omega_{ISCO}(\chi)$ is spin-dependent. This also means
that resonance at the ISCO sets the maximal value of the Love number that can
be detected $k_{2}^{max}=\mathcal{N}\omega_{ISCO}^{2}R^{2}$.

We consider moderate to extreme mass-ratio binaries with
$q=[10^{-3},10^{-4},10^{-5}]$, where the central object mass is
$M_{1}=10^{6}M_{\odot}$, and small to moderate Kerr spin parameters
$\chi^{i}=[0,0.1,0.2,0.3,0.4,0.5]$, at a luminosity distance $D_{l}=2$Gpc. We
also average over the sky location parameters [72]. We assume equal spins
$\chi_{1}=\chi_{2}$ that are aligned with the orbital angular velocity vector.
For the model-dependent order unity coefficient $\mathcal{N}$, we use the
estimation derived in [29], and consider $\mathcal{N}\in[0.1,1]$.

In Fig.1, the purple region shows the analytical Love-resonance-spin relation
described in Eq. (18) that is determined by our model, where a given Love
number corresponds to a specific resonance frequency and a spin parameter.
This region describes the parameter space accessible to our model and is
independent of the detector properties. In our analysis, the largest
accessible $k_{2}$ is reached for $\mathcal{N}=1$, $\alpha=1$ and $\chi=0.5$,
resulting in $k_{2}^{max}\approx 0.159$, larger values are inaccessible to our
model. The gray region is the parameter space region that our model cannot
describe.

###  4.1 Comparison to Tidal-deformability

We now turn to estimate the relative magnitude of the resonance phase shift
effects compared to the magnitude of tidal deformation effects on the phase
evolution. To leading PN order, the tidal deformability contribution to the
phase for $q\ll 1$ takes the form $\Phi_{TD}(f)\sim k_{2}v^{5}/q$, where
$v=(\pi Mf)^{1/3}$ is the orbital velocity. The accumulated phase throughout
the inspiral is given by

$\displaystyle\Delta\Phi_{TD}~{}=~{}\int_{f_{min}}^{f_{ISCO}}f\dfrac{d^{2}\Phi_{TD}(f)}{df^{2}}df\sim\dfrac{k_{2}}{q}v_{ISCO}^{5}~{}.$
(24)

Figure 2: The figure displays the relative statistical errors in measuring
$k_{2}$ with resonance excitations and tidal deformation. For a given $k_{2}$
and $\chi$ we calculate $\sigma_{k_{2}}^{\text{res}}$ without tidal
deformation effects and $\sigma_{k_{2}}^{TD}$ without resonances. The results
show a preference for detecting $k_{2}$ with resonance effects. The preference
is more apparent for a smaller mass ratio $q$. The colored regions enclosed by
the solid and the dashed lines mark the additional parameter space that
resonances can probe compared to tidal deformation.

For a case for which the central object mass is $M_{1}=10^{6}M_{\odot}$ and
for small to moderate spin parameters, we find $v_{ISCO}^{5}\sim 0.01$.
Comparing to the instantaneous resonance phase jump Eq. (11),
$\Delta\Phi_{TD}/\Delta\Phi_{Res}\sim v_{ISCO}^{5}$. Therefore, we would
expect to have a larger error in the measurement of the Love number relying on
tidal deformability.

We calculated the statistical error in measuring the Love number through tidal
deformability and compared it to a measurement via resonance effects and found
that the previous estimate is indeed correct. We repeated the statistical
evaluation performed above, excluding resonance effects and including the
leading tidal deformation terms entering the phase at 5PN and 6PN order [73,
74, 35].

The results of the calculation of the ratio of the relative errors in
measuring the Love numbers, denoted by
$\sigma^{Res}_{{k_{2}}}/\sigma^{TD}_{{k_{2}}}$ for different spin parameters
$0\leq\chi\leq 0.5$ are presented in Fig. 2.

##  5 Summary and conclusion

The future measurement of GWs produced during BBH inspirals by the planned GW
detector LISA will present an unprecedented opportunity to test GR.
Hypothetical tidal interactions between the inspiraling objects would affect
the waveform of the emitted GWs in a way that could only be possible if
astrophysical BHs were actually ultra-compact objects possesing an internal
structure rather than the structureless objects predicted by GR.

We discussed how the resonant excitation of the hypothetical non-relativistic
interior modes of astrophysical BHs changes the phase of the emitted GW
waveform when compared to the phase predicted by GR. The non-relativistic
nature of the modes was crucial to the possibility of resonantly exciting
them, because in this case they could be excited when the two objects are
still far apart. In this case, the resonance occurs a long time before the
ISCO is reached and leads to a significant dephasing. We find that regardless
the specific details of the primary’s interior composition, the phase shift is
governed by a single intrinsic quantity - the dimensionless tidal Love number
$k_{2}$.

We evaluated the statistical error in measuring the Love number $k_{2}$ by
LISA using the resonance effect. We concluded that the smallness of the
resulting statistical error indicates that $k_{2}$ could actually be detected
by LISA with impressive accuracy by observing intermediate and extreme mass-
ratio inspirals. We compared the statistical error for detection of the Love
number relying on tidal deformation effects with the error when using
resonance effects and concluded that prospects of measuring $k_{2}$ using
resonance effects are much better. The results reveal additional sensitivity-
enhancement factors whose origin is the Love-resonance-spin relation. First,
the statistical error in measuring the Love number reduces for BHs with higher
spin, because for such BHs, the inspiral duration is longer. Second, the
statistical error in measuring the Love number reduces if the inspiral
includes a range of higher orbital velocities, which could lead to excitation
of higher internal frequencies, which, in turn, correspond to the BH having a
larger Love number.

Our conclusion is that the effects of resonant excitation of astrophysical BHs
during intermediate and extreme mass-ratio inspirals provide the best
opportunity for putting bounds on, or detecting, the tidal Love number of
astrophysical BHs and thus providing evidence of physics beyond GR.
Nevertheless, we stress that the results of our statistical analysis should be
viewed as preliminary estimates for the detection prospects. A comprehensive
statistical treatment requires more accurate waveform modeling and should
consider LISA’s ability to track and discriminate several EMRIs simultaneously
[9].

Our analysis is based on a general theoretical framework which only requires
the existence of a set of non-relativistic internal modes, and does not
require specifying the detailed properties of the central object. The entire
dependence on the interior composition is parameterized in terms of the
dimensionless tidal Love numbers. Therefore our results can be applied to a
wide range of ultra-compact objects or BHs mimickers.

## Acknowledgments

We thank Vitor Cardoso, Tanja Hinderer and Ely Kovetz for useful discussions.
The research is supported by the German Research Foundation through a German-
Israeli Project Cooperation (DIP) grant “Holography and the Swampland” and by
VATAT (Israel planning and budgeting committee) grant for supporting
theoretical high energy physics.

## References

  * [1] J. D. Bekenstein, “Nonexistence of baryon number for static black holes,” Phys. Rev. D 5 (1972), 1239-1246
  * [2] B. P. Abbott et al. [LIGO Scientific and Virgo], “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116 (2016) no.6, 061102 [arXiv:1602.03837 [gr-qc]].
  * [3] E. E. Flanagan and T. Hinderer, “Constraining neutron star tidal Love numbers with gravitational wave detectors,” Phys. Rev. D 77 (2008), 021502 [arXiv:0709.1915 [astro-ph]].
  * [4] T. Hinderer, B. D. Lackey, R. N. Lang and J. S. Read, “Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral,” Phys. Rev. D 81 (2010), 123016 [arXiv:0911.3535 [astro-ph.HE]].
  * [5] S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger and C. M. Biwer, “Tidal Deformabilities and Radii of Neutron Stars from the Observation of GW170817,” Phys. Rev. Lett. 121 (2018) no.9, 091102 [erratum: Phys. Rev. Lett. 121 (2018) no.25, 259902] [arXiv:1804.08583 [astro-ph.HE]].
  * [6] A. Maselli, V. Cardoso, V. Ferrari, L. Gualtieri and P. Pani, “Equation-of-state-independent relations in neutron stars,” Phys. Rev. D 88 (2013) no.2, 023007 [arXiv:1304.2052 [gr-qc]].
  * [7] S. S. Yazadjiev, D. D. Doneva and K. D. Kokkotas, “Tidal Love numbers of neutron stars in $f(R)$ gravity,” Eur. Phys. J. C 78 (2018) no.10, 818 [arXiv:1803.09534 [gr-qc]].
  * [8] C. Pacilio, A. Maselli, M. Fasano and P. Pani, “Ranking Love Numbers for the Neutron Star Equation of State: The Need for Third-Generation Detectors,” Phys. Rev. Lett. 128 (2022) no.10, 101101 [arXiv:2104.10035 [gr-qc]].
  * [9] P. Amaro-Seoane et al. [LISA], “Laser Interferometer Space Antenna,” [arXiv:1702.00786 [astro-ph.IM]].
  * [10] K. G. Arun et al. [LISA], “New horizons for fundamental physics with LISA,” Living Rev. Rel. 25 (2022) no.1, 4 [arXiv:2205.01597 [gr-qc]].
  * [11] R. Penrose, “Gravitational Collapse and Space-Time Singularities,” Phys. Rev. Lett. 14, 57 (1965).
  * [12] S. W. Hawking and R. Penrose, “The singularities of gravitational collapse and cosmology,” Proc. R. Soc. Lond. A 314, 529 (1970).
  * [13] H. Buchdahl, “General Relativistic Fluid Spheres,” Phys. Rev. 116, 1027 (1959).
  * [14] S. Chandrasekhar “Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity,” Phys. Rev. Lett. 12, 114 (1964).
  * [15] S. Chandrasekhar, “The Dynamical Instability of Gaseous Masses Approaching the Schwarzschild Limit in General Relativity,” Astrophys. J. 140, 417 (1964).
  * [16] H. Bondi, “Massive spheres in general relativity,” Proc. Roy. Soc. Lond. A 282, 303 (1964).
  * [17] R. Brustein and A. J. M. Medved, “Resisting collapse: How matter inside a black hole can withstand gravity,” Phys. Rev. D 99, no.6, 064019 (2019) [arXiv:1805.11667 [hep-th]].
  * [18] V. P. Frolov and A. Zelnikov, “Quantum radiation from an evaporating nonsingular black hole,” Phys. Rev. D 95, no.12, 124028 (2017) [arXiv:1704.03043 [hep-th]].
  * [19] R. Carballo-Rubio, F. Di Filippo, S. Liberati, C. Pacilio and M. Visser, “On the viability of regular black holes,” JHEP 07, 023 (2018) [arXiv:1805.02675 [gr-qc]].
  * [20] T. Damour and A. Nagar, “Relativistic tidal properties of neutron stars,” Phys. Rev. D 80 (2009), 084035 [arXiv:0906.0096 [gr-qc]].
  * [21] T. Binnington and E. Poisson, “Relativistic theory of tidal Love numbers,” Phys. Rev. D 80 (2009), 084018 [arXiv:0906.1366 [gr-qc]].
  * [22] H. S. Chia, “Tidal deformation and dissipation of rotating black holes,” Phys. Rev. D 104 (2021) no.2, 024013 [arXiv:2010.07300 [gr-qc]].
  * [23] P. Charalambous, S. Dubovsky and M. M. Ivanov, “Hidden Symmetry of Vanishing Love Numbers,” Phys. Rev. Lett. 127 (2021) no.10, 101101 [arXiv:2103.01234 [hep-th]].
  * [24] P. Charalambous, S. Dubovsky and M. M. Ivanov, “On the Vanishing of Love Numbers for Kerr Black Holes,” JHEP 05 (2021), 038 [arXiv:2102.08917 [hep-th]].
  * [25] L. Hui, A. Joyce, R. Penco, L. Santoni and A. R. Solomon, “Ladder symmetries of black holes. Implications for love numbers and no-hair theorems,” JCAP 01 (2022) no.01, 032 [arXiv:2105.01069 [hep-th]].
  * [26] J. Ben Achour, E. R. Livine, S. Mukohyama and J. P. Uzan, “Hidden symmetry of the static response of black holes: applications to Love numbers,” JHEP 07 (2022), 112 [arXiv:2202.12828 [gr-qc]].
  * [27] A. Maselli, P. Pani, V. Cardoso, T. Abdelsalhin, L. Gualtieri and V. Ferrari, “Probing Planckian corrections at the horizon scale with LISA binaries,” Phys. Rev. Lett. 120 (2018) no.8, 081101 [arXiv:1703.10612 [gr-qc]].
  * [28] V. Cardoso, E. Franzin, A. Maselli, P. Pani and G. Raposo, “Testing strong-field gravity with tidal Love numbers,” Phys. Rev. D 95 (2017) no.8, 084014 doi:10.1103/PhysRevD.95.084014 [arXiv:1701.01116 [gr-qc]].
  * [29] R. Brustein and Y. Sherf, “Quantum Love numbers,” Phys. Rev. D 105, no.2, 024043 (2022) [arXiv:2008.02738 [gr-qc]].
  * [30] R. Brustein and Y. Sherf, “Classical Love number for quantum black holes,” Phys. Rev. D 105 (2022) no.2, 024044 [arXiv:2104.06013 [gr-qc]].
  * [31] R. Brustein and Y. Sherf, “Tidal deformation of quantum black holes,” Int. J. Mod. Phys. D 30 (2021) no.14, 2142011
  * [32] S. Datta, “Probing horizon scale quantum effects with Love,” Class. Quant. Grav. 39 (2022) no.22, 225016 [arXiv:2107.07258 [gr-qc]].
  * [33] T. Damour, A. Nagar and L. Villain, “Measurability of the tidal polarizability of neutron stars in late-inspiral gravitational-wave signals,” Phys. Rev. D 85 (2012), 123007 [arXiv:1203.4352 [gr-qc]].
  * [34] B. P. Abbott et al. [LIGO Scientific and Virgo], “GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral,” Phys. Rev. Lett. 119 (2017) no.16, 161101 [arXiv:1710.05832 [gr-qc]].
  * [35] N. Sennett, T. Hinderer, J. Steinhoff, A. Buonanno and S. Ossokine, “Distinguishing Boson Stars from Black Holes and Neutron Stars from Tidal Interactions in Inspiraling Binary Systems,” Phys. Rev. D 96 (2017) no.2, 024002 [arXiv:1704.08651 [gr-qc]].
  * [36] G. A. Piovano, A. Maselli and P. Pani, “Constraining the tidal deformability of supermassive objects with extreme mass ratio inspirals and semianalytical frequency-domain waveforms,” Phys. Rev. D 107 (2023) no.2, 024021 [arXiv:2207.07452 [gr-qc]].
  * [37] T. G. Cowling, M. A., D. Phil, Monthly Notices of the Royal Astronomical Society, Volume 101, Issue 8 (1941)
  * [38] D. Lai, “Resonant oscillations and tidal heating in coalescing binary neutron stars,” Mon. Not. Roy. Astron. Soc. 270 (1994), 611 [arXiv:astro-ph/9404062 [astro-ph]].
  * [39] A. Reisenegger and P. Goldreich, Astrophys. J. 426, 688 (1994).
  * [40] M. Shibata, Effects of Tidal Resonances in Coalescing Compact Binary Systems, Progress of Theoretical Physics, Volume 91, Issue 5 (1994)
  * [41] K. D. Kokkotas and G. Schaefer, “Tidal and tidal resonant effects in coalescing binaries,” Mon. Not. Roy. Astron. Soc. 275 (1995), 301 [arXiv:gr-qc/9502034 [gr-qc]].
  * [42] W. C. G. Ho and D. Lai, “Resonant tidal excitations of rotating neutron stars in coalescing binaries,” Mon. Not. Roy. Astron. Soc. 308 (1999), 153 [arXiv:astro-ph/9812116 [astro-ph]].
  * [43] E. E. Flanagan and E. Racine, “Gravitomagnetic resonant excitation of Rossby modes in coalescing neutron star binaries,” Phys. Rev. D 75 (2007), 044001 [arXiv:gr-qc/0601029 [gr-qc]].
  * [44] D. Lai and Y. Wu, “Resonant Tidal Excitations of Inertial Modes in Coalescing Neutron Star Binaries,” Phys. Rev. D 74 (2006), 024007 [arXiv:astro-ph/0604163 [astro-ph]].
  * [45] T. Hinderer, A. Taracchini, F. Foucart, A. Buonanno, J. Steinhoff, M. Duez, L. E. Kidder, H. P. Pfeiffer, M. A. Scheel and B. Szilagyi, et al. “Effects of neutron-star dynamic tides on gravitational waveforms within the effective-one-body approach,” Phys. Rev. Lett. 116 (2016) no.18, 181101 [arXiv:1602.00599 [gr-qc]].
  * [46] C. P. L. Berry, R. H. Cole, P. Cañizares and J. R. Gair, “Importance of transient resonances in extreme-mass-ratio inspirals,” Phys. Rev. D 94 (2016) no.12, 124042 [arXiv:1608.08951 [gr-qc]].
  * [47] H. Yu and N. N. Weinberg, “Resonant tidal excitation of superfluid neutron stars in coalescing binaries,” Mon. Not. Roy. Astron. Soc. 464 (2017) no.3, 2622-2637 [arXiv:1610.00745 [astro-ph.HE]].
  * [48] G. Pratten, P. Schmidt and T. Hinderer, “Gravitational-Wave Asteroseismology with Fundamental Modes from Compact Binary Inspirals,” Nature Commun. 11 (2020) no.1, 2553 [arXiv:1905.00817 [gr-qc]].
  * [49] B. Bonga, H. Yang and S. A. Hughes, “Tidal resonance in extreme mass-ratio inspirals,” Phys. Rev. Lett. 123 (2019) no.10, 101103 [arXiv:1905.00030 [gr-qc]].
  * [50] P. Gupta, B. Bonga, A. J. K. Chua and T. Tanaka, “Importance of tidal resonances in extreme-mass-ratio inspirals,” Phys. Rev. D 104 (2021) no.4, 044056 [arXiv:2104.03422 [gr-qc]].
  * [51] P. K. Gupta, J. Steinhoff and T. Hinderer, “Effect of dynamical gravitomagnetic tides on measurability of tidal parameters for binary neutron stars using gravitational waves,” [arXiv:2302.11274 [gr-qc]].
  * [52] P. K. Gupta, A. Puecher, P. T. H. Pang, J. Janquart, G. Koekoek and C. Broeck Van Den, “Determining the equation of state of neutron stars with Einstein Telescope using tidal effects and r-mode excitations from a population of binary inspirals,” [arXiv:2205.01182 [gr-qc]].
  * [53] J. Steinhoff, T. Hinderer, A. Buonanno and A. Taracchini, “Dynamical Tides in General Relativity: Effective Action and Effective-One-Body Hamiltonian,” Phys. Rev. D 94 (2016) no.10, 104028 [arXiv:1608.01907 [gr-qc]].
  * [54] S. Chakrabarti, T. Delsate and J. Steinhoff, “Effective action and linear response of compact objects in Newtonian gravity,” Phys. Rev. D 88, 084038 (2013) [arXiv:1306.5820 [gr-qc]].
  * [55] S. Chakrabarti, T. Delsate and J. Steinhoff, “New perspectives on neutron star and black hole spectroscopy and dynamic tides,” [arXiv:1304.2228 [gr-qc]].
  * [56] N. Andersson and P. Pnigouras, “Exploring the effective tidal deformability of neutron stars,” Phys. Rev. D 101, no.8, 083001 (2020) [arXiv:1906.08982 [astro-ph.HE]].
  * [57] V. Cardoso, A. del Rio and M. Kimura, “Distinguishing black holes from horizonless objects through the excitation of resonances during inspiral,” Phys. Rev. D 100 (2019), 084046 [erratum: Phys. Rev. D 101 (2020) no.6, 069902] [arXiv:1907.01561 [gr-qc]].
  * [58] Y. Asali, P. T. H. Pang, A. Samajdar and C. Van Den Broeck, “Probing resonant excitations in exotic compact objects via gravitational waves,” Phys. Rev. D 102, no.2, 024016 (2020) [arXiv:2004.05128 [gr-qc]].
  * [59] K. Fransen, G. Koekoek, R. Tielemans and B. Vercnocke, Phys. Rev. D 104, no.4, 044044 (2021) [arXiv:2005.12286 [gr-qc]].
  * [60] E. Maggio, M. van de Meent and P. Pani, “Extreme mass-ratio inspirals around a spinning horizonless compact object,” Phys. Rev. D 104 (2021) no.10, 104026 [arXiv:2106.07195 [gr-qc]].
  * [61] R. Brustein, A. J. M. Medved and K. Yagi, “When black holes collide: Probing the interior composition by the spectrum of ringdown modes and emitted gravitational waves,” Phys. Rev. D 96, no.6, 064033 (2017) [arXiv:1704.05789 [gr-qc]].
  * [62] R. Brustein, A. J. M. Medved and T. Simhon, “Black holes as frozen stars,” Phys. Rev. D 105, no.2, 024019 (2022) [arXiv:2109.10017 [gr-qc]].
  * [63] R. Brustein, A. J. M. Medved, T. Shindelman and T. Simhon, “Black holes as frozen stars: Regular interior geometry,” [arXiv:2301.09712 [gr-qc]].
  * [64] R. Brustein, A. J. M. Medved and T. Shindelman, ‘Defrosting frozen stars: spectrum of internal fluid modes,” [arXiv:2304.04984 [gr-qc]].
  * [65] T. Robson, N. J. Cornish and C. Liu, “The construction and use of LISA sensitivity curves,” Class. Quant. Grav. 36 (2019) no.10, 105011 [arXiv:1803.01944 [astro-ph.HE]].
  * [66] E. Poisson, “Tidal interaction of black holes and Newtonian viscous bodies,” Phys. Rev. D 80 (2009), 064029 [arXiv:0907.0874 [gr-qc]].
  * [67] Y. Sherf, “Tidal-heating and viscous dissipation correspondence in black holes and viscous compact objects,” Phys. Rev. D 103 (2021) no.10, 104003 [arXiv:2104.03766 [gr-qc]].
  * [68] K. S. Thorne, R. H. Price and D. A. Macdonald, “Black Holes: The Membrane Paradigm,” New Haven: Yale University Press, 1986.
  * [69] L. Lindblom, B. J. Owen and D. A. Brown, “Model Waveform Accuracy Standards for Gravitational Wave Data Analysis,” Phys. Rev. D 78 (2008), 124020 [arXiv:0809.3844 [gr-qc]].
  * [70] T. Damour, B. R. Iyer and B. S. Sathyaprakash, “A Comparison of search templates for gravitational waves from binary inspiral,” Phys. Rev. D 63 (2001), 044023 [erratum: Phys. Rev. D 72 (2005), 029902] [arXiv:gr-qc/0010009 [gr-qc]].
  * [71] K. G. Arun, B. R. Iyer, B. S. Sathyaprakash and P. A. Sundararajan, “Parameter estimation of inspiralling compact binaries using 3.5 post-Newtonian gravitational wave phasing: The Non-spinning case,” Phys. Rev. D 71 (2005), 084008 [erratum: Phys. Rev. D 72 (2005), 069903] [arXiv:gr-qc/0411146 [gr-qc]].
  * [72] A. Buonanno, B. Iyer, E. Ochsner, Y. Pan and B. S. Sathyaprakash, “Comparison of post-Newtonian templates for compact binary inspiral signals in gravitational-wave detectors,” Phys. Rev. D 80 (2009), 084043 [arXiv:0907.0700 [gr-qc]].
  * [73] D. Bini, T. Damour and G. Faye, “Effective action approach to higher-order relativistic tidal interactions in binary systems and their effective one body description,” Phys. Rev. D 85 (2012), 124034 [arXiv:1202.3565 [gr-qc]].
  * [74] K. Hotokezaka, K. Kyutoku, Y. i. Sekiguchi and M. Shibata, “Measurability of the tidal deformability by gravitational waves from coalescing binary neutron stars,” Phys. Rev. D 93 (2016) no.6, 064082 [arXiv:1603.01286 [gr-qc]].


# Dynamic Computerized Tomography using Inexact Models and Motion Estimation

Gesa Sarnighausen Institute for Numerical and Applied Mathematics, University
of Göttingen, Germany (g.sarnighausen@math.uni-goettingen.de). Anne Wald
Institute for Numerical and Applied Mathematics, University of Göttingen,
Germany<EMAIL_ADDRESS>Alexander Meaney Department of
Mathematics and Statistics, University of Helsinki, Finland
(alexander.meaney@helsinki.fi).

Abstract. Reconstructing a dynamic object with affine motion in computerized
tomography (CT) leads to motion artifacts if the motion is not taken into
account. In most cases, the actual motion is neither known nor can be
determined easily. As a consequence, the respective model that describes CT is
incomplete. The iterative RESESOP-Kaczmarz method can - under certain
conditions and by exploiting the modeling error - reconstruct dynamic objects
at different time points even if the exact motion is unknown. However, the
method is very time-consuming. To speed the reconstruction process up and
obtain better results, we combine the following three steps:

  1. 1.

RESESOP-Kacmarz with only a few iterations is implemented to reconstruct the
object at different time points.

  2. 2.

The motion is estimated via landmark detection, e.g. using deep learning.

  3. 3.

The estimated motion is integrated into the reconstruction process, allowing
the use of dynamic filtered backprojection.

We give a short review of all methods involved and present numerical results
as a proof of principle.

Keywords: dynamic computerized tomography, affine motion, regularization,
sequential subspace optimization, landmark detection, dynamic filtered
backprojection, hybrid reconstruction

##  1 Introduction

Computerized tomography (CT) is an imaging technique that exploits the
characteristic attenuation of X-radiation by different materials: By
reconstructing the density inside an unknown object from measurements of the
intensity loss, it is possible to visualize the space-dependent density of the
object. CT thus requires the stable solution of an ill-posed inverse problem,
see [23].

The Radon operator $R$ \- the forward operator of computerized tomography -
explains the relation between the scanned object’s density $f$ and the
measured data $g$. This is usually formulated as an operator equation

$\displaystyle Rf=g,$

where the impact $g$ is called a sinogram, which is measured by the CT
scanner. The inverse problem thus consists in computing the cause $f$ from the
knowledge of $R$ and $g$. The most popular method to solve this problem is
filtered backprojection (FBP), which works well for static objects, see [23,
24, 14].

However, the scanned object may move or deform during the scanning process
since it takes time to scan the object from all necessary angles. An overview
about dynamic reconstruction methods can be found in [4].

Let the map $\Gamma:\Omega\to\Omega$ on the domain $\Omega$ of $f$ describe
the motion. Then, instead of $R$, we consider the dynamic Radon operator
$R^{\Gamma}:=R\circ\Gamma$. If the motion is not taken into account, e.g. in
FBP, this typically leads to severe motion artifacts in the reconstruction. In
most applications, the motion is not a quantity of interest, such that one
usually aims at reconstructing a specific state of the object.

In practice, compensating for the motion can be important in medical imaging,
e.g. since the patient is breathing and the heart is beating. On the nano
scale (nano-CT), manufacturing tolerances, small vibrations and object drift
lead to a relative motion between the object and the scanner [13]. If the
motion is known, it is possible to incorporate it into the reconstruction,
see, e.g., [10, 12, 16] and the references therein. Otherwise, one either has
to estimate the motion first, for example by tracking markers, and use the
dynamic model for the inversion as in [11], or by correcting the motion in the
measured sinogram and use the standard model as in [19], or both the motion
and the object have to be reconstructed jointly [2, 7, 6], which is
particularly demanding if there is only a single projection per time step
available. The problem gets even more involved if only sparse data is
available [5].

In this article, we consider the standard CT setup in two dimensions with a
parallel beam geometry. The angle of the tomograph and, by consequence, the
position of the X-rays source, is given by $\varphi\in[0,2\pi]$ and the offset
of an X-ray beam by $s\in[-1,1]$. The goal is to reconstruct a fixed state of
the object from measurements affected by a rigid body motion. We combine two
ways of including the motion in the model:

  1. 1.

Directly, by using the dynamic model

$g^{\Gamma}(\varphi,s)=R^{\Gamma}f(\varphi,s)\coloneqq\int_{\mathbb{R}^{2}}f(\Gamma_{\varphi}x)\delta(s-x^{T}\theta(\varphi))\,\mathrm{d}x,$

for $s\in\mathbb{R},\varphi\in[0,2\pi]$,
$\theta(\varphi)=(\cos(\varphi),\sin(\varphi))^{\top}$ and the Dirac’s delta
distribution $\delta$.

  2. 2.

Indirectly, by interpreting the static operator as an inexact version of the
dynamic operator and adding the inexactness $\eta$ to the static model:

$g^{\Gamma}(\varphi,s)=Rf(\varphi,s)+\eta(\varphi,s)\coloneqq\int_{\mathbb{R}^{2}}f(x)\delta(s-x^{T}\theta(\varphi))\,\mathrm{d}x+\eta(\varphi,s).$

Figure 1: Combined Method

A class of adapted RESESOP-Kaczmarz algorithms [3], RESESOP being short for
REgularizing SEquential Subspace OPtimization, can be used for the second
approach. First, our dynamic CT problem is considered a semi-discrete inverse
problem $R^{\Gamma}_{(\varphi,s)}f=g(\varphi,s)$ for all discrete tuples
$(\varphi,s)$. The motion is interpreted as an inexactness in the forward
model with respect to the static model, which needs to be estimated and is
then compensated when computing the reconstruction. Hence, the subproblems are
combined similarly as in Kaczmarz’ iteration, but with an adapted RESESOP
projection step instead of the usual projection onto the solution set of the
respective subproblem.
Instead of knowing the motion, it is thus sufficient to estimate the
inexactness of the model. The estimation of the inexactness is still an open
task, but has certain benefits: the estimate needs not to be exact, but a
small overestimation often leads to reasonable results, it can incorporate
temporal and local information [3], and the respective reconstructions seem to
yield a good basis for data-driven image processing [20]. However, RESESOP-
Kaczmarz is relatively slow due to the Kaczmarz loops which can not be
parallelized here.

To obtain a more efficient reconstruction technique under the assumption that
an estimate of the inexactness is available, we propose a hybrid approach
consisting of three steps for dynamic CT as illustrated in Figure 1:

  1. 1.

RESESOP-Kaczmarz

First, the RESESOP-Kaczmarz algorithm is used with a very small amount of
iterations and in a small resolution to calculate rather rough reconstructions
of certain states of the object.

  2. 2.

Landmark detection

In each reconstruction, we identify landmarks of the object and track them to
get an estimate of the actual motion. This task can be done by applying a
trained deep neural network, but also manually.

  3. 3.

Dynamic Filtered Backprojection

As the last step, the dynamic filtered backprojection algorithm is used with
the estimated motion to get a better reconstruction of the object in high
resolution.

For a numerical proof of concept, we consider a simple directed rigid body
motion of the motion, i.e., the object shifts in a fixed direction with a
constant speed during data acquisition.

The article is structured as follows. Section 2 explains the three methods we
used. We start by presenting the RESESOP-Kaczmarz method in Section 2.1.
Section 2.2 addresses the derivation of the underlying motion from landmark
detection via deep learning or by hand. and in Section 2.3, the dynamic
filtered backprojection algorithm for affine deformations is reviewed. Section
3 contains the numerical results of the combined method introduced in this
article for generated and real measured data.

##  2 Methods

###  2.1 RESESOP-Kaczmarz

The regularizing sequential subspace optimization (RESESOP) method is a
regularizing iterative method to solve linear or nonlinear inverse problems in
Hilbert and Banach spaces [26, 25, 30] and has been extended to problems with
an inexact forward operator [3]. It can especially be used for dynamic inverse
problems if the static operator is considered the inexact version of the
dynamic forward operator. The method combines the projection step from RESESOP
with the semi-discrete setting from the Kaczmarz method and can be considered
a kind of regularized gradient method with several search directions and a
step size regulation by projecting onto suitable subsets of the source space.

Let ($X,\lVert\cdot\rVert_{X}$) and ($Y,\lVert\cdot\rVert_{Y}$) be real
Hilbert spaces with their corresponding norms and inner products and
$\mathcal{L}(X,Y)$ be the space of bounded linear operators from $X$ to $Y$.
The inverse problem we want to solve is then

(1) $\displaystyle\mathcal{A}f=g,\quad\mathcal{A}\in\mathcal{L}(X,Y)$

for measured data $g\in Y$ and the solution set
$\mathcal{M}_{\mathcal{A},g}\coloneqq\\{f\in X:\mathcal{A}f=g\\}$.

##### Sequential subspace optimization (SESOP)

The SESOP method iteratively solves inverse problems by projecting onto
suitable convex affine subspaces in a continuous setting with noise-free data
and an exact forward operator. To construct such subspaces, it is first
necessary to define hyperplanes, half-spaces and stripes.

######  Definition 1.

(Hyperplane, half-space, stripe)

For $u\in X\setminus\\{0\\}$ and $\alpha\in\mathbb{R}$ the set

$\displaystyle H(u,\alpha)\coloneqq\\{f\in X:\langle u,f\rangle_{X}=\alpha\\}$

is called a hyperplane in $X$ and the set

$\displaystyle H_{\leq}(u,\alpha)\coloneqq\\{f\in X:\langle
u,f\rangle_{X}\leq\alpha\\}$

is called a half-space. Analogously the half-spaces $H_{\geq}(u,\alpha)$,
$H_{<}(u,\alpha)$ and $H_{>}(u,\alpha)$ are defined.

Now let additionally be $\xi\geq 0$. Then a stripe can be defined as the set

$\displaystyle H(u,\alpha,\xi)\coloneqq\\{f\in X:\lvert\langle
u,f\rangle_{X}-\alpha\rvert\leq\xi\\}.$

Each stripe is bounded by the hyperplanes $H(u,\alpha+\xi)$ and
$H(u,\alpha-\xi)$.

With this notation, we are able to formulate

######  Algorithm 2.

(SESOP)

For $i\in I_{n}\subset\mathbb{N}_{0}$ a finite index set, the adjoint operator
$\mathcal{A}^{*}$ of $\mathcal{A}$ and search directions
$\mathcal{A}^{*}w_{n,i}$ for $w_{n,i}\in Y$, the $(n+1)$-th iterate $f_{n+1}$
can be calculated by

(2) $\displaystyle f_{n+1}=P_{H_{n}}(f_{n})$

as the metric projection $P_{H_{n}}$ of the $n$-th iterate $f_{n}$ onto the
intersection of hyperplanes

$\displaystyle H_{n}=\bigcap_{i\in I_{n}}H(\mathcal{A}^{*}w_{n,i},\langle
w_{n,i},g\rangle_{Y}).$

######  Remark 3.

For arbitrary elements $w\in Y$ the hyperplane $H(\mathcal{A}^{*}w,\langle
w,g\rangle_{Y})=\\{f\in X:\langle\mathcal{A}^{*}w,f\rangle_{X}=\langle
w,g\rangle_{Y}\\}$ includes $\mathcal{M}_{\mathcal{A},g}$, since for every
$z\in\mathcal{M}_{\mathcal{A},g}$ holds

(3) $\displaystyle\langle\mathcal{A}^{*}w,z\rangle_{X}=\langle
w,\mathcal{A}z\rangle_{Y}=\langle w,g\rangle_{Y}.$

This gives an intuition as to why the hyperplanes are chosen as above.

The calculation of the next iterate $f_{n+1}$ can also be considered an update
to the current iterate $f_{n}$ with multiple search directions per iteration
and regulation of the step size:

######  Lemma 4.

(Calculation of projection onto intersection of hyperplanes)

Let $H=\bigcap_{i=1}^{n}H(u_{i},\alpha_{i})$ be the intersection of
hyperplanes for $u_{i}\in X$ and $\alpha_{i}\in\mathbb{R}$ for
$i=\\{1,\dots,n\\},n\in\mathbb{N}$. Then the projection of $f\in X$ onto $H$
can be determined as

(4) $\displaystyle P_{H}(f)=f-\sum_{i=1}^{n}\tilde{t_{i}}u_{i},$

where $\tilde{t}=(\tilde{t_{1}},\dots,\tilde{t_{n}})$ minimizes the function

$\displaystyle
h(t)=\frac{1}{2}\big{\lVert}f-\sum_{i=1}^{n}t_{i}u_{i}\big{\rVert}_{X}^{2}+\sum_{i=1}^{n}t_{i}\alpha_{i}.$

If the elements $u_{1},\dots,u_{n}$ are linearly independent, $h$ has a unique
minimal solution $\tilde{t}$.

A proof of Lemma 4 can be found in Proposition 3.12 in [26]. If the search
directions $\mathcal{A}^{*}w_{n,i}$ for $i\in I_{n}$ and the finite index set
$I_{n}$ are chosen correctly, the iterates $f_{n}$ converge to a solution $f$
with $\mathcal{A}f=g$, see e.g. Section 6.2.2 of [28].

##### RESESOP

RESESOP is the regularizing version of SESOP and can take into account noisy
data and an inexact forward operator.

Let the global modelling error $\eta>0$ of the inexact forward operator
$\mathcal{A}^{\eta}$ and the global noise level $\delta$ of the noisy data
$g^{\delta}$ be given by

$\displaystyle\lVert\mathcal{A}^{\eta}f-\mathcal{A}f\rVert\leq\eta\rho\text{
for all }f\in B_{\rho}(0)\subset X,\quad\lVert
g-g^{\delta}\rVert_{Y}\leq\delta$

for $\rho>0$.

To regularize the SESOP method, the projection onto hyperplanes is then
changed to a projection onto stripes

$\displaystyle H_{n,i}^{\eta,\delta}:=\big{\\{}f\in X\ :\
\left\lvert\left\langle(\mathcal{A}^{\eta})^{*}w_{n,i}^{\eta,\delta},f\right\rangle-\left\langle
w_{n,i}^{\eta,\delta},g^{\delta}\right\rangle\right\rvert\leq(\delta+\eta\rho)\lVert
w_{n,i}^{\eta,\delta}\rVert\big{\\}}$

As in the case of exact data, the stripes contain the solution set
$\mathcal{M}_{\mathcal{A}f=g}$.

The discrepancy principle [22] assumes that in the inexact case the solution
cannot get better than the noise level and modeling error and yields a finite
stopping index $n^{*}$ as soon as

$\displaystyle\lVert\mathcal{A}^{\eta}f_{n^{*}}-g^{\delta}\rVert_{Y}\leq\tau(\eta\rho+\delta)$

for a fixed $\tau>1$. A proof of the convergence of this version of RESESOP
can be found in [3]. For the case of error-free forward operators, see, e.g.,
[28].

##### RESESOP-Kaczmarz

We summarize some elementary results from [3]. Let
$\mathcal{A}_{j,m}\in\mathcal{L}(X,Y_{j,m})$ and
$\mathcal{A}_{j,m}^{\eta}\in\mathcal{L}(X,Y_{j,m})$ be the inexact version of
$\mathcal{A}_{j,m}$.

In a semi-discrete setting, we can use the RESESOP-Kaczmarz method to solve a
system of equations

$\displaystyle\mathcal{A}_{j,m}f$ $\displaystyle=g_{j,m},$
$\displaystyle\lVert\mathcal{A}_{j,m}^{\eta}-\mathcal{A}_{j,m}\rVert_{\mathcal{L}(X,Y_{j,m})}$
$\displaystyle\leq\eta_{j,m},$ $\displaystyle\lVert
g^{\delta}_{j,m}-g_{j,m},\rVert_{Y_{j,m}}$ $\displaystyle\leq\delta_{j,m}$

for $j\in\mathcal{J}$, $m\in\mathcal{M}$ and finite index sets $\mathcal{J}$
and $\mathcal{M}$.

Combining RESESOP with Kaczmarz’s method allows the inclusion of local
modeling errors and local noise in the reconstruction process. This is
especially an advantage for dynamic computer tomography, since the model error
of the dynamic forward operator changes with time.

######  Remark 5.

In case of computerized tomography with $K\in\mathbb{N}$ source positions
$\varphi_{k}\in[0,2\pi]$, $k\in\mathcal{K}\coloneqq\\{0,\dots K-1\\}$, and
$L\in\mathbb{N}$ detector points $s_{l}\in[-1,1]$,
$l\in\mathcal{L}\coloneqq\\{0,\dots L-1\\}$, the algorithm allows four
different modalities:

  1. V1:

$\mathcal{J}=\mathcal{K}\times\mathcal{L}$ and $\mathcal{M}=\\{0\\}$ results
in Kaczmarz projections with respect to both source positions and detector
points.

  2. V2:

$\mathcal{J}=\mathcal{K}$ and $\mathcal{M}=\mathcal{L}$ results in Kaczmarz
projections with respect to the source positions and averaging with respect to
the detector points.

  3. V3:

$\mathcal{J}=\mathcal{L}$ and $\mathcal{M}=\mathcal{K}$ results in Kaczmarz
projections with respect to the detector points and averaging with respect to
the source directions.

  4. V4:

$\mathcal{J}=\\{0\\}$ and $\mathcal{M}=\mathcal{K}\times\mathcal{L}$ results
in averaging with respect to both detector points and search directions.

In this contribution, we use the version V1, which is the most precise
algorothm in the sense that each subproblem is considered individually. We
also emphasize that the modalities V2, V3 and V4 allow for (partial)
parallelization of the code in contrast to V1.

######  Algorithm 6.

(RESESOP-Kaczmarz)

Let $[n]=n\text{ mod }J$. Fix an initial value
$f_{0}^{\eta,\delta}\in\mathcal{B}_{\rho}(0)$ and precision constants
$\tau_{j,m}>1$ for $j\in\mathcal{J}$ and $m\in\mathcal{M}$.

While the discrepancy principle is not yet fulfilled for all $m\in\mathcal{M}$
and the current iterate is not close enough to the solution of the current
subproblem, i.e., $D_{n}^{\eta,\delta}\neq\emptyset$ with

$\displaystyle
D_{n}^{\eta,\delta}\coloneqq\\{m\in\mathcal{M}:\lVert\mathcal{A}^{\eta}_{[n],m}f_{n}^{\eta,\delta}-g^{\delta}_{[n],m}\rVert_{Y_{[n],m}}>\tau_{[n],m}(\eta_{[n],m}\rho+\delta_{[n],m})\\}$

for each $n\in\mathbb{N}_{0}$, perform the next iteration step

$\displaystyle f_{n+1}^{\eta,\delta}\coloneqq
P_{H_{n}^{\eta,\delta,\rho}}\left(f_{n}^{\eta,\delta}\right),\quad
H_{n}^{\eta,\delta,\rho}\coloneqq\bigcap_{i\in
I_{n}}H\left(u_{n,i}^{\eta,\delta},\alpha_{n,i}^{\eta,\delta},\xi_{n,i}^{\eta,\delta}\right).$

We choose for $n\in\mathbb{N}$ and $i\in
I_{n}^{\eta,\delta}\subset\\{n-N+1,\dots,n\\}\cap\mathbb{N}_{0}$,
$N\in\mathbb{N}\setminus\\{0\\}$

$\displaystyle u_{n,i}^{\eta,\delta}$
$\displaystyle\coloneqq\sum_{m\in\mathcal{M}}\left(A_{[i],m}^{\eta}\right)^{*}w_{n,i}^{\eta,\delta,m}$
$\displaystyle\alpha_{n,i}^{\eta,\delta}$
$\displaystyle\coloneqq\sum_{m\in\mathcal{M}}\langle
w_{n,i}^{\eta,\delta,m},g^{\delta}_{[i],m}\rangle_{Y}$
$\displaystyle\xi_{n,i}^{\eta,\delta}$
$\displaystyle\coloneqq\sum_{m\in\mathcal{M}}(\delta+\eta\rho)\lVert
w_{n,i}^{\eta,\delta,m}\rVert_{Y_{[i],m}}$ $\displaystyle
w_{n,i}^{\eta,\delta,m}$
$\displaystyle\coloneqq\begin{cases}\mathcal{A}^{\eta}_{[i],m}f_{n}^{\eta,\delta}-g^{\delta}_{[i],m}&\text{if
}m\in D_{n}^{\eta,\delta}\\\ 0&\text{if }m\notin
D_{n}^{\eta,\delta}\end{cases}.$

One full iteration over all $J$ subproblems is completed if $[n]=0$. Stop the
algorithm after a full iteration at the index $[n^{*}]=0$ when
$f_{n^{*}}^{\eta,\delta}=f_{n^{*}-J}^{\eta,\delta}$. This means that the set
$D^{\eta,\delta}_{n}$ has been empty at least for
$n\in\\{n^{*},\dots,n^{*}-J\\}$ and therefore the reconstruction cannot get
better.

######  Theorem 7.

(Convergence of SESOP-Kaczmarz)

Let $\\{f_{n}\\}_{n\in\mathbb{N}}$ be the sequence of iterates generated by
the RESESOP-Kaczmarz algorithm with an initial value $f_{0}\in X$ and exact
data, i.e., $\eta=\delta=0$.

If the optimization parameters $\lvert t_{n,i}\rvert\leq t$ for $i\in I_{n}$
are bounded by a constant $t>0$, then $\\{f_{n}\\}_{n\in\mathbb{N}}$ converges
strongly to a solution $f$ of the semi-discretized problem
$\mathcal{A}_{j,m}f=g_{j,m}$ for $j\in\mathcal{J}$, $m\in\mathcal{M}$.

######  Theorem 8.

(Convergence and Regularization of RESESOP-Kaczmarz)

Let all assumptions from the previous theorem be true and the set of search
directions be linearly independent. Then

  1. 1.

the RESESOP-Kaczmarz algorithm with the adapted discrepancy principle as a
stopping criterion leads to a finite stopping index $n_{*}$,

  2. 2.

for each iteration $n\in\mathbb{N}$ holds
$\lim_{(\eta,\delta)\to(0,0)}f_{n}^{\eta,\delta}=f_{n}$,

  3. 3.

$\lim_{\eta,\delta\to(0,0)}f^{\eta,\delta}_{r_{*}}=P_{M_{\mathcal{A}_{j,m}f=g_{j,m}}\cap\mathcal{B}_{\rho}(0)}(f_{0})$.

Proofs for both theorems can be found in [3].

We give a pseudo-code for implementing RESESOP-Kaczmarz, discretized with
respect to both source positions and detector points:

Listing 1: Pseudo code for V1 of RESESOP-Kaczmarz with two search directions

⬇

1$\mathbf{Input}$: $\mathcal{K}\text{ set of source positions,
}\mathcal{L}\text{ set of detector points, }\tau_{k,l}>1\text{ for }k\in
K,l\in L,$

2$g^{\delta}\text{ perturbed sinogram, }\eta,\delta\text{ model/data error for
each }k\in\mathcal{K},l\in\mathcal{L},\text{ max\\_iterations},$

3

4$\text{stopped}_{k,l}$ = $1\text{ for all }k\in\mathcal{K},l\in\mathcal{L}$

5$f,u_{\text{old}}$ = $0$

6$\alpha_{\text{old}},\xi_{\text{old}},\text{full\\_iter}$ = $0$

7while $\exists k,l:\text{stopped}_{k,l}\neq 0\text{ and
}\text{full\\_iter}<\text{max\\_iterations}$ do

8 for $k\in K$ do

9 for $l\in L$ do

10 res = $\mathcal{A}^{\eta}_{k,l}f-g^{\delta}_{k,l}$

11 if $\lvert\text{res}\rvert$ $\leq$
$\tau_{k,l}\cdot(\eta_{k,l}+\delta_{k,l})$ then

12 $\text{stopped}_{k,l}$ = $0$

13 else

14 $u_{\text{new}}$ = $(\mathcal{A}^{\eta}_{k,l})^{*}\cdot\text{res}$

15 if $\lVert u_{\text{new}}\rVert^{2}$ > $0$ then

16 $\alpha_{\text{new}}$ = $\text{res}\cdot g^{\delta}_{k,l}$

17 $\xi_{\text{new}}$ = $\lvert\text{res}\rvert\cdot(\eta_{k,l}+\delta_{k,l})$

18 # determine
$\tilde{f}=P_{H(u_{\text{new}},\alpha_{\text{new}}+\xi_{\text{new}})}(f)$

19 $\tilde{f}$ =
$f-\frac{\lvert\text{res}\rvert(\lvert\text{res}\rvert-(\eta_{k,l}+\delta_{k,l}))}{\lVert
u_{\text{new}}\rVert^{2}}u_{\text{new}}$

20 # determine
$f=P_{H(u_{\text{new}},\alpha_{\text{new}}+\xi_{\text{new}})\cap
H_{\gtrless}(u_{\text{old}},\alpha_{\text{old}}\pm\xi_{\text{old}})}(\tilde{f})$

21 if $\tilde{f}\in H(u_{\text{old}},\alpha_{\text{old}}+\xi_{\text{old}})$
then

22 $f$ = $\tilde{f}$

23 else if $\lVert u_{\text{new}}\rVert^{2}\cdot\lVert
u_{\text{old}}\rVert^{2}-\langle u_{\text{new}},u_{\text{old}}\rangle^{2}\neq
0$

24 Decide whether

25 1. $\tilde{f}\in
H_{>}(u_{\text{old}},\alpha_{\text{old}}+\xi_{\text{old}})\text{ or}$

26 2. $\tilde{f}\in
H_{<}(u_{\text{old}},\alpha_{\text{old}}-\xi_{\text{old}})$

27 $t$ = $\frac{\langle
u_{\text{old}},\tilde{f}\rangle-(\alpha_{\text{old}}\pm\xi_{\text{old}})}{\lVert
u_{\text{new}}\rVert^{2}\cdot\lVert u_{\text{old}}\rVert^{2}-\langle
u_{\text{new}},u_{\text{old}}\rangle^{2}}$

28 $f$ = $\tilde{f}+\langle u_{\text{new}},u_{\text{old}}\rangle\cdot t\cdot
u_{\text{new}}-\lVert u_{\text{new}}\rVert^{2}\cdot t\cdot u_{\text{old}}$

29 else:

30 $f$ = $\tilde{f}$

31 $\text{set negative values in }f\text{ to zero}$

32 $\alpha_{\text{old}}$ = $\alpha_{\text{new}}$

33 $\xi_{\text{old}}$ = $\xi_{\text{new}}$

34 $u_{\text{old}}$ = $u_{\text{new}}$

35 end

36 end

37 full_iter = $\text{full\\_iter}+1$

38end

39$\mathbf{Output}$: $f$

###  2.2 Motion Estimation via Landmark Detection

For the purpose of a proof of concept, we consider homogeneous rectangular
objects and restrict ourselves to the rather simple case of an affine linear
motion $\Gamma_{\theta}$ with constant speed of the form

(5)
$\displaystyle\begin{split}\Gamma_{\theta}x&=A_{\theta}x+b_{\theta},\quad\theta\in
S_{1},\\\ A_{\theta}&=A(t_{\theta})=I+\frac{t_{\theta}}{N-1}(A-I),\\\
b_{\theta}&=b(t_{\theta})=\frac{t_{\theta}}{N-1}b\end{split}$

for the unit circle $S_{1}$, the identity matrix $I\in\mathbb{R}^{2\times 2}$,
$A\in\mathbb{R}^{2\times 2}$, $b\in\mathbb{R}^{2}$, $N\in\mathbb{N}$ the
number of time points (i.e., the number of scanning angles $\varphi$), and
$t_{\theta}\in\\{0,\ldots,N-1\\}$.

In this case, it is sufficient to know the coordinates of four landmarks of
the object, each at the start and end time. By plugging these into (5) we get
a linear system of eight equations containing six unknowns for determining the
motion. The motion estimation then consists in computing the least-squares
solution of this overdetermined system.

Figure 2: CNN6 architecture

####  2.2.1 Data-driven landmark detection

As a proof-of-concept, we trained a small network for the purpose of landmark
detection for randomly generated rectangles, i.e., it is trained to determine
the coordinates of the four corners of a rectangle.

##### Architecture

The implemented architecture, denoted by CNN6 (see Figure 2), was inspired by
[8] where it was originally used for facial landmark detection.

##### Training Data Set

The network was trained in two steps with two different data sets, each
containing 500 pairs of images of size $128\times 128$ as an input for the
network and the corresponding coordinates of the landmarks as the output: For
pre-training, a data set with uncorrupted rectangular phantoms was used.

After pre-training, the net was trained with a data set containing the
reconstructions of the inital state of moving rectangles by the RESESOP-
Kaczmarz method with three iterations.

##### Hyperparameters

As a regularization technique, we used early stopping with a patience of 200
epochs and allowed a maximum of 10 000 epochs.
We remark that the network in the original article [8] was trained for 120 000
epochs, leading to our choice of this relatively high patience.

We used Adam [17] as an optimizer with an initial step size of $10^{-6}$, a
batch size of 8, and the mean squared error as a loss function.

##### Training Results

The pre-training took 10 000 epochs and already reached a good precision as
can be seen in figure 3.

Figure 3: Pre-training (left) and training loss (right)

The second phase of the training with the reconstructed data set only took 500
epochs until the early stopping ended the training. Since the reconstructions
are not as sharp and precise as the constructed phantoms, it is not surprising
that the validation loss does not get as small as on the pre-training data
set.

###  2.3 Dynamic Filtered Backprojection

To derive a dynamic filtered backprojection algorithm for affine deformations,
the method of the approximate inverse can be used. A detailed introduction to
the method of the approximate inverse can be found in [27]. This method can be
used to solve operator equations

(6) $\displaystyle\mathcal{A}f=g$

for a linear invertible operator $\mathcal{A}:X\to Y$ with
$X=L^{2}(\Omega_{1},\mu_{1})$ and $Y=L^{2}(\Omega_{2},\mu_{2})$ for
$\Omega_{i}\subset\mathbb{R}^{n_{i}}$ open, bounded domains, $\mu_{i}$
measures on $\Omega_{i}$, $i\in\\{1,2\\}$ and $L^{2}$ the space of Lebesgue
square-integrable functions.

To reconstruct $f$, a smoothed version $f^{\gamma}$ with

(7) $\displaystyle f(x)\approx f^{\gamma}(x)=\langle
f,\delta^{\gamma}(\cdot,x)\rangle_{X}$

is calculated for a suitable mollifier $\delta^{\gamma}(y,x)\in
L^{2}(\Omega_{1}\times\Omega_{1},\mu_{1}\times\mu_{1})$.

######  Definition 9.

(Mollifier)

A function $\delta^{\gamma}(\cdot,x)\in L^{2}(\mathbb{R}^{n_{1}},\mu_{1})$ for
$x\in\mathbb{R}^{n_{1}}$ is called a mollifier if

  * •

for all $\gamma>0$ it holds

$\displaystyle\int_{\mathbb{R}^{n_{1}}}\delta^{\gamma}(y,x)d\mu_{1}(y)=1.$

  * •

for $\gamma\to 0$ $f_{\gamma}$ converges to $f$ in
$L^{2}(\Omega_{1},\mu_{1})$.

With equation (7) it is not yet possible to calculate $f_{\gamma}$ since the
exact solution $f$ appears in the expression. This problem can be overcome:

For a mollifier $\delta^{\gamma}(\cdot,x)$ and a solution $\Psi^{\gamma}(x)\in
Y$ of the auxiliary problem

(8) $\displaystyle\mathcal{A}^{*}\Psi^{\gamma}(x)=\delta^{\gamma}(\cdot,x)$

we find that

$\displaystyle f^{\gamma}(x)=\langle
f,\delta^{\gamma}(\cdot,x)\rangle_{X}=\langle
f,\mathcal{A}^{*}\Psi^{\gamma}(x)\rangle_{X}=\langle\mathcal{A}f,\Psi^{\gamma}(x)\rangle_{Y}=\langle
g,\Psi^{\gamma}(x)\rangle_{Y}.$

######  Remark 10.

Since the operator $\mathcal{A}$ is invertible, $\delta^{\gamma}(\cdot,x)$ is
in the range of the adjoint operator $\mathcal{A}^{*}$ and therefore the
auxiliary problem (8) has a solution $\Psi^{\gamma}(x)\in Y$.

Now the approximate inverse depending on the chosen mollifier can be defined
as follows:

######  Definition 11.

(Approximate Inverse)

Let $\mathcal{A}:X\to Y$ be a linear, bounded and invertible operator and
$\delta^{\gamma}$ be a mollifier. Then the linear mapping
$\tilde{\mathcal{A}}_{\gamma}:Y\to X$ with

$\displaystyle\tilde{\mathcal{A}}_{\gamma}g(x)=\langle
g,\Psi^{\gamma}(x)\rangle_{Y}$

is called the approximate inverse of $\mathcal{A}$. The solution
$\Psi^{\gamma}(x)$ of (8) is called the reconstruction kernel associated to
$\delta^{\gamma}$.

Under certain conditions on $\Psi^{\gamma}(x)$ the approximate inverse
$\tilde{\mathcal{A}}_{\gamma}$ is continuous and
$\tilde{\mathcal{A}}_{\gamma}g$ converges pointwise to
$\mathcal{A}^{-1}_{\gamma}g$ for $\gamma\to 0$. Further details can be found
in the second chapter of [27].

An advantage of the computation of reconstruction kernels is that the
auxiliary problem (8) does not depend on the noisy measured data $g^{\delta}$.
Therefore the reconstruction kernel is independent of noise.

On the other hand, the computation of the reconstruction kernel is time-
consuming since it needs to be calculated for each reconstruction point $x$ in
the domain of $f$. However, it is possible to simplify the calculation under
certain invariance conditions of operator $\mathcal{A}$. This is shown in the
next lemma (see Lemma 3.2 in [10] and Theorem 2.6 in [27]).

######  Lemma 12.

Let $T_{1}^{x}:X\to X$ and $T_{2}^{x}:Y\to Y$ for $x\in\Omega_{1}$ be linear
operators with

(9) $\displaystyle T_{1}^{x}\mathcal{A}^{*}=\mathcal{A}^{*}T_{2}^{x}$

and let the mollifier $\delta^{\gamma}$ be generated by $T_{1}^{x}$, meaning
that there exists $x^{*}\in\Omega_{1}$ with

$\delta^{\gamma}(\cdot,x)=T_{1}^{x}\delta^{\gamma}(\cdot,x^{*}).$

If $\Psi^{\gamma}(x^{*})$ solves
$\mathcal{A}^{*}\Psi^{\gamma}(x^{*})=\delta^{\gamma}(\cdot,x^{*})$, then
$\Psi^{\gamma}(x)=T_{2}^{x}\Psi^{\gamma}(x^{*})$ solves

$\displaystyle\mathcal{A}^{*}\Psi^{\gamma}(x)=\delta^{\gamma}(\cdot,x)$

for any $x\in\Omega_{1}$.

###### Proof.

It holds

$\displaystyle\delta^{\gamma}(\cdot,x)=T_{1}^{x}\delta^{\gamma}(\cdot,x^{*})=T_{1}^{x}\mathcal{A}^{*}\Psi^{\gamma}(x^{*})=\mathcal{A}^{*}T_{2}^{x}\Psi^{\gamma}(x^{*}).$

Therefore, $T_{2}^{x}\Psi^{\gamma}(x^{*})$ solves the auxiliary problem
$\mathcal{A}^{*}\Psi^{\gamma}(x)=\delta^{\gamma}(\cdot,x)$. ∎

If Lemma (12) can be applied, only one auxiliary problem has to be solved. In
this case, the method of the approximate inverse can be formulated as

(10) $\displaystyle\tilde{\mathcal{A}}_{\gamma}g(x)=\langle
g,T_{2}^{x}\Psi_{\gamma}(x^{*})\rangle_{Y}$

Now we can apply the method of the approximate inverse to dynamic computerized
tomography as it is done in [10]:

######  Lemma 13 (Application to Affine Deformations in Dynamic CT).

Let the motion be given by

$\displaystyle\Gamma_{\theta}x=C_{\theta}x+b_{\theta},\quad\theta\in S_{1},$

for $b_{\theta}\in\mathbb{R}^{2}$ and $C_{\theta}\in\mathbb{R}^{2\times 2}$
with continuously differentiable entries, $\det C_{\theta}\neq 0$ for all
measured angles $\theta=(\cos(\varphi),\sin(\varphi))^{T}$ and

$\displaystyle
h(\theta):=(C_{\theta}^{-*}\theta)_{1}\cdot\frac{\partial}{\partial\varphi}(C_{\theta}^{-*}\theta)_{2}-(C_{\theta}^{-*}\theta)_{2}\cdot\frac{\partial}{\partial\varphi}(C_{\theta}^{-*}\theta)_{1}\neq
0,$

where $C_{\theta}^{-*}:=\big{(}C_{\theta}^{-1}\big{)}^{*}$.

Then the auxiliary problem can be solved analytically by using the inverse
dynamic Radon operator

$\displaystyle(R^{\Gamma})^{*}\Psi_{x}^{\gamma}=\delta_{x}^{\gamma}=(R^{\Gamma})^{-1}(R^{\Gamma})\delta_{x}^{\gamma}=(R^{\Gamma})^{*}\left(\frac{1}{4\pi}\lvert\text{det}C_{\theta}\rvert^{2}\lvert
h(\theta)\rvert\ell^{-1}(R^{\Gamma})\delta_{x}^{\gamma}\right)$

yielding the reconstruction kernel

(11)
$\displaystyle\Psi_{x}^{\gamma}(\theta,s)=\frac{1}{4\pi}\lvert\text{det}C_{\theta}\rvert^{2}\lvert
h(\theta)\rvert\ell^{-1}(R^{\Gamma})\delta_{x}^{\gamma}$

with the Riesz potential $\ell$ defined by
$\mathcal{F}\ell^{-1}g(\sigma)=\lvert\sigma\rvert\mathcal{F}g(\sigma)$ and the
two-dimensional Fourier transform $\mathcal{F}$ given by
$\mathcal{F}f(\xi)=\frac{1}{2\pi}\int_{\mathbb{R}^{2}}f(x)e^{-ix^{T}\xi}\,\mathrm{d}x$.

A derivation of this result can be found in Section 4 of [10].

Since the dynamic Radon operator fulfills the invariance condition

$T^{x}(R^{\Gamma})^{*}=(R^{\Gamma})^{*}T^{(0,(C_{\theta}^{-1}x)^{T}\theta)}$

from Lemma 12 for the shift operators

$\displaystyle T^{x}f(y)$ $\displaystyle=f(y-x)$ $\displaystyle
T^{\left(0,(C_{\theta}^{-1}x)^{T}\theta\right)}g(\theta,s)$
$\displaystyle=g\left(\theta,s-(C_{\theta}^{-1}x)^{T}\theta\right),$

it is enough to compute the reconstruction kernel
$\Psi^{\gamma}_{0}=\Psi^{\gamma}$ only for the point $x^{*}=0$ for a mollifier
with
$\delta^{\gamma}(\cdot,x)=T^{x}\delta^{\gamma}(\cdot,0)=\delta^{\gamma}(\cdot-x,0)$,
e.g., the Gaussian $\delta_{G}^{\gamma}$ with

$\displaystyle\delta_{G}^{\gamma}(y,x)=\frac{1}{2\pi\gamma^{2}}\exp\left(-\frac{\lVert
y-x\rVert^{2}}{(2\gamma^{2})}\right).$

This leads to the dynamic filtered backprojection formula

(12) $\displaystyle f^{\gamma}(x)$ $\displaystyle=\langle
g,T^{\left(0,(C_{\theta}^{-1}x)^{T}\theta\right)}\Psi^{\gamma}\rangle_{L^{2}(S^{1}\times\mathbb{R})}$
(13)
$\displaystyle=\int_{S^{1}}\int_{\mathbb{R}}g(\theta,s)\Psi^{\gamma}\left(\theta,s-(C_{\theta}^{-1}x)^{T}\theta\right)\,\mathrm{d}s\,\mathrm{d}\theta.$

######  Theorem 14.

(Reconstruction Kernel for the Gaussian Mollifier) The reconstruction kernel
for the Gaussian mollifier is given by

$\displaystyle\psi^{\gamma}(\theta,s)=\frac{\lvert\det C_{\theta}\rvert\lvert
h(\theta)\rvert}{4\pi^{2}\gamma^{2}\lVert
C_{\theta}^{-*}\theta\rVert^{2}}\left(1-\frac{\sqrt{2}\left[s+(C_{\theta}^{-1}b_{\theta})^{T}\theta\right]}{\gamma\lVert
C_{\theta}^{-*}\theta\rVert}D\left(\frac{s+(C_{\theta}^{-1}b_{\theta})^{T}\theta}{\sqrt{2}\gamma\lVert
C_{\theta}^{-*}\theta\rVert}\right)\right)$

with the Dawson integral

$\displaystyle D(s)=e^{-s^{2}}\int_{0}^{s}e^{t^{2}}\,\mathrm{d}t.$

###### Proof.

The computation can be done as in the static case [18] by using (11) and
equations 7.4.6 and 7.4.7 from [1]. The result can also be found in [9]. ∎

By applying the same kind of discretization as in the static case, we
altogether obtain the following pseudo code:

Listing 2: pseudo code for the dynamic filtered backprojection

⬇

1$\mathbf{Input}$: $p\text{: number of angles, }\tilde{q}\text{: uneven number
of offset-values, }\Psi^{\gamma}\text{: dynamic filter,}$

2$\text{numpix: number of pixels, }g(\theta_{j},s_{k})\text{: sinogram for
}\theta_{i}=(\cos(\varphi_{i}),\sin(\varphi_{i}))^{T},$

3$\varphi_{i}=\frac{i\pi}{p},s_{k}=hk\text{ for
}h=1/q,k=0,\ldots,\tilde{q}-1,q=(\tilde{q}-1)/2,i=0,...,p-1$

4

51.for $i$ = $0,\ldots,p-1$ do

6 for $k$ = $0,\ldots,\tilde{q}-1$ do

7 $v_{i,k}$ =
$h\sum_{j=0}^{\tilde{q}-1}\Psi^{\gamma}(s_{k}-s_{j})g(\theta_{i},s_{j})$

8 end

9 end

10

112.for $i$ = $0,\ldots,\text{numpix}-1$ do

12 $x_{1}$ = $1-\frac{2}{\tilde{q}}(i+1/2)$

13 for $j$ = $0,\ldots,\text{numpix}-1$ do

14 $x_{2}$ = $-1+\frac{2}{\tilde{q}}(j+1/2)$

15 sum = 0

16 if $x_{1}^{2}+x_{2}^{2}>1$ then

17 pass

18 else

19 for $l$ = $0,\ldots,p-1$ do

20 $s$ = $(C_{\theta}^{-1}x)^{T}\theta_{l}$

21 $k$ = $\lfloor sq\rfloor+q$

22 $\mu$ = $sq-k+q$

23 sum += $(1-\mu)v_{l,k}+\mu v_{l,k+1}$

24 end

25 $f^{\mathrm{FBI}}_{i,j}$ = $\frac{2\pi}{p}\cdot\text{sum}$

26 end

27 end

28 $\mathbf{Output}$: $f^{\mathrm{FBI}}$

##  3 Numerical Results

In the following, we test our method first for two simulated phantoms, the
first one performing a shift and the second one a stretch. In both cases,
landmark detection is performed via deep learning with the network described
in Section 2.2.1. Then we test our method on real data, where we perform
landmark detection by hand.

###  3.1 Simulated Data

##### Synthetic Measurements

Both simulated phantoms are of size $512\times 512$, contain one rectangle
with greyscale values in $[0,1]$, and perform an affine linear motion with
constant speed as in (5). In the first case, the rectangle moves with a simple
shift (rigid body motion) defined in the pixel domain by

$\displaystyle C_{1}(t_{\theta})=I,\quad
b_{1}(t_{\theta})_{pix}=\frac{t_{\theta}}{N-1}\begin{pmatrix}51\\\
51\end{pmatrix}$

In the second case, the rectangle is stretched as described by the motion

$\displaystyle
C_{2,\theta}=C_{2}(t_{\theta})=I+\frac{t_{\theta}}{N-1}(C-I),\quad
b_{2}(t_{\theta})=0$

for $t_{\theta}\in\\{0,\ldots,N-1\\}$, $N=450$ and $C=\begin{pmatrix}2&0\\\
0&1\end{pmatrix}.$

Next, we compute the dynamic sinogram as well as the two static sinograms for
both the initial and the final position to calculate the inexactness
$\eta(s,\varphi)$ for 300 offset values and 450 angles, and add uniformly
distributed noise in $[-0.02,0.02]$. These yield the input to our hybrid
three-step reconstruction scheme.

Figure 4: Landmark detection on reconstructions from RESESOP-Kaczmarz without
noise

##### Reconstructions

We reconstruct the object twice, once at the start time and again at the final
time point, using three iterations of the RESESOP-Kaczmarz algorithm with the
calculated inexactness. Note that we use an image size of only $128\times 128$
pixels for these reconstructions to reduce the computation time. For all our
reconstructions with RESESOP-Kaczmarz we choose $\tau=1.00001$. As can be seen
in Figure 4, the reconstructions are not satisfactory, but the quality is
sufficient to see the shift of the object. These reconstructions can now be
passed on to the trained CNN6 network described in Section 2.2.1 to
automatically detect the landmarks (here: corners), the results are shown in
Figure 4.

Figure 5: Shifted rectangle without noise Figure 6: Shifted rectangle with
noise

In the first case the shift is computed as the mean of the difference between
the corner coordinates from the end point and the corner coordinates from the
start point. In the second case, the entries of matrix $C$ are computed as a
least-squares solution of the respective system of linear equations.

To avoid committing an inverse crime we choose the size of all reconstruction
as $487\times 487$, which differs from the original phantom size of $512\times
512$. Some basic examples of inverse crimes and how to avoid them can be found
in [22, Chapter II].

For comparison, we computed RESESOP-Kaczmarz reconstructions with $30$
iterations. While the reconstruction of the RESESOP-Kaczmarz algorithm still
has a small motion artifact under the lower right corner as well as inside the
rectangle, the reconstruction of the combined method turns out to be a
slightly blurred version of the ground truth (see Figures 5 and 6). Our hybrid
method is thus able to recover the object without using the actual motion as
an input.

In the case of the stretched rectangle (figure 7), the output of the combined
method has stronger artifacts since the landmarks were not detected precisely,
and therefore also the motion estimation was faulty.

Considering the run time, our combined method leads to a speed-up and only
needs a little more than half of the computational time that is needed for 30
iterations of RESESOP-Kaczmarz.

Figure 7: Stretched Rectangle with Noise

###  3.2 Real Data

The data was measured by a cone-beam CT scanner designed and constructed in-
house in the Industrial Mathematics Computed Tomography Laboratory at the
University of Helsinki and was preprocessed using the Helsinki Tomography
Toolbox for Matlab [21]. The measurement settings were almost the same as in
the open access CT dataset [15] of a scanned seashell, where the following
measurement description is based on.

The scanner consists of a motorized rotation stage (Thorlabs CR1-Z7), a
molybdenum target X-ray tube (Oxford Instruments XTF5011), and a 12-bit,
2240x2368 pixel, energy-integrating flat panel detector (Hamatsu Photonics
C7942CA-22).

360 X-ray projections were acquired using an angle increment of 1 degree. The
X-ray source voltage and tube current were set at 45 kV and 1 mA,
respectively. The exposure time of the flat panel detector was set to 1200 ms.

Two correction images were acquired before scanning the sample. A dark current
image was created by averaging 255 images taken with the X-ray source off. A
flat-field image was created by averaging 255 images taken with the X-ray
source switched on with no sample placed in the scanner. After the scan, dark
current and flat-field corrections [29] were applied to each projection image
using the Hamamatsu HiPic imaging software version 9.3.

Fan-beam sinogram data was created by using only the measured intensity values
from the central plane of the X-ray cone beam. Due to a slightly misaligned
center of rotation in the scanner, the CT reconstructions can appear blurry.
It was empirically observed that this problem can be compensated for quite
well by shifting each projection to the left by five pixels, using circular
boundary conditions, before performing any other operations on the
projections.

Figure 8: RESESOP-Kaczmarz Reconstructions for the first and last time step
after three iterations

We scanned a composite Lego brick, shown in figure 8, with these settings for
four different time steps. For each measurement, we moved the brick one knob
forward such that it reached the other end at the last time step. Then we
computed the parallel-beam sinograms ranging from 0 to 180 degrees with the
inbuild Matlab binning algorithm ’fan2para’ and built the dynamic sinogram by
taking the first 45 degrees from the first time step, 46 to 90 degrees from
the second time step, and so on.

Figure 9: Comparison of the 3 Reconstruction Methods in relation to the FBP
computed with the static sinogram

As inexactness for the RESESOP-Kaczmarz algorithm, we took the absolute
differences between the dynamic sinogram and the sinogram of the first resp.
last position of the lego brick. We estimated the noise by taking the mean of
the left blank part of the dynamic sinogram.

Three iterations of RESESOP-Kaczmarz to determine the initial and final state,
and landmark detection by hand lead to the reconstructions of size $128\times
128$, shown in figure 8.

Afterwards, the mean of each coordinate shift yield the total movement of the
Lego brick. To compute the reconstruction of the combined algorithm with
dynamic filtered backprojection, we assumed that the movement took place in
four steps with the same amount of projections per step. The results of size
$512\times 512$ are displayed in Figure 9. As in the experiments with the
simulated data, the reconstruction by the combined method shows fewer
artifacts and is about twice as fast as a pure RESESOP-Kaczmarz
reconstruction.

##  4 Conclusion

In this article, we proposed a proof of concept for a combination of RESESOP-
Kaczmarz, landmark detection and filtered backprojection, i.e., a hybrid
reconstruction technique that incorporates classical as well as data-driven
methods for computerized tomography, to image an object undergoing an affine
linear motion with constant speed during data acquisition. The motion is not
known a priori, but estimated during the reconstruction. In addition, our
combined method leads to a faster and more accurate solution of dynamic
computerized tomography provided the respective model inexactness, stemming
from working with the static model, is accessible and if the landmark
detection works well.

By consequence, future research is dedicated to estimating the inexactness
directly from the data to provide a realistic initial setting for our
pipeline. Also, the landmark detection can be generalized towards more complex
objects

In the future, it will be necessary to find ways to estimate the modeling
inexactness given only the dynamic sinogram. Furthermore, it would be
interesting to train neural networks performing landmark detection for a
broader class of objects.

Our method could also be extended to more complex motions by interpolating the
motion piecewise by affine linear motions with constant speed. To achieve
this, it would be possible to estimate the inexactness at multiple different
time points.

Acknowledgements: G.S. and A.W. acknowledge support by the German Research
Foundation (DFG; Project-ID 432680300 ? SFB 1456), and thank Bernadette Hahn-
Rigaud and Alice Oberacker for fruitful discussions and help with the
implementation.

## References

  * [1] M. Abramowith and I. Stegun, Handbook of Mathematical Functions, Dover, New York, 1972.
  * [2] Simon Arridge, Pascal Fernsel, and Andreas Hauptmann, Joint reconstruction and low-rank decomposition for dynamic inverse problems, Inverse Problems and Imaging, 16 (2022), pp. 483–523.
  * [3] S. Blanke, B. Hahn, and A. Wald, Inverse problems with inexact forward operator: iterative regularization and application in dynmaic imaging, (2020).
  * [4] S. Bonnet, A. Koenig, S. Roux, P. Hugonnard, R. Guillemaud, and P. Grangeat, Dynamic x-ray computed tomography, Proceedings of the IEEE, 91 (2003), pp. 1574–1587.
  * [5] Tatiana A. Bubba, Maximilian März, Z. Purisha, Matti Lassas, and Samuli Siltanen, Shearlet-based regularization in sparse dynamic tomography, in Optical Engineering + Applications, 2017.
  * [6] Martin Burger, Hendrik Dirks, Lena Frerking, Andreas Hauptmann, Tapio Helin, and Samuli Siltanen, A variational reconstruction method for undersampled dynamic x-ray tomography based on physical motion models, Inverse Problems, 33 (2017), p. 124008.
  * [7] Chong Chen, Barbara Gris, and Ozan Öktem, A new variational model for joint image reconstruction and motion estimation in spatiotemporal imaging, SIAM Journal on Imaging Sciences, 12 (2019), pp. 1686–1719.
  * [8] Zhen-Hua Feng, Josef Kittler, Muhammad Awais, Patrik Huber, and Xiao-Jun Wu, Wing loss for robust facial landmark localisation with convolutional neural networks, 06 2018.
  * [9] Bernadette Hahn, Rekonstruktionsverfahren in der dynamischen Computertomographie, PhD thesis, 01 2013.
  * [10] Bernadette Hahn, Reconstruction of dynamic objects with affine deformations in computerized tomography, Journal of Inverse and Ill-posed Problems, 22 (2014), pp. 323–339.
  * [11] Bernadette Hahn, Motion estimation and compensation strategies in dynamic computerized tomography, Sensing and Imaging, 18 (2017).
  * [12] B. N. Hahn, Motion Compensation Strategies in Tomography, Springer International Publishing, Cham, 2021, pp. 51–83.
  * [13] A. Hölzing, Philipp Stahlhut, K. Dremel, Jonas Graetz Dittmann, J. Engel, Simon Zabler, and Randolf Hanke, First results on laboratory nano-ct with a needle reflection target and an adapted toolchain, 10 2016.
  * [14] A. Kak and M. Slaney, The Principles of Computerized Tomographic Imaging, 2001.
  * [15] Emma Kamutta, Sofia Mäkinen, and Alexander Meaney, Cone-beam computed tomography dataset of a seashell, Aug. 2022.
  * [16] A Katsevich, An accurate approximate algorithm for motion compensation in two-dimensional tomography, Inverse Problems, 26 (2010), p. 065007.
  * [17] D. Kingma and Jimmy Ba, Adam: A method for stochastic optimization, CoRR, abs/1412.6980 (2015).
  * [18] Alfred Louis, Diffusion reconstruction from very noisy tomographic data, Inverse Problems and Imaging, 4 (2010).
  * [19] Weiguo Lu and Thomas R Mackie, Tomographic motion detection and correction directly in sinogram space, Physics in Medicine & Biology, 47 (2002), p. 1267.
  * [20] T. Lütjen, Fabian Schönfeld, A. Oberacker, J. Leuschner, A. Wald, and T. Kluth, Learning-based approaches for reconstructions with inexact operators in nanoct applications, 2023.
  * [21] A. Meaney, HelTomo - Helsinki Tomography Toolbox, Matlab Central File Exchange, 2023-09-25.
  * [22] Jennifer Mueller and Samuli Siltanen, Linear and nonlinear inverse problems with practical applications, 08 2013.
  * [23] F. Natterer, The Mathematics of Computerized Tomography, Vieweg+Teubner Verlag, 1986.
  * [24] F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction, siam, 2001.
  * [25] F. Schöpfer and T. Schuster, Fast regularizing sequential subspace optimization in banach spaces, Inverse Problems 25, (2008).
  * [26] F. Schöpfer, T. Schuster, and A. Louis, Metric and bregman projections onto affine subspaces and their computation via sequential subspace optimization methods, J. Inverse ILL-Posed Problems 15, (2007).
  * [27] Thomas Schuster, The Method of Approximate Inverse: Theory and Applications, vol. 1906, 01 2007.
  * [28] T. Schuster, B. Kaltenbacher, B. Hofmann, and K. Kazimierski, Regularization Methods in Banach Spaces, Berlin: de Gruyter & Co, 2012.
  * [29] J. Anthony Seibert, John M. Boone, and Karen K. Lindfors, Flat-field correction technique for digital detectors, in Medical Imaging, 1998.
  * [30] A. Wald, A fast subspace optimization method for nonlinear inverse problems in Banach spaces with an application in parameter identification, Inverse Problems, 34 (2018), p. 085008.