# Shear and interface modes in neutron stars with pasta structures

Hajime Sotani<EMAIL_ADDRESS>Astrophysical Big Bang Laboratory,
RIKEN, Saitama 351-0198, Japan  Interdisciplinary Theoretical & Mathematical
Science Program (iTHEMS), RIKEN, Saitama 351-0198, Japan

###### Abstract

We carefully examine the shear and interface modes, which are excited due to
the presence of crust elasticity, in neutron stars with pasta structures,
adopting the relativistic Cowling approximation. We find that the shear modes
are independent of the presence of the cylindrical-hole and spherical-hole
nuclei at least up to a few kilohertz, while the interface modes strongly
depend on the presence of the cylindrical-hole and spherical-hole nuclei. In
addition, we find empirical relations for the interface mode frequencies
multiplied by the stellar mass and for the shear mode frequencies multiplied
by the stellar radius. These relations are expressed as a function of the
stellar compactness almost independently of the stiffness in a higher-density
region inside the neutron star, once one selects the crust equation of state.
Thus, if one would simultaneously observe the shear and interface modes from a
neutron star, one might extract the neutron star mass and radius with the help
of the constraint on the crust stiffness obtained from terrestrial
experiments.

###### pacs:

04.40.Dg, 97.10.Sj, 04.30.-w, 26.60.Gj

##  I Introduction

Neutron stars are one of the most suitable natural laboratories for probing
physics under extreme states. This object is provided as a massive remnant
remaining after the supernova explosion, which happens at the last moment of
the star’s life [1]. It is considered that the density inside the neutron
stars easily exceeds the standard nuclear density, while the gravitational and
magnetic fields inside/around the neutron stars become much stronger than
those observed in our Solar system. Thus, one could inversely see the aspect
of such extreme conditions by carefully observing the neutron stars and their
phenomena.

One of the most important observational pieces of evidence is the discovery of
the $2M_{\odot}$ neutron stars [2, 3, 4, 5]. It is theoretically known that
the neutron star has a maximum mass, which depends on the equation of state
(EOS) for neutron star matter. Owing to the discovery of the $2M_{\odot}$
neutron stars, now one can exclude the EOSs, with which the expected maximum
mass does not reach the observations. Meanwhile, through careful observations
of the pulsar light curve, one may extract the neutron star properties,
especially the stellar compactness. This is because the light emitted from the
neutron star’s surface can bend due to the strong gravitational field induced
by the neutron star, which is one of the relativistic effects (e.g., Refs. [6,
7, 8, 9, 10, 11]). In fact, the observations with the Neutron Star Interior
Composition ExploreR (NICER), which is a NASA telescope on the International
Space Station, give us the constraints on the neutron star mass and radius,
i.e., PSR J0030+0451 [12, 13] and PSR J0740+6620 [14, 15]. Moreover, the
observations of gravitational waves from the binary neutron star merger,
GW170817 [16], enables us to constrain the tidal deformability of the neutron
stars, which leads to the constraint on the $1.4M_{\odot}$ neutron star radius
[17]. In addition to these astronomical observations, the neutron star
properties in a lower-density region can be gradually constrained through
terrestrial nuclear experiments. It may be reasonable that one discusses both
constraints simultaneously in the neutron star mass and radius plane [18, 19,
20].

The oscillation frequencies of the neutron stars must be another important
piece of information for extracting the neutron star properties. Since the
oscillation frequencies strongly depend on the interior properties of the
object, one may be able to know the interior properties as an inverse problem
by observing the oscillation frequencies. This technique is known as
asteroseismology, which is similar to seismology on Earth and helioseismology
on Sun. In practice, by identifying the quasi-periodic oscillations observed
in the afterglow following the magnetar giant flares [21, 22] with the crustal
torsional oscillations, one can constrain the crust properties (e.g., Refs.
[23, 24, 25, 26, 27, 28]). This is a good example of how asteroseismology
works well. Furthermore, using the gravitational waves from neutron stars, one
could get information on neutron star mass, radius, and EOS (e.g., Refs. [29,
30, 31, 32, 33, 34, 35, 36, 37, 38, 39]).

In the observation of the neutron star oscillations, several eigenfrequencies
may be simultaneously excited, where each eigenfrequency corresponds to a
different physical process. By identifying an observed frequency with a
specific eigenfrequency one by one, one can see the physics corresponding to
the eigenmode. For example, if one would observe the fundamental oscillations
of a neutron star, one may extract the average density of the neutron star
[29, 30]. In a similar way, if one would detect the stellar oscillations
associated with the crust elasticity, such as torsional, shear, or interface
modes, one may glimpse crust properties. Up to now, there are a few studies
about the non-radial oscillations of neutron stars with an elastic crust,
e.g., Refs. [40, 41, 42, 36]. However, in these previous studies, the shear
and interface modes, which are polar-type oscillations exited due to the crust
elasticity, have been discussed only on a specific stellar model, where the
systematical study has never been done. So, in this study, we will
systematically examine the shear and interface modes, adopting the realistic
neutron star model with the crust elasticity, where not only the phase
composed of spherical nuclei but also the phase composed of non-spherical
nuclei, the so-called pasta phase, are considered as crust equilibrium models.
As the density increases, the shape of the nuclei changes from spherical (SP)
to cylindrical (C), slab-like (S), cylindrical-hole (CH), and spherical-hole
(SH) before the matter becomes uniform (U).

This manuscript is organized as follows. In Sec. II, we mention the
equilibrium models together with the EOS adopted in this study and also the
shear modulus in the pasta phase. In Sec. III, we briefly mention the
perturbation equations governing the non-radial oscillations in neutron stars
with elastic crust. Then, we discuss the eigenfrequencies in Sec. IV and
systematically examine their dependence on the neutron star properties in Sec.
V. Finally, we conclude this study in Sec. VI. Unless otherwise mentioned, we
adopt geometric units in the following, $c=G=1$, where $c$ and $G$ denote the
speed of light and the gravitational constant, respectively.

##  II EOS and Equilibrium models

In this study, we simply consider a non-rotating, strain-free, and spherically
symmetric neutron star as an equilibrium model. The metric describing such an
object is given by

$ds^{2}=-e^{2\Phi}dt^{2}+e^{2\Lambda}dr^{2}+r^{2}\left(d\theta^{2}+\sin^{2}\theta
d\phi^{2}\right),$ (1)

where $\Phi$ and $\Lambda$ are the metric functions depending on only $r$. In
particular, $\Lambda$ is directly connected to the enclosed gravitational
mass, $m$, inside the radial position, $r$, through $e^{-2\Lambda}=1-2m/r$.
The stellar models are constructed by integrating the Tolman-Oppenheimer-
Volkoff equation together with an appropriate EOS for neutron star matter.
Since the Fermi temperature of a neutron star is generally much higher than
the real temperature of a neutron star, one can neglect the thermal effect on
the neutron star structure. But, one may have to carefully handle the surface
envelope, where the density becomes low enough for appearing the thermal
effect [43]. Nevertheless, in this study, we simply consider the neutron star
models, assuming that the surface density is $10^{6}$ g/cm3 and the density of
the surface of the outer crust is $10^{10}$ g/cm3. In addition, in this study,
we will consider the shear oscillations, which can be excited due to the non-
zero crustal elasticity. Since the shear oscillations are a kind of polar-type
oscillations, we have to consider not only the crustal region but also the
core region, even though the shear oscillations are confined only inside the
crust region. This is a completely different situation from the study of the
torsional oscillations (axial-type oscillations), which can be discussed only
inside the neutron star crust, e.g., Refs. [25, 26], separately from the
neutron star core.

The bulk energy per nucleon for the zero temperature uniform nuclear matter
expected for any EOS models is expressed as a function of the baryon number
density, $n_{\rm b}$, and an asymmetry parameter, $\alpha$:

$\frac{E}{A}=w_{s}(n_{\rm b})+\alpha^{2}S(n_{\rm b})+{\cal O}(\alpha^{3}),$
(2)

where $n_{b}$ and $\alpha$ are given by $n_{\rm b}=n_{n}+n_{p}$ and
$\alpha=(n_{n}-n_{p})/n_{\rm b}$ with the neutron number density, $n_{n}$, and
the proton number density, $n_{p}$. In this expression, $w_{s}$ corresponds to
the energy per nucleon of symmetric nuclear matter ($\alpha=0$), while $S$
denotes the density-dependent symmetry energy. In addition, $w_{s}$ and $S$
can be expanded in the vicinity of the saturation density, $n_{0}$, for the
symmetric nuclear matter as a function of $u=(n_{\rm b}-n_{0})/(3n_{0})$:

$\displaystyle w_{s}(n_{\rm b})=w_{0}+\frac{K_{0}}{2}u^{2}+{\cal O}(u^{3}),$
(3) $\displaystyle S(n_{\rm b})=S_{0}+Lu+{\cal O}(u^{2}).$ (4)

The coefficients in these expansions are the nuclear saturation parameters,
which characterize each EOS. That is, each EOS has its own set of nuclear
saturation parameters. In this study, to examine the dependence on the
symmetry energy, we adopt the phenomenological EOSs constructed by Oyamatsu
and Iida [44, 45] (hereafter referred to as OI-EOSs). The OI-EOSs are
specially constructed for the given values of $K_{0}$ and $L$ so that the
other saturation parameters are tuned in such a way as to reproduce the
empirical nuclear data for stable nuclei, adopting a simplified version of the
extended Thomas-Fermi theory [44, 45]. This is because $n_{0}$, $w_{0}$, and
$S_{0}$ are well constrained from the terrestrial experiments, while the other
parameters, i.e., $K_{0}$ and $L$ (and the additional saturation parameters
associated with the higher terms), are relatively more difficult to be
constrained experimentally. Nevertheless, the constraint on $K_{0}$ is
gradually becoming more severe, i.e., $K_{0}=240\pm 20$ MeV [46], while $L$ is
still less constrained [47, 18]. So, in this study, we especially adopt the
EOS models with $K_{0}=230$ MeV. In Table 1, we list the EOS parameters
adopted in this study, where we also show the transition density from a
specific pasta phase to the next pasta phase. We note that the crust thickness
strongly depends on $L$ and stellar compactness, $M/R$, [45, 48], i.e., the
thickness decreases as $L$ and $M/R$ increase, and the pasta structures almost
disappear, using the EOS model with
$L\,\raisebox{-3.44444pt}{$\stackrel{{\scriptstyle\textstyle>}}{{\sim}}$ }100$
MeV [45].

