capbtabboxtable[][]

Aiming for Tops of ALPs
with a Muon Collider

So Chigusa1,2, Sudhakantha Girmohanta3,4,
Yuichiro Nakai3,4 and Yufei Zhang3,4
1Berkeley Center for Theoretical Physics, Department of Physics,
University of California, Berkeley, CA 94720, USA
2Theoretical Physics Group,
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA
3Tsung-Dao Lee Institute, Shanghai Jiao Tong University,
520 Shengrong Road, Shanghai 201210, China
4School of Physics and Astronomy, Shanghai Jiao Tong University,
800 Dongchuan Road, Shanghai 200240, China

Future muon colliders with center-of-mass energy of $\mathcal{O}(1-10)$ TeV
can provide a clean high-energy environment with advantages in searches for
TeV-scale axion-like particles (ALPs), pseudo-Nambu–Goldstone bosons
associated with spontaneously broken global symmetries, which are widely
predicted in physics beyond the Standard Model (SM). We exploit ALP couplings
to SM fermions, and guided by unitarity constraints, build a search strategy
focusing on the ALP decay to top quark pairs at muon colliders. It is found
that a large parameter space of TeV-scale ALPs with TeV-scale decay constants
can be probed by utilizing the ALP-top quark coupling.

##  1 Introduction

Axion-like particles (ALPs) are pseudo-Nambu–Goldstone bosons (pNGBs)
associated with spontaneously broken global symmetries. The exploration of
such particles is inspired by the QCD axion [1, 2, 3] to solve the strong CP
problem [4, 5], but unlike the QCD axion, their masses and decay constants are
independent parameters and can vary over a wide range. ALPs are predicted by
many extensions of the Standard Model (SM) and naturally arise in string
theory [6, 7]. They have various phenomenological, cosmological, and
astrophysical implications that have been extensively studied (for reviews,
see $e.g.$ refs. [8, 9]). Therefore, ALPs are considered to be one of the most
important targets in modern particle physics, and their discoveries will make
tremendous progress in the understanding of our universe or multiverse.

Among a wide range of ALP models, our focus in the present paper is on ALPs
whose masses and decay constants are around the TeV scale. Such TeV-scale ALPs
are well-motivated by solutions to the so-called axion quality problem that is
a challenge for the Peccei-Quinn (PQ) mechanism to solve the strong CP problem
[10, 11, 12, 13, 14]. A spontaneously broken global symmetry, resulting in the
QCD axion, must be very exact, otherwise, the axion potential would be
minimized at an inappropriate value to cancel the strong CP phase. However,
any global symmetry is expected to be explicitly broken by quantum gravity
effects. Among various possibilities, one approach to the axion quality
problem is to increase the axion mass without destroying the PQ solution to
the strong CP problem [15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27,
28]. If the axion is heavy enough, quantum gravity effects are no longer
problematic, and the axion decay constant can be set to the TeV scale. Another
interesting possibility is to consider a warped extra dimension model with
three 3-branes where the PQ symmetry is spontaneously broken on the
intermediate 3-brane and the electroweak symmetry is broken by Higgs fields
localized on the IR brane [29] (for recent developments of multi-brane models,
see $e.g.$ refs. [30, 31]). This model can address the electroweak naturalness
problem and the strong CP problem with a high-quality axion at the same time.
A unique prediction of the model is the presence of Kaluza-Klein resonances
associated with the invisible QCD axion, which are TeV-scale ALPs with TeV-
scale decay constants. Since the parameter space of the ALP mass and decay
constant that the Large Hadron Collider (LHC) can explore is limited, new
high-energy machines are needed to probe such TeV-scale ALPs.

A future muon collider with center-of-mass (CM) energy of $\mathcal{O}(1-10)$
TeV can furnish a clean high-energy environment with advantages in searches
for TeV-scale ALPs. Since the muon is much heavier than the electron and loses
much less energy by the synchrotron radiation, a multi-TeV muon beam with a
high luminosity can be realized in a relatively small circular ring. A muon
collider is also able to provide clean and precise collisions of point-like
particles, unlike hadron colliders that involve composite particles. Due to
the short muon lifetime of $2.2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$
at rest, it is required to produce a large number of muons, cool them down to
reduce their emittance, and accelerate them quickly before they decay. Recent
technological developments [32, 33, 34] make it possible to solve those issues
and move forward to the realization of a high-energy muon collider. Therefore,
a muon collider has currently received strong support and interest from the
particle physics community. For recent studies on the physics reach for a
future muon collider, see refs. [35, 36, 37, 38] and references therein.

Depending on the structure of ALP couplings to SM particles, ALP production
and decay at a high-energy muon collider can vary significantly. The authors
of refs. [39, 40] have initiated studies on searches for TeV-scale ALPs
coupling to the electroweak gauge bosons, demonstrating that a muon collider
can substantially expand the ALP mass coverage. In the present paper, we focus
on ALP couplings to SM fermions, which are naturally expected to be present in
the axion effective theory. Since the ALP-fermion coupling is proportional to
the fermion mass, TeV-scale ALPs predominantly decay into top quark pairs in
this scenario. We establish a search strategy for TeV-scale ALPs accordingly;
it is demonstrated that the number of jets and reconstructed top candidates
can serve as filters to distinguish the signal events we are interested in. We
elaborate on the detailed procedure for reducing the principal backgrounds,
including beam-induced fake jets. Furthermore, we reconstruct the dijet
invariant mass distribution for the outgoing top quark pairs. Searching for a
peak structure in the distribution helps to improve the sensitivity, and also
provides information about the ALP mass. We employ a likelihood analysis to
evaluate the significance of the peak observed in the dijet invariant mass
distribution for numerous ALP masses. It is shown that a large parameter space
can be probed by utilizing the ALP-top quark coupling.

The rest of the paper is organized as follows. Section 2 introduces the ALP
model, presents our setup of ALP couplings to SM particles, and investigates
the ALP signal at a muon collider. Then, in section 3, we explain our collider
simulations, discuss the event selection strategy, and explain our statistical
treatment. The results are presented in section 4. Section 5 is devoted to
conclusions and discussions.

##  2 Setup

This section establishes the effective field theoretic framework for the ALP
Lagrangian while elucidating the approximate unitarity constraints associated
with its couplings. Subsequently, we consider the production channels within a
prospective muon collider, focusing on the direct couplings of ALPs with
fermions.

###  2.1 The effective field theory framework

We define an effective field theory (EFT) by extending the SM with a CP-odd
pNGB $a$, which acts as the ALP and are associated with the breaking of a
$U(1)$ global symmetry at a scale $f_{a}$. We will work in the SM effective
field theory framework, where the electroweak symmetry is realized linearly.
The $U(1)$ breaking is assumed to be above the scale of the electroweak
symmetry breaking, i.e., $f_{a}>250$ GeV so that the effective operators are
full $SU(3)_{c}\otimes SU(2)_{L}\otimes U(1)_{Y}$ invariant and ALP couplings
to the EW gauge bosons become important for a future muon collider where the
vector boson fusion (VBF) processes can be substantial. The Lagrangian can be
expanded in powers of $a/f_{a}$, to the leading order,

${\cal L}_{\rm eff}={\cal L}_{\rm
SM}+\frac{1}{2}(\partial_{\mu}a)(\partial^{\mu}a)-\frac{1}{2}m_{a}^{2}a^{2}+\delta{\cal
L}^{(a)}\ .$ ( 2.1)

Here, ${\cal L}_{\rm SM}$ denotes the ordinary SM Lagrangian, $m_{a}$ is the
ALP mass, and

$\delta{\cal L}^{(a)}=c_{\widetilde{W}}{\cal
A}_{\widetilde{W}}+c_{\widetilde{B}}{\cal
A}_{\widetilde{B}}+c_{\widetilde{G}}{\cal
A}_{\widetilde{G}}+c_{a\Phi}O_{a\Phi}^{\psi}\ ,$ ( 2.2)

with

$\displaystyle O_{a\Phi}^{\psi}$
$\displaystyle=i\frac{a}{f_{a}}\sum_{\psi=Q,L}\bar{\psi}_{L}\mathbf{Y_{\psi}}\mathbf{\Phi}\sigma_{3}\psi_{R}+{\rm
h.c.}\ ,$ ( 2.3)

where $Q_{R}\equiv(u_{R},d_{R})$, $L_{R}\equiv(0,e_{R})$, $\sigma_{3}$ acts on
the weak isospin, ${\mathbf{Y}_{Q}}\equiv{\rm diag}(Y_{U},Y_{D})$,
$\mathbf{Y}_{L}\equiv{\rm diag}(0,Y_{E})$, and
$\mathbf{\Phi}\equiv(\tilde{\Phi},\Phi)$ with $\Phi$ being the SM Higgs
doublet and $\tilde{\Phi}=i\sigma_{2}\Phi^{\ast}$. Furthermore, $Y_{U}$,
$Y_{D}$, and $Y_{E}$ represent the Yukawa coupling matrices in the flavor
space for the up-type quarks, down-type quarks, and charged leptons,
respectively. In eq. (2.2),
$c_{\widetilde{W},\widetilde{B},\widetilde{G},a\Phi}$ are dimensionless
coefficients and we have used the notation,

${\cal A}_{\widetilde{X}}=-X_{\mu\nu}\widetilde{X}^{\mu\nu}\frac{a}{f_{a}}\ ,$
( 2.4)

for a gauge field $X\in\\{B,W,G\\}$ with the field strength tensor
$X^{\mu\nu}$ and the dual tensor
$\widetilde{X}^{\mu\nu}=\epsilon^{\mu\nu\rho\sigma}X_{\rho\sigma}$. For the
corresponding Feynman rules, see Appendix B of ref. [41].

ALP couplings to the SM vector bosons have played a central role in building
ALP search strategies at a future muon collider [39, 40]. On the other hand,
the couplings of the ALP to SM fermions were not considered. However, these
couplings are natural, and not prohibited by any symmetry principles.
Furthermore, the couplings are constrained less stringently by the unitarity
bound from partial wave analysis [42]. The helicity amplitude for fermion
scattering processes such as $u_{k}u_{k^{\prime}}\to Za$, where $k,k^{\prime}$
are flavor indices, places partial-wave unitarity bounds on the direct fermion
coupling, the most restrictive of which is [42]

$\left|c_{a\Phi}\right|\lesssim 6\left(\frac{f_{a}}{\rm
TeV}\right)\left(\frac{5\ \rm TeV}{\sqrt{s}}\right)\ .$ ( 2.5)

The couplings to vector bosons ($c_{\tilde{G}},c_{\tilde{B}},c_{\tilde{W}}$)
are constrained more stringently from the partial wave unitarity bounds of the
scatterings of vector bosons via the ALP (e.g., $ZZ\to ZZ$).

Figure 1:  Variation of ALP production cross-sections with CM energy for
different channels. We only have non-zero $c_{a\Phi}$, where all the other
couplings are set to zero. The plot demonstrates the best channel for the
search, namely $\mu^{+}\mu^{-}\to at\bar{t}$.

In the present paper, we initiate the study of ALP search strategies based on
the ALP couplings to the SM fermions at a muon collider by focusing on a
simplified limit in which the ALP-gauge boson couplings are dialed down to
zero, while the effect of the direct couplings to the SM fermions remains. A
generic feature of the ALP couplings to the SM fermions, owing to the shift
symmetry, is that these couplings are proportional to the corresponding
fermion masses. Therefore, an ALP in the TeV mass range dominantly decays to a
$t\bar{t}$ pair. This fact drives the search strategy we build in the current
study.

###  2.2 Collider phenomenology

