# Axion-like particles from primordial black hole evaporation and their
detection in neutrino experiments

Tong Li<EMAIL_ADDRESS>School of Physics, Nankai University, Tianjin
300071, China  Rui-Jia Zhang<EMAIL_ADDRESS>School of
Physics, Nankai University, Tianjin 300071, China

###### Abstract

The primordial black holes (PBHs) play as a novel source to radiate light
elementary particles of energies in the region of a few hundred MeV. We
explore the possibility that the light axion-like particles (ALPs) are
produced from PBH evaporation. The absorption of light ALPs in the underground
detector targets then induces energetic photoelectron signatures in current
and future neutrino experiments. Utilizing the PBH ALP event rate, we place
general exclusion limits on the axion couplings at Super-K and Hyper-K. We
also translate these limits into the upper bound on the fraction of DM
composed of PBHs $f_{\rm PBH}$.

##  I Introduction

Dark Matter (DM) is one of the most evidential facets of possible new physics
beyond the Standard Model (SM). Despite plenty of astrophysical evidence of DM
existence in the Universe, however, the constitution and the properties of DM
are still mystery. Up to now, there is no conclusive evidence of the weakly
interacting massive particle (WIMP) in DM direct detection (DD) experiments
Aalbers  _et al._ (2022); Aprile  _et al._ (2022). Much more effort has
recently been paid to the theoretical hypotheses beyond the WIMP and their
detection methods in experiments. The axion and primordial black hole (PBH)
are two of the most well-motivated alternative DM candidates.

Axions as a solution to the strong CP problem Baluni (1979); Crewther  _et
al._ (1979); Kim (1979); Shifman  _et al._ (1980); Dine  _et al._ (1981);
Zhitnitsky (1980); Baker  _et al._ (2006); Pendlebury  _et al._ (2015) receive
a wide interest in phenomenology and cosmology. Both the QCD axion Peccei and
Quinn (1977a, b); Weinberg (1978); Wilczek (1978) (for a recent review see
Ref. Di Luzio  _et al._ (2020a)) and axion-like particles (ALPs) Kim (1987);
Kuster  _et al._ (2008) can play as cold DM through the misalignment
mechanism. The general ALPs span a wide range of the mass and the coupling
constants in a variety of theories, from sub-eV wavelike axion to MeV or even
heavier point particle.

The PBHs formed from the collapse of local overdensities in the early Universe
can serve as hypothetical macroscopic objects in DM halos Zel’dovich (1967);
Carr and Hawking (1974); Carr (1975); Khlopov (2010) (see Ref. Carr and Kuhnel
(2021) for a recent review). Hawking’s famous discovery tells that PBHs with
the mass $\sim 10^{15}$ g would thermally radiate elementary particles Hawking
(1974) and thus they cannot provide all the observed abundance of DM Barrau
_et al._ (2004). The emitted particles such as gamma-rays, neutrinos and
$e^{\pm}$ in the evaporation process suffer from a variety of cosmological
constraints Carr  _et al._ (2020); Laha (2019); Laha  _et al._ (2020); Saha
and Laha (2022); Ray  _et al._ (2021); Dasgupta  _et al._ (2020); Wang  _et
al._ (2021); Calabrese  _et al._ (2021); De Romeri  _et al._ (2021); Ghosh
_et al._ (2021); Capanema  _et al._ (2021); Chao  _et al._ (2021); Bernal  _et
al._ (2022). The elementary DM particles emitted from the present PBHs were
also studied in very recent works Calabrese  _et al._ (2022a); Li and Liao
(2022); Calabrese  _et al._ (2022b). The emitted particles can acquire
energies to a few hundred MeV from the evaporation of PBHs with the mass
around $10^{14}-10^{15}$ g. The light particles then gain enough kinetic
energies to travel through the Earth and reach the terrestrial detectors.

An interesting scenario is that the PBHs play as a novel source to produce
light boosted ALPs Calzà  _et al._ (2021); Schiavone  _et al._ (2021). The
ALPs reach the Earth and interact with the targets in current and future
terrestrial facilities. The absorption of ALPs in the underground detector
targets then induces photoelectron signatures in the time-projection chambers
(TPCs). We refer the light ALPs from the PBH evaporation as “PBH ALP” below.
Unlike the solar axion residing mostly in the keV energy region, PBH ALPs
possess energies up to $\mathcal{O}(100)$ MeV and are boosted to
(semi-)relativistic velocities. As a result, the neutrino scattering
experiments are sensitive to the PBH ALP in light of the energetic
photoelectron signature. Moreover, the solar axions are produced by the
Primakoff process, the “ABC” processes and the nuclear transitions governed by
the axion-photon coupling, axion-electron coupling and axion-nucleon coupling,
respectively. By contrast, the calculation of PBH radiation into elementary
particles are based on quantum mechanics in curved spacetime. The conversion
between the vacuum quantum states close to the black hole and far away from
the black hole leads to a net thermal flux of particles with an emission rate
(see Refs. Hawking (1974); Carr and Hawking (1974); Carr (1975) and recent
reviews Carr  _et al._ (2020); Carr and Kuhnel (2021); Auffinger (2022)).
Thus, the PBH ALP flux is independent of the axion coupling constants. Only
the scattering processes in the terrestrial detectors rely on the couplings
and we thus expect the relevant constraints would be stronger in this sense.

In this work we explore the light ALPs from the PBH evaporation in light of
the neutrino experiments. Super-Kamiokande (Super-K) with 161.9
kiloton$\cdot$years exposure Kachulis  _et al._ (2018) provides a strong
constraint on sub-GeV boosted DM Ema  _et al._ (2019). The typical kinetic
energy of the recoiled electron is above 100 MeV at Super-K. One of us studied
the DM flux from PBH and the consequent DM-electron scattering at Super-K Li
and Liao (2022). Here we will evaluate the ALP scattering induced by the axio-
electric (AE) effect and the inverse Primakoff (IP) process. In addition, a
hypothetical axion-dependent nucleon electric dipole moment (nEDM) Lucente
_et al._ (2022) is introduced and we also examine the relevant PBH ALP event
rate. Utilizing the PBH ALP event rates and the Super-K IV data, we expect to
construct a method to place general exclusion limits on the axion couplings.
The next-generation water Cherenkov detector Hyper-Kamiokande (Hyper-K) Abe
_et al._ (2018), currently under construction, has larger fiducial volume of
the detector and can improve the current sensitivity of T2K. We will also
obtain the prospects of axion couplings for PBH ALP in future Hyper-K
experiment. On the other hand, these limits can also be translated into the
upper bound on the fraction of DM composed of PBHs $f_{\rm PBH}$. We expect
the bound on $f_{\rm PBH}$ can be improved from neutrino experiments, compared
with the current evaporation constraints from extragalactic gamma-rays.

This paper is organized as follows. In Sec. II we evaluate the ALP emission
rate and flux from PBH evaporation. The detection event rate in the
terrestrial facilities are then calculated in Sec. III. We discuss the axio-
electric effect, the inverse Primakoff process as well as a hypothetical
axion-dependent nucleon dipole portal. The resultant constraints on axion
couplings and $f_{\rm PBH}$ from Super-K and the future prospects in Hyper-K
are displayed in Sec. IV. Our conclusions are drawn in Sec. V.

##  II ALP flux from PBH evaporation

The famous discovery by Hawking shows that the PBHs have quantum properties
and thermally radiate with a temperature $T_{\rm PBH}$ given by Hawking
(1975); Page (1976a, 1977); MacGibbon and Webber (1990); MacGibbon (1991)

$\displaystyle k_{\rm B}T_{\rm PBH}={\hbar c^{3}\over 8\pi GM_{\rm
PBH}}\approx 1.06\Big{(}{10^{16}~{}{\rm g}\over M_{\rm PBH}}\Big{)}~{}{\rm
MeV}\;,$ (1)

where $k_{\rm B}$ is the Boltzmann constant, $G$ is the Newtonian constant of
gravitation, and $M_{\rm PBH}$ denotes the PBH mass. New particles of a mass
less than $T_{\rm PBH}$ can be emitted during the PBH evaporation. The PBHs
with $M_{\rm PBH}\lesssim 10^{16}$ g would produce elementary particles
lighter than 1 MeV which can have ultra-relativistic velocities. The general
emission rate can be written as Hawking (1971); Page (1976a, b)

$\frac{d^{2}N}{dEdt}=\frac{g}{2\pi}\frac{\Gamma(E,M_{\rm
PBH},a^{\ast})}{e^{E/T_{\rm PBH}}-(-1)^{2s}}\;,$ (2)

where $E$ is the energy of emitted particle, $a^{\ast}$ is the spin of PBH,
$\Gamma$ is the so-called “graybody” factor describing the probability of
elementary particles escaping the PBH gravitational well, and $g$ and $s$ are
the degrees of freedom and spin of the particle, respectively. Here we
consider the ALP emitted from PBH through the Hawking radiation and thus set
$g=1$, $s=0$. In principle, spinning PBHs evaporate faster than the non-
spinning PBHs and thus contribute to the content of elementary particles
nowadays for $M_{\rm PBH}\gtrsim 7\times 10^{15}~{}{\rm g}$ Dasgupta  _et al._
(2020). We use the public code BlackHawk v2.0 Arbey and Auffinger (2019, 2021)
to calculate the above number of ALP per units time and energy emitted by PBHs
Schiavone  _et al._ (2021). 111BlackHawk v2.0 provides the additional emission
of massive particle and we directly apply it to our massive ALP case.

The ALP flux emitted from PBH can then be divided into the contributions of
both the PBHs in galactic halo and the extragalactic PBHs

$\frac{d^{2}\varphi_{a}}{dEd\Omega}=\frac{d^{2}\varphi_{a}^{\rm{gal}}}{dEd\Omega}+\frac{d^{2}\varphi_{a}^{\rm{egal}}}{dEd\Omega}\;,$
(3)

where $\varphi_{a}^{\rm gal}$ and $\varphi_{a}^{\rm egal}$ correspond to the
ALP flux from galactic (gal) PBHs and extragalactic (egal) PBHs, respectively,
and $\Omega$ is the solid angle. The differential galactic ALP flux is given
by

$\displaystyle{d^{2}\varphi_{a}^{\rm gal}\over dEd\Omega}={f_{\rm PBH}\over
4\pi M_{\rm PBH}}{d^{2}N\over dEdt}\int{d\Omega_{s}\over 4\pi}\int dl\rho_{\rm
MW}[r(l,\psi)]\;,$ (4)

where $f_{\rm PBH}$ is the fraction of DM composed of PBHs, $\rho_{\rm
MW}[r(l,\psi)]$ denotes the DM density in the Milky Way (MW) halo,
$r(l,\psi)=\sqrt{r_{\odot}^{2}+l^{2}-2lr_{\odot}\cos\psi}$ is the
galactocentric distance with $r_{\odot}$ the solar distance from the galactic
center, $l$ the line-of-sight distance to the PBH and $\psi$ the angle between
these two directions. The maximum limit of distance $l$ is taken to be
$(r_{h}^{2}-r_{\odot}^{2}\sin^{2}\psi)^{1/2}+r_{\odot}\cos\psi$ with the halo
radius $r_{h}=200$ kpc Wang  _et al._ (2021). The angular integration is
defined as $\int d\Omega_{s}=\int_{0}^{2\pi}d\phi\int_{0}^{\pi}d\psi\sin\psi$
with the azimuthal angle $\phi$. We adopt the generalized Navarro-Frenk-White
(NFW) DM profile Navarro  _et al._ (1997)

$\displaystyle\rho_{\rm MW}(r)=\rho_{\odot}\Big{(}{r\over
r_{\odot}}\Big{)}^{-\gamma}\Big{(}{1+r_{\odot}/r_{s}\over
1+r/r_{s}}\Big{)}^{3-\gamma}\;,$ (5)

where $\rho_{\odot}=0.4~{}{\rm GeV}/{\rm cm}^{3}$ is the local DM density
Salucci  _et al._ (2010), $r_{\odot}=8.5$ kpc is the distance between the Sun
and the Galactic center Yuan  _et al._ (2017), $r_{s}=20$ kpc is the radius of
the galactic diffusion disk, and the inner slope of the NFW halo profile is
fixed as $\gamma=1$. For the extragalactic contribution, the differential ALP
flux over the full sky is

$\displaystyle{d^{2}\varphi_{a}^{\rm egal}\over dEd\Omega}={f_{\rm
PBH}\rho_{\rm DM}\over 4\pi M_{\rm PBH}}\int^{t_{\rm max}}_{t_{\rm
min}}dt[1+z(t)]{d^{2}N\over dEdt}\Big{|}_{E_{s}}\;,$ (6)

