# Distribution and Kinematics of H I through Raman He II Spectroscopy of NGC
6302

Seok-Jun Chang Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-
Stra$\beta$e 1, 85748 Garching b. München, Germany  Hee-Won Lee Sejong
University, 209 Neungdong-ro, Gwangjin-gu, Seoul 05006, Korea  Jiyu Kim
Sejong University, 209 Neungdong-ro, Gwangjin-gu, Seoul 05006, Korea  Yeon-Ho
Choi  Korea Astronomy and Space Science Institute, 776 Daedeokdae-ro, Yuseong-
gu, Daejeon 34055, Korea  University of Science and Technology, 217, Gajeong-
ro, Yuseong-gu, Daejeon 34113, Korea

###### Abstract

The young planetary nebula NGC 6302 is known to exhibit Raman-scattered He II
features at 6545 Å and 4851 Å. These features are formed through inelastic
scattering of He II$\lambda\lambda$ 1025 and 972 with hydrogen atoms in the
ground state, for which the cross sections are $1.2\times 10^{-21}$ and
$1.4\times 10^{-22}{\rm\ cm^{2}}$, respectively. We investigate the spectrum
of NGC 6302 archived in the ESO Science Portal. Our Gaussian line fitting
analysis shows that the Raman-scattered He II features are broader and more
redshifted than the hypothetical model Raman features that would be formed in
a cold static H I medium. We adopt a simple scattering geometry consisting of
a compact He II emission region surrounded by a H I medium to perform Monte
Carlo simulations using the radiative transfer code STaRS. Our simulations
show that the H I region is characterized by the H I column density $N_{\rm
HI}=3\times 10^{21}{\rm\ cm^{-2}}$ with the random speed component $v_{\rm
ran}=10{\rm\ km\ s^{-1}}$ expanding with a speed $v_{\rm exp}=13{\rm\ km\
s^{-1}}$ from the He II emission region. Based on our best fit parameters, we
estimate the H I mass of the neutral medium $M_{\rm HI}\simeq 1.0\times
10^{-2}\ {\rm M_{\odot}}$, pointing out the usefulness of Raman He II
spectroscopy as a tool to trace H I components.

Radiative transfer — Planetary nebulae — Scattering — Individual NGC 6302

††facilities: ESO Science Archive Facility

##  1 Introduction

NGC 6302 is a young planetary nebula exhibiting a well-known butterfly
morphology. The two main lobes are divided by an equatorial torus composed of
atomic, molecular, and dusty material (Matsuura et al., 2005; Kastner et al.,
2022). With high helium and nitrogen abundances, NGC 6302 is classified as a
Type I planetary nebula according to the classification scheme proposed by
Peimbert (1978). It belongs to the highest excitation class with prominent
emission lines including N V$\lambda\lambda$1238, 1243, C
IV$\lambda\lambda$1548, 1551 and Ne VI at $7.7{\rm\ \mu m}$ (Feibelman, 2001;
Pottasch et al., 1985). Many researchers investigated NGC 6302 for the
internal kinematics and various components, including the ionized, atomic,
molecular and dust components.

Meaburn et al. (2008) investigated the proper motions of the outflowing knots
to propose that the distance to NGC 6302 is 1.17 kpc. They also derived a
kinematic age of 2200 years from their analysis of the Hubble-type expansion
(e.g., Szyszka et al., 2011). Dinh-V-Trung et al. (2008) proposed that the
molecular torus is expanding with a speed of $\sim 15{\rm\ km\ s^{-1}}$ from
their measurement using the Submillimeter Array. Santander-García et al.
(2017) conducted a kinematical analysis using ALMA data to confirm the Hubble
type expansion.

With recent history of significant mass loss, a young planetary nebula is
expected to harbor abundant H I behind the ionization front, which is
expanding with respect to the central hot source. The hyperfine structure 21
cm line is regarded as currently the most effective spectroscopic tracer of
atomic hydrogen. The first successful detection of a neutral component in
planetary nebulae was made by Rodriguez & Moran (1982), who conducted radio
observations of NGC 6302 using the Very Large Array. CO and H I components
were also detected from radio observations for a number of planetary nebulae
(e.g., Gussie & Taylor, 1995). However, severe confusion from the Galactic
emission prevents one from investigating the distribution and the kinematics
of the H I components in planetary nebulae.

For young planetary nebulae, a very unique and useful spectroscopic probe is
provided by the Raman scattering process of far UV line radiation with atomic
hydrogen. In young planetary nebulae, far UV He II lines can be Raman
scattered by atomic hydrogen to form broad features blueward of hydrogen
Balmer lines (Nussbaumer et al., 1989). The first report of Raman-scattered He
II features was made for the young planetary nebula NGC 7027 by Péquignot et
al. (1997), who identified the broad feature at 4852 Å as Raman-scattered He
II. Groves et al. (2002) found the same feature in NGC 6302 while they
investigated the extinction in the nebula. Subsequently, Raman-scattered He II
at 6545 Å was detected in the young planetary nebulae IC 5117, NGC 6790, NGC
6881, and NGC 6886 (Lee et al., 2006; Kang et al., 2009; Choi & Lee, 2020).

Raman-scattered He II features are clearly detected in the high resolution
optical spectrum of NGC 6302 provided by the ESO Science Archive Facility. In
this paper, we investigate the physical properties of H I in NGC 6302 using
these features. The paper is organized as follows. In Section 2, we briefly
explain the basic atomic physics of Raman scattering with atomic hydrogen. In
Section 3, we analyze Raman-scattered He II features at 4851 Å as well as at
6545 Å. We present the results from our Monte Carlo simulations in Section 4.
A brief summary and discussions are presented in the final section.

Figure 1: A schematic illustration to show the formation of Raman-scattered He
II features blueward of H I Balmer lines in the H I region neighboring the He
II emission region. Far UV He II lines at 1025, 972 and 949 Å are slightly
more energetic than H I Ly$\beta$, Ly$\gamma$, and Ly$\delta$, resulting in
optical lines at 6545 Å, 4851 Å, and 4332 Å blueward of H$\alpha$, H$\beta$
and H$\gamma$.

##  2 Atomic Physics

Schmid (1989) proposed that the broad emission features at 6825 Å and 7082 Å
that appear in about a half of symbiotic stars are formed through inelastic
scattering processes of the far UV resonance doublet lines O
VI$\lambda\lambda$1032 and 1038. When a far UV photon blueward of Ly$\alpha$
incident on a hydrogen atom in the ground state can be converted into a lower
energy photon by an amount of 10.2 eV, the energy of a Ly$\alpha$ photon, then
the hydrogen atom makes a final de-excitation into the excited $2s$ state
instead of the initial ground state.

Additional examples are provided from far UV He II lines. He II $\lambda$1025
arising from transitions of $6\to 2$ is Raman scattered with atomic hydrogen
to form an optical spectral feature at 6545 Å, blueward of He II
$\lambda$6560, which is associated with transitions $6\to 4$. Similarly, Raman
scattering of He II $\lambda$972 and $\lambda 949$ yields spectral features at
4851 Å and 4332 Å, respectively (Nussbaumer et al., 1989; Lee, 2012).

Figure 1 shows a schematic illustration of the Raman scattering process
expected to operate in young planetary nebulae, when far UV He II line photons
enter the neutral region behind the ionization front. The atomic line center
wavelengths for far UV He II emission lines and their optical Raman lines are
shown in Table 1, where the cross sections $\sigma_{1s}^{\rm Ray}$ and
$\sigma_{2s}^{\rm Ram}$ for Rayleigh and Raman scattering are also available
(e.g., Lee, 2012; Chang et al., 2015).

The energy conservation requires

$\nu_{i}=\nu_{o}+\nu_{\rm Ly\alpha},$ (1)

where $\nu_{i},\ \nu_{o}$ and $\nu_{\rm Ly\alpha}$ are frequencies of the
incident, Raman-scattered and Ly$\alpha$ photons, respectively. One may
immediately note that the line widths of the incident and the Raman-scattered
radiation are related by

${\Delta\nu_{i}\over\nu_{i}}=\left({\nu_{o}\over\nu_{i}}\right){\Delta\nu_{o}\over\nu_{o}},$
(2)

which shows that the Raman He II features blueward of Balmer lines are
broadened by the factor $\nu_{i}/\nu_{o}$. The line broadening effect
associated with the inelasticity of Raman scattering allows one to readily
identify Raman-scattered features (Schmid, 1989).

Table 1: Atomic line center wavelength $\lambda_{0}$ of He II emission and Raman He II (first and second columns) and cross sections for Rayleigh scattering and Raman scattering into $2s$ (third and fourth columns) Transition |  $\lambda_{0}$ of He II Emissiona [Å] |  $\lambda_{0}$ of Raman He II b [Å] |  $\sigma^{\rm Ram}_{2s}$ [$\rm cm^{2}$] |  $\sigma^{\rm Ray}_{1s}$ [$\rm cm^{2}$]
---|---|---|---|---
$n=6\to 2$ | 1025.28 | 6544.70 | $1.2\times 10^{-21}$ | $6.2\times 10^{-21}$
$n=8\to 2$ | 972.13 | 4851.30 | $1.4\times 10^{-22}$ | $8.3\times 10^{-22}$
$n=10\to 2$ | 949.32 | 4331.74 | $2.9\times 10^{-23}$ | $1.9\times 10^{-22}$
a vacuum wavelength, b air wavelength |  |  |  |

##  3 Observation

###  3.1 Data

