DAMIC-M Collaboration

# Search for Daily Modulation of MeV Dark Matter Signals with DAMIC-M

I. Arnquist  Pacific Northwest National Laboratory (PNNL), Richland, WA,
United States  N. Avalos  Centro Atómico Bariloche and Instituto Balseiro,
Comisión Nacional de Energía Atómica (CNEA), Consejo Nacional de
Investigaciones Científicas y Técnicas (CONICET), Universidad Nacional de Cuyo
(UNCUYO), San Carlos de Bariloche, Argentina  D. Baxter  Kavli Institute for
Cosmological Physics and The Enrico Fermi Institute, The University of
Chicago, Chicago, IL, United States  X. Bertou  Centro Atómico Bariloche and
Instituto Balseiro, Comisión Nacional de Energía Atómica (CNEA), Consejo
Nacional de Investigaciones Científicas y Técnicas (CONICET), Universidad
Nacional de Cuyo (UNCUYO), San Carlos de Bariloche, Argentina  N. Castelló-Mor
Instituto de Física de Cantabria (IFCA), CSIC - Universidad de Cantabria,
Santander, Spain  A.E. Chavarria  Center for Experimental Nuclear Physics and
Astrophysics, University of Washington, Seattle, WA, United States  J. Cuevas-
Zepeda  Kavli Institute for Cosmological Physics and The Enrico Fermi
Institute, The University of Chicago, Chicago, IL, United States  A.
Dastgheibi-Fard  LPSC LSM, CNRS/IN2P3, Université Grenoble-Alpes, Grenoble,
France  C. De Dominicis  Laboratoire de physique nucléaire et des hautes
énergies (LPNHE), Sorbonne Université, Université Paris Cité, CNRS/IN2P3,
Paris, France  O. Deligny  CNRS/IN2P3, IJCLab, Université Paris-Saclay, Orsay,
France  J. Duarte-Campderros  Instituto de Física de Cantabria (IFCA), CSIC -
Universidad de Cantabria, Santander, Spain  E. Estrada  Centro Atómico
Bariloche and Instituto Balseiro, Comisión Nacional de Energía Atómica (CNEA),
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET),
Universidad Nacional de Cuyo (UNCUYO), San Carlos de Bariloche, Argentina  N.
Gadola  Universität Zürich Physik Institut, Zürich, Switzerland  R. Gaïor
Laboratoire de physique nucléaire et des hautes énergies (LPNHE), Sorbonne
Université, Université Paris Cité, CNRS/IN2P3, Paris, France  T. Hossbach
Pacific Northwest National Laboratory (PNNL), Richland, WA, United States  L.
Iddir  Laboratoire de physique nucléaire et des hautes énergies (LPNHE),
Sorbonne Université, Université Paris Cité, CNRS/IN2P3, Paris, France  B. J.
Kavanagh  Instituto de Física de Cantabria (IFCA), CSIC - Universidad de
Cantabria, Santander, Spain  B. Kilminster  Universität Zürich Physik
Institut, Zürich, Switzerland  A. Lantero-Barreda  Instituto de Física de
Cantabria (IFCA), CSIC - Universidad de Cantabria, Santander, Spain  I. Lawson
SNOLAB, Lively, ON, Canada  S. Lee  Universität Zürich Physik Institut,
Zürich, Switzerland  A. Letessier-Selvon  Laboratoire de physique nucléaire et
des hautes énergies (LPNHE), Sorbonne Université, Université Paris Cité,
CNRS/IN2P3, Paris, France  P. Loaiza  CNRS/IN2P3, IJCLab, Université Paris-
Saclay, Orsay, France  A. Lopez-Virto  Instituto de Física de Cantabria
(IFCA), CSIC - Universidad de Cantabria, Santander, Spain  K. J. McGuire
Center for Experimental Nuclear Physics and Astrophysics, University of
Washington, Seattle, WA, United States  P. Mitra  Center for Experimental
Nuclear Physics and Astrophysics, University of Washington, Seattle, WA,
United States  S. Munagavalasa  Kavli Institute for Cosmological Physics and
The Enrico Fermi Institute, The University of Chicago, Chicago, IL, United
States  D. Norcini  Kavli Institute for Cosmological Physics and The Enrico
Fermi Institute, The University of Chicago, Chicago, IL, United States
Department of Physics and Astronomy, Johns Hopkins University, Baltimore, MD,
United States  S. Paul  Kavli Institute for Cosmological Physics and The
Enrico Fermi Institute, The University of Chicago, Chicago, IL, United States
A. Piers  Center for Experimental Nuclear Physics and Astrophysics, University
of Washington, Seattle, WA, United States  P. Privitera  Kavli Institute for
Cosmological Physics and The Enrico Fermi Institute, The University of
Chicago, Chicago, IL, United States  Laboratoire de physique nucléaire et des
hautes énergies (LPNHE), Sorbonne Université, Université Paris Cité,
CNRS/IN2P3, Paris, France  P. Robmann  Universität Zürich Physik Institut,
Zürich, Switzerland  S. Scorza  LPSC LSM, CNRS/IN2P3, Université Grenoble-
Alpes, Grenoble, France  M. Settimo  SUBATECH, Nantes Université, IMT
Atlantique, CNRS/IN2P3, Nantes, France  R. Smida  Kavli Institute for
Cosmological Physics and The Enrico Fermi Institute, The University of
Chicago, Chicago, IL, United States  M. Traina  Center for Experimental
Nuclear Physics and Astrophysics, University of Washington, Seattle, WA,
United States  Laboratoire de physique nucléaire et des hautes énergies
(LPNHE), Sorbonne Université, Université Paris Cité, CNRS/IN2P3, Paris, France
R. Vilar  Instituto de Física de Cantabria (IFCA), CSIC - Universidad de
Cantabria, Santander, Spain  G. Warot  LPSC LSM, CNRS/IN2P3, Université
Grenoble-Alpes, Grenoble, France  R. Yajur  Kavli Institute for Cosmological
Physics and The Enrico Fermi Institute, The University of Chicago, Chicago,
IL, United States  J-P. Zopounidis  Laboratoire de physique nucléaire et des
hautes énergies (LPNHE), Sorbonne Université, Université Paris Cité,
CNRS/IN2P3, Paris, France

###### Abstract

Dark Matter (DM) particles with sufficiently large cross sections may scatter
as they travel through Earth’s bulk. The corresponding changes in the DM flux
give rise to a characteristic daily modulation signal in detectors sensitive
to DM-electron interactions. Here, we report results obtained from the first
underground operation of the DAMIC-M prototype detector searching for such a
signal from DM with MeV-scale mass. A model-independent analysis finds no
modulation in the rate of 1$e^{-}$ events with periods in the range 1-48 h. We
then use these data to place exclusion limits on DM in the mass range [0.53,
2.7] MeV/c2 interacting with electrons via a dark photon mediator. Taking
advantage of the time-dependent signal we improve by $\sim$2 orders of
magnitude on our previous limit obtained from the total rate of 1$e^{-}$
events, using the same data set. This daily modulation search represents the
current strongest limit on DM-electron scattering via ultralight mediators for
DM masses around 1 MeV/c2.

DAMIC-M, CCD, Dark Matter, Daily Modulation, Dark Current, DM-electron
scattering, Dark Photon

Introduction.— There are numerous astrophysical and cosmological observations
supporting the existence of a non-baryonic and non-luminous form of matter
[1], known as dark matter (DM). Despite a wealth of experimental efforts to
detect DM through its non-gravitational interactions [2, 3, 4], its identity
remains unknown. DM particles with MeV-scale masses are viable candidates that
appear in a range of natural, well-motivated scenarios [5, 6, 7, 6, 8, 9].
However, their discovery remains challenging since such light particles do not
have sufficient kinetic energy to produce detectable nuclear recoils in
direct-detection experiments [10]. Instead, DM-electron scattering would
produce eV-scale electronic recoils, which would be observed in semiconductor-
based detectors with $\mathcal{O}$(1 eV) ionization thresholds [11]. To date,
the sensitivity of semiconductor searches [12, 13, 14, 15, 16] has been partly
limited by background levels [17]. A time-dependent DM signal detected above a
time-independent background can significantly increase the sensitivity of
these searches.

Various phenomena give rise to daily-modulated DM signals in direct searches,
including the variation of the velocity of the lab frame due to the Earth’s
rotation [18], gravitational focusing of the DM flux by the Earth [19, 20,
18], and anisotropy in the response of the target material [21, 22, 23, 24,
25]. Furthermore, if the DM scattering cross-section is sufficiently large,
the velocity distribution of DM particles at the detector may be distorted by
their interactions in the Earth [26, 27, 28]. Over the course of a sidereal
day, the position of the detector rotates with respect to the incoming DM
flux. Thus, the DM particles will travel a greater or smaller distance across
the Earth at different times, leading to a daily modulation.

Most experimental searches for daily modulations [29, 30, 31, 32] have
targeted DM particles interacting with nuclei. Here, we focus on models with
MeV-scale DM, which may be coupled to Standard Model particles via a
kinetically-mixed dark photon [33, 34], and thus interact with electrons. For
DM-electron scattering cross-sections currently within reach of semiconductor
searches, interactions in the Earth may be substantial [35], which motivate
the search for a daily modulated signal on top of a constant background [36,
37].

In this letter, we present a search for a daily modulated signal with a
prototype of the DAMIC-M (Dark Matter in CCDs at Modane) experiment [38]. We
first perform a model-independent search for modulations over a wide range of
periods. We then report limits on DM particles with MeV-scale masses
interacting with electrons, which incorporate, for the first time with an
underground detector, a time-dependent analysis to provide a powerful
discriminant between signal and background.

Setup and data.—The DAMIC-M experiment will feature $\sim$200 skipper charge-
coupled devices (CCDs), for a target silicon mass of $\sim$700 g, installed
under the French Alps at the Laboratoire Souterrain de Modane (LSM). DAMIC-M
is sensitive to ionization signals of just a few electrons thanks to the
extremely low level of dark current and the sub-electron readout noise of the
devices. By operating prototype CCDs at LSM in the Low Background Chamber
(LBC) we obtained world-leading limits on sub-GeV DM particles interacting
with electrons [16]. The daily modulation results reported here are obtained
with the same apparatus and data set, which are summarized in the following.

