# Superfluid properties of a honeycomb dipolar supersolid

Albert Gallemí  Institut für Theoretische Physik, Leibniz Universität
Hannover, 30167 Hannover, Germany  Luis Santos  Institut für Theoretische
Physik, Leibniz Universität Hannover, 30167 Hannover, Germany

###### Abstract

Recent breakthrough experiments on dipolar condensates have reported the
creation of supersolids, including two-dimensional arrays of quantum droplets.
Droplet arrays are, however, not the only possible non-trivial density
arrangement resulting from the interplay of mean-field instability and quantum
stabilization. Several other possible density patterns may occur in trapped
condensates at higher densities, including the so-called honeycomb supersolid,
a phase that exists, as it is also the case of a triangular droplet
supersolid, in the thermodynamic limit. We show that compared to droplet
supersolids, honeycomb supersolids have a much-enhanced superfluid fraction
while keeping a large density contrast, and constitute in this sense a much
better dipolar supersolid. However, in contrast to droplet supersolids,
quantized vortices cannot be created in a honeycomb supersolid without driving
a transition into a so-called labyrinthic phase. We show that the reduced
moment of inertia, and with it the superfluid fraction, can be however
reliably probed by studying the dynamics following a scissors-like
perturbation.

††preprint: APS/123-QED

##  I Introduction

The recent realization of the supersolid phase in ultracold dipolar quantum
gases Tanzi  _et al._ (2019a); Böttcher  _et al._ (2019); Chomaz  _et al._
(2019) opens intriguing questions about the superfluidity in these systems
Leggett (1970). The interplay between mean-field interactions, quantum
fluctuations, and external confinement results in a rich landscape of possible
density patterns, depending on the atom-atom interactions, the density, and
the external trap Hertkorn  _et al._ (2021); Zhang  _et al._ (2019, 2021);
Poli  _et al._ (2021). Of these possible patterns, droplet arrays have been
the most intensively explored phase. Dipolar droplet supersolids (D-SS) have
been realized in recent years both in one- Böttcher  _et al._ (2019); Chomaz
_et al._ (2019); Tanzi  _et al._ (2019a) and two-dimensional Norcia  _et al._
(2021); Bland  _et al._ (2022) arrangements. In two dimensions, quantum
droplets arrange in a triangular crystalline-like order, surrounded by a low-
density region that ensures phase coherence amongst the different droplets.

Whereas most of the predicted density patterns demand an external confinement,
in addition to the D-SS, only another modulated phase may occur in homogeneous
space Zhang  _et al._ (2019), namely the so-called honeycomb supersolid
(H-SS). This pattern consists of a condensate containing hexagonally-arranged
density holes, constituting the complementary pattern to the triangular
droplet array. The observation of the H-SS demands significantly larger
densities than the D-SS. Although this is still an important limitation in
current experiments based on magnetic atoms, the problem could be circumvented
with a tight confinement Zhang  _et al._ (2021); Poli  _et al._ (2021), and
could be much less relevant in future experiments with condensates of polar
molecules Schmidt  _et al._ (2022).

Although recent experiments have revealed the expected presence of two
different Goldstone modes associated with the phase coherence and the
crystalline order Tanzi  _et al._ (2019b); Guo  _et al._ (2019); Natale  _et
al._ (2019), a direct clear proof of superfluid effects in dipolar supersolids
is still lacking. Whereas in standard condensates the study of the scissors
mode frequency provides a clear information about the reduction of the moment
of inertia, and with it about superfluidity Guéry-Odelin and Stringari (1999);
Maragò  _et al._ (2000), recent studies have revealed that a more subtle
analysis of the scissors response is necessary in D-SSs Tanzi  _et al._
(2021); Norcia  _et al._ (2022); Roccuzzo  _et al._ (2022). An alternative
would be provided by the observation of quantized vortices. Recent studies
have shown that vortices may be robustly created in D-SSs, where they occupy
the inter-droplet interstitials Roccuzzo  _et al._ (2020); Gallemí  _et al._
(2020); Šindik  _et al._ (2022). Vortices have not yet been observed in a
D-SS, but the recent experimental creation of vortices in an unmodulated
dipolar condensate Klaus  _et al._ (2022) opens optimistic perspectives in
this direction.

These recent works on dipolar density patterns and on superfluid effects in
D-SSs lead to intriguing questions about the superfluidity of other density
patterns, and very especially of the H-SS, and about how the superfluidity of
this phase may be experimentally revealed. This paper is devoted to these
questions. We show that the H-SS has for a large contrast of the density
modulation a much larger superfluid fraction than the D-SS. However, contrary
to the D-SS, the H-SS is not robust when vortices nucleate, rather
transitioning into a labyrinthic phase. Fortunately, the large superfluid
fraction of the H-SS could be revealed under realistic conditions by a careful
analysis of the scissors response.

The paper is organized as follows. In Sec. II, we introduce the formalism
employed, and review the possible ground-state density patterns. Section III
is devoted to the study of the density patterns under rotation. In Sec. IV, we
analyze the moment of inertia of the different patterns, paying special
attention to the H-SS. Section V discusses in detail how the response to a
scissors perturbation may be employed to provide a good estimation of the
reduction of the moment of inertia in the H-SS, revealing its large superfluid
fraction. Our conclusions are summarized in Sec. VI.

Figure 1: (Top left) Ground-state phase diagram as a function of the
scattering length and the atom number, for 162 Dy atoms confined in a harmonic
trap with frequencies $\omega_{\rm x,y,z}=2\pi\times(200,200,400)$. The panels
(a)-(e) show the density pattern in the $xy$ plane for $N=2\times 10^{5}$ and
$a/a_{0}=70,72,78,82$ and $84$, which correspond to the droplet, labyrinthic,
stripe, H-SS, and unmodulated patterns, respectively.

##  II Ground-state patterns

In the following, we consider a condensate of $N$ bosons with mass $m$ and
magnetic dipole moment $\mu$ oriented along the $z$ direction. A similar
physics is expected for the case of electric dipoles, as it is in particular
the case of polar molecules Schmidt  _et al._ (2022). The physics of a dipolar
condensate is given by the interplay between contact interactions,
characterized by the $s$-wave scattering length $a$, dipolar interactions, and
(despite of its weakly-interacting character) quantum fluctuations. The latter
provides the stabilization mechanism against mean-field collapse Petrov
(2015). The physics of a quantum-stabilized dipolar condensate is well
described by the so-called extended Gross-Pitaevskii equation Wächtler and
Santos (2016),

$\displaystyle i\hbar\dot{\Psi}({\bf r},t)$ $\displaystyle=$
$\displaystyle\left[-\frac{\hbar^{2}}{2m}\nabla^{2}+V_{\rm trap}({\bf
r})+g|\Psi({\bf r},t)|^{2}\right\delimiter 00$ (1) $\displaystyle+$
$\displaystyle\int V_{\rm dd}({\bf r}-{\bf r}^{\prime})|\Psi({\bf
r}^{\prime},t)|^{2}d{\bf r}^{\prime}$ $\displaystyle+$
$\displaystyle\left\delimiter 00\gamma|\Psi({\bf r},t)|^{3}\right]\Psi({\bf
r},t)\,.$

where $\Psi({\bf r},t)$ is the condensate wave function. In Eq. (1) $V_{\rm
trap}({\bf r})$ is the external trapping of the system, $g=4\pi\hbar^{2}a/m$
is the coupling constant, $V_{\rm
dd}=\frac{\mu_{0}\mu^{2}}{4\pi}\frac{1-3\cos^{2}\theta}{|{\bf r}-{\bf
r}^{\prime}|^{3}}$ is the dipole-dipole interaction potential, with $\mu_{0}$
the vacuum magnetic permeability, and $\theta$ the angle between ${\bf r}-{\bf
r}^{\prime}$ and the polarization direction. The last term of the equation
accounts for the effects of quantum fluctuations, the so-called Lee-Huang-Yang
correction Lima and Pelster (2012), with

$\\!\\!\\!\gamma\\!=\\!\frac{16}{3\sqrt{\pi}}ga^{3/2}\mbox{Re}\left[\int_{0}^{\pi}\\!\\!\\!d\phi\sin\phi[1\\!+\\!\varepsilon_{\rm
dd}(3\cos^{2}\phi-1)]^{\frac{5}{2}}\\!\right]\\!,$ (2)

where $\varepsilon_{\rm dd}=g_{\rm dd}/g$ is the ratio between the dipolar
strength $g_{\rm dd}=\mu_{0}\mu^{2}/3$ and the contact strength $g$. The
ground state is obtained by setting in Eq. (1) $\Psi({\bf r},t)=e^{-i\mu
t/\hbar}\Psi({\bf r})$, with $\mu$ the chemical potential. Normalizing $\int
d^{3}r|\Psi({\bf r})|^{2}=N$, we obtain $\Psi({\bf r})$ using standard
imaginary-time evolution techniques.

