# Unveiling the Fingerprint of Eccentric Binary Black Hole Mergers

Hao Wang<EMAIL_ADDRESS>Department of Astronomy, School of Physics,
Huazhong University of Science and Technology, Wuhan 430074, China  Yuan-Chuan
Zou<EMAIL_ADDRESS>Department of Astronomy, School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China  Qing-Wen Wu
<EMAIL_ADDRESS>Department of Astronomy, School of Physics, Huazhong
University of Science and Technology, Wuhan 430074, China  Yu Liu
<EMAIL_ADDRESS>State Key Laboratory of Public Big Data, Guizhou University,
Guiyang 550025, China

###### Abstract

The orbital eccentricity plays a crucial role in shaping the dynamics of
binary black hole (BBH) mergers. Remarkably, our recent findings reveal a
universal oscillation in essential dynamic quantities: peak luminosity
$L_{\text{peak}}$, masses $M_{f}$, spins $\alpha_{f}$, and recoil velocity
$V_{f}$ of the final remnant black hole, as the initial eccentricity $e_{0}$
undergoes variation. In this letter, by leveraging RIT’s extensive numerical
relativistic simulations of nonspinning eccentric orbital BBH mergers, we not
only confirm the universal oscillation in peak amplitudes (including harmonic
modes), similar to the oscillations observed in $L_{\text{peak}}$, $M_{f}$,
$\alpha_{f}$, and $V_{f}$, but also make the first discovery of a ubiquitous
spiral-like internal fine structure that correlates $L_{\text{peak}}$,
$M_{f}$, $\alpha_{f}$, $V_{f}$, and peak amplitudes. This distinctive feature,
which we term the “fingerprint” of eccentric orbital BBH mergers, carries
important implications for unraveling the intricate dynamics and astrophysics
associated with eccentric orbital BBH mergers.

††preprint: APS/123-QED††preprint: APS/123-QED

##  I Introduction

Following the groundbreaking detection of the gravitational wave (GW) event
GW150914 [1], the field of GW detection has witnessed a remarkable evolution,
transforming into routine practice. Concurrently, numerical relativity (NR),
our trusted tool for investigating the dynamics of BBH mergers, has made
significant strides in exploring the vast parameter space of BBH systems since
its initial breakthrough [2, 3, 4], encompassing configurations spinless
systems, spin alignment, spin precession, eccentric orbits and even extreme
mass ratios.

Although most of the research on NR and GW detection has focused predominantly
on circular orbits, attributed to the circularizing effect of GW radiation [5,
6], it is important to recognize various mechanisms through which BBHs can
exhibit nonzero eccentricity before their merger. These mechanisms include
double-single interactions [7, 8], double-double interactions [9, 10], and
gravitational capture [11, 12] within dense stellar environments such as
globular clusters [13, 14, 15, 16, 17, 18] and galactic nuclei [11, 19, 20,
21, 22]. Notably, in three-body systems [23] involving binary objects orbiting
a supermassive black hole, the eccentricity of the inner binary can experience
oscillations due to the Kozai-Lidov mechanism [24, 25, 26, 27]. These
eccentric BBH systems become detectable once they enter the frequency band of
ground-based GW detectors such as LIGO [28], VIRGO [29] and KAGRA [30]. A
notable example is GW190521 [31], considered a potential BBH merger with a
high eccentricity of $e=0.69_{-0.22}^{+0.17}$ [32, 33]. With continuous
advancements in detector sensitivity, future ground-based GW detectors such as
the Einstein Telescope [34] or Cosmic Explorer [35] are anticipated to observe
an increasing number of eccentric BBH mergers.

Over the past decades, several collaborations in NR have conducted extensive
simulations of binary compact object mergers, including SXS [36, 37], RIT [38,
39, 40, 41], and Georgia Tech. [42, 43]. These simulations have yielded
significant progress in modeling dynamic quantities that hold great
astrophysical significance, such as peak luminosity, recoil velocity, remnant
mass and spin. Various methods have been employed for modeling these dynamic
quantities, including Gaussian Process Regression [44, 45, 46], post-Newtonian
(PN) approaches [47, 48, 49], effective one body methods [50], and direct
fitting of formulas with NR data [51, 52, 53, 54, 55, 56, 57, 58]. These
approaches enable the modeling of dynamic quantities for both quasi-circular
and eccentric orbits, based on initial parameters such as mass ratio and spin.
There have been some investigations of eccentric BBH mergers in recent
decades, including studies on the influence of eccentricity on recoil velocity
from a PN perspective [59], the transition from inspiral to plunge in
eccentric orbits [60], orbital circularization for eccentric orbits [61],
remnant properties in low eccentricity orbits using NR [62], kick enhancement
caused by eccentricity [63], and anomalies in recoil due to eccentricity [64].
Recently, a few articles have explored correlations between the dynamic
quantities. For instance, Ref. [65] examined the correlation between the peak
amplitude and the remnant spin, while Ref. [41] utilized analytical formulas
to directly fit NR data for correlations. Furthermore, Ref. [66] obtained a
gauge-independent correlation fitted by polynomials, etc. However, these
methods did not fully consider the potential fine structures present in NR
data. In our recent work [67], incorporating extensive simulation results of
BBH mergers in eccentric orbits from RIT’s fourth release, we first identified
universal oscillations in variations of the aforementioned dynamic quantities
as a function of initial eccentricity, which indicate the presence of internal
structures in these quantities. In this Letter, we present additional findings
of oscillations in peak amplitudes (including higher-order harmonic modes). By
summarizing the correlations between quantities such as peak luminosity
$L_{\text{peak}}$ (maximum value of dimensionalized radiation energy [51]),
masses $M_{f}$, spins $\alpha_{f}$, recoil velocity $V_{f}$ (magnitude of
dimensional recoil velocity), and peak amplitudes (maximum amplitude values)
$\mathcal{A}_{22,\text{peak}}$, $\mathcal{A}_{32,\text{peak}}$,
$\mathcal{A}_{44,\text{peak}}$, we unveil, for the first time, the existence
of a spiral-like internal fine structure in these correlations. We refer to
this structure as the “fingerprint” of eccentric orbital BBH mergers according
to its characteristics.

Throughout this letter, we adopt geometric units where $G=c=1$. The component
masses of BBH are represented as $m_{1}$ and $m_{2}$, while the total mass is
denoted by $M=m_{1}+m_{2}$. For simplicity, we set the total mass $M$ at unity
(sometimes explicitly writing it for clarity). The mass ratio $q$ is defined
as $q=m_{1}/m_{2}$, and $m_{1}<m_{2}$.

##  II Method

We utilize NR simulations of BBH systems in eccentric orbits obtained from the
Rochester Institute of Technology (RIT) catalog [68]. These simulations were
conducted using the LazEv code [69], implemented within the Einstein Toolkit
[70] alongside the CACTUS/CARPET infrastructure [71]. The LazEv code employs
the moving puncture approach [3] and utilizes the BSSNOK formalism for
evolution systems [72, 73, 74]. In the initial stages, RIT employed
AHFinderDirect [75] to locate apparent horizons and employed the isolated
horizon algorithm to measure the amplitude of the horizon spins [76], denoted
as $S_{H}$. Subsequently, they calculated the horizon mass using the
Christodoulou formula:
$m_{H}=\sqrt{m_{\mathrm{irr}}^{2}+S_{H}^{2}/\left(4m_{\mathrm{irr}}^{2}\right)}$.
Here, $m_{\mathrm{irr}}$ represents the irreducible mass, defined as
$m_{\mathrm{irr}}=\sqrt{A_{H}/(16\pi)}$, where $A_{H}$ corresponds to the
surface area of the horizon [76].

In generating the initial data, RIT adopts the puncture approach [77] in
combination with the TwoPunctures code [78]. To enable continuous-eccentricity
simulations, RIT initially employs PN techniques, as outlined in Ref. [79], to
generate initial data for quasi-circular orbits. Subsequently, by introducing
a new parameter $\epsilon$ within the range of 0 to 1, the tangential linear
momentum is modified according to $p_{t}=p_{t,qc}(1-\epsilon)$. Within this
framework, the initial positions of the BBHs are fixed at the apocenter, and
the eccentricity of the orbit gradually increases throughout the simulations,
spanning from the quasi-circular orbit ($e_{0}=0$) to the head-on collision
limit ($e_{0}=1$). The initial eccentricity $e_{0}$ of the orbit can be
approximated by $e_{0}=2\epsilon-\epsilon^{2}$, offering a second-order
approximation in terms of $\epsilon$ that accurately captures the limits of
$e_{0}=0$ and $e_{0}=1$ at $\epsilon=0$ and $\epsilon=1$, respectively [41].

The RIT catalog provides a comprehensive dataset comprising waveform data and
accompanying metadata [68]. In our research, we utilize the gravitational wave
strain $h$ as the waveform data, which can be obtained as the harmonic mode
$rh=\sum_{l,m}rh_{lm}{}_{-2}Y_{l,m}(\theta,\phi)$, where $r$ represents the
extracted radius, and ${}_{-2}Y_{l,m}(\theta,\phi)$ denotes the spin-weighted
spherical harmonic functions. The harmonic mode $h_{lm}$ can be decomposed
into amplitude and phase as
$h_{lm}=\mathcal{A}_{lm}(t)\exp\left[-i\Phi_{lm}(t)\right]$. The metadata
within the catalog provides essential information regarding the initial data
of the simulations, which includes details such as the mass ratio, initial
distance, initial linear momentum, and more. Additionally, the metadata
contain significant simulation results, such as the peak luminosity
$L_{\text{peak}}$, masses $M_{f}$, spins $\alpha_{f}$, and recoil velocity
$V_{f}$ of the final remnant black hole. In all simulations conducted by RIT,
it has been ensured that the waveforms, at the resolutions provided in the
catalog, have achieved convergence up to fourth order with resolution [41].
The evaluation of quantities related to the black hole horizon, such as the
final masses $M_{f}$ and spins $\alpha_{f}$ of the remnant, yields errors of
the order of 0.1% using the isolated horizon algorithm [76]. Furthermore, the
radiative computed quantities, including the recoil velocities $V_{f}$ and the
peak luminosities $L_{\text{peak}}$, are evaluated with a typical error of 5%
[41]. Therefore, all of the data provided in the catalog meet the precision
requirements for research purposes.