Figure 2: The optical spectrum of NGC 6302 retrieved from the UVES ESO public
data archive. Strong emission lines, including H$\alpha$, H$\beta$, [O
III]$\lambda\lambda$4959, 5007, and [N II]$\lambda\lambda$6548, 6583 are
heavily saturated and He II emission lines at 4859 Å and 6560 Å are barely
noticeable. Due to the heavy saturation, [N II]$\lambda$6548 appears
comparable to [N II]$\lambda$6583, even though [N II]$\lambda$6548 is in fact
3 times weaker than [N II]$\lambda$6583\. The insets of the upper and lower
panels zoom in spectral regions around H$\beta$ and H$\alpha$, respectively.
In the insets, the two broad features at 4851 Å and 6545 Å are clearly seen,
which are Raman-scattered He II features.

We retrieved the spectra of NGC 6302 taken with UV-Visual Echelle Spectrograph
(UVES) attached on ESO Very Large Telescope (VLT) from ESO Science
Archive111http://archive.eso.org/scienceportal/home under the ESO programme
65.I-0465 (P.I. Casassus). NGC 6302 was observed on May 24, 2000 through a 0.3
arcsec slit, resulting in spectral resolution of $R\sim 107,200$. A total of 4
spectra were obtained. The exposure times are 1200 seconds for the first 2
spectra and 2400 seconds for the other 2 spectra. The median signal-to-noise
ratio of each exposure ranges from 7.1 to 11.5.

The 1-D and wavelength-calibrated spectrum covers the spectral range from 4727
Å to 6835 Å. In Figure 2, we show the full UVES spectrum of NGC 6302. The
strong lines of H$\alpha$, H$\beta$, [O III]$\lambda\lambda$4959, 5007 and [N
II]$\lambda\lambda$6548, 6583 are all severely saturated in this spectrum. The
presence of He II emission lines at 6560 Å and 4859 Å is barely noticed in the
scale.

In the inset of the upper panel of Figure 2, we show a part around H$\beta$ to
show clearly the broad feature at 4851 Å, which is identified as Raman-
scattered He II. In a similar way, the inset of the lower panel shows a part
of the spectrum near H$\alpha$, where we find a broad feature blended with [N
II]$\lambda$6548\. The intrinsic flux ratios $F_{6583}/F_{6548}$ of [N
II]$\lambda\lambda 6548,6583$ as well as $F_{5007}/F_{4959}$ of [O
III]$\lambda\lambda 4959,5007$ is fixed to be 3 due to the fact that they
arise from the same excited state. Hence, if the broad weak emission feature
near [N II]$\lambda$6548 is due to N II, then it implies the presence of a
more conspicuous broad feature by a factor 3 near [N II]$\lambda$6583\. No
such broad feature is apparent near [N II]$\lambda$6583, which confirms that
the broad feature around [N II]$\lambda$6548 is not associated with [N II] but
is to be identified with Raman-scattered He II.

###  3.2 Gaussian Fitting

In this subsection, we provide a quantitative profile analysis of He II
emission lines and Raman-scattered He II lines. In addition to He II
$\lambda$6560 and He II $\lambda$4859, we also consider He II $\lambda$6527,
which is formed from transitions between $n=14$ and $n=5$. This emission line
is very convenient because Raman He II $\lambda$6545 is located redward of He
II $\lambda$6527 with comparable strength. The atomic line centers of these He
II $\lambda$4859, $\lambda$6527 and $\lambda$6560 are 4859.32 Å, 6527.10 Å,
6560.10 Å, respectively (e.g., Lee et al., 2006; Hyung & Feibelman, 2004).

The He II emission lines and Raman-scattered features of NGC 6302 appear to be
well-fitted using a single Gaussian function in the form given by

$G(\lambda)={F\over\sqrt{2\pi}\
\sigma_{\lambda}}\exp\left(-{(\lambda-\lambda_{c})^{2}\over{2\sigma_{\lambda}^{2}}}\right),$
(3)

where $\lambda_{c}$ and $\sigma_{\lambda}$ are the wavelength of line center
and the line width, respectively.

In Figure 3, we show the result of our single Gaussian fits. The observational
data are shown with the solid gray lines, and the Gaussian fitting results are
shown with the black dashed lines. In the two left panels, the results for He
II $\lambda$6527 and $\lambda$6560 and Raman He II $\lambda$6545 are shown,
and the right panels show the results for the counterparts blueward of
H$\beta$. The left lower panel shows the Raman He II $\lambda$6545, which is
severely blended with [N II]$\lambda$6548\. The right lower panel shows the
Raman He II $\lambda$4851.

Table 2 provides the fitting parameters used for the results that are shown
with the black dashed lines in Figure 3. The widths $\sigma_{\lambda}$ of the
two He II emission lines $\lambda$6560 and $\lambda$6527 near H$\alpha$ are
0.283 Å, and that of He II $\lambda$4859 is 0.212 Å. In contrast, the line
width $\sigma_{\lambda}=3.25{\rm\ \AA\ }$ of Raman He II $\lambda$6545 is 11.5
times wider than those of He II $\lambda$6527 and $\lambda$6560\. In the case
of Raman He II $\lambda$4851, the line width $\sigma_{\lambda}=1.70$ Å is
wider than that of He II $\lambda$4859 by a factor of 8.02. These factors
exceed the line broadening factor $\lambda_{o}/\lambda_{i}$ due to the
inelasticity of Raman scattering given in Eq. (2). This is consistent with the
proposal that the Raman-scattered He II features of NGC 6302 are formed in
expanding neutral regions, where the expansion provides additional line
broadening (e.g. Jung & Lee, 2004; Choi et al., 2020).

Figure 3:  Gaussian line profile fits of NGC 6302 near H$\alpha$ (left panels)
and H$\beta$ (right panels). The observed data are represented by the gray
solid lines. The black dashed lines show the results of our single Gaussian
fit adopting parameters in Table 2. In the lower two panels, the red dashed
lines show hypothetical Raman-scattered features that would be formed in a H I
region stationary with respect to the He II emission region. The observed
Raman He II are broader and more red-sfhited than the hypothetical profiles.
Table 2: Parameters of Gaussian fitting in Figure 3

Line | He II 6560 | Raman He II 6545 | He II 6527 | He II 4859 | Raman He II 4851
---|---|---|---|---|---
$\lambda_{c}$ [Å] | 6559.00 | 6546.06 | 6525.98 | 4858.45 | 4851.45
$\sigma_{\lambda}$ [Å] | 0.283 | 3.25 | 0.283 | 0.212 | 1.70
F [$10^{-14}\rm erg\,s^{-1}\,cm^{-2}$] | 7.08 | 1.17 | 0.292 | 1.16 | 0.121
$\sigma_{v}$ [ km s-1 ] | 12.9 | 23.4a | 12.9 | 13.09 | 21.03a
$\Delta V_{c}$ [ km s-1 ] | - | 17.86b1 | - | - | 12.88b2

a: corrected for Raman broadening, b1 & b2: velocity offset from He II
$\lambda$6560 & $\lambda$4859

###  3.3 Line Ratio of He II $\lambda$6560 & $\lambda$4859

The observed flux ratio of H$\alpha$/H$\beta$ is regarded as an extinction
indicator, because the Case B recombination theory predicts the flux ratio of
$\sim 2.8$ (e.g. Osterbrock, 1989; Storey & Hummer, 1995). However, H$\alpha$
and H$\beta$ in the UVES spectra are saturated near the line center,
preventing one from estimating the extinction. Instead, the nearby emission
lines He II $\lambda$6560 and $\lambda$4859 can be used for the same purpose.

According to Case B recombination theory for He II, the intrinsic line ratio
of He II $\lambda$6560 and $\lambda$4859 is $\sim 2.5$, whereas it is observed
to be $\sim 6.12$ from the UVES spectra (Storey & Hummer, 1995). By adopting
the dust model of our Galactic interstellar medium provided by Draine (2003),
the observed line ratio of He II $\lambda$6560 and $\lambda$4859 is consistent
with $N_{\rm H}=1.35\times 10^{22}{\rm\,cm^{-2}}$ corresponding to $A_{V}\sim
7.5$.

###  3.4 Line Broadening

Raman He II lines are broadened significantly due to inelasticity of Raman
scattering, which is illustrated in Eq. (2). Therefore, a correction factor of
$\lambda_{\rm UV}/\lambda_{\rm Ram}$ is required in order to yield the
velocity width in the parent UV spectral space, where $\lambda_{\rm UV}$ and
$\lambda_{\rm Ram}$ are the wavelengths of the incident UV He II emission and
the corresponding Raman scattered feature, respectively. The velocity width
$\sigma_{v}$ of Raman-scattered He II is calculated by

$\sigma_{v}={{c\sigma_{\lambda}}\over{\lambda_{c}}}{\left(\lambda_{\rm
UV}\over\lambda_{\rm Ram}\right)}=\sigma_{v}^{\rm app}{\left(\lambda_{\rm
UV}\over\lambda_{\rm Ram}\right)},$ (4)

where $\sigma_{v}^{\rm app}=c\sigma_{\lambda}/\lambda_{c}$ is an apparent
velocity width of Raman He II. For example, the apparent width
$\sigma_{\lambda}=3.25{\rm\ \AA}$ for Raman He II $\lambda$6545 is converted
to $\sigma_{v}=23.4{\rm\ km\ s^{-1}}$ instead of $\sigma_{v}^{\rm
app}=150{\rm\ km\ s^{-1}}$.

Our single Gaussian fit analysis shows that $\sigma_{v}=23.4{\rm\ km\ s^{-1}}$
and $21.0{\rm\ km\ s^{-1}}$ for Raman He II $\lambda$6545 and $\lambda$4851,
respectively. These values of $\sigma_{v}$ for Raman He II lines are broader
than those of He II emission lines by $8-10\ {\rm\,km\,s^{-1}}$. In the lower
two panels of Figure 3, we show the hypothetical profiles of Raman scattered
features, which would be formed in a static H I medium with
$\sigma_{v}=13{\rm\,km\,s^{-1}}$. In Section 4, we use a Monte Carlo approach
to find that the difference of $\sigma_{v}$ between Raman He II and He II
emission lines originates from the random motion of scattering medium.