Table 1:  The EOS parameters adopted in this study. SP-C, C-S, S-CH, CH-SH, and SH-U denote the transition densities for the OI-EOSs characterized by $K_{0}$ and $L$. In addition, for the $1.4M_{\odot}$ neutron star model constructed with each EOS, the ratio of the thickness of the elastic region composed of spherical and cylindrical nuclei, $\Delta R_{\rm SpCy}$, and that composed of cylindrical-hole and spherical-hole nuclei, $\Delta R_{\rm CHSH}$, to the stellar radius, $R$, is also listed.  $K_{0}$ (MeV) |  $L$ (MeV) | SP-C (fm-3) | C-S (fm-3) | S-CH (fm-3) | CH-SH (fm-3) | SH-U (fm-3) | $\Delta R_{\rm SpCy}/R$ | $\Delta R_{\rm CHSH}/R$
---|---|---|---|---|---|---|---|---
230 | 42.6 | 0.06238 | 0.07671 | 0.08411 | 0.08604 | 0.08637 | 0.0007731 | 0.05406
230 | 73.4 | 0.06421 | 0.07099 | 0.07284 | 0.07344 | 0.07345 | 0.0003773 | 0.06740

The shear modulus, $\mu$, is an additional integrant to discuss the shear
oscillations. The shear modulus, $\mu_{\rm sp}$, in the body-centered cubic
(bcc) lattice composed of the spherical nuclei has been formulated as a
function of the ion number density, $n_{i}$, the charge number of the ion,
$Z$, and a Wigner-Seitz cell radius, $a$, i.e., $4\pi a^{3}/3=1/n_{i}$ [49]:

$\mu_{\rm sp}=0.1194\frac{n_{i}(Ze)^{2}}{a}.$ (5)

This expression of the shear modulus should be modified a little due to the
phonon contribution [50], the electron screening effect [51], the
polycrystalline effect [52], and the effect of finite-sizes of atomic nuclei
[53], but in this study, we simply adopt the standard expression given as Eq.
(5). On the other hand, the shear modulus, $\mu_{\rm cy}$, in the phase
composed of the cylindrical nuclei is expressed as a function of the Coulomb
energy per volume of a Wigner-Seitz cell, $E_{\rm Coul}$, and the volume
fraction of cylindrical nuclei, $w_{2}$, as

$\mu_{\rm cy}=\frac{2}{3}E_{\rm Coul}\times 10^{2.1(w_{2}-0.3)},$ (6)

and the shear modulus, $\mu_{\rm sl}$, in the phase composed of the slab-like
nuclei can be considered as

$\mu_{\rm sl}=0$ (7)

against the linear perturbations, i.e., the deformation energy due to the
distortion becomes of higher order contribution [54]. We note that the elastic
properties of phases with nonspherical nuclei, i.e., pasta phase, in a neutron
star have been also discussed in Refs. [55, 56], which suggested the
possibility that the tiny but non-zero elastic constant may appear in the
polycrystalline lasagna (slab-like nuclei). By taking into account this
feature, the results shown in this study may be changed. Additionally, the
shear modulus, $\mu_{\rm ch}$ ($\mu_{\rm sh}$), in the phase composed of the
cylindrical-hole (spherical-hole) nuclei can be derived in the same way as
$\mu_{\rm cy}$ ($\mu_{\rm sp}$) because the liquid crystalline structure of
cylindrical-hole (spherical-hole) nuclei is the same as that of cylindrical
(spherical) nuclei (see Ref. [27] for details). In Figs. 1 and 2, as an
example, we show the radial profile of energy density and shear modulus for
the stellar model with $1.4M_{\odot}$ and $12.4$ km constructed using the OI-
EOS with $K_{0}=230$ and $L=73.4$ MeV.

Finally, we should mention the effect of superfluidity, although we simply
neglect such an effect in this study. In general, the enthalpy density
effectively decreases because of the fact that the superfluid neutron does not
contribute to the oscillations, which leads to an increase in the frequencies
(at least for the torsional oscillations) [25, 57]. We will see how the
superfluidity can change the frequencies of the shear and interface modes in
the future.

Figure 1: Radial profile of energy density for the neutron star model with
$1.4M_{\odot}$ and 12.4 km, using the EOS with $L=73.4$ MeV. In the top panel,
the vertical dashed lines from left to right denote the boundary between the
core and crust, the boundary between the inner and outer crust, and the
boundary between the crust and envelope.  Figure 2: The effective shear
modulus for various phases in neutron star curst. The bottom panel is an
enlarged view of the shaded region shown in the top panel, where Sp, Cy, CH,
and SH denote the phases composed of spherical, cylindrical, cylindrical-hole,
and spherical-hole nuclei, respectively. We note that the shear modulus in the
phase of slab-like nuclei becomes zero against the linear response [54]. For
reference, we show the boundary between the core and crust and the boundary
between the crust and envelope in the top panel, and the boundary between the
core and SH; the boundary between CH and the phase of slab-like nuclei; and
the boundary between the phase of slab-like nuclei and Cy in the bottom panel
with the dashed lines.

##  III Perturbation equations

In this study, we simply adopt the relativistic Cowling approximation, i.e.,
the metric is fixed during the fluid oscillations. Even with this
approximation, one can qualitatively discuss the behavior of eigenfrequencies
[58]. The Lagrangian displacement, $\xi^{i}$, for the polar-type oscillations
is generally given with the spherical harmonics, $Y_{\ell m}(\theta,\phi)$, by

$\xi^{i}=\left(rW,V\frac{\partial}{\partial\theta},V\frac{1}{\sin^{2}\theta}\frac{\partial}{\partial\phi}\right)Y_{\ell
m}(\theta,\phi)e^{i\sigma t},$ (8)

where $W$ and $V$ are the functions of $r$ and $\sigma$ is the eigenvalue,
while the pressure perturbation, $\delta p$, is expressed with the energy
density, $\epsilon$, and pressure, $p$, for the equilibrium models as

$\displaystyle\delta p=(\epsilon+p)H(r)Y_{\ell m}(\theta,\phi)e^{i\sigma t}.$
(9)

We note that the polar-type oscillations can be discussed completely apart
from the axial-type oscillations, because of the nature of the spherically
symmetric background. In this study, we also consider the adiabatic
oscillations, i.e.,

$\Delta p=\frac{p\Gamma}{\epsilon+p}\Delta\epsilon,$ (10)

where $\Gamma$ is the adiabatic index and $\Delta Q$ denotes the Lagrangian
perturbation of a quantity $Q$, which is associated with the Eulerian
perturbation, $\delta Q$, through

$\Delta Q=\delta Q+\xi^{r}\frac{dQ}{dr}.$ (11)

So, one can derive the relation between $\delta p$ and $\delta\epsilon$ as

$\delta p=c_{s}^{2}\delta\epsilon+p\Gamma\xi^{r}A_{r},$ (12)

where $c_{s}$ is the sound velocity and $A_{r}$ is the relativistic
Schwarzschild discriminant given by

$\displaystyle c_{s}^{2}\equiv\left(\frac{\partial
p}{\partial\epsilon}\right)_{s}=\frac{\Delta
p}{\Delta\epsilon}=\frac{p\Gamma}{\epsilon+p},$ (13) $\displaystyle
A_{r}=\frac{1}{\epsilon+p}\frac{d\epsilon}{dr}-\frac{1}{p\Gamma}\frac{dp}{dr}.$
(14)

On the other hand, the shear strain tensor, $\Sigma_{\mu\nu}$, is described
through the relation of

${\cal
L}_{u}\Sigma_{\mu\nu}=\frac{2}{3}\Sigma_{\mu\nu}\nabla_{\alpha}u^{\alpha}+\sigma_{\mu\nu},$
(15)

where ${\cal L}_{u}$ is the Lie derivative along the direction of fluid four-
velocity, $u^{\mu}$, and $\sigma_{\mu\nu}$ is the rate of shear tenor [59],
which are respectively expressed as

$\displaystyle{\cal
L}_{u}\Sigma_{\mu\nu}=u^{\alpha}\nabla_{\alpha}\Sigma_{\mu\nu}+\Sigma_{\alpha\nu}\nabla_{\mu}u^{\alpha}+\Sigma_{\mu\alpha}\nabla_{\nu}u^{\alpha},$
(16) $\displaystyle\sigma_{\mu\nu}=\frac{1}{2}\left(P^{\alpha}_{\
\nu}\nabla_{\alpha}u_{\mu}+P^{\alpha}_{\
\mu}\nabla_{\alpha}u_{\nu}\right)-\frac{1}{3}P_{\mu\nu}\nabla_{\alpha}u^{\alpha}.$
(17)

Here, $P_{\mu\nu}$ is the projection tensor given by

$P_{\mu\nu}=g_{\mu\nu}+u_{\mu}u_{\nu}.$ (18)

Using these quantities, the contribution from the shear strain in perturbation
of the energy-momentum tensor is given by

$\delta T_{\mu\nu}^{(s)}=-2\mu\delta\Sigma_{\mu\nu},$ (19)

assuming a Hookean relationship [60], where $\mu$ is the shear modulus
discussed in the previous section.

With the relativistic Cowling approximation, the perturbation equations are
derived from the linearized energy-momentum conservation laws, i.e.,
$\nabla_{\mu}\delta T^{\mu\nu}=0$. The concrete system of equation is shown in
Appendix A for the elastic region and in Appendix B for the fluid region
($\mu=0$), while the boundary and junction conditions, which should be
imposed, are shown in Appendix C [41]. In practice, the perturbation equations
are integrated outward from the center and inward from the stellar surface
with the appropriate boundary conditions, where the corresponding solutions
are named as $(y_{1}^{\rm in},y_{2}^{\rm in})$ and $(y_{1}^{\rm
out},y_{2}^{\rm out})$, respectively. Then, the eigenfrequencies are
determined via the condition of