The interplay between dipole-induced mean-field instability, quantum
stabilization, and external confinement results in a variety of possible
ground state phases. If the condensate is mean-field stable, the density
profile remains unmodulated, presenting (for sufficiently large interactions)
a standard Thomas-Fermi parabolic form. We denote this phase in the following
as the unmodulated phase. For a sufficiently low $a$, and small-enough $N$,
the condensate is mean-field unstable, breaking into quantum droplets, which
may remain coherently linked (supersolid regime) or eventually disconnected
(independent-droplet regime). These phases have been experimentally observed
both in quasi-one- (Tanzi  _et al._ (2019a); Böttcher  _et al._ (2019); Chomaz
_et al._ (2019)) and quasi-2D geometries (Norcia  _et al._ (2021); Bland  _et
al._ (2022)). In the latter case, in the thermodynamic limit, the droplet
pattern is expected to acquire a triangular geometry Zhang  _et al._ (2019,
2021); Bland  _et al._ (2022). Recent works Hertkorn  _et al._ (2021); Zhang
_et al._ (2021) have shown, however, that the ground state phase diagram may
be significantly richer for larger atoms numbers (for fixed external
confinement). Particularly interesting is the so-called honeycomb phase, which
as for the case of the triangular droplet array, extends all the way to the
thermodynamic limit Zhang  _et al._ (2019, 2021). This phase is characterized
by the formation of empty regions, which are arranged in an hexagonal pattern.
In the presence of harmonic confinement, a rich variety of additional patterns
may occur, including the so-called labyrinthic, stripe, and pumpkin phases
Hertkorn  _et al._ (2021).

Figure 1 (top left panel) illustrates these possible ground-state phases for
the particular case of a condensate of 162Dy atoms (which possess a dipole
moment $\mu=10\mu_{B}$, $\mu_{B}$ being the Bohr magneton). We consider
throughout this paper an external harmonic confinement

$V_{\rm trap}=\frac{1}{2}m(\omega_{\rm x}x^{2}+\omega_{\rm y}y^{2}+\omega_{\rm
z}z^{2})\,,$ (3)

with frequencies $\omega_{\rm x,y}=\omega=2\pi\times 200\,\mathrm{Hz}$ and
$\omega_{\rm z}=2\pi\times 400\,\mathrm{Hz}$. For $N$ up to $3\times 10^{5}$
atoms and $70<a/a_{0}<84$, with $a_{0}$ the Bohr radius, the system presents
different density patterns: unmodulated (panel (a)), H-SS (panel (b)), stripe
phase (panel (c)) (formed by a pattern of concentric rings), labyrinthic
(panel (d)) and droplet array (panel (e)). The so-called pumpkin pattern only
appears in this case for much larger atom numbers Hertkorn  _et al._ (2021).

Figure 2: Ground-state phase diagram of the system under rotation for
$N=2\times 10^{5}$ as a function of the rotation frequency $\Omega$, for the
same parameters as in Fig. 1. Above a certain critical value of $\Omega$
(represented by the line in the plot and the red squares), a vortex nucleates.
Note that for the stripe and the H-SS phases this results in the breaking of
the pattern into a labyrinthic arrangement (shaded region).

##  III Robustness of the density pattern under rotation

The creation of quantized vortices constitute a possible way to directly probe
superfluidity. We consider at this point that the condensate is rotated around
$z$ with a rotation frequency $\Omega$. We obtain the ground state by
evaluating the eGPE in the rotating frame, minimizing $H-\Omega\hat{L}_{z}$,
with $\hat{L}_{z}$ the angular momentum operator along $z$. In the unmodulated
phase, beyond a critical $\Omega>\Omega_{c}$, the solution with a vortex at
the trap center becomes energetically favorable. In the presence of a density
modulation, rotation does not only eventually result in the presence of
vortices in the ground state, but may compromise as well the density pattern
of the non-rotating ground state. In the D-SS phase, the ground state under
rotation remains a droplet array (although the number of droplets may change),
presenting, beyond a critical rotation frequency, vortices in the inter-
droplet interstitials Roccuzzo  _et al._ (2020); Gallemí  _et al._ (2020);
Šindik  _et al._ (2022). Interestingly, the situation is different for other
density patterns.

We focus our attention on the particular case of $N=2\times 10^{5}$ atoms,
which for the range of scattering lengths considered comprises all possible
density patterns. Figure 2 depicts the ground state phase diagram as a
function of $\Omega$ and $a$. As mentioned above, the droplet array is robust
under rotation. In contrast, both the stripe phase and the H-SS are not.
Increasing the rotation frequency results in a transition of the ground-state
from a non-rotating stripe or H-SS pattern into a labyrinthic phase with a
vortex. There is hence no ground state with a vortex with a stripe or H-SS
pattern.

This has important consequences for the actual nucleation of vortices in
experiments, which is typically realized by stirring the condensate by means
of a slightly anisotropic rotating confinement on the $xy$ plane. Although the
solution without vortex is not the ground state, it remains metastable for
$\Omega>\Omega_{c}$ because moving a vortex to the center of the trap demands
overcoming a potential barrier Sinha and Castin (2001). Such a barrier is
eventually circumvented for a sufficiently large rotational frequency,
typically much larger $\Omega_{c}$, at which the quadrupole surface mode is
destabilized. This remains true also for a D-SS, although the critical
rotation frequency for dynamical vortex nucleation may be significantly
reduced compared to the unmodulated case Gallemí  _et al._ (2020).

In the D-SS phase, since the ground state under rotation remains a droplet
array, the surface deformation is followed by the penetration of a vortex or
vortices through the interstitials while keeping the droplet structure
Roccuzzo  _et al._ (2020); Gallemí  _et al._ (2020). The situation is
radically different in the stripe and H-SS patterns, due to the absence of a
ground state with a vortex in those cases. When the condensate is stirred at a
rotation frequency high-enough to destabilize the quadrupole surface mode, the
stripe and the honeycomb patterns are destroyed while vortices nucleate.
Therefore, as expected from the ground-state properties discussed above, a
labyrinthic pattern is formed with vortices in its interstitials. This is
illustrated in Fig. 3. Starting with the non-rotating H-SS ground state for
$a=82a_{0}$ (left panel), we increase adiabatically $\Omega$. The eventual
surface deformation for rotational frequencies above the quadrupole frequency
leads to the creation of a vortex or vortices, but, prior to that, also to the
destruction of the honeycomb pattern (right panel).

As discussed below, the H-SS constitutes a much better supersolid than the
droplet array. However, our results show that the lack of robustness of the
density pattern under rotation is expected to prevent under typical
experimental conditions the direct probe of superfluidity of that phase by
means of the creation of quantized vortices. Fortunately, as shown in the
following, the study of the moment of inertia constitutes a feasible
alternative to directly probe the superfluidity of the H-SS.

Figure 3: (left panel) Initial condensate (without rotation) in the H-SS phase
at $t=0$ for $a=82a_{0}$ (and the same parameters as Fig. 4). When the trap is
rotated with a trap anisotropy $\delta=0.025$ (see text) increasing
adiabatically with a linear ramp of $0.03$ ms the rotation frequency up to
$\Omega=0.9\omega_{x}$ (over the threshold frequency for the dynamical
nucleation of vortices), the honeycomb pattern is destroyed into a labyrinthic
arrangement, shown in the right panel at $t=43$ ms.

##  IV Moment of inertia

The moment of inertia may be theoretically evaluated from the response to an
infinitesimally slow rotation, as $I_{\rm exact}=\lim_{\Omega\to
0}\frac{d\langle L_{z}\rangle}{d\Omega}$, where $\langle\hat{L}_{z}\rangle$ is
the expected value of the angular momentum operator in the ground state
calculated in the rotating frame. The reduction of the moment of inertia with
respect to its classical value in a cylindrically-symmetric quantum gas,
$I_{\rm C}=m\int n(\mathbf{r})(x^{2}+y^{2})d{\bf r}$, is a clear signature of
a finite superfluid fraction, which together with the density modulation,
becomes a direct probe for supersolidity Leggett (1970); Roccuzzo  _et al._
(2020, 2022).

Figure 4(a) shows, for $N=2\times 10^{5}$, $I_{\rm exact}/I_{\rm C}$ (black
dots) as a function of $a$ . The moment of inertia provides crucial
information about the superfluid properties, although the lack of cylindrical
symmetry may partially mask the true superfluid nature in some configurations.
Due to the (quasi-)cylindrical symmetry of the unmodulated and stripe phases,
the moment of inertia of these configurations turns out to be very close to
zero. In both the labyrinth and the droplet patterns, the moment of inertia
increases dramatically and approaches the classical rigid value, indicating
the rigid response of the solid part. Most relevantly, the H-SS presents a
very low moment of inertia, showing that, compared to the droplet array, it
has typically a much larger superfluid fraction, while presenting a large
density contrast.

##  V Scissors-like perturbation

The reduction of the moment of inertia associated to superfluidity may be
probed by monitoring the response of the system against a scissors-like
perturbation. In the unmodulated phase, the frequency of the scissors mode may
be directly linked to the superfluid fraction Guéry-Odelin and Stringari
(1999); Maragò  _et al._ (2000). As shown by recent experiments on dipolar
condensates, the situation is more subtle in the case of a D-SS Tanzi  _et
al._ (2021); Norcia  _et al._ (2022); Roccuzzo  _et al._ (2022), since
determining the dominant frequency of the scissors mode does not generally
allow for a clear proof of global superfluidity. However, a careful study of
the overall spectrum of the system response against the scissors-like
perturbation may provide an experimental estimation of the reduction of the
moment of inertia under proper conditions Roccuzzo  _et al._ (2022). A similar
procedure may be employed for the study of superfluidity in the H-SS.