###  3.5 Line Center Shift

In the bottom row of Table 2, we show the velocity offset $\Delta V_{c}$ of
Raman He II lines relative to nearby optical He II emission lines. The
apparent velocity line shift $V_{\rm app}$ is given by

$V_{\rm app}=\left({\lambda_{c}-\lambda_{0}}\over\lambda_{0}\right)c$ (5)

where $\lambda_{0}$ is the atomic line center wavelength. In turn, the
velocity offset $\Delta V_{c}$ of Raman He II $\lambda$6545 is calculated by

$\Delta V_{c}={\lambda_{0,1025}\over\lambda_{0,6545}}V_{\rm app,6545}-V_{\rm
app,6560}.$ (6)

The wavelengths $\lambda_{0,1025},\lambda_{0,6545}$, and $\lambda_{c,6545}$
are defined analogously. In the case of Raman He II $\lambda$4851, the
velocity offset $\Delta V_{c}$ is obtained with $\lambda_{0,972}$ and
$\lambda_{0,4859}$ replacing $\lambda_{0,1025}$ and $\lambda_{0,6560}$,
respectively.

From our line fit analysis, $\Delta V_{c}=10.24$ km s-1 and 13.10 km s-1 for
Raman He II $\lambda$6545 and $\lambda$4851, respectively. In the bottom
panels of Figure 3, the Raman He II features are clearly displayed redward of
the hypothetical features shown in the red dashed lines. This redward line
shift of Raman He II indicates that the H I medium is moving away from the He
II emission region (e.g. Choi et al., 2020).

###  3.6 Raman Conversion Efficiency

In this subsection, we present the Raman conversion efficiency (RCE), which is
defined as the ratio of the number of Raman-scattered photons divided by that
of the incident far UV photons. Explicitly, RCE of Raman He II $\lambda$6545
is given by

$\displaystyle{\rm RCE}_{6545}$ $\displaystyle=$
$\displaystyle{{F_{6545}/E_{6545}}\over{F_{1025}/E_{1025}}}$ (7)
$\displaystyle=$
$\displaystyle\left(\lambda_{0,6545}\over\lambda_{0,1025}\right){\left({F_{6545}}\over{F_{6560}}\right)}\left(F_{6560}\over
F_{1025}\right),$

where $E_{\lambda}=hc/\lambda$ is the energy of a photon with wavelength
$\lambda$. RCE of Raman He II $\lambda$4851 is also given in an analogous way
by replacing 1025 and 6560 by 972 and 4859.

Firstly, the ratio of Raman and neighboring optical He II line fluxes,
$({F_{6545}}/{F_{6560}})$ or $({F_{4851}}/{F_{4859}})$, is obtained using the
observed values presented in Table 2. Secondly, the flux ratio of the
neighboring optical He II and directly unavailable incident far UV He II
lines, $(F_{6560}/F_{1025})$ or $(F_{4859}/F_{972})$, is deduced, assuming
that Case B recombination is valid. According to Storey & Hummer (1995),
$(F_{6560}/F_{1025})\sim 0.19$ and also $(F_{4859}/F_{972})\sim 0.19$.
Combining these results, we obtain RCE${}_{6545}=0.21$ and
RCE${}_{4851}=0.10$.

Figure 4: Schematic illustration of the scattering geometry considered in this
work for Monte Carlo radiative transfer. The scattering geometry is composed
of a point-like central He II emission source (orange), a disk-like H I region
(green), and an outer dusty and molecular region (gray). Closely related to
the line shift and broadening of Raman He II, the H I region has two kinematic
components: one is a radial expansion with a velocity ${v_{\rm exp}}$, and the
other is a random motion with a representative speed of $v_{\rm ran}$. The
scattering geometry is specified by the inner and outer radii $R_{i}$, $R_{o}$
and the half opening angle $\theta_{o}$ of the H I region with respect to the
central source.

##  4 Monte Carlo Approach

Figure 5: Raman conversion efficiency (RCE) of Raman He II 4851 (right) and
6545 (center) lines and the ratio of RCE of those two Raman He II lines in our
model of Figure 4 for various half opening angle $\theta_{\rm o}$. The black
vertical lines represent the values from Gaussian fitting in Table 2.

Raman He II lines carry important physical information of the H I region. The
line center, profile width and the strength of Raman-scattered He II can be
used to put strong constraints on the distribution and kinematics of the H I
component near the He II emission nebula (e.g., Nussbaumer et al., 1989). Due
to the difference in scattering cross sections of He II $\lambda$ 1025 and He
II $\lambda$972, the Raman conversion efficiencies also differ so that much
detailed information can be obtained if both Raman-scattered features are
secured with sufficient data quality. In particular, Choi et al. (2020)
investigated the line formation of Raman-scattered He II in an expanding H I
region to show that the redward shift of the line center of Raman-scattered He
II features is conspicuously enhanced due to the sharp rise of the cross
section toward H I Lyman line centers (see also Jung & Lee, 2004). In this
section, we adopt a Monte Carlo approach to propose a simple scattering
geometry consistent with the observed spectrum considered in this work.

The Monte Carlo code ’STaRS’ developed for radiative transfer in thick neutral
regions by Chang & Lee (2020) is used to simulate the formation of Raman-
scattered He II and find the best fitting parameters. In this simulation, we
set the numbers of the incident far UV He II 1025 and 972 photons to be
$10^{8}$. The ratio of the two UV lines is $\sim 2.5$ in accordance with Case
B recombination. Because optical Raman-scattered He II features are subject to
dust extinction after leaving the H I region, the Raman conversion efficiency
is computed using the optical Raman photons before dust extinction. However,
the final line fitting of Raman-scattered He II is carried out after
correction for dust extinction. The incident far UV He II photons are assumed
to be described by a single Gaussian profile with a velocity width
$\sigma_{v}=13{\rm\ km\ s^{-1}}$ corresponding to the widths of He II
$\lambda$6560 and $\lambda$4859 as illustrated in Table 2.

###  4.1 Scattering Geometry

The central region of NGC 6302 is highly obscured and can be probed with high
angular resolution observations achievable by radio interferometry. Peretto et
al. (2007) investigated the kinematics of the molecular torus to report the
expansion velocity of $8{\rm\ km\ s^{-1}}$. Wright et al. (2011) presented
their 3D photoionization computation to propose that the inner and outer radii
of the circumstellar disk are $r_{\rm in}=1.2\times 10^{16}{\rm\ cm}$ and
$r_{\rm out}=3.0\times 10^{17}{\rm\ cm}$, respectively, based on their best-
fitting model result. Furthermore, the circumstellar disk is geometrically
thin with a half-opening angle $\sim 10^{\circ}$.

In our Monte Carlo simulation, the He II emission region is assumed to be an
unresolved compact source surrounded by a circumstellar disk, where Raman
scattering takes place. Figure 4 shows a schematic illustration of the
scattering geometry considered in this work. For the sake of simplicity, the
neutral region is assumed to take the form of a disk with the inner and outer
radii $R_{\rm i}=10^{16}{\rm\ cm}$ and $R_{o}=5R_{\rm i}$ and also
characterized by the half opening angle $\theta_{o}$ of the H I disk with
respect to the point-like central He II region. The neutral region is of
uniform H I density $n_{\rm HI}={N_{\rm HI}}/(R_{o}-R_{i})$, where $N_{\rm
HI}$ is H I column density.

The H I medium is assumed to move away from the central He II emission region
in the radial direction. In addition, we denote by $v_{\rm ran}$ the random
motion component contributed by the thermal motion. A Hubble-type outflow is
chosen in this work in accordance with the observations of the ionized and
molecular components by Szyszka et al. (2011) and Santander-García et al.
(2017), respectively. Specifically, the radial velocity $v(r)$ at a distance
$r$ from the He II source is chosen to follow

$v(r)={v_{\rm exp}}\left({r\over R_{o}}\right),$ (8)

where the parameter $v_{\rm exp}$ is the expansion velocity at the outer
radius $R_{o}$.

In our simulation, we collect Raman scattered photons escaping along the line
of sight, which coincides with the direction specified by the polar angle
$\theta$= $90^{\circ}$, in view of the fact that the central star of NGC 6302
is highly obscured by dust (Kastner et al., 2022). Thus, we consider dust
extinction in the line of sight, which is coincident with the equatorial
direction of the scattering medium. In Figure 4, a dust component is added
outside the neutral region so that optical Raman He II photons are subject to
dust extinction before reaching the detector. The dust optical depth is chosen
to be consistent with the reddening found in the line ratio of He II
$\lambda$6560 and $\lambda$4859 discussed in Section 3.3. The central He II
emission region is assumed to inject far UV He II line photons with a profile
described by a single Gaussian and strengths in accordance with the Case B
recombination theory.

###  4.2 Simulated Raman Conversion Efficiency

Figure 6: Best fit simulation profiles of Raman He II 4851 (left) and 6545
(right) superimposed on the observed spectrum (gray) and the single Gaussian
fits (black). The simulated profiles are shown by red and blue solid lines in
the left and right panels, respectively. The model parameters are
$\theta_{0}=20^{\circ}$, $N_{\rm HI}=3\times 10^{21}{\rm\ cm^{-2}}$, ${v_{\rm
exp}}=13{\rm\,km\,s^{-1}}$, and ${v_{\rm ran}}=10{\rm\,km\,s^{-1}}$.