Two high-resistivity, ($>$10 k$\Omega$ cm) $n$-type silicon CCDs [39, 40, 41]
are installed in the LBC. Each CCD has $6144\times 4128$ pixels (pixel area
$15\times 15$ $\mu$m2) and is $670$ $\mu$m thick, for a mass of $\sim$$9$ g.
The detectors are mounted in a high-purity copper box, which also shields
infrared radiation. A lead shield of $\sim$$7.5$ cm thickness, with the
innermost 2 cm of ancient origin, encloses the copper box inside the LBC
cryostat. An additional 15 cm of lead and 20 cm of high density polyethylene
surround the cryostat to shield from environmental radiation. The measured
background is $\sim$$10$ events/keV/kg/day, consistent with expectations. The
CCDs are kept at a temperature of $\sim 130$ K inside the vacuum cryostat at
$\sim$$5\times 10^{-6}$ mbar.

Charge carriers produced by ionizing radiation in the silicon bulk of the CCD
drift towards the $x$-$y$ plane where pixels are located. Thermal diffusion
while drifting results in a spatial variance, $\mathrm{\sigma^{2}_{xy}}$, of
the charge collected in contiguous pixels [42]. In the readout process, the
pixel charge is moved across the array by appropriate voltage clocking and
read out serially by two amplifiers (referred as U and L) located at the
corners of the CCD at the end of the serial register. These skipper amplifiers
[43, 44, 45] allow for multiple non-destructive measurements (NDCMs) of the
pixel charge. Sub-electron resolution is achieved by averaging a sufficient
number ($N_{\mathrm{skips}}$) of NDCMs, with the resolution improving as
$1/\sqrt{N_{\mathrm{skips}}}$. The performance of DAMIC-M CCDs is detailed in
Ref. [46].

For this analysis we use 63 days of uninterrupted data-taking beginning on
June 8th 2022 (the SR2 data set in Ref. [16]). Relevant CCD operating
parameters include $N_{\mathrm{skips}}=650$, corresponding to a pixel charge
resolution of $\sim$$0.2\,e^{-}$, and a $\mathrm{10\times 10}$ pixel binning.
In subsequent text, the term pixel is used to describe a $10\times 10$ bin of
pixels from which charge is summed before being read out. To minimize dark
current counts, only a fraction of the CCD is read out ($640\times 110$
pixels), and its charge is cleared between consecutive images.

Image data reduction and pixel selection criteria are detailed in Ref. [16].
In brief, images are first calibrated from the fitted position of the 0$e^{-}$
and 1$e^{-}$ peaks in the pixel charge distribution. Then contiguous pixels
with charge are joined into clusters. Those with total charge $>$7$e^{-}$ are
excluded from further analysis since they are unlikely to be produced by sub-
GeV DM interactions with electrons. For each cluster we also exclude 10
trailing pixels in the horizontal and vertical directions to account for
charge transfer inefficiencies. Defects in the CCDs may release charge during
the readout, appearing in the images as “hot” pixels and columns [47]. These
are identified by their increased rate of pixels with 1$e^{-}$, and rejected.
Lastly, we exclude one of the prototype CCDs due to the presence of several
charge traps in the serial register (identified by a decreased rate of pixels
with 1$e^{-}$). A portion of the remaining CCD is also excluded for the same
reason. The final integrated exposure after applying these selection criteria
is 39.97 g-days, distributed over $N_{\mathrm{im}}=8779$ images.

For each image $i$ the charge distribution of the selected pixels is fitted
with the sum of two Gaussian functions corresponding to the 0$e^{-}$ and
1$e^{-}$ peaks. The number of pixels with 1$e^{-}$ charge, $N_{1}^{i}$, and
its uncertainty $\sigma_{N_{1}}^{i}$ are then obtained from the fitted
parameters. The U and L sides of the image are fitted independently to account
for small differences in the calibration and resolution of the amplifiers. We
find the value of $N_{1}^{i}$ to be $\sim$80 and $\sim$30 pixels/image for the
U and L amplifier, respectively. The difference in counts is due to a large
portion of the L side being excluded by the serial register trap criteria.
Corresponding rates of 1$e^{-}$, $R_{1}^{i}$, are then obtained dividing
$N_{1}^{i}$ by the effective image exposure time ($t_{\mathrm{exp}}=0.00356$
days) and target mass after pixel selection. For each image the relative
uncertainty $\sigma_{R_{1}}^{i}/R_{1}^{i}$ is $\sim$$13(21)\%$ for the U (L)
side. This set of 1$e^{-}$ rates, measured every $\sim$10 min for 63
consecutive days, is used to search for a daily modulation.

The rate $R_{1}^{i}$ for the U side is shown as a function of time in Fig. 1.
The rate is well-parameterized by a constant plus an exponential with time
decay constant $\tau\sim 33.5$ days. The L side has the same behavior. Since
$\tau\gg$ 1 day, a daily modulation analysis is largely unaffected. The
observed time dependence is consistent with the dark current stabilizing over
time, a characteristic of these types of devices.

Figure 1: Measured rate of 1$e^{-}$, $R_{1}^{i}$, as a function of time for
the U-side. Each data point is the rate averaged over 50 images. The red line
is the result of the fit of a constant plus an exponential.