$\Delta\equiv y_{1}^{\rm in}y_{2}^{\rm out}-y_{1}^{\rm out}y_{2}^{\rm in}=0$
(20)

at some position inside the star, e.g., the boundary at the crust and
envelope. In this study, we especially focus on the $\ell=2$ modes.

##  IV Eigenfrequencies

In order to understand the dependence of eigenfrequencies excited in the
neutron stars on the presence of elasticity, first we consider (i) the neutron
star composed of fully zero-elastic “fluid”; (ii) the stellar model with
elastic phase composed of spherical nuclei, “Sp”; (iii) the stellar model with
elastic phase composed of spherical and cylindrical nuclei, “Sp+Cy”; and (iv)
the “realistic” stellar model with elastic phase composed of spherical,
cylindrical, cylindrical-hole, and spherical-hole nuclei as shown in Fig. 2.
That is, we focus on a specific neutron star model with $1.4M_{\odot}$ and
$12.4$ km constructed with the EOS with $L=73.4$ MeV, but the shear moduli in
some elastic phases are artificially put to zero except for the realistic
stellar model. In Fig. 3, one can see how the eigenfrequencies depend on the
elastic phases, where the absolute value of $\Delta$ at the boundary between
the crust and envelope is shown as a function of the frequency, i.e., the
eigenfrequencies correspond to the frequency where abs$(\Delta)=0$. The
resultant eigenfrequencies excited in the stellar models shown in Fig. 3 are
listed in Table 2.

Figure 3: The absolute values of $\Delta$ given by Eq. (20) are shown as a function of the frequencies. The eigenfrequencies correspond to the specific frequencies, with which the absolute value of $\Delta$ at some position inside the star becomes zero. In the left pane, we show the results for the stellar model composed of only the fluid with the solid line and that including a non-zero elastic region composed of spherical nuclei with the dotted line. In the middle panel, we show the results for the stellar model including a non-zero elastic region composed of spherical nuclei with the dotted line and that including a non-zero elastic region composed of spherical and cylindrical nuclei with the solid line. The right panel is just an enlarged view of the middle panel, where we also show the result for the “realistic” stellar model with the dashed line. The neutron star model is the same as in Fig. 2. Table 2:  Eigenfrequencies excited in the stellar models shown in Fig. 3 in the unit of kHz.  | fluid | Sp | SpCy | SpCy+CHSH
---|---|---|---|---
$f$ | 2.237 | 2.237 | 2.237 | 2.237
$p_{1}$ | 6.075 | 6.074 | 6.074 | 6.074
$i_{1}$ | — | 0.039 | 0.036 | 0.036
$i_{2}$ | — | 0.030 | 0.030 | 0.030
$i_{3}$ | — | — | — | 0.027
$s_{1}$ | — | 0.730 | 0.668 | 0.668
$s_{2}$ | — | 1.201 | 1.111 | 1.111
$s_{3}$ | — | 1.552 | 1.477 | 1.477
$s_{4}$ | — | 1.877 | 1.752 | 1.752
$s_{5}$ | — | 2.316 | 2.139 | 2.139

In the left panel of Fig. 3, we show the results for the “fluid” model with
the solid line and for the “Sp” model with the dotted line. From this figure,
one can observe the excitation of the shear ($s_{i}$-) and interface
($i_{i}$-) modes together with the fundamental ($f$-) and pressure ($p_{i}$-)
modes in the “Sp” model (see also the right panel), where $f$\- and
$p_{1}$-mode frequencies excited in the “Sp” model are almost the same as
those excited in the “fluid” model. That is, the presence of the elasticity
hardly affects the acoustic oscillations. In the middle panel, we show the
results for the “Sp” model with the dotted line and for the “Sp+Cy” model with
the solid line. From this result, one can observe that the $s_{i}$-mode
frequencies in the “Sp+Cy” model become smaller than those in the “Sp” model.
In the right panel, we show an enlarged view of the middle panel, where we
also show the results for the “realistic” model (SpCy+CHSH). From the right
panel, one can observe that the $i_{2}$-mode frequency in the “Sp+Cy” model is
the same as that in the “Sp” model, while the $i_{1}$-mode frequency strongly
depends on the presence of the phase composed of cylindrical nuclei. In
addition, for the “realistic” model one can observe an additional mode, i.e.,
the $i_{3}$-mode, together with the $i_{1}$\- and $i_{2}$-modes excited in the
“Sp+Cy” model. In particular, except for the excitation of the $i_{3}$-mode,
we find that the eigenfrequencies excited in the “realistic” model are the
same as those in the “Sp+Cy” model at least in the frequency domain shown in
Fig. 3. This may be because the phase composed of cylindrical-hole and
spherical-hole nuclei is quite narrow and the effect can not appear in the
frequency domain considered here, as discussed below.

Figure 4: For a realistic stellar model with $1.4M_{\odot}$ and $12.4$ km
using the EOS with $L=73.4$ MeV, the eigenfunctions of $W$ (solid lines) and
$V$ (dashed lines) for the $i_{2}$-mode ($i_{1}$-mode) are shown in the top
(bottom) panel. The vertical dotted lines from right to left denote the
boundary between the crust and envelop (Sp/envelope), the boundary between the
phases composed of cylindrical and slab-like nuclei (Slab/Cy), the boundary
between the phases of slab-like and cylindrical-hole nuclei (CH/Slab), and the
boundary between the phase composed of spherical-hole nuclei and core
(core/SH).  Figure 5: For a realistic stellar model with $1.4M_{\odot}$ and
$12.4$ km using the EOS with $L=73.4$ MeV, the eigenfunctions of $W$ (solid
lines) and $V$ (dashed lines) for the $i_{3}$-mode are shown in the top panel,
while an enlarged view of the top panel is shown in the bottom panel. The
meaning of the vertical lines is the same as in Fig. 4.

The $i_{i}$-modes are the eigenmodes excited due to the presence of the
interface between the phases with zero and non-zero elasticity [42]. In this
paper, we simply assign the $i_{i}$-modes in order from the highest to the
lowest frequencies. That is, we have only two $i$-modes in the “Sp+Cy” model
as shown in the right panel of Fig. 3, because there are two interfaces, i.e.,
the interface between the envelope and crust and the interface between the
phases composed of cylindrical and slab-like nuclei. In Fig. 4 we show the
amplitude of eigenfunctions, $W$ and $V$, for $i_{2}$-mode ($i_{1}$-mode) in
the top (bottom) panel excited in the “realistic” model, where the vertical
dotted lines denote the boundary between the envelope and crust (Sp/envelope);
the boundary between the phases of cylindrical and slab-like nuclei (Slab/Cy);
the boundary between the phases of slab-like and cylindrical-hole nuclei
(CH/Slab); and the boundary between the phase of spherical-hole nuclei and
core (core/SH). From this figure, one can see that the $i_{2}$-mode is
associated with the interface between the envelope and crust, while the
$i_{1}$-mode is with the interface between the phases composed of cylindrical
and slab-like nuclei. This may be a reason why the $i_{2}$-mode frequency in
the “Sp” model is the same as that in the “Sp+Cy” model, as shown in the right
panel of Fig. 3. Anyway, in the “realistic” model, one can see the tiny effect
of the non-zero elasticity in the phase composed of cylindrical-hole and
spherical-hole nuclei in the amplitude of eigenfunctions, even for the
$i_{1}$\- and $i_{2}$-modes (e.g., see the inset in the top panel of Fig. 4).
Moreover, in Fig. 5, we show the amplitude of the $i_{3}$-mode in the top
panel, while an enlarged view of the top panel is shown in the bottom panel.
One can observe that the amplitude of the $i_{3}$-mode becomes dominant inside
the region composed of the cylindrical-hole and spherical-hole nuclei. We also
find that only three interface modes are excited in the “realistic” model,
even though four interfaces exist in the “realistic” model, i.e., Sp/envelope,
Slab/Cy, CH/Slab, and core/SH. This may come from the fact that the region
composed of cylindrical-hole and spherical-hole nuclei is too narrow. In fact,
if the elastic region becomes too narrow, the number of excited interface
modes can become less than the number of interfaces, as shown in Appendix D.
Furthermore, if one considers the neutron star model using the EOS with
$L=42.6$ MeV, where the ratio of the thickness of the elastic region composed
of cylindrical-hole and spherical-hole nuclei to the stellar radius is
relatively larger than that considered in Fig. 3 as shown in Table 1, one can
observe the $i_{4}$-mode together with the $i_{3}$-mode by introducing the
elastic region composed of cylindrical-hole and spherical-hole nuclei (see
Sec. V for details).

Figure 6: For a realistic stellar model with $1.4M_{\odot}$ and $12.4$ km
using the EOS with $L=73.4$ MeV, the eigenfunctions of $W$ (solid lines) and
$V$ (dashed lines) for the $s_{1}$-mode ($s_{2}$-mode) are shown in the top
(bottom) panel. The meaning of the vertical lines is the same as in Fig. 4.
Figure 7: Eigenfunctions of $V$ for the $s_{i}$-modes with $i=1-5$.

The $s_{i}$-modes are also the eigenmodes excited due to the presence of
elasticity. Unlike the $i_{i}$-modes, the $s_{i}$-modes are basically confined
inside the elastic region. In Fig. 6, we show the amplitude of the
eigenfunction, $W$ and $V$, for the $s_{1}$-mode ($s_{2}$-mode) in the top
(bottom) panel. From this figure, one can see that $W$ is continuous even at
the boundaries between the anelastic and elastic regions owing to the junction
condition (see in Appendix C), while $V$ is discontinuous at the boundaries.
In addition, one can observe that the nodal number in the eigenfunction of the
$s_{i}$-modes is equivalent to the subscript $i$. This behavior is easily
observed in the amplitude of $V$, as shown in Fig. 7, where we show the
amplitude of V for the $s_{i}$-modes with $i=1-5$. So, the wavelength of the
$s_{i}$-mode, $\lambda_{i}$, is roughly estimated as

$\lambda_{i}\simeq 2\Delta R/i,$ (21)