In Figure 5, the Raman conversion efficiencies for Raman-scattered He II at
6545 Å and 4851Å are shown in the middle and right panels for three values of
the half opening angle $\theta_{\rm o}$ = 10, 20, and 30∘ and for a range of H
I column densities $10^{21}-10^{22}{\rm\,cm^{-2}}$. Here, we fix the expansion
speed $v_{\rm exp}=13{\rm\,km\,s^{-1}}$ and the random speed ${v_{\rm
ran}}=10{\rm\,km\,s^{-1}}$ according to the best fitting result (see Appendix
A). The horizontal dotted lines indicate RCEs of 0.21 and 0.10 for Raman-
scattered He II at 6545 and 4851, respectively, presented in Section 3.6. In
the left panel, the ratio of the two Raman conversion efficiencies is shown.
The horizontal dotted line represents the ratio of the two RCEs $\sim 2.1$.

In the range of H I column density $N_{\rm HI}$ considered in Figure 5, the
Raman conversion efficiency is nearly proportional to the half opening angle
$\theta_{0}$. However, in contrast, the ratio of the two Raman conversion
efficiencies is relatively insensitive to $\theta_{0}$ and decreases as
$N_{\rm HI}$ increases. In this range of $N_{\rm HI}$, too small RCEs are
obtained from the cases for $\theta_{\rm o}$ = $10^{\circ}$ and those with
$\theta_{\rm o}$ = $30^{\circ}$ lead to much larger RCEs to account for the
measured values. In view of this, we propose that the measured RCEs are
consistent with the scattering geometry characterized by ${N_{\rm HI}}=3\times
10^{21}{\rm\,cm^{-2}}$ and $\theta_{\rm o}$ = $20^{\circ}$.

###  4.3 Best Fit Profiles

In the left panel of Figure 6, we show the best line fit to the observed
spectrum of Raman He II $\lambda$4851, for which the parameters adopted are
$v_{\rm exp}$ = 13 km s-1 and $v_{\rm ran}$ = 10 km s-1 in addition to $N_{\rm
HI}=3\times 10^{21}{\rm\ cm^{-2}}$ and $\theta_{\rm o}=20^{\circ}$. The red
line shows the best fit and the black dotted line is the Gaussian fit to the
observed data.

Using our best fit parameters, the H I number density is estimated to $n_{\rm
HI}=7.5\times 10^{4}{\rm\ cm^{-3}}$, from which we estimate the total mass of
the H I disk $\simeq 1.0\times 10^{-2}{\rm\ M_{\odot}}$. Taylor et al. (1990)
investigated the H I mass of several planetary nebulae from 21 cm radio
observations. For example, they propose that the planetary nebula BD + 30⋅3639
has the H I mass of $0.028{\rm\ M_{\odot}}$. The random speed of H I medium
$\sim 10{\rm\,km\,s^{-1}}$ corresponds to the thermal speed at $T=5000\rm\ K$,
which is also consistent with the excitation temperature deduced from 21 cm
radio observations.

We use the same set of parameters to show the fitting result for the Raman He
II $\lambda$6545 by the blue line in the right panel. The fit quality of Raman
He II $\lambda$6545 is slightly poorer than of Raman He II at 4851 Å. More
specifically, the simulated profile is shifted blueward of the observed data
by an amount of $\sim 6{\rm\ km\ s^{-1}}$. Because of the severe blending of
Raman He II $\lambda$6545 with N II$\lambda$6548, we consider it a more
appropriate strategy to focus on the Raman He II $\lambda$4851 in the
determination of the best fit parameters pertinent to NGC 6302. In Appendix A,
we show the dependence of spectral profile on parameters, demonstrating the
best fit convincingly.

##  5 Summary and Discussion

Using the archived UVES spectrum of NGC 6302, we have carried out profile
analyses of He II emission lines and Raman-scattered He II at 6545 Å and 4851
Å. Raman He II features are found to be broader and more redshifted than the
hypothetical Raman features that would be formed in a static H I medium. The
Monte Carlo simulation code ’STaRS’ (Chang & Lee, 2020) is used to produce
satisfactory line fits to the Raman-scattered He II features. For our Monte
Carlo approach, the He II emission region is assumed to be compact near the
hot central star surrounded by the H I region with H I column density $N_{\rm
HI}$ = $3\times 10^{21}{\rm\,cm^{-2}}$ and a half opening angle
$\theta_{o}=20^{\circ}$. The kinematics of the H I region is characterized by
the expansion speed $v_{\rm exp}=13{\rm\ km\ s^{-1}}$ and a random speed of
${v_{\rm ran}}=10{\rm\ km\ s^{-1}}$. The physical properties of H I region are
imprinted in Raman He II features via scattering.

In this work, dust particles are assumed to be distributed outside the neutral
medium, presuming that those inside the neutral region are completely
destroyed by strong UV radiation from the central star. However, additional
complications are expected if we introduce a considerable amount of dust in
the H I medium. In a sophisticated model involving a dusty neutral medium, the
formation of Raman He II $\lambda$4851 would be more suppressed than that of
Raman He II $\lambda$6545, because He II $\lambda$972 has the Raman cross
section smaller almost by an order of magnitude than He II $\lambda$1025\. A
line photon of He II $\lambda$972 has to traverse a longer dusty path in order
to yield a Raman optical photon than He II $\lambda$1025 so that dust
extinction is more effective for He II $\lambda$972 than He II $\lambda$1025\.
We defer the formation of Raman features in a dusty neutral medium for future
work.

On the observational side, we expect that Raman-scattered He II feature at
4332 Å can be obtained from very deep spectroscopic observations. Raman He II
$\lambda$4332 was reported in the symbiotic stars, including V1016 Cygni and
RR Telescopii, and also in young planetary nebulae including NGC 7027 (e.g.,
Lee, 2012; van Groningen, 1993; Péquignot et al., 1997). With the future
availability of Raman He II $\lambda$4332, strong constraints on the
scattering geometry and the amount of dust extinction will be placed.

We are grateful to an anonymous referee for constructive comments. This
research has made use of the services of the ESO Science Archive Facility.
H.L. and J.K. were supported by the National Research Foundation of Korea
(NRF) grants funded by the Korea government (NRF-2023R1A2C1006984).

##  Appendix A Monte Carlo Profile Fitting

Figure 7:  Fitting results using Raman He II 4851 (top) & 6545 (bottom)
spectra from STaRS. The gray and black solid lines are observed spectra and
the best-fitting results, respectively. The red-blue solid dashed lines
represent the spectra for smaller-larger parameters ($N_{\rm HI}$,
$\theta_{\rm o}$, $v_{\rm exp}$, and $v_{\rm ran}$). Figure 8:  Simulated
spectra of Raman He II $\lambda$4851 (left) and $\lambda$6545 (right) for
various escaping directions without dust extinction. The colors represent the
different directions $\mu_{z}=\cos\theta_{z}$, where $\theta_{z}$ is the polar
angle measured from the symmetry axis, $z$-axis. Raman-scattered features are
observed to be weaker for $\mu_{z}=0.8-1.0$ than for $\mu_{z}=0.0-0.2$ mainly
due to the anisotropic scattering phase function favoring the forward and
backward scattering. The black lines are identical to the best-fit spectra
shown in Figure 7.

The expanding velocity and random speed of H I are mainly responsible for the
center shift and line broadening of Raman He II, respectively. However, the
line profiles of Raman-scattered He II are determined in a complicated way
involving the multiple scattering effect and sharply varying cross sections as
a function of wavelength in addition to the scattering geometry and the
kinematics. For example, when the covering factor of the neutral scattering
region is significant, far UV He II photons that escape through Rayleigh
scattering may re-enter the scattering region, introducing an additional line
broadening effect to the final Raman-scattered He II (e.g., Choi et al.,
2020). However, the half opening angle of $20^{\circ}$ deduced from Figure 5
is not big enough to enhance the line broadening through the re-entry effect.
For this reason, we consider the random speed component of the H I medium in
our simulation as the main factor affecting the broadened Raman-scattered He
II.

In Figure 7, we show line profiles obtained from our Monte Carlo simulations
adopting the scattering geometry illustrated in Figure 4. In particular, we
show the dependence of the Raman line profiles on the parameters a H I column
density $N_{\rm HI}$, an half opening angle $\theta_{\rm o}$, an expansion
velocity $v_{\rm exp}$, and a random speed $v_{\rm ran}$. The upper and lower
panels show the results for Raman He II $\lambda$4851 and Raman He II
$\lambda$6545, respectively.

In the two left panels, we investigate the effect of $N_{\rm HI}$, where the
three cases for ${N_{\rm HI}}/(10^{21}{\rm\ cm^{-2}})=1,3$ and 5 are shown.
The other parameters are fixed to the best fit values, i.e., ${\theta_{\rm
o}}=20^{\circ}$, ${v_{\rm exp}}=13{\rm\,km\,s^{-1}}$ and ${v_{\rm
ran}}=10{\rm\,km\,s^{-1}}$. Because of the best fit value ${N_{\rm
HI}}=3\times 10^{21}{\rm\ cm^{-2}}$ corresponding to Raman optical depth
exceeding unity for He II $\lambda$1025, the Raman conversion efficiency
increases only slightly compared to that for He II $\lambda$972.

The next two panels show the effect of $\theta_{\rm o}$, where the three
values of ${\theta_{\rm o}}=10^{\circ},20^{\circ}$ and $30^{\circ}$ are
considered. The Raman conversion efficiency is nearly proportional to
$\theta_{\rm o}$, which determines the covering factor of the scattering
region with respect to the He II emission region.

The third two panels show the dependence of $v_{\rm exp}$, which mainly
affects the line center location of the Raman-scattered He II. The line
centers move redward as $v_{\rm exp}$ increases. It is particularly notable
that Raman He II $\lambda$4851 becomes stronger as $v_{\rm exp}$ increases.
This is due to the increase of Raman scattering cross section as He II photons
get redshifted toward the hydrogen resonance, while the Raman optical depth is
less than unity (e.g., Jung & Lee, 2004; Choi et al., 2020). No such
conspicuous increase is seen in the case of Raman He II $\lambda$6545 because
the Raman optical depth exceeds unity.