where the average DM density of the Universe at the present epoch is
$\rho_{\rm DM}=2.35\times 10^{-30}~{}{\rm g}/{\rm cm}^{3}$ determined by
Planck Aghanim  _et al._ (2020), and
$E_{s}=\sqrt{(E^{2}-m_{a}^{2})(1+z(t))^{2}+m_{a}^{2}}$ denotes the energy at
the source which is related to the energy $E$ in the observer’s frame by the
redshift $z(t)$ for massive ALP. In the massless limit, it simplifies to
$E_{s}=E(1+z(t))$. The ALP emitted from the PBHs in the very early Universe is
sufficiently redshifted. As a result, their energies are very low today and
cannot be detected. Thus, for the integral limits, we fix $t_{\rm
min}=10^{11}$ s being close to the era of matter-radiation equality Wang  _et
al._ (2021) and $t_{\rm max}$ as the age of Universe. In fact, it turns out
that changing the lower limit has no impact on the results De Romeri  _et al._
(2021).

In the left panel of Fig. 1, we show the differential ALP flux from PBHs
$d^{2}\varphi_{a}/dEd\Omega$ for $m_{a}=1$ keV, $M_{\rm PBH}=10^{15}$ g and
$f_{\rm PBH}=3.9\times 10^{-7}$ allowed by extragalactic gammay-ray
constraints Carr  _et al._ (2020). The galactic and extragalactic
contributions are denoted by solid and dashed lines, respectively. One can see
that the ALP produced from PBH gains energies peaked at order of
$\mathcal{O}(10)-\mathcal{O}(100)$ MeV. The light ALPs would thus receive
ultra-relativistic velocities and are highly boosted. As seen in Eqs. (4) and
(6), smaller PBHs would exhibit harder spectra of the emitted ALP. We also
demonstrate that spinning PBHs make the ALPs gain higher energies along with
increasing PBH spin, because the black hole with non-zero angular moment can
accelerate the Hawking evaporation. In the following study, as being customary
in the literature Wang  _et al._ (2021); Calabrese  _et al._ (2021, 2022a), we
focus on $M_{\rm PBH}<1\times 10^{16}$ g and only consider non-spinning PBHs.
To compare with the ALP, we also display the neutrino flux from PBH
evaporation. There are two contributions to the neutrino emission rate
generated by the Hawking radiation Dasgupta  _et al._ (2020); Wang  _et al._
(2021); Calabrese  _et al._ (2021); De Romeri  _et al._ (2021); Ghosh  _et
al._ (2021); Capanema  _et al._ (2021); Chao  _et al._ (2021); Bernal  _et
al._ (2022). One of them is the primary contribution consisting of neutrinos
directly emitted in the evaporation. The other one is the secondary
contribution which origins from the hadronization and the subsequent decay of
the primary particles. It turns out that the inclusion of secondary production
injects more neutrinos at low energies.

Figure 1: Left: The ALP flux from PBH evaporation, for different PBH spins
(spinless: black line, $a^{\ast}=0.5$: red line and $a^{\ast}=0.9$: blue line)
with fixed ALP mass $m_{a}=1$ keV, PBH mass $M_{\rm{PBH}}=10^{15}$ g and
$f_{\rm PBH}=3.9\times 10^{-7}$. The galactic and extragalactic contributions
are denoted by solid and dashed lines, respectively. The neutrino flux is also
shown for comparison (suppose spinless PBH with the same mass and $f_{\rm
PBH}$). Right: The flux of solar axion, which consists of Primakoff process
(dashed) and ABC production (solid) with the coupling constants
$g_{a\gamma}=10^{-11}~{}\rm{GeV}^{-1}$ and $g_{ae}=10^{-13}$. We set three
finite values of axion mass: $m_{a}=$ 1 keV (red), 5 keV (orange) and 10 keV
(brown).

To further demonstrate the properties of ALPs from PBH evaporation, we need to
compare with the conventional source of axions produced in the Sun. The
components of solar axion mainly originate from the following two mechanisms:
(i) the Primakoff production process (i.e. $\gamma\mathbb{Z}\to a\mathbb{Z}$
with $\mathbb{Z}$ denoting the atomic system $\mathbb{Z}=(e^{-},N)$) which
depends on the axion-photon coupling $g_{a\gamma}$, and (ii) the so-called
“ABC” process Redondo (2013). The “ABC” production consists of three
processes: the Atomic de-excitation and recombination, Bremsstrahlung, and
Compton scattering process. These production channels are all relevant with
axion-electron coupling $g_{ae}$ via the axion-lepton interaction term. Thus,
the relevant interaction terms of solar axion are

$\mathcal{L}\supset-\frac{1}{4}g_{a\gamma}aF_{\mu\nu}\tilde{F}^{\mu\nu}+ig_{ae}a\bar{\psi}_{e}\gamma_{5}\psi_{e}\;.$
(7)

In previous studies of solar axion flux, people usually utilized the property
of axion being a nearly massless pseudo-scalar boson whose mass is far less
than $\mathcal{O}(1)$ eV. Thus, the calculations of solar axion flux can be
simplified under the ultra-relativistic approximation. Nevertheless, computing
the flux of ALP with finite mass requires to generally extend the method. In
more detail, we should recompute the matrix element and phase space integral
without the approximation of massless axion. In fact, some computations of
solar axion flux depending on arbitrary $m_{a}$ have been done in previous
works Raffelt (1986); Gondolo and Raffelt (2009); Caputo  _et al._ (2020);
Carenza and Lucente (2021). We follow them to explicitly calculate the
primakoff, bremsstrahlung and compton scattering process. For the flux
produced in axio-deexcitation and axio-recombination (i.e. free-bound and
bound-bound electron transitions), we simply use the data in Ref. Redondo
(2013) and interpolate them to energies above 10 keV. We also multiply every
flux component (even though it is not rigorous for “A” flux) by a rescaling
factor $(1-m_{a}^{2}/E_{a}^{2})$ to compensate the variation from phase space
integrals of energy-loss rate for non-zero $m_{a}$ case.

The solar ALP flux is shown in the right panel of Fig. 1 for different choices
of $m_{a}=1,5,10$ keV. It turns out that the PBH ALP flux is nearly 16 orders
of magnitude smaller than the solar axion with the coupling constants
$g_{a\gamma}=10^{-11}~{}{\rm GeV}^{-1}$ and $g_{ae}=10^{-13}$ favored by the
stellar cooling Giannotti  _et al._ (2017). However, the spectrum of solar
axion drops rapidly after the peak right above the axion mass threshold for
$m_{a}\sim\mathcal{O}(1)$ keV. Thus, the solar axion is usually used to
account for the measurements from DM direct detection experiments with low
recoil energy, such as the former XENON1T excess in the range of 1-5 keV
Aprile  _et al._ (2020). By contrast, the ALP flux from PBH gets enhanced
along with increasing energies and peaks well above the threshold around
$\mathcal{O}(10)-\mathcal{O}(100)$ MeV. In spite of much low flux, the PBHs
emit more energetic ALPs. The neutrino experiments such as Super-K and Hyper-K
are sensitive to MeV-scale energy regime and can provide compelling bounds for
ALP from PBH evaporation. We will discuss the event rate in these neutrino
experiments and the bounds on axion couplings and PBHs in the following
sections.

##  III Event rates at terrestrial detectors

After being captured by detector, there are two main types of scattering
between the target atom and the ALP: (i) axio-electric effect, referring to
the process that ALP is absorbed by bound-state electron of target atom and
leads to an ionization signal, which depends on the axion-electron coupling
$g_{ae}$, and (ii) inverse Primakoff process
$a+\mathbb{Z}\to\gamma+\mathbb{Z}$. In the latter case (ii), one ALP can be
converted into two photons through the axion-photon coupling $g_{a\gamma}$.
Then, one virtual photon makes atom gain energy (or be excited) while the
other is emitted in the form of real photon. In addition, we will also
consider the detection channel $a+p\to p+\gamma$ on water Cherenkov detector
through the nucleon dipole portal Lucente  _et al._ (2022), which also
produces a visible photon signal. Next, we individually discuss the scattering
cross-sections and event rates corresponding to these types of process.

###  III.1 The axio-electric effect

For the AE effect, the total cross-section is given by Pospelov  _et al._
(2008); Bloch  _et al._ (2017)

$\sigma_{AE}(E_{a})=\sigma_{PE}(E_{a})\frac{3g_{ae}^{2}}{16\pi\alpha_{EM}v_{a}}\frac{E_{a}^{2}}{m_{e}^{2}}(1-\frac{1}{3}v_{a}^{2/3})\;,$
(8)

where $\alpha_{EM}$ is the fine-structure constant, $v_{a}$ is the axion
velocity and $\sigma_{PE}$ is the photoabsorption cross-section of atoms
measured by low-energy x-ray experiment Henke  _et al._ (1993). The cross-
section gives a narrow jagged peak which helps to separate signal from
background as mentioned in Ref. Bloch  _et al._ (2017). In Fig. 2, we show
$\sigma_{AE}$ for ${}^{131}{\rm Xe}$ with the axion-electric coupling
$g_{ae}=10^{-13}$. The ALP absorption cross-section $\sigma_{AE}$ is limited
within the energy range of $0.01\sim 30$ keV because $\sigma_{PE}$ is
suppressed by high photon energies. Here we assume that the incoming velocity
of axion should be relativistic, as the axion radiated from PBH has a velocity
near the speed of light, i.e., $v_{a}\approx c$. Taking into account the ALP
flux $d\varphi_{a}/dE_{a}$ as well as the cross-section $\sigma_{AE}$, we can
further obtain the differential event rate of the scattering process

$\frac{dR}{dE_{R}}={N_{A}\over
m_{A}}\frac{d\varphi_{a}}{dE_{a}}\sigma_{AE}(E_{a})\;,$ (9)

where $N_{A}$ is the Avogadro constant and $m_{A}$ is the molar mass of atom
or molecule. Since the axion mass is much smaller than that of electron or
nucleus, the recoil energy delivered by axion is approximately equal to the
incoming axion energy $E_{a}$. Thus, we can directly rewrite the right-handed
side of above formula which depends on $E_{a}$ into the form relying on the
recoil energy $E_{R}$. Furthermore, to compare with the results in Refs.
Aprile  _et al._ (2020); Gao  _et al._ (2020); Dent  _et al._ (2020); Bloch
_et al._ (2021), the above event rate $dR/dE_{R}$ is smeared by a Gaussian
distribution with the variance of energy as $\sigma_{E}~{}[{\rm
keV}]=0.3171~{}[{\rm keV}^{1/2}]\sqrt{E_{R}}+0.0015E_{R}~{}[{\rm keV}]$ Bloch
_et al._ (2021) which depends on the detector energy and energy resolution.

Combining Eq. (8) and Eq. (9), we can get the event rate of AE effect as shown
in Fig. 3. The left (right) panel is for PBH ALP (solar axion). For PBH we
take $M_{\rm PBH}=10^{15}$ g and $f_{\rm PBH}=3.9\times 10^{-7}$ for
illustration. As a common detection channel for solar axion, due to the
sensitive recoil energy below $\mathcal{O}(10)$ keV, AE effect is usually
considered to account for the results from DM direct detection. But for the
PBH ALP that we focus in this paper, the ALP flux is much small in the working
range of AE effect and its contribution is negligible when fitting to the
results from XENON1T/XENONnT Aprile  _et al._ (2020, 2022). In addition, the
AE effect loses viability at high energies above several tens of keV. To
explore the PBH ALP in the energy range above 100 MeV, we quest for some other
detection channels viable at higher energy below.

Figure 2: Cross-sections for different detection processes including AE effect
for ${}^{131}{\rm Xe}$ (red), inverse Primakoff for xenon (blue), inverse
Primakoff for oxygen (green), and nucleon EDM for free proton (purple). The
axion mass is assumed to 1 keV for inverse Primakoff process (only inverse
Primakoff process dependent on axion mass). The axion coupling constants are
set as $g_{a\gamma}=10^{-11}~{}\rm{GeV}^{-1}$, $g_{ae}=10^{-13}$ favored by
the stellar cooling Giannotti  _et al._ (2017); Di Luzio  _et al._ (2020b) and
$g_{d}=10^{-11}~{}\rm{GeV}^{-2}$, respectively.

Figure 3: Left: Event rates of PBH ALP from AE absorption for ${}^{131}{\rm
Xe}$ (red), inverse Primakoff for ${}^{131}{\rm Xe}$ (blue), inverse Primakoff
for ${\rm H_{2}O}$ (green) and the nucleon EDM for free proton (purple). Note
that the event rate from nEDM have been scaled by a factor $10^{-8}$ to make
them visible. The coupling constants are set as
$g_{a\gamma}=10^{-11}~{}\rm{GeV}^{-1}$, $g_{ae}=10^{-13}$ and
$g_{d}=10^{-11}~{}\rm{GeV}^{-2}$, respectively. We take $M_{\rm PBH}=10^{15}$
g, $f_{\rm PBH}=3.9\times 10^{-7}$ and $m_{a}=1$ keV. Right: Event rates of
solar axion from AE absorption (red) and inverse Primakoff (blue) for
${}^{131}{\rm Xe}$, in which the solar axion flux is obtained by summing over
ABC flux and Primakoff flux. Three different axion masses are assumed as
$m_{a}=1$ keV (solid), 5 keV (dashed) and 10 keV (dotted).