Figure 2: A few representative Feynman diagrams for the signal event
$\mu^{+}\mu^{-}\to at\bar{t}$. The final state ALP predominantly decays to a
$t\bar{t}$ pair.

In figure 1, we show the production cross-sections of the ALP for a few
representative channels and their variation with CM energy $\sqrt{s}$ of a
muon collider. As can be seen in the figure, the top-associated production
$\mu^{-}\mu^{+}\to t\bar{t}a(+g)$, with the matching of jets [43] to avoid
double counting, is the dominant production channel for the ALP (see figure 2
for a few representative Feynman diagrams). Note that there are also VBF
contributions to the final states including $t\bar{t}a(+g)$, which has a
cross-section as large as $\mathcal{O}(10)\,\mathrm{\%}$ of that of
$\mu^{-}\mu^{+}\to t\bar{t}a(+g)$ in spite of the more severe phase space
suppression. This can be interpreted as a consequence of the large log factor
$\mathop{\rm ln}(s/\Lambda_{W}^{2})$ with
$\Lambda_{W}=\mathcal{O}(100)\,\mathrm{GeV}$ being the weak scale, which
exhibits another role of the muon collider as a vector boson collider. In
addition, we have checked that the processes with more than one gluon emission
only give $\lesssim 10\,\mathrm{\%}$ correction to the cross-section for
$\sqrt{s}=5\,\mathrm{TeV}$. Note, however, that these additional QCD processes
can be important for larger $\sqrt{s}$. In the present work, we focus only on
the leading production channel $at\bar{t}(+g)$.

Since the ALP couples to the SM fermion with a coupling strength proportional
to the Yukawa coupling as shown in eq. 2.3, the ALP most likely decays into a
$t\bar{t}$ pair dominantly, when it is kinematically allowed. Moreover, the
ALP decays promptly for models with large fermionic couplings of our interest.
For example, an illustrative choice for the parameter set,
$m_{a}=1\,\mathrm{TeV}$ and $|c_{a\Phi}|/f_{a}=1/(50\,\mathrm{TeV})$, leads to
the branching ratio $\mathrm{Br}\,(a\to t\bar{t})>99.9\,\%$ with the partial
decay width $\Gamma(a\to t\bar{t})\sim 10^{-3}\,\mathrm{GeV}$, which
corresponds to the typical flight length $c\tau_{a}\sim 200\,\mathrm{fm}$ with
$c$ and $\tau_{a}$ being the speed of light and the ALP lifetime,
respectively. Thus, combined with the production process $at\bar{t}$ of our
interest, we need to consider a way to look for events with four top quarks.

As explained in detail below, we focus on the hadronic decays of top quarks
and identify the boosted top-jet candidate by using the reconstructed jet
mass. This method effectively selects the top-initiated jets, in which all
three hadronic decay products of a single top are within the jet cone. The
method works well because of the high collision energy of a muon collider, for
example, with $\sqrt{s}=$5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$; in
contrast, for four-top events at the LHC with
$\sqrt{s}=$14\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$, the boosted top-jet
reconstruction is not very efficient since the parton-level collision energy
is typically $\sqrt{s^{\prime}}\lesssim$ a few
$\mathrm{T}\mathrm{e}\mathrm{V}$, and the analysis with leptonic decays is
more preferable [44, 45, 46]. In our analysis, we require at least three top-
jet candidates in the event. Therefore, most of the background events come
from the $t\bar{t}t\bar{t}$ production within the SM, which does not have an
axion as an intermediate state. As we will show later in section 3.2, some
background events also come from $t\bar{t}+$jets channels, where a non-top-jet
can be misidentified as a top candidate. For reference, table 2 shows the
cross-sections of processes relevant to our analysis for the benchmark model
of table 1.

$c_{\tilde{W}}$ | $c_{\tilde{B}}$ | $c_{\tilde{G}}$ |  $|c_{a\Phi}|/f_{a}$ [TeV-1] |  $m_{a}$ [TeV]
---|---|---|---|---
$0$ | $0$ | $0$ | $6$ | $1$
Table 1: A benchmark model parameters we use in our analysis. The coupling
value is the maximum of the allowed values for
$\sqrt{s}=$5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$ from the unitarity bound
(2.5).

##  3 Analysis

We now describe the event generation process for the signal together with some
comments on the reduction of fake jets sourced from the beam-induced
backgrounds (BIBs). Then we outline our event selection strategy, considering
various potential backgrounds that may arise there. Principal backgrounds are
identified and the way to reduce them is explained. Finally, we outline the
statistical method to estimate the signal significance.

###  3.1 Event generation

In the following discussion, we choose $\sqrt{s}=5\,\mathrm{TeV}$,
$\mathcal{L}=$100\text{\,}\mathrm{a}\mathrm{b}^{-1}$$, and benchmark
parameters shown in table 1 for an illustrative purpose. The events are first
generated by MadGraph5 [47], which is further processed by Pythia8 [48] for
matching of jets, showering and hadronization. The jet clustering is performed
by FastJet [49] with the default clustering parameter $R=0.5$ and the detector
simulation by Delphes [50] using the muon collider detector card [51]. To
generate the signal events with an axion in an intermediate state, the model
Lagrangian is implemented in FeynRules [52], which then generates a model file
compatible with MadGraph5.

| $t\bar{t}a$ | $t\bar{t}t\bar{t}$ | $t\bar{t}W^{+}W^{-}$ | $t\bar{t}h$ | $t\bar{t}Z$
---|---|---|---|---|---
$\sigma$[ab] | 90 | 23 | 465 | 324 | 158
Table 2:  Cross-sections for the signal and principal backgrounds for the
chosen benchmark parameters in table 1.

As the muon lifetime is $2.2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{s}$
[53] in its rest frame, the decay products of the muon beam and their
secondary interactions give rise to the BIBs that pose a major challenge in
distinguishing true jets from the beam-induced fake jets. The fake jets have
different characteristics from the true jets, which can be exploited to
separate them from each other. For example, most of the BIBs have low
transverse momenta, asynchronous time of arrival, displaced origin, and high
absolute value of pseudo-rapidity. Therefore, the following filtering
procedure during the jet reconstruction is shown to efficiently reduce the
effect of the BIBs [54]:

  * •

Tracks are filtered by requiring a number of Vertex Detector hits greater than
$3$ and a number of Inner Tracker hits greater than $2$.

  * •

The normalized calorimeter hit time $t_{N}$ should be within a time window
$|t_{N}|<$250\text{\,}\mathrm{p}\mathrm{s}$$. Here, $t_{N}\equiv t-t_{0}-cD$,
with $t$ being the absolute hit time, $t_{0}$ the collision time, and $D$ the
hit distance from the collision point.

  * •

Each jet should contain at least one track.

Furthermore, the following kinematic cuts are useful to reduce the number of
fake jets:

  * •

The pseudo-rapidity cut $|\eta|<1.5$ to reject most of the fake jets that lie
along the beam axis.

  * •

The transverse momentum cut with $p_{\rm
T}>$50\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$.

According to ref. [54], the filtering procedure has only a mild effect on the
real jets with $|\eta|<1.5$ and $p_{\rm
T}>$50\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$, leading to the reconstruction
efficiency of roughly $\sim 90\%$, while combined with the kinematic cuts it
reduces the number of fake jets by orders of magnitude. Thus, in our analysis,
we simply assume that the above procedure reduces the effect of the BIBs to a
negligible level, and adopt $90\%$ of real jets that pass the kinematic
cuts.111 Our approach can be regarded as an approximation that the filtering
procedure does not drastically change the kinematic distribution of
reconstructed jets. For a more precise treatment, including the effect of the
remnant fake jets, we need a full simulation of the BIBs and detector
response, which is beyond the scope of the present paper.

###  3.2 Event selection

Our general strategy for the event selection is as follows. Firstly, we look
for four-top events that contain both contributions from the axion signal and
the SM processes. In our analysis, we focus on the hadronic decay of top
quarks and identify candidates of boosted top jets using the jet mass
information. Then, we utilize the number of top candidates and the number of
jets to trim the events with relatively smaller number of top quarks. Finally,
we search for a resonance peak structure in the dijet invariant mass
distribution of top candidates. Note that the use of boosted top jets under a
sufficiently high $\sqrt{s}$ allows us to reconstruct the resonance peak,
which is distinct from the search strategy using the leptonic decays in, e.g.,
ref. [44].

(a) $n_{\mathrm{top}}$

(b) $n_{\mathrm{jet}}$

Figure 3:  Distribution of signal and background events with respect to the
number of top candidates $n_{\mathrm{top}}$ (left) and jets $n_{\mathrm{jet}}$
(right) defined in the text.

To demonstrate what kind of cuts are useful to reduce the number of background
events, we first define the number of top candidates $n_{\mathrm{top}}$ and
jets $n_{\mathrm{jet}}$ for each event. Firstly, $n_{\mathrm{jet}}$
corresponds to jets that have passed $|\eta_{\rm jet}|<1.5$ and $p_{\rm T,\rm
jet}>$50\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$ criteria and have been
chosen randomly with $90\%$ efficiency as discussed above. Then, among
$n_{\mathrm{jet}}$ jets in an event, we identify those with mass $m_{j}$ close
to the top mass $m_{t}\simeq$173\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$ [53]
as top candidates, which counts to be $n_{\mathrm{top}}$. We use the mass
window
$$140\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$<m_{j}<$220\text{\,}\mathrm{G}\mathrm{e}\mathrm{V}$$
for our analysis. In fig. 3, we show the histograms of $n_{\mathrm{top}}$
(left) and $n_{\mathrm{jet}}$ (right) from the ALP signal as well as some of
top-rich SM processes. From the left panel, it is clear that all the SM
processes with only two top quarks tend to have $n_{\mathrm{top}}\leq 2$ as
expected, so we can significantly reduce the number of background events with
the cut $n_{\mathrm{top}}\geq 3$. Note that even four-top processes tend to
have $n_{\mathrm{top}}<4$ due to the reconstruction efficiency of boosted top
jets smaller than $1$. From the right panel, it can be seen that the number-
of-jet cut, e.g., $n_{\mathrm{jet}}\geq 4$ also helps us to reduce backgrounds
while maintaining a large part of the signal events.

cut processes | $t\bar{t}a$ | $t\bar{t}t\bar{t}$ | $t\bar{t}W^{+}W^{-}$ | $t\bar{t}h$ | $t\bar{t}Z$
---|---|---|---|---|---
origin | 9015 | 2285 | 46510 | 32380 | 15790
$n_{\rm top}\geq 3$ | 495 | 82 | 37 | 26 | 3
$n_{\rm jet}\geq 4$ | 417 | 65 | 19 | 0 | 0
Table 3:  Cutflow for the signal process with benchmark parameters in table 1
and top-rich SM processes for
$\sqrt{s}=$5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$.

Tables 3 and 4 show our cutflow for a simulation of the ALP signal and top-
rich SM processes.222 It turns out that contributions from the VBF processes
are negligibly small after our cutflow for both the signal and background
events. Thus, we neglect these contributions in the tables and the following
analysis.  Here, we adopt the ALP parameters
$m_{a}=$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$ and $|c_{a\Phi}|/f_{a}=6$
TeV-1 as shown in table 1. From tables 3 and 4, we can see that the
$n_{\mathrm{top}}\geq 3$ cut significantly reduces the number of events from
processes with only two top quarks. The preceding cut $n_{\mathrm{jet}}\geq 4$
further reduces contaminations from $t\bar{t}W^{+}W^{-}$ and $t\bar{t}h$
processes, slightly making the signal significance larger. After all cuts are
imposed, most of the events are originated from the four-top processes either
through the ALP generation or purely SM processes. Note that, from comparison
between two tables, it is clear that the top-jet reconstruction efficiency
becomes higher when we increase $\sqrt{s}$ from
$5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$ to
$10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$.