In the right two panels, the simulated line profiles for three values of
${v_{\rm ran}}=0,10$ and $20{\rm\ km\ s^{-1}}$ are illustrated. The line
profile becomes broader with increasing $v_{\rm ran}$. Because the He II
emission region is assumed to be characterized by a random velocity
$\sigma_{v}=13{\rm\ km\ s^{-1}}$ corresponding to the width of He II emission
lines in Table 2, the resultant profiles are significantly affected by choice
of $v_{\rm ran}$ in the range 0-20 km s-1 . In view of the profile fit
quality, we may safely conclude that the random velocity component in the
neutral region is ${v_{\rm ran}}\simeq 10{\rm\,km\,s^{-1}}$.

We conclude that the best fitting parameters are obtained from the fitting of
Raman He II $\lambda$4851 and the blueward of Raman He II $\lambda$6545\. In
the third two panels, the simulated spectra of Raman He II $\lambda$6545 for
$v_{\rm exp}$ = 13 and 26 km s-1 provide well-fitted features in the blueward
and redward, respectively. We set $v_{\rm exp}$ = 13 km s-1 since the strong
[N II] emission above Raman He II $\lambda$6545 can affect the spectral
profile in the redward of 6545 Å.

Figure 8 shows the simulated line spectra for various values of escaping
direction $\mu_{z}=\cos\theta_{z}$, where $\theta_{z}$ is the angle of the
emergent photon making with the symmetry axis. The spectra for the line of
sight ($\mu_{z}=0$) in our simulation are stronger than those for other
directions with $\mu_{z}\geq 0.8$. This behavior is explained by the fact that
both Raman and Rayleigh processes prefer scattering in the forward and
backward directions to lateral directions, where the phase function of the two
scattering processes is given by

$\Phi(\mu_{s})={3\over 8}(1+\mu_{s}^{2}).$ (A1)

Here, $\mu_{s}={\bf\hat{k}_{i}}\cdot{\bf\hat{k}_{s}}$ is the cosine of the
angle between the wavevectors ${\hat{k}_{i}}$ and ${\hat{k}_{s}}$ of the
incident and scattered photons, respectively (Chang & Lee, 2020).

## References

  * Chang et al. (2015) Chang, S.-J., Heo, J.-E., Di Mille, F., et al. 2015, ApJ, 814, 98, doi: 10.1088/0004-637X/814/2/98
  * Chang & Lee (2020) Chang, S.-J., & Lee, H.-W. 2020, Journal of Korean Astronomical Society, 53, 169, doi: 10.5303/JKAS.2020.53.6.169
  * Choi et al. (2020) Choi, B.-E., Chang, S.-J., Lee, H.-G., & Lee, H.-W. 2020, ApJ, 889, 2, doi: 10.3847/1538-4357/ab61f9
  * Choi & Lee (2020) Choi, B.-E., & Lee, H.-W. 2020, ApJ, 903, L39, doi: 10.3847/2041-8213/abc559
  * Dinh-V-Trung et al. (2008) Dinh-V-Trung, Bujarrabal, V., Castro-Carrizo, A., Lim, J., & Kwok, S. 2008, ApJ, 673, 934, doi: 10.1086/524373
  * Draine (2003) Draine, B. T. 2003, ARA&A, 41, 241, doi: 10.1146/annurev.astro.41.011802.094840
  * Feibelman (2001) Feibelman, W. A. 2001, ApJ, 550, 785, doi: 10.1086/319798
  * Groves et al. (2002) Groves, B., Dopita, M. A., Williams, R. E., & Hua, C.-T. 2002, PASA, 19, 425, doi: 10.1071/AS02010
  * Gussie & Taylor (1995) Gussie, G. T., & Taylor, A. R. 1995, MNRAS, 273, 801, doi: 10.1093/mnras/273.3.801
  * Hyung & Feibelman (2004) Hyung, S., & Feibelman, W. A. 2004, ApJ, 614, 745, doi: 10.1086/423660
  * Jung & Lee (2004) Jung, Y.-C., & Lee, H.-W. 2004, MNRAS, 355, 221, doi: 10.1111/j.1365-2966.2004.08311.x
  * Kang et al. (2009) Kang, E.-H., Lee, B.-C., & Lee, H.-W. 2009, ApJ, 695, 542, doi: 10.1088/0004-637X/695/1/542
  * Kastner et al. (2022) Kastner, J. H., Moraga Baez, P., Balick, B., et al. 2022, ApJ, 927, 100, doi: 10.3847/1538-4357/ac51cd
  * Lee (2012) Lee, H.-W. 2012, ApJ, 750, 127, doi: 10.1088/0004-637X/750/2/127
  * Lee et al. (2006) Lee, H.-W., Jung, Y.-C., Song, I.-O., & Ahn, S.-H. 2006, ApJ, 636, 1045, doi: 10.1086/498143
  * Matsuura et al. (2005) Matsuura, M., Zijlstra, A. A., Molster, F. J., et al. 2005, MNRAS, 359, 383, doi: 10.1111/j.1365-2966.2005.08903.x
  * Meaburn et al. (2008) Meaburn, J., Lloyd, M., Vaytet, N. M. H., & López, J. A. 2008, MNRAS, 385, 269, doi: 10.1111/j.1365-2966.2007.12782.x
  * Nussbaumer et al. (1989) Nussbaumer, H., Schmid, H. M., & Vogel, M. 1989, A&A, 211, L27
  * Osterbrock (1989) Osterbrock, D. E. 1989, Astrophysics of gaseous nebulae and active galactic nuclei
  * Peimbert (1978) Peimbert, M. 1978, in Planetary Nebulae, ed. Y. Terzian, Vol. 76, 215–224
  * Péquignot et al. (1997) Péquignot, D., Baluteau, J. P., Morisset, C., & Boisson, C. 1997, A&A, 323, 217
  * Peretto et al. (2007) Peretto, N., Fuller, G., Zijlstra, A., & Patel, N. 2007, A&A, 473, 207, doi: 10.1051/0004-6361:20066973
  * Pottasch et al. (1985) Pottasch, S. R., Preite-Martinez, A., Olnon, F. M., et al. 1985, A&A, 143, L11
  * Rodriguez & Moran (1982) Rodriguez, L. F., & Moran, J. M. 1982, Nature, 299, 323, doi: 10.1038/299323a0
  * Santander-García et al. (2017) Santander-García, M., Bujarrabal, V., Alcolea, J., et al. 2017, A&A, 597, A27, doi: 10.1051/0004-6361/201629288
  * Schmid (1989) Schmid, H. M. 1989, A&A, 211, L31
  * Storey & Hummer (1995) Storey, P. J., & Hummer, D. G. 1995, MNRAS, 272, 41, doi: 10.1093/mnras/272.1.41
  * Szyszka et al. (2011) Szyszka, C., Zijlstra, A. A., & Walsh, J. R. 2011, MNRAS, 416, 715, doi: 10.1111/j.1365-2966.2011.19087.x
  * Taylor et al. (1990) Taylor, A. R., Gussie, G. T., & Pottasch, S. R. 1990, ApJ, 351, 515, doi: 10.1086/168489
  * van Groningen (1993) van Groningen, E. 1993, MNRAS, 264, 975, doi: 10.1093/mnras/264.4.975
  * Wright et al. (2011) Wright, N. J., Barlow, M. J., Ercolano, B., & Rauch, T. 2011, MNRAS, 418, 370, doi: 10.1111/j.1365-2966.2011.19490.x


11institutetext:  Department of Mechanical and Aerospace Engineering,
University of Central Florida, Orlando FL
<EMAIL_ADDRESS>
Beams, plates, and shells Static elasticity Tribology and mechanical contacts
Structural acoustics and vibration

# Frictional Damping in Biomimetic Scale Beam Oscillations

Hessein Ali 11 Hossein Ebrahimi 11 Ranajay Ghosh 11aa11aa

###### Abstract

Stiff scales adorn the exterior surfaces of fishes, snakes, and many reptiles.
They provide protection from external piercing attacks and control over global
deformation behavior to aid locomotion, slithering, and swimming across a wide
range of environmental condition. In this letter, we investigate the dynamic
behavior of biomimetic scale substrates for further understanding the origins
of the nonlinearity that involve various aspect of scales interaction, sliding
kinematics, interfacial friction, and their combination. Particularly, we
study the vibrational characteristics through an analytical model and
numerical investigations for the case of a simply supported scale covered
beam. Our results reveal for the first time that biomimetic scale beams
exhibit viscous damping behavior even when only Coulomb friction is postulated
for free vibrations. We anticipate and quantify the anisotropy in the damping
behavior with respect to curvature. We also find that unlike static pure
bending where friction increases bending stiffness, a corresponding increase
in natural frequency for the dynamic case does not arise for simply supported
beam. Since both scale geometry, distribution and interfacial properties can
be easily tailored, our study indicates a biomimetic strategy to design
exceptional synthetic materials with tailorable damping behavior.

###### pacs:

46.70.De

###### pacs:

46.25.-y

###### pacs:

46.55.+d

###### pacs:

43.40.+s

Dermal scales are one of the oldest evolutionary adaptation in animals and
found in many organisms of the kingdom Animalia including fishes [1, 2],
reptiles [3, 4] and mammals [5, 6]. Due to their ubiquity and sheer longevity
in evolutionary history, scales have been garnering intense scrutiny from
materials and structures community [7, 8, 9]. From a materials perspective,
characterization studies of various scales have revealed highly intricate
microstructures. Such scales can be hybrid [10, 11, 12], hierarchical [13, 14,
15], and composite in nature [16, 17, 18] capable of engaging multiple length
scales depending on the load [19, 20, 21]. Scales are also naturally
multifunctional [22, 23, 24, 25]. The most apparent function is the protective
properties of scales for the underlying substrate, which has been an
inspiration of ancient armor designs [26]. In addition to protection against
external objects, which is mechanically an indentation type local deformation
problem, global deformation modes such as bending and twisting of a substrate
reveal equally interesting properties. These include reversible nonlinear
stiffening and locking behavior due to the sliding kinematics of the scales in
one-dimensional substrates [27, 28, 29, 30, 31, 32, 33]. In spite of rapidly
accumulating literature on the mechanics of scales and scale covered systems,
their influence on the dynamic behavior of substrates has not been revealed
beyond characterizing strain rate effects under transient loading [34, 35,
36].