Model-independent search for a daily modulation.—We first perform a model-
independent analysis by fitting the measured rates $R_{1}^{i}$ to a function
$F(t)=Be^{-t/\tau}+C+A\cos((2\pi(t-\phi)/T)$ where $B,C~{}$ and $\tau$
parameterize the exponential time dependence observed in the data, and
$A,~{}T$ and $\phi$ are the modulation amplitude, period and phase. This
simplified model approximates a periodic signal with the first term of its
Fourier expansion. We use a binned likelihood method to test the null
hypothesis, $i.e.$ absence of modulation, corresponding to $A=0$. The
likelihood function is

$\mathcal{L}(\theta)=\prod_{i=1}^{N_{\mathrm{im}}}\sqrt{2\pi}\sigma_{R_{1}}^{i}\exp\left\\{-\frac{1}{2}\left(\frac{R_{1}^{i}-F(t_{i}|\theta)}{\sigma_{R_{1}}^{i}}\right)^{2}\right\\}\,,$
(1)

where $\theta=\\{A,B,C,\tau,T,\phi\\}$ is the set of model parameters and
$t_{i}$ is the solar time when image $i$ was taken. The spread in the arrival
time of the events within each image has negligible impact for the range of
periods explored. The likelihood of the no-modulation hypothesis $H_{0}$ is
found by maximizing Eq. (1) under the constraint $A=0$,
$\mathcal{L}_{H_{0}}=\sup_{A=0}\mathcal{L(\theta)}$, while the global maximum
for a given period $T$ is
$\mathcal{L}_{H_{1}}=\sup_{\theta}\mathcal{L}(\theta)$. The likelihood-ratio
test statistic for the $H_{0}$ hypothesis is given by
$t_{q}=-2\ln\left(\mathcal{L}_{H_{0}}/\mathcal{L}_{H_{1}}\right)$. The local
significance of a departure from the $H_{0}$ hypothesis for any given period
$T$ can then be quantified by the value of $t_{q}$ [48].

Figure 2:  Model-independent search for daily modulation: the test statistic
$t_{q}$ and corresponding local significance as a function of the modulation
period $T$ in a range around 24 h; the global significance levels for 1, 2 and
3 $\sigma$ are also shown as horizontal lines. Both local and global
significances are one-sided.

In Fig. 2 we show the value of the test statistic and the corresponding local
significance as a function of $T$ in a range around 24 h. Global
significances, calculated by examining the distribution of the maximum value
of the test statistic under $H_{0}$ [49], are also shown in Fig. 2. We find
that the measured rate is consistent with the null hypothesis (no-modulation)
for periods close to 24 h, as well as over the wider range 1-48 h (not shown
in Fig. 2). This null result indicates the lack of a discernible DM signal,
and demonstrates the stability of our detectors over a time scale of days.

DM-$e^{-}$ scattering.—In this section, we introduce a new approach for
constraining DM-$e^{-}$ scattering, which exploits the expected daily
modulation due to DM interactions within the Earth. The rate of DM-$e^{-}$
scattering in a semiconductor is obtained by convolving the scattering cross-
section with the DM velocity distribution at the detector $f(\mathbf{v},t)$
and with the crystal form factor $f_{\mathrm{c}}$ [11, 50, 51]:

$\displaystyle\begin{split}\frac{\mathrm{d}R}{\mathrm{d}E_{e}}\propto\bar{\sigma}_{e}\int\frac{\mathrm{d}q}{q^{2}}\left[\int\frac{f(\mathbf{v},t)}{\mathrm{v}}\,\mathrm{d^{3}v}\right]\left|F_{\mathrm{DM}}(q)\right|^{2}\left|f_{\mathrm{c}}(q,E_{e})\right|^{2}\,.\end{split}$
(2)

Here, $\bar{\sigma}_{e}$ is a reference DM-$e^{-}$ scattering cross-section,
$E_{e}$ is the energy deposited as ionization and $q$ is the momentum
transfer. We consider a model of DM which couples to Standard Model particles
via a kinetically mixed dark photon. In this case, the DM form factor
$F_{\mathrm{DM}}=(\alpha m_{e}/q)^{n}$ depends on the mass of the mediator
$m_{A^{\prime}}$, with $n=0$ for a heavy mediator
($m_{A^{\prime}}\gg\alpha\,m_{e}$) and $n=2$ for ultralight mediators
($m_{A^{\prime}}\ll\alpha\,m_{e}$) [35].

The velocity distribution of DM is expected to follow a Maxwell-Boltzmann
distribution in the Galactic frame [52]111We use the Standard Halo Model [76]
parameters with the Earth’s velocity fixed to its average value during the
data-taking period, $v_{E}=263$ km/s.. The velocity distribution at the
detector is obtained by boosting into the rest-frame of the laboratory, which
is moving at a velocity
$\mathbf{v}_{\mathrm{lab}}(t)=\mathbf{v}_{\odot}+\mathbf{v}_{\mathrm{E}}(t)+\mathbf{v}_{\mathrm{rot}}(t)$.
Here, $\mathbf{v}_{\odot}$ is the velocity of the Sun in the Galactic frame,
$\mathbf{v}_{\mathrm{E}}(t)$ the Earth’s velocity in the Solar frame [54], and
$\mathbf{v}_{\mathrm{rot}}(t)$ is the rotational velocity of the lab in an
Earth-centered inertial frame [55, App. A]. The velocity distribution
$f(\mathbf{v},t)$ at the detector is therefore always time-dependent due to
the Earth’s motion.

If the DM-scattering cross section is sufficiently large, however, the flux of
DM may also be distorted by scattering of particles in the Earth before they
reach the detector. The velocity distribution then depends on the isodetection
angle $\Theta$ [55], defined as the angle between the local zenith and the
direction of the mean DM flux
$\langle\mathbf{v}_{\chi}\rangle=-\mathbf{v}_{\mathrm{lab}}$. For
$\Theta=0^{\circ}$, the mean DM flux comes from directly above the detector,
while at $\Theta=180^{\circ}$, it comes from directly below. The isodetection
angle oscillates over a sidereal day with the Earth’s rotation, and thus the
expected DM flux. During the data-taking period the isodetection angle at LSM
($45.2^{\circ}$ N) varied in the range $\Theta\in[0^{\circ},92^{\circ}]$.

Depending on the regime of interest, a number of formalisms have been
developed to estimate the DM flux at the detector as a function of $\Theta$,
taking into account Earth-scattering (e.g. [56, 36, 57, 58, 55, 59, 60]).
Here, we employ a modified version of the Verne code [61, 62]222Note that
Verne uses an alternative convention for the isodetection angle
$\gamma=180^{\circ}-\Theta$., which assumes that light DM particles travel
along straight-line trajectories and either continue unaffected or are
reflected back along their incoming path when they scatter [64]. DM scattering
with Earth nuclei, implemented in Verne with charge-screening [35], is the
dominant process for daily modulation effects. Since each $t_{i}$ corresponds
to a $\Theta(t_{i})$ 333Note that the isodetection angle corresponding to a
given time in a sidereal day changes by several degrees during the data taking
period, this code allows us to efficiently calculate $f(\mathbf{v},t_{i})$
over a wide range of parameter space.

We search for a daily modulation of the rate $R_{1}^{i}$ by performing an
unbinned likelihood fit with Eq. (1), where $F$ includes both the background
and signal model. For the background we use a Poisson model,
$\textup{Pois}(n_{q}|\lambda)$, which gives the probability of measuring
$n_{q}$ charges in a pixel given a dark current $\lambda$. As shown in Fig. 1,
the dark current in our data stabilizes with an exponential decrease in time.
Thus, we model the time-dependent dark current as
$\lambda(t_{i})=\lambda_{A}\cdot\exp(-t_{i}/\tau)+\lambda_{\mathrm{eq}}$,
where $\lambda_{A}$ and $\tau$ are the amplitude and time constant of the
exponential decay of the dark current, and $\lambda_{\mathrm{eq}}$ is its
value once stabilized. The signal model,
$S(n_{q}|m_{\chi},\bar{\sigma}_{e},t_{i})$, gives the probability that $n_{q}$
charges are measured in a pixel for DM particles of mass $m_{\chi}$ and cross
section $\bar{\sigma}_{e}$; the time dependence in $S$ accounts for the
expected daily modulation. The fitting function in Eq. (1) is then:

$F(t_{i}|\theta)=\frac{1}{t_{\mathrm{exp}}m_{\mathrm{pix}}}\sum_{j=0}^{1}\textup{Pois}(1-j|\lambda(t_{i}))S(j|m_{\chi},\bar{\sigma}_{e},t_{i})\,,$
(3)

where $m_{\mathrm{pix}}$ is the mass of one CCD pixel and
$\theta=\\{\lambda_{A},\tau,\lambda_{\mathrm{eq}},\bar{\sigma}_{e}\\}$ are
free parameters in the fit.

Figure 3:  U-side residual rate after the subtraction of the background model,
binned as a function of local apparent sidereal time. As a reference, the
upper $x$-axis gives the isodetection angle $\Theta(t)$ for the first day of
data taking. Each data point is obtained from the average of $\sim$ 300
images. The expected signal for a DM particle of mass 1 MeV/c2,
$\bar{\sigma}_{e}=10^{-32}$ cm2 interacting via an ultralight dark photon
mediator is superimposed (light red; the band reflects the time shift and
amplitude change of the signal during the data taking period.)

Figure 4: DAMIC-M 90% C.L. upper limits (solid thick black) on DM-electron
interactions through an ultralight (left) and heavy (right) dark photon
mediator obtained from the daily modulation analysis. Also shown are previous
limits from DAMIC-M [16] (solid black) and other experiments: DAMIC-SNOLAB
[13] (dashed black); SENSEI [14] (solid gray); EDELWEISS [15] (dashed gray);
SuperCDMS [12] (dotted gray); XENON1T combined result from [66, 67] (dashed
violet); PandaX-II [68] (dotted violet); a limit obtained from XENON10 data in
Ref. [69] (dash-dotted violet); and a limit obtained from XENON1T data
considering “solar reflected DM” (dashed orange) from Ref. [70] (left) and
Ref. [71] (right). Theoretical expectations assuming a DM relic abundance from
freeze-in and freeze-out mechanisms are also shown in light blue [72].

The DM signal $S$ is computed with Eq. (2) using QEDark [11] for the crystal
form factor $f_{\mathrm{c}}$, which encodes the electronic response of the
target. To convert the energy $E_{e}$ deposited in the DM interaction to
electron-hole pairs, we use a semi-empirical model based on probabilities
$P_{\mathrm{pair}}(n_{q}|E_{e})$ from the charge yield model in Ref. [73]. A
Monte Carlo simulation of the charge diffusion process is then applied to
include the detector response, with $\sigma_{xy}^{2}$ obtained from a surface
laboratory calibration with cosmic rays [74, 16]. Lastly, a $10\times 10$
binning of the simulated pixel array is applied to match the data-taking
conditions.

We minimize the joint log-likelihood for the U and L sides at fixed DM masses
between 0.53 and 2.7 MeV/c2. To illustrate the sensitivity of the method, we
show in Fig. 3 the residual rate after background subtraction, obtained by
fitting the U-side data to the background model only. Also shown is the
expectation for a DM particle of mass 1 MeV/c2 interacting via an ultralight
mediator with $\bar{\sigma}_{e}=10^{-32}$ cm2 (close to the current best limit
on the cross section at this mass).

The fit finds no preference for a DM signal at any mass. Consequently, we
derive exclusion limits using the likelihood-ratio test statistic $t_{q}$,
with null hypothesis $H_{0}$ given by $\bar{\sigma}_{e}=0$. The corresponding
90% C.L. exclusion limits for ultralight (left) and heavy (right) mediators
are shown in Fig. 4. These limits fall within the expected 68% sensitivity
band as estimated by Monte Carlo simulations. Also shown are previous results
from DAMIC-M and other direct detection experiments. Notably, the daily
modulation analysis improves up to $\sim$2 orders of magnitude the previous
DAMIC-M limits [16] obtained with the same data set. Cross-checks of the
analysis are included in the Supplemental Material.

For DM masses $\leq 2.7$ MeV/c2 the signal is overwhelmingly comprised of
1$e^{-}$. At higher masses the 2$e^{-}$ signal becomes relevant in
establishing the exclusion limit. However, a daily modulation analysis of the
$R_{2}^{i}$ rate will not improve existing DAMIC-M limits at higher masses,
since they already constrain the cross section to a level where Earth
interactions are negligible. For this reason, we do not perform a full daily
modulation analysis of $R_{2}^{i}$.

Conclusions.— This DAMIC-M search for DM particles interacting with electrons
excludes unexplored regions of parameter space between 0.53 and 2.7 MeV/c2.
This is the first time that the daily modulation due to scattering of DM
particles in the Earth, before they reach the detector, constrains DM-electron
interactions. When combined with our previous limits, DAMIC-M provides the
current best constraints from searches for a non-relativistic flux of DM
particles incident on Earth, over the mass ranges $[0.53,1000]$ MeV/c2 and
$[0.53,15.1]$ MeV/c2 for ultralight and heavy mediator interactions,
respectively 444Complementary searches for semi-relativistic DM fluxes from
Solar reflection [70, 71], and cosmological [77, 78] and stellar evolution
constraints [79, 80] also place limits in this mass range.. In addition, a
model-independent search for modulations with period up to 48 h demonstrates
the stability of our detector and allows for further improvements in time-
dependent searches for a DM signal.

Acknowledgments.—We thank the LSM for their support in the installation and
operation of the detector underground. The DAMIC-M project has received
funding from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation programme Grant Agreement No. 788137, and
from NSF through Grant No. NSF PHY-1812654. The work at University of Chicago
and University of Washington was supported through Grant No. NSF PHY-2110585.
This work was supported by the Kavli Institute for Cosmological Physics at the
University of Chicago through an endowment from the Kavli Foundation. We thank
the College of Arts and Sciences at University of Washington for contributing
the first CCDs to the DAMIC-M project. IFCA was supported by project
PID2019-109829GB-I00 funded by MCIN/ AEI /10.13039/501100011033. BJK
acknowledges funding from the Ramón y Cajal Grant RYC2021-034757-I, financed
by MCIN/AEI/10.13039/501100011033 and by the European Union
“NextGenerationEU”/PRTR. The Centro Atómico Bariloche group is supported by
ANPCyT grant PICT-2018-03069. The University of Zürich was supported by the
Swiss National Science Foundation. The CCD development work at Lawrence
Berkeley National Laboratory MicroSystems Lab was supported in part by the
Director, Office of Science, of the U.S. Department of Energy under Contract
No. DE-AC02-05CH11231. We thank Teledyne DALSA Semiconductor for CCD
fabrication. We acknowledge Santander Supercomputing support group at the
University of Cantabria who provided access to the supercomputer Altamira at
the Institute of Physics of Cantabria (IFCA-CSIC), member of the Spanish
Supercomputing Network, for performing simulations.

## References

  * [1] G. Bertone, D. Hooper, and J. Silk,  Phys. Rept. 405, 279 (2005), hep-ph/0404175.
  * [2] J. M. Gaskins,  Contemp. Phys. 57, 496 (2016), 1604.00014.
  * [3] O. Buchmueller, C. Doglioni, and L. T. Wang,  Nature Phys. 13, 217 (2017), 1912.12739.
  * [4] J. Billard et al.,  Rept. Prog. Phys. 85, 056201 (2022), 2104.07634.
  * [5] C. Boehm, P. Fayet, and J. Silk,  Phys. Rev. D 69, 101302 (2004), hep-ph/0311143.
  * [6] D. Hooper and K. M. Zurek,  Phys. Rev. D 77, 087302 (2008), 0801.3686.
  * [7] M. Pospelov, A. Ritz, and M. B. Voloshin,  Phys. Lett. B 662, 53 (2008), 0711.4866.
  * [8] X. Chu, T. Hambye, and M. H. G. Tytgat,  JCAP 05, 034 (2012), 1112.0493.
  * [9] S. Knapen, T. Lin, and K. M. Zurek,  Phys. Rev. D 96, 115021 (2017), 1709.07882.
  * [10] R. Essig, J. Mardon, and T. Volansky,  Phys. Rev. D 85, 076007 (2012), 1108.5383.
  * [11] R. Essig et al.,  JHEP 05, 046 (2016), 1509.01598.
  * [12] SuperCDMS Collaboration, R. Agnese et al.,  Phys. Rev. Lett. 121, 051301 (2018), 1804.10697,  [Erratum: Phys.Rev.Lett. 122, 069901 (2019)].
  * [13] DAMIC Collaboration, A. Aguilar-Arevalo et al.,  Phys. Rev. Lett. 123, 181802 (2019), 1907.12628.
  * [14] SENSEI Collaboration, L. Barak et al.,  Phys. Rev. Lett. 125, 171802 (2020), 2004.11378.
  * [15] EDELWEISS Collaboration, Q. Arnaud et al.,  Phys. Rev. Lett. 125, 141301 (2020), 2003.01046.
  * [16] DAMIC-M Collaboration, I. Arnquist et al.,  Phys. Rev. Lett. 130, 171003 (2023), 2302.02372.
  * [17] P. Adari et al.,  SciPost Phys. Proc. 9, 001 (2022), 2202.05097.
  * [18] C. Kouvaris and N. G. Nielsen,  Phys. Rev. D 92, 075016 (2015), 1505.02615.
  * [19] P. Sikivie and S. Wick,  Phys. Rev. D 66, 023504 (2002), astro-ph/0203448.
  * [20] M. S. Alenazi and P. Gondolo,  Phys. Rev. D 74, 083518 (2006), astro-ph/0608390.
  * [21] N. Bozorgnia, G. B. Gelmini, and P. Gondolo,  Nuclear Instruments and Methods in Physics Research A 654, 162 (2011), 1011.6006.
  * [22] N. Bozorgnia, G. B. Gelmini, and P. Gondolo,  Phys. Rev. D84, 023516 (2011), 1101.2876.
  * [23] S. Knapen, T. Lin, M. Pyle, and K. M. Zurek,  Phys. Lett. B 785, 386 (2018), 1712.06598.
  * [24] S. Griffin, S. Knapen, T. Lin, and K. M. Zurek,  Phys. Rev. D 98, 115034 (2018), 1807.10291.
  * [25] C. Boyd et al.,  (2022), 2212.04505.
  * [26] J. I. Collar and F. T. Avignone,  Phys. Lett. B 275, 181 (1992).
  * [27] J. I. Collar and F. T. Avignone, III,  Phys. Rev. D 47, 5238 (1993).
  * [28] F. Hasenbalg et al.,  Phys. Rev. D 55, 7350 (1997), astro-ph/9702165.
  * [29] R. Bernabei et al.,  The European Physical Journal C 74, 2827 (2014).
  * [30] LUX Collaboration, D. S. Akerib et al.,  Phys. Rev. D 98, 062005 (2018).
  * [31] PandaX-II Collaboration, X. Cui et al.,  Phys. Rev. Lett. 128, 171801 (2022).
  * [32] PROSPECT Collaboration, M. Andriamirado et al.,  Phys. Rev. D 104, 012009 (2021).
  * [33] B. Holdom,  Physics Letters B 166, 196 (1986).
  * [34] P. Galison and A. Manohar,  Physics Letters B 136, 279 (1984).
  * [35] T. Emken, R. Essig, C. Kouvaris, and M. Sholapurkar,  JCAP 09, 070 (2019), 1905.06348.
  * [36] C. Kouvaris and I. M. Shoemaker,  Phys. Rev. D 90, 095011 (2014), 1405.1729.
  * [37] N. Ávalos et al.,  Journal of Physics: Conference Series 2156, 012074 (2021).
  * [38] I. Arnquist et al.,  The DAMIC-M Experiment: Status and First Results, 2022, arXiv:2210.12070.
  * [39] S. E. Holland,  Exper. Astron. 14, 83 (2002).
  * [40] S. E. Holland, D. E. Groom, N. P. Palaio, R. J. Stover, and M. Wei,  IEEE Transactions on Electron Devices 50, 225 (2003).
  * [41] S. E. Holland, W. F. Kolbe, and C. J. Bebek,  IEEE Transactions on Electron Devices 56, 2612 (2009).
  * [42] DAMIC Collaboration, A. Aguilar-Arevalo et al.,  Phys. Rev. D 94, 082006 (2016).
  * [43] J. Janesick et al.,  Proc. SPIE 1242, 223 (1990).
  * [44] C. E. Chandler et al.,  Proc. SPIE 1242, 238 (1990).
  * [45] J. Tiffenberg et al.,  Phys. Rev. Lett. 119, 131802 (2017).
  * [46] DAMIC-M Collaboration, D. Norcini et al.,  Phys. Rev. D 106, 092001 (2022).
  * [47] J. R. Janesick,  Scientific charge-coupled devices (SPIE Optical Engineering Press, 2001).
  * [48] G. Cowan, K. Cranmer, E. Gross, and O. Vitells,  Eur. Phys. J. C 71, 1554 (2011), 1007.1727,  [Erratum: Eur.Phys.J.C 73, 2501 (2013)].
  * [49] E. Gross and O. Vitells,  Eur. Phys. J. C 70, 525 (2010), 1005.1891.
  * [50] S. K. Lee, M. Lisanti, S. Mishra-Sharma, and B. R. Safdi,  Phys. Rev. D 92, 083517 (2015), 1508.07361.
  * [51] S. M. Griffin, K. Inzani, T. Trickle, Z. Zhang, and K. M. Zurek,  Phys. Rev. D 104, 095015 (2021), 2105.05253.
  * [52] A. M. Green,  J. Phys. G 44, 084001 (2017), 1703.10102.
  * [53] We use the Standard Halo Model [76] parameters with the Earth’s velocity fixed to its average value during the data-taking period, $v_{E}=263$ km/s.
  * [54] C. McCabe,  JCAP 02, 027 (2014), 1312.1355.
  * [55] T. Emken and C. Kouvaris,  JCAP 10, 031 (2017), 1706.02249.
  * [56] G. D. Starkman, A. Gould, R. Esmailzadeh, and S. K. Dimopoulos,  Phys. Rev. D 41, 3594 (1990).
  * [57] B. J. Kavanagh, R. Catena, and C. Kouvaris,  JCAP 01, 012 (2017), 1611.05453.
  * [58] M. S. Mahdawi and G. R. Farrar,  JCAP 12, 004 (2017), 1709.00430.
  * [59] D. Hooper and S. D. McDermott,  Phys. Rev. D 97, 115006 (2018), 1802.03025.
  * [60] M. S. Mahdawi and G. R. Farrar,  JCAP 10, 007 (2018), 1804.03073.
  * [61] B. J. Kavanagh,  Phys. Rev. D 97, 123013 (2018), 1712.04901.
  * [62] B. J. Kavanagh,  Verne [Code, v1.3],  https://github.com/bradkav/verne, DOI:10.5281/zenodo.7193430, 2017-2023.
  * [63] Note that Verne uses an alternative convention for the isodetection angle $\gamma=180^{\circ}-\Theta$.
  * [64] A. Lantero-Barreda and B. J. Kavanagh,  in preparation, 2023.
  * [65] Note that the isodetection angle corresponding to a given time in a sidereal day changes by several degrees during the data taking period.
  * [66] XENON Collaboration, E. Aprile et al.,  Phys. Rev. Lett. 123, 251801 (2019).
  * [67] XENON Collaboration, E. Aprile et al.,  Phys. Rev. D 106, 022001 (2022).
  * [68] PandaX-II Collaboration, C. Cheng et al.,  Phys. Rev. Lett. 126, 211803 (2021).
  * [69] R. Essig, T. Volansky, and T.-T. Yu,  Phys. Rev. D 96, 043017 (2017).
  * [70] H. An, H. Nie, M. Pospelov, J. Pradler, and A. Ritz,  Phys. Rev. D 104, 103026 (2021).
  * [71] T. Emken,  Phys. Rev. D 105, 063020 (2022), 2102.12483.
  * [72] M. Battaglieri et al.,  US Cosmic Visions: New Ideas in Dark Matter 2017: Community Report, arXiv:1707.04591.
  * [73] K. Ramanathan and N. Kurinsky,  Phys. Rev. D 102, 063026 (2020).
  * [74] DAMIC Collaboration, A. Aguilar-Arevalo et al.,  Phys. Rev. D 105, 062003 (2022).
  * [75] Complementary searches for semi-relativistic DM fluxes from Solar reflection [70, 71], and cosmological [77, 78] and stellar evolution constraints [79, 80] also place limits in this mass range.
  * [76] D. Baxter et al.,  The European Physical Journal C 81, 907 (2021).
  * [77] G. Krnjaic and S. D. McDermott,  Phys. Rev. D 101, 123022 (2020), 1908.00007.
  * [78] C. Giovanetti, M. Lisanti, H. Liu, and J. T. Ruderman,  Phys. Rev. Lett. 129, 021302 (2022), 2109.03246.
  * [79] S. Davidson, S. Hannestad, and G. Raffelt,  JHEP 05, 003 (2000), hep-ph/0001179.
  * [80] J. H. Chang, R. Essig, and S. D. McDermott,  JHEP 09, 051 (2018), 1803.00993.
  * [81] T. Emken and C. Kouvaris,  Phys. Rev. D 97, 115047 (2018), 1802.04764.
  * [82] T. Emken and C. Kouvaris,  Dark Matter Simulation Code for Underground Scatterings (DaMaSCUS) [Code, v1.1],  https://github.com/temken/damascus, DOI:10.5281/zenodo.3726878, 2017-2020.
  * [83] S. Knapen, J. Kozaczuk, and T. Lin,  Phys. Rev. D 105, 015014 (2022).
  * [84] S. M. Griffin, K. Inzani, T. Trickle, Z. Zhang, and K. M. Zurek,  Phys. Rev. D 104, 095015 (2021).
  * [85] T. Trickle,  Phys. Rev. D 107, 035035 (2023).

## Supplemental Material

Independent cross-checks have been performed at every step in the analysis,
starting from the low-level image processing to the generation of the data
pixel distribution, the identification of defects, the modeling of the DM
signal, and the extraction of the DM signal upper limit. All of these checks
yield consistent results.

The observed time dependence of the rate $R_{1}^{i}$ is consistent with dark
current stabilizing over time, characteristic in these types of devices. To
confirm this interpretation we check the time dependence of the 2$e^{-}$ rate,
$R_{2}^{i}$, in the same data set. We sum the pixel charge distribution of 100
consecutive images to guarantee a sufficient number of pixels with 2$e^{-}$
and then fit it with the sum of 3 Gaussians to obtain $R_{2}^{i}$. The
exponential time decay constant for $R_{2}^{i}$ when compared with that for
$R_{1}^{i}$ is consistent with a Poisson process, as expected for dark
current.

In Fig. 5 we compare the Verne DM velocity distributions at the detector with
those obtained when the Earth-scattering effects are calculated using DaMaSCUS
[55, 81, 82], a more accurate but computationally demanding 3D Monte Carlo
simulation. The simplified approach adopted in Verne reproduces to better than
20% the DaMaSCUS predictions. We have repeated our analysis at a DM mass of 1
MeV/c2 using the DaMaSCUS code; the corresponding limits differ from those
obtained with Verne by less than 30%. The use of Verne has also allowed for a
number of additional checks. In the Verne calculation of the DM velocity
distribution the Earth’s speed is fixed to $263$ km/s. To check the effect of
this assumption we repeat the analysis using the largest and smallest value of
$\mathrm{v}_{\mathrm{E}}$, $267$ and $258$ km/s, corresponding to the
beginning and the end of the data-taking period. The exclusion limits change
by at most 15%.

We also perform a cross-check of the exclusion limits with a less sensitive
but simpler method by computing the correlation

$\rho=\frac{\sqrt{2}}{N_{\rm{im}}||R_{1}||}\sum_{i=1}^{N_{\rm{im}}}R_{1}^{i}\cos\left[\frac{2\pi(t_{i}-\phi)}{T}\right]\,,$
(4)

where $||R_{1}||=\sqrt{\sum_{i=1}^{N_{\rm{im}}}(R_{1}^{i})^{2}/N_{\rm{im}}}$
and $T$ and $\phi$ are fixed to the period and phase predicted by the
DM-$e^{-}$ scattering model described in the main text. The probability
distribution of $\rho$ is obtained from toy Monte Carlo simulations, which
include only the background model, and then used to derive a 90% limit on the
amplitude of a putative daily modulation signal. We convert the rate to cross
section using Eq. (2). The corresponding exclusion limits are a factor of 1.2
weaker than the nominal limits obtained from the likelihood.

For the model-independent analysis, we have performed toy-MC simulations under
$H_{0}$ (no modulation), finding that the test statistic $t_{q}$
asymptotically follows a $\chi^{2}$ distribution with two degrees of freedom.
We have also verified that our limits in Fig. 4 have the correct coverage by
performing the analysis on toy Monte Carlo simulation, which include both a
signal and background. We find that the signal is properly rejected or
discovered according to our 90% C.L.

We use QEDark as the reference theoretical model for proper comparison with
previous and forthcoming results from other experiments. We note that
theoretical uncertainties are significant for DM masses below $2.7$ MeV/c2,
with QEDark-based limits almost a factor 100 better than those obtained with
DarkELF [83] and EXCEED-DM [84, 85] models (see Supplemental Material of Ref.
[16].)

Figure 5:  Speed distributions at different isodetection angles for DM
particles of mass 1 MeV/c2 interacting via an ultralight dark photon mediator
(Verne, solid line; DaMaSCUS, dots). The distributions are calculated for
$\bar{\sigma}_{e}=10^{-33}$ cm2, close to the DAMIC-M limit at $m_{\chi}$=1
MeV/c2 reported in this letter.


# Humans prefer interacting with slow, less realistic butterfly simulations

Paige L. Reiter 1 & Talia Y. Moore1,2 1 Mechanical Engineering, Robotics
Institute, University of Michigan, Ann Arbor, MI<EMAIL_ADDRESS>Robotics, Ecology and Evolutionary Biology, Museum of Zoology, University of
Michigan, Ann Arbor, MI<EMAIL_ADDRESS>ORCID: 0000-0003-0867-4512

###### Abstract

How should zoomorphic, or bio-inspired, robots indicate to humans that
interactions will be safe and fun? Here, a survey is used to measure how human
willingness to interact with a simulated butterfly robot is affected by
different flight patterns. Flapping frequency, flap to glide ratio, and
flapping pattern were independently varied based on a literature review of
butterfly and moth flight. Human willingness to interact with these
simulations and demographic information were self-reported via an online
survey. Low flapping frequency and greater proportion of gliding were
preferred, and prior experience with butterflies strongly predicted greater
interaction willingness. The preferred flight parameters correspond to
migrating butterfly flight patterns that are rarely directly observed by
humans and do not correspond to the species that inspired the wing shape of
the robot model. The most realistic butterfly simulations were among the least
preferred. An analysis of animated butterflies in popular media revealed a
convergence on slower, less realistic flight parameters. This iterative and
interactive artistic process provides a model for determining human
preferences and identifying functional requirements of robots for human
interaction. Thus, the robotic design process can be streamlined by leveraging
animated models and surveys prior to construction.

###### Index Terms:

Human-Robot Interaction, Human-Animal Interaction, Zoomorphic Robot, Bio-
Inspired Design, Lepidoptera

##  I Introduction

Human interaction with animals and green spaces can have psychological and
emotional benefits [1]. For example, human perception of butterflies in
outdoor environments is strongly associated with restorative and relaxing
experiences [2, 3, 4]. Because the frequency and outcomes of these situations
are subject to chance, humans have developed lures — birdfeeders, salt licks,
nest boxes, etc. — to increase the likelihood of a positive interaction.
However, lures are less effective in urbanized areas, where wild plants and
animals are more rare [5].

Figure 1: Comparison of flight characteristics for real butterflies (purple
and images), animations of butterflies in popular media (orange), and
butterfly robots simulated for this study (green). The letters indicate the
popular media source, corresponding to rows in Tab. III. Note that in a) the
character analyzed is a moth; all other examples are butterflies. The numbers
indicate the video simulation, corresponding to trajectories in Fig. 3. The
preferred video of survey participants is indicated with a black arrow. Note
that for simulated dataset, each specified combination of wingbeat frequency
and proportion of gliding includes one aperiodic and one periodic datapoint.
Images of the biological butterflies obtained from Yale Peabody Museum
(_Danaus plexippus_), Muséum d’Histoire Naturelle de Toulouse (_Pieris
rapae_), and London (_Manduca sexta_). The shape of each plotted data point
indicates periodicity of the wingbeats.