Figure 4: (a) Moment of inertia $I_{\rm exact}/I_{\rm C}$ (black dots) for
$N=2\times 10^{5}$ atoms, as a function of the scattering lengths (all other
parameters as in Fig. 4). Panels (b)-(d) show the Fourier transform ${\cal
F}(\omega)$ of the signal $\langle xy\rangle(t)$ of the scissors mode (see
text), measured after a sudden stop of a rotation of $1$ Hz for
$a/a_{0}=81.5$, $82$, and $82.5$, respectively. The estimation of the moment
of inertia using Eq. (4) for those cases is indicated with red stars in panel
(a).

We consider in the following a slightly deformed trap on the $xy$ plane, with
$\omega_{x,y}=\omega(1\pm\delta)$, with $\delta=0.025$. Two possible
experimental procedures may be employed for the study of the moment of
inertia. The scissors-like perturbation may be induced either by starting with
the ground state and tilting the trap on the $xy$ plane, or by rotating the
condensate around $z$ and then suddenly stopping the rotation. In both cases,
we may monitor the scissors-mode response $\langle xy\rangle(t)$ as a function
of time and obtain the Fourier transform within a given time window $0<t<T$:
${\cal F}(\omega)=\int_{0}^{T}dt\langle xy\rangle(t)e^{i\omega t}$. Using sum
rules, it is possible to show that the moment of inertia may be estimated as
Roccuzzo  _et al._ (2022):

$I_{\rm p}=m(\omega_{y}^{2}-\omega_{x}^{2})\langle
x^{2}-y^{2}\rangle\frac{\int d\omega\mathcal{F}(\omega)\omega^{p-2}}{\int
d\omega\mathcal{F}(\omega)\omega^{p}},$ (4)

with $p=0$ ($p=1$) for the tilting (rotation) protocol. Recent studies of the
D-SS, have shown that both probing methods may in principle provide a
reasonably good estimation of the moment of inertia. However, the tilting
mechanism is particularly sensitive to low-frequency modes. As a result, that
procedure demands prohibitively long integration times in order to provide
reliable estimations of the moment of inertia. In contrast, the study of the
scissors mode after suddenly stopping the rotation is much less sensitive to
low frequencies, and hence provides a much better alternative.

Similar procedures may be employed to study experimentally the moment of
inertia in the H-SS. We have evaluated for $N=2\times 10^{5}$ atoms the system
response during $T=1$s (the results are very consistent down to integration
times as short as $200$ ms, well within realistic experimental lifetimes). The
tilting procedure is very much affected by the very low-frequency modes
characteristic of the polar quasi-symmetry of the H-SS pattern. As a result,
that procedure provides a very bad estimation of $I_{\rm exact}/I_{\rm C}$;
for $a=82a_{0}$, the estimation is $5.36$ versus the correct $0.08$ value. A
much better result is obtained using the rotation technique (see the red-star
symbols, and panels (b) to (d) in Fig. 4). Note that, although there is still
a deviation from the actual value, the absolute difference with the exact
result is small. Hence, the rotation method can be experimentally applied to
reveal the much larger superfluid fraction of the H-SS compared to that of the
D-SS.

##  VI Conclusions

Although research on dipolar supersolids has been limited up to now to the
case of droplet supersolid arrays, experiments using larger number of atoms or
tighter traps, or working with polar molecules, may open in the near future
the possibility for studying other intriguing two-dimensional density patterns
in quantum-stabilized dipolar condensates. Of those, honeycomb supersolids are
particularly interesting, since, in contrast to droplet supersolids, their
superfluid fraction remains very large even for a large density contrast.
However, as discussed in this paper, the superfluidity of the honeycomb
supersolid cannot be probed by the observation of quantized vortices due to
the lack of robustness of the density pattern. Under rotation, vortex
nucleation results in the destruction of the honeycomb pattern into a
labyrinthic one. We have shown, however, that the large superfluid fraction of
the honeycomb supersolid may be probed by a proper monitoring of the response
after a scissors-like perturbation.

## Acknowledgments

We acknowledge support of the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) under Germany’s Excellence Strategy– EXC-2123
QuantumFrontiers–390837967, and FOR 2247.

## References

  * Tanzi  _et al._ (2019a) L. Tanzi, E. Lucioni, F. Famà, J. Catani, A. Fioretti, C. Gabbanini, R. N. Bisset, L. Santos, and G. Modugno, Phys. Rev. Lett. 122, 130405 (2019a).
  * Böttcher  _et al._ (2019) F. Böttcher, J.-N. Schmidt, M. Wenzel, J. Hertkorn, M. Guo, T. Langen, and T. Pfau, Phys. Rev. X 9, 011051 (2019).
  * Chomaz  _et al._ (2019) L. Chomaz, D. Petter, P. Ilzhöfer, G. Natale, A. Trautmann, C. Politi, G. Durastante, R. M. W. van Bijnen, A. Patscheider, M. Sohmen, M. J. Mark, and F. Ferlaino, Phys. Rev. X 9, 021012 (2019).
  * Leggett (1970) A. J. Leggett, Phys. Rev. Lett. 25, 1543 (1970).
  * Hertkorn  _et al._ (2021) J. Hertkorn, J.-N. Schmidt, M. Guo, F. Böttcher, K. S. H. Ng, S. D. Graham, P. Uerlings, T. Langen, M. Zwierlein, and T. Pfau, Phys. Rev. Research 3, 033125 (2021).
  * Zhang  _et al._ (2019) Y.-C. Zhang, F. Maucher, and T. Pohl, Phys. Rev. Lett. 123, 015301 (2019).
  * Zhang  _et al._ (2021) Y.-C. Zhang, T. Pohl, and F. Maucher, Phys. Rev. A 104, 013310 (2021).
  * Poli  _et al._ (2021) E. Poli, T. Bland, C. Politi, L. Klaus, M. A. Norcia, F. Ferlaino, R. N. Bisset, and L. Santos, Phys. Rev. A 104, 063307 (2021).
  * Norcia  _et al._ (2021) M. A. Norcia, C. Politi, L. Klaus, E. Poli, M. Sohmen, M. J. Mark, R. N. Bisset, L. Santos, and F. Ferlaino, Nature 596, 357 (2021).
  * Bland  _et al._ (2022) T. Bland, E. Poli, C. Politi, L. Klaus, M. A. Norcia, F. Ferlaino, L. Santos, and R. N. Bisset, Phys. Rev. Lett. 128, 195302 (2022).
  * Schmidt  _et al._ (2022) M. Schmidt, L. Lassablière, G. Quéméner, and T. Langen, Phys. Rev. Research 4, 013235 (2022).
  * Tanzi  _et al._ (2019b) L. Tanzi, S. M. Roccuzzo, E. Lucioni, F. Famà, A. Fioretti, C. Gabbanini, G. Modugno, A. Recati, and S. Stringari, Nature 574, 382 (2019b).
  * Guo  _et al._ (2019) M. Guo, F. Böttcher, J. Hertkorn, J.-N. Schmidt, M. Wenzel, H. P. Büchler, T. Langen, and T. Pfau, Nature 574, 386 (2019).
  * Natale  _et al._ (2019) G. Natale, R. M. W. van Bijnen, A. Patscheider, D. Petter, M. J. Mark, L. Chomaz, and F. Ferlaino, Phys. Rev. Lett. 123, 050402 (2019).
  * Guéry-Odelin and Stringari (1999) D. Guéry-Odelin and S. Stringari, Phys. Rev. Lett. 83, 4452 (1999).
  * Maragò  _et al._ (2000) O. M. Maragò, S. A. Hopkins, J. Arlt, E. Hodby, G. Hechenblaikner, and C. J. Foot, Phys. Rev. Lett. 84, 2056 (2000).
  * Tanzi  _et al._ (2021) L. Tanzi, J. G. Maloberti, G. Biagioni, A. Fioretti, C. Gabbanini, and G. Modugno, Science 371, 1162 (2021), https://www.science.org/doi/pdf/10.1126/science.aba4309 .
  * Norcia  _et al._ (2022) M. A. Norcia, E. Poli, C. Politi, L. Klaus, T. Bland, M. J. Mark, L. Santos, R. N. Bisset, and F. Ferlaino, Phys. Rev. Lett. 129, 040403 (2022).
  * Roccuzzo  _et al._ (2022) S. M. Roccuzzo, A. Recati, and S. Stringari, Phys. Rev. A 105, 023316 (2022).
  * Roccuzzo  _et al._ (2020) S. M. Roccuzzo, A. Gallemí, A. Recati, and S. Stringari, Phys. Rev. Lett. 124, 045702 (2020).
  * Gallemí  _et al._ (2020) A. Gallemí, S. M. Roccuzzo, S. Stringari, and A. Recati, Phys. Rev. A 102, 023322 (2020).
  * Šindik  _et al._ (2022) M. Šindik, A. Recati, S. M. Roccuzzo, L. Santos, and S. Stringari, “Creation and robustness of quantized vortices in a dipolar supersolid when crossing the superfluid-to-supersolid transition,”  (2022).
  * Klaus  _et al._ (2022) L. Klaus, T. Bland, E. Poli, C. Politi, G. Lamporesi, E. Casotti, R. N. Bisset, M. J. Mark, and F. Ferlaino, “Observation of vortices and vortex stripes in a dipolar bose-einstein condensate,”  (2022).
  * Petrov (2015) D. S. Petrov, Phys. Rev. Lett. 115, 155302 (2015).
  * Wächtler and Santos (2016) F. Wächtler and L. Santos, Phys. Rev. A 93, 061603 (2016).
  * Lima and Pelster (2012) A. R. P. Lima and A. Pelster, Phys. Rev. A 86, 063609 (2012).
  * Sinha and Castin (2001) S. Sinha and Y. Castin, Phys. Rev. Lett. 87, 190402 (2001).