###  III.2 The inverse Primakoff process

In this subsection, we consider the inverse Primakoff process
$a+\mathbb{Z}\to\gamma+\mathbb{Z}$. The outgoing photon would mimic electronic
signal as the current detectors are not able to distinguish the signal caused
by scattered photons from the electron recoil. Its contribution was taken into
account to accommodate the electron recoil events in XENON1T/nT Gao  _et al._
(2020); Dent  _et al._ (2020); Aprile  _et al._ (2022).

Since the typical mass of the incoming ALP we consider is less than
$\sim\rm{MeV}$, the energy transfer is not enough to excite an atom. Thus, we
assume that the atom target does not change its state after being scattered,
i.e., only the elastic inverse Primakoff process is considered. The
differential elastic cross-section per solid angle for $m_{a}\ll E_{a}\approx
E_{\gamma}$ is given by Abe  _et al._ (2021)

$\frac{d\sigma_{\rm{el}}}{d\Omega}\simeq\frac{\alpha_{EM}g_{a\gamma}^{2}\sin\theta}{16\pi(1-\cos\theta)^{2}}|Z-F(\bm{q})|^{2}\;,$
(10)

where $Z$ is the atomic number, and $\theta$ is the angle between the momentum
of axion and outgoing photon. The atomic form factor $F(\bm{q})$ contains the
information about screened effect from electric charge of nucleus. In the
previous studies, one usually used the cross-section induced by the screened
Coulomb potential: $F(q,r_{0})=Z/(1+q^{2}r_{0}^{2})$ Abe  _et al._ (2021),
where $r_{0}$ is the screening length. It gives the cross-section of inverse
Primakoff as follows Creswick  _et al._ (1998); Dent  _et al._ (2020)

$\sigma_{IP}=\frac{\alpha_{EM}Z^{2}g_{a\gamma}^{2}}{2}\left[\frac{2r_{0}^{2}E_{a}^{2}+1}{4r_{0}^{2}E_{a}^{2}}\ln{(1+4r_{0}^{2}E_{a}^{2})}-1\right]\;,$
(11)

where $E_{a}\approx p_{a}$ in the massless limit of axion. However, the above
formula is only valid for ideal atom. In fact, the most commonly used method
in numerical calculation is relativistic-Hartree-Fock (RHF) approximation.
Here we take the total cross-section from Ref. Abe  _et al._ (2021)

$\sigma_{IP}^{\rm
RHF}=\frac{\alpha_{EM}g_{a\gamma}^{2}}{8}\int_{E_{a}-p_{a}}^{E_{a}+p_{a}}{dq\frac{[q^{2}-(E_{a}-p_{a})^{2}][(E_{a}+p_{a})^{2}-q^{2}]}{p_{a}^{2}q^{3}}|Z-F(q)|^{2}}\;,$
(12)

where $F(q)$ is the atomic form factor approximated in a simple analytic form,
which agrees well with the RHF approximation and improves the efficiency of
calculation. Apparently, this cross-section is quadratically dependent on
$g_{a\gamma}$. We display the cross-sections for xenon and oxygen atoms in
Fig. 2 with the coupling $g_{a\gamma}=10^{-11}~{}\rm{GeV}^{-1}$. It turns out
that the inverse Primakoff cross-sections get enhanced for large energies of
axion. Replacing $\sigma_{AE}$ in Eq. (9) by $\sigma_{IP}^{\rm RHF}$, one can
obtain the event rate of inverse Primakoff process as shown in the left panel
of Fig. 3. We can see that the event rate of inverse Primakoff from PBH ALP
(blue and green lines) extends to energies as large as 100 MeV and is 1-2
orders of magnitude greater than that of AE absorption due to the energy
distribution of ALP flux from PBH evaporation.

For solar axion, we show the event rate of both AE effect and inverse
Primakoff process for xenon in the right panel of Fig. 3. The energy range of
event rate is limited within $1\sim 20$ keV due to the profile of solar axion
flux. The event rate for these two detection channels decreases as the axion
mass increases. However, for the PBH ALP in the left panel of Fig. 3, we find
that there is obvious separation between these two kinds of processes. As we
mentioned before, AE effect mainly occurs in low-energy region while at high
energies the inverse Primakoff process dominates and starts to drop around
$\mathcal{O}(10^{4})$ keV. There is a small overlap for the two processes only
at the energy of $\mathcal{O}(10)$ keV. It is thus not difficult to
distinguish the signals caused by these two types of scatterings. Moreover, it
should be emphasized that the signal of PBH ALP is less sensitive to the given
coupling constants of axion because its flux is independent of the couplings.
By contrast, both flux and cross-section of solar axion are proportional to
the coupling squared.

###  III.3 The nucleon dipole portal

Finally, we consider a dynamic interaction for axion inducing an analogous
signal. The so-called axion-dependent nucleon electric dipole moment (EDM) is
given by the following Lagrangian Lucente  _et al._ (2022)

$\mathcal{L}_{a}^{\rm{nEDM}}=-\frac{i}{2}g_{d,N}a\bar{N}\gamma_{5}\sigma_{\mu\nu}NF^{\mu\nu}\;,$
(13)

where $N$ denotes proton or neutron. There is a model-independent feature for
coupling $g_{d,N}$ similar to QCD axion, i.e., $g_{d}\equiv g_{d,n}=-g_{d,p}$.
Consequently, we have a new detection channel, that is $a+p\to p+\gamma$ which
implies that an axion can be absorbed by free proton and then emits a free
photon. This reaction produces a visible flux of photons within detector. This
hypothesis was used to predict the signal of axions emitted from supernova in
free-proton-rich water Cherenkov detector such as Hyper-K Lucente  _et al._
(2022). After integrating the phase space, the total cross-section of
scattering process is given by the neat formula

$\displaystyle\sigma_{\rm{nEDM}}=g_{d}^{2}E_{a}^{2}/2\pi$ (14)

under the limit of small axion mass, which does not depend on axion mass as
shown in Fig. 2. Considering PBH as the source of incoming axion, the total
event rate will be only proportional to $g_{d}^{2}$. In the left panel of Fig.
3 (green line), we show the event rate of nEDM induced scattering for free
proton. The scattered target in this case is the hydrogen atom in water
molecule in which we take into account two free protons. One can see that the
energy range sensitive to nEDM overlaps that of inverse Primakoff for
$E_{a}>1$ MeV and peaks around several tens of MeV. Thus, this detection
channel can also probe the PBH ALPs with energies around 100 MeV.

###  III.4 Relevant constraints

####  III.4.1 The ALP lifetime

Given the decay of ALP to diphoton, the lifetime of ALP becomes

$\displaystyle\tau_{a}={64\pi\over m_{a}^{3}g_{a\gamma}^{2}}\;.$ (15)

We require the ALP lifetime to be longer than the era of matter-radiation
equality, i.e. $10^{11}~{}{\rm s}\sim 1~{}{\rm kpc}$, to give a conservative
constraint on $g_{a\gamma}$. Note that taking into account the relativistic
Lorentz factor $\gamma\simeq E_{a}/m_{a}$ would make the constraint much
weaker when $m_{a}\ll E_{a}$. Given $E_{a}\simeq 10$ MeV and $m_{a}=1$ keV for
illustration, the limit on $g_{a\gamma}$ will be weakened by two orders of
magnitude.

####  III.4.2 The attenuation of ALPs

The PBH ALPs can travel a distance in the atmosphere and underground before
reaching the terrestrial detectors and would be stopped due to the process
$a+p\to p+\gamma$ therein. Their path $\ell$ is thus required to satisfy
Lucente  _et al._ (2022)

$\displaystyle\ell=\Big{(}n_{N}\sigma_{\rm nEDM}\Big{)}^{-1}>d^{\oplus}~{}{\rm
or}~{}d^{\rm atm}\;,$ (16)

where $n_{N}$ is the nucleon density, $d^{\oplus}\simeq 1$ km is the depth of
the detector from the Earth surface and $d^{\rm atm}=86$ km is the altitude of
atmosphere in the US Standard Atmosphere model Emken and Kouvaris (2018).
Taking a constant mass density of Earth crust $\rho^{\oplus}=2.7~{}{\rm
g}~{}{\rm cm}^{-3}$ Emken and Kouvaris (2018) and running over the most
abundant elements in the Earth, we obtain $n_{N}^{\oplus}=1.6\times
10^{24}~{}{\rm cm}^{-3}$. The atmosphere was modeled as 4 layers of falling
mass density as a function of the altitude. We take a constant density
$\rho^{\rm atm}=1.1~{}{\rm kg}~{}{\rm m}^{-3}$ in the first layer Emken and
Kouvaris (2018) and the most abundant elements in the atmosphere (nitrogen and
oxygen) to obtain the number density of nucleons as $n_{N}^{\rm atm}=6.5\times
10^{20}~{}{\rm cm}^{-3}$. Taking $E_{a}>100$ MeV as a reasonable spectrum
limit, we find the upper limit of $g_{d}$ as $3~{}(17)~{}{\rm GeV}^{-2}$ from
the attenuation underground (in the atmosphere) and thus the constraint due to
the attenuation is very weak.

####  III.4.3 The supernova bounds and others

The supernova is also a novel source of ALPs and the relevant observations of
supernova (SN) permit to place constraints on the ALP couplings. There exists
study of supernova population with low explosion energies placing the
constraint on $a\to\gamma\gamma$ and the coupling $g_{a\gamma}$ Caputo  _et
al._ (2022a). The observation of the SN 1987A neutrino burst also allows to
constrain the axion coupling $g_{d}$ by evaluating the axion luminosity
Lucente  _et al._ (2022). There are also other constraints on $g_{a\gamma}$
from visible decays of ALPs produced in SN 1987A and beam dum experiments. We
place the relevant bounds in the figures below.

##  IV Detection and constraints in neutrino experiments

Next we examine the detection of PBH ALP events and the relevant constraints
in neutrino experiments. The event number per recoil energy is given by

$\displaystyle{dN\over dE_{R}}=M_{\rm det}t{dR\over dE_{R}}\;,$ (17)

where $M_{\rm det}$ is the fiducial mass of the detector and $t$ is the
exposure time. As the neutrino experiments are only sensitive to visible
energy above 100 MeV, we ignore the negligible AE absorption in $dR/dE_{R}$
and only consider the inverse Primakoff process for ALP coupling $g_{a\gamma}$
with water molecule or for nEDM hypothesis with free protons.

The Super-K detector performed 2628.1 days of running in a single tank
consisting of 22.5 kiloton (kton) fiducial volume pure water. This corresponds
to an exposure of 162 kiloton$\cdot$years. In the Super-K analysis, there are
three energy bins with $100~{}{\rm MeV}<E<1.33~{}{\rm GeV}$, $1.33~{}{\rm
GeV}<E<20~{}{\rm GeV}$ and $E>20~{}{\rm GeV}$ for the “electron elastic
scatter-like” events Kachulis  _et al._ (2018). The numbers of measured event
and the estimated background for each energy bin are listed in Table I of Ref.
Kachulis  _et al._ (2018). To set bounds on the PBH ALP parameter space from
Super-K, we define the $\chi^{2}$ with regard to the three energy bins

$\displaystyle\chi^{2}=\sum_{i=1}^{3}\min_{\alpha_{i}}\left[2\left(N_{\text{th}}^{i}-N_{\text{obs}}^{i}+N_{\text{obs}}^{i}\ln\frac{N_{\text{obs}}^{i}}{N_{\text{th}}^{i}}\right)+\left(\frac{\alpha_{i}}{\sigma_{i}}\right)^{2}\right]\;,$
(18)

where $N_{\text{th}}^{i}=N_{\text{bkg}}^{i}(1+\alpha_{i})+N_{\text{sig}}^{i}$.
Here $N_{\text{bkg}}^{i}$ ($N_{\text{obs}}^{i}$) [$N_{\text{sig}}^{i}$] is the
number of background (observed) [signal] events in each bin, $\alpha_{i}$ is
the background normalization factor that is minimized over, and $\sigma_{i}$
is the uncertainty for each bin. Follow Ref. Kachulis  _et al._ (2018), we
take a conservative value of $19\%$ for the uncertainties in the two bins with
lower energies, and 23% for the uncertainty in the third energy bin. When
calculating our signal event number, we also convolute the signal efficiency
taken from the Fig. 1 in Ref. Kachulis  _et al._ (2018). As the event rate of
inverse Primakoff decreases rapidly for energies above $\sim 100$ MeV due to
the ALP flux from PBH, only the signal event for the first energy bin in fact
matters for our analysis. The background-only hypothesis yields a
$\chi^{2}_{0}=2.7$, and the 2$\sigma$ exclusion bound is set with
$\Delta\chi^{2}=\chi^{2}-\chi^{2}_{0}=4.0$ for one degree of freedom Li and
Liao (2022).