In this letter, our focus is solely on the nonspinning configuration, as it
provides a larger dataset compared to spin alignment and spin precession
simulations [67]. In addition to the previously mentioned dynamic quantities
$L_{\text{peak}}$, $M_{f}$, $\alpha_{f}$, and $V_{f}$, we also introduce the
peak amplitudes of different harmonic modes, specifically
$\mathcal{A}_{22,\text{peak}}$, $\mathcal{A}_{32,\text{peak}}$, and
$\mathcal{A}_{44,\text{peak}}$. We have chosen not to investigate other
harmonic modes and peak frequencies due to issues with the former in RIT’s
catalog (although we speculate that other harmonics behave similarly) and the
fact that the latter does not accurately represent the merger structure of BBH
systems. Our analysis focuses on characterizing the variations of the peak
amplitudes as functions of $e_{0}$ and examining the correlations between the
dynamic quantities and the peak amplitudes. Fig. 1 presents the parameter
space of the eccentric nonspinning BBH NR simulations used in our study.
Specifically, we analyze a total of 192 data sets for initial coordinate
separation $D_{\text{ini}}=11.3M$ and mass ratios $q=1,0.75,0.5,0.25$, along
with 48 data sets for $D_{\text{ini}}=24.6M$ and mass ratio $q=1$. Due to the
limited availability of simulated data points, nonspinning simulations for
other mass ratios, as well as simulations involving spin alignment and spin
precession, cannot provide sufficient information for our research [67].

Figure 1: Parameter space of the eccentric nonspinning BBH NR simulations used
in our study. There are 192 sets simulations with $D_{\text{ini}}=11.3M$ and
$q=1$ (43 sets), $q=0.25$ (67 sets), $q=0.5$ (41 sets), $q=0.75$ (41 sets),
and 48 sets simulations with $D_{\text{ini}}=24.6M$ and $q=1$.

##  III Results

In our previous study, described in Ref. [67], we investigate the universal
oscillations observed in dynamical quantities $L_{\text{peak}}$, $M_{f}$,
$\alpha_{f}$, and $V_{f}$ as a function of $e_{0}$. Supplementary to these
findings, Fig. 2 presents the relationship between the peak amplitude of the
2-2 mode $\mathcal{A}_{22,\text{peak}}$ and $e_{0}$ for both
$D_{\text{ini}}=11.3M$ and $D_{\text{ini}}=24.6M$ (for
$\mathcal{A}_{32,\text{peak}}$, $\mathcal{A}_{44,\text{peak}}$, refer to Fig.
S1 of supplementary materials). Each point in Fig. 2 corresponds to a
simulation result. Fig. 2 reveal the same universal oscillatory behavior
observed in the peak amplitudes. While these oscillations share similarities
with peak luminosity $L_{\text{peak}}$ oscillations described in Ref. [67],
they also demonstrate distinct characteristics influenced by mass ratio $q$
when considering higher harmonic modes. In our previous work [67], we proposed
that these peculiar oscillations may arise from orbital transitions based on
the integer orbital cycles $N_{\text{waves }}=\frac{\Delta\Phi}{4\pi}$ of GW.
Here, the phase difference $\Delta\Phi$ is calculated as
$\Delta\Phi=\Phi\left(t_{\mathrm{merger}}\right)-\Phi\left(t_{0}+t_{\mathrm{relax}}\right)$
from the 2-2 mode of GW, where $t_{\mathrm{merger}}$, $t_{0}$, and
$t_{\mathrm{relax}}$ represent the merger moment, initial moment and duration
of the junk radiation, respectively. These fine oscillatory structures can
only be discerned when there is a sufficient number of data points from
eccentric numerical simulations [67]. In Fig. S2 of supplementary material, we
provide further evidence to emphasize their relevance by illustrating the
relationship between peak amplitudes and the integer $N_{\text{waves }}$ as
described in Ref. [67]. The shift and enhancement of the oscillations for
higher $e_{0}$ for $D_{\text{ini}}=24.6M$ compared to $D_{\text{ini}}=11.3M$,
along with the correspondence between integer $N_{\text{waves }}$ and the peak
and valley points, indicates that these oscillations originate from strong
field dynamics, as highlighted in Ref. [67]. Consequently, we aim to explore
the impact of this strong field effect due to the existence of eccentricity on
the correlations between the dynamic quantities $L_{\text{peak}}$, $M_{f}$,
$\alpha_{f}$ and $V_{f}$, and peak amplitudes $\mathcal{A}_{22,\text{peak}}$,
$\mathcal{A}_{32,\text{peak}}$, $\mathcal{A}_{44,\text{peak}}$. Fig. 3
showcases three representative correlations between the aforementioned dynamic
quantities and peak amplitudes. To maintain brevity, we include the remaining
18 correlations in Figs. S3, S4 and S5 of supplementary material.

Figure 2: Variations of 2-2 mode peak amplitude $\mathcal{A}_{22,\text{peak}}$
as a function of initial eccentricity $e_{0}$ at $D_{\text{ini}}=11.3M$ and
$D_{\text{ini}}=24.6M$ for nonspinning configuration with different mass
ratios. Figure 3: Correlations between quantities $L_{\text{peak}}$, $M_{f}$,
and $\alpha_{f}$, for four mass ratios and two initial coordinate separations
$D_{\text{ini}}=11.3M$ and $D_{\text{ini}}=24.6M$. The “fingerprints” are
drawn in insets as partial enlargement for details.

In Fig. 3, we present the functional relationships between the quantities
$L_{\text{peak}}$, $M_{f}$, and $\alpha_{f}$ for four mass ratios and two
initial coordinate separations. Partial enlargements are included as insets
for a detailed examination. Taking panels (a) and (b) in Fig. 3 as examples,
we observe a captivating spiral structure in the curves. The center of the
spiral corresponds to the quasi-circular orbit counterpart ($e_{0}=0$) of this
family of eccentric orbital BBH merger simulations. As the initial
eccentricity $e_{0}$ gradually increases, the spiral rotates outwards. When
the spiral reaches the transition point from inspiral to plunge (the maximum
or minimum value of oscillations, see Ref. [67]), it terminates, and the curve
gradually reaches the specific value corresponding to the maximum initial
eccentricity (i.e., $e_{0}=1$ for the head-on collision limit). The insets (1)
and (2) in panel (a) represent the mass ratios $q=0.5$ (40 data points) and
$q=0.25$ (67 data points), respectively, for $D_{\text{ini}}=11.3M$. From
panel (a) of Fig. 3, we observe that a smaller mass ratio $q$ corresponds to a
smaller spiral, originating weaker oscillations of the corresponding quantity
with changes in $e_{0}$, as observed in Ref. [67]. The right panel (b)
corresponds to the case of the same mass ratio $q=1$ and different initial
separations $D_{\text{ini}}=11.3M$ and $D_{\text{ini}}=24.6M$. We can discern
that a larger $D_{\text{ini}}$ corresponds to a larger spiral, originating
stronger oscillations of the corresponding quantity with changes in $e_{0}$.
In panels (c) and (d), the spiral becomes more elliptical; however, we still
observe a similar spiral structure. The elliptical shape is due to the weaker
oscillation of $\alpha_{f}$ compared to $L_{\text{peak}}$ and $M_{f}$ (see
Ref. [67]). In Ref. [67], we obtained irregular oscillations in $V_{f}$ and
provided an explanation. When any other quantity is combined with $V_{f}$, an
irregular spiral pattern emerges (please refer to the supplementary material).
These curves of $V_{f}$ exhibit various irregular behaviors and overlaps,
distinguishing them from panels (a), (b), (c), (d). However, it is important
to acknowledge that, despite their irregularity, their internal structure
remains intact. Based on the supplementary material and Fig. 3, we can
conclude that a ubiquitous spiral structure is present in the correlations
between the dynamic quantities $L_{\text{peak}}$, $M_{f}$, $\alpha_{f}$,
$V_{f}$, and peak amplitudes. The position and size of a spiral are
characterized by the mass ratio $q$ and initial coordinate separation
$D_{\text{ini}}$. Two types of spiral structures are identified: one is
irregular, corresponding to the irregularly oscillating $V_{f}$, while the
other is regular and not combined with $V_{f}$. Due to the spiral structure
observed in these correlations, we refer to the prevalent correlations in
eccentric orbital BBH mergers as “fingerprint”, which uniquely represents a
family of eccentric BBH merger simulations, analogous to how a fingerprint
represents an individual.

In some cases, the orbital transition points of the integer $N_{\text{waves}}$
may not precisely align with the peaks and valleys of oscillation. Several
factors contribute to this discrepancy, including errors and pericenter
precession [67]. To investigate the relationship between the spirals and the
peaks and valleys of oscillation, we mark the points corresponding to the
peaks and valleys in the relationship between the dynamic quantities
$L_{\text{peak}}$, $M_{f}$, $\alpha_{f}$, and $e_{0}$ on the spirals in Fig.
S6 of supplementary material. Notably, we find that the peaks and valleys of
the oscillation are positioned on opposite sides of the spiral and are
concentrated near a straight line. This phenomenon suggests that the spiral
structure arises from the oscillations of $L_{\text{peak}}$, $M_{f}$, and
$\alpha_{f}$, which in turn originate from the orbital transitions of the
periodic process.

##  IV Discussion

Previously, many studies have approached modeling dynamic quantities by
directly fitting scatter points using methods such as polynomials. However,
the existence of internal fine structures within these dynamic quantities,
particularly in relation to eccentricity, was not recognized. As we discussed
previously, these internal structures originate from the oscillations of
dynamic quantities as a function of eccentricity, which are driven by the
strong-field dynamics of eccentric BBH mergers, specifically, orbital
transitions that correspond to different orbital cycle numbers. This unique
characteristic distinguishes BBH mergers in eccentric orbits from those in
quasi-circular orbits, where the spiral structure dissipates and returns to
the central point.

In summary, by leveraging RIT’s extensive numerical relativity simulations of
nonspinning eccentric orbital binary black hole mergers, we make the first
discovery of a universal spiral-like internal fine structure in the
correlations between various quantities, including $L_{\text{peak}}$, $M_{f}$,
$\alpha_{f}$, and peak amplitudes such as $\mathcal{A}_{22,\text{peak}}$,
$\mathcal{A}_{32,\text{peak}}$, and $\mathcal{A}_{44,\text{peak}}$, arising
due to eccentricity. This spiral structure presents a fresh perspective on the
investigation of eccentric orbital BBH mergers and their strong-field
dynamics. Furthermore, it holds significant implications for the study of
eccentric orbital BBH mergers within dense stellar environments like globular
star clusters and galactic nuclei [8, 15]. Moreover, it enables a more precise
understanding and modeling of various astrophysical quantities, including
remnant black hole mass, spin and recoil velocity.

It is important to note that our examination of nonspinning eccentric BBHs
represents only a subset of the parameter space explored by NR simulations. In
the future, we anticipate the generation of additional simulations that
encompass spin-aligned and spin-precessing configurations, which will offer
crucial insights into the interplay between spin and eccentricity in strong-
field dynamics.

###### Acknowledgements.

The authors are very grateful to the RIT collaboration for the numerical
simulation of eccentric BBH mergers, and thanks to Yan-Fang Huang, Zhou-Jian
Cao for their helpful discussions. The computation is partially completed in
the HPC Platform of Huazhong University of Science and Technology. The
languages was polished by ChatGPT during the revision of the draft. This work
is supported by the National Key R&D Program of China (2021YFA0718504).