In sustained dynamic behavior such as oscillations, the role of scale
interfacial friction can be critical in affecting the long-term response. This
type of behavior is especially apparent for applications such as swimming,
locomotion and structural vibration. Recent research for the static pure
bending case showed that dry friction adds resistance to scale sliding
inevitably increases bending stiffness [28]. Thus there seems to be an
inherent interplay between the stiffness, geometry, range of motion and
friction behavior in these type of systems.

In this letter, we study for the first time the damped oscillations of
biomimetic scale beam using a combination of numerical and analytical
modeling. The work addresses the first mode of bending vibrations.

We first simplify the complex structure of biomimetic scale beam by assuming a
linear-elastic substrate with partially embedded rigid scales on its top
surface Fig. 1a. This is valid because the scales used in this study are
several orders of magnitude stiffer than the substrate material [27, 29] and
typical substrate strains are moderate, with excessive deformations arrested
by additional stiffness induced by the stiff scales. This assumption also
leads to isolation of nonlinearity stemming purely from scales interaction,
rather than material effects. We assume that the substrate is of length
$L_{B}$ and height $h$ and with Young’s modulus $E_{B}=1.5$ MPa, Poisson’s
ratio $\nu=0.42$ [31, 33], and density $\rho_{B}=854kg/m^{3}$. Although these
values are not critical to the outcomes of the study, these are typical of
silicone polymer (Vinylpolysiloxane) used for fabricating many soft
substrates. Each individual scale is of thickness $D$, density $\rho_{s}$ ,
and total length $l_{s}$ that consists of an exposed part $l$ and embedded
part $L$, or $l_{s}=l+L$, Fig. 1a. In addition, we impose $D\ll l_{s}$ and
$L\ll h$, (shallow embedding). This allows each individual scale to be modeled
as a linear torsional spring rotating about fixed point with constant
$\tilde{K}_{s}$ [27, 29, 31]. This constant was analytically computed as
${K}_{s}=E_{B}D^{2}C_{B}(L/D)^{n}$ where $C_{B}$,$n$ are constants with
corresponding values $0.66$ and $1.75$, respectively for same materials
earlier [27]. The ratio of scale length to spacing is denoted by
$\eta=l/\tilde{d}$, where $\tilde{d}$ is the distance between any two
neighboring scales [27]. The scales are embedded onto the substrate at an
initial angle $\theta_{0}$, measured with respect to the beam’s centerline,
and increase nonlinearly once the engagement of scales commences [27, 29, 31].
Due to small strains in the beam substrate, we use Euler-Bernoulli beam
assumptions.

(a)

(b)

Figure 1: (a) Natural history fish image illustrating a typical scale covered
skin [37] and a biomimetic scale covered beam. (b) Schematic diagram of simply
supported biomimetic scale covered beam and RVE geometry.

We also consider uniform scale distribution over the substrate. In this
configuration, scales preserve periodicity post-engagement under pure bending.
However, beyond this limited loading regime, such periodicity would not hold
[31]. This prevents us from using the closed form of kinematic relationship
existing in current literature [27, 28, 29]. However, in the limit of
sufficiently dense distribution of scales we can reasonably impose local
periodicity, Fig. 1b. This assumption leads to a locally uniform pure bending
mode thereby allowing previously developed kinematic relationships to be valid
at least in a local sense. The global deformation of the scaly beam can be
envisioned as a combination of substrate bending and scales rotation [27, 31].
Accordingly, the energy stored in the system is an additive combination of the
two deformation modes. Friction is modeled using Coulomb dry friction with a
coefficient of friction $\mu$.

The local periodicity assumption introduced earlier allows the isolation of a
representative volume element (RVE), Fig. 1b. Thus, locally the kinematic
relationship can be given by
$\theta_{RVE}=\sin^{-1}{(\eta\psi_{RVE}\cos{\psi_{RVE}/2})-\psi_{RVE}/2}$ [28]
where $\psi_{RVE}=\tilde{\kappa}_{RVE}\tilde{d}$ is the angular rotation of
each RVE measured with respect to the beam’s instantaneous curvature
$\tilde{\kappa}_{RVE}$, Fig. 1b.

We derive the dynamic equation of motion by invoking the Hamilton’s principle
of least action leading to the variational energy equation
$\delta\int_{\tilde{t}_{1}}^{\tilde{t}_{2}}{(\hat{T}-\hat{V}+W})d\tilde{t}=0$
for the entire beam. Here $\hat{T}$ is the total kinetic energy per unit
length,$\hat{V}$ is the total strain energy per unit length, while the term
$W$ denotes the work done by the applied traction (see supplemental material
(SM) for their formulas). We apply a harmonic load on the smooth (side
opposite to scales) of form
$\tilde{f}(\tilde{x},\tilde{t})=\tilde{f}_{0}\phi(\tilde{x})\cos{\tilde{\Omega}\tilde{t}}$
to ensure that only the first mode of vibration is activated [38] with simply
supported boundary conditions. Here $\tilde{f}_{0}$ is the load amplitude,
$\phi(\tilde{x})$ is the shape function of the first eigenmode of simply
supported beam [39], while $\tilde{\Omega}$ is the load frequency. This work
neglects the mass of scales (See SM for justification).

$\displaystyle{\rho_{B}}A{{{\partial^{2}}\tilde{y}}\over{\partial{{\tilde{t}}^{2}}}}\;+{E_{B}}{I_{B}}{{{\partial^{4}}\tilde{y}}\over{\partial{{\tilde{x}}^{4}}}}\;+\tilde{C}{{\partial\tilde{y}}\over{\partial\tilde{t}}}+{{{\partial^{2}}}\over{\partial{{\tilde{x}}^{2}}}}{1\over{{N_{RVE}}}}\bigg{[}\mathop{\sum}\limits_{RVE}{{\tilde{K}}_{s}}\left({{\theta_{RVE}}-{\theta_{0}}}\right){{\partial{\theta_{RVE}}}\over{\partial{\psi_{RVE}}}}+$
(1)
$\displaystyle\mathop{\sum}\limits_{RVE}{{\sin\left(\beta\right){{\tilde{K}}_{s}}\left({{\theta_{RVE}}-{\theta_{0}}}\right)\;{\rm{sgn}}\left({\mathop{\dot{\tilde{y}}}\limits}\right)}\over{\cos\left({{\psi_{RVE}}+\beta}\right)-\bar{r}\cos\left(\beta\right)}}\bigg{]}H\left({{{\tilde{\kappa}}_{RVE}}-{{\tilde{\kappa}}_{e}}}\right)=f\left({\tilde{x},\tilde{t}}\right).\vspace{-10pt}$

The Hamilton’s principle leads to the differential equation expressed in Eq.1
(See SM for derivation). Where $A$ refers to the cross sectional area of the
substrate, $I_{B}$ is the second moment of area, $\beta=\tan^{-1}{\mu}$ is the
scales friction angle, and $\bar{r}$ is the distance per scales length
measured from scales base to the interaction with neighboring scale (Fig. 1b).
The Heaviside step function ensures that scales contribution is accounted for
only in the case when the beam is deflected downward and when
$\theta_{RVE}\geq\theta_{0}$. Note that $\tilde{\kappa}_{e}$ refers to the RVE
curvature at the instant $\theta_{RVE}=\theta_{0}$. A small artificial viscous
damping term with constant $\tilde{C}$ is included in the equation to model
any other source of damping and to stabilize numerical integrators for the
forced oscillation cases. The variable $N_{RVE}$ refers to the number of RVEs
utilized in the solution of the system. Note that the results would not depend
on the actual number of RVEs for a sufficiently large number of RVEs. In our
case, we found $N_{RVE}=10$ to be sufficient. In this equation, the terms
inside the summations agree with the quasi-static moment equations expressions
derived earlier in literature [27, 28].

The governing Equation 1 can be written in a non-dimensional form by first
noting that two natural length scales emerge from the geometry, the length of
beam $L_{B}$ (horizontal length scales) and the radius of gyration of the
cross section $R=\sqrt{I_{B}/A}$ (vertical length scales). We denote the ratio
of the two as $\gamma=R/L_{B}$. In addition, a natural time scale for the
problem can be extracted from the natural frequency expression of the
underlying elastic substrate, $\zeta=\sqrt{\rho_{B}AL_{B}^{4}/E_{B}I_{B}}$.
This allows us to non-dimensionalize the following quantities appearing in Eq.
1 $x=\tilde{x}/L_{B}$, $d=\tilde{d}/L_{B}$, $y=\tilde{y}/R$,
$t=\tilde{t}/\zeta$ and $\Omega=\tilde{\Omega}\zeta$. These intrinsic length
and time scales can also be used to obtain the vertical force normalization
unit by noting that Force $\sim MLT^{-2}=\rho_{B}AL_{B}R\zeta^{-2}$. This
leads to a vertical traction normalizer given by
$F^{\prime}=E_{B}I_{B}R/L_{B}^{4}$ and thus
$f_{0}=\tilde{f}_{0}/{F^{\prime}}$. Similar arguments also lead to
normalization of other variables
$C={\tilde{C}L_{B}^{2}\over\sqrt{\rho_{B}AE_{B}I_{B}}}$ and
$K_{s}={\tilde{K}_{s}L_{B}\over E_{B}I_{B}}$. With these non-
dimensionalizations, the non-dimensionalized equation of motion (EOM) is now