For the Hyper-K with underground water Cherenkov detector, the fiducial mass
is 374 kton and we refer the exposure as 1870 kton$\cdot$years Abe  _et al._
(2018) with an ideal detection efficiency of signal as $\epsilon=1$. Following
Ref. Lucente  _et al._ (2022), we place a conservative bound on PBH ALP from
Hyper-K as

$\epsilon\times N_{s}<N_{\rm HK}\;,$ (19)

where $N_{s}$ is the signal event number for $E>100$ MeV and $N_{\rm HK}=2$ is
the assumed number of observation.

We first show the bounds on the ALP couplings versus $m_{a}$ from Super-K
(solid line) and Hyper-K (dashed line) in Fig. 4 (left panel for $g_{a\gamma}$
and right panel for $g_{d}$) for different choices of PBH masses with
corresponding maximum $f_{\rm PBH}$ allowed by the extragalactic gamma-ray
constraints. It turns out that Super-K (Hyper-K) can confine the axion-photon
coupling $g_{a\gamma}$ down to $8\times 10^{-4}~{}{\rm GeV}^{-1}$ ($8\times
10^{-6}~{}{\rm GeV}^{-1}$), $3\times 10^{-3}~{}{\rm GeV}^{-1}$ ($3\times
10^{-5}~{}{\rm GeV}^{-1}$) and $6\times 10^{-6}~{}{\rm GeV}^{-1}$ ($6\times
10^{-8}~{}{\rm GeV}^{-1}$) for $M_{\rm PBH}=1.0\times 10^{15},~{}0.5\times
10^{15},~{}0.1\times 10^{15}$ g, respectively. We also display the constraints
from ALP lifetime, supernova 1987A, the cooling of horizontal branch (HB)
stars as well as beam dump experiments. It turns out that the above sensitive
region of $g_{a\gamma}$ from PBH ALP has already been excluded by these
constraints. For axion-nEDM coupling $g_{d}$, the corresponding limits become
$2\times 10^{-8}~{}{\rm GeV}^{-2}$ ($2\times 10^{-10}~{}{\rm GeV}^{-2}$),
$5\times 10^{-8}~{}{\rm GeV}^{-2}$ ($5\times 10^{-10}~{}{\rm GeV}^{-2}$) and
$6\times 10^{-11}~{}{\rm GeV}^{-2}$ ($6\times 10^{-13}~{}{\rm GeV}^{-2}$) from
Super-K (Hyper-K).

These bounds can be translated into the upper limit on $f_{\rm PBH}\cdot
g^{2}_{a\gamma(d)}$ as a function of $M_{\rm PBH}$, as shown in Fig. 5 for
$f_{\rm PBH}\cdot g^{2}_{d}$ and $m_{a}=1$ keV for illustration. The combined
excluded region of $f_{\rm PBH}\cdot g^{2}_{d}$ from extragalactic gamma-rays
for $f_{\rm PBH}$ Carr  _et al._ (2020) and SN 1987A for $g_{d}$ Lucente  _et
al._ (2022) is shown in gray. One can see that the neutrino experiments can
provide viable search for the axion-nEDM coupling of PBH ALP and are more
sensitive to low PBH masses due to the visible energy range. As shown above,
$g_{a\gamma}$ is highly constrained and we do not show the limit on $f_{\rm
PBH}\cdot g^{2}_{a\gamma}$ here.

Figure 4: The bounds from Super-K (solid) and projected sensitivity at Hyper-K
(dashed) on $g_{a\gamma}$ (left) and $g_{d}$ (right) as a function of $m_{a}$.
The PBH masses are assumed as $M_{\rm PBH}=1.0\times 10^{15}$ g (black),
$0.5\times 10^{15}$ g (red) and $0.1\times 10^{15}$ g (blue) with
corresponding maximum $f_{\rm PBH}$ allowed by the extragalactic gamma-ray
constraints. We also present other constraints on $g_{a\gamma}$ from low-
energy supernovae radiation Caputo  _et al._ (2022a) (pink), beam dump
experiments Dolan  _et al._ (2017) (purple), SN 1987A gamma-rays Caputo  _et
al._ (2022b); Jaeckel  _et al._ (2018) (sky blue), LEP (yellow green),
horizontal branch stars Jaeckel  _et al._ (2018) (orange), SN 1987A neutrinos
Caputo  _et al._ (2022b, a) (teal) and ALP lifetime (green). The upper limit
of $g_{d}=6.7\times 10^{-9}~{}\rm{GeV}^{-2}$ Lucente  _et al._ (2022) obtained
by neutrino from SN1987A absorbed by proton is also shown in gray.  Figure 5:
The bounds on $f_{\rm PBH}\cdot g_{d}^{2}$ as a function of $M_{\rm PBH}$ from
Super-K (solid) and Hyper-K (dashed) for $m_{a}=1$ keV as an illustration. The
current evaporation constraint on $f_{\rm PBH}$ from extragalactic gamma-rays
Carr  _et al._ (2020) and the SN 1987A bound on $g_{d}$ Lucente  _et al._
(2022) are combined and the excluded region is shown in gray for comparison.

##  V Conclusion

We explore the possibility that the light axion-like particles are produced
from PBH evaporation. The ALPs receive energies in the region of a few hundred
MeV. The absorption of light ALPs induces energetic photoelectron signatures
in the underground detector targets of neutrino experiments. Given the
galactic and extragalactic PBH ALP flux, we calculate the PBH ALP event rate
produced by the inverse Primakoff process with axion-photon coupling
$g_{a\gamma}$ and with axion-nucleon coupling $g_{d}$ for the hypothesis of
axion-dependent nucleon electric dipole moment. We then place general
exclusion limits on the axion couplings utilizing the Super-K data and predict
the prospects for the future Hyper-K. We also translate these limits into the
upper bound on the fraction of DM composed of PBHs $f_{\rm PBH}$.

We find the following conclusions.

  * •

Super-K can constrain the axion-photon coupling $g_{a\gamma}$ as small as
$8\times 10^{-4}~{}{\rm GeV}^{-1}$, $3\times 10^{-3}~{}{\rm GeV}^{-1}$ and
$6\times 10^{-6}~{}{\rm GeV}^{-1}$ for $M_{\rm PBH}=1.0\times
10^{15},~{}0.5\times 10^{15},~{}0.1\times 10^{15}$ g, respectively. The above
sensitive region of $g_{a\gamma}$ from PBH ALP has already been excluded by
other observations and constraints. For axion-nEDM coupling $g_{d}$, the
corresponding limits become $2\times 10^{-8}~{}{\rm GeV}^{-2}$, $5\times
10^{-8}~{}{\rm GeV}^{-2}$ and $6\times 10^{-11}~{}{\rm GeV}^{-2}$ from
Super-K.

  * •

Hyper-K can improve the limits of axion coupling constants by two orders of
magnitude.

  * •

The neutrino experiments provide improved bounds on $f_{\rm PBH}\cdot
g_{d}^{2}$ for small PBH masses, compared with the constraints on $f_{\rm
PBH}$ from extragalactic gamma-rays and $g_{d}$ from SN 1987A.

Note that these bounds on axion couplings and $f_{\rm PBH}$ are obtained
simultaneously and are fully correlated.

Note Added: During the completion of this work, a study Cui  _et al._ (2022)
appeared and investigated the boosted ALPs from dark matter decay in neutrino
experiments.

## ACKNOWLEDGMENTS

We would like to thank Andrea Caputo and Wei Chao for useful communication.
T.L. is supported by the National Natural Science Foundation of China (Grant
No. 11975129, 12035008) and “the Fundamental Research Funds for the Central
Universities”, Nankai University (Grants No. 63196013).