# An interpretable deep learning method for bearing fault diagnosis

Hao Lu  Austin M. Bray  Chao Hu  Andrew T. Zimmerman  Hongyi Xu

###### Abstract

Deep learning (DL) has gained popularity in recent years as an effective tool
for classifying the current health and predicting the future of industrial
equipment. However, most DL models have black-box components with an
underlying structure that is too complex to be interpreted and explained to
human users. This presents significant challenges when deploying these models
for safety-critical maintenance tasks, where non-technical personnel often
need to have complete trust in the recommendations these models give. To
address these challenges, we utilize a convolutional neural network (CNN) with
Gradient-weighted Class Activation Mapping (Grad-CAM) activation map
visualizations to form an interpretable DL method for classifying bearing
faults. After the model training process, we apply Grad-CAM to identify a
training sample’s feature importance and to form a library of diagnosis
knowledge (or health library) containing training samples with annotated
feature maps. During the model evaluation process, the proposed approach
retrieves prediction basis samples from the health library according to the
similarity of the feature importance. The proposed method can be easily
applied to any CNN model without modifying the model architecture, and our
experimental results show that this method can select prediction basis samples
that are intuitively and physically meaningful, improving the model’s
trustworthiness for human users.

###### keywords:

Interpretable deep learning , Health monitoring , Fault diagnosis , Bearing ,
Physical knowledge

††journal: Journal of Manufacturing Systems

[ISU-ME]organization=Department of Mechanical Engineering, Iowa State
University,city=Ames, postcode=50011, state=IA, country=USA

[ISU-ECE]organization=Department of Electrical and Computer Engineering, Iowa
State University,city=Ames, postcode=50011, state=IA, country=USA

[UCONN]organization=Department of Mechanical Engineering, University of
Connecticut,city=Storrs, state=CT, postcode=06269, country=USA

[Percev]organization=Perc$\bar{e}v$, LLC,city=Davenport, postcode=52807,
state=IA, country=USA

[Grace]organization=Grace Technologies,city=Davenport, postcode=52807,
state=IA, country=USA

An interpretable DL method allows human users to understand model predictions
of bearing health

Two-step procedure: health library creation and prediction basis retrieval

Representative training samples are selected as prediction basis samples for
each test sample

Prediction basis samples have similar feature importance distributions to the
test sample

##  1 Introduction

Recent years have ushered in a continuous flow of new machine learning
research contributions and a steady increase in the number of implementations
of machine learning techniques in commercial applications. These advances have
been enabled in part by ever-faster computing equipment becoming available
that enable cost-effective machine learning deployments. Engineering
applications for machine learning techniques are no exception to this trend,
as machine learning-enabled health management and fault classification systems
are rapidly expanding to monitor and supervise physical assets such as
rotating equipment (e.g., bearings and gearboxes) [1, 2], structural systems
(e.g., bridges, buildings, and planes) [3, 4, 5], energy infrastructure (e.g.,
batteries and turbines) [6, 7, 8], and manufacturing processes [9].

In this study, we focus on using machine learning to supervise rotating
industrial equipment (e.g., pumps, motors, fans, or generators), as these
types of assets represent a specific mission-critical class of mechanical
system, the failure or which can lead to substantial unplanned downtime,
damage to other elements of the equipment, and a significant safety risk [8,
10]. The criticality of these rotating systems have provided the stimulus for
extensive research dedicated to monitoring, classifying, and predicting the
health of rotating equipment [11, 12, 13].

In recent years, deep learning (DL) - a subset of machine learning that
employs layered artificial neural networks to model complex data abstractions
- has gained popularity because of its excellent performance across different
fields [14, 15]. Despite its promising performance, DL models build the input-
output knowledge by optimizing model parameters via the gradient
backpropagation algorithm. However, the huge amount of parameters inside the
DL models makes it too complex to understand the model’s decision-making
process, which creates challenges for end users of the models interpreting
model outputs and determining if results can be trusted in critical
applications [16, 17]. This challenge is amplified in fields including
biomedical engineering and medicine, where scientific ethics demands a high
level of trust in machine learning models prior to making medical decisions
informed by model outputs [18, 19]. Similarly, for safety-critical production
tasks or maintenance decisions where rotating equipment failure could lead to
millions of dollars in lost revenue, an operator requires far more information
from a trustworthy model than a simple prediction result.

Many efforts have been made to improve the performance and subsequent trust
that can be placed in the outputs of a DL model. One common approach to
improving model trust involves increasing the quality and quantity of training
dataset. Unfortunately, generating a sufficient quantity of training data is
either cost-prohibitive or physically impossible [20]. Interpretable DL
modeling techniques were developed to improve the model’s trustworthiness
without increasing the amount of training data. By generating an improved
model with outputs that are interpretable by the end user, it becomes possible
to trust a model that was trained on a smaller but obtainable dataset.

Multiple techniques can be utilized to develop interpretable DL methods, such
as guided backpropagation, deconvolution, and class activation maps. Class
Activation Mapping (CAM) enhance model explainability by generating activation
maps via the weighted sum of the convolutional feature maps before the final
output layer [21]. The CAM is applicable for CNNs which do not contain any
fully-connected layers. Gradient-weighted Class Activation Mapping (Grad-CAM)
was developed as a generalization of CAM that can provide visual explanations
to CNNs with fewer restrictions regarding model architecture [22]. This
visualization enables the identification of the regions of a sample most
important to a trained model’s prediction (classification or regression).

Selvaraju et al. implement Grad-CAM as a gradient-based visual explanation
which provides a localization map identifying important regions for predicting
the sample’s class [23]. This results in improvements to model training
procedures and makes the model output more explainable to end users. Camalan
et al. (2021) utilize CAMs to generate heat maps overlaid on test images to
identify the location of oral cancer and calculate a patch around the region
[24]. Lee et al. (2018) propose a novel Pyramid Gradient-based Class
Activation Mapping (PG-CAM) method that overcomes a lack of highly annotated
data by training in a weakly supervised state [20].

The aforementioned studies are limited to generating activation heatmaps
visualizing important regions of an input sample and still leave challenges
for end users looking to fully understand and trust a model’s decision-making
process. For example, an activation map generated by Grad-CAM can provide end
users with an indication of important regions for a particular sample, but
this leaves the end user with little context when comparing other samples of
the same classification (i.e., how the activation map differs for samples of
each classification or how activations may differ between similarly classified
samples).

A user-friendly interpretable DL method aims to provide training samples that
are similar to the test sample during the model evaluation process. Kim et al.
(2022) utilize CAM to identify regions of an image most important for
classifying chest X-rays [25]. Then, the model’s decision is supported by
comparing an activated region of the test X-ray to a library of activation
patches generated from training data. From this library a prediction basis is
selected that provides images of similarly activated image patches of the same
classification. In a similar fashion, Lee et al. (2019) utilize CAM activation
maps as an explainable element to a CNN in order to classify hemorrhage
conditions within human brain CT scans [18]. The authors generate a library of
training activation maps selected based on a relevance count metric determined
for each classification and for activation maps within each layer of the CNN.
For a given test sample, prediction basis samples can be selected by
calculating the Euclidean distance between the normalized test activation and
each library entry. From these distances, the most similar entries are
selected to justify the model classification for the end user.

In this paper, we investigate a novel interpretable DL method applied for the
first time in the field of bearing fault diagnosis. First, a CNN model is
trained using the training bearing dataset. Then, we calculate each sample’s
Grad-CAM activation map to generate an activation vector. Each sample’s
activation vector and its predicted health condition are utilized to create
the health library. During the model evaluation process, the proposed method
compares the activation vector’s similarity of testing samples with each
library entry. The entries most similar to the test activation vector are
selected as the prediction basis samples.

##  2 Background

This section discusses approaches implemented for data pre-processing, DL
model design, and Grad-CAM calculation and comparisons. Section 2.1 introduces
the calculation of physics-informed fault frequencies. Section 2.2 discusses
the signal processing steps that are required to convert the vibration signal
from time domain to order domain via envelope analysis and resampling. The CNN
architecture is then overviewed in section 2.3 with the subsequent calculation
of Grad-CAM activation maps discussed in section 2.4.

###  2.1 Bearing Fault Characteristic Frequencies

Rolling element bearings are key components in rotating mechanical systems. A
rolling element bearing contains four main components: an inner race, an outer
race, rolling elements (e.g., balls or rollers), and a cage that maintains
equal spacing between rolling elements [26]. If a defect presents in any
bearing component during bearing operation, the contact between the defect
area and other components will generate periodic vibrations. Studies have
shown that the frequency of that periodic vibration, denoted as bearing fault
characteristic frequency, can be calculated based on the geometry of the
bearing [27, 28].