$\displaystyle{{{\partial^{2}}y}\over{\partial{t^{2}}}}\;+{{{\partial^{4}}y}\over{\partial{x^{4}}}}\;+C{{\partial
y}\over{\partial
t}}+{{{\partial^{2}}}\over{\partial{x^{2}}\;}}{1\over{\gamma{N_{RVE}}}}\bigg{[}\mathop{\sum}\limits_{RVE}{K_{s}}({{\theta_{RVE}}-{\theta_{0}}})$
$\displaystyle{{\partial{\theta_{RVE}}}\over{\partial{\psi_{RVE}}}}+\mathop{\sum}\limits_{RVE}{{\sin\left(\beta\right){K_{s}}\left({{\theta_{RVE}}-{\theta_{0}}}\right)\;{\rm{sgn}}\left({\mathop{\dot{\tilde{y}}}\limits}\right)}\over{\cos\left({{\psi_{RVE}}+\beta}\right)-\bar{r}\cos\left(\beta\right)}}\bigg{]}$
$\displaystyle
H\left({{\kappa_{RVE}}-{\kappa_{e}}}\right)=f\left({x,t}\right).$ (2)

For simply supported beam, the boundary conditions are $y=0$ and
$\partial^{2}y/\partial x^{2}=0$ at $x=0$ and $1$. Two types of initial
conditions would be studied, displacement and velocity. The displacement type
would be denoted by $y(x,0)=AA(\phi(x))$ where $AA$ is a constant and
$\phi(x)$ is the non-dimensional counterpart of the shape function introduced
earlier. The velocity initial condition is given by $v(x,0)=BB(\phi(x))$ where
$BB$ is a constant. Equation 2 can be solved by utilizing separation of
variables $y(x,t)=\phi(x)T(t)$. Here $T(t)$ is a time-dependent function while
$\phi(x)$ is the first eignmode of a simply supported beam, which is known to
be $\sin{\pi x}$ [39]. Using Bubnov-Galerkin weighted residual method we
obtain the following second order nonlinear ordinary differential equation for
the time-dependent part (See Supplemental Information for detailed derivation)

$\displaystyle\ddot{T}+C\dot{T}+{\pi^{4}}T+{K_{s}}{1\over{{N_{RVE}}}}\bigg{[}\mathop{\sum}\limits_{RVE}{f_{1}}\left({T,\eta,d,{T_{e}},\gamma}\right)+$
(3)
$\displaystyle\mathop{\sum}\limits_{RVE}{f_{2}}\left({T,\eta,d,{T_{e}},\gamma,\mu}\right)\;{\rm{sgn}}\left({{{\dot{T}}_{RVE}}}\right)\bigg{]}H\left({{\kappa_{RVE}}-{\kappa_{e}}}\right)$
$\displaystyle={f_{0}}\cos{\rm{\Omega}}t$

Here, $f_{1}$ and $f_{2}$ are complex nonlinear functions (see SM for their
form). This nonlinear ordinary differential equation is solved using a direct
numerical integrator such as the Newmark beta scheme which nests a Newton-
Raphson nonlinear equation solver [40].

The results of our mathematical model are verified numerically via FE
simulations using a commercially available FE software ABAQUS (Dassault
systems). In the FE model, an assembly of 2D deformable shell parts (for
scales and substrate) was made with rigidity imposed on the embedded scales
and linear-elastic properties were assigned for the substrate. An implicit
dynamic step was then utilized. The initial conditions were imposed through a
predefined velocity in the form $\dot{T}(0)=BB(\phi(x))$ where BB is the
magnitude of velocity. The imposition of contact between scales, mesh
convergence, and boundary condition were done consistent with literature [31].
Interfacial friction was introduced using penalty function formulation with
coefficient $\mu$ [41].

In the small deformation regime, two important geometrical parameters are of
interest in these simulations, $\eta$ and initial angle $\theta_{0}$. We
compare and contrast the effect of these two sources in free vibration regime
$f_{0}=0$ in Eq. 3 in Figs. 2a and 2b. In these simulations, an initial
velocity $BB=10.4$ (chosen to assure engagement of scales for the assigned
$\theta_{0}$) is applied to the beam and the resulting normalized mid-point
deflections are plotted with normalized time. Here the time is normalized by
the undamped period of a plain beam $\tau_{n}=2\pi/w_{n}$ , $w_{n}=\pi^{2}$.
We do not use any artificial viscous damping for free vibration simulations.

We first investigate the effect of geometric parameters on the free vibration
of the biomimetic beam without any friction. We denote the side with scale
engagements as the concave side and the one where scales do not engage as
convex, see Fig. 2a inset. In Fig. 2a, we plot the midpoint deflection of the
biomimetic beam without taking the effect of friction with increasing overlap
ratio $\eta$. There are two visible effects of increasing $\eta$ on the free
vibration regime. The first is the expected increase in the natural frequency
of vibration due to the added stiffness. However, increasing $\eta$ also
decreases the amplitudes of vibration on the concave side, with little effect
in the deflection on the convex side. We expect this behavior since higher
$\eta$ leads to greater energy stored in scales rotation compared to the
substrate. Note that the change in $\eta$ corresponds to varying d and fixing
$l$ for example $\eta=5$ corresponds to $d=0.05$ and $\eta=0$ refers to a
plain beam. Our results are an excellent match with FE simulations, which are
depicted by black dots for $\eta=5$. Therefore, overall, higher overlap ratio
leads to a higher frequency and lower mean amplitude of vibration for same
initial velocity profile. The other geometrical parameter of consequence is
the initial inclination angle $\theta_{0}$. We fix the overlap ratio to
$\eta=5$ and simulate the free vibration response under same initial
conditions as before. As inclination angle increases, the scales engagement is
delayed. Therefore, the overall stiffness of the system is lesser at a given
displacement. This leads to decrease in frequency of response as seen in Fig.
2b. The plots also reveal that the asymmetry in deflection between the scale
side and plain side is made less severe with an increase in scale angle. For
the case with $\theta_{0}=10^{\circ}$, no engagement takes place for the given
initial condition and the vibrational symmetry is fully restored. In summary,
these simulations highlight the critical role played by geometrical parameters
of the scale $\eta$ and $\theta_{0}$ in tailoring the frequency and symmetry
of vibrating biomimetic beam, analogous to deformation in the static regime
[27, 31].

(a)

(b)

Figure 2: (a) The response of biomimetic scale beam under velocity initial
condition for $\theta_{0}=5^{\circ}$ and various $\eta$ (The block dots
represent FE simulations). (b) The effect of $\theta_{0}$ on the behavior of
biomimetic scale beam with $\eta=5$ and under velocity initial condition. For
these plots and all others, the number of RVE was $10$ in the model. Insets
show the concave, neutral, and convex side of deformation.

We now study the interplay of friction with the critical geometric parameters
for free vibrations under same initial conditions. In Fig. 3a, we fix $\eta=5$
, $\theta_{0}=5^{\circ}$ and vary the coefficient of friction. The black dots
for $\eta=5$ in the figure indicate FE simulations, which are in excellent
agreement with our model results. These results highlight several interesting
phenomena. First, as expected we find a pronounced damping behavior with
increasing $\mu$. However, unlike a dry friction spring-mass damper system,
here the amplitude decay is not linear but exponential resembling a viscous
damping behavior. This indicates a geometrically induced regime transformation
of damping from dry friction to viscous drag. We also observe a clear increase
in the time-period of oscillations with increasing $\mu$ as expected in
viscous damped system. The similarities with conventional viscously damped
system ends here. In typical conventional damped oscillators, the decaying
behavior continues till motion is altogether arrested. In contrast, for the
biomimetic system, the decay would not end in a complete stop but rather till
deflection is small enough to preclude scales engagement, returning it to the
frictionless conservative system, vibrating indefinitely. Increasing friction
also decreases the amplitudes of vibration, mirroring the stiffening effect of
scales overlap observed in Fig. 3a, but here, the amplitude reduction occurs
on both sides (concave and convex) of beam and not just on one side. This
apparent distinction between the two stiffening mechanisms has to do with
continuous energy loss brought about in every cycle in the damped oscillation
case, lowering the overall amplitude in contrast to the frictionless case
explored earlier. Furthermore, for this initial velocity oscillation, the
natural frequency decreases with increasnig friction. Interesting this is in
contrast to the stiffening role of friction found in the static case [28].
Thus, in the dynamic case, no additional stiffening effect is apparent. In
spite of friction affecting the amplitude of both sides of vibration, there is
still an observed asymmetry with more rapid reduction occurring on the convex
side. To investigate this asymmetry, we use a logarithmic decrement defined as
$\delta={1\over\Delta}\ln(AA_{n+1}/AA_{n})$ where $\Delta$ is the time between
$n$ and $n+1$ peaks for both sides of oscillations. We then take a ratio of
the two to quantify the asymmetry, with
$\alpha=\delta_{convex}/\delta_{concave}$. This quantity is used to generate
an asymmetry phase map, Fig. 3b which is mapped by $\mu$ and $\eta$. This
phase diagram indicates that anisotropy is considerably accentuated by
increasing both $\eta$ and $\mu$. Specifically, the decays are larger on the
convex side. Also, increasing the friction at lower $\eta$ hardly contributes
to the asymmetry, presumably due to late stage scale engagement. On the other
hand, higher $\eta$ amplifies even small friction coefficients appreciably.
Overall, this study indicates that both distribution of the scales and its
surface properties can be engineered to control the nature of oscillation of
these systems. Another distinguishing feature of the system is the long-term
behavior of damped oscillators. In classical damped systems, the oscillators
eventually cease motion. However, in this case, a damped system first reduces
the amplitude of motion. However, when the displacement is small enough to
prevent further scales engagement, friction ceases to act and the system
continues to oscillate linearly, indefinitely in the ideal case.