## References

  * Aalbers  _et al._ (2022) J. Aalbers _et al._ (LUX-ZEPLIN), (2022), arXiv:2207.03764 [hep-ex] .
  * Aprile  _et al._ (2022) E. Aprile _et al._ , (2022), arXiv:2207.11330 [hep-ex] .
  * Baluni (1979) V. Baluni, Phys. Rev. D 19, 2227 (1979).
  * Crewther  _et al._ (1979) R. J. Crewther, P. Di Vecchia, G. Veneziano, and E. Witten, Phys. Lett. B 88, 123 (1979), [Erratum: Phys.Lett.B 91, 487 (1980)].
  * Kim (1979) J. E. Kim, Phys. Rev. Lett. 43, 103 (1979).
  * Shifman  _et al._ (1980) M. A. Shifman, A. I. Vainshtein, and V. I. Zakharov, Nucl. Phys. B 166, 493 (1980).
  * Dine  _et al._ (1981) M. Dine, W. Fischler, and M. Srednicki, Phys. Lett. B 104, 199 (1981).
  * Zhitnitsky (1980) A. R. Zhitnitsky, Sov. J. Nucl. Phys. 31, 260 (1980).
  * Baker  _et al._ (2006) C. A. Baker _et al._ , Phys. Rev. Lett. 97, 131801 (2006), arXiv:hep-ex/0602020 .
  * Pendlebury  _et al._ (2015) J. M. Pendlebury _et al._ , Phys. Rev. D 92, 092003 (2015), arXiv:1509.04411 [hep-ex] .
  * Peccei and Quinn (1977a) R. D. Peccei and H. R. Quinn, Phys. Rev. Lett. 38, 1440 (1977a).
  * Peccei and Quinn (1977b) R. D. Peccei and H. R. Quinn, Phys. Rev. D 16, 1791 (1977b).
  * Weinberg (1978) S. Weinberg, Phys. Rev. Lett. 40, 223 (1978).
  * Wilczek (1978) F. Wilczek, Phys. Rev. Lett. 40, 279 (1978).
  * Di Luzio  _et al._ (2020a) L. Di Luzio, M. Giannotti, E. Nardi, and L. Visinelli, Phys. Rept. 870, 1 (2020a), arXiv:2003.01100 [hep-ph] .
  * Kim (1987) J. E. Kim, Phys. Rept. 150, 1 (1987).
  * Kuster  _et al._ (2008) M. Kuster, G. Raffelt, and B. Beltran, eds.,  _Axions: Theory, cosmology, and experimental searches. Proceedings, 1st Joint ILIAS-CERN-CAST axion training, Geneva, Switzerland, November 30-December 2, 2005_ , Vol. 741 (2008).
  * Zel’dovich (1967) Y. B. . N. Zel’dovich, I. D., Soviet Astron. AJ (Engl. Transl. ), 10, 602 (1967).
  * Carr and Hawking (1974) B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * Carr (1975) B. J. Carr, Astrophys. J. 201, 1 (1975).
  * Khlopov (2010) M. Y. Khlopov, Res. Astron. Astrophys. 10, 495 (2010), arXiv:0801.0116 [astro-ph] .
  * Carr and Kuhnel (2021) B. Carr and F. Kuhnel, in  _Les Houches summer school on Dark Matter_ (2021) arXiv:2110.02821 [astro-ph.CO] .
  * Hawking (1974) S. W. Hawking, Nature 248, 30 (1974).
  * Barrau  _et al._ (2004) A. Barrau, D. Blais, G. Boudoul, and D. Polarski, Annalen Phys. 13, 115 (2004), arXiv:astro-ph/0303330 .
  * Carr  _et al._ (2020) B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, (2020), arXiv:2002.12778 [astro-ph.CO] .
  * Laha (2019) R. Laha, Phys. Rev. Lett. 123, 251101 (2019), arXiv:1906.09994 [astro-ph.HE] .
  * Laha  _et al._ (2020) R. Laha, J. B. Muñoz, and T. R. Slatyer, Phys. Rev. D 101, 123514 (2020), arXiv:2004.00627 [astro-ph.CO] .
  * Saha and Laha (2022) A. K. Saha and R. Laha, Phys. Rev. D 105, 103026 (2022), arXiv:2112.10794 [astro-ph.CO] .
  * Ray  _et al._ (2021) A. Ray, R. Laha, J. B. Muñoz, and R. Caputo, Phys. Rev. D 104, 023516 (2021), arXiv:2102.06714 [astro-ph.CO] .
  * Dasgupta  _et al._ (2020) B. Dasgupta, R. Laha, and A. Ray, Phys. Rev. Lett. 125, 101101 (2020), arXiv:1912.01014 [hep-ph] .
  * Wang  _et al._ (2021) S. Wang, D.-M. Xia, X. Zhang, S. Zhou, and Z. Chang, Phys. Rev. D 103, 043010 (2021), arXiv:2010.16053 [hep-ph] .
  * Calabrese  _et al._ (2021) R. Calabrese, D. F. G. Fiorillo, G. Miele, S. Morisi, and A. Palazzo, (2021), arXiv:2106.02492 [hep-ph] .
  * De Romeri  _et al._ (2021) V. De Romeri, P. Martínez-Miravé, and M. Tórtola, (2021), arXiv:2106.05013 [hep-ph] .
  * Ghosh  _et al._ (2021) D. Ghosh, A. Guha, and D. Sachdeva, (2021), arXiv:2110.00025 [hep-ph] .
  * Capanema  _et al._ (2021) A. Capanema, A. Esmaeili, and A. Esmaili, (2021), arXiv:2110.05637 [hep-ph] .
  * Chao  _et al._ (2021) W. Chao, T. Li, and J. Liao, (2021), arXiv:2108.05608 [hep-ph] .
  * Bernal  _et al._ (2022) N. Bernal, V. Muñoz Albornoz, S. Palomares-Ruiz, and P. Villanueva-Domingo, (2022), arXiv:2203.14979 [hep-ph] .
  * Calabrese  _et al._ (2022a) R. Calabrese, M. Chianese, D. F. G. Fiorillo, and N. Saviano, Phys. Rev. D 105, L021302 (2022a), arXiv:2107.13001 [hep-ph] .
  * Li and Liao (2022) T. Li and J. Liao, Phys. Rev. D 106, 055043 (2022), arXiv:2203.14443 [hep-ph] .
  * Calabrese  _et al._ (2022b) R. Calabrese, M. Chianese, D. F. G. Fiorillo, and N. Saviano, Phys. Rev. D 105, 103024 (2022b), arXiv:2203.17093 [hep-ph] .
  * Calzà  _et al._ (2021) M. Calzà, J. March-Russell, and J. a. G. Rosa, (2021), arXiv:2110.13602 [astro-ph.CO] .
  * Schiavone  _et al._ (2021) F. Schiavone, D. Montanino, A. Mirizzi, and F. Capozzi, JCAP 08, 063 (2021), arXiv:2107.03420 [hep-ph] .
  * Auffinger (2022) J. Auffinger, (2022), arXiv:2206.02672 [astro-ph.CO] .
  * Kachulis  _et al._ (2018) C. Kachulis _et al._ (Super-Kamiokande), Phys. Rev. Lett. 120, 221301 (2018), arXiv:1711.05278 [hep-ex] .
  * Ema  _et al._ (2019) Y. Ema, F. Sala, and R. Sato, Phys. Rev. Lett. 122, 181802 (2019), arXiv:1811.00520 [hep-ph] .
  * Lucente  _et al._ (2022) G. Lucente, L. Mastrototaro, P. Carenza, L. Di Luzio, M. Giannotti, and A. Mirizzi, Phys. Rev. D 105, 123020 (2022), arXiv:2203.15812 [hep-ph] .
  * Abe  _et al._ (2018) K. Abe _et al._ (Hyper-Kamiokande), (2018), arXiv:1805.04163 [physics.ins-det] .
  * Hawking (1975) S. W. Hawking, Commun. Math. Phys. 43, 199 (1975), [Erratum: Commun.Math.Phys. 46, 206 (1976)].
  * Page (1976a) D. N. Page, Phys. Rev. D 13, 198 (1976a).
  * Page (1977) D. N. Page, Phys. Rev. D 16, 2402 (1977).
  * MacGibbon and Webber (1990) J. H. MacGibbon and B. R. Webber, Phys. Rev. D 41, 3052 (1990).
  * MacGibbon (1991) J. H. MacGibbon, Phys. Rev. D 44, 376 (1991).
  * Hawking (1971) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * Page (1976b) D. N. Page, Phys. Rev. D 14, 3260 (1976b).
  * Arbey and Auffinger (2019) A. Arbey and J. Auffinger, Eur. Phys. J. C 79, 693 (2019), arXiv:1905.04268 [gr-qc] .
  * Arbey and Auffinger (2021) A. Arbey and J. Auffinger, Eur. Phys. J. C 81, 910 (2021), arXiv:2108.02737 [gr-qc] .
  * Navarro  _et al._ (1997) J. F. Navarro, C. S. Frenk, and S. D. M. White, Astrophys. J. 490, 493 (1997), arXiv:astro-ph/9611107 .
  * Salucci  _et al._ (2010) P. Salucci, F. Nesti, G. Gentile, and C. F. Martins, Astron. Astrophys. 523, A83 (2010), arXiv:1003.3101 [astro-ph.GA] .
  * Yuan  _et al._ (2017) Q. Yuan, S.-J. Lin, K. Fang, and X.-J. Bi, Phys. Rev. D 95, 083007 (2017), arXiv:1701.06149 [astro-ph.HE] .
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), arXiv:1807.06209 [astro-ph.CO] .
  * Redondo (2013) J. Redondo, JCAP 12, 008 (2013), arXiv:1310.0823 [hep-ph] .
  * Raffelt (1986) G. G. Raffelt, Phys. Rev. D 33, 897 (1986).
  * Gondolo and Raffelt (2009) P. Gondolo and G. G. Raffelt, Phys. Rev. D 79, 107301 (2009), arXiv:0807.2926 [astro-ph] .
  * Caputo  _et al._ (2020) A. Caputo, A. J. Millar, and E. Vitagliano, Phys. Rev. D 101, 123004 (2020), arXiv:2005.00078 [hep-ph] .
  * Carenza and Lucente (2021) P. Carenza and G. Lucente, Phys. Rev. D 103, 123024 (2021), arXiv:2104.09524 [hep-ph] .
  * Giannotti  _et al._ (2017) M. Giannotti, I. G. Irastorza, J. Redondo, A. Ringwald, and K. Saikawa, JCAP 10, 010 (2017), arXiv:1708.02111 [hep-ph] .
  * Aprile  _et al._ (2020) E. Aprile _et al._ (XENON), Phys. Rev. D 102, 072004 (2020), arXiv:2006.09721 [hep-ex] .
  * Pospelov  _et al._ (2008) M. Pospelov, A. Ritz, and M. B. Voloshin, Phys. Rev. D 78, 115012 (2008), arXiv:0807.3279 [hep-ph] .
  * Bloch  _et al._ (2017) I. M. Bloch, R. Essig, K. Tobioka, T. Volansky, and T.-T. Yu, JHEP 06, 087 (2017), arXiv:1608.02123 [hep-ph] .
  * Henke  _et al._ (1993) B. L. Henke, E. M. Gullikson, and J. C. Davis, Atom. Data Nucl. Data Tabl. 54, 181 (1993).
  * Gao  _et al._ (2020) C. Gao, J. Liu, L.-T. Wang, X.-P. Wang, W. Xue, and Y.-M. Zhong, Phys. Rev. Lett. 125, 131806 (2020), arXiv:2006.14598 [hep-ph] .
  * Dent  _et al._ (2020) J. B. Dent, B. Dutta, J. L. Newstead, and A. Thompson, Phys. Rev. Lett. 125, 131805 (2020), arXiv:2006.15118 [hep-ph] .
  * Bloch  _et al._ (2021) I. M. Bloch, A. Caputo, R. Essig, D. Redigolo, M. Sholapurkar, and T. Volansky, JHEP 01, 178 (2021), arXiv:2006.14521 [hep-ph] .
  * Di Luzio  _et al._ (2020b) L. Di Luzio, M. Fedele, M. Giannotti, F. Mescia, and E. Nardi, Phys. Rev. Lett. 125, 131804 (2020b), arXiv:2006.12487 [hep-ph] .
  * Abe  _et al._ (2021) T. Abe, K. Hamaguchi, and N. Nagata, Phys. Lett. B 815, 136174 (2021), arXiv:2012.02508 [hep-ph] .
  * Creswick  _et al._ (1998) R. J. Creswick, F. T. Avignone, III, H. A. Farach, J. I. Collar, A. O. Gattone, S. Nussinov, and K. Zioutas, Phys. Lett. B 427, 235 (1998), arXiv:hep-ph/9708210 .
  * Emken and Kouvaris (2018) T. Emken and C. Kouvaris, Phys. Rev. D 97, 115047 (2018), arXiv:1802.04764 [hep-ph] .
  * Caputo  _et al._ (2022a) A. Caputo, H.-T. Janka, G. Raffelt, and E. Vitagliano, Phys. Rev. Lett. 128, 221103 (2022a), arXiv:2201.09890 [astro-ph.HE] .
  * Dolan  _et al._ (2017) M. J. Dolan, T. Ferber, C. Hearty, F. Kahlhoefer, and K. Schmidt-Hoberg, Journal of High Energy Physics 2017 (2017), 10.1007/jhep12(2017)094.
  * Caputo  _et al._ (2022b) A. Caputo, G. Raffelt, and E. Vitagliano, Physical Review D 105 (2022b), 10.1103/physrevd.105.035022.
  * Jaeckel  _et al._ (2018) J. Jaeckel, P. Malta, and J. Redondo, Physical Review D 98 (2018), 10.1103/physrevd.98.055032.
  * Cui  _et al._ (2022) Y. Cui, J.-L. Kuo, J. Pradler, and Y.-D. Tsai, (2022), arXiv:2207.13107 [hep-ph] .


# Learning Material Parameters and Hydrodynamics of Soft Robotic Fish via
Differentiable Simulation

John Z. Zhang1,2, Yu Zhang1, Pingchuan Ma3, Elvis Nava1,4, Tao Du3, Philip
Arm1,
Wojciech Matusik3, Robert K. Katzschmann1,4 *We are grateful for funding
received by the ETH AI Center, and the Defense Advanced Research Projects
Agency.1Soft Robotics Lab, ETH Zurich, Switzerland2Mechanical Engineering,
MIT, Cambridge, MA, USA3Computer Science and AI Lab, MIT, Cambridge, MA,
USA4ETH AI Center, ETH Zurich<EMAIL_ADDRESS>
###### Abstract

The high dimensionality of soft mechanisms and the complex physics of fluid-
structure interactions render the sim2real gap for soft robots particularly
challenging. Our framework allows high fidelity prediction of dynamic behavior
for composite bi-morph bending structures in real hardware to millimeter-
accuracy. We address this gap with our differentiable simulation tool by
learning the material parameters and hydrodynamics of our robots. We
demonstrate an experimentally-verified, fast optimization pipeline for
learning the material parameters and hydrodynamics from quasi-static and
dynamic data via differentiable simulation. Our method identifies physically
plausible Young’s moduli for various soft silicone elastomers and stiff acetal
copolymers used in creation of our three different fish robot designs and is
compatible with varying internal geometry of the actuators, such as number of
air chambers. For these robots we provide a differentiable and more robust
estimate of the thrust force than analytical models and we successfully
predict deformation to millimeter accuracy in dynamic experiments under
various actuation signals. Although we focus on a specific application for
underwater soft robots, our framework is applicable to any pneumatically
actuated soft mechanism. This work presents a prototypical hardware and
simulation problem solved using our framework that can be extended
straightforwardly to higher dimensional parameter inference, learning control
policies, and computational design enabled by its differentiability.

##  I Introduction

Soft robots are favored over traditional hard robots in a growing list of
scenarios, including bio-inspired design, conformal gripping, co-bot
situations, and more. While the design and simulation of hard robots are well
understood and extensively shown thanks to decades of development, the same
simulation task for compliant structures and their interaction with the
environment is an ongoing challenge for soft roboticists.

Although many simulation schemes exist, none have demonstrated sufficient
sim2real matching through verifiable experiments for soft robots in a
hydrodynamic situation such as a robotic fish. The high dimensionality of soft
mechanisms and the complex physics of fluid-structure interactions (FSI) make
this problem particularly difficult. The problem has recently become more
tractable because of advances in physics-based simulation, such as
differentiability and data-driven approaches for modeling complex phenomena.
These recent advances show promise in efficiently bridging the sim2real gap
from soft body simulation to real scenarios for soft robots. Straightforward
sim2real methods for roboticists as users are within reach: such methods could
be used confidently for design optimization if they are verifiable by real
data.

Figure 1: Left: A model of a pneumatic fish tail (_Nemo_). Middle: The same
fish tail inflated to $250\text{\,}\mathrm{m}\mathrm{b}\mathrm{a}\mathrm{r}$
in simulation. Right: The actual hardware of the fish tail in an experimental
setup for measuring the deformation and the thrust of the fish in air and
underwater. Our framework is capable of reproducing the buckling effect in
simulation, a valuable result.