where $\Delta R$ denotes the thickness of the elastic region, in which the
shear modes are confined. Then, the corresponding frequencies are also
estimated as

$f_{i}\approx v_{s}/\lambda_{i},$ (22)

where $v_{s}\approx(\mu/\epsilon)^{1/2}$ denotes the shear velocity [36]. With
this simple estimation, one may understand why the $s_{i}$-mode frequencies
excited in the “realistic” model are the same as those in the “Sp+Cy” model as
discussed in Fig. 3. That is, the $s_{i}$-mode frequencies excited in the
phase of cylindrical-hole and spherical-hole nuclei must be much higher (maybe
more than 100 times higher) than the $s_{i}$-mode frequencies excited in the
phase of spherical and cylindrical nuclei, because $\Delta R$ for the phase of
cylindrical-hole and spherical-hole nuclei, $\Delta R_{\rm CHSH}$, is much
thinner than $\Delta R$ for the phase of spherical and cylindrical nuclei,
$\Delta R_{\rm SpCy}$, as shown in Table 1, i.e., $\Delta R_{\rm SpCy}/\Delta
R_{\rm CHSH}=179$.

##  V Dependence on the neutron star properties

First, in Fig. 8, we show the eigenfrequencies of the $i$-, $s$-, and
$f$-modes as a function of the stellar compactness for the neutron star models
constructed using the EOSs with $L=42.6$ and $73.4$ MeV. From this figure, one
can observe that the $i$-mode frequencies weakly depend on the stellar
compactness, while the $s$-mode frequencies monotonically increase with the
stellar compactness. This is because the ratio of the thickness of the elastic
region to the stellar radius decreases as the stellar compactness increases
[48], which leads to the increase of the $s$-mode frequencies as discussed
with Eqs. (21) and (22). In addition, as mentioned before, since the thickness
of the elastic region composed of the cylindrical-hole and spherical-hole
nuclei for the neutron star model constructed with $L=42.6$ MeV is relatively
larger than that with $L=73.4$ MeV, one can observe the $i_{4}$-mode together
with $i_{1}$-, $i_{2}$-, and $i_{3}$-modes in the stellar model with $L=42.6$
MeV.

Figure 8: Eigenfrequencies of the $i_{i}$-modes (top panels) and the $s_{i}$\-
and $f$-modes (bottom panels) are shown as a function of the stellar
compactness, $M/R$, for the neutron star models constructed using the EOSs
with $L=42.6$ (left panels) and $73.4$ MeV (right panels).

Next, we examine how the frequencies of the $i$\- and $s$-modes depend on the
neutron star properties. In particular, since the nuclear properties in the
core region (or in a higher-density region) are quite uncertain, we examine
the frequencies of the $i$\- and $s$-modes by changing the stiffness of the
EOS in a higher-density region. For this purpose, in addition to the original
OI-EOSs listed in Table 1, we simply consider that the EOS for a lower density
region of $\varepsilon\leq\varepsilon_{t}$, i.e., OI-EOSs, is connected to the
one-parameter EOS characterized by $\alpha$ for a higher density region of
$\varepsilon\geq\varepsilon_{t}$, i.e.,

$p=\alpha(\varepsilon-\varepsilon_{t})+p_{t},$ (23)

where $p_{t}$ is given from the EOS for a lower-density region with
$\varepsilon=\varepsilon_{t}$ and $\alpha$ is associated with the sound
velocity, $c_{s}$, as $c_{s}^{2}=\alpha$ [61]. In this study, we especially
adopt that $\epsilon_{t}$ is equivalent to twice the saturation density,
focusing on the value of $\alpha$ in the range of $1/3\leq\alpha\leq 1$.

In practice, if one calculates the frequencies of $i$\- and $s$-mode with this
type of EOS, the frequencies depend on the value of $\alpha$. However, we find
that the $i$-mode frequencies multiplied by the stellar mass, $fM$, can be
expressed as a function of the stellar compactness almost independently of the
value of $\alpha$ (or the stiffness in a higher density region inside the
neutron stars), only depending on the stiffness of the curst EOS, as shown in
Fig. 9. In this figure, the solid lines are the fitting lines given by the
functional form as

$fM\ ({\rm kHz}/M_{\odot})=a_{0}+a_{1}(x/0.1)+a_{2}(x/0.1)^{2},$ (24)

where $x$ denotes the stellar compactness, $M/R$, and $a_{0}$, $a_{1}$, and
$a_{2}$ are adjusted coefficients.

Figure 9: The $i_{i}$-mode frequencies multiplied by the stellar mass are
shown as a function of the stellar compactness. The filled and open marks
denote the results for the neutron star model using the EOSs with $L=42.6$ and
$73.4$ MeV, respectively. The circles denote the results for the neutron star
models with the original OI-EOSs, while the diamonds, squares, and triangles
denote the results for the neutron star models with the OI-EOS connected to
the EOS characterized by $\alpha=1/3$, 0.6, and 1, respectively. The solid
lines denote the fitting lines using the formula given by Eq. (24).

In a similar way, we also find that the $s$-mode frequencies multiplied by the
stellar radius, $fR$, can be expressed as a function of the stellar
compactness almost independently of the value of $\alpha$, which depends only
on the crust stiffness, as shown in Fig. 10. In this figure, the solid lines
are fitting lines given by the functional form as

$fR\ ({\rm kHz\ km})=b_{0}+b_{1}(x/0.1),$ (25)

where $x$ denotes the stellar compactness, $M/R$, and $b_{0}$ and $b_{1}$ are
adjusted coefficients. Now, we find two different types of fitting formulae
for the $i$\- and $s$-mode frequencies. Thus, if one would simultaneously
observe the $i$\- and $s$-modes, one might extract the stellar mass and radius
with the help of the constraint on the crust stiffness from the terrestrial
experiments.

Figure 10: The $s$-mode frequencies multiplied by the stellar radius are shown
as a function of the stellar compactness. The solid lines denote the fitting
lines using the formula given by Eq. (25).

##  VI Conclusion

We carefully examine the frequencies of the interface and shear oscillations,
which are excited due to the presence of the elasticity, by considering the
neutron star models with the pasta structures, i.e., cylindrical, slab-like,
cylindrical-hole, and spherical-hole nuclei at the basis of the crust. We find
that the shear mode frequencies excited in a realistic stellar model are
basically the same as those in the neutron star model composed of only
spherical and cylindrical nuclei, if we focus only on the frequency range up
to a few kHz. This is because the shear modes are only excited inside the
elastic region, which leads to the feature that the frequencies are inversely
proportional to the thickness of the elastic region, and the thickness of the
elastic region composed of cylindrical-hole and spherical-hole nuclei is
extremely thin. On the other hand, the interface mode frequencies strongly
depend on the elastic region composed of cylindrical-hole and spherical-hole
nuclei. We find that the number of the interface mode frequencies depends on
the thickness of the elastic region composed of cylindrical-hole and
spherical-hole nuclei (or the value of the slop parameter $L$). In addition,
we find the empirical relations for the interface mode frequencies multiplied
by the stellar mass and for the shear mode frequencies multiplied by the
stellar radius as a function of the stellar compactness, which is almost
independent of the stiffness in a higher-density region inside the neutron
stars, once one selects the crust equation of state. Via our empirical
relations, if one would simultaneously observe the interface and shear mode
oscillations from a neutron star, one might extract the stellar mass and
radius with the help of the constraint on the crust stiffness obtained from
the terrestrial experiments.

###### Acknowledgements.

We are grateful to Shijun Yoshida, Christian J. Krüger, and Akira Harada for
their valuable comments. This work is supported in part by Japan Society for
the Promotion of Science (JSPS) KAKENHI Grant Numbers JP19KK0354 and
JP21H01088, and by Pioneering Program of RIKEN for Evolution of Matter in the
Universe (r-EMU).

##  Appendix A Perturbation equations inside the elastic region

The perturbation equations are derived from the linearized energy-momentum
conservation laws:

$\displaystyle r\frac{dz_{1}}{dr}=$
$\displaystyle-\left(1+\frac{2\alpha_{2}}{\alpha_{3}}+U_{2}\right)z_{1}+\frac{1}{\alpha_{3}}z_{2}+\frac{\alpha_{2}}{\alpha_{3}}\ell(\ell+1)z_{3},$
(26) $\displaystyle r\frac{dz_{2}}{dr}=$
$\displaystyle\left[\left(-3-U_{2}+U_{1}-e^{2\Lambda}c_{1}\bar{\sigma}^{2}\right)U_{3}+\frac{4\alpha_{1}}{\alpha_{3}}\left(3\alpha_{2}+2\alpha_{1}\right)\right]z_{1}$
$\displaystyle+\left(U_{4}-\frac{4\alpha_{1}}{\alpha_{3}}\right)z_{2}+\left[U_{3}-2\alpha_{1}\left(1+\frac{2\alpha_{2}}{\alpha_{3}}\right)\right]\ell(\ell+1)z_{3}$
$\displaystyle+e^{2\Lambda}\ell(\ell+1)z_{4},$ (27) $\displaystyle
r\frac{dz_{3}}{dr}=$
$\displaystyle-e^{2\Lambda}z_{1}+\frac{1}{\alpha_{1}}e^{2\Lambda}z_{4},$ (28)
$\displaystyle r\frac{dz_{4}}{dr}=$
$\displaystyle\left(U_{3}-6\Gamma\frac{\alpha_{1}}{\alpha_{3}}\right)z_{1}-\frac{\alpha_{2}}{\alpha_{3}}z_{2}$
$\displaystyle-\left[c_{1}\bar{\sigma}^{2}U_{3}+2\alpha_{1}-\frac{2\alpha_{1}}{\alpha_{3}}\left(\alpha_{2}+\alpha_{3}\right)\ell(\ell+1)\right]z_{3}$
$\displaystyle-\left(3+U_{2}-U_{4}\right)z_{4},$ (29)

where the variables, $z_{i}$ for $i=1-4$, are defined as