This study considers two common bearing fault types: outer race fault and
inner race fault. Their related fault characteristic frequencies are
calculated and marked as $f_{\mathrm{BPFO}}$ and $f_{\mathrm{BPFI}}$,
respectively. Bearing characteristic frequencies serve as invaluable features
for identifying fault conditions within roller bearings [29]. With known
parameters, including bearing shaft speed $f_{\mathrm{r}}$, bearing inner race
diameter $d$, outer race diameter $D$, load angle from the radial plane
$\phi$, and the number of rollers $n$, $f_{\mathrm{BPFO}}$ and
$f_{\mathrm{BPFI}}$ can be calculated using the following equations:

$f_{\mathrm{BPFO}}=\frac{(n\times
f_{\mathrm{r}})}{2}[1-\frac{d}{D}\times\cos(\phi)]$ (1)

$f_{\mathrm{BPFI}}=\frac{(n\times
f_{\mathrm{r}})}{2}[1+\frac{d}{D}\times\cos(\phi)]$ (2)

###  2.2 Signal Processing

Localized bearing faults creating fault characteristic frequencies cause
amplitude modulation which results from a combination of the bearing fault
frequency as well as the resonance of the entire system [27]. This reality
makes fault frequencies challenging to inspect from raw frequency spectrum
data. Data processing is employed to address this challenge and improve
usability of frequency data.

Envelope analysis is a technique capable of demodulating bearing vibration
data [30, 31]. The first step of this process is to transform each time series
sample $x(t)$ via a Hilbert transform shown in Eq. (3). This process converts
$x(\tau)$, which is the sample at time $\tau$, to generate Hilbert-transformed
$\hat{x}(t)$ [32]. Envelope analysis is then completed by taking the absolute
value of the Hilbert-transformed sample.

$\hat{x}\left(t\right)=\frac{1}{\pi}\int_{-\infty}^{\infty}{\frac{x\left(\tau\right)}{t-\tau}d\tau}$
(3)

To obtain data in the frequency domain, $x(\omega)$, a Fourier transform is
completed as shown in Eq. (4). This step takes the Hilbert-transformed time
domain sample $\hat{x}(t)$ and results in a Fourier-transformed sample in the
frequency domain $f(\omega)$.

$f\left(\omega\right)=\int_{-\infty}^{\infty}{\hat{x}\left(t\right)e^{-j\omega
t}dt}$ (4)

As each bearing sample’s characteristic frequency is a function of the
constant shaft speed, $f_{r}$, order frequency normalization is utilized to
transform the data to be measured in orders of shaft angle. The transformation
from frequency to order follows Eq. 5 where $f$ is a given frequency prior to
being transformed to order domain [33].

$o=\frac{f}{f_{r}}$ (5)

Figure 1: Processed envelope order spectrum.

This pre-processed vibration data is now suitable for model training. By
working in the frequency domain and completing an envelope analysis, bearing
characteristic frequencies can be overlaid on the bearing data. This will also
facilitate the identification of class activation values near these
characteristic frequencies, which will be discussed later in Sec. 3.1.

###  2.3 CNN Architecture for Bearing Health Classification

We utilize a 1D-CNN model to classify bearing health based on pre-processed
vibration signals. Each pre-processed vibration signal is in the form of an
envelope order spectrum (see an example in Fig. 1). The 1D-CNN model comprises
a feature extractor, a global max pooling layer, and a classifier. The feature
extractor consists of a series of three convolutional blocks. Each
convolutional block is composed of a 1D convolutional layer, a batch
normalization layer, and a max pooling layer. The input and output of the
$l^{\mathrm{th}}$ convolutional layer is denoted as $\mathbf{x}^{l}$ and
$\mathbf{y}^{l}$, respectively.

As the first layer inside the $l^{\mathrm{th}}$ convolutional block, the 1D
convolutional layer applies a convolution operation to the input feature map
$\mathbf{x}$, generating an output feature map $\mathbf{y}$, expressed as:

$\mathbf{y}^{l}_{b}=\mathrm{\phi}\left(\sum_{a}\mathbf{x}^{l}_{a}\ast\mathbf{w}^{l}_{ab}\right)$
(6)

where $\mathbf{x}^{l}_{a}$ denotes the $a^{\mathrm{th}}$ channel of the input
feature map, $\mathbf{y}^{l}_{b}$ denotes the $b^{\mathrm{th}}$ channel of the
output feature map, $\mathbf{w}^{l}_{ab}$ denotes the 1D convolutional kernel
weights (from the $a^{\mathrm{th}}$ channel of the input feature map to the
$b^{\mathrm{th}}$ channel of the output feature map), $\ast$ denotes a
convolution operation, and $\mathrm{\phi}$ is an activation function. Note
that the input feature map of the first convolutional layer, $\mathbf{x}^{1}$,
is a vector of vibration amplitude readings from an envelope order spectrum,
which can be viewed as a feature vector.

Nonlinear functions are utilized as activation functions in CNNs to give the
network the ability to model nonlinear input-output relationships. [34].
Popular nonlinear activation functions, including sigmoid, tanh, and rectified
linear unit (ReLU), can be used within CNNs [35]. The ReLU activation function
shown in Eq. (7) is adopted in our work. Note that $z$ in Eq. (7) is the
scalar input to the ReLU activation function.

$\mathrm{ReLU}(z)=\max{(0,z)}$ (7)

In each convolutional block, the batch normalization layer is implemented
after the 1D convolutional layer to improve model training stability. Then the
output of the batch normalization layer is fed as the input to the 1D max
pooling layer. The 1D max pooling layer extracts the maximum feature value
within a sub-region by sliding a fixed-length window over its input feature
map.

In this study, three convolutional blocks are stacked one by one to extract
lower-dimensional features. Following the feature extractor, the global max
pooling layer calculates the maximum value within each channel to generate a
downsampled feature map [36].

As the last component of the CNN, the classifier is designed as a single fully
connected layer plus a softmax function. The fully connected layer calculates
the weighted sum across all the inputs and adds a bias. Then, the softmax
function converts the output of the fully connected layer to class
probabilities, i.e., a probability distribution over predicted health classes.
Specifically, the softmax function outputs a probability for each health class
$c$ with $c=1,...,C$, where $C$ denotes the number of possible health classes.

###  2.4 Grad-CAM Activation Map Visualization

Grad-CAM visualizations of class activation maps rely on average gradients as
weights applied to a map. To calculate a Grad-CAM map for a specific class $c$
out of $C$ possible class labels, the score for class $c$, $y^{c}$, and the
activation map $A_{i}^{k}$ for each $i^{\mathrm{th}}$ element and
$k^{\mathrm{th}}$ filter are retrieved. These values are calculated within a
CNN during the back-propagation step as $\frac{\partial y^{c}}{\partial
A_{i}}$ and are subsequently retrievable without significant additional
computational overhead [22].

$\alpha_{k}^{c}=\frac{1}{C}\sum_{i}\frac{\partial y^{c}}{\partial A_{i}^{k}}$
(8)

Utilizing the activation map $A_{i}$ retrieved from the back-propagation step
as well as the importance-weight $\alpha_{k}^{c}$ for class $c$ and filter
$k$, the overall class $c$ importance-weight, $L_{\mathrm{Grad-CAM}}^{c}$, can
be calculated using Eq. (9).

$L_{\mathrm{Grad-CAM}}^{c}=\sum_{k}\alpha_{k}^{c}A^{k}$ (9)

Activation maps can be generated for each block within a CNN model. Displaying
activation maps for different blocks results in varied resolution and detail
of the activation map. Later block’s feature maps include less noise as a
result of passing through additional max pooling layers in each block [37]. As
a clearer, readable activation map is desired, the feature map from the last
convolutional layer is selected.

##  3 Methodology

The proposed method is a post hoc interpretable DL method applicable to any
CNN-based model without modifying the model architecture. An overview of the
proposed method is provided in Fig. 2. The proposed interpretable DL method
has two major steps: (1) generating a library of candidate prediction basis
samples (Sec. 3.1) and (2) retrieving prediction basis samples from the
library based on the similarity between the feature importance of a test
sample and all training samples in the library (Sec. 3.2).

Figure 2: The overview of proposed interpretable DL method

###  3.1 Creation of Health Library

Step 1 of our proposed interpretable DL method constructs a health library
based on the class activation maps generated by the last convolutional layer
of a trained CNN. The health library is designed as a set with
$N_{\mathrm{train}}$ entries, where $N_{\mathrm{train}}$ denotes the number of
training samples. Each entry consists of an activation vector that reflects
the importance of input features and the training sample’s health class. Two
interpretable DL algorithms are investigated in this study: CAM-Full and CAM-
Sub. These two algorithms vary in the activation vector for the health library
creation and prediction basis retrieval.

CAM-Full creates a health library and retrieves prediction basis samples using
activation vectors associated with input features (i.e., vibration amplitudes
on envelop spectra) across the entire order range $[0,30]$. CAM-Sub differs
from CAM-Full by only considering the feature importance values within the
pre-defined frequency sub-bands corresponding to the sample’s fault type. The
remainder of this subsection describes the design of fault frequency sub-
bands.