There are a few notable related works that have addressed parts of the
challenge. Drake is a rigid multi-body simulator that supports auto-
differentiation and allows the user to tune the parameters of a multi-link
hard robot to maximize an objective [1]. Hahn et al. propose a method for
sim2real matching that is limited to an open-loop controlled tendon-driven
crawler [2]. SOFA is an non-differentiable simulator that supports pneumatic
and tendon-driven actuation, but requires extensive computation time and setup
overhead for complicated systems [3]. Additionally, SOFA is not capable of
supporting FSI, and there is a lack of other soft actuation methods.

In this work, we present a simulation tool of general use to soft roboticists
by demonstrating an experimentally-verified simulation framework that can be
used to accurately predict the deformation of a soft pneumatic actuator and
identify its material properties. Specifically, we focus on the case of a
pneumatically-actuated fish tail. The deformable body of the fish tail is cast
out of silicone elastomer and has a flexible spine made of an acetal plastic
sheet in the center (see Figure 1). This prediction and identification problem
is of interest since composite bi-morph structures are difficult to simulate
due to the disparate Young’s moduli and aspect ratios of the deformable body
and flexible but stiff spine.

We demonstrate an experimentally-verified, fast pipeline for automatic system
identification of material parameters using a gradient search on a quasistatic
data set. The fast system identification is achieved using our finite element
method (FEM) [4, 5], which combines projective dynamics [6] and differentiable
simulation [7]. The parameters we identify from quasi-static experiments fall
within the range of traditionally measured Young’s moduli for our tested
materials. In dynamic flapping experiments for various amplitudes and
frequencies of pressure actuation, we successfully predict the deformation to
within millimeter accuracy.

To verify our simulation results, we present an experimental setup suitable
for collecting position, pressure, and force data using marker tracking, a
pneumatic valve array, and a load cell that is suitable for both air and
water. The raw data is synchronized and analyzed in our MATLAB-based software
pipeline to enable data export to a convenient format. In total, we design,
build and characterize three different fish designs for use in our
identification and prediction pipeline (see Table I).

Using the data sets collected underwater, we construct a neural network
surrogate for predicting the hydrodynamic thrust produced by the motion of the
fish tail. We show that our model is able to generalize to actuation signals
not seen in the training set and that it brings advantages over simple
analytical models such as the large amplitude elongated body theory (EBT) [8].
We chose to use a neural network so that our simulation framework remains
fully differentiable.

In this paper, we contribute:

  * •

a simulation framework that solves the multi-physics problem for a soft
pneumatic bending actuator consisting of disparate materials and geometries;

  * •

a system identification method that uses differentiable simulation and
gradient-based optimization to accurately learn material parameters, i.e., the
Young’s moduli of two corotated materials;

  * •

verification of our simulation results using an experimental setup suitable
for data acquisition in both air and water; and

  * •

a data-driven hydrodynamics model using a neural network as a simple predictor
of the thrust force generated by the fish, preserving differentiability of our
simulator and providing greater robustness than analytical models.

##  II Related Work

###  II-A Soft underwater robots

Soft robots are difficult to optimally design and control when compared to
their rigid counterparts due to the infinite dimensionality of their compliant
structures. Due to this modeling complexity, an experienced designer must hand
craft each design guided by intuition, experiments, and approximate models.
Marchese et al. offer approaches to designing and fabricating soft fluidic
elastomer robots, the type of robot we are also using in this work [9].
Katzschmann et al. present the design, fabrication, control, and testing of a
soft robotic fish with interior cavities that is hydraulically actuated. Their
manually designed robot can swim at multiple depths and record aquatic life in
the ocean [10, 11]. Zhu et al. manually optimize the swimming performance of
their robotic fish called Tunabot [12]. The authors measured kinematics,
speed, and power at increasing flapping frequencies to quantify swimming
performance and find agreement in performance between real fish and their
Tunabot over a wide range of frequencies. Zheng et al. propose to design soft
robots by pre-checking controllability during the numerical design phase [13].
FEM is used to model the dynamics of cable-driven parallel soft robot and a
differential geometric method is applied to analyze the controllability of the
points of interest. Katzschmann et al. [14] manually tweak the material
parameters of their reduced-order FEM [15] with an experimental soft robotic
arm to perform dynamic closed-loop control. Van et al. present a DC motor
driven soft robotic fish which is optimized for speed and efficiency based on
experimental, numerical and theoretical investigation into oscillating
propulsion [16]. Wolf et al. use a pneumatically-actuated fish-like stationary
model to investigate how parameters like stiffness, strength, and frequency
affect thrust force generation [17]. Wolf et al. measure thrust, side forces,
and torques generated during propulsion and use a statistical linear model to
examine the effects of parameter combinations on thrust generation; they show
that both stiffness and frequency substantially affect swimming kinematics. We
are not aware of any work that uses experimental observations to
systematically identify the material parameters of a soft robotic design
through the use of a differentiable FEM simulation environment.

###  II-B Differentiable soft-body simulators

Our work is also relevant to the recent developments of robotic simulators,
particularly for soft robots. [18] present a differentiable multi-body
dynamics solver that is able to handle frictional contact for rigid and
deformable objects. [3] notably present a non-differentiable framework for
modeling, simulation, and control of soft-bodied robots using continuum
mechanics for modeling the robotic components and using Lagrange multipliers
for boundary conditions like actuators and contacts. Most related to our work
are the recent works on differentiable soft-body and fluid simulators [19, 4,
2, 7, 20, 21, 5]. These papers develop numerical methods for computing
gradients in a traditional soft-body or fluid simulator. Furthermore, they
demonstrate the power of gradient information in robotics applications, e.g.,
system identification or trajectory optimization. Most of the works present
simulation results only, with ChainQueen [20] and Real2Sim [2] the two notable
exceptions that discuss real-world soft-robot applications. Notably, [2]
optimizes visco-elastic material parameters of a finite element simulation to
approximate the dynamic deformations of real-world soft objects, such as an
open-loop controlled tendon-driven crawling robot.

###  II-C Hydrodynamic Surrogates

For underwater soft robots, the challenge of simulation is exacerbated by the
hydrodynamic interaction with the soft body. Several previous works tackle the
fluid-structure interaction problem through different methods, including
heuristic hydrodynamics [5, 22], physically-informed neural network approaches
[23], and data-driven learning approaches [24]. Note that many of these
hydrodynamic surrogates involve only the fluid and neglect the interaction
with the soft body. Compared to these previous methods, our work is the first
to apply differentiable simulation tools to a systematic approach for system
identification of real-world robots actuated by a pneumatically actuated
artificial muscles. In this paper, we focus on identification of a underwater
robot, however, our framework is applicable to other pneumatically actuated
soft mechanism.

##  III System Overview

Our pipeline (Figure 2) automatically learns the material parameters and
hydrodynamic thrust of a soft robotic fish using differentiable simulation and
measurements from quasistatic experiments in air and dynamic experiments
underwater. We examine a soft robotic fish in the form of a pneumatically
actuated soft bending actuator (see Figure 1). Despite this actuator choice,
our proposed method for system identification and simulation is agnostic to
any chosen type of soft robot. Matching will be reasonably good for any soft
fluidic actuator provided that stretches are not exceedingly large and the
material behaviour is approximately linear, which we demonstrate to be the
case for our soft robotic fish prototypes even under significant bending.

Figure 2: Flow diagram of our differentiable simulation and learning pipeline.
To learn the material parameters, we compute the loss as the Euclidean
distance between the measured positions and the simulated position data and
then minimize this loss using gradient-search enabled by our differentiable
FEM simulation. To learn the thrust force, our neural predictor network takes
as input the positions and pressure at a previous time step and outputs the
thrust force for the next time step.

In Section IV-C, forward simulation is accomplished using implicit Euler time
stepping and FEM spatial discretization. We leverage projective dynamics for a
substantial speed up in computation [4]. Our experimental setup needed to
gather the data for system identification and dynamic verification of our
simulation is described in Section V. In Section VI-A, the material properties
of the body and spine are specified by two different corotational materials
stitched together at the boundary. The whole structure of body and spine is
represented by a tetrahedral mesh. The pneumatic actuation of the interior
chambers of the body is accomplished via a traction boundary condition
specified as an external force normal to the deformed surface of the interior
pressure chambers. Our differentiable simulator enables system identification
of the two unknown Young’s moduli using a gradient-based descent algorithm
such as the Adam optimizer used for our system identification method [25]. We
verify the results (Figure 5) of our system identification by comparing
simulations to measurements from a gamut of dynamic experiments at varying
amplitudes and frequencies. The experiments are done for three different
designs of the fish tail including modifications of the material stiffness and
the number of air chambers within each pneumatic actuator. Finally, in Section
IV-B we describe the architecture and training of our thrust predictor
network.

##  IV Physical Modeling

###  IV-A Solid mechanics

The deformation of a soft material can be modeled by the equations of
elasticity. For the constitutive model we choose a corotated material [26]
since it is straightforward to implement for projective dynamics and can
predict bending to within millimeter precision for our application.

A corotated material is defined by its energy density function using the
Frobenius Norm and the trace operator,

$\Psi=\mu\left\lVert\mathbf{U}-\mathbf{I}\right\rVert^{2}_{F}+\frac{\lambda}{2}\text{tr}^{2}(\mathbf{U}-\mathbf{I}),$
(1)

where $\mathbf{U}$ is the stretch tensor, $\mathbf{I}$ is the identity, and
$\mu$ and $\lambda$ are the Lamé parameters, which have the following relation
to the Young’s modulus and Poisson ratio for an isotropic material,

$\mu=\frac{E}{2(1+\nu)},\quad\lambda=\frac{E\nu}{(1+\nu)(1-2\nu)}.$ (2)

To simulate a pneumatically actuated fish tail, we use a spatially-constant
time-varying pressure boundary condition on the interior surface of each
chamber of the tail.

###  IV-B Hydrodynamics

The hydrodynamic effects that produce thrust and drag on a swimmer are
important, however, the physics governing the fluid can be complicated due to
turbulent effects and two-way coupling between the solid and fluid domains.
One may choose to model water using the Navier-Stokes equations, however, in
practice simulating the full FSI problem is under time constraints
computationally unfeasible. For this reason, many workers in the field of
computer graphics and simulation choose to use simplified heuristic
hydrodynamic models such as the one presented by Min et al. [22]. Analytical
models for fish propulsion have also been studied extensively in the
literature [8, 27, 28], but these approximations require severe simplifying
assumptions necessary for analytical tractability. Using the elongated-body
theory (EBT) described in Section IX-A [8], the average thrust produced by a
fish is approximated by $f_{\textrm{EBT}}=m\hat{g}(\dot{x},\dot{y},y)$, where
$m$ is the virtual mass and $\hat{g}$ is a function of the velocity
$\dot{x},\dot{y}$ and deformation $y$ of the fish tail. We use EBT for
comparison to our learned hydrodynamic model Section IV-B.

###  IV-C Simulation

####  IV-C1 Discretization

We use our finite element method with tetrahedra to discretize the continuous
spatial domain occupied by the soft fish tail. Furthermore, we implement the
implicit Euler time-stepping scheme for time discretization. Formally, let
$\mathbf{x}_{i}$ and $\mathbf{v}_{i}$ be vectors of size $3N$ stacking up the
nodal positions and velocities from the tetrahedra at the $i$-th time step. We
then rewrite the governing equations in the following discretized form

$\displaystyle\mathbf{x}_{i+1}=$ $\displaystyle\mathbf{x}+h\mathbf{v}_{i+1},$
(3) $\displaystyle\mathbf{v}_{i+1}=$
$\displaystyle\mathbf{v}_{i}+h\mathbf{M}^{-1}[\mathbf{f}_{ela}(\mathbf{x}_{i+1})+\mathbf{f}_{act}].$
(4)

Here, $h$ is the time step, $M$ is the mass matrix, and $\mathbf{f}_{ela}$ and
$\mathbf{f}_{act}$ represent the elastic force and the actuator force
(computed by the pressure from the pneumatic actuator and the surface area of
the chambers), respectively. The time integration is implicit because we make
the elastic force depend on the nodal position $\mathbf{x}_{i+1}$ at the
beginning of the next time step instead of the current one.

####  IV-C2 Solver

Solving Equation 3 typically requires expensive numerical computation due to
its implicit time-stepping scheme. However, mature numerical tools exist since
this problem has been extensively studied in computer graphics and physics
simulation. In this work, we choose to use DiffPD [4], a recent differentiable
simulator with projective dynamics that is suitable for simulating terrestrial
and underwater soft robots. We choose differentiable projective dynamics over
other numerical solvers, for example traditional Newton’s methods, due to its
simultaneous accommodation of speed, robustness, and differentiability.

##  V Experimental Setup

In Figure 3, we show renderings of the three types of fishtail actuators we
investigate in this work. We vary material parameters as well as geometry.

Figure 3: Three different actuator designs. _Nemo_ and _Bruce_ share the same
geometry but _Bruce_ has a stiffer DS20 material for its body. _Nemo_ and
_Dory_ share the same materials but _Dory_ has a greater number of air
chambers.