## References

  * Abbott  _et al._ [2016] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Observation of Gravitational Waves from a Binary Black Hole Merger, Phys. Rev. Lett. 116, 061102 (2016), arXiv:1602.03837 [gr-qc] .
  * Pretorius [2005] F. Pretorius, Evolution of binary black hole spacetimes, Phys. Rev. Lett. 95, 121101 (2005), arXiv:gr-qc/0507014 .
  * Campanelli  _et al._ [2006] M. Campanelli, C. O. Lousto, P. Marronetti, and Y. Zlochower, Accurate evolutions of orbiting black-hole binaries without excision, Phys. Rev. Lett. 96, 111101 (2006), arXiv:gr-qc/0511048 .
  * Baker  _et al._ [2006] J. G. Baker, J. Centrella, D.-I. Choi, M. Koppitz, and J. van Meter, Gravitational wave extraction from an inspiraling configuration of merging black holes, Phys. Rev. Lett. 96, 111102 (2006), arXiv:gr-qc/0511103 .
  * Peters and Mathews [1963] P. C. Peters and J. Mathews, Gravitational radiation from point masses in a Keplerian orbit, Phys. Rev. 131, 435 (1963).
  * Peters [1964] P. C. Peters, Gravitational Radiation and the Motion of Two Point Masses, Phys. Rev. 136, B1224 (1964).
  * Samsing  _et al._ [2014] J. Samsing, M. MacLeod, and E. Ramirez-Ruiz, The Formation of Eccentric Compact Binary Inspirals and the Role of Gravitational Wave Emission in Binary-Single Stellar Encounters, Astrophys. J. 784, 71 (2014), arXiv:1308.2964 [astro-ph.HE] .
  * Samsing  _et al._ [2018] J. Samsing, A. Askar, and M. Giersz, MOCCA-SURVEY Database. I. Eccentric Black Hole Mergers during Binary–Single Interactions in Globular Clusters, Astrophys. J. 855, 124 (2018), arXiv:1712.06186 [astro-ph.HE] .
  * Zevin  _et al._ [2019] M. Zevin, J. Samsing, C. Rodriguez, C.-J. Haster, and E. Ramirez-Ruiz, Eccentric Black Hole Mergers in Dense Star Clusters: The Role of Binary–Binary Encounters, Astrophys. J. 871, 91 (2019), arXiv:1810.00901 [astro-ph.HE] .
  * Arca-Sedda  _et al._ [2021] M. Arca-Sedda, G. Li, and B. Kocsis, Order in the chaos - Eccentric black hole binary mergers in triples formed via strong binary–binary scatterings, Astron. Astrophys. 650, A189 (2021), arXiv:1805.06458 [astro-ph.HE] .
  * Gondán and Kocsis [2021] L. Gondán and B. Kocsis, High eccentricities and high masses characterize gravitational-wave captures in galactic nuclei as seen by Earth-based detectors, Mon. Not. Roy. Astron. Soc. 506, 1665 (2021), arXiv:2011.02507 [astro-ph.HE] .
  * East  _et al._ [2013] W. E. East, S. T. McWilliams, J. Levin, and F. Pretorius, Observing complete gravitational wave signals from dynamical capture binaries, Phys. Rev. D 87, 043004 (2013), arXiv:1212.0837 [gr-qc] .
  * O’Leary  _et al._ [2006] R. M. O’Leary, F. A. Rasio, J. M. Fregeau, N. Ivanova, and R. W. O’Shaughnessy, Binary mergers and growth of black holes in dense star clusters, Astrophys. J. 637, 937 (2006), arXiv:astro-ph/0508224 .
  * Rodriguez  _et al._ [2015] C. L. Rodriguez, M. Morscher, B. Pattabiraman, S. Chatterjee, C.-J. Haster, and F. A. Rasio, Binary Black Hole Mergers from Globular Clusters: Implications for Advanced LIGO, Phys. Rev. Lett. 115, 051101 (2015), [Erratum: Phys.Rev.Lett. 116, 029901 (2016)], arXiv:1505.00792 [astro-ph.HE] .
  * Samsing [2018] J. Samsing, Eccentric Black Hole Mergers Forming in Globular Clusters, Phys. Rev. D 97, 103014 (2018), arXiv:1711.07452 [astro-ph.HE] .
  * Rodriguez  _et al._ [2018a] C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, and F. A. Rasio, Post-Newtonian Dynamics in Dense Star Clusters: Highly-Eccentric, Highly-Spinning, and Repeated Binary Black Hole Mergers, Phys. Rev. Lett. 120, 151101 (2018a), arXiv:1712.04937 [astro-ph.HE] .
  * Rodriguez  _et al._ [2018b] C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, K. Kremer, F. A. Rasio, J. Samsing, C. S. Ye, and M. Zevin, Post-Newtonian Dynamics in Dense Star Clusters: Formation, Masses, and Merger Rates of Highly-Eccentric Black Hole Binaries, Phys. Rev. D 98, 123005 (2018b), arXiv:1811.04926 [astro-ph.HE] .
  * Park  _et al._ [2017] D. Park, C. Kim, H. M. Lee, Y.-B. Bae, and K. Belczynski, Black Hole Binaries Dynamically Formed in Globular Clusters, Mon. Not. Roy. Astron. Soc. 469, 4665 (2017), arXiv:1703.01568 [astro-ph.HE] .
  * Hoang  _et al._ [2018] B.-M. Hoang, S. Naoz, B. Kocsis, F. A. Rasio, and F. Dosopoulou, Black Hole Mergers in Galactic Nuclei Induced by the Eccentric Kozai–Lidov Effect, Astrophys. J. 856, 140 (2018), arXiv:1706.09896 [astro-ph.HE] .
  * Gondán  _et al._ [2018] L. Gondán, B. Kocsis, P. Raffai, and Z. Frei, Eccentric Black Hole Gravitational-Wave Capture Sources in Galactic Nuclei: Distribution of Binary Parameters, Astrophys. J. 860, 5 (2018), arXiv:1711.09989 [astro-ph.HE] .
  * Samsing  _et al._ [2022] J. Samsing, I. Bartos, D. J. D’Orazio, Z. Haiman, B. Kocsis, N. W. C. Leigh, B. Liu, M. E. Pessah, and H. Tagawa, AGN as potential factories for eccentric black hole mergers, Nature 603, 237 (2022), arXiv:2010.09765 [astro-ph.HE] .
  * Tagawa  _et al._ [2021] H. Tagawa, B. Kocsis, Z. Haiman, I. Bartos, K. Omukai, and J. Samsing, Eccentric Black Hole Mergers in Active Galactic Nuclei, Astrophys. J. Lett. 907, L20 (2021), arXiv:2010.10526 [astro-ph.HE] .
  * Naoz  _et al._ [2013] S. Naoz, B. Kocsis, A. Loeb, and N. Yunes, Resonant Post-Newtonian Eccentricity Excitation in Hierarchical Three-body Systems, Astrophys. J. 773, 187 (2013), arXiv:1206.4316 [astro-ph.SR] .
  * Silsbee and Tremaine [2017] K. Silsbee and S. Tremaine, Lidov-Kozai Cycles with Gravitational Radiation: Merging Black Holes in Isolated Triple Systems, Astrophys. J. 836, 39 (2017), arXiv:1608.07642 [astro-ph.HE] .
  * Blaes  _et al._ [2002] O. Blaes, M. H. Lee, and A. Socrates, The kozai mechanism and the evolution of binary supermassive black holes, Astrophys. J. 578, 775 (2002), arXiv:astro-ph/0203370 .
  * Antognini  _et al._ [2014] J. M. Antognini, B. J. Shappee, T. A. Thompson, and P. Amaro-Seoane, Rapid Eccentricity Oscillations and the Mergers of Compact Objects in Hierarchical Triples, Mon. Not. Roy. Astron. Soc. 439, 1079 (2014), arXiv:1308.5682 [astro-ph.HE] .
  * Katz  _et al._ [2011] B. Katz, S. Dong, and R. Malhotra, Long-Term Cycling of Kozai-Lidov Cycles: Extreme Eccentricities and Inclinations Excited by a Distant Eccentric Perturber, Phys. Rev. Lett. 107, 181101 (2011), arXiv:1106.3340 [astro-ph.EP] .
  * Aasi  _et al._ [2015] J. Aasi _et al._ (LIGO Scientific, VIRGO), Characterization of the LIGO detectors during their sixth science run, Class. Quant. Grav. 32, 115012 (2015), arXiv:1410.7764 [gr-qc] .
  * Acernese  _et al._ [2015] F. Acernese _et al._ (VIRGO), Advanced Virgo: a second-generation interferometric gravitational wave detector, Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 [gr-qc] .
  * Akutsu  _et al._ [2019] T. Akutsu _et al._ (KAGRA), KAGRA: 2.5 Generation Interferometric Gravitational Wave Detector, Nature Astron. 3, 35 (2019), arXiv:1811.08079 [gr-qc] .
  * Abbott  _et al._ [2020] R. Abbott _et al._ (LIGO Scientific, Virgo), GW190521: A Binary Black Hole Merger with a Total Mass of $150M_{\odot}$, Phys. Rev. Lett. 125, 101102 (2020), arXiv:2009.01075 [gr-qc] .
  * Gayathri  _et al._ [2022] V. Gayathri, J. Healy, J. Lange, B. O’Brien, M. Szczepanczyk, I. Bartos, M. Campanelli, S. Klimenko, C. O. Lousto, and R. O’Shaughnessy, Eccentricity estimate for black hole mergers with numerical relativity simulations, Nature Astron. 6, 344 (2022), arXiv:2009.05461 [astro-ph.HE] .
  * Romero-Shaw  _et al._ [2020] I. M. Romero-Shaw, P. D. Lasky, E. Thrane, and J. C. Bustillo, GW190521: orbital eccentricity and signatures of dynamical formation in a binary black hole merger signal, Astrophys. J. Lett. 903, L5 (2020), arXiv:2009.04771 [astro-ph.HE] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , The Einstein Telescope: A third-generation gravitational wave observatory, Class. Quant. Grav. 27, 194002 (2010).
  * Reitze  _et al._ [2019] D. Reitze _et al._ , Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO, Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Mroue  _et al._ [2013] A. H. Mroue _et al._ , Catalog of 174 Binary Black Hole Simulations for Gravitational Wave Astronomy, Phys. Rev. Lett. 111, 241104 (2013), arXiv:1304.6077 [gr-qc] .
  * Boyle  _et al._ [2019] M. Boyle _et al._ , The SXS Collaboration catalog of binary black hole simulations, Class. Quant. Grav. 36, 195006 (2019), arXiv:1904.04831 [gr-qc] .
  * Healy  _et al._ [2017a] J. Healy, C. O. Lousto, Y. Zlochower, and M. Campanelli, The RIT binary black hole simulations catalog, Class. Quant. Grav. 34, 224001 (2017a), arXiv:1703.03423 [gr-qc] .
  * Healy  _et al._ [2019] J. Healy, C. O. Lousto, J. Lange, R. O’Shaughnessy, Y. Zlochower, and M. Campanelli, Second RIT binary black hole simulations catalog and its application to gravitational waves parameter estimation, Phys. Rev. D 100, 024021 (2019), arXiv:1901.02553 [gr-qc] .
  * Healy and Lousto [2020] J. Healy and C. O. Lousto, Third RIT binary black hole simulations catalog, Phys. Rev. D 102, 104018 (2020), arXiv:2007.07910 [gr-qc] .
  * Healy and Lousto [2022] J. Healy and C. O. Lousto, Fourth RIT binary black hole simulations catalog: Extension to eccentric orbits, Phys. Rev. D 105, 124010 (2022), arXiv:2202.00018 [gr-qc] .
  * Jani  _et al._ [2016] K. Jani, J. Healy, J. A. Clark, L. London, P. Laguna, and D. Shoemaker, Georgia Tech Catalog of Gravitational Waveforms, Class. Quant. Grav. 33, 204001 (2016), arXiv:1605.03204 [gr-qc] .
  * Ferguson  _et al._ [2023] D. Ferguson _et al._ , Second MAYA Catalog of Binary Black Hole Numerical Relativity Waveforms, (2023), arXiv:2309.00262 [gr-qc] .
  * Taylor and Varma [2020] A. Taylor and V. Varma, Gravitational wave peak luminosity model for precessing binary black holes, Phys. Rev. D 102, 104047 (2020), arXiv:2010.00120 [gr-qc] .
  * Varma  _et al._ [2019a] V. Varma, S. E. Field, M. A. Scheel, J. Blackman, D. Gerosa, L. C. Stein, L. E. Kidder, and H. P. Pfeiffer, Surrogate models for precessing binary black hole simulations with unequal masses, Phys. Rev. Research. 1, 033015 (2019a), arXiv:1905.09300 [gr-qc] .
  * Varma  _et al._ [2019b] V. Varma, D. Gerosa, L. C. Stein, F. Hébert, and H. Zhang, High-accuracy mass, spin, and recoil predictions of generic black-hole merger remnants, Phys. Rev. Lett. 122, 011101 (2019b), arXiv:1809.09125 [gr-qc] .
  * Blanchet  _et al._ [2005] L. Blanchet, M. S. S. Qusailah, and C. M. Will, Gravitational recoil of inspiralling black-hole binaries to second post-Newtonian order, Astrophys. J. 635, 508 (2005), arXiv:astro-ph/0507692 .
  * Racine  _et al._ [2009] E. Racine, A. Buonanno, and L. E. Kidder, Recoil velocity at 2PN order for spinning black hole binaries, Phys. Rev. D 80, 044010 (2009), arXiv:0812.4413 [gr-qc] .
  * Buonanno  _et al._ [2008] A. Buonanno, L. E. Kidder, and L. Lehner, Estimating the final spin of a binary black hole coalescence, Phys. Rev. D 77, 026004 (2008), arXiv:0709.3839 [astro-ph] .
  * Damour and Gopakumar [2006] T. Damour and A. Gopakumar, Gravitational recoil during binary black hole coalescence using the effective one body approach, Phys. Rev. D 73, 124006 (2006), arXiv:gr-qc/0602117 .
  * Healy and Lousto [2017] J. Healy and C. O. Lousto, Remnant of binary black-hole mergers: New simulations and peak luminosity studies, Phys. Rev. D 95, 024037 (2017), arXiv:1610.09713 [gr-qc] .
  * Healy  _et al._ [2014] J. Healy, C. O. Lousto, and Y. Zlochower, Remnant mass, spin, and recoil from spin aligned black-hole binaries, Phys. Rev. D 90, 104004 (2014), arXiv:1406.7295 [gr-qc] .
  * Lousto  _et al._ [2010] C. O. Lousto, M. Campanelli, Y. Zlochower, and H. Nakano, Remnant Masses, Spins and Recoils from the Merger of Generic Black-Hole Binaries, Class. Quant. Grav. 27, 114006 (2010), arXiv:0904.3541 [gr-qc] .
  * Hemberger  _et al._ [2013] D. A. Hemberger, G. Lovelace, T. J. Loredo, L. E. Kidder, M. A. Scheel, B. Szilágyi, N. W. Taylor, and S. A. Teukolsky, Final spin and radiated energy in numerical simulations of binary black holes with equal masses and equal, aligned or anti-aligned spins, Phys. Rev. D 88, 064014 (2013), arXiv:1305.5991 [gr-qc] .
  * Lousto and Zlochower [2014] C. O. Lousto and Y. Zlochower, Black hole binary remnant mass and spin: A new phenomenological formula, Phys. Rev. D 89, 104052 (2014), arXiv:1312.5775 [gr-qc] .
  * Zlochower and Lousto [2015] Y. Zlochower and C. O. Lousto, Modeling the remnant mass, spin, and recoil from unequal-mass, precessing black-hole binaries: The Intermediate Mass Ratio Regime, Phys. Rev. D 92, 024022 (2015), [Erratum: Phys.Rev.D 94, 029901 (2016)], arXiv:1503.07536 [gr-qc] .
  * Lousto and Zlochower [2011] C. O. Lousto and Y. Zlochower, Modeling maximum astrophysical gravitational recoil velocities, Phys. Rev. D 83, 024003 (2011), arXiv:1011.0593 [gr-qc] .
  * Keitel  _et al._ [2017] D. Keitel _et al._ , The most powerful astrophysical events: Gravitational-wave peak luminosity of binary black holes as predicted by numerical relativity, Phys. Rev. D 96, 024006 (2017), arXiv:1612.09566 [gr-qc] .
  * Sopuerta  _et al._ [2007] C. F. Sopuerta, N. Yunes, and P. Laguna, Gravitational recoil velocities from eccentric binary black hole mergers, Astrophys. J. Lett. 656, L9 (2007), arXiv:astro-ph/0611110 .
  * Sperhake  _et al._ [2008] U. Sperhake, E. Berti, V. Cardoso, J. A. Gonzalez, B. Bruegmann, and M. Ansorg, Eccentric binary black-hole mergers: The Transition from inspiral to plunge in general relativity, Phys. Rev. D 78, 064069 (2008), arXiv:0710.3823 [gr-qc] .
  * Hinder  _et al._ [2008] I. Hinder, B. Vaishnav, F. Herrmann, D. Shoemaker, and P. Laguna, Universality and final spin in eccentric binary black hole inspirals, Phys. Rev. D 77, 081502 (2008), arXiv:0710.5167 [gr-qc] .
  * Huerta  _et al._ [2019] E. A. Huerta _et al._ , Physics of eccentric binary black hole mergers: A numerical relativity perspective, Phys. Rev. D 100, 064003 (2019), arXiv:1901.07038 [gr-qc] .
  * Sperhake  _et al._ [2020] U. Sperhake, R. Rosca-Mead, D. Gerosa, and E. Berti, Amplification of superkicks in black-hole binaries through orbital eccentricity, Phys. Rev. D 101, 024044 (2020), arXiv:1910.01598 [gr-qc] .
  * Radia  _et al._ [2021] M. Radia, U. Sperhake, E. Berti, and R. Croft, Anomalies in the gravitational recoil of eccentric black-hole mergers with unequal mass ratios, Phys. Rev. D 103, 104006 (2021), arXiv:2101.11015 [gr-qc] .
  * Ferguson  _et al._ [2019] D. Ferguson, S. Ghonge, J. A. Clark, J. Calderon Bustillo, P. Laguna, D. Shoemaker, and J. Calderon Bustillo, Measuring Spin of the Remnant Black Hole from Maximum Amplitude, Phys. Rev. Lett. 123, 151101 (2019), arXiv:1905.03756 [gr-qc] .
  * Carullo  _et al._ [2023] G. Carullo, S. Albanesi, A. Nagar, R. Gamba, S. Bernuzzi, T. Andrade, and J. Trenado, Unveiling the merger structure of black hole binaries in generic planar orbits, (2023), arXiv:2309.07228 [gr-qc] .
  * Wang  _et al._ [2023] H. Wang, Y.-C. Zou, Q.-W. Wu, Y. Liu, and X. Liu, Characterizing the effect of eccentricity on the dynamics of binary black hole mergers in numerical relativity, (2023), arXiv:2310.04777 [gr-qc] .
  * [68] Caltech-Cornell-CITA, Rit binary black hole simulation results, https://ccrg.rit.edu/numerical-simulations.
  * Zlochower  _et al._ [2005] Y. Zlochower, J. G. Baker, M. Campanelli, and C. O. Lousto, Accurate black hole evolutions by fourth-order numerical relativity, Phys. Rev. D 72, 024021 (2005), arXiv:gr-qc/0505055 .
  * Loffler  _et al._ [2012] F. Loffler _et al._ , The Einstein Toolkit: A Community Computational Infrastructure for Relativistic Astrophysics, Class. Quant. Grav. 29, 115001 (2012), arXiv:1111.3344 [gr-qc] .
  * Schnetter  _et al._ [2004] E. Schnetter, S. H. Hawley, and I. Hawke, Evolutions in 3-D numerical relativity using fixed mesh refinement, Class. Quant. Grav. 21, 1465 (2004), arXiv:gr-qc/0310042 .
  * Nakamura  _et al._ [1987] T. Nakamura, K. Oohara, and Y. Kojima, General Relativistic Collapse to Black Holes and Gravitational Waves from Black Holes, Prog. Theor. Phys. Suppl. 90, 1 (1987).
  * Shibata and Nakamura [1995] M. Shibata and T. Nakamura, Evolution of three-dimensional gravitational waves: Harmonic slicing case, Phys. Rev. D 52, 5428 (1995).
  * Baumgarte and Shapiro [1998] T. W. Baumgarte and S. L. Shapiro, On the numerical integration of Einstein’s field equations, Phys. Rev. D 59, 024007 (1998), arXiv:gr-qc/9810065 .
  * Thornburg [2004] J. Thornburg, A Fast apparent horizon finder for three-dimensional Cartesian grids in numerical relativity, Class. Quant. Grav. 21, 743 (2004), arXiv:gr-qc/0306056 .
  * Campanelli  _et al._ [2007] M. Campanelli, C. O. Lousto, Y. Zlochower, B. Krishnan, and D. Merritt, Spin Flips and Precession in Black-Hole-Binary Mergers, Phys. Rev. D 75, 064030 (2007), arXiv:gr-qc/0612076 .
  * Brandt and Bruegmann [1997] S. Brandt and B. Bruegmann, A Simple construction of initial data for multiple black holes, Phys. Rev. Lett. 78, 3606 (1997), arXiv:gr-qc/9703066 .
  * Ansorg  _et al._ [2004] M. Ansorg, B. Bruegmann, and W. Tichy, A Single-domain spectral method for black hole puncture data, Phys. Rev. D 70, 064011 (2004), arXiv:gr-qc/0404056 .
  * Healy  _et al._ [2017b] J. Healy, C. O. Lousto, H. Nakano, and Y. Zlochower, Post-Newtonian Quasicircular Initial Orbits for Numerical Relativity, Class. Quant. Grav. 34, 145011 (2017b), arXiv:1702.00872 [gr-qc] .