Zoomorphic, or bio-inspired, robots can be used to increase the likelihood of
an interaction that is known to have positive emotional associations, such as
a butterfly landing on a hand. However, the positive emotional and
psychological benefits of interacting with a butterfly-mimicking robot would
be nullified if the combination of the motion and the appearance evoke an
uncanny or repulsive effect. Humans are repulsed by some robotic depictions of
living organisms, an effect which is often referred to as the ‘Uncanny Valley’
[6]. Such effects can arise from appearance and motion having different levels
of abstraction, stylization, or anthropomorphism from their biological
inspiration. For humanoid animations, naturalistic movement tends to enhance
acceptability of characters that are unacceptable when static, but distorted
movement decreases the acceptability of realistic characters [7]. Animal-
mimicking, or zoomorphic robots and animations can also evoke an uncanny
effect [8], especially when anthropomorphism is combined with more realistic
depictions [7]. The extent to which previous research on the Uncanny Valley
can be used to predict responses to invertebrate-mimicking robots is unknown
because invertebrate animals have less natural resemblance to human
appearances and are less familiar to humans.

Butterflies are frequently depicted in television shows, movies, and video
games to enhance the naturalistic ambiance of an environment [9, 10, 11, 4,
12]. Because the success of artistic depictions depends on human emotional
responses, the depictions of butterflies in popular media are likely to
indicate human preference. However, the extent to which popular animations are
informed by biological data varies greatly [13], and there are no guidelines
for assessing the likeability of zoomorphic animations based on the biological
accuracy of their movement.