$\displaystyle z_{1}=$ $\displaystyle W,$ (30) $\displaystyle z_{2}=$
$\displaystyle
2\alpha_{1}e^{-\Lambda}\frac{d}{dr}\left(re^{\Lambda}W\right)+\left(\Gamma-\frac{2\alpha_{1}}{3}\right)\left[\frac{1}{r^{2}}e^{-\Lambda}\frac{d}{dr}\left(r^{3}e^{\Lambda}W\right)-\ell(\ell+1)V\right],$
(31) $\displaystyle z_{3}=$ $\displaystyle V,$ (32) $\displaystyle z_{4}=$
$\displaystyle\alpha_{1}\left(e^{-2\Lambda}r\frac{dV}{dr}+W\right).$ (33)

We note that the variables $z_{2}$ and $z_{4}$ are proportional to the radial
and transverse tractions [62]. The various quantities in the coefficients are
defined as

$\displaystyle\alpha_{1}=$ $\displaystyle\frac{\mu}{p},$ (34)
$\displaystyle\alpha_{2}=$ $\displaystyle\Gamma-\frac{2\alpha_{1}}{3},$ (35)
$\displaystyle\alpha_{3}=$ $\displaystyle\Gamma+\frac{4\alpha_{1}}{3},$ (36)
$\displaystyle U_{1}=$
$\displaystyle\left(\frac{d\Phi}{dr}\right)^{-1}\frac{d}{dr}\left(r\frac{d\Phi}{dr}\right),$
(37) $\displaystyle U_{2}=$ $\displaystyle r\frac{d\Lambda}{dr},$ (38)
$\displaystyle U_{3}=$
$\displaystyle\left(1+\frac{\epsilon}{p}\right)r\frac{d\Phi}{dr},$ (39)
$\displaystyle U_{4}=$ $\displaystyle\frac{r\epsilon}{p}\frac{d\Phi}{dr},$
(40) $\displaystyle c_{1}=$
$\displaystyle\frac{M}{R^{3}}re^{-2\Phi}\left(\frac{d\Phi}{dr}\right)^{-1},$
(41)

and $M$, $R$, and $\bar{\sigma}$ are the stellar mass, radius, and
$\bar{\sigma}\equiv\sigma(R^{3}/M)^{1/2}$.

##  Appendix B Perturbation equations in the fluid region

One can derive the perturbation equations in the fluid region, where $\mu=0$:

$\displaystyle r\frac{dy_{1}}{dr}=$
$\displaystyle-\left(3-\frac{U_{3}}{\Gamma}+U_{2}\right)y_{1}-\left(\frac{U_{3}}{\Gamma}-\frac{\ell(\ell+1)}{c_{1}\bar{\sigma}^{2}}\right)y_{2},$
(42) $\displaystyle r\frac{dy_{2}}{dr}=$
$\displaystyle\left(e^{2\Lambda}c_{1}\bar{\sigma}^{2}+rA_{r}\right)y_{1}-\left(U_{1}+rA_{r}\right)y_{2},$
(43)

where the variables, $y_{1}$ and $y_{2}$, are defined as

$\displaystyle y_{1}=$ $\displaystyle W,$ (44) $\displaystyle y_{2}=$
$\displaystyle\left(r\frac{d\Phi}{dr}\right)^{-1}H=c_{1}\bar{\sigma}^{2}V.$
(45)

##  Appendix C Boundary and junction conditions

The boundary condition at the stellar surface is that the Lagrangian
perturbation of pressure should be zero, i.e., $\Delta p=0$, which is
expressed as

$y_{1}-y_{2}=0,$ (46)

while the boundary condition at the center is the regularity condition, such
as

$c_{1}\bar{\sigma}^{2}y_{1}-\ell y_{2}=0.$ (47)

On the other hand, the junction conditions at the interface are the continuity
of the radial displacement, the radial and transverse tractions, and the
Lagrangian perturbation of pressure. In practice, the junction conditions at
the interface between the elastic and fluid regions [40, 41] are given as

$\displaystyle z_{1}=$ $\displaystyle y_{1},$ (48) $\displaystyle z_{2}=$
$\displaystyle U_{3}(y_{1}-y_{2}),$ (49) $\displaystyle z_{4}=$ $\displaystyle
0.$ (50)

In a similar way, the junction conditions at the interface, where the non-zero
shear modulus becomes discontinuity, e.g., the interface between the phases
composed of spherical and cylindrical nuclei, are the continuity of the
variables $z_{i}$ for $i=1-4$. Finally, one can set that $y_{1}=1$ at the
stellar surface (or at the stellar center), as a normalization condition in
the linear perturbation system.

##  Appendix D Dependence on the thickness of elastic region

In this appendix, we show how the frequencies of the interface modes depend on
the thickness of an elastic region. To see this behavior, we especially
consider the “Sp+Cy” model as discussed in Sec. IV by artificially increasing
the density between the phase composed of spherical nuclei and envelope from
$10^{12}$ g/cm3 up to $10^{14}$ g/cm3. As this transition density ($\rho_{\rm
Sp/En}$) increases, the thickness of the elastic region, $\Delta R$, decreases
from $\Delta R=0.503$ to $0.075$ km. In Fig. 11, the $i$-mode frequencies are
shown as a function of $\Delta R$, where the $i_{1}$\- and $i_{2}$-modes are
the same meaning as discussed in Sec. IV. From this figure, one can observe
that the frequency of the $i_{1}$-mode is almost independent of $\Delta R$ (or
$\rho_{\rm Sp/En}$) with a lower value of $\Delta R$, while that of the
$i_{2}$-mode strongly depends on $\Delta R$ (or $\rho_{\rm Sp/En}$). That is,
as $\Delta R$ decreases (or as $\rho_{\rm Sp/En}$ increases), the frequency of
the $i_{2}$-mode decreases and eventually disappears (or at least too small to
be determined numerically). That is, even if two interfaces exist, only one
interface mode can be excited in the stellar model with a much narrower
elastic region.

Figure 11: The frequencies of the $i_{1}$\- and $i_{2}$-modes excited due to
the presence of the elastic region composed of spherical and cylindrical
nuclei are shown as a function of the thickness of the elastic region, $\Delta
R$, by artificially increasing the density between the phase composed of
spherical nuclei and envelope, for the neutron star model with $1.4M_{\odot}$
and 12.4 km, using the EOS with $L=73.4$ MeV.

## References

  * [1] S. L. Shapiro and S. A. Teukolsky, Black Holes, White Dwarfs, and Neutron Stars: The Physics of Compact Objects (Wiley-Interscience, New York, 1983).
  * [2] P. Demorest, T. Pennucci, S. Ransom, M. Roberts, and J. Hessels, Nature 467, 1081 (2010).
  * [3] J. Antoniadis et al., Science 340, 6131 (2013).
  * [4] H. T. Cromartie et al., Nature Astronomy 4, 72 (2020).
  * [5] E. Fonseca et al., Astrophys. J. 915, L12 (2021).
  * [6] K. R. Pechenick, C. Ftaclas, and J. M. Cohen, Astrophys. J. 274, 846 (1983).
  * [7] D. A. Leahy and L. Li, Mon. Not. R. Astron. Soc. 277, 1177 (1995).
  * [8] J. Poutanen and M. Gierlinski, Mon. Not. R. Astron. Soc. 343, 1301 (2003).
  * [9] D. Psaltis and F. Özel, Astrophys. J. 792, 87 (2014).
  * [10] H. Sotani and U. Miyamoto, Phys. Rev. D 98, 044017 (2018); 98, 103019 (2018).
  * [11] H. Sotani, Phys. Rev. D 101, 063013 (2020).
  * [12] T. E. Riley et al., Astrophys. J. 887, L21 (2019).
  * [13] M. C. Miller et al., Astrophys. J. 887, L24 (2019).
  * [14] T. E. Riley et al., Astrophys. J. 918, L27 (2021).
  * [15] M. C. Miller et al., Astrophys. J. 918, L28 (2021).
  * [16] B. P. Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * [17] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018).
  * [18] H. Sotani, N. Nishimura, and T. Naito, Prog. Theor. Exp. Phys. 2022, 041D01 (2022).
  * [19] H. Sotani and S. Ota, Phys. Rev. D 106, 103005 (2022).
  * [20] H. Sotani and T. Naito, Phys. Rev. C 107, 035802 (2023).
  * [21] T. E. Strohmayer and A. L. Watts, Astrophys. J. 632, L111 (2005).
  * [22] T. E. Strohmayer and A. L. Watts, Astrophys. J. 653, 593 (2006).
  * [23] A. W. Steiner and A. L. Watts, Phys. Rev. Lett. 103, 181101 (2009).
  * [24] M. Gearheart, W. G. Newton, J. Hooker, and B.A. Li, Mon. Not. R. Astron. Soc. 418, 2343 (2011).
  * [25] H. Sotani, K. Nakazato, K. Iida, and K. Oyamatsu, Phys. Rev. Lett. 108, 201101 (2012); Mon. Not. R. Astron. Soc. 428, L21 (2013); 434, 2060 (2013).
  * [26] H. Sotani, K. Iida, and K. Oyamatsu, New Astron. 43, 80 (2016); Mon. Not. R. Astron. Soc. 464, 3101 (2017); 479, 4735 (2018).
  * [27] H. Sotani, K. Iida, and K. Oyamatsu, Mon. Not. R. Astron. Soc. 489, 3022 (2019).
  * [28] H. Sotani, K. D. Kokkotas, and N. Stergioulas, arXiv:2303.03150.
  * [29] N. Andersson and K. D. Kokkotas, Phys. Rev. Lett. 77, 4134 (1996).
  * [30] N. Andersson and K. D. Kokkotas, Mon. Not. R. Astron. Soc. 299, 1059 (1998).
  * [31] H. Sotani, K. Tominaga, and K. I. Maeda, Phys. Rev. D 65, 024010 (2001).
  * [32] H. Sotani and T. Harada, Phys. Rev. D 68, 024019 (2003); H. Sotani, K. Kohri, and T. Harada, ibid. 69, 084008 (2004).
  * [33] H. Sotani, N. Yasutake, T. Maruyama, and T. Tatsumi, Phys. Rev. D 83 024014 (2011).
  * [34] A. Passamonti and N. Andersson, Mon. Not. R. Astron. Soc. 419, 638 (2012).
  * [35] D. D. Doneva, E. Gaertig, K. D. Kokkotas, and C. Krüger, Phys. Rev. D 88, 044052 (2013).
  * [36] C. J. Krüger, W. C. G. Ho, and N. Andersson, Phys. Rev. D 92, 063009 (2015)
  * [37] H. Sotani, Phys. Rev. D 102, 063023 (2020); 102, 103021 (2020); 103, 123015 (2021).
  * [38] H. Sotani and B. Kumar, Phys. Rev. D 104, 123002 (2021).
  * [39] H. Sotani and A. Dohi, Phys. Rev. D 105, 023007 (2022).
  * [40] L. S. Finn, Mon. Not. R. Astron. Soc. 245, 82 (1990).
  * [41] S. Yoshida and U. Lee, Astron. Astrophys. 395, 201 (2002).
  * [42] A. L. Piro and L. Bildsten, Astrophys. J. 619, 1054 (2005).
  * [43] E. H. Gudmundsson, C. J. Pethick, and R. I. Epstein, Astrophys. J. 272, 286 (1983).
  * [44] K. Oyamatsu and K. Iida, Prog. Theor. Phys. 109, 631 (2003).
  * [45] K. Oyamatsu and K. Iida, Phys. Rev. C 75, 015801 (2007).
  * [46] S. Shlomo, V. M. Kolomietz, and G. Colò, Eur. Phys. J. A 30, 23 (2006).
  * [47] B. A. Li, P. G. Krastev, D. H. Wen, and N. B. Zhang, Euro. Phys. J. A 55, 117 (2019).
  * [48] H. Sotani, K. Iida, and K. Oyamatsu, Mon. Not. R. Astron. Soc. 470, 4397 (2017).
  * [49] T. Strohmayer, H. M. van Horn, S. Ogata, H. Iyetomi, and S. Ichimaru, Astrophys. J. 375, 679 (1991).
  * [50] D. A. Baiko, Mon. Not. R. Astron. Soc. 416, 22 (2011).
  * [51] D. Kobyakov and C. J. Pethick, Phys. Rev. C 87, 055803 (2013).
  * [52] D. Kobyakov and C. J. Pethick, Mon. Not. R. Astron. Soc. 449, L110 (2015).
  * [53] H. Sotani, H. Togashi, M. Takano, Mon. Not. R. Astron. Soc. 516, 5440 (2022).
  * [54] C. J. Pethick and A. Y. Potekhin, Phys. Lett. B 427, 7 (1998).
  * [55] M. E. Caplan, A. S. Schneider, and C. J. Horowitz, Phys. Rev. Lett. 121, 132701 (2018).
  * [56] C. J. Pethick, Z.-W. Zhang, and D. N. Kobyakov, Phys. Rev. C 101, 055802 (2020).
  * [57] D. N. Kobyakov and C. J. Pethick, Sov. Phys. JETP 127, 851 (2018).
  * [58] S. Yoshida and Y. Kojima, Mon. Not. R. Astron. Soc. 289, 117 (1997).
  * [59] B. Carter and H. Quintana, Proc. R. Soc. Lond. A 331, 57 (1972).
  * [60] B. L. Schumaker and K. S. Thorne, Mon. Not. R. Astron. Soc. 203, 457 (1983).
  * [61] H. Sotani, Phys. Rev. C 95, 025802 (2017).
  * [62] P. N. McDermott, H. M. van Horn, and C. J. Hansen, Astrophys. J. 325, 725 (1988).