##  V Supplementary material

Supplemental materials

##  VI Oscillations in peak amplitudes

In this section, in Fig. S1 we show the variations of higher order harmonics
peak amplitudes, $\mathcal{A}_{32,\text{peak}}$ (panel (a)),
$\mathcal{A}_{44,\text{peak}}$ (panel (b)) as a function of the initial
eccentricity $e_{0}$ at the initial coordinate separation of
$D_{\text{ini}}=11.3M$ and $D_{\text{ini}}=24.6M$ for nonspinning
configuration with different mass ratio $q=1,0.75,0.5,0.25$. We can see that
there is a universal oscillation in higher order harmonics peak amplitudes
similar to $L_{\text{peak}}$, masses $M_{f}$, spins $\alpha_{f}$, recoil
velocity $V_{f}$, and $\mathcal{A}_{22,\text{peak}}$ [67].

Figure S1: Variations of higher order harmonics peak amplitudes
$\mathcal{A}_{32,\text{peak}}$ (panel (a)), $\mathcal{A}_{44,\text{peak}}$
(panel (b)) as a function of the initial eccentricity $e_{0}$ at the initial
coordinate separation of $D_{\text{ini}}=11.3M$ and $D_{\text{ini}}=24.6M$ for
nonspinning configuration with different mass ratio $q=1,0.75,0.5,0.25$.