From the cutflow, we can order estimate the reconstruction efficiency of a top
quark as a top candidate to be of $\mathcal{O}(10)\,\%$, while the
misidentification rate of the normal jets to be $\mathcal{O}(1)\,\%$. Using
these values, we estimate how many events survive the cut
$n_{\mathrm{top}}\geq 3$, focusing on the SM processes with less than two top
quarks in the final state, and conclude that contaminations from these
processes are negligible for our analysis. Thus, the top-rich processes
considered in tables 3 and 4 are the dominant background processes in our
analysis. Then, we can see from the tables that the number of ALP signal
events can be significantly larger than that of the background events after
the cutflow, indicating the efficacy of our approach.

cut processes | $t\bar{t}a$ | $t\bar{t}t\bar{t}$ | $t\bar{t}W^{+}W^{-}$ | $t\bar{t}h$ | $t\bar{t}Z$
---|---|---|---|---|---
origin | 6625 | 1509 | 34970 | 12300 | 67730
$n_{\rm top}\geq 3$ | 784 | 146 | 49 | 9 | 61
$n_{\rm jet}\geq 4$ | 757 | 136 | 24 | 0 | 7
Table 4:  Same as table 3, but for
$\sqrt{s}=$10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$.

Next, we search for a Gaussian peak in the dijet invariant mass distribution
of top jets produced by the ALP decay. Since there are more than two top
candidates in the events that passed our cutflow, we need to select two
candidates that are most likely sourced from the ALP decay. We have found that
the pair of two top quarks closest to each other in the angular space is a
good choice. More specifically, we define the variable $\Delta R$, which is
the angular separation of the two jets calculated as

$\Delta R=\sqrt{(\Delta\eta)^{2}+(\Delta\phi)^{2}}\ ,$ ( 3.1)

where $\Delta\eta$ denotes the pseudo-rapidity difference of the two jets and
$\Delta\phi$ is the azimuthal angle between them. We calculate $\Delta R$ for
every two-jet pair and select the one with the least $\Delta R$. Fig. 4 shows
the distribution of the dijet invariant mass $m_{jj}$ of the pair of jets
selected as such. We clearly see a peak at
$m_{jj}\sim$1\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$ for the signal events,
corresponding to the ALP mass of our choice, while the smooth distribution for
$m_{jj}\neq m_{a}$ is expected to be due to the failure in the choice of two
top jets originated from the ALP decay. We find a systematic shift of the peak
position slightly lower than $m_{a}$, which might come from the missing
$p_{T}$ associated with invisible partons inside jets. However, we comment
that the peak position and height could be used to reconstruct the ALP mass
$m_{a}$ and its coupling with a further understanding of the systematic shift
by, e.g., a more profound Monte Carlo simulation of the process.

Figure 4:  Dijet invariant mass distribution of the two top candidates
selected as described in the text.

###  3.3 Statistical treatment

In order to quantify the signal significance taking into account the Gaussian
peak structure in the dijet mass distribution shown in fig. 4, we perform
fitting of the distribution. Let $\lambda_{i}$ be the expected number of
events in the $i$-th bin of the histogram. For the current setup, we decompose
$\lambda_{i}$ as

$\displaystyle\lambda_{i}=b_{i}+s_{i}\;,$ ( 3.2)

where $b_{i}$ denotes the smooth distribution that we model with a cubic
function,

$\displaystyle b_{i}=\sum_{k=0}^{3}\beta_{k}(m_{jj}^{(i)})^{k}\;,$ ( 3.3)

with $m_{jj}^{(i)}$ being a representative value of $m_{jj}$ for the $i$-th
bin and $\beta_{k}$ ($k=0\sim 3$) being coefficients, while $s_{i}$ denotes
the peak structure on top of it that we model as

$\displaystyle
s_{i}=A\exp\left(-\frac{\left(m_{jj}^{(i)}-m_{0}\right)^{2}}{2\sigma^{2}}\right)\;.$
( 3.4)

Here, $A$, $m_{0}$, and $\sigma$ are fitting parameters.

Practically, we first prepare the simulation data for each choice of $m_{a}$
and the coupling. We fit the data by $\lambda_{i}$ with fitting parameters
being $\beta_{k}$ ($k=0\sim 3$), $A$, $m_{0}$, and $\sigma$. Let
$\tilde{b}_{i}$ and $\tilde{s}_{i}$ be the values of $b_{i}$ and $s_{i}$,
respectively, evaluated with the best-fit values of these parameters. Defining
their combination,

$\displaystyle\tilde{\lambda}_{i}(\mu)=\tilde{b}_{i}+\mu\tilde{s}_{i}\;,$ (
3.5)

we treat $\lambda_{i}(\mu=1)$ to be the correct expected number of events for
the corresponding ALP model to reduce statistical uncertainties of the
simulation. Note that $\tilde{b}_{i}$ not only contains a contribution from
the SM processes, but also contains a smooth contribution from the process
with an intermediate ALP. This reduces the contribution of the ALP signal to
$\tilde{s}_{i}$, but is useful to make the whole procedure robust to any
systematic uncertainties that could modify the di-jet mass distribution
smoothly. Experimentally, $\tilde{b}_{i}$ and $\tilde{s}_{i}$ correspond to
the fit result of the sideband region with the smooth distribution (3.3) and
that of the whole region with the full expression (3.2).

Let $o_{i}$ be the histogram data obtained from data. In our analysis, we use
$o_{i}=\tilde{\lambda}_{i}(\mu=1)$ as a representative data set. We define the
likelihood function as a function of $\mu$ as

$\displaystyle
L(\bm{o};\mu)=\prod_{i=1}^{N_{B}}\frac{e^{-\tilde{\lambda}_{i}(\mu)}\tilde{\lambda}_{i}(\mu)^{o_{i}}}{\Gamma(o_{i}+1)}\;,$
( 3.6)

with $N_{B}$ being the number of bins. We also define a log-likelihood test-
statistic,

$\displaystyle q_{0}=-2\mathop{\rm
ln}\frac{L(\bm{o};\mu=0)}{L(\bm{o};\mu=1)}\;,$ ( 3.7)

where the choice of $\mu=1$ maximizes the denominator. According to Wilk’s
theorem [55], $q_{0}$ asymptotically obeys a chi-squared distribution with one
degree of freedom. Thus, we identify $\sqrt{q_{0}}=5$ as a criteria for the
discovery of the ALP model at the $5\sigma$ confidence level.

Note that the procedure described so far gives a local significance of the
signal for a fixed axion mass $m_{a}$ used throughout the analysis. For the
global $p$-value, one can define the test-statistic $q$ by making the $m_{a}$
dependence of the likelihood function explicit and maximizing the denominator
of eq. 3.7 by varying $\mu$ and $m_{a}$ at the same time. Since the
probability distribution of $q$ is not known a priori, one has to perform a
large number of pseudo-experiments to estimate the probability distribution to
evaluate the $p$-value corresponding to the expected value of $q$. This way to
estimate a global $p$-value is beyond the scope of the current paper and we
leave it as a future task.

##  4 Results

Figure 5: The test-statistic $\sqrt{q_{0}}$ as a function of the ALP mass
$m_{a}$ for $\sqrt{s}=5$ TeV (red solid), and $\sqrt{s}=10$ TeV (blue dashed)
and $\mathcal{L}=$100\text{\,}\mathrm{a}\mathrm{b}^{-1}$$ with other model
parameters chosen from table 1, while the gray band incorporates the
uncertainties in the fitting procedure. The black dot-dashed line denotes
$\sqrt{q_{0}}=5$ corresponding to the $5\sigma$ reach.

In fig. 5, we show the variation of the test-statistic $\sqrt{q_{0}}$ as a
function of different ALP masses while keeping the coupling fixed to
$|c_{a\Phi}|/f_{a}=$6\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}^{-1}$$, which is
the maximally allowed value by unitarity for
$\sqrt{s}=$5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$ as shown in table 1. The
assumed collider setup is
$\sqrt{s}=$5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$
($\sqrt{s}=$10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$) and
$\mathcal{L}=$100\text{\,}\mathrm{a}\mathrm{b}^{-1}$$ for the red solid (blue
dashed) line. The black dot-dashed line denotes the criteria for the $5\sigma$
discovery, which shows that we can probe up to $m_{a}\sim 1.3$ TeV ($2.7$ TeV)
for $\sqrt{s}=5$ TeV (10 TeV).

The nature of the curve is understood as follows. Firstly, for $m_{a}$ barely
above the threshold $2m_{t}$ of having the $a\to t\bar{t}$ decay channel, top
quarks generated from the ALP decay have small momenta in the ALP-rest frame,
and they are highly collimated in the lab frame. Thus, there is a non-
negligible probability that hadronic decay products of two top quarks are
clustered as a single jet. We confirm this by verifying the existence of a
peak at the ALP mass in the jet mass histogram together with a drastic
reduction of $n_{\rm top}=4$ events in our analysis. This explains the
sensitivity decrease for smaller $m_{a}$, which is more severe for
$\sqrt{s}=$10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$ with more collimated
top jets and more background events. Note that this sensitivity loss can be
well compensated by fitting the jet mass distribution since all the signal
events we lose in our analysis should make a more clear peak in the jet mass
distribution.

Secondly, for a larger $m_{a}$ the production cross-section decreases due to
the reduction of the available phase space, which basically sets the upper
limit on $m_{a}$ to which our analysis is sensitive. Now, it is useful to
understand how this compares with a larger CM energy, as although we do not
gain much in the production cross-section as can be inferred from fig. 1, the
top-jet reconstruction efficiency increases thanks to higher boost factors. To
verify this, and to compare with the case of
$\sqrt{s}=$5\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$, we set aside the
unitarity constraint for the moment, and keep the coupling fixed while
increasing CM energy to
$\sqrt{s}=$10\text{\,}\mathrm{T}\mathrm{e}\mathrm{V}$$. Clearly, due to
increased phase space in the larger ALP mass regime, we can probe to higher
values of $m_{a}$, and it is also confirmed that one effectively gets more
sensitivity, thanks to the enhanced top-jet reconstruction efficiency.

Figure 6: Projected $5\sigma$ reach contours in the normalized coupling
$|c_{a\Phi}|/f_{a}$ and $m_{a}$ plane using the four-top channel search for
two different CM energies, namely $\sqrt{s}=5$ TeV (red) and $10$ TeV (blue)
for integrated luminosities ${\cal L}=10$ ab-1 and $100$ ab-1. The approximate
unitarity bound [42] is shown as the red dashed (blue dot-dashed) line for
$\sqrt{s}=5$ TeV (10 TeV).