#  Benchmark tasks for Quality-Diversity applied to Uncertain domains

Manon Flageat<EMAIL_ADDRESS>Imperial College LondonLondonU.K. ,
Luca Grillotti<EMAIL_ADDRESS>Imperial College
LondonLondonU.K. and  Antoine Cully<EMAIL_ADDRESS>Imperial College
LondonLondonU.K.

(2018)

###### Abstract.

While standard approaches to optimisation focus on producing a single high-
performing solution, Quality-Diversity (QD) algorithms allow large diverse
collections of such solutions to be found. If QD has proven promising across a
large variety of domains, it still struggles when faced with uncertain
domains, where quantification of performance and diversity are non-
deterministic. Previous work in Uncertain Quality-Diversity (UQD) has proposed
methods and metrics designed for such uncertain domains. In this paper, we
propose a first set of benchmark tasks to analyse and estimate the performance
of UQD algorithms. We identify the key uncertainty properties to easily define
UQD benchmark tasks: the uncertainty location, the type of distribution and
its parameters. By varying the nature of those key UQD components, we
introduce a set of $8$ easy-to-implement and lightweight tasks, split into $3$
main categories. All our tasks build on the Redundant Arm: a common QD
environment that is lightweight and easily replicable. Each one of these tasks
highlights one specific limitation that arises when considering UQD domains.
With this first benchmark, we hope to facilitate later advances in UQD.

Quality-Diversity optimisation, Uncertain domains, MAP-Elites, Behavioural
diversity.

††copyright: acmcopyright††journalyear: 2018††doi:
10.1145/1122445.1122456††conference: Woodstock ’18: ACM Symposium on Neural
Gaze Detection; June 03–05, 2018; Woodstock, NY††booktitle: Woodstock ’18: ACM
Symposium on Neural Gaze Detection, June 03–05, 2018, Woodstock, NY††price:
15.00††isbn: 978-1-4503-XXXX-X/18/06††ccs: Computational Methodologies
Evolutionary Robotics

##  1\. Introduction

Quality-Diversity (QD) optimisation (Cully and Demiris, 2017; Chatzilygeroudis
et al., 2021) has proven in recent years to be a promising approach across
multiple application domains (Cully et al., 2015; Gaier et al., 2017; Fontaine
et al., 2021). QD aims to find a collection $\mathcal{A}$ of diverse and high-
performing solutions to an optimisation problem. It does so by storing in
$\mathcal{A}$ solutions $i$ with high-fitness $f_{i}$ that induce novel
behaviours, quantified using behaviour descriptors $d_{i}$. In other words, QD
aims to optimise:

(1)
$\displaystyle\begin{split}\max_{\mathcal{A}}\left(\textrm{QDScore}(\mathcal{A})\right)=\max_{\mathcal{A}}\sum_{i\in\mathcal{A}}{f_{i}}\\\
\quad\text{w.r.t}\quad\forall
i\in\mathcal{A},d_{i}\in\textrm{cell}_{i}\end{split}$

Figure 1.  Illustration of the effect of uncertainty in the Redundant Arm QD
task; in this case, the uncertainty is applied to the parameters of the arm.
The descriptor is given as the end-effector $(x,y)$ position and the fitness
as the variance of the joint angles. The background figures the corresponding
QD Illusory archive (left) and Corrected archive (right) after multiple
reevaluations of each solution.

Despite their encouraging results, multiple works have high-lightened the
limitations of QD algorithms when considering uncertain domains (Justesen et
al., 2019; Flageat and Cully, 2020, 2023; Grillotti et al., 2023). A QD domain
is seen as uncertain if a solution $i$ can get different fitness $f_{i}$ and
descriptor $d_{i}$ values from one evaluation to another. In other words, in
Uncertain Quality-Diversity (UQD) (Flageat and Cully, 2023) fitness and
descriptor associated with one solution are no more fixed values but
distributions: $\tilde{f_{i}}\sim\mathcal{D}_{f_{i}}$,
$\tilde{d_{i}}\sim\mathcal{D}_{d_{i}}$. When applied to such domains, usual QD
algorithms such as MAP-Elites (Mouret and Clune, 2015) tend to keep solutions
that have been lucky during evaluation and get outlier fitness or descriptor
values. This leads to a loss in the overall performance of their final
collection. Contrary to standard QD approaches, UQD algorithms take into
consideration the distributional nature of their variables, and adapt the
problem formulation accordingly (Flageat and Cully, 2023):

(2)
$\begin{split}\max_{\mathcal{A}}\left(\textrm{QDScore}(\mathcal{A})\right)=\max_{\mathcal{A}}\sum_{i\in\mathcal{A}}{\mathbb{E}_{\tilde{f_{i}}\sim\mathcal{D}_{f_{i}}}\left[\tilde{f_{i}}\right]}\\\
\text{w.r.t}\quad\forall
i\in\mathcal{A},\mathbb{E}_{\tilde{d_{i}}\sim\mathcal{D}_{d_{i}}}\left[\tilde{d_{i}}\right]\in\textrm{cell}_{i}\end{split}$

Table 1. Performance-estimation Redundant Arm benchmark tasks. | Uncertainty location | Distribution | Parameters
---|---|---|---
Gaussian fitness noise |  Fitness after evaluation: $\begin{cases}\tilde{f}=f+\epsilon\\\ \tilde{d}=d\end{cases}$ | $\epsilon\sim\mathcal{N}(0,\sigma^{2})$ | $\sigma>0$
Bimodal fitness noise | $\epsilon\sim\alpha\mathcal{N}(0,\sigma_{1}^{2})+(1-\alpha)\mathcal{N}(-1,\sigma_{2}^{2})$ | $\sigma_{1},\sigma_{2}>0\quad\alpha\in(0,1)$
Small-Gaussian descriptor noise |  Descriptor after evaluation: $\begin{cases}\tilde{f}=f\\\ \tilde{d}=d+\epsilon\end{cases}$ | $\epsilon\sim\mathcal{N}({\bm{{0}}},\sigma^{2}{\bm{{I}}})$ | $\sigma\in(0,0.02)$
Large-Gaussian descriptor noise | $\epsilon\sim\mathcal{N}({\bm{{0}}},\sigma^{2}{\bm{{I}}})$ | $\sigma>0.1$
Bimodal descriptor noise | $\epsilon\sim\alpha\mathcal{N}({\bm{{0}}},\sigma_{1}^{2}{\bm{{I}}})+(1-\alpha)\mathcal{N}(\left[\begin{smallmatrix}1\\\ 1\end{smallmatrix}\right],\sigma_{2}^{2}{\bm{{I}}})$ | $\sigma_{1},\sigma_{2}>0\quad\alpha\in(0,1)$
Table 2. Reproducibility-maximisation Redundant Arm benchmark tasks. | Uncertainty location | Distribution | Parameters
---|---|---|---
Gaussian descriptor noise with $2$ variance values | Descriptor after evaluation: $\begin{cases}\tilde{f}=0\\\ \tilde{d}=d+\epsilon\end{cases}$ | $\begin{cases}\epsilon\sim\mathcal{N}({\bm{{0}}},\sigma_{1}^{2}{\bm{{I}}})\mbox{ if }sign(\prod_{j}{\theta^{j}})\geq 0\\\ \epsilon\sim\mathcal{N}({\bm{{0}}},\sigma_{2}^{2}{\bm{{I}}})\mbox{ otherwise}\end{cases}$ | $\sigma_{2}\gg\sigma_{1}>0$
Gaussian descriptor noise with continuous variance values | $\epsilon\sim\mathcal{N}({\bm{{0}}},(\eta*\mathbb{V}(\theta))^{2}{\bm{{I}}})$ | $\eta>0$
Table 3. Realistic Redundant Arm benchmark tasks. | Uncertainty location | Distribution | Parameters
---|---|---|---
Gaussian Phenotype $J$ noise |  Dimension $J$ of the phenotype: $\tilde{\theta}=\begin{pmatrix}\theta^{1}&\cdots&\theta^{J-1}&\theta^{J}+\epsilon&\theta^{J+1}&\cdots&\theta^{N}\end{pmatrix}$ | $\epsilon\sim\mathcal{N}\left(0,\sigma^{2}\right)$ | $\sigma>0$