Figure 2:  A) Female _Pieris rapae_ specimen (Photo by Sarefo (Wikimedia
Commons) licensed under CC BY-SA 3.0). B) Simple butterfly model in the
Blender environment, with the armature used for animations outlined in blue.

With limited information regarding human preferences with respect to
invertebrate motion, a simulated and interactive approach can direct the
design of robots for positive interactions prior to fabrication. Just as
simulations can aid in optimizing terrestrial robotic locomotion for stability
[14], human preference can be substituted as the goal [15].

In this paper, aspects of real butterfly flight behaviors are reviewed (Sec.
II) to inform flight simulations of a butterfly-mimicking robot (Sec. III).
The simulations are then used in a survey (Sec. IV) to determine the effect of
flight parameters on human willingness to interact with such a butterfly-
mimicking robot. The results from the survey (Sec. V) are compared to
depictions of butterflies in popular media (Sec. VI). The results of this
survey demonstrate that precise biomimicry is not always the most effective
strategy for designing zoomorphic robots, and further conclusions are drawn in
Section VII. Instead, by analyzing user feedback to simulations, human
preferences for specific motions, and the corresponding functional
capabilities of a robot, can be identified prior to fabrication.

##  II Biological Template

Butterflies and moths are flying insects of the order Lepidoptera, which
contains over $175,000$ species. Within Lepidoptera, wing shape and size, body
size, and flight characteristics all vary greatly. In fact, if we consider a
theoretical morphospace in which wing shape, wing size with respect to body,
and body size represent independent axes, lepidopteran insects have evolved
representatives across the entire available morphospace. Because wing shape
determines biomechanical capabilities, this morphological diversity implies
that there are many different types of flight characteristics that would be
considered realistic. Thus, the design of butterfly-inspired robots should
likely be informed by the intended function (in this case, human preference)
and the limitations of the robotic system, rather than a narrow range of
realistic movements.