We collect tail deformation, thrust force, and air chamber pressure data of
pneumatically actuated silicone fishtail both in air and in water. To this
end, the fishtail is mounted in a water tank and the fish head is rigidly
connected to a load cell (TAL2210, $1\text{\,}\mathrm{k}\mathrm{g}$ model) to
measure the thrust force in the heading direction (Figure 4). The load cell
data is obtained using an amplifier board and a microcontroller.111Amplifier:
HX711, Microcontroller: Arduino Leonardo Black markers are painted onto the
fish’s back and tracked using Kanade-Lucas-Tomasi (KLT) feature tracking in a
video captured from above using a GoPro Hero 6 camera [29, 30]. We found in
practice that three markers are sufficient for capturing the deformation of
the fishtail and can be used as input to our learning pipeline. An LED is used
to synchronize the pressure and the force data with the recorded video. The
fish is actuated with a Festo proportional valve manifold with controlled
pressure222MBA-FB-VI,
$0\text{\,}\mathrm{-}$$2\text{\,}\mathrm{b}\mathrm{a}\mathrm{r}$ range, 1%
accuracy, which outputs desired and actual pressure data for both air chamber
sides of the fishtail.

Figure 4: Thrust experiment setup actuated with a Festo proportional valve manifold. A 3D-printed fixture (black) connects the load cell (grey) to a fish adapter (opaque). The fish adapter is mounted within a casted soft robotic fish tail made of silicone elastomer. The deformation of the tail is captured using a GoPro positioned above the tail. Markers on the back of the fish are tracked using Kanade-Lucas-Tomasi (KLT) feature tracking. The thrust generated by the fish is measured using a load cell. TABLE I: Material, Young’s moduli, air chambers, and maximum applied pressure for each prototype. Reported material parameters are nominal values. Property | _Nemo_ | _Dory_ | _Bruce_
---|---|---|---
Body Material | DS10 | DS10 | DS20
Body Young’s Modulus [$\mathrm{MPa}$] | $0.1-0.25$ | $0.1-0.25$ | $1.1$
# of Chambers | 9 | 12 | 9
Spine Young’s Modulus [$\mathrm{GPa}$] | $2.5-5$ | $2.5-5$ | $2.5-5$
Max. Pressure [$\mathrm{bar}$] | 0.2 | 0.35 | 0.5

Two types of experiments are conducted using the same experimental setup:
Firstly, a quasistatic experiment is carried out, where the fishtail is
deformed by actuating one side of the fish with constant pressure. Secondly,
flapping experiments are conducted where the fishtail is actuated using a
square pressure wave, which leads to a natural motion of the tail. Both the
quasistatic and the flapping experiment are carried out with a variety of fish
designs using silicone rubber with different Young’s moduli.333Smooth-On
Dragon Skin 10/20 slow with shore hardness 10A/20A The fish materials,
geometries, and actuation pressures are shown in Table I. No pre-straining was
done on the fishtails prior to experiments and negligible hysteresis was found
when actuating under periodic pressure signals. The maximum observed
displacements in the fishtail were on the order of $0.2$, normalized to tail
length.

Since we are not fully constraining the fish to one-dimensional translation,
there exists bending moments that may result in spurious forces on the load
cell. However, we assume these to be negligible compared to the peak of the
thrust force. Furthermore, due to the rather small size of the tank, a
standing wave persists in the tank even quite long after interacting with it.
This can lead to a disturbance force on the load cell. However, since the
frequency of this disturbance is known, it can be filtered out during post-
processing. We also note that the measurement system has its own compliance
and damping, which account for the recoil and negative force readings.

##  VI System Identification

###  VI-A Learning Material Parameters

In this section, we explain our method for system identification. The process
begins with exporting two STL files from a CAD file of the soft actuator, one
for the body and one for the spine. To streamline the tetrahedral mesh
generation, we simplify the geometry by patching the small mounting holes used
during the casting process to stabilize the fish. Otherwise, the geometry of
the body and spine are imported as designed without modifications. In
practice, we have found that a correct chamber and spine geometry are
essential for adequate system identification of material parameters that are
physically plausible.

We tetrahedralize the surface triangle mesh by modeling it as a whole piece of
surface, a few holes representing the air chambers, and a set of internal
points defining the separation of spine and body. The tetrahedralization is
done according to the method by Hu et al. [31]. We control the target length
of edges to be $\nicefrac{{1}}{{50}}$ of the whole body length to get an
expressive and representative tetrahedral mesh.

After the tetrahedralization, we split the spine and body elements and assign
different Young’s moduli to them. The internal surfaces, which are identified
as air chambers, are modeled as actuators to drive the tail. In physical
experiments, we actuate the tail in the air and record its deformation through
the tracked markers. This information provides supervision to material
parameter search for both body and spine. We measure the error of matching by
the Euclidean distance between measurement data and simulation results from
the sensors.

We conduct our simulation experiments based on a linear elastic corotated
material. This material energy model is popular for physical simulation and
animation, but it is not the most accurate choice for closing the sim2real gap
for large deformations [26, 32]. The model can be improved by introducing a
more sophisticated elastic energy model, for example Neo-Hookean material,
which was not yet supported in our simulation framework. Furthermore, we chose
not to include dissipation explicitly in the material model since we learn
material parameters using static deformation data. For the dynamic
experiments, we observed that the deformations were sufficiently small for
accurate characterization and that the effect of material damping is small
thus having a negligible effect on the forced response.

###  VI-B Learning Hydrodynamics

To learn a simplified model of the hydrodynamic thrust force predictor
$f_{\textrm{thrust}}$, we use a feedforward neural network trained on data
from the experimental setup (see appendix). The inputs of the neural network
consist of the positions and velocities of the tracking points at time $t$,
together with the measured pressure and its time derivative, while the output
is a one-dimensional value for $f_{\textrm{thrust}}$ at time $t$. The
feedforward neural network consists of three hidden layers with RELU
activation functions and respectively a dimension of 200, 300, and 200 hidden
units. The network is trained with the Adam optimizer using a learning rate of
0.001 and a mean squared error loss. Note that our hydrodynamic modeling
method does not estimate fluid properties such as density or viscosity, but
rather the thrust force prediction, which is of more immediate use for
designers of swimming robots.

TABLE II: Specification of the training and test sets for the thrust predictor network. In the Nemo experiment, we generalize to unseen trials with seen parameters. In the Bruce experiment we generalize to unseen higher actuation frequency. Name | Mode | N. trials | Pressure | Frequency
---|---|---|---|---
Nemo | Training | 32 | 200,300 mbar | 1,2,3,4 Hz
Validation | 8 | 200,300 mbar | 1,2,3,4 Hz
Bruce | Training | 30 | 500,750 mbar | 1,2,3 Hz
Validation | 10 | 500,750 mbar | 4 Hz

In Table II, we provide a description of the training and the test set used
for the thrust predictor network. The table clarifies how we tested the
generalizability of our predictor. We split the experimental data to show the
ability of the network to generalize to unseen trials. For the fishtail named
_Nemo_ , while we keep the parameters the same, we generalize in the
validation set to unseen trials under various actuation amplitudes and
frequencies. In the experiments for _Bruce_ , we generalize to an unseen
higher actuation frequency.

##  VII Results

###  VII-A Learned Material Parameters

Figure 5: Top: Exhaustive search of the Young’s moduli pair for the spine and
body. The moduli pair with the lowest Euclidean loss is located at the red
$\times$. The red dot indicates the start of a gradient search and the white
line shows the progress towards convergence. Middle: Convergence comparison
between a gradient-based (Adam) and gradient-free search (CMA-ES). Bottom:
Comparison of static deformation of the _Nemo_ fish for increasing pressure in
experiment and simulation. Maximum displacement error increases with pressure,
but remains within the measurement error (on the order of the diameter of
markers). The grid in white under the real robot has
$10\text{\,}\mathrm{m}\mathrm{m}$ spacing. Figure 6: Learning four material
parameters from deformation data take with the Nemo fish prototype using a
gradient-based approach. The final values for Young’s moduli and Poisson’s
ratios depicted as an asterisk agree with plausible material parameter values
with lower loss.

For the _Nemo_ tail actuator described in Table I, we perform a grid search of
the loss landscape centered around the ground truth datasheet values for the
silicone and acetal materials of the body and spine. We report that if the
exact geometry of the actuator is reproduced with high fidelity in the
simulation, we converge to values within the range of typical measured values
for the material Young’s moduli (see Figure 5). Further, we see that there is
a unique minimum value that is within the acceptable range of measured moduli
for both parameters. Note that we assume a priori that the Poisson ratio for
silicone is $\nu\rightarrow 0.5$ or nearly incompressible, a standard
assumption for silicone, and the acetal sheet Poisson ratio is $\nu=0.37$ as
reported by the manufacturer. Typical values for the Young’s modulus of Dragon
Skin 10 range from $0.1\text{\,}\mathrm{M}\mathrm{P}\mathrm{a}$ to
$0.25\text{\,}\mathrm{M}\mathrm{P}\mathrm{a}$ and the value for Dragon Skin 20
was measured to be in the range of
$1.1\text{\,}\mathrm{M}\mathrm{P}\mathrm{a}$ [33] and typical values for the
Young’s modulus of acetal sheets range from
$2.5\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$ to
$5\text{\,}\mathrm{G}\mathrm{P}\mathrm{a}$.444https://dielectricmfg.com/knowledge-
base/acetal/ In Figure 6, we demonstrate that our method can be extended to
higher dimensional parameter spaces such as a search over both Young’s moduli
and Poisson’s ratios.

###  VII-B Gradient-based and Gradient-free Solver Methods

We compare the runtime of the gradient-free method _CMA-ES_ against the
runtime of the gradient-based method _Adam_ in Table III. The comparison shows
that although _CMA-ES_ [34] is slightly faster in runtime per iteration, the
use of the gradient-based method Adam [25] is significantly more effective for
convergence. These comparative experiments shown in Figure 5 were carried out
on a computer with an Intel Core i9-9900K @ 3.60GHz with 16 cores processor
and 64.0 GB of memory.

TABLE III: Comparison of Adam, CMA-ES, and Grid search. The forward simulation time is 318.9 seconds equivalent to grid search. Method |  Total Iterations |  Time Per Iteration |  Total Time |  Loss After 4 Iterations
---|---|---|---|---
Adam | 4 | 334.4 s | 22 m | 0.0027
CMA-ES | 40 | 321.0 s | 214 m | 0.021
Grid search | 25 | 318.9 s | 132.9 m | 0.023

###  VII-C Dynamic Experiments

For dynamic experiments, we compare the results of our simulation output with
the measured data of the furthest tracked dot on the tail. In Fig. 7, we
report our findings for sim2real performance in both the standard fish
(_Nemo_) and two other data sets for a tail with different Young’s moduli for
the body (_Bruce_) and a tail with a greater number of air chambers (_Dory_).
We demonstrate that if system identification is done correctly our simulation
results can predict the performance of a novel actuator design to within
millimeter precision. The dynamic experiments revealed that a larger actuation
pressure is required for similar amplitude displacements for the Bruce
prototype in comparison to Nemo since the material is stiffer. The prototype
Dory with a greater number of actuation chambers resulted in larger
deflections requiring greater pressure. For videos of the dynamic experiments
and simulation, we ask readers to refer to our supplemental material.

Figure 7: Simulation and measurement data for the _Nemo_ fish at (a) 200 mbar
and 2 Hz, (b) 200 mbar and 4 Hz, (c) the _Bruce_ fish at 500 mbar and 2 Hz,
and (d) the _Dory_ fish at 350 mbar and 2 Hz. The color bands indicate the
variance from $N=5$ trials. For both actuation signals, we are capable of
achieving sub-millimeter accuracy between experiment and simulation. The same
method is used to identify the parameters of Bruce and Dory, exhibiting
accuracy still within $3\text{\,}\mathrm{m}\mathrm{m}$.

###  VII-D Learned Hydrodynamics

We compare our simple predictor of thrust force with the measurement data and
the theoretical thrust from EBT in Figure 8. This serves as a benchmark of our
learning-based method against a classic, non-learning-based heuristic from the
literature. It is evident that the neural network model is capable of over
fitting to the training data. However, the model is also capable of
generalizing to actuation signals not used in the training for the same fish,
including at frequencies previously unobserved. However, we find that
generalization between fishes is difficult. In comparing to the thrust
prediction from EBT, we see that the analytical thrust is strictly positive.
The negative measured force is due to recoil of the measurement system. The
advantage of our approach is the differentiability and flexibility of the
neural network provided enough data is collected since estimating the virtual
mass of the fish necessary for EBT is generally imprecise. We note that the
network is capable of learning the measurement dynamics due to the compliance
of the load cell and the damping of the water though these effects can be
considered separately and corrected for if desired (see Section IX-B). We
conjecture that the asymmetry of the measured thrust may be due to
imperfections in the fabrication process of the fish tails favoring one
direction more strongly.