##  VII Integer orbital cycles in peak amplitudes

In this section, we follow the ideas of Ref. [67] to consider the orbital
transition. The number of orbital cycles $N$ can be determined through the
phase of the gravitational waveform. In our analysis, we specifically focus on
the 2-2 mode. To calculate the phase difference, we evaluate the expression:

$\Delta\Phi=\Phi\left(t_{\mathrm{merger}}\right)-\Phi\left(t_{0}+t_{\mathrm{relax}}\right).$
(S1)

where $t_{\mathrm{merger}}$ represents the time of BBH merger, $t_{0}$ denotes
the initial moment of the waveform, and $t_{\mathrm{relax}}$ signifies the
time required to the transition from the initial moment to a physically stable
state. The number of orbital cycles accomplished by the BBH system can be
obtained as:

$N_{\text{waves }}=\frac{\Delta\Phi}{4\pi}.$ (S2)

Here, we divide the phase difference $\Delta\Phi$ by $4\pi$ since the waveform
phase of 2-2 mode is twice that of the orbital phase. Fig. S2 displays the
relationship between the integer orbital cycle number $N_{\text{waves}}$ and
peak amplitudes $\mathcal{A}_{22,\text{peak}}$ (panel (a)),
$\mathcal{A}_{32,\text{peak}}$ (panel (b)), $\mathcal{A}_{44,\text{peak}}$
(panel (c)) at initial coordinate separations $D_{\text{ini}}=11.3M$ and
$D_{\text{ini}}=24.6M$ for nonspinning configuration with different mass ratio
$q=1,0.75,0.5,0.25$. These points, denoted by red “x” markers, correspond to
either an integer multiple or are in close proximity to an integer multiple of
the orbital cycles. Moving from right to left, each red “x” corresponds to
successive orbital cycles, starting from cycle 1 and continuing indefinitely.

Figure S2: Relationship between the integer orbital cycle number
$N_{\text{waves}}$ and peak amplitudes $\mathcal{A}_{22,\text{peak}}$ (panel
(a)), $\mathcal{A}_{32,\text{peak}}$ (panel (b)),
$\mathcal{A}_{44,\text{peak}}$ (panel (c)) at initial coordinate separations
$D_{\text{ini}}=11.3M$ and $D_{\text{ini}}=24.6M$ for nonspinning
configuration with different mass ratio $q=1,0.75,0.5,0.25$. These points,
denoted by red “x” markers, correspond to either an integer multiple or are in
close proximity to an integer multiple of the orbital cycles. Moving from
right to left, each red “x” corresponds to successive orbital cycles, starting
from cycle 1 and continuing indefinitely.

##  VIII Other correlations

In this section, we show the other 18 correlations between dynamic quantities
$L_{\text{peak}}$, masses $M_{f}$, spins $\alpha_{f}$, and recoil velocity
$V_{f}$ and peak amplitudes $\mathcal{A}_{22,\text{peak}}$ ,
$\mathcal{A}_{32,\text{peak}}$, $\mathcal{A}_{44,\text{peak}}$ in Fig. S3,
Fig. S4 and Fig. S5.

Figure S3: Correlations between dynamic quantities $L_{\text{peak}}$, masses
$M_{f}$, spins $\alpha_{f}$, and recoil velocity $V_{f}$ and peak amplitudes
$\mathcal{A}_{22,\text{peak}}$ , $\mathcal{A}_{32,\text{peak}}$,
$\mathcal{A}_{44,\text{peak}}$. Figure S4: Correlations between dynamic
quantities $L_{\text{peak}}$, masses $M_{f}$, spins $\alpha_{f}$, and recoil
velocity $V_{f}$ and peak amplitudes $\mathcal{A}_{22,\text{peak}}$ ,
$\mathcal{A}_{32,\text{peak}}$, $\mathcal{A}_{44,\text{peak}}$. Figure S5:
Correlations between dynamic quantities $L_{\text{peak}}$, masses $M_{f}$,
spins $\alpha_{f}$, and recoil velocity $V_{f}$ and peak amplitudes
$\mathcal{A}_{22,\text{peak}}$ , $\mathcal{A}_{32,\text{peak}}$,
$\mathcal{A}_{44,\text{peak}}$.

##  IX The origin of fingerprints

In this section, we select two typical cases to analyze the relationship
between oscillation and this spiral structure. Other situations are the same
as it. We mark the positions of the peaks and valleys of oscillations of
dynamic quantities $L_{\text{peak}}$, masses $M_{f}$, spins $\alpha_{f}$ in
the spiral structure in Fig. S6. Notably, we find that the peaks and valleys
of oscillation are positioned on opposite sides of the spiral and are
concentrated near a straight line. For $M_{f}$ and $L_{\text{peak}}$, there
appears to be an approximate phase difference (see panel (a), (c), (e), (g) of
Fig. S6), while $L_{\text{peak}}$ and $\alpha_{f}$ nearly coincide (see panel
(b), (d), (f), (h) of Fig. S6).

Figure S6:  Points corresponding to the peaks and valleys of oscillations in
the relationship between dynamic quantities $L_{\text{peak}}$, masses $M_{f}$,
spins $\alpha_{f}$ and $e_{0}$ marked with red “x”.


[datatype=bibtex] [overwrite]
[fieldsource=shortjournal,fieldtarget=journaltitle] [bib]nametitledelim,

# Analysis of Iterative Deblurring:
No Explicit Noise

###### Abstract.

Iterative deblurring, notably the Richardson-Lucy algorithm with and without
regularization, is analyzed in the context of nuclear and high-energy physics
applications. In these applications, probability distributions may be
discretized into a few bins, measurement statistics can be high, and
instrument performance can be well understood. In such circumstances, it is
essential to understand the deblurring first without any explicit noise
considerations. We employ singular value decomposition for the blurring matrix
in a low-count pixel system. A strong blurring may yield a null space for the
blurring matrix. Yet, a nonnegativity constraint for images built into the
deblurring may help restore null-space content in a high-contrast image with
zero or low intensity for a sufficient number of pixels. For low-contrast
images, the control over null-space content may be gained through
regularization. When the regularization is applied, the blurred image is, in
practice, restored to an image that is still blurred but less than the
starting one.

###### Key words and phrases:

Richardson-Lucy algorithm, singular value decomposition, null space, image
contrast, regularization

###### 1991 Mathematics Subject Classification:

Primary: 65R30, 65R32; Secondary: 65Z05.

∗Corresponding author

Paweł Danielewicz was supported by the U.S. Department of Energy Office of
Science Grant No. DE-SC0019209.