Let us now turn to the following question: what is the $5\sigma$ reach for the
coupling-mass parameter space, i.e., in the $|c_{a\Phi}|/f_{a}$-$m_{a}$ plane?
Note that it is difficult to answer this analytically as the test-statistic
does not scale with the coupling, and there are non-trivial effects from the
background dependence and top-jet reconstruction efficiencies. Therefore, we
utilize the following argument, namely that although the test-statistic does
not scale with the coupling, the signal cross-section scales as
$\left(c_{a\Phi}/f_{a}\right)^{2}$. We have verified this numerically within
1% accuracy. This is because once the ALP is produced on-shell, for
$m_{a}>2m_{t}$, it immediately decays to two top pairs with $\sim 100\%$
branching fraction. Therefore, we can follow the same technique for evaluating
$\sqrt{q_{0}}$, with the same background as before, while re-scaling the
signal events, i.e., a reduction in the coupling will reduce the signal peak,
and will result in a reduced $\sqrt{q_{0}}$. Now we view this procedure as a
continuous function of the coupling and extract the value of the coupling at
which one gets $\sqrt{q_{0}}=5$ for each fixed value of $m_{a}$. Hence, this
gives the $5\sigma$ reach for the normalized coupling for the chosen ALP mass.
We repeat this procedure for all the different events generated for distinct
$m_{a}$, and thereby get the $5\sigma$ projected reach in the coupling-mass
plane. This is depicted in fig. 6 for two different $\sqrt{s}$ and two
different $\mathcal{L}$. The estimated unitarity bounds are shown for each of
the CM energy as dashed lines, and we can see that one is able to probe a
large parameter space including the region where the EFT description is
expected to breakdown with possibly the underlying physics signature
appearing. For the extra-dimensional KK mode ALP, for example, this could mean
that the contribution from the next KK resonance becomes important, while for
other UV completion, this could mean the fields that had been integrated out
for obtaining the EFT Lagrangian can be excited. We show both the cases for
the integrated luminosities $10\text{\,}\mathrm{a}\mathrm{b}^{-1}$ and
$100\text{\,}\mathrm{a}\mathrm{b}^{-1}$ that help to elucidate the dependence
of the parameter space probe as a function of the luminosity.

##  5 Discussions

Future muon colliders with CM energy of $\mathcal{O}(1-10)$ TeV are able to
provide a clean high-energy environment with advantages in searches for TeV-
scale ALPs. We have exploited ALP couplings to the SM fermions, and guided by
unitarity constraints, built a search strategy focusing on the ALP decay to SM
top quark pairs. We analyzed the extraction of the signal from principal SM
backgrounds, and provided a projected sensitivity reach for the parameter
space of the ALP.

In our current work, we were agnostic regarding the UV completion, but our
search strategy for the four-top channel that we have elucidated here can be
applied to specific models to find the reach of future collider machines.
Furthermore, the ALP-fermion coupling was considered in isolation to
illustrate its effects on the production and decay channel. However, in the
future, it will be interesting to consider the gauge boson couplings together
with the direct fermion couplings we considered here as this will mimic a more
realistic scenario. Regarding the region $2m_{b}<m_{a}<2m_{t}$, although we do
not take up this task in the current work, note that here the $b\bar{b}$ decay
channel for the ALP decay will dominate, and a similar search strategy could
be followed to project a constraint in this regime. We also note that High-
Luminosity (HL)-LHC will have a comparable or better handle at this mass
range, and one can do a comparative analysis of the efficacy of HL-LHC and a
future lepton collider for this regime, and we relegate it to a future
work.333 Four-top events have been detected recently at both ATLAS and CMS
[44, 45, 46].  To conclude, our results show that the ALP-top coupling can be
crucial to probe the ALP in four-top channels, as this coupling is natural
from a UV theory perspective and can provide a complementary or even better
chance at probing it in a future muon collider machine to the more
investigated channels involving gauge boson couplings primarily [39, 40].

## Acknowledgements

SC is supported by the Director, Office of Science, Office of High Energy
Physics of the U.S. Department of Energy under the Contract No. DE-
AC02-05CH1123. YN is supported by Natural Science Foundation of China under
grant No. 12150610465.

## References

  * [1] R. D. Peccei and H. R. Quinn, _CP Conservation in the Presence of Instantons_ , _Phys. Rev. Lett._ 38 (1977) 1440.
  * [2] S. Weinberg, _A New Light Boson?_ , _Phys. Rev. Lett._ 40 (1978) 223.
  * [3] F. Wilczek, _Problem of Strong $P$ and $T$ Invariance in the Presence of Instantons_, _Phys. Rev. Lett._ 40 (1978) 279.
  * [4] C. A. Baker et al., _An Improved experimental limit on the electric dipole moment of the neutron_ , _Phys. Rev. Lett._ 97 (2006) 131801 [hep-ex/0602020].
  * [5] J. M. Pendlebury et al., _Revised experimental upper limit on the electric dipole moment of the neutron_ , _Phys. Rev. D_ 92 (2015) 092003 [1509.04411].
  * [6] P. Svrcek and E. Witten, _Axions In String Theory_ , _JHEP_ 06 (2006) 051 [hep-th/0605206].
  * [7] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper and J. March-Russell, _String Axiverse_ , _Phys. Rev. D_ 81 (2010) 123530 [0905.4720].
  * [8] J. Beacham et al., _Physics Beyond Colliders at CERN: Beyond the Standard Model Working Group Report_ , _J. Phys. G_ 47 (2020) 010501 [1901.09966].
  * [9] K. Choi, S. H. Im and C. Sub Shin, _Recent Progress in the Physics of Axions and Axion-Like Particles_ , _Ann. Rev. Nucl. Part. Sci._ 71 (2021) 225 [2012.05029].
  * [10] R. Holman, S. D. H. Hsu, T. W. Kephart, E. W. Kolb, R. Watkins and L. M. Widrow, _Solutions to the strong CP problem in a world with gravity_ , _Phys. Lett. B_ 282 (1992) 132 [hep-ph/9203206].
  * [11] M. Kamionkowski and J. March-Russell, _Planck scale physics and the Peccei-Quinn mechanism_ , _Phys. Lett. B_ 282 (1992) 137 [hep-th/9202003].
  * [12] S. M. Barr and D. Seckel, _Planck scale corrections to axion models_ , _Phys. Rev. D_ 46 (1992) 539.
  * [13] S. Ghigna, M. Lusignoli and M. Roncadelli, _Instability of the invisible axion_ , _Phys. Lett. B_ 283 (1992) 278.
  * [14] L. M. Carpenter, M. Dine and G. Festuccia, _Dynamics of the Peccei Quinn Scale_ , _Phys. Rev. D_ 80 (2009) 125017 [0906.1273].
  * [15] S. Dimopoulos, _A Solution of the Strong CP Problem in Models With Scalars_ , _Phys. Lett. B_ 84 (1979) 435.
  * [16] B. Holdom and M. E. Peskin, _Raising the Axion Mass_ , _Nucl. Phys. B_ 208 (1982) 397.
  * [17] M. Dine and N. Seiberg, _String Theory and the Strong CP Problem_ , _Nucl. Phys. B_ 273 (1986) 109.
  * [18] J. M. Flynn and L. Randall, _A Computation of the Small Instanton Contribution to the Axion Potential_ , _Nucl. Phys. B_ 293 (1987) 731.
  * [19] V. A. Rubakov, _Grand unification and heavy axion_ , _JETP Lett._ 65 (1997) 621 [hep-ph/9703409].
  * [20] Z. Berezhiani, L. Gianfagna and M. Giannotti, _Strong CP problem and mirror world: The Weinberg-Wilczek axion revisited_ , _Phys. Lett. B_ 500 (2001) 286 [hep-ph/0009290].
  * [21] H. Fukuda, K. Harigaya, M. Ibe and T. T. Yanagida, _Model of visible QCD axion_ , _Phys. Rev. D_ 92 (2015) 015021 [1504.06084].
  * [22] T. Gherghetta, N. Nagata and M. Shifman, _A Visible QCD Axion from an Enlarged Color Group_ , _Phys. Rev. D_ 93 (2016) 115010 [1604.01127].
  * [23] S. Dimopoulos, A. Hook, J. Huang and G. Marques-Tavares, _A collider observable QCD axion_ , _JHEP_ 11 (2016) 052 [1606.03097].
  * [24] P. Agrawal and K. Howe, _Factoring the Strong CP Problem_ , _JHEP_ 12 (2018) 029 [1710.04213].
  * [25] M. K. Gaillard, M. B. Gavela, R. Houtz, P. Quilez and R. Del Rey, _Color unified dynamical axion_ , _Eur. Phys. J. C_ 78 (2018) 972 [1805.06465].
  * [26] A. Hook, S. Kumar, Z. Liu and R. Sundrum, _High Quality QCD Axion and the LHC_ , _Phys. Rev. Lett._ 124 (2020) 221801 [1911.12364].
  * [27] T. Gherghetta, V. V. Khoze, A. Pomarol and Y. Shirman, _The Axion Mass from 5D Small Instantons_ , _JHEP_ 03 (2020) 063 [2001.05610].
  * [28] T. Gherghetta and M. D. Nguyen, _A Composite Higgs with a Heavy Composite Axion_ , _JHEP_ 12 (2020) 094 [2007.10875].
  * [29] S. J. Lee, Y. Nakai and M. Suzuki, _High quality axion via a doubly composite dynamics_ , _JHEP_ 03 (2022) 038 [2112.08083].
  * [30] S. J. Lee, Y. Nakai and M. Suzuki, _Multiple hierarchies from a warped extra dimension_ , _JHEP_ 02 (2022) 050 [2109.10938].
  * [31] S. Girmohanta, S. J. Lee, Y. Nakai and M. Suzuki, _Multi-brane cosmology_ , 2304.05586.
  * [32] MICE collaboration, _Demonstration of cooling by the Muon Ionization Cooling Experiment_ , _Nature_ 578 (2020) 53 [1907.08562].
  * [33] J. P. Delahaye, M. Diemoz, K. Long, B. Mansoulié, N. Pastrone, L. Rivkin et al., _Muon Colliders_ , 1901.06150.
  * [34] N. Bartosik et al., _Detector and Physics Performance at a Muon Collider_ , _JINST_ 15 (2020) P05001 [2001.04431].
  * [35] H. Al Ali et al., _The muon Smasher’s guide_ , _Rept. Prog. Phys._ 85 (2022) 084201 [2103.14043].
  * [36] K. M. Black et al., _Muon Collider Forum Report_ , 2209.01318.
  * [37] C. Aime et al., _Muon Collider Physics Summary_ , 2203.07256.
  * [38] Y. Hamada, R. Kitano, R. Matsudo, H. Takaura and M. Yoshida, _$\mu$ TRISTAN_, _PTEP_ 2022 (2022) 053B02 [2201.06664].
  * [39] T. Han, T. Li and X. Wang, _Axion-Like Particles at High Energy Muon Colliders – A White paper for Snowmass 2021_ , in _2022 Snowmass Summer Study_ , 3, 2022, 2203.05484.
  * [40] Y. Bao, J. Fan and L. Li, _Electroweak ALP Searches at a Muon Collider_ , 2203.04328.
  * [41] I. Brivio, M. B. Gavela, L. Merlo, K. Mimasu, J. M. No, R. del Rey et al., _ALPs Effective Field Theory and Collider Signatures_ , _Eur. Phys. J. C_ 77 (2017) 572 [1701.05379].
  * [42] I. Brivio, O. J. P. Éboli and M. C. Gonzalez-Garcia, _Unitarity constraints on ALP interactions_ , _Phys. Rev. D_ 104 (2021) 035027 [2106.05977].
  * [43] S. Frixione and B. R. Webber, _Matching NLO QCD computations and parton shower simulations_ , _JHEP_ 06 (2002) 029 [hep-ph/0204244].
  * [44] ATLAS collaboration, _Measurement of the t $\overline{t}$t$\overline{t}$ production cross section in $pp$ collisions at $\sqrt{s}$ = 13 TeV with the ATLAS detector_, _JHEP_ 11 (2021) 118 [2106.11683].
  * [45] ATLAS collaboration, _Observation of four-top-quark production in the multilepton final state with the ATLAS detector_ , _Eur. Phys. J. C_ 83 (2023) 496 [2303.15061].
  * [46] CMS collaboration, _Observation of four top quark production in proton-proton collisions at $\sqrt{s}$ = 13 TeV_, 2305.13439.
  * [47] J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer et al., _The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations_ , _JHEP_ 07 (2014) 079 [1405.0301].
  * [48] C. Bierlich et al., _A comprehensive guide to the physics and usage of PYTHIA 8.3_ , 2203.11601.
  * [49] M. Cacciari, G. P. Salam and G. Soyez, _FastJet User Manual_ , _Eur. Phys. J. C_ 72 (2012) 1896 [1111.6097].
  * [50] DELPHES 3 collaboration, _DELPHES 3, A modular framework for fast simulation of a generic collider experiment_ , _JHEP_ 02 (2014) 057 [1307.6346].
  * [51] M. Selvaggi, “Talk at mdi studies meeting of the muon collider collaboration.” https://indico.cern.ch/event/957299/contributions/4023467/attachments/2106044/3541874/delphes_card_mucol_mdi%20.pdf.
  * [52] A. Alloul, N. D. Christensen, C. Degrande, C. Duhr and B. Fuks, _FeynRules 2.0 - A complete toolbox for tree-level phenomenology_ , _Comput. Phys. Commun._ 185 (2014) 2250 [1310.1921].
  * [53] Particle Data Group collaboration, _Review of Particle Physics_ , _PTEP_ 2022 (2022) 083C01.
  * [54] Muon Collider collaboration, _Simulated Detector Performance at the Muon Collider_ , 2203.07964.
  * [55] S. S. Wilks, _The Large-Sample Distribution of the Likelihood Ratio for Testing Composite Hypotheses_ , _The Annals of Mathematical Statistics_ 9 (1938) 60 .