Butterflies are some of the most popular insects in the world, likely due to
their bright coloration, role in pollination, and lack of stinging or biting
features [16, 17]. On the other hand, moths are often considered undesirable
pests [16]. Non-experts are easily able to distinguish between these groups
and quickly decide whether to attempt an interaction [18].

Butterflies and moths may be distinguished by examining morphology and flight
characteristics (Fig. 1). Butterflies generally have large forewings and
hindwings with respect to their small body size. For example, the large
forewings of monarch butterflies (_Danaus plexippus_) enable their long-range,
low-energy $3000$ mile gliding migrations between Canada and Mexico. During
their cross-continental journey, monarchs spend over $80\%$ of their time
gliding [19]. Most humans generally do not observe the high-altitude migratory
flight of monarch butterflies. Instead, we are much more likely to encounter
butterflies like the cabbage white (_Pieris rapae_), which exists in over 32
countries on every continent except South America and Antarctica [20]. Cabbage
white butterflies have average wingbeat frequencies between $10$ Hz and $15$
Hz and spend about $20\%$ of their time gliding [21]. On the other end of the
spectrum, hawk-moth (_Manduca sexta_) wings have a very high aspect ratio and
are relatively small with respect to their large body size. Hawk-moths are
capable of incredibly fast, high-energy flapping flight, sometimes reaching
speeds over $5.3$ m/s. Moths hardly ever glide, instead hovering in place at
low speeds for precise tasks, such as drinking nectar from flowers.
Furthermore, migratory monarchs can flap with wingbeat frequencies as low as
$3$ Hz to conserve energy, while many moths flap at frequencies over $24$ Hz
as they travel with rapid horizontal speeds [19], [22]. These flapping
frequencies are within the range of frequencies that can be distinguished by
the human eye [23, 24]. Data regarding the wingbeat frequency, approximate
flap to glide ratio, and wingbeat pattern were gathered from previous
biomechanical analyses and video examination. Because wing shape may be
difficult to distinguish during flapping flight, the proportion of gliding and
flapping and the wingbeat frequency are simple metrics that are likely used to
distinguish butterflies from moths.

##  III Simulated Model

Figure 3: A screenshot from a butterfly simulation that was presented to
participants, with the trajectories of all butterfly trials overlaid.

###  III-A Appearance of butterfly model

The basic shape of a butterfly was modeled in Blender [25] based on the female
_Pieris rapae_ (Fig. 2A), which was selected due to its ubiquity across
several continents [20]. Splines were used to trace the basic outline of the
wings, and these splines were extruded to make thin 3D meshes. Basic spheres
were stretched into the shapes of the head, abdomen, and thorax to match the
proportions of the butterfly template. Because anthropomorphism has been
demonstrated to produce diverging responses within a participant group [26],
human-like stylization was not explored in this study.

These shapes were attached to a simple armature to generate the animations
(Fig. 2 B). Because the motion of the wings and the overall flight trajectory
were the focus of this particular study, the head, abdomen, and thorax were
held fixed throughout the animated flights. Simple keyframes were used to
create a looping 1 second animation of a 3Hz flapping motion, a 5Hz flapping
motion, and a wings-extended gliding motion. These simple animations and the
model of simple shapes were exported from Blender and imported into Unity.

###  III-B Simulation of butterfly flight

The model and animation described in Section III-A were imported into Unity
[27]. C# scripts were written to control the motion of the game objects in
Unity. The butterfly is rendered in purple and shown against a simple blue and
green background representing the sky and ground, respectively (Fig. 3). The
simple visuals were used to provide an abstracted version of a naturalistic
background so as not to distract from the simulated butterfly.

###  III-C Simulated parameters

Based on the biomechanical features described in Section II, the following
parameters were independently varied in this study: wingbeat frequency,
proportion of gliding, and the pattern of wingbeats. For each parameter, two
values within the natural ranges of butterfly movement were selected to
present during the videos for comparison (Tab. I). The low option for wingbeat
frequency was set at $3$ Hz, and the high option was $5$ Hz, which is the
highest frequency renderable before the video outputs fail to generate smooth
flapping motion. The difference between the frequencies is noticeable by
humans [23, 24].

Random motions are often added to repetitive movements to make animations
appear more realistic and believable [28] and likely would have an effect on
human willingness to interact with the robot. Thus, periodic patterns were
compared with aperiodic patterns of flapping and gliding. Periodic simulations
consisted of repeating patterns of either $4:1$ or $1:1$ ratios of flapping to
gliding durations. For aperiodic flapping, although the flapping to gliding
ratios of the entire videos were either $4:1$ or $1:1$, the duration of each
bout of flapping was determined by randomly obtaining $50$ samples, with
replacement, from the set of numbers from $0$ to $9$, inclusive. Because only
a small number of samples were required to produce these simulations,
resampling the numbers for each video would result in trajectories that likely
differed greatly. Therefore, the same set of sampled numbers was used for all
trials with aperiodic wingbeat patterns. To alter flap to glide ratio of
trials with aperiodic wingbeat patterns, the sampled numbers were assigned to
different categories. To produce a $4:1$ flap to glide ratio, the gliding
motion was associated with values of $0$ to $1$ and the flapping motion was
associated with values of $2$ to $9$. For the $1:1$ flap to glide ratio, the
gliding motion was associated with values of $0$ to $4$, while the flapping
motion was associated with values of $5$ to $9$. The resulting trajectories
are shown in Figure 3.

Parameters | Frequency | Pattern | Flap to Glide Ratio
---|---|---|---
Video 1 |  $3$ Hz | Periodic | $4:1$
Video 2 |  $3$ Hz | Periodic | $1:1$
Video 3 |  $3$ Hz | Aperiodic | $4:1$
Video 4 |  $3$ Hz | Aperiodic | $1:1$
Video 5 |  $5$ Hz | Periodic | $4:1$
Video 6 |  $5$ Hz | Periodic | $1:1$
Video 7 |  $5$ Hz | Aperiodic | $4:1$
Video 8 |  $5$ Hz | Aperiodic | $1:1$
TABLE I: The flight parameters specified for each animated video.

While the aforementioned factors were varied, other parameters that might vary
in nature were held constant in all videos. The horizontal speed was held
constant at $3.89$ body lengths per second. To prevent distracting violations
of physics, a vertical speed was specified for each wing flap ($1.55$ body
lengths per second) and glide period ($-1.30$ body lengths per second). The
butterfly moved within a plane perpendicular to the viewer, thereby
maintaining a constant distance from the viewer. The butterfly did not land at
any point during any video, and it was consistently illuminated with diffuse
lighting as it traveled across the scene. No audio tracks were included in any
videos. Videos are posted in a playlist at https://bit.ly/33ZSYn3.

##  IV Experimental Methods

A survey was created in Qualtrics [29] to determine whether biomechanical
parameters influence the perceived friendliness of butterflies. The University
of Michigan Institutional Review Board determined that the study was exempt
from review (HUM$00184754$). The survey was primarily distributed via one
email to the faculty, staff, and students associated with the University of
Michigan’s Robotics Institute and the UM Ecology and Evolutionary Biology
department. A few responses came from secondary contacts to members within
those groups. Though the survey was anonymous, participants had the option of
listing their affiliation, age, familiarity with butterflies, and global
regional locations where they had spent at least five years of their life. The
core portion of the survey involved participants watching eight video
animations and responding to the following question for each: “Would you reach
your hand out to let this butterfly land on you?” Responses were reported on a
scale from $0$ to $7$, in which an answer of $0$ corresponds to “No thank
you,” while an answer of $7$ corresponds to “I hope it lands on me!” These
values are reported hereafter as Interaction Willingness (IW) scores. The
order in which the videos were presented for each participant was randomized,
but all videos were available at the same time, so participants could view the
videos as many times and in whatever order they chose. Lastly, participants
had the option to leave a short answer suggestion for what they would change
to make the butterfly appear friendlier. There was no time limit to complete
the survey. The survey was open for seven days and received $131$ responses.
The survey is posted online at https://bit.ly/2H2n68L.

##  V Experimental Results and Discussion

###  V-A Human preferences

Overall, the eight butterfly simulations were perceived as friendly, receiving
an average IW rating of $5.07$ overall on a scale from $0$ to $7$. The overall
favorite butterfly motion was featured in Video 2, which corresponded to $3$
Hz frequency wingbeat, a prescribed wingbeat pattern, and a $1:1$ ratio
between time spent gliding and time spent flying (Tab. II). In fact, the three
most highly ranked videos share a wingbeat frequency of $3$ Hz. The difference
in IW between the highest ranked video (2) and the lowest ranked video (5) was
highly significant (p = $1.67e^{-4}$)

| Mean(IW) | SD(IW) | Median(IW)
---|---|---|---
Video 1 | $5.16$ | $1.85$ | $6$
Video 2 | $5.36$ | $1.74$ | $6$
Video 3 | $5.23$ | $1.77$ | $5$
Video 4 | $5.05$ | $1.88$ | $5$
Video 5 | $4.78$ | $1.88$ | $5$
Video 6 | $5.05$ | $1.84$ | $5$
Video 7 | $4.85$ | $1.90$ | $5$
Video 8 | $5.06$ | $1.86$ | $5$
TABLE II: A summary of Interaction Willingness (IW) scores for each video

While all of the parameters shown in the videos were possible for real
butterflies, participants demonstrated a strong preference for the very slow
wingbeat frequency of $3$ Hz (p = $4.97e^{-5}$). Participants also
demonstrated a weak preference for more gliding (p = $0.0730$). Unless a
butterfly is actively descending, riding an air current like they might during
a long migratory event, or taking advantage of ground effects, butterflies are
very unlikely to have an even distribution of flapping and gliding time.
Unexpectedly, wingbeat pattern did not have a significant effect on IW scores
(p = $0.489$). Overall, the videos with the highest IW scores display flight
parameters corresponding to migratory flights above the tree canopy, rather
than flying through gardens at approximately human height.

These surprising results prompted further analysis to investigate whether
demographic factors informed these preferences.

###  V-B Demographic Factors