As described in Sec. 2.1, features located close to the bearing fault
characteristic frequencies are important in diagnosing the bearing health. To
leverage that domain knowledge, CAM-Sub focuses on the feature importance
values associated with input features (vibration amplitudes) near the bearing
fault characteristic frequencies. For a given fault type whose health class is
${c}$, the selected input features fall within three frequency sub-bands that
are defined using the fault characteristic frequency $f^{c}$ and sub-band
width $\epsilon$. These three sub-bands are centered at $f^{c}$ and it’s first
two harmonics ($2f^{c}$ and $3f^{c}$), respectively. The width of each sub-
band is defined as $2\epsilon$ of the sub-band center frequency, and each sub-
band extends symmetrically in both directions by $\epsilon$ times the
corresponding center frequency. After identifying the order range of the three
sub-bands for each fault type, the feature importance (or class activation)
values in the sub-bands are chosen to create the CAM-Sub health library.

Figure 3: Fault sub-bands (green) overlaid on a pre-processed outer race
sample for the bearing outer race frequency $f_{\mathrm{BPFO}}$ and its first
two harmonics.

Figure 3 shows the envelope spectrum of an outer race fault sample, annotated
with the three fault frequency sub-bands. The sub-bands include the first
outer race fault characteristic frequency and its following two harmonics.
Note that the frequency sub-band increases in width (by $2f^{c}\epsilon$) for
each harmonic. The increasing sub-band width ensures that the first two
harmonics of the fault characteristic frequency fall within the respective
sub-bands in the presence of errors in the shaft speed measurement.

###  3.2 Retrieval of Prediction Basis Samples

Step 2 of our proposed interpretable DL method takes place after a DL model
has been trained. Given a test sample, this step calculates the sample’s class
activation (or feature importance) vector using the CAM-Full or CAM-Sub
algorithm described in Sec. 3.1, then evaluates the similarity between the
class activation vector of this test sample and the class activation vector of
each training entry in the health library that has the same health class as
the test sample’s predicted class. The entries most similar to the test
activation vector are selected as the prediction basis samples.

Given a test dataset with $N_{\mathrm{test}}$ samples, Grad-CAM is used to
determine each test sample’s class activation vector
$\mathbf{TestCA}^{c_{j}}_{j}$, where $c_{j},j=1,...,N_{\mathrm{test}}$,
denotes the $j^{\mathrm{th}}$ test sample’s predicted class. Then, we select
entries in the health library with the same health class $c_{j}$, and their
training sample indices form a set of $N_{\mathrm{train}}^{c_{j}}$ indices,
denoted as $\mathrm{I}^{c_{j}}$, which satisfies the following condition:
$\forall i\in\mathrm{I}^{c_{j}},c_{i}=c_{j}$. Here,
$N_{\mathrm{train}}^{c_{j}}$ is the number of training numbers whose health
classes are ${c_{j}}$. The class activation vectors of these
$N_{\mathrm{train}}^{c_{j}}$ library entries also form a set, denoted as
$\mathrm{S}^{c_{j}}=\\{\mathbf{TrainCA}^{c_{i}}_{i}\\}_{i\in\mathrm{I}^{c_{j}}}$.
Training samples with the health class $c_{j}$ provide critical information
backing the classification of the test sample into $c_{j}$, assuming that we
do not have test samples that fall outside the training data distribution.
Considering the health library may contain training data collected from
different bearings under different operating conditions, the class activation
vector may vary drastically from one training sample to another, even for
those with the same health class; thus, it is expected that not all training
samples are important to the model prediction on a test sample. To identify
prediction basis samples that are highly relevant on a per-test-sample basis,
we calculate the similarity between $\mathbf{TestCA}^{c_{j}}_{j}$ and
$N_{\mathrm{train}}^{c_{j}}$ class activation vectors in $\mathrm{S}^{c_{j}}$.
Entries most similar to $\mathbf{TestCA}^{c_{j}}_{j}$ are selected as the
prediction basis samples. Note again that Sec. 3.1 describes the derivation of
class activation vectors with two algorithms: CAM-Full and CAM-Sub. The key
difference between these two algorithms is that CAM-Full selects class
activation values associated with vibration amplitudes across the entire order
range $[0,30]$, while CAM-Sub concatenates the class activation values within
the three pre-defined frequency sub-bands for a given fault type.

Comparisons between two samples (a test sample and an entry in the health
library) are achieved by calculating the Euclidean distance between their
normalized class activation vectors. The vector normalization is completed
using Eqs. (10) and (11), where the class activation vector of the
$j^{\mathrm{th}}$ test sample, $\mathbf{TestCANorm}_{j}^{c_{j}}$, and the
class activation vector of each library entry,
$\mathbf{TrainCANorm}_{i}^{c_{i}}\in\mathrm{S}^{c_{j}}$, are normalized by
their respective L2 norms (or Euclidean norms). Utilizing vector normalization
minimizes the effect of vector length when measuring the similarity between
two vectors. This treatment may be beneficial as our focus is comparing the
distributions of relative feature importance rather than the absolute feature
importance values.

$\mathbf{TestCANorm}_{j}^{c_{j}}=\frac{\mathbf{TestCA}_{j}^{c_{j}}}{||\mathbf{TestCA}_{j}^{c_{j}}||}$
(10)

$\mathbf{TrainCANorm}_{i}^{c_{i}}=\frac{\mathbf{TrainCA}_{i}^{c_{i}}}{||\mathbf{TrainCA}_{i}^{c_{i}}||}$
(11)

After the two class activation vectors are L2-normalized, the Euclidean
distance between these normalized vectors measures the class activation
similarity between the test sample and training entry. This similarity measure
takes the following form:

$\mathbf{Dis}_{i,j}=||\mathbf{TestCANorm}_{i}^{c_{i}}-\mathbf{TrainCANorm}_{j}^{c_{j}}||$
(12)

The Euclidean distance $\mathbf{Dis}_{i,j}$ between a given test sample $j$
and training sample $i$ is calculated for each of the
$N_{\mathrm{train}}^{c_{j}}$ training entries in the set $\mathrm{S}^{c_{j}}$.
Then, these $N_{\mathrm{train}}^{c_{j}}$ selected entries,
{$\mathbf{TrainCA}_{i}^{c_{i}}\\}_{i\in\mathrm{I}^{c_{j}}}$, are sorted in
ascending order of $\mathbf{Dis}_{i,j}$. Finally, the prediction basis samples
are retrieved by selecting the top $K$ $(\ll N_{\mathrm{train}}^{c_{j}})$
entries with the smallest Euclidean distances to the test sample in the
normalized class activation vector space. These representative training
samples from the health library with similar distributions of feature
importance can be used to justify the prediction of the CNN on the test
sample.

##  4 Case Study

###  4.1 Design of Case Study

To assess the effectiveness of the proposed interpretable DL method in
retrieving trustworthy prediction basis samples, we designed experiments using
the open-source dataset from the Paderborn University (PU) Bearing DataCenter.
As one of the standard references in the bearing diagnosis field, the PU
bearing dataset has been widely used to evaluate bearing fault diagnosis
algorithms. The bearing health condition can be classified as healthy, inner
race fault, outer race fault, roller fault, or a combination of faults. The PU
dataset generates damaged bearings in two ways: (1) using drilling/electrical
discharge machines to introduce artificial damage and (2) performing
accelerated life tests to introduce real damage [38].

The dataset used in our experiments included bearings subjected to accelerated
life tests, shown in Table 4. The data collected from each bearing consisted
of three vibration time series samples corresponding to three operating
conditions. Each time series sample was cut into smaller samples, with each
smaller sample containing 6,400 measurement points. Our signal processing step
then converted these vibration amplitude measurements from the time domain
into the order domain. We then selected the features within the order range
$[0,30]$ as the model input. After signal processing, 80% of the samples were
used as the training dataset, and the remainder were used as the test dataset.

Figure 4: A summary of PU dataset samples selected for this case study.

###  4.2 Evaluation of Prediction Basis Selection

The proposed method aims to provide end users with the model’s prediction for
each test sample and intuitively show the prediction basis samples, i.e., a
set of training samples similar to the test sample that support the model’s
prediction. The prediction basis samples are expected to reveal information
most relevant to the prediction while being easily comprehensible to end
users. We conducted experiments to assess the performance of the proposed
interpretable DL method in two aspects: (1) the ability to identify prediction
basis samples that are important with respect to making accurate predictions
on the entire test dataset (Sec. 4.2.1) and (2) the intuitiveness of
identified prediction basis samples (Sec. 4.2.2).

####  4.2.1 Importance of Prediction Basis Samples

We first evaluated the effectiveness of the proposed interpretable DL method
in selecting the _most relevant_ training samples (prediction basis) that
contribute, on average, the most to model predictions on the entire test
dataset. The basic idea of our evaluation was to compare the originally
trained model with models trained after removing some selected samples from
the original training dataset. Specifically, the evaluation was carried out in
the following three steps:

  1. 1.

Given each test sample, we used Grad-CAM to compute its feature importance. We
then calculated the Euclidean distance between the feature importance vector
of the selected test sample and the saved importance vector of each training
sample in the health library. We repeated the above process for each test
sample to obtain the Euclidean distance between each pair of test and training
samples, generating a distance matrix $\mathbf{Dis}$ of size
$N_{\mathrm{train}}\times N_{\mathrm{test}}$, where $\mathbf{Dis}_{i,j}$
denotes the Euclidean distance between the $i^{\mathrm{th}}$ training sample
and $j^{\mathrm{th}}$ test sample, and $N_{\mathrm{train}}$ and
$N_{\mathrm{test}}$, respectively, denote the numbers of training and test
samples.

  2. 2.