# Deep network series for large-scale high-dynamic range imaging

###### Abstract

We propose a new approach for large-scale high-dynamic range computational
imaging. Deep Neural Networks (DNNs) trained end-to-end can solve linear
inverse imaging problems almost instantaneously. While unfolded architectures
provide necessary robustness to variations of the measurement setting,
embedding large-scale measurement operators in DNN architectures is
impractical. Alternative Plug-and-Play (PnP) approaches, where the denoising
DNNs are blind to the measurement setting, have proven effective to address
scalability and high-dynamic range challenges, but rely on highly iterative
algorithms. We propose a residual DNN series approach, where the reconstructed
image is built as a sum of residual images progressively increasing the
dynamic range, and estimated iteratively by DNNs taking the back-projected
data residual of the previous iteration as input. We demonstrate on
simulations for radio-astronomical imaging that a series of only few terms
provides a high-dynamic range reconstruction of similar quality to state-of-
the-art PnP approaches, at a fraction of the cost.

Index Terms—  computational imaging, deep neural networks, unfolded
architectures, plug-and-play, astronomical imaging.

##  1 Introduction

Designing algorithms for solving ill-posed imaging inverse problems that both
deliver high precision reconstructions and scale to large data volumes is a
significant challenge of imaging sciences. A problem where such considerations
are paramount is aperture synthesis by Radio Interferometry (RI) in astronomy,
an iconic example of linear inverse Fourier imaging problems [1]. In this
context, the acquisition of an unknown target radio image
$\overline{x}\in\mathbb{R}^{N}$ from observed data $z\in\mathbb{C}^{M}$
follows the observation model

$z=H\overline{x}+e,$ (1)

where $e\sim\mathcal{N}(0,\sigma^{2})$ is the realisation of a complex
Gaussian random noise, and where the linear measurement operator
$H\colon\mathbb{R}^{N}\to\mathbb{C}^{M}$ represents (in its simplest form, see
[2] for details) a non-uniform Fourier transform. More precisely, $H=UFZ$,
where $Z$ is a zero-padding operator, $F$ denotes the 2D discrete Fourier
transform, and $U$ is an interpolation matrix mapping the discrete Fourier
coefficients onto the measured values of the continuous Fourier transform. The
set of sampled values is characterised by a so-called sampling pattern. With
the advent of a whole new generation of radio telescopes arrays aiming to
probe the sky with much higher resolution and sensitivity, such as the
flagship MeerKAT (sarao.ac.za) and SKA (skatelescope.org), extreme image
dimensions and data volumes are to be expected, along with dynamic ranges
(_i.e._ the ratio between the brightest and faintest sources in
$\overline{x}$) spanning multiple orders of magnitude.

Deep Neural Networks (DNNs) have shown outstanding results in solving linear
inverse imaging problems. On the one hand, end-to-end approaches provide
extremely fast reconstruction. They are widespread in other imaging
communities [3, 4, 5], but their use remains limited in RI imaging [2, 6, 7].
This is mainly due to the lack of ground-truth datasets, combined with a wide
intrinsic variability of the RI observation model, leading to generalisation
issues. While unfolded architectures [8, 9] provide necessary robustness to
variations of the measurement setting, embedding large-scale measurement
operators in DNN architectures is impractical, both for training and
inference. On the other hand, Plug-and-Play (PnP) algorithms [10, 11, 12],
substituting learned DNN denoisers in lieu of proximal regularisation
operators in optimisation algorithms, have shown outstanding performance and
robustness, including for high-dynamic range imaging [2, 13]. However, PnP
approaches remain highly iterative and will still struggle to scale to the
image sizes and data volumes of interest in applications such as RI imaging.

In this paper, a new residual DNN series approach to large-scale high-dynamic
range computational imaging is devised, where the reconstructed image is built
as a sum of few residual images progressively increasing the dynamic range,
and estimated iteratively as output of DNNs taking the back-projected data
residual of the previous iteration as input. We show on preliminary
simulations for RI imaging that such a strategy may enable to achieve an
imaging quality matching that delivered by the most advanced optimisation and
PnP algorithms, while being orders of magnitude faster thanks to its very
limited iterative nature.

##  2 Methodology

###  2.1 Proposed approach

We propose to train a series of DNNs $(G_{k})_{1\leq k\leq K}$ yielding a
sequence of estimated reconstructions $(x_{k})_{1\leq k\leq K}$ of
$\overline{x}$, following the update rule