Although butterflies are globally distributed and generally considered
positively by many cultures, there is a chance that participants may be more
familiar with static representations of butterflies as decorative motifs if
they have not had personal encounters. There may also be differences among
individuals who generally like or dislike all butterflies, versus those who
are more discerning among butterflies. Thus, participant responses to several
optional questions regarding demographic information (e.g. age, experience
with butterflies, academic affiliation) were analyzed to determine their
effect on IW scores. A participant’s level of discerning among butterfly
simulations was evaluated by measuring the range of their IW scores, while
overall affinity was evaluated by measuring the mean.

Figure 4: Survey participant willingness to interact varied with respect to
different flight parameters. Figure 5: Comparison of Interaction Willingness
(IW) with respect to demographic factors. Asterisks indicate significance
level $p<0.05$. A period indicates significance level $0.05<p<0.1$. Academic
affiliation and previous experience with butterflies significantly predicted
mean IW scores for each participant. The range of IW scores for a participant
decreased with age, but this trend did not explain the majority of the
variation in the data. The y-axes were truncated at the maximum available
score of seven.

Participants who interact with butterflies as part of their job (n = $9$)
demonstrated highest overall affinity for butterfly videos and the widest
range of discerning. Participants who have not meaningfully interacted with
butterflies (n = $25$) scored the videos significantly lower than those who
have interacted with butterflies for fun (n = $97$, p = $0.0302$), and those
who interact with butterflies as a part of their job (p = $0.0450$). The mean
level of discerning was not significantly determined by affiliation. Although
this option was available, no participants reported zero familiarity with
butterflies.

Participants affiliated with the Ecology and Evolutionary Biology department
(n = $53$) showed the highest overall and most consistent IW scores regardless
of previous experience with butterflies. Biologists rated videos significantly
higher than participants affiliated with the Robotics Institute (n = $51$, p =
$5.22e^{-3}$) and somewhat higher than those with Other affiliations (n =
$27$, p = $0.0759$). All participants who interact with butterflies as a part
of their jobs were affiliated with the Biology department, but biologists who
do not study butterflies were the least discerning participants. Roboticists
who had interacted with butterflies for fun ranked videos higher on average
than those who have had insignificant interaction with butterflies (p =
$0.0256$), but all roboticists had similar levels of discerning. Participants
with Other affiliations who have interacted with butterflies for fun were
almost as discerning as those who study butterflies as their occupation, but
differences in discerning according to experience were not significant within
this group.

Participant ages ranged from 7 to 71111Children were not solicited for
participation in the survey, but multiple survey entries per household were
accepted., with $50$% of participants within the range of 24 to 34. Age of
participant did not significantly affect the overall affinity for butterflies
(p = $0.2$). Discerning among butterfly videos significantly decreased with
age (p = $0.0218$), but this trend did not explain the variation in most of
the data (Adjusted $R^{2}$ = $0.0337$).

There were no interactions among the biomechanical (fixed effects) and
demographic (random effects) when tested as a linear mixed-effects model using
the lmer function in the lmerTest package for R [30].

###  V-C Open responses

During the survey period, multiple participants directly contacted the
researchers with positive messages, such as “… thinking about robot
butterflies made my day a little brighter.” An analysis of the verbatim
responses at the end of the survey brought to light some interesting topics.
$66$% of the $131$ participants chose to leave a note at the end of the survey
in the box asking for suggestions on how to make the butterflies appear even
more friendly. $26$ participants explicitly suggested alterations to the
static visual appearance of the butterfly, asking researchers to make the
robot look like a familiar monarch species or to add anthropomorphic eyes. The
butterfly in this study was intentionally designed to be nondescript so that
positive associations with a specific species due to previous interactions
would be equally unlikely for all participants, regardless of regional origin.
The effects of stylization and anthropomorphization can be examined in future
studies.

$54$ participants made suggestions related to behavior, which was broken down
further into several categories. $15$ participants explicitly described the
flapping motion, including requests to increase or decrease the overall speed
or adding variability in the wingbeat frequency. $15$ participants described
the gliding behavior, with some calling for more gliding and others calling
for less. A few participants questioned the physics of such a rapid descent
associated with gliding. Lastly, $10$ participants addressed the horizontal
speed of the butterfly, which was a fixed variable across all trials. Most
participants suggested slowing butterflies to speeds that correspond to values
below the range of realistic speeds.

$30$ suggestions described the overall flight trajectory. These suggestions
ranged from choosing a more simple trajectory to suggestions of more erratic
trajectories. $16$ participants requested more interaction with the butterfly.
Several suggestions included alterations in trajectories and in interaction,
and were counted in both. $16$ participants described wanting to know whether
the butterfly wanted to interact with them, suggesting either that the flight
path come near or circle the observer, or that the flapping frequency decrease
with proximity to the observer. Although such indications of intent are
potentially informed by interactions with humans, dogs, or cats instead of
butterflies, the underlying concept of the robot giving the appearance of
informed consent is a fascinating avenue for future research.

##  VI Depictions of butterflies in popular media

###  VI-A Data collection

| Popular Media Name | Type | Year | Studio | Director/Creator
---|---|---|---|---|---
a | _A Bug’s Life_ | Feature-length film | $1998$ | Pixar Animation Studios | John Lasseter
b | _The Road to El Dorado_ | Feature-length film | $2000$ | DreamWorks Animation | Eric Bergeron and Don Paul
c | _Big Hero 6_ | Feature-length film | $2014$ | Walt Disney Animation Studios | Don Hall and Chris Williams
d | _Miraculous: Tales of Ladybug and Cat Noir_ | Television show | $2015$ | Zagtoon and Method Animation | Thomas Astruc, et al.
e | _Lovecraft Country_ | Television show | $2020$ | Monkeypaw Productions, et al. | Cheryl Dunye
f | _Second Life_ | Video game | $2003$ | Linden Lab | Hana Kozlowski
g | _Legend of Zelda: Breath of the Wild_ | Video game | $2017$ | Nintedo | Hidemaro Fujibayashi
TABLE III: The popular media sources of the animated butterflies analyzed for
this study.

The authors recorded animated butterflies encountered during their normal
interaction with video games, animated television shows, live-action
television shows, and animated feature-length films during a four-month
period, from June 11 to October 11, 2020 (see Tab. III). Wingbeat frequencies
and patterns were determined by counting the number of video frames per wing
beat, and the proportion of gliding was measured over the entire duration of
the recorded flight.

###  VI-B Trends in animated butterfly motion

Butterflies depicted in popular media had wingbeat frequencies between 3 and 6
Hz, much like the animations created for the survey (Fig. 1). This is on the
extremely low end for real-life butterflies, similar to the frequencies used
by migrating butterflies that are rarely observed by casual viewers. The
depictions of butterflies in feature-length films and live-action television
shows featured periods of gliding between the flaps, though the examples from
video games and animated television shows did not include any gliding. In the
characterized films, the butterflies spent about $20\%$ to $25\%$ of their
time gliding. For the other media types, the lack of gliding may be a
consequence of simplifying the animation process by limiting butterfly flight
to a fully periodic pattern so that the same animations can be reused in
multiple episodes of a television show or in different locations in a video
game. Feature-length films and live-action television shows were also more
likely to feature aperiodic flapping cycles.

##  VII Conclusion

Based on the overall positive rankings of simulations, the results of this
study suggest that butterfly-mimicking robots can facilitate positive
interactions with humans. Although the results were generally positive, they
were not uniform, indicating that humans do have preferences for slower forms
of butterfly motions. Therefore, butterflies are an appropriate inspiration
for zoomorphic robots, and tuning the flight behaviors is likely to have an
effect on user acceptance and interaction. This preference may be explained by
an aversion to moths, which display a higher wingbeat frequency and greater
proportion of flapping to gliding. Alternatively, a preference for slower and
more predictable motion may be driving this pattern. Such ‘sensory bias’ for a
specific feature has been revealed to drive the evolution of extreme features
in animals undergoing intersexual selection [31]. Repeating this experiment
using more abstract images, such as floating spheres, can help determine
whether moth aversion or sensory bias drives this preference, and how these
preferences compare to predictions based on the Uncanny Valley.

Our results provide foundational information for examining human affinity for
representations of invertebrate animals as a function of their similarity to a
real organism. Unlike previous examinations of static vertebrate animal
representations that revealed Uncanny Valley patterns of likeability with
respect to verisimilitude [7, 8], humans in this study preferred motions that
are slower than observed in real butterflies. To determine whether the
preference patterns found in this study can be extended to all invertebrates,
a range that spans human likeability, from spiders and centipedes to
butterflies and ladybugs, could be examined.

Anthropomorphism may also influence the Uncanny Valley effect with respect to
vertebrate animal representations [7]. Because invertebrate animals have less
morphological resemblance to humans, the effect anthropomorphism has on human
preference involves more suspension of disbelief and is therefore difficult to
predict. Furthermore, examining the response of children to picture books with
or without anthropomorphism demonstrated diverging effects in boys and girls
[26]. Thus, future studies should carefully consider whether and when
anthropomorphism should be applied to invertebrate-mimicking robots.

This study also adds to a growing body of work examining the effect of
movement on the Uncanny Valley. In humanoid representations, realistic
movement has been shown to rehabilitate statically unfamiliar representations,
but distorted movement decreases the familiarity of statically realistic
representations [32, 33]. Here, although static and dynamic representations
were not compared, variation in movement patterns associated with different
proportions of gliding had only a weak effect on human willingness to interact
with the simulations. Humans are likely attuned to sensing changes in
trajectory randomness, as evidenced by targeting accuracy decreasing for more
erratically moving targets [34, 35]. Animals may also experience strong
selective pressures favoring trajectory randomness for predator evasion [36,
37, 38]. Although the trajectories simulated in this study varied as a result
of altering the flapping and gliding patterns, the trajectory randomness was
not directly varied. In future research, studying the effect of random
movements by varying trajectory entropy can precisely characterize how humans
respond to erratic motion of zoomorphic robots.

The findings from this study converge with the representations of butterflies
in animated media. Interestingly, a survey of popular media indirectly
indicates that slower and more predictable butterflies are aesthetically
pleasing, even if they share less resemblance to real animal movement (Fig.
1). This phenomenon could potentially be explained by the fact that artistic
animations often go through rounds of feedback and revision to shape their
final product. Here, we demonstrate how a similar interactive and iterative
process can be used to design robotic systems for user interface. Leveraging
simulated motion is a useful strategy to decrease the effort required for each
iteration, thereby streamlining the manufacturing process by determining
functional requirements prior to robot construction [15].