Previous work in UQD (Justesen et al., 2019; Flageat and Cully, 2023, 2020;
Grillotti et al., 2023) have proposed solutions tailored for these domains, as
well as new metrics that better estimate the performance of QD algorithms when
applied to UQD domains. In this paper, we aim to propose a set of benchmark
tasks to help future developments of UQD algorithms. As first contribution, we
propose a more systematic way to define UQD benchmark tasks for a given
environment, as the ensemble of (1) an uncertainty location, (2) an
uncertainty distribution type and (3) its corresponding parameters. These
properties also constitute initial tools to analyse inherent uncertainty in
more realistic UQD tasks, allowing to better select tailored approaches. As
second contribution, we define, based on this definition, $8$ new tasks, split
into $3$ main categories. These tasks aim to highlight some of the main
limitations that arise in UQD settings. All our tasks are based on the
Redundant Arm task, a simple-to-implement and lightweight QD task that has
been widely studied across the existing literature (Cully et al., 2015;
Flageat et al., 2022). To facilitate the use of our benchmark, we make our
code available at URL https://github.com/adaptive-intelligent-robotics/QDax .

##  2\. Properties of UQD tasks

One can define multiple UQD benchmark tasks for a given environment by varying
the properties of the uncertainty induced in the environment. We identified 3
main uncertainty properties: the uncertainty location, the type of
distribution, and the parameters of this distribution. In this section, we
analyse those properties and explain how they can be used to define UQD
benchmark tasks. They can also be used as tools to analyse more realistic UQD
tasks.

###  2.1. Uncertainty location

As explained above, in UQD tasks, fitnesses and descriptors are no more fixed
values but distributions. These distributions can arise from uncertainty at
different levels within the environment. The clear definition of this
uncertainty location is critical to define benchmark tasks and study UQD
tasks. Previous work in uncertain evolutionary algorithms proposed a
classification that separates prior noise, applied to the phenotype before the
fitness evaluation, and posterior noise, applied to the fitness after the
evaluation (Neumann et al., 2020). In the case of UQD, we propose to decompose
further the prior noise category to encompass dynamic and applied problems as
considered in QD. In simple QD environments such as Rastrigin or the Redundant
Arm, one could consider adding (1) noise on the fitness and descriptor
reading, or (2) noise on the phenotype (i.e. the parameters in Rastrigin and
the joint angles in the Arm), which corresponds directly to the prior and
posterior noise categories. However, in more complex settings, such as robot
control environments (Cully et al., 2015), one could consider adding (1) noise
on the fitness and descriptor reading, (2) noise on the phenotype, i.e. the
controller parameters, (3) noise on the initial position of the robot, and (4)
stochasticity in the dynamics of the environment. Similar locations could be
considered in closed-loop QD-RL environments (Flageat et al., 2022; Chalumeau
et al., 2022).

###  2.2. Distribution Type

Another key property of UQD tasks is the distribution of the noise or
stochasticity. Most previous work in UQD used Gaussian or Uniform
distributions. However, more complex distributions such as skewed, multi-
modal, or non-continuous distributions can lead to interesting properties as
will be discussed later in this paper. Depending on the uncertainty location
within the environment, the structure of this distribution will then be shaped
and modified to give the final fitness and descriptor distributions, leading
to complex uncertainty structures.

###  2.3. Distribution Parameters

Varying the parameters of a distribution may also have a significant impact on
the performance of UQD algorithms. For example, gaussian descriptor
distributions with a very low variance may have a negligible impact on an
algorithm’s results; indeed, if an individual’s descriptor has a low variance,
then most re-evaluations of that descriptor will land in the same behaviour
niche (i.e. MAP-Elites cell). Conversely, a really large fitness variance can
make the comparison between two fitness values almost non-meaningful and thus
prevent convergence.

##  3\. Redundant Arm UQD benchmark

###  3.1. Redundant Arm

Figure 2.  Distributions in the Performance-estimation tasks.

The Redundant Arm (Cully et al., 2015; Cully and Demiris, 2017;
Chatzilygeroudis et al., 2021), illustrated in Figure 1, is one of the most
commonly used QD benchmark task. It aims to learn how to achieve every
reachable point with a planar $N$-DoF robotic Arm:

  * •

Controller:  real-valued angles of the $N$-DoF:
$\theta_{i}=(\theta_{i}^{j})_{1\leq j\leq N}$.

  * •

Descriptor:  defined as the $(x,y)$ position of the end-effector, after
applying the control values in the joints.

  * •

Fitness:  in this paper, we consider two versions of this task:

    1. (1)

fitness $0$ for all individuals.

    2. (2)

fitness as negative variance of the $N$ joint angles $(\theta_{i}^{j})_{1\leq
j\leq N}$:
$f_{i}=-\mathbb{V}(\theta_{i})=-\frac{1}{N}\sum_{j=1}^{N}{(\theta_{i}^{j}-\overline{\theta_{i}})^{2}}$.

###  3.2. UQD benchmark tasks

In the following, we propose a set of benchmark UQD tasks in the Redundant Arm
environment, chosen to illustrate some of the common challenges in UQD. As
highlighted by Flageat and Cully (2023), UQD approaches usually tackle two
main problems:

  * •

Performance estimation: correctly estimating the expected fitness
$\mbox{E}_{f_{i}\sim\mathcal{D}_{f}}\left[f_{i}\right]$ and expected
descriptor $\mbox{E}_{d_{i}\sim\mathcal{D}_{d}}\left[d_{i}\right]$ of
solutions from their evaluations.

  * •

Reproducibility maximisation: encouraging solutions that are more likely to
get similar fitness and descriptor values from one evaluation to another. This
corresponds to maintaining solutions with lower variance of fitness
$\mbox{V}_{f_{i}\sim\mathcal{D}_{f}}\left[f_{i}\right]$ and variance of each
descriptor dimension
$\mbox{V}_{d^{j}_{i}\sim\mathcal{D}_{d^{j}}}\left[d^{j}_{i}\right]$. This
issue only arises in tasks where solutions can have different variances.

In Section 4, we present $5$ benchmark tasks designed to assess the
Performance estimation ability of UQD algorithms; and in Section 5, we
introduce $2$ other tasks to assess their ability to select reproducible
solutions. We also propose in Section 6 a more realistic task that combines
both problems. For each section, we give some preliminary results obtained by
MAP-Elites and MAP-Elites-sampling (with $30$ samples) with $100\times 100$
grid. These approaches are used as baselines in existing UQD works (Flageat
and Cully, 2020; Justesen et al., 2019; Flageat and Cully, 2023; Grillotti et
al., 2023).

##  4\. Performance Estimation tasks

In this section, we propose benchmark tasks that assess the ability of UQD
approaches to estimate performance. We consider two cases: tasks where only
the fitness is a distribution and tasks where only the descriptor is a
distribution. One can easily combine these cases to get complex tasks with
noise on both. These tasks all build on the Redundant Arm with non-zero
fitness from Section 3.1. We introduce each task below, give corresponding
equations in Table. 1 and Figure 2, and some results in Figures 4 and 5.

Figure 3.  Illustration of the effect of noise applied on the $6^{\text{th}}$
joint of 2 distinct arm configurations $\theta_{1}$ (left) and $\theta_{2}$
(right). The red crosses indicate the sampled descriptors, and the blue cross
their average. The descriptor reproducibility achieved by $\theta_{2}$ is way
higher than the one achieved by $\theta_{1}$.

###  4.1. Gaussian fitness noise

Gaussian noise on the fitness is the most commonly used noise distribution in
previous UQD works. Here, the variance parameter $\sigma^{2}$ for this
distribution can take almost any value. The results shows the benefit of
sampling-based approaches in this case (Fig. 4).

Figure 4.  Results on the Performance-estimation benchmark tasks. Each column
corresponds to one of the tasks from Table 1. For each baseline, we give the
Illusory (left - dashed line) and Corrected (left - plain line) QD-Score and
Coverage (Flageat and Cully, 2023; Grillotti et al., 2023). Solid lines are
the median across 5 replications, and shaded area the quartiles. We also give
the Loss QD-Score (right - top) and Loss Coverage (right - bottom) from
Flageat et al. (2023), that quantify the estimation capability of each
baseline.

###  4.2. Bimodal fitness noise

One of the particularity of Gaussian distributions, as proposed in the
previous section, is that sampling is always beneficial. Here, we propose a
task where too low number of samples can actually be detrimental. We define a
bimodal distribution with a distant low-probability second mode. This second
mode is really detrimental for low number of samples, leading to an
interesting phenomenon where MAP-Elites outperforms MAP-Elites-sampling (Fig.
4).

###  4.3. Small-Gaussian descriptor noise