For each training sample, we computed the average of the distances of the
training samples to the test samples with the same label
$\mathbf{AvgDis}_{j}=\frac{1}{N_{\mathrm{train}}}\sum_{i=1}^{N_{\mathrm{train}}}\mathbf{Dis}_{i,j}$.
This average distance measures the importance of a training sample to model
predictions on test data, on average, across all test samples combined. Then,
we ranked the training samples according to their average distances to the
test data (i.e., their overall importance to the predictions on the test
samples), removed the top $K\%$ ($K=10,15,20,...,45$) training samples, and
trained a new DL model with this newly formed, reduced training dataset.

  3. 3.

The third and last step evaluated the newly trained model using the test
dataset (original, unchanged) and calculated the accuracy difference between
the newly trained and original models.

Our intuition was that removing training samples of higher importance is
expected to cause a more significant decrease in model accuracy. If the
Euclidean distance metric in Eq. (12) allows for effectively identifying high-
importance training samples, we expect a more rapid accuracy decrease than
randomly removing training samples. Note that the proposed interpretable DL
method selects prediction basis samples according to this Euclidean distance
metric. If this metric proves effective, the selected prediction basis samples
will likely contain more relevant information that contributes to correctly
classifying the test sample’s health. The baseline in our comparison was a
random selection method, denoted as “Random”, which randomly selects a certain
fraction of training samples as the prediction basis for sample removal and
model re-training. To account for the run-to-run variation in the model
training process, we independently re-trained ten new models for each of the
three approaches (i.e., Random, CAM-Full, and CAM-Sub).

The comparison results are summarized in Fig. 5. These two plots visualize how
the classification accuracy (Fig. 5(a)) and loss (Fig. 5(b)) on the test
dataset change as the fraction of removed training samples increases. Suppose
the Euclidean distance metric in Eq. (12) is a good measure of training sample
importance. In that case, the classification accuracy (loss) should decrease
(increase) monotonically and rapidly (at least initially) with the fraction of
removed high-importance samples. This seems to be the case for CAM-Full and
CAM-Sub, as shown in Fig. 5. The decreasing (increasing) trend of the
classification accuracy (loss) by CAM-Sub seems slightly more rapid than that
by CAM-Full, suggesting that CAM-Sub performed marginally better than CAM-Full
in identifying important training samples. This better performance by CAM-Sub
is also evidenced by another observation that model “re-training” after
removing important training samples identified by CAM-Sub exhibited larger
run-to-run variation than CAM-Full and Random. In other words, removing
prediction basis samples caused increases in the instability of model training
with respect to the test dataset, indicating these removed training samples
are highly important to ensure consistent accuracy on the test samples in the
presence of randomness in the model training process. In contrast, the random
selection method that does not consider sample importance produced an almost
flat curve for accuracy and loss. This observation provides evidence that
looking at the similarity between a test sample’s class activation and a
training sample’s class activation offers a meaningful way to measure the
importance of the training sample to the prediction for the test sample.

Figure 5: Test accuracy and loss of re-trained models vs. fraction of training
samples removed. A more rapidly decreasing trend of accuracy or increasing
trend of loss indicates a more effective ranking of training samples according
to their overall importance. The error bars for Random and CAM-Sub at each
fraction number of interest are offset slightly and horizontally for ease of
visualization.

So far, we have analyzed how the overall classification accuracy and loss
change as functions of the fraction of training samples removed. Additionally,
we examined a detailed summary of classification results in the form of a
confusion matrix. Given a confusion matrix $\mathbf{C}$, the entry
$\mathbf{C}_{i,j}$ is equal to the number of samples with true label $i$
classified as label $j$.

Figure 6: Confusion matrices produced by DL models trained using varying
quantities of data. Three methods (Random, CAM-Full, and CAM-Sub) were
implemented to remove 20% and 40% of all training samples to produce the
confusion matrices on the right-hand side.

The confusion matrices produced by the three methods when removing 20% and 40%
of all training samples are summarized in Fig. 6. No significant changes in
the models’ performance were observed after removing the top 20% of samples
from the training dataset.

After removing the top 40% of samples from the prediction basis pool, the
models developed by CAM-Sub and CAM-Full yielded substantially lower
prediction accuracy on inner race fault samples. Compared with the CAM-Full
models, the CAM-Sub models exhibited more noticeable decreases in the
classification accuracy on the healthy and outer race fault samples.
Considering the training samples removed consisted of healthy, inner race
fault, and outer race fault samples, this performance degradation indicates
CAM-Sub outperformed CAM-Full in identifying training samples of high
importance with respect to the other two health classes (i.e., healthy and
outer race fault).

####  4.2.2 Intuitiveness of Prediction Basis Samples

The second aspect of our method evaluation focused on how intuitive the
prediction basis samples were and how well these samples could help end users
interpret the model’s decision-making process. We selected a representative
outer race fault test sample and an inner race fault test sample, and used
CAM-Full and CAM-Sub to generate the top four prediction basis samples for
each selected test sample.

Let us first look at the outer race fault test sample. Figure 7 shows the
envelope order spectrum of the test sample (top row) and the prediction basis
samples CAM-Full and CAM-Sub retrieved from the health library (four bottom
rows), all overlaid with the feature importance by Grad-CAM. High amplitude
peaks can be observed at orders 3.05, 6.10, and 9.15, which are the outer race
fault frequency and its first two harmonics. The DL model correctly classified
the bearing health condition for this selected test sample. Both CAM-Full and
CAM-Sub retrieved prediction basis samples with similar activation maps to the
test sample. These prediction basis samples would allow human users to
understand the model’s decision-making process, particularly its justification
for this prediction. Such understanding can be obtained by comparing the
fault-signature similarities between the test and prediction basis samples in
the order ranges of interest (e.g., the frequency sub-bands for the outer race
fault).

Figure 7: Envelop spectra of an outer race fault test sample (a) and its four
prediction basis samples retrieved by CAM-Full (b) and CAM-Sub (c)

We can make a similar conclusion on the intuitiveness of prediction basis
samples when looking at the inner race fault test sample, as shown in Figure
8. The main difference is that compared to the outer race fault sample shown
in Fig. 7, the inner race fault sample contained more complex fault
signatures, for example, peaks outside of the inner race fault-related sub-
bands and weaker peaks. This observation is consistent with prior literature
and may be attributed to the fact that the inner race is not stationary like
the outer race [39]. CAM-Full selected prediction basis samples according to
the Euclidean norm of the vector difference of the entire activation map.
Although the four selected training samples were supposed to be the closest to
the test sample in terms of the normalized activation (or feature importance)
vector, the importance weight distributions of these training samples appeared
different from that of the test sample. In particular, the DL model assigned
high weights to vibration amplitudes at around order 4.9 (the first fault
frequency sub-band), i.e., the inner race fault frequency, for the test sample
(see Fig. 8(a)), while in contrast, very low importance weights were
associated with features in the first fault frequency sub-band for the
prediction basis samples retrieved by CAM-Full. Looking at the Grad-CAM
outputs of the top four prediction basis samples by CAM-Sub, we observe strong
activation (high feature importance) at around order 4.9. This observation
suggests that CAM-Sub selected training samples with similar class activation
maps to the test sample’s class activation map in the order ranges of high
relevance to the inner race fault. The results are expected as CAM-Sub selects
prediction basis samples according to the similarity of (normalized)
importance weight vectors within fault frequency sub-bands. Overall, the
prediction basis samples retrieved by CAM-Sub make more intuitive and physical
sense than those by CAM-Full, at least for this test sample.

Figure 8: Envelop spectra of an inner race fault test sample (a) and the
prediction basis samples retrieved by CAM-Full (b) and CAM-Sub (c)

##  5 Conclusion

In this study, we have proposed an interpretable DL technique capable of
creating a prediction basis for each test sample by ranking the 1D activation
maps of training samples. These prediction basis samples are a subset of the
training dataset deemed most relevant to the DL model’s prediction. Our work
furthers our understanding of how to leverage Grad-CAM activation maps to
enhance the interpretability of convolutional neural networks (CNNs) in
bearing fault classification tasks. A unique benefit offered by our method is
an extension of model interpretability from evaluating feature importance
(Grad-CAM) to evaluating training sample importance. This extension was made
possible by comparing the (normalized) class activation vector of a test
sample with those of training samples across the entire frequency range (CAM-
Full) or pre-defined fault frequency sub-bands (CAM-Sub). This comparison
results in a prediction basis set of training samples (or simply prediction
basis samples) that are expected to contribute the most to the mode’s
prediction on the test sample. The effectiveness of both CAM-Full and CAM-Sub
was evaluated using the Paderborn University (PU) Bearing DataCenter’s
experimental dataset and validated against randomized sample selection. Our
future work will examine ways to account for model uncertainty in retrieving
prediction basis samples and investigate model interpretability when applying
deep learning to predict bearing remaining useful life.

## Acknowledgments

This work was supported in part by the U.S. National Science Foundation under
Grant IIP-2222630. Any opinions, findings, or conclusions in this paper are
those of the authors and do not necessarily reflect the sponsor’s views. We
also thank Dr. Venkat Nemani and Dr. Mohammad Behtash for the helpful
discussions on interpretable DL for bearing health monitoring.