The feedback from the participants, especially the open responses, suggest how
future interactive and simulated studies can provide additional guidance for
the design and deployment of such a robotic system. An important takeaway from
this study is that demographic factors strongly determined human willingness
to interact with a butterfly robot. While the generic visual appearance of the
butterfly in this study was designed to ensure consistent levels of
familiarity across cultural backgrounds, examination of regional preferences
or interaction with the local setting can inform the design of distinct models
for different target audiences. As animals are also frequently represented in
animation and robotics, examining the effect of their diverse movement
patterns on human acceptance will likely benefit many different research
fields and industries.

##  VIII Acknowledgements

The authors thank E.L. Moore for providing footage and description of the
animated butterfly in _Second Life_ , A. Reiter for assistance with
simulations, members of the EMBiR lab for providing feedback on a pilot
version of the survey, and members of the ROAHM Lab and J. Socha for providing
feedback on the manuscript.

## References

  * [1] Thea Cameron-Faulkner, Joanna Melville and Merideth Gattis  “Responding to nature: Natural environments improve parent-child communication”  In _Journal of Environmental Psychology_ 59.February Elsevier Ltd, 2018, pp. 9–15  DOI: 10.1016/j.jenvp.2018.08.008
  * [2] Richard A. Fuller et al.  “Psychological benefits of greenspace increase with biodiversity”  In _Biology Letters_ 3.4, 2007, pp. 390–394  DOI: 10.1098/rsbl.2007.0149
  * [3] Erminia Pedretti and Barbara J. Soren  “Reconnecting to the natural world through an immersive environment”  In _Canadian Journal of Science, Mathematics and Technology Education_ 6.1, 2006, pp. 83–96  DOI: 10.1080/14926150609556689
  * [4] Masashi Soga et al.  “Urban residents’ perceptions of neighbourhood nature: Does the extinction of experience matter?”  In _Biological Conservation_ 203 Elsevier Ltd, 2016, pp. 143–150  DOI: 10.1016/j.biocon.2016.09.020
  * [5] Robert B. Blair and Alan E. Launer  “Butterfly diversity and human land use: Species assemblages along an urban gradient”  In _Biological Conservation_ 3207.96, 1997, pp. 113–125
  * [6] Masahiro Mori  “Bukimi no tani [The Uncanny Valley]”  In _Energy_ 7, 1970, pp. 33–35
  * [7] V. Schwind et al.  “Is there an uncanny valley of virtual animals? A quantitative and qualitative investigation”  In _International Journal of Human Computer Studies_ 111.November 2016 Elsevier Ltd, 2018, pp. 49–61  DOI: 10.1016/j.ijhcs.2017.11.003
  * [8] Diana Loffler, Judith Dorrenbacher and Marc Hassenzahl  “The uncanny valley effect in zoomorphic robots: The U-shaped relation between animal likeness and likeability”  In _ACM/IEEE International Conference on Human-Robot Interaction_ , 2020  DOI: 10.1145/3319502.3374788
  * [9] Richard J. Leskosky and May R. Berenbaum  “Insects in Animated Films: Not All ”Bugs” Are Bunnies”  In _Bulletin of the Entomological Society of America_ 34.2, 1988, pp. 55–63
  * [10] Penesta Dika  “The Research of Motifs in Interactive Media Art Concerning the Visual Aspect”, 2008
  * [11] Meike Uhrig  “Emotion in Animated Films”  New York, NY: Routledge, 2018
  * [12] James N. Hogue  “Cultural Entomology”  In _Annual Review of Entomology_ 32, 1987, pp. 181–99  DOI: 10.1016/B978-0-12-374144-8.00072-2
  * [13] Bruce Blumberg  “Biomimetic Animated Creatures”  In _Biologically Inspired Intelligent Robots_ SPIE Press, 2003, pp. 47–73
  * [14] C. Paul and J.C. Bongard  “The road less travelled: morphology in the optimization of biped robot locomotion”  In _Proceedings 2001 IEEE/RSJ International Conference on Intelligent Robots and Systems. Expanding the Societal Role of Robotics in the the Next Millennium (Cat. No.01CH37180)_ 1, 2001, pp. 226–232  DOI: 10.1109/IROS.2001.973363
  * [15] James M. Bern, Kai Hung Chang and Stelian Coros  “Interactive design of animated plushies”  In _ACM Transactions on Graphics_ 36.4, 2017  DOI: 10.1145/3072959.3073700
  * [16] Stephen R. Kellert  “Values and Perceptions of Invertebrates”  In _Conservation Biology_ 7.4, 1993, pp. 845–855  DOI: 10.1046/j.1523-1739.1993.740845.x
  * [17] Raynald H. Lemelin et al.  “Humans, Insects and Their Interaction: A Multi-faceted Analysis”  In _Animal Studies Journal_ 5.1, 2016, pp. 65–79
  * [18] Jürg Schlegel, Gabriele Breuer and Reto Rupf  “Local insects as flagship species to promote nature conservation? A survey among primary school children on their attitudes toward invertebrates”  In _Anthrozoos_ 28.2, 2015, pp. 229–245  DOI: 10.1080/08927936.2015.11435399
  * [19] Alejandro Ortega Ancel et al.  “Aerodynamic evaluation of wing shape and wing orientation in four butterfly species using numerical simulations and a low-speed wind tunnel, and its implications for the design of flying micro-robots”  In _Interface Focus_ , 2017  DOI: 10.1098/rsfs.2016.0087
  * [20] Sean F. Ryan et al.  “Global invasion history of the agricultural pest butterfly _Pieris rapae_ revealed with genomics and citizen science”  In _Proceedings of the National Academy of Sciences of the United States of America_ 116.40, 2019, pp. 20015–20024  DOI: 10.1073/pnas.1907492116
  * [21] Ngoc San Ha and Nam Seo Goo  “Flapping frequency and resonant frequency of insect wings”  In _2013 10th International Conference on Ubiquitous Robots and Ambient Intelligence, URAI 2013_ , 2013  DOI: 10.1109/URAI.2013.6677463
  * [22] Alexander P. Willmott and Charles P. Ellington  “The mechanics of flight in the hawkmoth Manduca sexta I. Kinematics of hovering and forward flight”  In _Journal of Experimental Biology_ , 1997
  * [23] J. Heck  “The Flicker Electroretinogram of the Human Eye”  In _Acta Physiologica Scandinavica_ 39.2-3, 1957, pp. 158–166  DOI: 10.1111/j.1748-1716.1957.tb01417.x
  * [24] Thomas L. Harrington and Marcia K. Harrington  “Perception of motion using blur pattern information in the moderate and high-velocity domains of vision”  In _Acta Psychologica_ 48.1-3, 1981, pp. 227–237  DOI: 10.1016/0001-6918(81)90064-0
  * [25] Blender Foundation  “Blender”, https://blender.org/, 1994–2020
  * [26] Megan Geerdts  “Anthropomorphic media and children’s biological knowledge”, 2014  DOI: 10.7282/T32Z13RV
  * [27] Unity Technologies  “Unity, Version 2019.3.8f1”, https://unity.com/, 2019
  * [28] Sooha Park Lee, Jeremy B. Badler and Norman I. Badler  “Eyes alive”  In _Proceedings of the 29th Annual Conference on Computer Graphics and Interactive Techniques, SIGGRAPH ’02_ , 2002, pp. 637–644  DOI: 10.1145/566570.566629
  * [29] “Qualtrics”, https://qualtrics.com/, 2002–2020
  * [30] Alexandra Kuznetsova, Per B. Brockhoff and Rune H. B. Christensen  “lmerTest Package: Tests in Linear Mixed Effects Models”  In _Journal of Statistical Software_ 82.13, 2017, pp. 1–26  DOI: doi:10.18637/jss.v082.i13
  * [31] Rebecca C. Fuller, David Houle and Joseph Travis  “Sensory bias as an explanation for the evolution of mate preferences”  In _American Naturalist_ 166.4, 2005, pp. 437–446  DOI: 10.1086/444443
  * [32] Álvaro Castro-González, Henny Admoni and Brian Scassellati  “Effects of form and motion on judgments of social robots’ animacy, likability, trustworthiness and unpleasantness”  In _International Journal of Human Computer Studies_ 90 Elsevier, 2016, pp. 27–38  DOI: 10.1016/j.ijhcs.2016.02.004
  * [33] Lukasz Piwek, Lawrie S. McKay and Frank E. Pollick  “Empirical evaluation of the uncanny valley hypothesis fails to confirm the predicted effect of motion”  In _Cognition_ 130.3 Elsevier B.V., 2014, pp. 271–277  DOI: 10.1016/j.cognition.2013.11.001
  * [34] Graham Richardson et al.  “Unpredictable movement as an anti-predator strategy”  In _Proceedings of the Royal Society B: Biological Sciences_ 218.20181112., 2018
  * [35] Christofer J. Clemente and Robbie S. Wilson  “Speed and maneuverability jointly determine escape success: Exploring the functional bases of escape performance using simulated games”  In _Behavioral Ecology_ 27.1, 2016, pp. 45–54  DOI: 10.1093/beheco/arv080
  * [36] Talia Yuki Moore, Kimberly L. Cooper, Andrew A. Biewener and Ramanarayan Vasudevan  “Unpredictability of escape trajectory explains predator evasion ability and microhabitat preference of desert rodents”  In _Nature Communications_ 8.440, 2017, pp. 1–9  DOI: 10.1038/s41467-017-00373-2
  * [37] Stacey A Combes, D. E. Rundle, J. M. Iwasaki and James D. Crall  “Linking biomechanics and ecology through predator-prey interactions: flight performance of dragonflies and their prey”  In _The Journal of Experimental Biology_ 215, 2012, pp. 903–913  DOI: 10.1242/jeb.059394
  * [38] Marta A. Skowron Volponi, Donald James McLean, Paolo Volponi and Robert Dudley  “Moving like a model: Mimicry of hymenopteran flight trajectories by clearwing moths of Southeast Asian rainforests”  In _Biology Letters_ 14.5, 2018  DOI: 10.1098/rsbl.2018.0152