Gaussian noise on the descriptor is the most commonly used noise distribution
in previous UQD works. The value of the variance $\sigma^{2}$ for this
distribution is critical. Here, we consider small values that lead individuals
to stay with high probability within the same behavioural niche (i.e. MAP-
Elites cell). The results show the benefit of sampling-based approaches in
this case (Fig. 4).

###  4.4. Large-Gaussian descriptor noise

In this second version, we choose the standard deviation $\sigma$ to be
sufficiently high to lead solutions to change behavioural niche with high
probability. This task makes the correct estimation of descriptor highly
complex. For instance, with a MAP-Elites grid of size $100\times 100$ and a
value of $\sigma$ higher than $0.1$, most descriptor estimates with fewer than
$100$ samples will not be placed in the appropriate cell. The results
highlight the limitation of MAP-Elites and MAP-Elites-sampling in this case
(Fig. 4 and 5).

###  4.5. Bimodal descriptor noise

The Bimodal distribution is the equivalent of the Bimodal Fitness Noise from
Section 4.2 for the descriptor. In this task, the noise takes the form of a
bimodal distribution in 2 dimensions with a distant low-probability second
mode (see Fig. 2). This distant second mode makes it difficult to estimate
expected descriptors based on a low amount of samples, which is why MAP-
Elites-sampling achieves the worst results in this task (Fig. 4 and 5).

##  5\. Reproducibility maximisation tasks

In this section, we propose $2$ benchmark tasks that assess the ability of UQD
approaches to keep solutions with high reproducibility. To simplify the
problem, we consider only the reproducibility in term of descriptor, and we
only focus on the Redundant Arm with zero-fitness from Section 3.1. We
introduce each task below, give corresponding equations in Table. 2, and some
results in Figures 6 and 7. To quantify the performance of approaches, we
display the Reproducibility-Score as defined by Flageat and Cully (2023).
Additionally, inspired by Grillotti et al. (2023), we add as baseline MAP-
Elites-sampling-Reproducibility, which directly optimises for the
reproducibility of solutions.

###  5.1. Gaussian descriptor noise with $2$ possible variance values

We first consider a simple task where there are two distinct possible values
of variance for the descriptor: one small $\sigma_{1}^{2}$ and one large
$\sigma_{2}^{2}$. Thus, solutions ending in the same cell can be either highly
reproducible or poorly reproducible. In this task, we attribute variance based
on phenotype: we give low variance $\sigma_{1}^{2}$ to solutions with positive
product of joint angles, and high variance $\sigma_{2}^{2}$ to the others.
Ideally, as the fitness is the same (zero) for all solutions, UQD approaches
are expected to prefer low-variance solutions to high-variance ones. However,
the results show that none of the two considered approaches manages to select
reproducible solutions across the descriptor space, unlike MAP-Elites-
sampling-Reproducibility (Grillotti et al., 2023) (Fig. 6 and 7). UQD
algorithms require more complex mechanisms to be able to effectively solve
this task (Justesen et al., 2019; Flageat and Cully, 2023; Grillotti et al.,
2023).

Figure 5.  Illusory (top) and Corrected (bottom) archives for the two
Performance-estimation descriptor noise tasks.  Figure 6.  Results on the
Reproducibility-maximisation and Realistic-combination benchmark tasks. Each
column corresponds to one of the task defined in Table 2 or 3. For each
baseline, we give the Illusory Coverage (dashed line - left) and Corrected
Coverage (plain line - left). We also display the Reproducible-Score from
Flageat and Cully (2023) that quantifies the ability of baselines to maintain
reproducible solutions. It is computed based on the normalised variance
$\sigma^{2}_{norm}$ as the sum over the Corrected archive of
$1-\sigma^{2}_{norm}$. Solid lines are the median across 5 replications, and
shaded area the quartiles.  Figure 7.  Reproducibility archives (Flageat and
Cully, 2023) for benchmark tasks of Reproducibility-maximisation and
Realistic-combination.

###  5.2. Gaussian descriptor noise with continuous variance values

We propose a more complex version of the previous task, with the aim to refine
the study of UQD algorithms in reproducibility-maximisation tasks. In this
task, the variance can take any values, continuously, within a range. Here, we
attribute variance proportionally to the variance of the joint angles.
Ideally, as the fitness is the same (zero) for all solutions, UQD approaches
are expected to select the most reproducible solutions for each cell. However,
as for the previous task, results in Figure 6 show that neither MAP-Elites nor
MAP-Elites-sampling favour reproducible solutions.

##  6\. Realistic combination task

In addition to the targeted tasks proposed in the previous sections, we
propose here one more realistic task that combines both the issues of
performance estimation and reproducibility maximisation: the Gaussian
Phenotype $J$ noise. We give the corresponding equations in Table 3 and some
results in Figures 6 and 7.

###  6.1. Gaussian Phenotype $J$ noise

In this task, we propose to change the uncertainty location to consider the
phenotype (i.e. the joint angles). Noise is applied on only one single pre-
defined joint angle, and all solutions have zero-fitness. UQD algorithms are
expected to maximise the reproducibility of the solutions while also providing
accurate estimates of their expected descriptors. In this task, the best
solutions are fully reproducible, i.e. have a descriptor variance of 0, as
depicted on Figure 3. The results show that MAP-Elites and MAP-Elites-sampling
do not manage to find these solutions (Fig. 6 and 7).

##  7\. Conclusion and Discussion

This paper presents the first set of Uncertain Quality-Diversity benchmark
tasks. Each of these tasks aims to highlight one limitation of existing QD
algorithms in uncertain domains and points toward interesting directions of
improvement. They are all defined for the Redundant Arm, a commonly used QD
environment. In total, we propose $8$ tasks that can be split into $3$ main
categories: Performance-estimation, Reproducibility-maximisation and Realistic
tasks. The first two correspond to the two main issues that arise when
considering Uncertain domains, as highlighted in previous work. We hope that
this paper will constitute a meaningful first benchmark for later advances in
UQD and that it will lead to more similar tasks to assess QD algorithms in
Uncertain domains.

## References

  * (1)
  * Chalumeau et al. (2022) Felix Chalumeau, Thomas Pierrot, Valentin Macé, Arthur Flajolet, Karim Beguir, Antoine Cully, and Nicolas Perrin-Gilbert. 2022\.  Assessing Quality-Diversity Neuro-Evolution Algorithms Performance in Hard Exploration Problems.  _arXiv preprint arXiv:2211.13742_ (2022).
  * Chatzilygeroudis et al. (2021) Konstantinos Chatzilygeroudis, Antoine Cully, Vassilis Vassiliades, and Jean-Baptiste Mouret. 2021\.  Quality-Diversity Optimization: a novel branch of stochastic optimization.  In _Black Box Optimization, Machine Learning, and No-Free Lunch Theorems_. Springer, 109–135.
  * Cully et al. (2015) Antoine Cully, Jeff Clune, Danesh Tarapore, and Jean-Baptiste Mouret. 2015.  Robots that can adapt like animals.  _Nature_ 521, 7553 (2015), 503–507.
  * Cully and Demiris (2017) Antoine Cully and Yiannis Demiris. 2017.  Quality and diversity optimization: A unifying modular framework.  _IEEE Transactions on Evolutionary Computation_ 22, 2 (2017), 245–259.
  * Flageat et al. (2023) Manon Flageat, Felix Chalumeau, and Antoine Cully. 2023\.  Empirical analysis of PGA-MAP-Elites for Neuroevolution in Uncertain Domains.  _ACM Transactions on Evolutionary Learning_ 3, 1 (2023), 1–32.
  * Flageat and Cully (2020) Manon Flageat and Antoine Cully. 2020.  Fast and stable MAP-Elites in noisy domains using deep grids. In _ALIFE 2020: The 2020 Conference on Artificial Life_. MIT Press, 273–282.
  * Flageat and Cully (2023) Manon Flageat and Antoine Cully. 2023.  Uncertain Quality-Diversity: Evaluation methodology and new methods for Quality-Diversity in Uncertain Domains.  _arXiv preprint arXiv:2302.00463_ (2023).
  * Flageat et al. (2022) Manon Flageat, Bryan Lim, Luca Grillotti, Maxime Allard, Simón C Smith, and Antoine Cully. 2022\.  Benchmarking Quality-Diversity Algorithms on Neuroevolution for Reinforcement Learning.  _arXiv preprint arXiv:2211.02193_ (2022).
  * Fontaine et al. (2021) Matthew C Fontaine, Ruilin Liu, Ahmed Khalifa, Jignesh Modi, Julian Togelius, Amy K Hoover, and Stefanos Nikolaidis. 2021.  Illuminating mario scenes in the latent space of a generative adversarial network. In _Proceedings of the AAAI Conference on Artificial Intelligence_ , Vol. 35. 5922–5930.
  * Gaier et al. (2017) Adam Gaier, Alexander Asteroth, and Jean-Baptiste Mouret. 2017\.  Aerodynamic design exploration through surrogate-assisted illumination. In _18th AIAA/ISSMO multidisciplinary analysis and optimization conference_. 3330\.
  * Grillotti et al. (2023) Luca Grillotti, Manon Flageat, Bryan Lim, and Antoine Cully. 2023.  Don’t Bet on Luck Alone: Enhancing Behavioral Reproducibility of Quality-Diversity Solutions in Uncertain Domains.  _arXiv preprint arXiv:2304.03672_ (2023).
  * Justesen et al. (2019) Niels Justesen, Sebastian Risi, and Jean-Baptiste Mouret. 2019\.  Map-elites for noisy domains by adaptive sampling. In _Proceedings of the Genetic and Evolutionary Computation Conference Companion_. 121–122.
  * Mouret and Clune (2015) Jean-Baptiste Mouret and Jeff Clune. 2015.  Illuminating search spaces by mapping elites.  _arXiv preprint arXiv:1504.04909_ (2015).
  * Neumann et al. (2020) Frank Neumann, Mojgan Pourhassan, and Vahid Roostapour. 2020\.  Analysis of evolutionary algorithms in dynamic and stochastic environments.  _Theory of evolutionary computation: recent developments in discrete optimization_ (2020), 323–357.