## References

  * [1] B. Wang, Y. Lei, N. Li, T. Yan, Deep separable convolutional network for remaining useful life prediction of machinery, Mechanical systems and signal processing 134 (2019) 106330.
  * [2] Z. Liu, L. Zhang, A review of failure modes, condition monitoring and fault diagnosis methods for large-scale wind turbine bearings, Measurement 149 (2020) 107002.
  * [3] Y. Wang, C. Gogu, N. Binaud, C. Bes, R. T. Haftka, N. H. Kim, A cost driven predictive maintenance policy for structural airframe maintenance, Chinese journal of aeronautics 30 (2017) 1242–1257.
  * [4] W. Yiwei, G. Christian, N. Binaud, B. Christian, R. T. Haftka, et al., A cost driven predictive maintenance policy for structural airframe maintenance, Chinese Journal of Aeronautics 30 (3) (2017) 1242–1257.
  * [5] Y.-W. Wang, Y.-Q. Ni, S.-M. Wang, Structural health monitoring of railway bridges using innovative sensing technologies and machine learning algorithms: A concise review, Intelligent Transportation Infrastructure 1 (2022).
  * [6] Y. Liu, B. Guo, X. Zou, Y. Li, S. Shi, Machine learning assisted materials design and discovery for rechargeable batteries, Energy Storage Materials 31 (2020) 434–450.
  * [7] A. Downey, Y.-H. Lui, C. Hu, S. Laflamme, S. Hu, Physics-based prognostics of lithium-ion battery using non-linear least squares with dynamic bounds, Reliability Engineering & System Safety 182 (2019) 1–12.
  * [8] K. L. Tsui, N. Chen, Q. Zhou, Y. Hai, W. Wang, et al., Prognostics and health management: A review on data driven approaches, Mathematical Problems in Engineering 2015 (2015).
  * [9] C. Wang, X. Tan, S. B. Tor, C. Lim, Machine learning in additive manufacturing: State-of-the-art and perspectives, Additive Manufacturing 36 (2020) 101538.
  * [10] M. Cerrada, R.-V. Sánchez, C. Li, F. Pacheco, D. Cabrera, J. V. de Oliveira, R. E. Vásquez, A review on data-driven fault severity assessment in rolling bearings, Mechanical Systems and Signal Processing 99 (2018) 169–196.
  * [11] H. Khorasgani, G. Biswas, S. Sankararaman, Methodologies for system-level remaining useful life prediction, Reliability Engineering & System Safety 154 (2016) 8–18.
  * [12] K. Medjaher, D. A. Tobon-Mejia, N. Zerhouni, Remaining useful life estimation of critical components with application to bearings, IEEE Transactions on Reliability 61 (2) (2012) 292–302.
  * [13] N. Li, Y. Lei, J. Lin, S. X. Ding, An improved exponential model for predicting remaining useful life of rolling element bearings, IEEE Transactions on Industrial Electronics 62 (12) (2015) 7762–7773.
  * [14] C. Lou, M. A. Atoui, X. Li, Recent deep learning models for diagnosis and health monitoring: A review of research works and future challenges, Transactions of the Institute of Measurement and Control (2023) 01423312231157118\.
  * [15] S. Shamshirband, M. Fathi, A. Dehzangi, A. T. Chronopoulos, H. Alinejad-Rokny, A review on deep learning approaches in healthcare systems: Taxonomies, challenges, and open issues, Journal of Biomedical Informatics 113 (2021) 103627\.
  * [16] E. M. Kenny, C. Ford, M. Quinn, M. T. Keane, Explaining black-box classifiers using post-hoc explanations-by-example: The effect of explanations and error-rates in xai user studies, Artificial Intelligence 294 (2021) 103459.
  * [17] K. Aas, M. Jullum, A. Løland, Explaining individual predictions when features are dependent: More accurate approximations to shapley values, Artificial Intelligence 298 (2021) 103502.
  * [18] H. Lee, S. Yune, M. Mansouri, M. Kim, S. H. Tajmir, C. E. Guerrier, S. A. Ebert, S. R. Pomerantz, J. M. Romero, S. Kamalian, et al., An explainable deep-learning algorithm for the detection of acute intracranial haemorrhage from small datasets, Nature biomedical engineering 3 (3) (2019) 173–182.
  * [19] B. H. Van der Velden, H. J. Kuijf, K. G. Gilhuijs, M. A. Viergever, Explainable artificial intelligence (xai) in deep learning-based medical image analysis, Medical Image Analysis 79 (2022) 102470.
  * [20] S. Lee, J. Lee, J. Lee, C.-K. Park, S. Yoon, Robust tumor localization with pyramid grad-cam, arXiv preprint arXiv:1805.11393 (2018).
  * [21] B. Zhou, A. Khosla, A. Lapedriza, A. Oliva, A. Torralba, Learning deep features for discriminative localization, in: Proceedings of the IEEE conference on computer vision and pattern recognition, 2016, pp. 2921–2929.
  * [22] J. Kim, J.-M. Kim, Bearing fault diagnosis using grad-cam and acoustic emission signals, Applied Sciences 10 (6) (2020) 2050.
  * [23] R. R. Selvaraju, M. Cogswell, A. Das, R. Vedantam, D. Parikh, D. Batra, Grad-cam: Visual explanations from deep networks via gradient-based localization, in: Proceedings of the IEEE international conference on computer vision, 2017, pp. 618–626.
  * [24] S. Camalan, H. Mahmood, H. Binol, A. L. D. Araujo, A. R. Santos-Silva, P. A. Vargas, M. A. Lopes, S. A. Khurram, M. N. Gurcan, Convolutional neural network-based clinical predictors of oral dysplasia: Class activation map analysis of deep learning results, Cancers 13 (6) (2021) 1291.
  * [25] D. Kim, J. Chung, J. Choi, M. D. Succi, J. Conklin, M. G. F. Longo, J. B. Ackman, B. P. Little, M. Petranovic, M. K. Kalra, et al., Accurate auto-labeling of chest x-ray images based on quantitative similarity to an explainable ai model, Nature communications 13 (1) (2022) 1867.
  * [26] S. A. McInerny, Y. Dai, Basic vibration signal processing for bearing fault detection, IEEE Transactions on education 46 (1) (2003) 149–156.
  * [27] S. Shen, H. Lu, M. Sadoughi, C. Hu, V. Nemani, A. Thelen, K. Webster, M. Darr, J. Sidon, S. Kenny, A physics-informed deep learning approach for bearing fault detection, Engineering Applications of Artificial Intelligence 103 (2021) 104295.
  * [28] T. Wang, M. Liang, J. Li, W. Cheng, Rolling element bearing fault diagnosis via fault characteristic order (fco) analysis, Mechanical Systems and Signal Processing 45 (1) (2014) 139–153.
  * [29] A. Rai, S. H. Upadhyay, A review on signal processing techniques utilized in the fault diagnosis of rolling element bearings, Tribology International 96 (2016) 289–306.
  * [30] R. B. Randall, J. Antoni, Rolling element bearing diagnostics-a tutorial, Mechanical Systems and Signal Processing 25 (2011) 485–520.  doi:10.1016/j.ymssp.2010.07.017.
  * [31] P. Borghesani, R. Ricci, S. Chatterton, P. Pennacchi, A new procedure for using envelope analysis for rolling element bearing diagnostics in variable operating conditions, Mechanical Systems and Signal Processing 38 (2013) 23–35.  doi:10.1016/j.ymssp.2012.09.014.
  * [32] N. E. Huang, Z. Shen, S. R. Long, A new view of nonlinear water waves: the hilbert spectrum, Annual review of fluid mechanics 31 (1) (1999) 417–457.
  * [33] K. Fyfe, E. Munck, Analysis of computed order tracking, Mechanical systems and signal processing 11 (2) (1997) 187–205.
  * [34] W. Hao, W. Yizhou, L. Yaqin, S. Zhili, The role of activation function in cnn, in: 2020 2nd International Conference on Information Technology and Computer Application (ITCA), IEEE, 2020, pp. 429–432.
  * [35] U. Michelucci, Applied Deep Learning A Case-Based Approach to Understanding Deep Neural Networks, 1st Edition, Apress, Berkeley, CA, 2018.
  * [36] Z. Ma, D. Chang, J. Xie, Y. Ding, S. Wen, X. Li, Z. Si, J. Guo, Fine-grained vehicle classification with channel max pooling modified cnns, IEEE Transactions on Vehicular Technology 68 (4) (2019) 3224–3233.
  * [37] S. Li, T. Li, C. Sun, R. Yan, X. Chen, Multilayer grad-cam: An effective tool towards explainable deep neural networks for intelligent fault diagnosis, Journal of Manufacturing Systems 69 (2023) 20–30.
  * [38] C. Lessmeier, J. K. Kimotho, D. Zimmer, W. Sextro, Condition monitoring of bearing damage in electromechanical drive systems by using motor current signals of electric motors: A benchmark data set for data-driven classification, in: PHM Society European Conference, Vol. 3, 2016.
  * [39] F. Cong, J. Chen, G. Dong, M. Pecht, Vibration model of rolling element bearings in a rotor-bearing system for fault diagnosis, Journal of sound and vibration 332 (8) (2013) 2081–2097.