Figure 8: Tail lateral displacement (blue) with thrust measurement and
prediction. We compare the measured thrust force $\mathbf{f}_{\textrm{m}}$
(black), analytical thrust force from EBT $\mathbf{f}_{\textrm{EBT}}$
(magenta), and the neural network thrust prediction
$\mathbf{f}_{\textrm{thrust}}$ (green). We show the thrust prediction for two
training sets (a) and (c) and we report the time-averaged thrust for the two
test cases (b) and (d). We note that EBT tends to over-predict the thrust
measurement while our neural network thrust prediction accurately reproduces
the frequency for a given actuation signal and matches amplitude more robustly
than EBT.

##  VIII Conclusion

We demonstrated an experimentally-verified simulation framework that can be
used to accurately predict the deformations of a soft pneumatic actuator in
air and underwater. We focused on the case of a pneumatically actuated fish
tail with a flexible spine, but the framework is applicable to other soft
pneumatic actuators. Our pipeline can accurately learn material parameters and
hydrodynamics given quasi-static and dynamic data sets. The parameters we
found are not only within typical range of measured material parameters for
our materials, but can be used to successfully predict the behavior of dynamic
experiments for different pressure actuation amplitudes and frequencies to
within measurement error (millimeter accuracy). We further show how a data-
driven approach can be used to capture complex physics such as the
hydrodynamic effects from our fish tail.

Since our thrust predictor is also differentiable, we have the ability to
leverage in the future our fully differentiable pipeline for physically-
verified computational design. Our aim is to further progress towards a
systematic method by which soft roboticists can simulate and optimize their
designs and controllers, whether they be soft fish, manipulators, or other
flavors of soft robots. A fast and physically-verified co-optimization method
of design and control is the goal.

##  IX Appendix

###  IX-A Large Amplitude Elongated-Body Theory

The thrust of a fish swimming in an inviscid fluid can be described by
Lighthill’s Elongated-body theory (EBT) [8]. The thrust based on EBT depends
upon the full motion of the fish $x(a,t)$ and $y(a,t)$, where $a$ is the
parameter denoting the location along the fishtail with $a=0$ being at the
head and $a=L$ being at the tail.

To calculate the thrust predicted by EBT, we make use of a simplified version
of equation (6) in Lighthill [8],

$f_{\textrm{EBT}}\approx\frac{1}{2}m\dot{x}^{2},$ (5)

where we have evaluated the expression for the marker closest to the tail and
we have assumed the unit vectors $\frac{\partial x}{\partial a}\approx 0$ and
$\frac{\partial y}{\partial a}\approx 1$ when considering small tail
deformations in dynamic actuation.

We also use Lighthill’s expression for the virtual mass $m=\frac{1}{4}\pi\rho
s^{2}$, where $\rho$ is the density of the fluid and $s$ is the cross-
sectional depth, by estimating the geometry from the CAD model of our robotic
fish.

###  IX-B Load Cell Measurement Dynamics

We choose to model the dynamics of the load cell measurement jig (Figure 4) as
a lumped element system shown in Figure 9. We assume that the stiffness of the
fixture is far greater than the stiffness of the connection between the load
cell and the fish. Solving for the transfer function describing the load cell
dynamics, we have

$\frac{f_{\textrm{m}}}{f_{\textrm{hydro}}}=\frac{bs+k}{ms^{2}+bs+k},$ (6)

where $s$ is the complex Laplace variable. From this model, it is clear how
the measured force $f_{\textrm{m}}$ can be negative even for strictly positive
$f_{\textrm{hydro}}$ as discussed in Section VII-D.

It is possible to measure each parameter in Equation 6 and compensate it by
inverting the transfer function, however, our proposed data-driven neural
network force predictor is agnostic to the dynamics of the measurement system.
We are able to learn the system dynamics just from the data without the need
for characterizing the measurement system.

$k$$b$$m$$f_{\textrm{hydro}}$$f_{\textrm{m}}$ Figure 9: Load cell dynamics.
The force $f_{\textrm{m}}$ is measured with a pre-calibrated load cell. The
force $f_{\textrm{hydro}}$ is the actual hydrodynamic force experienced by the
fish. The stiffness $k$ and the damping coefficient $b$ describe the lumped
mechanical impedance of the measurement system. The lumped mass $m$ is
dominated by the fish. The wall is the connection point to the load cell.

## References

  * [1] R. Tedrake and the Drake Development Team, “Drake: Model-based design and verification for robotics,” 2019. [Online]. Available: https://drake.mit.edu
  * [2] D. Hahn, P. Banzet, J. Bern, and S. Coros, “Real2sim: Visco-elastic parameter estimation from dynamic motion,” _ACM Trans. Graph._ , vol. 38, no. 6, Nov. 2019.
  * [3] E. Coevoet, T. Morales-Bieze, F. Largilliere, Z. Zhang, M. Thieffry, M. Sanz-Lopez, B. Carrez, D. Marchal, O. Goury, J. Dequidt, _et al._ , “Software toolkit for modeling, simulation, and control of soft robots,” _Advanced Robotics_ , vol. 31, no. 22, pp. 1208–1224, 2017.
  * [4] T. Du, K. Wu, P. Ma, S. Wah, A. Spielberg, D. Rus, and W. Matusik, “Diffpd: Differentiable projective dynamics with contact,” _arXiv preprint arXiv:2101.05917_ , 2021.
  * [5] P. Ma, T. Du, J. Z. Zhang, K. Wu, A. Spielberg, R. K. Katzschmann, and W. Matusik, “Diffaqua: A differentiable computational design pipeline for soft underwater swimmers with shape interpolation,” _arXiv preprint arXiv:2104.00837_ , 2021.
  * [6] S. Bouaziz, S. Martin, T. Liu, L. Kavan, and M. Pauly, “Projective dynamics: Fusing constraint projections for fast simulation,” _ACM transactions on graphics (TOG)_ , vol. 33, no. 4, pp. 1–11, 2014.
  * [7] Y. Hu, L. Anderson, T.-M. Li, Q. Sun, N. Carr, J. Ragan-Kelley, and F. Durand, “DiffTaichi: Differentiable programming for physical simulation,” _International Conference on Learning Representations_ , 2020.
  * [8] M. J. Lighthill, “Large-amplitude elongated-body theory of fish locomotion,” _Proceedings of the Royal Society of London. Series B. Biological Sciences_ , vol. 179, no. 1055, pp. 125–138, 1971.
  * [9] A. D. Marchese, R. K. Katzschmann, and D. Rus, “A recipe for soft fluidic elastomer robots,” _Soft robotics_ , vol. 2, no. 1, pp. 7–25, 2015.
  * [10] R. K. Katzschmann, A. D. Marchese, and D. Rus, “Hydraulic autonomous soft robotic fish for 3d swimming,” in _Experimental Robotics_. Springer, 2016, pp. 405–420.
  * [11] R. K. Katzschmann, J. DelPreto, R. MacCurdy, and D. Rus, “Exploration of underwater life with an acoustically controlled soft robotic fish,” _Science Robotics_ , vol. 3, no. 16, 2018.
  * [12] J. Zhu, C. White, D. K. Wainwright, V. Di Santo, G. V. Lauder, and H. Bart-Smith, “Tuna robotics: A high-frequency experimental platform exploring the performance space of swimming fishes,” _Science Robotics_ , vol. 4, no. 34, 2019.
  * [13] G. Zheng, O. Goury, M. Thieffry, A. Kruszewski, and C. Duriez, “Controllability pre-verification of silicone soft robots based on finite-element method,” in _2019 International Conference on Robotics and Automation (ICRA)_. IEEE, 2019, pp. 7395–7400.
  * [14] R. K. Katzschmann, M. Thieffry, O. Goury, A. Kruszewski, T.-M. Guerra, C. Duriez, and D. Rus, “Dynamically closed-loop controlled soft robotic arm using a reduced order finite element model with state observer,” in _2019 2nd IEEE International Conference on Soft Robotics (RoboSoft)_. IEEE, 2019, pp. 717–724.
  * [15] M. Thieffry, A. Kruszewski, C. Duriez, and T.-M. Guerra, “Control design for soft robots based on reduced-order model,” _IEEE Robotics and Automation Letters_ , vol. 4, no. 1, pp. 25–32, 2018.
  * [16] S. C. van den Berg, R. B. Scharff, Z. Rusák, and J. Wu, “Biomimetic design of a soft robotic fish for high speed locomotion,” in _Conference on Biomimetic and Biohybrid Systems_. Springer, 2020, pp. 366–377.
  * [17] Z. Wolf, A. Jusufi, D. Vogt, and G. Lauder, “Fish-like aquatic propulsion studied using a pneumatically-actuated soft-robotic model,” _Bioinspiration & biomimetics_, vol. 15, no. 4, p. 046008, 2020.
  * [18] M. Geilinger, D. Hahn, J. Zehnder, M. Bächer, B. Thomaszewski, and S. Coros, “Add: analytically differentiable dynamics for multi-body systems with frictional contact,” _ACM Transactions on Graphics (TOG)_ , vol. 39, no. 6, pp. 1–15, 2020.
  * [19] T. Du, K. Wu, A. Spielberg, W. Matusik, B. Zhu, and E. Sifakis, “Functional optimization of fluidic devices with differentiable stokes flow,” _ACM Trans. Graph._ , vol. 39, no. 6, Nov. 2020. [Online]. Available: https://doi.org/10.1145/3414685.3417795
  * [20] Y. Hu, J. Liu, A. Spielberg, J. B. Tenenbaum, W. T. Freeman, J. Wu, D. Rus, and W. Matusik, “ChainQueen: A real-time differentiable physical simulator for soft robotics,” _International Conference on Robotics and Automation_ , 2019\.
  * [21] Z. Huang, Y. Hu, T. Du, S. Zhou, H. Su, J. Tenenbaum, and C. Gan, “Plasticinelab: A soft-body manipulation benchmark with differentiable physics,” in _ICLR 2021_ , 2021.
  * [22] S. Min, J. Won, S. Lee, J. Park, and J. Lee, “Softcon: simulation and control of soft-bodied animals with biomimetic actuators,” _ACM Transactions on Graphics (TOG)_ , vol. 38, no. 6, pp. 1–12, 2019.
  * [23] N. Wandel, M. Weinmann, and R. Klein, “Learning Incompressible Fluid Dynamics from Scratch – Towards Fast, Differentiable Fluid Models that Generalize,” _arXiv:2006.08762 [cs, stat]_ , Mar. 2021, arXiv: 2006.08762. [Online]. Available: http://arxiv.org/abs/2006.08762
  * [24] R. T. Chen, Y. Rubanova, J. Bettencourt, and D. Duvenaud, “Neural ordinary differential equations,” in _Proceedings of the 32nd International Conference on Neural Information Processing Systems_ , 2018, pp. 6572–6583.
  * [25] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” in _International Conference on Learning Representations_ , 2015.
  * [26] E. Sifakis and J. Barbic, “Fem simulation of 3d deformable solids: a practitioner’s guide to theory, discretization and model reduction,” in _Acm siggraph 2012 courses_ , 2012, pp. 1–50.
  * [27] T. Y.-T. Wu, “Swimming of a waving plate,” _Journal of Fluid Mechanics_ , vol. 10, no. 3, pp. 321–344, 1961.
  * [28] M. S. Triantafyllou, G. Triantafyllou, and D. Yue, “Hydrodynamics of fishlike swimming,” _Annual review of fluid mechanics_ , vol. 32, no. 1, pp. 33–53, 2000.
  * [29] B. D. Lucas, T. Kanade, _et al._ , “An iterative image registration technique with an application to stereo vision.” Vancouver, British Columbia, 1981.
  * [30] C. Tomasi and T. K. Detection, “Tracking of point features,” _Int. J. Comput. Vis_ , pp. 137–154, 1991.
  * [31] Y. Hu, Q. Zhou, X. Gao, A. Jacobson, D. Zorin, and D. Panozzo, “Tetrahedral meshing in the wild,” _ACM Transactions on Graphics (TOG)_ , vol. 37, no. 4, pp. 1–14, 2018.
  * [32] S. Adeeb, _Introduction to solid mechanics and finite element analysis using Mathematica_. Kendall Hunt, 2011\.
  * [33] L. Marechal, P. Balland, L. Lindenroth, F. Petrou, C. Kontovounisios, and F. Bello, “Toward a common framework and database of materials for soft robotics,” _Soft Robotics_ , 2020.
  * [34] C. Igel, N. Hansen, and S. Roth, “Covariance matrix adaptation for multi-objective optimization,” _Evolutionary computation_ , vol. 15, no. 1, pp. 1–28, 2007.