$(\forall 1\leq k\leq K)\quad\left\\{\begin{aligned} r_{k-1}&=\gamma
H^{*}(z-Hx_{k-1}),\\\ x_{k}&=x_{k-1}+G_{k}(r_{k-1}),\end{aligned}\right.$ (2)

where $x_{0}=0$ and where $(r_{k-1})_{1\leq k\leq K}$ corresponds to the data
residuals back-projected in the image domain. In (2), we set
$\gamma=1/\operatorname{max}(H^{*}H\delta)$ to ensure the normalisation of
back-projected data residuals at the input of the networks ($\delta$ being the
Dirac delta). The last element of the sequence and final reconstruction
$x_{K}$ also reads

$x_{K}=\sum_{k=1}^{K}G_{k}(r_{k-1}),$ (3)

hence the name “deep network series”. In practice, we aim at keeping $K$ as
small as possible. We also note that the first iteration is equivalent to a
standard end-to-end learning approach, as $x_{1}=G_{1}(r_{0})$ with
$r_{0}=H^{*}z$ the back-projected data. The approach is summarised in Figure
1.

Fig. 1: Summary of the proposed approach. Top row shows the successive image
estimates $x_{k}$ for $1\leq k\leq K$ while the bottom row shows the back-
projected data residuals $r_{k-1}$, input to $G_{k}$, and the predicted image
residuals $x_{k}-x_{k-1}$, output of $G_{k}$. Each network is shown in a
different colour to stress that it is trained separately. In blue, the
computation of $r_{k-1}$. Note that $x_{0}=0$ and $r_{0}=H^{*}z$.

In order to train the sequence of DNNs $(G_{k})_{1\leq k\leq K}$, we start
from a dataset of $L$ synthetic ground-truth images $(\overline{x}_{l})_{1\leq
l\leq L}$ from which we simulate measurements $(z_{l})_{1\leq l\leq L}$ as per
(1). Then, for every $1\leq k\leq K$, we generate a dataset of residuals
$(r_{k-1,l})_{1\leq l\leq L}$ as per (2) and train $G_{k}$ to

$\underset{\theta_{k}}{\text{minimise}}\,\,\frac{1}{L}\sum_{l=1}^{L}\|G_{k_{,}\theta_{k}}(r_{k-1,l})+x_{k-1,l}-\overline{x}_{l}\|_{1},$
(4)

where $\theta_{k}$ denotes the learnable parameters of $G_{k}$.

###  2.2 Relation to unfolding and PnP

The proposed approach is reminiscent of both unfolded and PnP algorithms based
on the Forward-Backward optimisation algorithm [14], as the back-projected
data residuals in (2) can be written as $r_{k-1}=-\gamma\nabla f(x_{k-1})$ for
$f(\cdot)=\frac{1}{2}\|H\cdot-z\|_{2}$. There are however important
differences with such schemes. Firstly, instead of applying $G_{k}$ to
$\operatorname{Id}-\gamma\nabla f$, $G_{k}$ is applied to the residual
$-\gamma\nabla f$ only. Secondly, unfolded architectures are trained in an
end-to-end fashion, whereas our approach proposes to train $G_{k}$
sequentially, with a dataset of back-projected data residuals
$(r_{k-1,l})_{1\leq l\leq L}$ as input. We stress that this procedure,
involving the training of $K$ DNNs, is dictated by the large-scale nature of
$H$, precluding its direct embedding into a DNN architecture. Thirdly, DNNs in
PnP algorithms are generally the same across iterations and trained as simple
denoisers.

In its series expression (3), our approach is reminiscent of Neumann networks
[15]. But the latter are unfolded architectures, embedding a single DNN
regulariser across iterations.

##  3 Experiments

###  3.1 Training RI dataset & training

To circumvent the absence of an appropriate dataset of target radio images, we
rely on a synthetic dataset from [2], where noisy, low-dynamic range optical
images have been processed through denoising and exponentiation procedures to
create a dataset of $L=2235$ clean images $(\overline{x}_{l})_{1\leq l\leq
L}$, of size $N=512\times 512$, and with nominal dynamic range slightly above
$10^{4}$. The full dataset is normalised so that the maximum pixel value
across all images is equal to 1.

For each image $\overline{x}_{l}$ we create data
$z_{l}=H_{l}\overline{x}_{l}+e_{l}$, with a measurement operator $H_{l}$
resulting from simulating a sampling pattern of the MeerKAT telescope, with
fixed integration and observation times, and randomly chosen pointing
direction, leading to data vectors of size $M=1.6\times 10^{6}$. The
measurement noise $e_{l}$ is generated with fixed standard deviation set to
ensure a back-projected noise level
$\sigma/\sqrt{2\|H_{l}^{*}H_{l}\|_{2}}\lesssim 10^{-4}$, thus leading to a
target dynamic range of the order of, or slightly above, $10^{4}$. Following
standard practice in the field, the data $z_{l}$, operator $H_{l}$, and noise
$e_{l}$ are weighted component-wise by the inverse square-root of the sampling
density at each point to mitigate the side-lobes of the point spread function
$H_{l}^{*}H_{l}\delta$ (see Figure 2.c).

We train the sequence of networks $(G_{k})_{1\leq k\leq K}$ following the
procedure detailed in Section 2.1. Each network $G_{k}$ has the same
architecture, namely a UNet [11] in which convolutional blocks were replaced
by WDSR blocks [16]. The networks also include an initial normalisation layer,
ensuring that any input image effectively has zero mean and unit standard
deviation, mitigating potential generalisation issues related to variations of
statistics between the training samples and any test sample. The output image
is de-normalised as part of a last layer, using the same mean and standard
deviation as computed at the input. Each network is trained with the Adam
algorithm, for approximately 200 epochs for each network. In practice, we
noticed that the reconstruction quality saturates for $K>4$ on the training
dataset, and therefore chose $K=4$ for all our experiments.

###  3.2 Test RI dataset & simulation set up

We validate the proposed methodology on a test dataset of 3 radio images
(3c353, Hercules A, and Centaurus A, see [2] for more details) of size
$N=512\times 512$, with peak value normalised to 1, and dynamic range slightly
above $10^{4}$. For each of them, 5 RI data vectors are simulated from 5
MeerKAT sampling patterns generated as in Section 3.1, using the same
integration and observation time, same noise standard deviation, and random
pointing directions. This leads to a total of 15 inverse problems with data
vectors of size $M=1.6\times 10^{6}$.

Reconstruction quality is evaluated with the signal-to-noise ratio metric, as
$\text{SNR}(x,\overline{x})=20\operatorname{log}_{10}(\|\overline{x}\|_{2}/\|\overline{x}-x\|_{2})$
for an estimate $x$ of $\overline{x}$. Since we aim at reconstructing high-
dynamic range images, we visualise reconstructions in logarithmic scale
(_i.e._ after applying the transform $\operatorname{rlog}\colon
x\mapsto\operatorname{log}_{10}(10^{3}x+1)/3$), which helps visualising very
faint features. We also compute the SNR in logarithmic scale as
$\text{logSNR}(x,\overline{x})=\text{SNR}(\operatorname{rlog}(x),\operatorname{rlog}(\overline{x}))$.
When evaluating the reconstruction quality of the proposed method across
iterations, we also use the relative norm of the back-projected data residuals
$\eta_{k}=\|r_{k}\|_{2}/\|r_{0}\|_{2}$, measuring the relative amount of data
left to be processed after each iteration.

###  3.3 Benchmark RI algorithms

We compare the proposed method to various RI imaging algorithms, from the
recent and advanced sparsity-promoting optimisation algorithm SARA [17, 18],
to its PnP counterpart AIRI [2, 13], to the traditional RI imaging algorithm
CLEAN [19], and to an end-to-end approach resulting from cutting the network
series to $K=1$ in the proposed approach.

Interestingly, CLEAN relies on a matching pursuit approach iteratively
identifying model components from a back-projected data residual. From this
perspective our approach is conceptually close to CLEAN. However, while the
simplicity of CLEAN’s implicit sparsity image model makes it scalable to large
image dimensions and data volumes, this comes at the cost of a now commonly
accepted severe sub-optimality in reconstruction quality. Our approach
leverages a fundamentally different, DNN-based, component identification
process from the back-projected data residual at each iteration.

###  3.4 Results

|  |  |  |
---|---|---|---|---
(a) $\overline{x}$ |  (b) $H^{*}z$ | (c) sampling |  |
|  |  |  |
(d) $r_{0}=H^{*}z$ |  (e) $r_{1}$ |  (f) $r_{2}$ |  (g) $r_{3}$ |
|  |  |  |
(h) $x_{1}$ |  (i) $x_{2}$ |  (j) $x_{3}$ |  (k) $x_{4}$ |
$(21.0,11.0)$ | $(27.9,22.3)$ | $(28.1,24.2)$ | $(28.2,25.1)$ |
Fig. 2: Different stages of the proposed approach illustrated on one of the three images (3c353) and for one of the five sampling patterns of the test dataset. Top row shows (a) the ground truth $\overline{x}$, (b) the back-projected data $H^{*}z$, and (c) the weighted Fourier sampling pattern (colorbar displaying weight values). The second row (panels (d)-(g)) and third row (panels (h)-(k)) show the back-projected data residuals $r_{k-1}$ and reconstructions $x_{k}$ for $1\leq k\leq 4$, below which associated (SNR, logSNR) metrics are displayed. |
---|---
Fig. 3: Evolution of average reconstruction quality across iterations $1\leq k\leq 4$ over the test dataset, with error bars indicating $95\%$ confidence intervals. Left: SNR and logSNR metrics for $x_{k}$. Right: Back-projected data residual norm $\eta_{k}$. (a) Ground truth | (b) Back-projected | (c) CLEAN | (d) SARA | (e) AIRI |  (f) $x_{1}$ |  (g) $x_{4}$ |
---|---|---|---|---|---|---|---
|  |  |  |  |  |  |
| $(-11.0,-6.9)$ | $(-2.7,0.2)$ | $(\mathbf{27.0},24.6)$ | $(25.0,24.7)$ | $(20.8,11.8)$ | $(26.0,\mathbf{25.1})$ |
|  |  |  |  |  |  |
| $(-13.5,-6.8)$ | $(-4.6,0.4)$ | $(26.8,\mathbf{24.7})$ | $(25.3,24.3)$ | $(18.0,16.9)$ | $(\mathbf{27.7},23.8)$ |
Fig. 4: Reconstruction results for the various algorithms considered, for the Hercules A (top) and Centaurus A (bottom) test images, each probed through one of the 5 test MeerKAT sampling patterns. (SNR, logSNR) metrics (in dB) are displayed below each reconstruction. Algorithm | SNR | logSNR | time (s.) | iteration #
---|---|---|---|---
CLEAN | $5.3\pm 0.1$ | $8.2\pm 0.4$ | $217\pm 5$ | $8.8\pm 0.2$
SARA | $\mathbf{26.6}\pm 0.8$ | $\mathbf{24.9}\pm 0.3$ | $11291\pm 563$ | $4825\pm 190$
AIRI | $25.1\pm 0.8$ | $\mathbf{25.0}\pm 0.4$ | $4258\pm 140$ | $4515\pm 103$
$x_{1}$ | $20.3\pm 1.1$ | $13.3\pm 1.2$ | $\mathbf{1.0}\pm 0.2$ | $1$
$x_{4}$ | $\mathbf{26.8}\pm 1.1$ | $24.5\pm 0.4$ | $3.9\pm 0.6$ | $4$
Table 1: Reconstruction SNR, logSNR, times, and number of iterations, for the
different imaging algorithms considered. Values reported are averages over the
test dataset, with associated 95% confidence intervals.

Figure 2 shows the different stages of the proposed reconstruction method, on
3c353 and for the sampling pattern shown in panel (c). Figure 3 shows the
average SNR and logSNR metrics and residual norms on the 15 inverse problems
of the test dataset. The first estimate $x_{1}=G_{1}(r_{0})$ provides good
results at high intensities but is not accurate at faint intensities. This
translates into relatively high SNR, but poor logSNR values. The associated
residual $r_{1}$ contains important information and is far from random noise.
The output $x_{2}$ is much more accurate than $x_{1}$, not only at high
intensities (improvement of the SNR), but also at fainter intensities
(improvement of the logSNR). It is also associated with a residual $r_{2}$
containing much less information than $r_{1}$, also evidenced by
$\eta_{2}<\eta_{1}$, supporting the observation that the $x_{2}$
reconstruction is more faithful to the measurements than $x_{1}$. While the
SNR does not evolve significantly between $x_{2}$ and $x_{3}$, with also
$\eta_{3}\simeq\eta_{2}$, it appears clearly that $r_{3}$ contains less
residual signal structure than $r_{2}$, suggesting that $x_{3}$ provides
further improvement at faint intensities. This is evidenced by a noticeable
improvement in logSNR. A final logSNR increase is brought in $x_{4}$.

Table 1 shows the average SNR, logSNR, reconstruction times, and number of
iterations, as well as associated 95% confidence intervals over the 15 test
inverse problems, with comparison to the benchmark algorithms. Figure 4
provides visual reconstruction results of Hercules A and Centaurus A, for 2
sampling patterns different from the one in Figure 2. In terms of imaging
quality, the reported statistics confirm that the proposed approach is on par
with SARA and AIRI, all three methods achieving high resolution and high
dynamic range. The end-to-end approach falls short from delivering a similar
resolution and dynamic range, but is still significantly better than CLEAN. In
terms of reconstruction times, the proposed approach is two orders of
magnitude faster than CLEAN, and three orders of magnitude faster than AIRI,
itself more than twice as fast as SARA. The computation time of each method
crudely evolves proportionally to the number of iterations. More precisely,
the proposed method, the end-to-end approach, and AIRI, share a very similar
computation time per iteration, due to a similar iteration structure
consisting in computing a residual followed by application of a DNN. SARA’s
iteration cost is heavier due to costly proximal operations associated with
its regularisation approach. CLEAN’s iteration cost is also heavier, mainly
due to the fact it does not keep the measurement operator in memory.

We finally note that preliminary results not reported here suggest that the
proposed hybrid DNN architecture used outperforms architectures such as UNet
and WDSR, both in the end-to-end and residual network series approaches.

##  4 Conclusion

In this paper, we proposed a novel residual DNN series approach for solving
large-scale high-dynamic range computational imaging problems, circumventing
scalability challenges encountered by both unfolded end-to-end architectures
and PnP approaches. Each network of the series, trained to transform a back-
projected data residual into an image residual, improves the dynamic range
from the previous iteration, with the sum of the output residual images across
iterations representing the final reconstruction.

Simulation results for RI imaging show that a series of only a few terms
achieves similar reconstruction quality to the state-of-the-art optimisation
algorithm SARA and PnP method AIRI, whilst being orders of magnitude faster.

Future work should investigate the robustness of the approach to the wide
variations of measurements settings characterising RI observation, and its
practical effectiveness on real large-scale high-dynamic range data. The
suggested superiority of the DNN architecture proposed over state-of-the-art
architectures should also be confirmed. The application of the proposed
methodology to other computational imaging modalities is naturally an open
question.

## References

  * [1] Yves Wiaux, Laurent Jacques, Gilles Puy, Anna MM Scaife, and Pierre Vandergheynst,  “Compressed sensing imaging techniques for radio interferometry,”  Monthly Notices of the Royal Astronomical Society, vol. 395, no. 3, pp. 1733–1742, 2009.
  * [2] Matthieu Terris, Arwa Dabbech, Chao Tang, and Yves Wiaux,  “Image reconstruction algorithms in radio interferometry: From handcrafted to learned regularization denoisers,”  Monthly Notices of the Royal Astronomical Society, 2022.
  * [3] Rizwan Ahmad, Charles A Bouman, Gregery T Buzzard, Stanley Chan, Sizhuo Liu, Edward T Reehorst, and Philip Schniter,  “Plug-and-play methods for magnetic resonance imaging: Using denoisers for image recovery,”  IEEE signal processing magazine, vol. 37, no. 1, pp. 105–116, 2020\.
  * [4] Matthew J Muckley, Bruno Riemenschneider, Alireza Radmanesh, Sunwoo Kim, Geunu Jeong, Jingyu Ko, Yohan Jun, Hyungseob Shin, Dosik Hwang, Mahmoud Mostapha, et al.,  “Results of the 2020 fastMRI challenge for machine learning MR image reconstruction,”  IEEE transactions on medical imaging, vol. 40, no. 9, pp. 2306–2317, 2021.
  * [5] Jingyun Liang, Jiezhang Cao, Guolei Sun, Kai Zhang, Luc Van Gool, and Radu Timofte,  “Swinir: Image restoration using swin transformer,”  in IEEE/CVF International Conference on Computer Vision, 2021, pp. 1833–1844.
  * [6] Claudio Gheller and Franco Vazza,  “Convolutional deep denoising autoencoders for radio astronomical images,”  Monthly Notices of the Royal Astronomical Society, 2021.
  * [7] Liam Connor, Katherine L Bouman, Vikram Ravi, and Gregg Hallinan,  “Deep radio-interferometric imaging with polish: Dsa-2000 and weak lensing,”  Monthly Notices of the Royal Astronomical Society, vol. 514, no. 2, pp. 2614–2626, 2022.
  * [8] Jonas Adler and Ozan Öktem,  “Learned primal-dual reconstruction,”  IEEE Transactions on Medical Imaging, vol. 37, no. 6, pp. 1322–1332, 2018.
  * [9] Sebastian Banert, Axel Ringh, Jonas Adler, Johan Karlsson, and Ozan Oktem,  “Data-driven nonsmooth optimization,”  SIAM Journal on Optimization, vol. 30, no. 1, pp. 102–131, 2020\.
  * [10] Singanallur V Venkatakrishnan, Charles A Bouman, and Brendt Wohlberg,  “Plug-and-play priors for model based reconstruction,”  in 2013 IEEE Global Conference on Signal and Information Processing. IEEE, 2013, pp. 945–948.
  * [11] Kai Zhang, Yawei Li, Wangmeng Zuo, Lei Zhang, Luc Van Gool, and Radu Timofte,  “Plug-and-play image restoration with deep denoiser prior,”  IEEE Transactions on Pattern Analysis and Machine Intelligence, 2021\.
  * [12] Jean-Christophe Pesquet, Audrey Repetti, Matthieu Terris, and Yves Wiaux,  “Learning maximally monotone operators for image recovery,”  SIAM Journal on Imaging Sciences, vol. 14, no. 3, pp. 1206–1237, 2021.
  * [13] Arwa Dabbech, Matthieu Terris, Adrian Jackson, Mpati Ramatsoku, Oleg M Smirnov, and Yves Wiaux,  “First AI for deep super-resolution wide-field imaging in radio astronomy: unveiling structure in ESO 137-006,”  Accepted for publication in The Astrophysical Journal Letters, 2022\.
  * [14] Heinz H Bauschke and Patrick L Combettes,  Convex analysis and monotone operator theory in Hilbert spaces,  Springer, 2017.
  * [15] Davis Gilton, Greg Ongie, and Rebecca Willett,  “Neumann networks for linear inverse problems in imaging,”  IEEE Transactions on Computational Imaging, vol. 6, pp. 328–343, 2019.
  * [16] Jiahui Yu, Yuchen Fan, Jianchao Yang, Ning Xu, Zhaowen Wang, Xinchao Wang, and Thomas Huang,  “Wide activation for efficient and accurate image super-resolution,”  arXiv preprint arXiv:1808.08718, 2018.
  * [17] Alexandru Onose, Rafael E Carrillo, Audrey Repetti, Jason D McEwen, Jean-Philippe Thiran, Jean-Christophe Pesquet, and Yves Wiaux,  “Scalable splitting algorithms for big-data interferometric imaging in the ska era,”  Monthly Notices of the Royal Astronomical Society, vol. 462, no. 4, pp. 4314–4335, 2016.
  * [18] Pierre-Antoine Thouvenin, Abdullah Abdulaziz, Arwa Dabbech, Audrey Repetti, and Yves Wiaux,  “Parallel faceted imaging in radio interferometry via proximal splitting (faceted hypersara): I. algorithm and simulations,”  Monthly Notices of the Royal Astronomical Society, 2022.
  * [19] AR Offringa, Benjamin McKinley, Natasha Hurley-Walker, FH Briggs, RB Wayth, DL Kaplan, ME Bell, Lu Feng, AR Neben, JD Hughes, et al.,  “WSCLEAN: an implementation of a fast, generic wide-field imager for radio astronomy,”  Monthly Notices of the Royal Astronomical Society, vol. 444, no. 1, pp. 606–619, 2014.


# Deep Slap Fingerprint Segmentation for Juveniles and Adults ††thanks: *
Authors contributed equally to this work. ** This material is based upon work
supported by the Center for Identification Technology Research and the
National Science Foundation (NSF) under Grant No.$650503$. Authors would also
like to thank the Chameleon cloud for providing the computational
infrastructure required for this work [1].

M. G. Sarwar Murshed *, Robert Kline*, Keivan Bahmani*, Faraz Hussain,
Stephanie Schuckers
Electrical and Computer Engineering
Clarkson University, Potsdam, NY
{murshem, kliner, bahmank, fhussain<EMAIL_ADDRESS>
###### Abstract

Many fingerprint recognition systems capture four fingerprints in one image.
In such systems, the fingerprint processing pipeline must first segment each
four-fingerprint slap into individual fingerprints. Note that most of the
current fingerprint segmentation algorithms have been designed and evaluated
using only adult fingerprint datasets. In this work, we have developed a
human-annotated in-house dataset of 15790 slaps of which 9084 are adult
samples and 6706 are samples drawn from children from ages 4 to 12.
Subsequently, the dataset is used to evaluate the matching performance of the
NFSEG, a slap fingerprint segmentation system developed by NIST, on slaps from
adults and juvenile subjects. Our results reveal the lower performance of
NFSEG on slaps from juvenile subjects. Finally, we utilized our novel dataset
to develop the Mask-RCNN based Clarkson Fingerprint Segmentation (CFSEG). Our
matching results using the Verifinger fingerprint matcher indicate that CFSEG
outperforms NFSEG for both adults and juvenile slaps. The CFSEG model is
publicly available at https://github.com/keivanB/Clarkson_Finger_Segment

###### Index Terms:

juvenile and adult fingerprints, deep slap Segmentation, Mask-RCNN

##  I Introduction

Due to their high accuracy and convenience, fingerprint-based recognition
systems are now in widespread use, e.g. in border crossings, cell-phone
authentication, and health care. Many fingerprint recognition systems use
high-end multi-finger scanners instead of single finger scanners in order to
achieve more accurate identification [2]. In such systems, the fingerprint
processing pipeline relies on a fingerprint segmentation model to localize
individual fingerprints [3] within each four-finger slap. As a result, slap
fingerprint segmentation models are an integral part of fingerprint
recognition systems as the loss of fingerprint ridge structure due to improper
segmentation can significantly degrade the overall matching performance [4].

To the best of our knowledge, all current slap fingerprint segmentation models
have been developed using adult datasets [3, 5, 6]. However, since aging is
known to affect the size and ridge-valley structure of the fingerprints, lower
the quality and diminish the matching performance [7, 8], careful
consideration and modeling are required to accurately segment and match slaps
from juvenile subjects. However, the lack of publicly available juvenile
datasets, as well as privacy concerns, have led to a dearth of work on this
topic.

In this work, we have developed a human-annotated in-house dataset of 15790
slaps (Children: 6706, Adult: 9084). Our new dataset allows us to accurately
evaluate the performance of the NFSEG model on both adult and juvenile
subjects. NFSEG is one of the most popular fingerprint segmentation system
published by NIST [5]. This segmentation system is used as a benchmark
algorithm to evaluate the performance of a newly developed slap segmentation
algorithm. Additionally, we utilized our human-annotated dataset to develop
Clarkson Fingerprint SEGmentation (CFSEG) – a novel deep learning-based
fingerprint segmentation model designed to accurately segment both adult and
juvenile slap fingerprints. Our test results using the Verifinger fingerprint
matcher (version 10, compliant with NIST MINEX [9]), show that CFSEG
outperforms NFSEG in both adult and juvenile cohorts of our dataset.

##  II Collection and Annotation of Slap Fingerprints

We combined several in-house adult and juvenile datasets to create a balanced
and representative dataset for training and evaluating CFSEG. Our dataset
contains a total of 15790 slap fingerprints (adults: 9084, juvenile: 6706)
from 203 adults and 242 juvenile subjects. All slaps are captured at 500 ppi
using the FBI certified Crossmatch L Scan Guardian (9000251) fingerprint
scanner [10].

Initially, all slaps are segmented using the NFSEG model. Subsequently, we
developed a Graphic User Interface (GUI) based on the open-source `labelImg`
[11] platform for fingerprint annotation. Our GUI allows us to utilize NFSEG
segmentations as a starting point for each capture. This process significantly
reduced our annotation time to a practical time frame. Additionally, this
pipeline allows us to utilize the estimated angle of rotation from the NFSEG
and rotate each slap to the upright position. In the first stage, human
annotators evaluated each slap for the correct angle of rotation. In the
second stage, annotators cycle through slaps and correct the fingerprints that
are incorrectly segmented by NFSEG while the rest of the bounding boxes remain
untouched. This process resulted in a cleaned and annotated dataset with a
total of 52674 localized fingerprints (adult: 30304, juvenile: 22370).
Finally, we used a rigorous two stage 10-fold cross-validation process based
on the unique identities in the dataset to derive our training, validation,
and test sets. At each fold, $80\%$ of the adult and $80\%$ of the juvenile
identities were used for training while the remaining $20\%$ of identities of
both types were used for validation ($10\%$) and testing ($10\%$).

##  III Clarkson Fingerprint Segmentation (CFSEG)

Our CFSEG model for slap segmentation is based on the Mask-RCNN architecture
[12]. In this section, we discuss the overall Mask R-CNN model, the loss
functions, and the training scheme used for training the CFSEG.

###  III-A Architecture of Mask R-CNN

sh Mask-RCNN is a simple and flexible two-stage deep architecture developed to
perform semantic segmentation, object localization, and object instance
segmentation of natural images [12]. As a result, we adopt Mask-RCNN
architecture for CFSEG.

The first stage of this architecture contains the Region Proposal Network
(RPN) which proposes candidate object bounding boxes. It consists of two
networks, a Convolution Neural Network (CNN) and a Region Proposal Network
(RPN). The CNN is the backbone of the Mask R-CNN architecture, and it is
responsible for extracting feature maps from the input images. Any CNN model
designed for image classification tasks (such as ResNet, MobileNet, or VGG)
could be used as the backbone network [13]. Previous research has shown that
the ResNet-FPN backbone provides better performance in terms of accuracy and
speed on object detection tasks [12]. Therefore, we have used ResNet-101 [14]
as the backbone network in our experiments. On the other hand, the RPN is
responsible for generating a set of ‘Region Proposals’, which signify regions
in the feature map that have a high probability of containing objects. In the
second stage, a small, Fully Connected Neural Network (FCNN) takes the
proposed regions from the first stage and predicts bounding boxes and object
classes for each of them. Note that the proposed regions generated by the RPN
can be of different sizes. However, fully connected layers in the networks
only take a fixed size vector to make predictions. The ROIAlign [12], which is
a modified version of Max-Pooling, is used to fix the size of these proposed
regions. ROIAlign also fixes the misalignment which preserves the exact
spatial location of an object and helps to improve overall accuracy.

A fully convolutions network is added (parallel with the existing branch for
classification and bounding box regression) which is responsible for
predicting segmentation masks in a pixel-to-pixel manner on each Region Of
Interest (ROI).

###  III-B Loss function

For training RPNs, a class label indicating whether it is an object or not is
assigned to each anchor. Anchor boxes are reference boxes placed at different
positions in the input image. A positive label is assigned to an anchor if
that anchor has an Intersection over Union (IoU) greater than 0.7 with any
ground truth objects. A negative label is assigned to an anchor if that anchor
has an IoU less than 0.3 with any ground truth objects. The anchors that have
an IoU value between 0.3 and 0.7 are considered neutral and excluded from the
training set. Shift and resizing operations are performed before training the
RPN which helps make the anchor cover the ground truth object completely. With
this definition, a multi-task loss function is used to train the Mask R-CNN.
This loss function is divided into three parts which combine the loss of
classification, localization, and segmentation as follows:

$L=L_{cls}+L_{box}+L_{mask}$ (1)

The class, bounding box, and mask losses are represented by $L_{cls}$,
$L_{box}$ and $L_{mask}$ respectively. $L_{cls}$ is the log loss function over
two classes carried out as a binary classification by predicting an object
being a target object or not. The smooth L1 loss is used for bounding box loss
$L_{box}$. Mask loss $L_{mask}$ is the average binary cross-entropy loss, and
it is calculated for each class separately. These calculations prevent
competition among the classes when generating masks. The details of different
loss functions are as follows:

$L(\\{p_{i}\\},\\{t_{i}\\})=\frac{1}{N_{cls}}\sum_{i}L_{cls}(p_{i},p_{i}^{*})+\frac{\lambda}{N_{box}}\\\
\sum_{i}p_{i}^{*}L_{box}^{smooth}(t_{i}-t_{i}^{*})+(-1*\frac{1}{m^{2}})\sum_{1\leq
i,j\leq m}[y_{ij}\log y_{ij}^{k}\\\ +(1-y_{ij})\log(1-y_{ij}^{k})]$ (2)

Here, $i$ represents the index number of an anchor and $p_{i}$ is the
predicted probability of an anchor being an object. The ground-truth binary
label is represented by $p_{i}^{*}$. $t_{i}$ and $t_{i}^{*}$ represent the
vector with four coordinates of the predicted and ground-truth bounding box
respectively. $m$ is the size of the mask output; $y_{ij}$ is the label of a
cell (i,j) in the true mask for the region of size $m\times m$; $y^{k}_{ij}$
is the k-th mask predicted value of the same cell in the mask learned for the
ground-truth class k.

###  III-C Training

An image-centric training approached [15] is utilized in this work and only
positive anchors are used to determine the loss of the network. Our training
approach is as follows:

  * •

N ROIs are generated from each image with the ratio of positive to negative
ROIs of 1:3, where N is selected as 64 and 512 for backbone and FPN stages
respectively.

  * •

When training RPN, the anchors aspect ratio is 3, and the span is 5 scales.
The RPN is trained individually and does not share features with other parts
of the Mask-RCNN network.

  * •

The Stochastic Gradient Descent (SGD) with a learning rate of 0.001, momentum
of 0.9, and decay of 0.0001 are used to train the network.

During inference, the region proposal number increases to 300 and 1000 for the
backbone and FPN stages, respectively. This increment helps the network not
miss any potential regions during inference time. The bounding box prediction
branch considers each proposal proposed by the backbone and FPN network. Then
the non-maximum suppression is performed to remove the low-scoring ROIs. To
keep computational overhead to a minimum during training and inference, the
mask branch only runs over the 50 highest scoring ROIs.

Figure 1: Example for calculating positive and negative errors between the
predicted and ground truth bounding box.

##  IV Evaluation and Results

In this work, we utilized both the Mean Absolute Error (MAE) of the predicted
bounding boxes and fingerprint matching performance as our evaluation metrics.
In this section, we describe our evaluation process and results.

###  IV-A Mean Absolute Error (MAE)

The MAE measures the distance between the detected bounding box and the
annotated ground-truth bounding box in terms of pixels. Slap fingerprint
segmentation models have to find the balance between over-extending the
bounding boxes (over-segmentation) vs. excessively reducing the predicted
bounding box (under-segmentation). Over-segmentation can lead to the ridge-
valley structure of other fingerprints leaking into the bounding box of other
fingerprints. This extra noise can potentially degrade the matching
performance. On the other hand, under-segmentation can cause the loss of the
valuable parts of the fingerprints.

Previous work has identified that under-segmentation beyond 32 pixels from the
sides and 64 pixels from the top or bottom of the fingerprint can negatively
affect the matching performance [4]. As a result, we calculate and report MAE
for _each side_ of the bounding box individually, i.e. the error is calculated
by measuring how far those four sides are from the four sides of the ground
truth bounding box. In Figure 1, the red dotted rectangle represents the
ground-truth bounding box generated by a human-annotator and the blue
rectangle represents the detected bounding box by a fingerprint detection
model. If any side of a detected bounding box captures more information than
the side of the ground-truth bounding box, we consider it as a positive error.
On the other hand, if any side of a detected bounding box captures less
information than the ground-truth bounding box, we consider it as a negative
error. Finally, the MAE for each side is calculated using the equation 3,
Where, $n$ is the total number of fingerprints in the test dataset. $X$
represents any side such as left, right, top, bottom, of the bounding box.

$MAE=\frac{\sum_{i=0}^{n}abs(X\>error_{i})}{Total\\_fingerprints\\_in\\_the\\_dataset}$
(3)

Figure 2: The histogram of the bottom error (pixels) from predicted bounding
boxes by the NFSEG and CFSEG models for adult (Left) and Juvenile (Right)
subjects.

Table I shows the mean and standard deviation of the MAE for NFSEG and CFSEG
models evaluated using our dataset. We can observe that, compared to NSFEG,
CFSEG has lower MAE in almost all directions and produces more precise
bounding boxes for both adults and juvenile subjects. Additionally, our
results confirm that compared to adults, NFSEG suffers from higher MAE in all
directions in the juvenile cohort of our dataset. Furthermore, we observe that
in both adult and juvenile subjects, NFSEG has difficulty in accurately
localizing the lower boundary of fingerprints. Figure 2 illustrates the
histogram of the bottom errors for adult and juvenile subjects. We can observe
that while CFSEG maintains the performance in Juvenile subjects, NFSEG is
particularly susceptible to this type of error in juvenile subjects. Our
failure analysis on such cases revealed that the high error observed in this
direction is because of the over-segmentation at the distal interphalangeal
joint. Figure 3 depicts an example of this problem, where the red bounding box
indicates the ground truth and the blue bounding box presents the bounding box
predicted by NFSEG. Additionally, the long tail of the error histogram of the
NFSEG model for the juvenile subjects suggests that NFSEG over-segmented many
juvenile samples below the Minimum Tolerance Limit (MTL) of -64 pixels
suggested by Slapseg-II [4]. This can potentially reduce the matching
performance of the system.

TABLE I: Mean Absolute Error (MAE) [Mean (Std.)] for NFSEG and CFSEG Age Group | Adults | Children
---|---|---
Model |  NFSEG |  CFSEG |  NFSEG |  CFSEG
Left |  07.13(28.59) |  08.21(14.66) |  12.09(38.19) |  08.36(16.25)
Top |  13.18(44.70) |  13.05(20.26) |  28.47(92.51) |  13.77(21.13)
Right |  09.46(31.52) |  07.73(13.78) |  13.60(34.55) |  07.21(14.16)
Bottom |  35.23(84.73) |  16.60(24.41) |  102.87(159.00) |  15.00(22.92)
Figure 3: Example of over-segmentation error in NFSEG. The red rectangle shows
the ground-truth bounding box while the blue rectangle shows the detected
bounding box by the NFSEG.

###  IV-B Fingerprint Matching

In the second stage of the evaluation process, we utilized the Verifinger
fingerprint matcher (version 10, compliant with NIST MINEX [9]) to compare the
matching accuracy using the fingerprints segmented with NFSEG and CFSEG. This
helps us better analyze the effects of segmentation accuracy on the overall
performance of the fingerprint recognition system. In order to evaluate the
matching accuracy, we used the annotated ground truth, NFSEG, and CFSEG
bounding boxes. For every fingerprint in the dataset, we evaluated all the
mated comparisons for our genuine distribution (172348), while randomly
selecting 20 non-mated fingerprints to construct an imposter distribution
(441322 imposter comparisons). In total, we performed $~{}600,000$ comparisons
using all 10 fingers.

Our results indicate that for both adults and juvenile cohorts the
fingerprints segmented with CFSEG provide higher accuracy across operation
points. Additionally, we observe lower overall matching performance in
juvenile subjects compared to adults. Finally, Table II depicts the True
Positive Rate (TPR) at False Positive Rate (FPR) of 0.1% for fingerprints
segmented using NFSEG, CFSEG, and human-annotators.

TABLE II: TPR [Mean (Std.)] at FPR of 0.001 for NFSEG, CFSEG, and Ground Truth. Model | Adults | Children
---|---|---
NFSEG | 0.9972 (0.0027) | 0.9675 (0.0135)
CFSEG | 0.9977 (0.0026) | 0.9687 (0.0135)
Ground-Truth | 0.9991 (0.0011) | 0.9716 (0.0137)

##  V Conclusion and Limitations

The experimental results indicate that our novel CFSEG model provides more
precise bounding box predictions with respect to NFSEG for both adult and
juvenile subjects. However, CFSEG currently depends on the estimated rotation
angle of the slap acquired from NFSEG. This will be addressed in future work
by the introduction of a regression sub-network for angle prediction and
artificially augmenting the training dataset with rotated slaps. Additionally,
our fingerprint matching experiment is conducted using the human-annotated
class labels (i.e. index, middle, or ring fingers) for all algorithms to
produce a meaningful comparison between the predicted bounding boxes. Future
work can focus on designing a class prediction sub-network to integrate this
task into the CFSEG and further improve the overall performance. Finally, the
trained CFSEG model is publicly available to other researchers at
https://github.com/keivanB/Clarkson_Finger_Segment.

## References

  * [1] K. Keahey, J. Anderson, Z. Zhen, P. Riteau, P. Ruth, D. Stanzione, M. Cevik, J. Colleran, H. S. Gunawi, and C. Hammock, “Lessons learned from the chameleon testbed,” in _Proceedings of the 2020 USENIX Annual Technical Conference (USENIX ATC ’20)_. USENIX Association, 2020, pp. 219–233.
  * [2] D. Maltoni, D. Maio, A. K. Jain, and S. Prabhakar, _Handbook of fingerprint recognition_. Springer Science & Business Media, 2009.
  * [3] P. Feulner, “Slap fingerprint segmentation evaluationII (slapsegII),” https://www.nist.gov/itl/iad/image-group/slap-fingerprint-segmentation-evaluation-ii-slapsegii, Feb. 2010.
  * [4] C. I. Watson, “SlapSegII Analysis: Matching Segmented Fingerprint Images,” pp. 1–53, Nov. 2010.
  * [5] Kenneth Ko, Wayne J. Salamon, “NIST Biometric Image Software (NBIS),” https://www.nist.gov/services-resources/software/nist-biometric-image-software-nbis, February 2010.
  * [6] P. Gupta and P. Gupta, “A slap fingerprint based verification system invariant to halo and sweat artifacts,” _Applied Mathematical Modelling_ , vol. 54, pp. 413–428, 2018.
  * [7] J. Galbally, R. Haraksim, and L. Beslay, “A study of age and ageing in fingerprint biometrics,” _IEEE Transactions on Information Forensics and Security_ , vol. 14, pp. 1351–1365, 2018.
  * [8] R. Haraksim, J. Galbally, and L. Beslay, “Fingerprint growth model for mitigating the ageing effect on children’s fingerprints matching,” _Pattern Recognition_ , vol. 88, pp. 614–628, 2019.
  * [9] C. Watson, G. Fiumara, E. Tabassi, S. Cheng, P. Flanagan, and W. Salamon, “Fingerprint vendor technology evaluation, nist interagency/internal report 8034: 2015,” 2014.
  * [10] “Cross Match L SCAN Patrol and Patrol ID fingerprint scanners,” https://neurotechnology.com/fingerprint-scanner-cross-match-l-scan-patrol.html.
  * [11] Tzutalin, “Labelimg,” https://github.com/tzutalin/labelImg, 2015.
  * [12] K. He, G. Gkioxari, P. Dollár, and R. Girshick, “Mask r-cnn,” in _Proceedings of the IEEE international conference on computer vision_ , 2017, pp. 2961–2969.
  * [13] S. Gonzalez, C. Arellano, and J. Tapia Farias, “Deepblueberry: Quantification of blueberries in the wild using instance segmentation,” _IEEE Access_ , vol. PP, pp. 1–1, 08 2019.
  * [14] K. He, X. Zhang, S. Ren, and J. Sun, “Deep Residual Learning for Image Recognition,” _arXiv:1512.03385 [cs]_ , Dec. 2015, arXiv: 1512.03385.
  * [15] S. Ren, K. He, R. Girshick, and J. Sun, “Faster r-cnn: Towards real-time object detection with region proposal networks,” in _Advances in Neural Information Processing Systems_ , vol. 28, 2015.