Sinethemba Neliswa Mamba${}^{{\href<EMAIL_ADDRESS>and Paweł
Danielewicz ${}^{{\href danielew@msu.edu}*1,2}$

1African Institute for Mathematical Sciences, Kigali, Rwanda

2Facility for Rare Isotope Beams and Department of Physics and Astronomy,

Michigan State University, East Lansing, Michigan, USA

(Communicated by Handling Editor)

##  1\. Introduction

Blurring is common in image acquisition, making the images less sharp and
clear. When the blurring process is understood, it may be possible to remove,
fully or partially, the blurring effects or to deblur the images. Here, we
will analyze what a standard iterative deblurring methodology can and cannot
do, with an eye on deblurring to improve quantitative results of particle
yield measurements in nuclear and high-energy physics. In the latter areas,
the deblurring gained attention relatively recently [2, 13, 3, 9, 14]. The
measurement statistics can be higher than in many optical applications, making
no-nominal-noise analysis relevant. At the same time, the number of bins in
these areas, equivalent to pixels, can be modest.

A convenient framework in which the deblurring can be understood is that of
the Singular Value Decomposition (SVD) for the blurring matrix, already used
in the past in the literature [5, 10, 15]. With some advantages that the
Richardson-Lucy (RL) [11, 8] and related iterative methods [6, 4, 3] have over
other iterative methods for the restoration of brightness [1, 12, 7], we will
concentrate here on the RL method. One novel aspect of our work will examine
the contrast’s impact on the objects imaged in the deblurring process.
Further, we will actively seek the cases for which the deblurring fails, try
to understand the failures, and predict when these may occur. Small or
vanishing eigenvalues from SVD will play a role there. In the present work, we
will concentrate on the deblurring when no explicit noise appears in the
problem, which already has some complexity. In future work, we shall consider
practicalities around the noise in the problem, for which the present work
will provide a base.

This paper is organized as follows. Section 2 discusses general blurring
concepts and formulates a schematic model for explorations. We decompose the
blurring matrix using SVD. Section 3 describes selected deblurring methods,
including the RL and Landweber (LW) algorithms, and illustrates their
practical operation. Section 4 explores the limits of the RL deblurring with
and without regularization. We conclude in Section 5.

##  2\. Blurring Matrices and Their Singular Value Decomposition

The blurring relation can be stated in the form of the equation [8]

$\displaystyle\mathbf{g}=T\mathbf{F}\,.$ (1)

Here, $\mathbf{g}\in\mathbb{R}^{n}$ is the vector of values for the blurred
image, $T\in\mathbb{R}^{n\times n}$ is the blurring matrix or transfer matrix
and $\mathbf{F}\in\mathbb{R}^{n}$ is the vector for the blur-free image.
Probability conservation for the blurring matrix implies that the sum of all
the elements of each column (index $i$) is precisely equal to 1:

$\sum_{j=0}^{n-1}T_{ji}=1\,.$

Without practical loss of generality, we shall consider that the indices
pertain to locations in one dimension, representing pixels or other
distribution bins. We shall adopt cyclic boundary conditions within that
space, where the vector and matrix indices are computed modulo $n$. To fix the
attention, we shall consider a couple of exemplary blurring (or point-
spreading) functions, three- and five-bin, out of which we shall consider the
five-bin one most often:

Three-bin

$\displaystyle T_{ji}=\begin{cases}\begin{array}[]{cl}0.550\,,&\text{for
}j=i\,,\\\ 0.225\,,&\text{for }j=(i\pm 1)\mod n\,,\\\
0\,,&\text{otherwise.}\\\ \end{array}\end{cases}$ (2)

Five-bin

$\displaystyle T_{ji}=\begin{cases}\begin{array}[]{cl}0.250\,,&\text{for
}j=i\,,\\\ 0.250\,,&\text{for }j=(i\pm 1)\mod n\,,\\\ 0.125\,,&\text{for
}j=(i\pm 2)\mod n\,,\\\ 0\,,&\text{otherwise.}\\\ \end{array}\end{cases}$ (3)

In the grayscale convention, the pixel brightness values span the interval
[0,1], with 1 representing white and 0 representing black. In deblurring, it
is important not to end up with values outside of that interval, as lacking
physical sense. We illustrate the three- and five-bin blurring functions in
Fig. 1(a) and Fig. 2(a), respectively, where we show the action of these
functions on an original image where only pixel 5 has a value of 1, and other
have 0. The overall size of the space is $n=12$.

(a)

(b)

Figure 1. Illustration of the three-bin blurring function. (a) Image from
blurring an original with pixel 5 at 1 and other at 0, for $n=12$. (b)
Singular values of the blurring function, ordered by their magnitude. Lines
join the values to guide the eye.

(a)

(b)

Figure 2. Illustration of the five-bin blurring function. (a) Image from
blurring an original with pixel 5 at 1 and other at 0, for $n=12$. (b)
Singular values of the blurring function, ordered by their magnitude. Lines
join the values to guide the eye.

Upon applying the transposed blurring matrix to both sides of Eq. (1), we find
the norm equation

$T^{\top}\mathbf{g}=T^{\top}T\,\mathbf{F}\equiv M\mathbf{F}\,.$ (4)

If we apply SVD to $T$, we get

$T=\sum_{\alpha=0}^{n-1}\sigma_{\alpha}\,\mathbf{u}_{\alpha}\,\mathbf{v}_{\alpha}^{\top}\,,$
(5)

where $\mathbf{u}$ and $\mathbf{v}$ are, respectively, the separately
orthonormal left and right singular vectors, and we adopt the convention of
the singular values $\sigma$ to be nonnegative and ordered in nonincreasing
order. Then for the positive definite matrix $M$ we arrive at the following
spectral decomposition

$M=T^{\top}\,T=\sum_{\alpha,\beta}\sigma_{\alpha}\,\sigma_{\beta}\,\mathbf{v}_{\alpha}\,\mathbf{u}_{\alpha}^{\top}\,\mathbf{u}_{\beta}\,\mathbf{v}_{\beta}^{\top}=\sum_{\alpha}\sigma_{\alpha}^{2}\,\mathbf{v}_{\alpha}\,\mathbf{v}_{\alpha}^{\top}\,.$
(6)

Generally, for original and deblurred images, we shall employ a decomposition
in terms of the right singular vectors of the blurring matrix and for the
blurred images - in terms of the left:

$\displaystyle\begin{cases}\mathbf{F}=\sum_{\alpha}F_{\alpha}\,\mathbf{v}_{\alpha}\,,\\\
\mathbf{g}=\sum_{\alpha}g_{\alpha}\,\mathbf{u}_{\alpha}\,.\end{cases}$ (7)

With (1), we then find

$\mathbf{g}=T\mathbf{F}=\sum_{\alpha}\sigma_{\alpha}\,F_{\alpha}\,\mathbf{u}_{\alpha}\,,$
(8)

or $g_{\alpha}=\sigma_{\alpha}\,F_{\alpha}$.

The three- and five-bin blurring matrices, Eqs. (2) and (3), which we use here
as examples, have some peculiar features. These features may often be
encountered in practice, possibly as approximate rather than exact. Thus,
these blurring matrices are symmetric, meaning their singular vectors are
matrix eigenvectors. Further, with the matrix elements of $T$ depending only
on the absolute difference of indices, that dependence extends to the
corresponding $M$ matrices. This implies two types of invariance for the
matrices $T$ and $M$, impacting the singular vectors. One is the invariance
under a shift of matrix indices by some integer. With this, a singular vector,
corresponding to some singular value, remains a singular vector, corresponding
to the same singular value, after a shift in its indices by an integer. The
second invariance is under a reflection of the indices around any chosen
index. With this, a singular vector, corresponding to some singular value,
remains a singular vector, corresponding to the same singular value, after a
reflection of its indices around any chosen index. With these two symmetries,
any singular value is at least doubly degenerate unless the corresponding
singular vector transforms onto itself, up to a factor of -1, under these
transformations.

The spectra of singular values for the two exemplary blurring matrices are
shown in Figs. 1(b) and 2(b), and the singular vectors for the five-bin
function are shown in Fig. 3. Double degeneracy is seen for most of the
singular values. For Fig. 3, we choose the corresponding pairs of singular
vectors in the degenerate subspaces so that some symmetry behind the
degeneracies is evident, either a reflection or translational one, or both. In
each of the spectra, in Figs. 1(b) and 2(b), a nondegenerate singular value of
1 is seen. This value corresponds to the vector that is constant across the
pixels. Under blurring, this vector transforms into itself. Also, this vector
is invariant under the symmetries. The spectrum of singular values falls off
faster for the five- than for the three-bin function, demonstrating a more
significant deterioration of the information under the blurring with a broader
function. In Fig. 3, it is evident that a faster variation of the vectors over
pixels is generally associated with lower singular values. The singular vector
with the most rapid variation, alternating between the same positive and
negative values from one pixel to another, is invariant under the symmetries.
It is the other one that may be nondegenerate. Also, under blurring, it
transforms into itself, up to a factor, because it is uniquely alternating and
reflection-symmetric about every pixel, which will be preserved by any
blurring that depends only on the difference of indices. With a broadening of
the blurring function, this vector will likely be the first to be nullified.
In Fig. 1(b), the singular value for that vector, with index 11, approaches
zero, and the vector is nullified when the blurring matrix evolves into
$T_{jj}=0.5$ and $T_{ji}=0.25$, for $j=(i\pm 1)\mod 12$. In Fig. 2(b), the
maximally alternating vector is accompanied in the null space by two more
vectors, cf., panels (a) and (f) in Fig. 3.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3. Singular vectors for the five-bin blurring function. Within
degenerate spaces, pairs of vectors are chosen to transform onto the other or
itself under reflection and/or translation. Panel (a) shows the two unique
vectors that transform onto themselves under the transformations.

##  3\. Deblurring Methods

We next turn to deblurring. We illustrate the deblurring in Fig. 4(a). We take
a $\delta$ distribution where only one bin gets intensity 1, and the others
get 0. We blur the original image with the 5-bin function (3) and deblur the
outcome with the RL method [11, 8], of which the details will be discussed.
The blurred and restored images are shown in Fig. 4(a) together with the
original. The eye cannot distinguish the RL-restored and original images, thus
represented with one set of graphics. In addition, intensity values for the
images are shown in Table 1, and within a precision nearing that of the
machine we employ, they cannot be distinguished. This may be perceived as
astounding, as the singular vectors with the fastest bin-to-bin variation span
the null space for the blurring matrix, and they most certainly contribute to
the SVD of the $\delta$ distribution.

(a)

(b)

Figure 4. Deblurring test with a single pixel with index 5 at the intensity 1 and all other pixels at the intensity of 0. (a) The original image $F$ (circles), image $g$ (triangles) from blurring $F$ with the 5-pixel blurring function, and the restored images $f$ from different deblurring methods are displayed. As the image restored following the RL method is not distinguishable by the eye from the original, the circles represent that image, too. The image restored using the LW method is represented by diamonds. Finally, the image that was restored with SVD, where the operational space was reduced ($m=6$ in Eqs. (21) and (22)) compared to LW, is represented by circles. The symbols corresponding to each case are connected with lines. (b) Coefficients of decomposition in the basis of singular vectors vs singular value index for the original and blurred images and images restored with the RL and LW methods. $i$ | $F_{i}$ | $g_{i}$ | $f_{i}^{RL}$ | $f_{i}^{LW}$ | $f_{i}^{SVD}$
---|---|---|---|---|---
0 | 0 | 0.000 | 0.000 | 0.083 | -0.061
1 | 0 | 0.000 | 0.000 | -0.250 | -0.167
2 | 0 | 0.000 | 0.000 | 0.083 | 0.083
3 | 0 | 0.125 | 0.000 | 0.083 | 0.000
4 | 0 | 0.250 | 0.000 | 0.083 | 0.228
5 | 1 | 0.250 | 1.000 | 0.750 | 0.583
6 | 0 | 0.250 | 0.000 | 0.083 | 0.228
7 | 0 | 0.125 | 0.000 | 0.083 | 0.000
8 | 0 | 0.000 | 0.000 | 0.083 | 0.083
9 | 0 | 0.000 | 0.000 | -0.250 | -0.167
10 | 0 | 0.000 | 0.000 | 0.083 | -0.061
11 | 0 | 0.000 | 0.000 | 0.083 | 0.250
Table 1. Bin intensities in deblurring the $\delta$-distribution centered at
index 5, illustrated in Fig. 4(a). From left to right, the columns represent
the bin index ($i$), and then the intensity for the original image ($F_{i}$),
for the 5-bin blurred image ($g_{i}$), and the image restored with three
different methods: RL ($f_{i}^{RL}$), LW ($f_{i}^{LW}$) and SVD
($f_{i}^{SVD}$).

In the following, we review exemplary deblurring methods within the context of
their potential utility in quantitative research analyses. After this section,
we shall examine under what circumstances the deblurring can yield faithful
results, such as in the example above, and we shall seek conditions and manner
under which the deblurring can fail. The iterative methods aim to restore a
blurred image to its original form by iteratively improving the estimate of
the blur-free image. These methods take as input some initial guess of the
blur-free image, the blurring matrix, and a criterion for the number of
iterations to be carried out. The typical initial guesses are the uniform
image or the blurred image. As the blur-free image is to improve with each
iteration, naively, a more accurate output is to be obtained with more
iterations - in practice, this may not be the case.

We use below $k$ as an iteration index while discussing the iterative methods.

###  3.1. Richardson-Lucy method, classical and regularized

The classical RL algorithm iteratively upgrades the restored image $f$, given
the blurred image $g$, according to [11, 8, 3]

$f_{i}^{(k+1)}=A_{i}^{(k)}\,f_{i}^{(k)}\,.$ (9)

Here, the amplification factor $A_{i}$ for the $i$’th pixel intensity is

$A_{i}^{(k)}=\frac{1}{q_{i}}\sum_{j}T_{ji}\frac{g_{j}}{g_{j}^{(k)}}\,,$ (10)

$i,j,m=0,1,2,\ldots,n-1$, $g_{j}^{(k)}=\sum_{m}T_{jm}f_{m}^{(k)}$ is the
blurred-image prediction for the restored image $f^{(k)}$ at the $k$’th step
of restoration, $T$ is the before introduced point spread function or blurring
matrix, and $q_{i}=\sum_{j}T_{ji}$ serves the normalization. The normalization
factors for the introduced blurring-function examples are unity, $q_{i}\equiv
1$.

The regularized RL algorithm modifies the iteration equation to

$\displaystyle f_{i}^{(k+1)}=A_{i}^{(k)}\,I^{k}_{i}\,f_{i}^{(k)}\quad
i,j,m=0,1,2,\ldots,n-1\,,$ (11)

where the regularization factor is [4, 3]

$\begin{split}I^{(k)}_{i}&=\frac{1}{1-\Lambda\operatorname{div}\left(\frac{\nabla
f_{i}^{(k)}}{|\nabla f_{i}^{(k)}|}\right)}\\\
&=\begin{cases}\frac{1}{1-\Lambda}\,,&\text{ if
}f^{(k)}_{i}<f^{(k)}_{[(i-1)\,\mathrm{mod}\,n],\,[(i+1)\,\mathrm{mod}\,n]}\,,\\\
\frac{1}{1+\Lambda}\,,&\text{ if
}f^{(k)}_{i}>f^{(k)}_{[(i-1)\,\mathrm{mod}\,n],\,[(i+1)\mathrm{mod}\,n]}\,,\\\
1\,,&\text{ }\text{ otherwise }\,.\end{cases}\end{split}$ (12)

for $i=0,1,2,\ldots,n-1$. Here, the r.h.s. is the one-dimensional version [3]
of the center expression [4] in which $\nabla$ is the gradient and $|\nabla
f_{i}^{(k)}|$ is $\ell_{1}$-norm; $\Lambda$ is a small positive regularization
coefficient. The role of the regularization factor is the suppression of
uncontrolled rise, during iterations, of contributions to the restored image
from the fastest varying singular vectors corresponding to vanishing or
minimal singular values. The typical signature of such an uncontrolled rise is
developing a bin-to-bin sea-saw pattern in the restored image. As SVD might
not be applied in parallel to the RL-restoration and the direct connection of
the regularization factor $I$ to SVD is not obvious, it is common to make
$\Lambda$ small enough that no pathological signs in the restored image
appear.

The restoration of the $\delta$-distribution of Fig. 4(a) with the RL method,
from the blurred image there, with or without regularization, yields a result
that is not distinguishable by eye from the original, so they share the
graphical representation in the figure. The numerical values for the
restoration are provided in Table 1. Another test to which we subject the
deblurring methods is presented in Fig. 5. The original image $F$ is a ramp
with pixels 5-8 at intensity 1 and other pixels at intensity 0. The blurring,
yielding $g$, is performed with the 5-bin blurring function. The results of RL
deblurring, without or with regularization, are not distinguishable by eye
from the original, but so are the results of deblurring with other methods
that we will discuss next. The numerical values from the ramp tests can be
found in Table 2.

Figure 5. Deblurring test with a 4-pixel ramp. In the original image, the intensity is 1 for the pixels 5-8 and 0 for the others. The operational space for SVD is here the full row space of $T$. $i$ | $F_{i}$ | $g_{i}$ | $f_{i}^{RL}$ | $f_{i}^{LW}$ | $f_{i}^{SVD}$
---|---|---|---|---|---
0 | 0 | 0.000 | 0.000 | 0.000 | 0.000
1 | 0 | 0.000 | 0.000 | 0.000 | 0.000
2 | 0 | 0.000 | 0.000 | 0.000 | 0.000
3 | 0 | 0.125 | 0.000 | 0.000 | 0.000
4 | 0 | 0.375 | 0.000 | 0.000 | 0.000
5 | 1 | 0.625 | 1.000 | 1.000 | 1.000
6 | 1 | 0.875 | 1.000 | 1.000 | 1.000
7 | 1 | 0.875 | 1.000 | 1.000 | 1.000
8 | 1 | 0.625 | 1.000 | 1.000 | 1.000
9 | 0 | 0.375 | 0.000 | 0.000 | 0.000
10 | 0 | 0.125 | 0.000 | 0.000 | 0.000
11 | 0 | 0.000 | 0.000 | 0.000 | 0.000
Table 2. The numerical values from the 4-pixel ramp deblurring test, with the
RL, SVD, and LW methods, illustrated in Fig. 5.

###  3.2. Landweber Method

The Landweber (LW) method [7] can be viewed as the steepest-descent method for
minimizing the square deviation between the brightness distribution for the
generated and blurred images,

$\left\|\mathbf{g}-T\mathbf{f}\right\|^{2}_{2}\,.$ (13)

Iterative adjustments in $\mathbf{f}$ take the form

$\mathbf{f}^{(k+1)}=\mathbf{f}^{(k)}+\beta
T^{T}(\mathbf{g}-T\mathbf{f}^{(k)})\,,$ (14)

where $0<\beta<2/\|T^{T}T\|$ for convergence. With the $T$ multiplication in
(14), it is apparent that the iterations do not change the null-space content
in $f$. The null-space content may be suppressed completely in $f$ by
replacing the minimized quantity (13) with

$\left\|\mathbf{g}-T\mathbf{f}\right\|^{2}_{2}+\epsilon\left\|\mathbf{f}\right\|^{2}_{2}\,,$
(15)

where $0<\epsilon\ll(\sigma_{\alpha})_{min}^{2}$ and $(\sigma_{\alpha})_{min}$
is the smallest singular value of which the contribution to $f$ should not be
disturbed. The iterations are then modified to

$\mathbf{f}^{(k+1)}=(1-\beta\epsilon)\mathbf{f}^{(k)}+\beta
T^{T}(\mathbf{g}-T\mathbf{f}^{(k)})\,.$ (16)

The results from the deblurring of the blurred $\delta$-distribution with the
LW method are provided in Fig. 4 and Table 1. The corresponding results from
processing the 4-pixel ramp are provided in Fig. 5 and Table 2. While the LW
method struggles to restore the $\delta$ distribution, yielding, in
particular, negative intensity values, it faithfully restores the 4-pixel
ramp. This may be surprising, and we will return to this issue later in this
Section.

###  3.3. Deblurring Using Singular Value Decomposition

The deblurring can be further carried out by employing SVD explicitly. In
general, the SVD deblurring [5, 15, 10] seeks the best approximation to the
blurred image in a reduced subspace of singular vectors, which in particular
should exclude the null space for $T$, minimizing the square deviation

$\left\|\mathbf{g}-T_{m}\mathbf{f}\right\|^{2}_{2}\,,$ (17)

where

$T_{m}=\sum_{\alpha=0}^{m-1}\sigma_{\alpha}\,\mathbf{u}_{\alpha}\,\mathbf{v}_{\alpha}^{\top}\,,$
(18)

and $m\leq n$. The minimization of (17) is equivalent to the minimization of
the square deviation within the reduced subspace only

$\left\|\mathbf{g}_{m}-T_{m}\mathbf{f}\right\|^{2}_{2}\,,$ (19)

where

$\mathbf{g}_{m}=\sum_{\alpha=0}^{m-1}g_{\alpha}\,\mathbf{u}_{\alpha}\,,$ (20)

with the minimum reached for

$f_{\alpha}=\frac{g_{\alpha}}{\sigma_{\alpha}}\,,\hskip
20.00003pt\text{for}\hskip 5.0pt\alpha<m\,.$ (21)

Outside of the considered subspace, the expansion coefficients for
$\mathbf{f}$ are set to zero,

$f_{\alpha}=0\,,\hskip 20.00003pt\text{for}\hskip 5.0ptm\leq\alpha<n\,.$ (22)

In the tests of this Section, the SVD method yields the same results as LW
when its operational space is maximal, i.e., identical to the row space of the
blurring matrix $T$. Regarding the 4-pixel ramp test with results in Fig. 5
and Table 2, we keep the operational space for SVD maximal. However, in the
$\delta$-function test with results in Fig. 4 and Table 1, we reduce the
operational space for SVD to $m=6$, see Fig. 3 for the corresponding singular
vectors. With the reduction, the quality of the restoration significantly
deteriorates. The location and strength of negative intensities in the
restored image are similar for $m=8$ (identical to LW within precision) and
$m=6$, cf. Eqs. (21) and (22).

###  3.4. Methods’ Assessment

We call the original images, such as in Figs. 4(a) and 5, where a substantial
fraction of the pixels has an intensity zero or low compared to the maximal,
as having high contrast. When such images are moderately blurred, the
deblurring methods with the nonnegativity constraint for intensity built into
the deblurring process, such as RL, generally perform far better than methods
that lack the constraint, such as LW or SVD, see the case of the
$\delta$-distribution of Fig. 4(a). In this case, the two latter methods yield
negative intensity values in the restored images. Amazingly, while the
$\delta$-distribution contains contributions from the $v_{10}$ and $v_{11}$
vectors in the null space, the nonnegativity constraint allows the RL method
to restore their contributions faithfully despite their singular values being
zero, see Fig. 4(b). Without the nonnegativity constraint, the content of the
null space, in the case of LW, or excluded space, in the case of SVD, stays at
zero.

One puzzling outcome of the tests so far is that the LW method and the SVD
method, when its excluded space is identical to the null space, perform as
well as RL in restoring the 4-pixel ramp in Fig. 5, even though they struggle
with the $\delta$-distribution. The mystery is solved when looking at the
null-space singular vectors in Fig. 3 and the original ramp in Fig. 5: a ramp
with an even number of pixels at the same intensity lacks any null-space
contribution for our blurring function of Eq. (3). When the ramp’s extension
is changed to an odd number of pixels or the ramp is sloped, the LW and SVD
methods can develop negative intensity values for the restored image, just
like in the case of the $\delta$ distribution.

(a)

(b)

Figure 6. Restored-image expansion coefficients in singular vectors, as a
function of the RL iteration step $k$, normalized to the coefficients for the
original $\delta$-distribution image, for even-index (a) and odd-index (b)
singular vectors.

Given the benefits of the built-in nonnegativity constraint, we will further
exclusively rely on the RL deblurring method. Other methods may perform
similarly when supplemented with this type of constraint. In employing the RL
method, the iteration progress is sometimes used as a tool. Figure 6 shows how
the expansion coefficients in singular vectors for the restored image behave
as a function of the RL iteration step. Here, the original image is the
$\delta$-distribution, and the starting image is the uniform distribution. One
can observe that the higher the index of the singular vector, the generally
slower the approach of the coefficient to its asymptotic value. This can be
understood in terms of the weakening impact of expansion coefficients on the
amplification factors (10) when the singular value decreases, as

$\frac{\text{d}A_{i}^{(k)}}{\text{d}f_{\alpha}^{(k)}}\propto\sigma_{\alpha}\,,$
(23)

for $g_{i}>0$. Null-space and low singular-value coefficients will be impacted
by the positivity constraint and/or regularization, and in Fig. 6, we can see
their pace of approach to the asymptotic values to be lumped.

##  4\. Exploring the limits of deblurring

In the context of the deblurring applications in quantitative research in
nuclear and particle physics, it is important to understand when and how the
deblurring can falsify the results. We have observed the role that positive
definiteness could play in restoring the null-space contributions to a
distribution. However, intuitively, one could expect the positive definiteness
to play a little role when the original distribution is uniformly far from
zero on the general distribution’s variation scale. Below, we will compare the
performance of the classic RL method in those different situations. In cases
when there is no positive definiteness that one can fall on, the
regularization can help improve the performance of a deblurring method.
However, the regularization must be a compromise with some adverse effects,
too, as may be evident from Eq. (15).

###  4.1. High vs Low Contrast

For high-contrast original images, such as in Fig. 4(a) and 5, we have
observed that the deblurring with the nonnegativity constraint can restore the
null-space components from the blurred images. Notably, for the nonnegativity
constraint to be impactful, the number of pixels at low intensity compared to
the rest must be comparable to the dimension of the null space. We now
confront that situation to one for low-contrast original images where the
intensity for most pixels varies within a small range relative to the maximal
intensity.

An exemplary low-contrast case of a $\delta$-function superimposed on a high
uniform background is illustrated in Fig. 7(a). When the classical RL method
is applied to the blurred image, the method fails, on the scale of the
$\delta$-function norm, in a similar manner as the methods without the
nonnegativity constraint in the high-contrast case of Fig. 4(a). Again, we
start the RL iterations in this illustration with a uniform image. Though the
nonnegativity is there for the RL method, it does not impact restoration when
the background is sufficiently elevated. Fig. 7(b) complements the results of
(a) with coefficients of decomposition in the basis of singular vectors,
demonstrating that the coefficients for the null space are not restored for
the classical RL method, unlike in Fig. 4(b). For impactful constraints, we
turn to the regularized RL method, where restoration is modified compared to
the classical method depending on the differences in the intensity of adjacent
pixels, suppressing local extrema. With a low $\Lambda=0.0001$, the original
image is satisfactorily restored in Fig. 7(a).

(a)

(b)

Figure 7. Deblurring test with a single pixel with index 5 at intensity
elevated by 0.25 over a high uniform background at 0.75. (a) Crosses represent
the original image $F$ and triangles - the image $g$ from blurring $F$ with
the 5-pixel blurring function. Moreover, the restored images $f$ from
different versions of the RL method are displayed. Circles represent the image
restored using the classical RL method. Diamonds, wide and slim, and stars
represent the images restored with the regularized RL method at, respectively,
$\Lambda=0.001$, 0.001, and 0.01. The symbols corresponding to each case are
connected with lines. (b) Coefficients of decomposition in the basis of
singular vectors for the cases in (a) vs singular-value index.

Additional examples of restoration of images at different contrast levels are
shown in Fig. 8.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 8. Deblurring tests such as in Fig. 7 for other exemplary low-contrast
images, with matching information for each of these presented in panels (a)
and (b), (c) and (d), and (e) and (f), respectively.

###  4.2. Impact of Initialization

For a low-contrast original image $F$, it is possible to add a significant
null-space component to that image, without violating intensity nonnegativity.
When such a combination is fed to the classical RL method as the starting
guess $f^{(0)}$, mathematically, the restoration iterations will never
suppress the added spurious null-space component, yielding that starting guess
as the fully restored image $f$. Even added components out of low-singular-
value vectors may be only slowly modified in iterations.

To get around excessive spurious null-space or low-singular value admixtures
to the restored image, it helps to start the classical RL iterations with an
image that lacks such components, such as the uniform image. In Fig. 9, we
show what happens within the classical RL restoration when the restoration of
a low-contrast image is started with a significant null-space component in the
first guess for the image instead of the uniform guess. The null-space
component’s strength persists for $\Lambda=0$ at significant strength in the
limit of many restoration steps. Notably, the null-space components are not
likely to adjust to the values for the original image in the absence of an
impact of the nonnegativity constraint, no matter how the iterations are
started, including with the uniform image. In Figs. 7 and 8, we can see that
classical restorations of low-contrast images started with uniform guesses,
tend to underestimate the magnitude of the null-space components in the
restored images. In particular, the panels (a) and (b) in Fig. 8 are analogous
to Fig. 9 but differ importantly in initialization. As may be seen in Figs. 7,
8, and 9, the null-space components get under some control after the
regularization is switched on.

(a)

(b)

Figure 9. Sensitivity of an RL restoration to null-space content at iteration
start and to the regularization strength. The RL iterations here are started
with the original image with an $\alpha=11$ component added. The classical RL
method ($\Lambda=0$) leaves the null-space content intact. An RL method with
regularization ($\Lambda>0$) suppresses the superimposed spurious null-space
content, but excessive regularization strength, such as for $\Lambda=0.01$,
can chip off the null-space content in the restored image compared to the
original.

###  4.3. Impact of Regularization

In examining the cases of restoration of interesting original images with few
significant intensity jumps in Figs. 7 and 8, we can see that the
regularization can dramatically improve image restoration. The particular
regularization of Eqs. (11) and (12) biases against multiple pronounced
extrema and, with that, suppresses significant spurious null-space
contributions that could be added to the restored low-contrast image for the
classical RL method.

While a modest regularization within the RL method can help restore the null-
space content of the original image, an excessive regularization can act to
suppress not only that content but also contributions from singular vectors
corresponding to low nonzero singular values. This is illustrated in Figs. 7
and 8, which include results from RL restorations with progressively
strengthening regularization. We can see that with the progressing strength of
regularization, the features of the restored images begin to evolve away from
the original image toward the blurred image. In a sense, the regularization
needed for low-contrast images amounts to trading stronger blurring for
weaker.

##  5\. Summary

We have examined the iterative RL deblurring method and partially other
methods in the context of the utility of deblurring in nuclear and high-energy
physics. We disregarded explicit noise in the images for the time. SVD of the
blurring matrix is an essential tool in understanding both the blurring and
deblurring processes. For decreasing singular values, the singular vectors
generally exhibit faster variation with pixel position. Symmetries of the
blurring matrix generally give rise to symmetries of the singular vectors and
may make the left- and right-hand vectors coincide. Blurring strength may be
characterized by the falloff of singular values with a singular-vector index.
A strong enough blurring can give rise to a null space in the SVD. A model
one-dimensional system of 12 pixels with a periodic boundary condition has
served to illustrate our points. For high-contrast images, where the number of
pixels at zero or near-zero intensity, compared to maximal in the image, is
comparable to the null-space dimension, the nonnegativity constraint built
into a deblurring can help restore the null-space content in the processed
image. For lower-contrast images, a regularization acting against multiple
extrema in the restored image helps control the null-space content. However, a
too-strong regularization may start affecting the restored image like the
blurring one tries to correct.

## References

  * [1] J. Biemond, R.L. Lagendijk and R.M. Mersereau  “Iterative methods for image deblurring” Conference Name: Proceedings of the IEEE  In _Proc. IEEE_ 78.5, 1990, pp. 856–883  DOI: 10.1109/5.53403
  * [2] G. D’Agostini  “A multidimensional unfolding method based on Bayes’ theorem”  In _Nucl. Instrum. Methods Phys. Res. A_ 362.2, 1995, pp. 487–498  DOI: 10.1016/0168-9002(95)00274-X
  * [3] Pawel Danielewicz and Mizuki Kurata-Nishimura  “Deblurring for nuclei: 3D characteristics of heavy-ion collisions” Publisher: American Physical Society  In _Phys. Rev. C_ 105.3, 2022, pp. 034608  DOI: 10.1103/PhysRevC.105.034608
  * [4] Nicolas Dey et al.  “Richardson–Lucy algorithm with total variation regularization for 3D confocal microscope deconvolution” _eprint: https://analyticalsciencejournals.onlinelibrary.wiley.com/doi/pdf/10.1002/jemt.20294  In _Microsc. Res. Tech._ 69.4, 2006, pp. 260–266  DOI: https://doi.org/10.1002/jemt.20294
  * [5] Per Christian Hansen, James G. Nagy and Dianne P. O’Leary  “Deblurring Images”, Fundamentals of Algorithms  Society for IndustrialApplied Mathematics, 2006  URL: https://epubs.siam.org/doi/book/10.1137/1.9780898718874
  * [6] M.. Khan  “Iterative Methods of Richardson-Lucy-Type for Image Deblurring”  In _NMTMA_ 6.1, 2013, pp. 262–275  DOI: 10.4208/nmtma.2013.mssvm14
  * [7] L. Landweber  “An Iteration Formula for Fredholm Integral Equations of the First Kind” Publisher: Johns Hopkins University Press  In _Am. J. Math._ 73.3, 1951, pp. 615–624  DOI: 10.2307/2372313
  * [8] L.. Lucy  “An iterative technique for the rectification of observed distributions”  In _Astron. J._ 79, 1974, pp. 745  DOI: 10.1086/111605
  * [9] Pierre Nzabahimana et al.  “Deconvoluting experimental decay energy spectra: The 26O case” Publisher: American Physical Society  In _Phys. Rev. C_ 107.6, 2023, pp. 064315  DOI: 10.1103/PhysRevC.107.064315
  * [10] Lothar Reichel  “Introduction to Numerical Computing I. Lecture Notes”, 2013  URL: https://www.math.kent.edu/~reichel/courses/intr.num.comp.1/syllabus.html
  * [11] William Hadley Richardson  “Bayesian-Based Iterative Method of Image Restoration” Publisher: Optical Society of America  In _J. Opt. Soc. Am., JOSA_ 62.1, 1972, pp. 55–59  DOI: 10.1364/JOSA.62.000055
  * [12] Fagun Vankawala, Amit Ganatra and Amit Patel  “A Survey on Different Image Deblurring Techniques” Publisher: Foundation of Computer Science (FCS)  In _Int. J. Comput. Appl._ 116.13, 2015, pp. 15–18  URL: https://www.ijcaonline.org/archives/volume116/number13/20396-2697
  * [13] J. Vargas, J. Benlliure and M. Caamaño  “Unfolding the response of a zero-degree magnetic spectrometer from measurements of the $\Delta$ resonance”  In _Nucl. Instrum. Methods Phys. Res. A_ 707, 2013, pp. 16–25  DOI: 10.1016/j.nima.2012.12.087
  * [14] JunHuai Xu et al.  “Reconstruction of Bremsstrahlung $\gamma$-rays Spectrum in Heavy Ion Reactions with Richardson-Lucy Algorithm”  arXiv, 2024  DOI: 10.48550/arXiv.2405.06711
  * [15] Christian D. Zuninga  “Singular Value Decomposition for Imaging Applications”  SPIE, 2021  URL: https://spie.org/Publications/Book/2611523

Received xxxx 20xx; revised xxxx 20xx; early access xxxx 20xx.