(a)

(b)

Figure 3: (a) The response of the scaly beam with $\eta=5$ under velocity
initial condition for various $\mu$ (The block dots represent FE simulations).
(b) phase map of the ratio between the convex and concave damping coefficient
($\alpha$) spanned by $\eta$ and $\mu$. For these plots and all others, the
number of RVE was $10$ in the model.

Next, we turn our attention to the role of initial conditions in influencing
the time-period of oscillations. We now give an initial displacement condition
to the beam which is varied and corresponding time-period of oscillation
calculated. We find that this variation leads to different time-periods for
the system. We undertake two different studies to investigate the role of
$\eta$ and $\mu$ on this phenomena. In Figure.4a, frictionless biomimetic
scale beams are studied for various values of $\eta$. For a plain beam
($\eta$=0), there is a complete insensitivity of time-period on initial
displacement, as expected. Note that we have ignored nonlinearity from large
deformation or material sources in this case. As scales with higher $\eta$ are
introduced, for a given initial displacement, the time period of oscillations
are smaller for higher $\eta$. At the same time, we also observe that as
$\eta$ increases, there is dramatic increase in time-period sensitivity to
initial amplitude for a given beam, leading to pronounced decrease of time-
period with increasing initial amplitude. However, this trend is gradually
arrested as $\eta$ increases further since the system now begins to oscillate
near its locking limit [28]. This leads to restoration of insensitivity.
Furthermore, this insensitivity restoration phenomena can also be introduced
by increasing the coefficient of friction for a given $\eta$. This can be seen
in Fig. 4b where we plot the dynamic response of a biomimetic beam with
$\eta=10$ for different $\mu$. This trend continues till $\mu$ is large enough
to transition the system into a really small deformation regime oscillation
due to proximity to frictional locking curvature ($AA>20$ and $\mu=2$ ). Due
to very high friction, the system quickly attenuates, returning to a linear
stage due to lack of scale engagements. This results in almost no sensitivity
to initial displacement. Overall, for a given initial condition, scales with
higher $\mu$ have higher time periods. This is consistent with the results
obtained for initial velocity simulations, 3b.

(a)

(b)

Figure 4: (a) Period of free vibration of frictionless biomimetic scale beam
for various initial displacements $AA$ and $\eta$. (b) Period of free
vibration of frictional biomimetic scale beam for various initial
displacements $AA$ and $\mu$. $\theta_{0}=0^{\circ}$ for both these
simulations.

At this point it is interesting to compare this behavior with that of
hardening Duffing oscillator which also exhibits this characteristic amplitude
dependent time period [42]. Although not exactly a Duffing system, we have in
effect obtained a tailorable hardening spring system.

Next, we probe the forced-vibration behavior of this system by obtaining the
amplitude-frequency response. We plot the maximum mid-point deflection of the
beam on the convex side with forcing frequency. We first isolate the effect of
$\eta$ by neglecting interfacial friction. However, a small artificial
viscosity $C=0.5$ (which amounts to a damping ratio of $\approx 0.025$) is
used to encourage convergence. Without loss of generality, the forcing
amplitude is taken to be $f_{0}=200$ and initial scale inclination
$\theta_{0}=5^{\circ}$. The results of this simulation are shown in Fig. 5a in
which we plot the amplitude versus the non-dimensional applied
frequency(normalized by the natural frequency of a plain beam). In this
figure, we get a pronounced backbone curve, similar to a hardening nonlinear
oscillator. We observe that increasing $\eta$ lead to a greater flip over
behavior indicating a stronger nonlinearity, taking the resonance farther away
from that of a plain beam. Next, in Fig. 5b we investigate the effect of
friction. We take $\eta=5$ and vary $\mu$. As in the case of free vibration,
the role of friction is found to be mainly that of dampener in forced
vibration response. The resonant frequency is therefore moved closer to the
natural frequency of the substrate with increasing friction.

(a)

(b)

Figure 5: (a) Amplitude-frequency behavior at steady state for frictionless
biomimetic scale beam with $\theta_{0}=5^{\circ}$ and various $\eta$. (b) The
amplitude-frequency response at steady state for a frictional biomimetic scale
beam with $\eta=5$ and various $\mu$.

In conclusion, although previous static deflection studies had highlighted the
dual role of friction as a dissipater as well as stiffness enhancer, the
dynamic studies revealed phenomena that are even more interesting. We found
that in spite of postulating only dry Coulomb friction, the oscillations
resemble viscous damping behavior. The scales provide significant anisotropy
of frictional behavior, which can be tailored using interfacial and geometric
parameters. We have also found various degrees of sensitivity to initial
displacements depending upon scale geometry and friction. Unlike the static
case, no appreciable stiffening effect of friction on the dynamics was
observed. This was also confirmed in forced vibration studies. Therefore, our
study indicates that biomimetic scale beams can be an excellent platform for
synthesizing substrates with tailorable damping indices.

###### Acknowledgements.

This work was supported by the United States National Science Foundation’s
Civil, Mechanical, and Manufacturing Innovation Grant No. 1825801.

## References

  * [1] Dou Z., Wang J. Chen D. Journal of bionic Engineering92012457.
  * [2] Onozato H. Watabe N. Cell and tissue research2011979409.
  * [3] Chang C., Wu P., Baker R. E., Maini P. K., Alibardi L. Chuong C.-M. The International journal of developmental biology532009813.
  * [4] Di-Poï N. Milinkovitch M. C. Science advances22016e1600708.
  * [5] Chen I. H., Kiang J. H., Correa V., Lopez M. I., Chen P.-Y., McKittrick J. Meyers M. A. Journal of the mechanical behavior of biomedical materials42011713.
  * [6] Wang B., Yang W., Sherman V. R. Meyers M. A. Acta biomaterialia41201660.
  * [7] Ehrlich H. Materials design principles of fish scales and armor in Biological Materials of Marine Origin (Springer) 2015 pp. 237–262.
  * [8] Funk N., Vera M., Szewciw L. J., Barthelat F., Stoykovich M. P. Vernerey F. J. ACS applied materials & interfaces720155972.
  * [9] Martini R. Barthelat F. Bioinspiration & biomimetics112016066001.
  * [10] Ning G., Li T., Yan J., Xu C., Wei T. Fan Z. Carbon542013241.
  * [11] Rudykh S., Ortiz C. Boyce M. C. Soft Matter1120152547.
  * [12] Drelich A. J., Monteiro S. N., Brookins J. Drelich J. W. Advanced Biosystems220181800055.
  * [13] Fratzl P. Weinkamer R. Progress in materials Science5220071263.
  * [14] Lakes R. Nature3611993511.
  * [15] Qing H. Mishnaevsky Jr L. Mechanics of Materials4120091034.
  * [16] Buehler M. J. Proceedings of the National Academy of Sciences103200612285.
  * [17] Espinosa H. D., Juster A. L., Latourte F. J., Loh O. Y., Gregoire D. Zavattieri P. D. Nature communications22011173.
  * [18] Ortiz C. Boyce M. C. Science31920081053.
  * [19] Song J., Ortiz C. Boyce M. C. Journal of the mechanical behavior of biomedical materials42011699.
  * [20] Yang W., Gludovatz B., Zimmermann E. A., Bale H. A., Ritchie R. O. Meyers M. A. Acta biomaterialia920135876.
  * [21] Zhu D., Szewciw L., Vernerey F. Barthelat F. Journal of the mechanical behavior of biomedical materials24201330.
  * [22] Colbert E. H. et al. Evolution of the vertebrates. A history of the backboned animals through time.1955.
  * [23] Sire J.-Y., Donoghue P. C. Vickaryous M. K. Journal of Anatomy2142009409.
  * [24] Vickaryous M. K. Sire J.-Y. Journal of Anatomy2142009441.
  * [25] Willmer P. Invertebrate relationships: patterns in animal evolution (Cambridge University Press) 1990.
  * [26] White Z. W. Vernerey F. J. Bioinspiration & biomimetics132018041004.
  * [27] Ghosh R., Ebrahimi H. Vaziri A. Applied Physics Letters1052014233701.
  * [28] Ghosh R., Ebrahimi H. Vaziri A. EPL (Europhysics Letters)113201634003.
  * [29] Vernerey F. J. Barthelat F. International Journal of Solids and Structures4720102268.
  * [30] Vernerey F. J., Musiket K. Barthelat F. International Journal of Solids and Structures512014274.
  * [31] Ali H., Ebrahimi H. Ghosh R. International Journal of Solids and Structures2019.
  * [32] Horton R., Ebrahimi H., Ali H., Galvez J., Gordon A. et al. 2018\.
  * [33] Ali H., Ebrahimi H. Ghosh R. arXiv preprint arXiv:1810.106292018.
  * [34] Lin Y., Wei C., Olevsky E. Meyers M. A. Journal of the mechanical behavior of biomedical materials420111145.
  * [35] Ghods S., Murcia S., Ossa E. Arola D. Journal of the mechanical behavior of biomedical materials902019451.
  * [36] Arola D., Ghods S., Son C., Murcia S. Ossa E. Journal of the Royal Society Interface16201920180775.
  * [37] See source at https://thegraphicsfairy.com/natural-history-fish-image/.
  * [38] Abhyankar N., Hall E. Hanagud S. Journal of Applied Mechanics601993167.
  * [39] Rao S. S. Vibration of continuous systems (Wiley) 2019.
  * [40] Ahmadian Saieni H. Non-linear vibrations of tensegrity structures (2012).
  * [41] Abaqus C. Analysis user’s manual, version 6.12 (2012).
  * [42] Kovacic I. N. Rand R. H. Duffing-type oscillators with amplitude-independent period in Applied Non-Linear Dynamical Systems (Springer) 2014 pp. 1–10.