# Robotic in-hand manipulation with relaxed optimization

Ali Hammoud1, Valerio Belcamino2, Quentin Huet1, Alessandro Carfi2, Mahdi
Khoramshahi1,
Veronique Perdereau1 and Fulvio Mastrogiovanni2 This work was supported by the
European project Index1A. Hammoud and V. Perdereau are with Sorbonne
Universite, Institut des Systemes Intelligents et de Robotique, ISIR, F-75005
Paris, France ali.hammoud, quentin.huet, mahdi.khoramshahi,
<EMAIL_ADDRESS>Belcamino, A. Carfì and F.
Mastrogiovanni are with TheEngineRoom, Department of Informatics,
Bioengineering, Robotics, and Systems Engineering, University of Genoa, Via
Opera Pia 13, 16145, Genoa, Italy<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>2024 IEEE.
Personal use of this material is permitted. Permission from IEEE must be
obtained for all other uses, in any current or future media, including
reprinting/republishing this material for advertising or promotional purposes,
creating new collective works, for resale or redistribution to servers or
lists, or reuse of any copyrighted component of this work in other works.

###### Abstract

Dexterous in-hand manipulation is a unique and valuable human skill requiring
sophisticated sensorimotor interaction with the environment while respecting
stability constraints. Satisfying these constraints with generated motions is
essential for a robotic platform to achieve reliable in-hand manipulation
skills. Explicitly modelling these constraints can be challenging, but they
can be implicitly modelled and learned through experience or human
demonstrations. We propose a learning and control approach based on
dictionaries of motion primitives generated from human demonstrations. To
achieve this, we defined an optimization process that combines motion
primitives to generate robot fingertip trajectories for moving an object from
an initial to a desired final pose. Based on our experiments, our approach
allows a robotic hand to handle objects like humans, adhering to stability
constraints without requiring explicit formalization. In other words, the
proposed motion primitive dictionaries learn and implicitly embed the
constraints crucial to the in-hand manipulation task.

##  I INTRODUCTION

Humans’ ability to manipulate objects with dexterity is essential for
interacting with their surroundings, allowing them to grasp, explore, and
reorient objects. These skills are developed through a lifelong learning
process that involves observing other people’s behaviour and personal attempts
and failures. A primary goal of human-robot interaction is to integrate robots
into human-centred environments. However, the effectiveness of this
integration depends on the robot’s ability to move and operate in a human-like
manner[1]. Therefore, robots should be trained to manipulate unfamiliar
objects and apply their prior knowledge to new situations (as displayed in
Fig. 1). Additionally, robots should be able to learn new manipulation
techniques by observing the actions of other agents. A robotic platform can
achieve this by incorporating advanced perception tools and adaptable learning
techniques to generate new smooth and dexterous manipulation
operations.Planning a dexterous manipulation requires defining appropriate
finger trajectories to reach a predetermined target configuration. Robotic
manipulation planning strategies can be categorized into two groups: data-
driven and analytical [2, 3, 4, 5, 6, 7, 8]. In data-driven approaches,
dexterous manipulation models are trained by robot trial and error or by
observing human demonstrations [2, 3, 4]. On the other hand, analytical
solutions are based on robotics principles; complex tasks are divided into
sets of elementary actions that solve sections of the more significant
challenge [5, 6, 7].

Figure 1: The Shadow Hand manipulating a cylinder. The reference frame for the
palm ($R_{P}$) is shown, with the $x$, $y$, and $z$ axes in orange, green, and
blue arrows, respectively.

Most data-driven in-hand manipulation path planning approaches rely on Dynamic
Movement Primitives (DMP) [2]. By modifying a simple linear dynamical system
with a non-linear component, DMP enables the creation of smooth movements of
any shape [9]. Therefore, DMP can learn from humans when data from human
demonstrations are used to identify the non-linear components of a DMP
system.On the other hand, analytical solutions for in-hand manipulation
planning typically require modelling the robotic hand and its surroundings.
Robot hand actions are decomposed into atomic processes, such as in-grasp
reorientation and finger reallocation. In-grasp reorientation involves moving
the object relative to the palm without changing the contact points between
the fingers and the surface, for example, turning a dial. Sundaralingam and
Hermans (2017) proposed an efficient solution to this problem with purely
kinematic trajectory optimization [5]. In finger reallocation, the robot moves
a single finger to a new contact location on the object while the remaining
fingers maintain a stable grasp. Sundaralingam and Hermans (2018) [6] and Fan
et al. (2017)[7] both used geometric approaches to offer solutions to finger
gaiting. Classical methods provide predictable planning outcomes because the
designer sets the constraints. However, fully modelling manipulation
restrictions is complex, so these solutions can only approach atomic
operations rather than the entire in-hand manipulation task. Both analytical
and data-driven approaches have their limitations. Data-driven approaches
require large training datasets and long training times. In contrast, the
analytical approaches must model complex constraints to generate an entire in-
hand manipulation trajectory.

In our previous work[8], we introduced a novel approach that aims to balance
the strengths and weaknesses of data-driven methods by extracting human
skills. This is achieved through the creation of a dictionary of motion
primitives, which are derived from human demonstrations. This dictionary
comprises sparse vectors representing manipulation primitives and can be
combined to generate human-like in-hand manipulation trajectories. In our
subsequent research [10], addressing the generation of trajectories emerged as
a pivotal concern. We leveraged a dictionary as the foundation for a more
flexible path-planning algorithm. Our path planner successfully identified the
necessary motion primitives to transition the hand from configuration A to
configuration B while adhering to constraints such as finger reachable space,
finger collisions, and the number of contact points with the object. Although
these constraints were not explicitly imposed during the optimization process,
they were implicitly extracted from human demonstrations. While this approach
achieved human-like fingertip trajectories respecting environmental
constraints, it did so without explicitly considering the object’s pose.
However, effective and realistic robotic manipulation requires considering the
object’s dynamics and kinematics.

In this work, we address the limitations of our previous solution by adopting
a human in-hand manipulation solution that considers the relationship between
the object pose and the fingertip positions. This solution can generate
consistent fingertip trajectories given a desired change in the object pose.
To accomplish this result, the data from human demonstrations must include the
manipulated object’s full pose (position and orientation) and fingertip
positions. By integrating this comprehensive data, we can harness the learned
primitives to plan in-hand manipulations that adhere to stability constraints
while effecting the desired transformations on the handled object.
Additionally, our proposed approach tackles the computational problem often
associated with data-driven approaches by building a primitive dictionary from
a small sample of demonstrations.

The rest of this paper is structured as follows. Section II discusses the
concept of a primitives dictionary and how it relates to in-hand manipulation.
Section III explains path planning using an in-hand manipulation dictionary.
Section IV outlines the actions needed to incorporate an in-hand manipulation
dictionary into the path-planning problem, including data collection, model
training, and model testing. In Section V, we analyze robot simulation and
experimental results. Finally, the last section presents discussions,
conclusions and future work.

##  II Primitives Dictionary

Dictionaries of primitives are commonly used tools in literature for
clustering and reducing the dimensionality of datasets. This method has broad
applications, including computer vision [11], document clustering [12],
astronomy [13], and motion planning [14]. When focusing on human actions, the
dictionary is considered to be composed of motion primitives. Although there
is no universal definition for motion primitives, they are often described as
simple motions whose concatenation leads to complex human actions [15][16].
Furthermore, primitives’ dictionaries can be learned from training data using
the non-negative matrix factorization method (NMF) [17].

Let us begin with a set of demonstrations organized in the matrix form as
$V\in\mathbb{R}^{n\times m}$, where $n$ and $m$ respectively represent the
length and the number of demonstrations and each element of $V$ is $\geq 0$.
Each column of $V$ denotes a single demonstration and we can apply the NMF
algorithm that leads to the following decomposition:

$V=WH$ (1)

Where $W\in\mathbb{R}^{n\times l}$ is a matrix representing the extracted
primitives and each column of $W$ represents a single primitive. Given the
ability of NMF to reduce the data dimensionality, the number of extracted
primitives $l$ is considerably lower than the number of demonstrations $m$
(i.e. $l<m$). Furthermore, $H\in\mathbb{R}^{l\times m}$ is also a matrix that
represents the corresponding activation matrix and contains the weights to
combine the primitives. For instance, each $i$th column of $H$ encodes the
weights that can be used to reconstruct the $i$th demonstration with a linear
combination of the primitives. These two matrices $W$ and $H$ only contain
non-negative elements

Using the extracted primitives, we can also generate new samples: given the
desired behaviour to represent ($v$), it is possible to reconstruct it with a
dictionary of primitives ($W$) by combining them with a set of weights ($h$):

$v=Wh$ (2)

Each weight $h_{j}$ belonging to the activation vector $h$ influences to which
extent every primitive $W_{j}$ contributes to the reconstruction of $v$.

###  II-A In-Hand Manipulation Primitives

Object manipulation aims to achieve a specific object pose or trajectory
through finger motions. Therefore, in-hand manipulations can be described as
the evolution of fingertips’ positions and the object’s pose over time,
starting from an initial configuration and ending at the desired pose.
Therefore, the problem can be decomposed into two main components: the
trajectories of the fingers and the object’s trajectory. With this assumption
and recalling Eq. 2, we can approach the problem by describing the overall
trajectory as:

$v=\left[{\begin{array}[]{c}P(1)\\\ \vdots\\\ P(k)\\\ \vdots\\\ P(N)\\\
\end{array}}\right]\\\ $ (3)

where $v\in\mathbb{R}^{21N}$, $N$ is the number of time steps and 21, as we
will see, is the number of features. In fact, each element $P(k)$ of $v$
represents the fingertips and object pose at time $k\Delta t$ and is defined
as follows:

$P(k)=\left[{\begin{array}[]{c}P_{1}(k\Delta t)\\\ P_{2}(k\Delta t)\\\
P_{3}(k\Delta t)\\\ P_{4}(k\Delta t)\\\ P_{5}(k\Delta t)\\\ P_{O}(k\Delta
t)\\\ \end{array}}\right]\\\ $ (4) $\textit{s.t.}\hskip
2.84544ptk\in[1,\dots,N]$

with $P(k)\in\mathbb{R}^{21}$ and where

$P_{i}(k\Delta t)=\left[{\begin{array}[]{c}x_{i}(k\Delta t)\\\ y_{i}(k\Delta
t)\\\ z_{i}(k\Delta t)\\\ \end{array}}\right]\\\ $ (5)

and

$P_{O}(k\Delta t)=\left[{\begin{array}[]{c}x_{O}(k\Delta t)\\\ y_{O}(k\Delta
t)\\\ z_{O}(k\Delta t)\\\ \psi_{O}(k\Delta t)\\\ \theta_{O}(k\Delta t)\\\
\phi_{O}(k\Delta t)\\\ \end{array}}\right]\\\ $ (6)

Eq. 5 and 6 represent the fingertip position in Cartesian space and the object
pose, respectively. Therefore, there are 21 features at each time instant.
Each of the five fingers is represented by x, y, and z coordinates, and the
object pose includes its x, y, and z position as well as roll ($\psi$), pitch
($\theta$), and yaw ($\phi$).

Given these definitions, the problem of in-hand manipulation involves finding
the fingertips and object trajectories from an initial pose $\hat{P}(1)$ to a
final desired pose $\hat{P}(N)$.

This formalization leads to the definition of $W$ as follows:

$W=\left[{\begin{array}[]{c}W(1)\\\ \vdots\\\ W(k)\\\ \vdots\\\ W(N)\\\
\end{array}}\right]\\\ $ (7)

where $W\in\mathbb{R}^{21N\times l}$ and $l$ is the number of primitives. Each
element $W(k)\in\mathbb{R}^{21\times l}$ represents the $l$ primitives at the
$k$th of the 21 features describing the object and fingertips pose. As a
result of this formalization, we can write Eq. 2 as below:

$\left[{\begin{array}[]{c}P(1)\\\ \vdots\\\ P(N)\\\
\end{array}}\right]=\left[{\begin{array}[]{c}W(1)\\\ \vdots\\\ W(N)\\\
\end{array}}\right]h$ (8)

We can see from this representation how each of the l values of $h$ acts as a
weight, assigning importance to each primitive of $W$ during the
reconstruction of $v$.

##  III Generation of Manipulations

Once the dictionary of primitives has been created from the dataset of
demonstrations, we can combine the primitives using weights to produce new in-
hand manipulation trajectories. However, the created trajectories must respect
the constraints of the robotic application. Before introducing the constraints
and describing the process of generating in-hand manipulation, it is important
to note that a few assumptions about the manipulation task limit the proposed
approach:

  * •

Only two factors can influence the pose of the object: the robot and gravity.

  * •

Both the robotic hand and the object are rigid.

  * •

The initial and final object grasps are stable.

  * •

Only the fingertips make contact with the object.

Now that the problem has been formalized and the initial assumption has been
established, we can introduce the process of determining the proper weights to
solve a particular manipulation. This is achieved through an optimization
process that aims to minimize the following cost function:

$h=\arg\underset{h}{\min}\lVert P(1)-\hat{P}(1)\rVert^{2}+\lambda\lVert
P(N)-\hat{P}(N)\rVert^{2}$ (9)

and by substituting Eq. 8 into Eq. 9, we get:

$h=\arg\underset{h}{\min}\lVert W(1)h-\hat{P}(1)\rVert^{2}+\lambda\lVert
W(N)h-\hat{P}(N)\rVert^{2}$ (10)

where the first and the second terms minimize the difference between the
desired and achieved fingertips positions and object pose at steps $1$ and
$N$, respectively. Lambda $\lambda$ is the scalar weight fine-tuning the
trade-off between the two cost components.

The optimization process must consider other constraints linked to the robotic
hand kinematics while solving the problem. The i-th fingertip should be in the
reachable workspace $\mathbb{P}_{i}$. That is mathematically defined as the
set of points in three-dimensional space that the fingertip can reach. This
set is typically constrained by the physical limitations of the robotic arm or
hand controlling the fingertip.

$P_{i}(k\Delta t)\in\mathbb{P}_{i},\forall i\in[1,5],\forall k\in[1,N]$ (11)

Additionally, each fingertip must comply with kinematic constraints on its
instantaneous velocity.

$\dot{P}_{i,min}\leq\dot{P}_{i}(k\Delta t)\leq\dot{P}_{i,max},\forall
i\in[1,5],\forall k\in[1,N]$ (12)

Other constraints must also be considered. For example, the fingertips should
not overlap during the motion. Additionally, during manipulation, the object
must maintain contact with at least two fingertips of the Shadow Hand to
ensure a stable grip, leveraging the compliant nature of the fingertips. This
soft compliance allows for slight deformation upon contact, enhancing grasping
capabilities, especially for irregular objects. Optimal stability is achieved
by distributing forces through two fingertips, creating a moment about the
object’s centre of mass to prevent slippage or rotation. While any two fingers
can establish contact, stability is maximized when forces oppose or act along
different axes, simplifying grasp planning while maintaining effective
manipulation [18]. We can express the first constraint by checking the
Euclidean distances between fingertips for all the manipulation time:

$\min_{\forall j\in(i,5]}\lVert P_{i}(k\Delta t)-P_{j}(k\Delta t)\rVert\neq
0,$ (13) $\forall i\in[1,4],\forall k\in[1,N]$

The constraint on the points of contact is represented by the distance between
the fingertip ($P_{i}(k\Delta t)$) and the object surface, represented as a
point-cloud ($O(k\Delta t)$) of 3D points ($o_{c}(k\Delta t)$, with
$c\in[1,|O|]$). This representation presupposes a process for computing the
point cloud influenced by the object’s shape, dimension, and spatial pose. In
light of this representation, the distance ($d_{i}(k\Delta t)$) between the
$i$th fingertip and the object is considered as the distance between the
fingertip and the closest point of the object point cloud:

$d_{i}(k\Delta t)=\min_{\forall c\in[1,|O|]}\lVert P_{i}(k\Delta
t)-o_{c}(k\Delta t))\rVert$

and the subset of fingertips in contact with the object ($P_{c}(k\Delta t)$)

$P_{c}(k\Delta t)=\\{P_{j}(k\Delta t)\hskip 2.84544pt|\hskip
2.84544ptd_{j}(k\Delta t)\leq\tau\\}$ (14)

where $\tau$ is a threshold that defines the maximum distance between a
fingertip and an object to be considered in contact. Based on this definition,
at each time step $k$, the cardinality of this subset must always be greater
than 2 for soft fingertips to hold the object during manipulation.

$|P_{c}(k\Delta t)|\geq 2,\hskip 2.84544pt\forall k\in[1,N]$ (15)

###  III-A Relaxed Optimization Problem

Humans naturally adhere to the constraints presented in the previous section.
Since our approach uses human data to train the motion primitives dictionary,
we hypothesize that trajectories generated from the human manipulation
dictionary will follow these constraints without explicitly including them in
the optimization problem. This hypothesis applies to all constraints, except
for the fingertips’ velocity, because of the differences in velocity ranges
between humans and robots. Therefore, the optimization process should find the
weights $h$ that minimize Eq. 10 while considering the fingertips’ velocity
constraint expressed in Eq. 12. The optimization process does not take into
account the other constraints (Eq. 11,13,15). However, the generated
trajectories will be evaluated to determine if they respect them.

##  IV Implementation

Based on the problem statement introduced in Section II and the optimization
procedure detailed in Section III, human demonstrations of in-hand
manipulations are required to generate in-hand manipulation. These
demonstrations should include both the 3D positions of the fingertips and the
object pose.

###  IV-A Data Acquisition

The acquisition of human demonstrations has been performed using a motion
capture (MoCap) system comprising six OptiTrack Flex-3
cameras111https://optitrack.com/cameras/flex-3/ (see Fig. 2). These infrared
cameras can track the position of small reflective objects within the field of
view of at least three cameras (see Fig. 3). The same concept can be extended
to determine orientation, but, in this case, it is necessary to define a rigid
body by attaching three or more markers to the object’s surface. The
Motive222https://www.optitrack.com/software/motive/ software acquires the
information mentioned above at a frequency of 100 Hz. In case of interruption
in the tracking, the user can fill in the missing parts of the sequences
through cubic interpolation.

Figure 2: Six Flex-3 cameras arranged on a one-meter cube frame.

We positioned markers to record the human demonstrations as shown in Fig. 3.
We used one reflective sphere for each finger on the distal phalanges and a
4-marker rigid body positioned on the hand back to track the palm pose.

For the experiments, we considered two differently shaped objects: a cube with
a 5-centimeter edge and a cylinder with a diameter and height of 5 cm. We used
the same 4-marker rigid body configuration to track the objects, as shown in
Fig. 3, ensuring that the markers’ centroid corresponded to their axis of
symmetry. This choice allowed us to determine the position of the centre of
the objects by applying a simple translation in the local system of reference.

Due to the limited number of cameras and the close positioning of the markers
during manipulation, the markers could easily occlude each other, resulting in
poor tracking results. To solve this issue, we restricted the manipulations to
a specific hand pose in space, i.e., the palm facing up in the middle of the
tracking area, and disallowed any significant wrist rotation. These
precautions allowed us to achieve high measurement precision by minimizing
occlusions and collisions.

We conducted six 5-minute trials for each object using the described setup.
Each trial included rotations and translations along multiple axes and
required the object’s entire weight to be held exclusively by the fingertips.
As previously mentioned, to ensure that the manipulations were solely
resulting from finger motions and to optimize the recordings’ quality, wrist
motions were not allowed. In addition, the rotations of the objects were
limited by the markers mounted on them and the need to avoid occlusion.

The recording phase resulted in two datasets of in-hand manipulations, one for
each object shape. Each of the two datasets contains 36 minutes of recordings
split into a 30-minute training set and a 6-minute test set.

Figure 3: On the left, the markers are placed on two different objects, while
on the right, the markers are placed on a hand. The back of the hand and the
two objects have multiple markers on rigid supports for tracking their
orientation.

###  IV-B Model Training and Testing

Before the training process can begin, the collected data must go through
three processing phases. At first, all the data is transformed into the hand-
palm reference. Based on the first assumption presented in Section III, we
will only focus on the manipulation components related to the fingertips.
Following this assumption, we should consider the hand palm static and only
focus on the motion of the fingertips and objects. In the second phase, the
data is filtered by a low pass filter with a cut-off frequency of 20Hz. The
cut-off frequency is based on a study by Xiong and Quek (2006) [19], which
indicates that a 10Hz sampling rate is sufficient to detect human hand motion.
The third and final step ensures that the dataset $V$ contains only positive
values. This is ensured by offsetting all position data by $0.8$ m and all
object orientation data by $2\pi$.

After the processing phase, we followed the training and segmentation steps as
presented in our previous research [8, 10]. The processed data was segmented
into sequences of 1 second each and stacked to build the columns of the $V$
matrix. We then applied the non-negative matrix factorization method (NMF) to
obtain the $W$ dictionary from the $V$ data matrix. This procedure was
performed separately for the cube and the cylinder, resulting in two distinct
dictionaries of primitives. The reason for using 1-second segments is to model
intermediate configurations of the in-hand manipulation and simplify the
optimization process. The final dictionaries each contain 200 motion
primitives. The processing and training phase took 50 minutes using a single
graphical process unit.

Finally, we tested the dictionary’s ability to recreate human in-hand
manipulation trajectories[8]. We collected errors between the original and
recreated trajectories (see Table I) by considering both the Euclidean
distances of the positions and the difference in object orientation. Both
dictionaries showed good accuracy, with a mean error on fingertips position
equal to $0.6$ mm and a standard deviation of $0.45$ mm. Additionally, the
object orientation and position error had similar results, with mean values of
$0$ rad and $0.4$ mm, and standard deviations of $0.55$ rad and $0.3$ mm,
respectively.

TABLE I: This table shows Euclidean distances between the real and recreated fingers/objects positions, along with object-orientation differences. | Cube Dictionary | Cylinder Dictionary
---|---|---
Fingers
Thumb |  $0.4755\pm 0.4470$ (mm) |  $0.6912\pm 0.6903$ (mm)
Index |  $0.4402\pm 0.4246$ (mm) |  $0.7874\pm 1.0882$ (mm)
Middle |  $0.4276\pm 0.3915$ (mm) |  $0.6846\pm 0.8708$ (mm)
Ring |  $0.4462\pm 0.4398$ (mm) | $0.5591\pm 0.7200\textit{(mm)}$
Little |  $0.4493\pm 0.4271$ (mm) | $0.5511\pm 0.7042\textit{(mm)}$
Object
Translation |  $0.5108\pm 0.5106$ (mm) |  $0.8176\pm 1.0774$ (mm)
Roll |  $0\pm 0.0132$ (rad) |  $0.0001\pm 0.0179$ (rad)
Pitch |  $0.0001\pm 0.0198$ (rad) |  $0\pm 0.0175$ (rad)
Yaw |  $0\pm 0.0017$ (rad) |  $0\pm 0.0039$ (rad)

After confirming that the two dictionaries can accurately recreate human
motion, the following step is integrating them into a pipeline for generating
in-hand manipulation for a robotic platform. To solve the optimization problem
introduced in Section III, we used IBM ILOG’s CPLEX optimization studio in
C++. This software enabled us to tackle complex optimization tasks based on
linear and mixed linear programming. The optimizer requires initial and final
desired stable grasp configurations to obtain fingertip trajectories in the
Cartesian space. Next, an Inverse Kinematics algorithm, based on dynamics-
based recursive linearization (RDBL) principles, converts the fingertip
trajectories from Cartesian to joint space. Finally, the trajectory is used to
control a Shadow Hand. We have created and tested various trajectories
representing different in-hand manipulation scenarios. The experimental setup
and results are discussed in the following section.

##  V Experiments and Results

Our method aims to overcome the limitations of analytical and data-driven in-
hand manipulation approaches. We achieve this by generating a non-complex
method that can cover in-hand manipulation and finger gaiting without
requiring a heavy training process or complex representation of the
constraints.

We tested the ability of our system to achieve the desired object pose by
measuring its position and rotation after each trajectory. Additionally, we
compared the performance of our method to analytical and data-driven
approaches.

###  V-A System overview

We implemented the proposed approach for generating in-hand manipulations
using the ROS middleware, and we tested it on an anthropomorphic robotic hand
(refer to Fig. 4). The experimental setup consists of a Shadow
Hand333https://www.shadowrobot.com/dexterous-hand-series/ connected to a
Shadow Arm and custom-designed fingertips with 6-axis ATI nano17 force and
torque sensors. A state-of-the-art algorithm processes data from the
fingertips’ sensors to estimate the contact of the fingers with the object
based on force and torque sensor measurement [20]. Additionally, the setup
comprised a goniometer and a calliper for measuring the orientations and
translations of the objects.

Figure 4: The picture shows the platform used for the experiments. This
includes the Shadow Hand, a mechanical arm, and 5 force sensors placed on the
fingertips (AT Nano17).

###  V-B Results

The training approach is based on extracting manipulation features unique to
each object. Therefore, for each object, we used the corresponding dictionary.

Since we only considered in-hand manipulations, we fixed the Shadow Hand,
adding stable support, and all object trajectories were executed through
finger actuation. Our experimental scenario was designed to consider the most
common actions among those described by Elliott et al. (1984)[21] in their
classification of in-hand manipulations. Thus, the three operations we
considered were rotation on the $x$ and $y$ axes, and translation on the $y$
axis. These axes correspond to the palm reference system $R_{P}$ in Fig. 1.

After completing each trajectory, we measured the object’s orientation using a
goniometer and its translation with a calliper. The performance of our
algorithm was tested on 21 in-hand manipulation trajectories, evenly
distributed as follows: seven rotations around the $x$ axis, seven rotations
around the $y$ axis, and seven translations along the $y$ axis. For rotation
actions, the objects were rotated within an angle range of 15 to 20 degrees,
while translation actions involved displacements of 5 to 10 cm. Table II shows
the errors associated with each movement, indicating their average, minimum,
and maximum values.

Figure 5: This series of images shows the manipulation of a cube over time.
First, the Shadow Hand rotates the object around the y-axis. Then, the Shadow
Hand translates the object along the y-axis. Finally, the Shadow Hand performs
a composite motion to return the object to its initial pose. (The system of
reference $R_{P}$ is the same considered in Fig. 1).

For the first constraint, we defined the reachable workspace $\mathbb{P}_{i}$
for each finger as the entire workspace of every fingertip. Next, we checked
if the fingertip positions satisfied the reachability constraint described in
Eq. 11. All the points generated with our approach met the reachability
constraint.

For the finger collision constraint, we computed the minimum finger-to-finger
distance as defined in Eq. 13. We detected five times instants with a
collision for the cube and six times instants for the cylinder over the
twenty-one trajectories generated for each object. Therefore instances in
which the collision constraint is not satisfied are rare and, in any case, it
never caused a task failure.

Finally, we recorded the contact data between hand fingertips and the object
using the fingertips sensor addressed in subsection V-A to satisfy the minimum
contact points between the fingertips and the object constraint. For all
generated trajectories, a minimum of two fingers were in contact with the
object at any given time. This is sufficient since the fingertips of the
Shadow Hand provide soft contact.

After testing the ability of the relaxed optimization solver to generate
trajectories that adhere to in-hand manipulation constraints, we evaluated our
method for generating complex tasks. To achieve this, we created trajectories
resulting from a composite of manipulations. First, the Shadow Hand rotated
the object, then translated it, and finally performed a complex motion in
which the object returned to its initial pose, involving a composite rotation
and translation. As shown in Fig. 5, the Shadow Hand successfully returned the
object to its initial pose after the series of three actions. A video
demonstrating the entire procedure is
available444https://www.youtube.com/watch?v=2uszlWyuYWw.

Our approach to the in-hand manipulation task focuses on re-grasping and
finger gaiting from a high-level view. We aim to verify whether the path
planner can execute both re-grasping and finger gaiting actions. By analyzing
recorded data on how the fingertips contact the object, we demonstrate that
the number of contacting fingers varies over time, which confirms finger
gaiting. Additionally, if finger gaiting is detected, any observed changes in
the object’s orientation can be considered strong evidence that re-grasping
has occurred.

In Fig. 6, we present two examples of the variation in the number of contacts
between the object and hand during object rotation and translation. This means
that the object is undergoing in-hand re-grasping since its pose is changing,
and finger gaiting is occurring since the number of fingers in contact with
the object changes during the action.

TABLE II: Performance of the algorithm on rotation and translation motions over the different in-hand manipulation trials. The mean error and the full error range are presented. | Cube | Cylinder
---|---|---
Translation error on y-axis [mm] | $-3.183$ | $0.280$
$[-9,-1]$ | $[+6,-2]$
Rotation error on x-axis [∘] | $-0.573$ | $0.785$
$[-1,+1]$ | $[-1.5,+2.5]$
Rotation error on y-axis [∘] | $-1.089$ | $-0.499$
$[-2,-1]$ | $[-2,+3]$

The median and maximum time for generating a trajectory in Cartesian space and
applying inverse kinematics to find it in joint space is $150$ and $350$ ms,
respectively. We compared our results for the same problem with a classical
approach, as presented in [6], which represents manipulations as a sequence of
finger gating and in-hand re-grasping. In this approach, an optimization
process finds the optimal trajectory for each atomic action. We tested two
different types of solvers for the optimization problem. The first solver aims
to reduce the signed distance between the reachable workspace and the desired
contact points, and it showed a median and maximum planning time of $729.05$
and $3513.96$ s, respectively. The second solver is based on the singular
value decomposition, and it showed a median and maximum planning time of $75$
and $134.275$ s. Therefore, our method can create a trajectory including in-
hand re-grasping and gating faster than this analytical solution alternating
between the two atomic actions.

Figure 6: The picture illustrates the number of fingers in contact with a cube
while executing two manipulations. The left and right graphs display
respectively data for pure rotation and translation. The over-time variation
in the number of contact points demonstrates the method’s ability to generate
finger gaiting.

Finally, it is important to evaluate whether the fast response in trajectory
generation comes at the expense of lengthy procedures for data acquisition and
dictionary training. In our study, collecting human demonstrations and
conducting dictionary training took approximately 50 minutes using a single
GPU. In contrast, using the data-driven strategy described in [22] takes up to
50 hours of training on the experimental set-up using 8 GPUs, roughly
equivalent to 291.6 hours using a single GPU.

##  VI Conclusions

This paper presents a method for generating fingertip trajectories for in-hand
manipulation using motion primitive dictionaries. The approach utilizes human
demonstrations to extract a manipulation dictionary that implicitly respects
in-hand manipulation constraints.

The proposed approach has been tested for in-hand manipulations of two objects
(a cube and a cylinder) and could find trajectories to reach the desired pose
while respecting three constraints that characterize in-hand manipulation
(reachability, finger collision, and minimum contact points). These results
show that the proposed approach retains the constraints from human
demonstrations without requiring a formal representation. The generated
trajectories solve the in-hand manipulation problem by including in-grasp re-
orienting and finger gaiting as actions to move the object to the desired
pose. Furthermore, our approach achieves competitive results in terms of
trajectory generation time and training time compared to analytical and data-
driven approaches, respectively.

In conclusion, we want to mention the possibility of adapting our method to
robotic hands that do not replicate human anatomy. Our approach, centred on
observing the poses of the fingertips and not relying on detailed information
about the kinematic chains of the fingers, suggests that it could be extended
to non-anthropomorphic robotic hands. However, this extension comes with
certain assumptions and limitations. We presuppose that the robotic fingers
possess the necessary dexterity for flexion-extension, adduction and abduction
movements. Furthermore, our method implicitly assumes the use of a five-finger
configuration. This is because the constraints we apply to ensure object
stability are derived from a human motion dictionary, which naturally
presupposes a hand with five fingers. While our method holds promise for
broader application, these considerations highlight the need for careful
adaptation when using robotic hands with configurations that deviate from the
human one.

Given these results, this approach’s trajectories can enable robotic hands to
perform dexterous manipulations. However, future work should extend the
proposed approach to more advanced manipulation tasks, implementing more
precise control strategies and relaxing our initial assumptions.

## ACKNOWLEDGMENT

This work is supported by the CHIST-ERA (2014-2020) project InDex and received
funding from Agence Nationale de la Recherche (ANR) under grant agreement No.
ANR-18-CHR3-0004 and the Italian Ministry of Education and Research (MIUR).

## References

  * [1] A. Carfì, T. Patten, Y. Kuang, A. Hammoud, M. Alameh, E. Maiettini, A. I. Weinberg, D. Faria, F. Mastrogiovanni, G. Alenyà, et al., “Hand-object interaction: From human demonstrations to robot manipulation,” Frontiers in Robotics and AI, vol. 8, p. 316, 2021.
  * [2] E. Theodorou, J. Buchli, and S. Schaal, “A generalized path integral control approach to reinforcement learning,” The Journal of Machine Learning Research, vol. 11, pp. 3137–3181, 2010.
  * [3] G. E. Monahan, “State of the art—a survey of partially observable markov decision processes: theory, models, and algorithms,” Management science, vol. 28, no. 1, pp. 1–16, 1982.
  * [4] U. Prieur, V. Perdereau, and A. Bernardino, “Modeling and planning high-level in-hand manipulation actions from human knowledge and active learning from demonstration,” in IEEE/RSJ International Conference on intelligent Robots and Systems, (Vilamoura-Algarve, Portugal), December 2012.
  * [5] B. Sundaralingam and T. Hermans, “Relaxed-rigidity constraints: kinematic trajectory optimization and collision avoidance for in-grasp manipulation,” Autonomous Robots, vol. 43, no. 2, pp. 469–483, 2019.
  * [6] B. Sundaralingam and T. Hermans, “Geometric in-hand regrasp planning: Alternating optimization of finger gaits and in-grasp manipulation,” in IEEE International Conference on Robotics and Automation (ICRA), (Brisbane Convention & Exhibition Centre, Brisbane, Australia), May 2018.
  * [7] Y. Fan, W. Gao, W. Chen, and M. Tomizuka, “Real-time finger gaits planning for dexterous manipulation,” IFAC-PapersOnLine, vol. 50, no. 1, pp. 12765–12772, 2017.
  * [8] A. Hammoud, A. Diouf, and V. Perdereau, “A robotic in-hand manipulation dictionary based on human data,” in 20th International Conference on Advanced Robotics (ICAR), (Ljubljana, Slovenia), December 2021.
  * [9] A. J. Ijspeert, J. Nakanishi, and S. Schaal, “Learning attractor landscapes for learning motor primitives,” in 15th International Conference on Neural Information Processing Systems (NIPS), (Vancouver, British Columbia, Canada), January 2002.
  * [10] A. Hammoud, V. Belcamino, A. Carfi, V. Perdereau, and F. Mastrogiovanni, “In-hand manipulation planning using human motion dictionary,” in 2022 31st IEEE International Conference on Robot and Human Interactive Communication (RO-MAN), pp. 927–933, 2022.
  * [11] E. Cao, K. Cao, K. Feng, and J. Wang, “Nmf based image sequence analysis and its application in gait recognition,” CCF Transactions on Pervasive Computing and Interaction, vol. 2, no. 2, pp. 86–96, 2020.
  * [12] Z. Akata, C. Thurau, and C. Bauckhage, “Non-negative matrix factorization in multimodality data for segmentation and label prediction,” in 16th Computer vision winter workshop, (Mitterberg, Austria), February 2011.
  * [13] O. Berne, C. Joblin, Y. Deville, J. Smith, M. Rapacioli, J. Bernard, J. Thomas, W. Reach, and A. Abergel, “Analysis of the emission of very small dust particles from spitzer spectro-imagery data using blind signal separation methods,” Astronomy & Astrophysics, vol. 469, no. 2, pp. 575–586, 2007.
  * [14] C. Vollmer, S. Hellbach, J. Eggert, and H.-M. Gross, “Sparse coding of human motion trajectories with non-negative matrix factorization,” Neurocomputing, vol. 124, pp. 22–32, 2014.
  * [15] A. Vignolo, N. Noceti, A. Sciutti, F. Odone, and G. Sandini, “Learning dictionaries of kinematic primitives for action classification,” in 2020 25th International Conference on Pattern Recognition (ICPR), pp. 5965–5972, 2021.
  * [16] J. R. Flanagan, M. C. Bowman, and R. S. Johansson, “Control strategies in object manipulation tasks,” Current Opinion in Neurobiology, vol. 16, pp. 650–659, 2006.
  * [17] S. Sra and I. Dhillon, “Generalized nonnegative matrix approximations with bregman divergences,” Advances in neural information processing systems, vol. 18, p. 283–290, 2005.
  * [18] D. Prattichizzo and J. C. Trinkle, Grasping, pp. 671–700.  Berlin, Heidelberg: Springer Berlin Heidelberg, 2008.
  * [19] Y. Xiong and F. Quek, “Hand motion gesture frequency properties and multimodal discourse analysis,” International Journal of Computer Vision, vol. 69, no. 3, pp. 353–371, 2006.
  * [20] H. Liu, X. Song, J. Bimbo, L. Seneviratne, and K. Althoefer, “Surface material recognition through haptic exploration using an intelligent contact sensing finger,” in 2012 IEEE/RSJ international conference on intelligent robots and systems, pp. 52–57, IEEE, 2012.
  * [21] J. M. Elliott and K. Connolly, “A classification of manipulative hand movements,” Developmental Medicine & Child Neurology, vol. 26, no. 3, pp. 283–296, 1984.
  * [22] O. M. Andrychowicz, B. Baker, M. Chociej, R. Jozefowicz, B. McGrew, J. Pachocki, A. Petron, M. Plappert, G. Powell, A. Ray, et al., “Learning dexterous in-hand manipulation,” The International Journal of Robotics Research, vol. 39, no. 1, pp. 3–20, 2020.


2022 xx Month 2022

# Applications of Gaussian Mutation for Self Adaptation in Evolutionary
Genetic Algorithms

Okezue Bell1,2  1Harvard Medical School-MIT<EMAIL_ADDRESS>2Moravian
Academy Upper School

(4 January 2022)

###### Abstract

In recent years, optimization problems have become increasingly more prevalent
due to the need for more powerful computational methods. With the more recent
advent of technology such as artificial intelligence, new metaheuristics are
needed that enhance the capabilities of classical algorithms. More recently,
researchers have been looking at Charles Darwin’s theory of natural selection
and evolution as a means of enhancing current approaches using machine
learning. In 1960, the first genetic algorithm was developed by John H.
Holland and his students (Holland, 1975). We explore the mathematical
intuition and implications of the genetic algorithm in developing systems
capable of evolving using Gaussian mutation.

###### keywords:

EvolutionGaussian MutationGenetic Algorithm

††journal: JMLFS-ID

##  1 Gaussian Mutations in Genetic Algorithms

  * •

Introduction to genetic algorithms

  * •

Self-adaptation

  * •

The Gaussian mutation operator

  * •

Implementation of gaussian mutation results

  * •

Conclusion

##  2 Introduction to genetic algorithms

A genetic algorithm is a type of stochastic search algorithm that functions
off of the fundamental principles of natural selection. Essentially, the
algorithm works off of given string structures, typically organisms, which
gradually evolve, and their fitness is characterized by survival and
reproduction. The genetic algorithm sustains the survival of the fittest rule,
where there is a randomized information exchange but with a structured method
(Gholizadeh, 2013). With every succeeding generation, new strings are created
and the fittest strings from the previous generations are used as a reference.
The ”winners” are considered the individuals whose genes enable them to
survive and reproduce. With each successive selection, more resources are
retained. The algorithm is capable of simulating this process and uses it to
calculate the optimum of objective functions, thus making the algorithm
suitable for solving optimization problems. In summary, the genetic algorithm
(GA) uses neuro-evolutionary concepts and artificial intelligence to design
better neural networks (Roetzel et al., 2020).

When most of the algorithm’s intricacies are removed, it essentially takes in
a set of something and evolves it to its optimal state. Figure 1 shows the
abbreviated process for a cluster of carbon atoms being evolved to organize in
the lowest energy carbon structure (Figure 1).

Figure 1: The GAs process of determining the buckminsterfullerene as the
lowest energy configuration of 60 carbon atoms. It has a carbon population and
knows it wants to solve for the lowest energy structure. It then generates a
bunch of random configurations, which are then analyzed for their fitness
based on their energy levels. The best ones mate and create a new generation,
and this process repeats until an ultralow energy level generation is created
(the buckminsterfullerene structure). It took scientists over 3 years to solve
this, while the GA can complete this process in about 4 hours (self-made).

GAs have been applied to a variety of fields, including inverse kinematic
problems (IKPs) for advanced robotic designs for transhumeral prosthetics
(Števo et al., 2019) and tumor tracking in combination with intermittent
Kalman predictors (IKPs) by using genetic algorithms to solve combinatorial
optimization problems to conclude the measurements of tumor motion (Aspeel et
al., 2019). Currently, the popular usage of GAs falls in the category of
generating an ANN capable of determining a very precise solution for an
infinite set of possibilities, more specifically a search or optimization
problem. Genetic algorithms are considered a subset of artificial intelligence
and are a part of the broader scope of evolutionary computing (Eiben, 2003).

##  3  Self-adaptation

What is currently considered one of the most exciting fields of research for
evolutionary computing and genetic algorithms (Hinterding et al., 2005) is the
concept of self-adaptation in GAs, especially for numeric functions. Self-
adaptation is the property of a genetic algorithm to adapt its algorithm to a
specific problem when calculating the solution for it. This can be done for
one or more facets of evolution, and approaches for self-adapting for both
mutation strength and population size, for example, have been successful on
initial testing (Hinterding et al., 2005). The mutation genetic operator
introduces genetic diversity with the genetic algorithm population of
chromosomes. It can be likened to biological mutation. There are many types of
mutation operators, including bit string, boundary, and Gaussian, most of
which deal with the introduction of a random value to a specific gene, which
influences the expression of the offspring.

The Gaussian mutation operator has proven to be an optimal and popular choice
for self-adaptation in genetic algorithms (Li-hui Zhang Ph.D. et all., 2015).
In Gaussian mutation, a random unit Gaussian distribution value is added to a
selected gene by the operator. The value is thus added to each element of the
gene carrier’s vector, which causes the development of a new offspring (Figure
2).

Figure 2: Gaussian Mutation of the a) parent to form the b) offspring, causing
variance (Dorronsoro, 2004).

Aside from Gaussian mutation, there are other popular operators, including the
one-point crossover reproduction operator (Figure 3), and the bit-flipping
mutation operator (Figure 4).

Figure 3: One point crossover of parents a) and b) to form the c) and d)
offspring from reproduction. A subsequence of the parents, represented as
strings, is swapped to create two new offspring (Dorronsoro, 2005).

Figure 4: The bit flipping mutation of the a) parent to form the b) offspring.
A singular binary bit in the original parent string is flipped to create a new
string for a different offspring (Dorronsoro, 2005).

However, Gaussian mutation operators are still widely used in genetic
algorithms, especially in a variety of object-oriented programming languages
and their use of Gaussian distribution values. This is because of the Gaussian
curve’s (Figure 5) usage in probability density functions and standard
deviation. Usage of Gaussian values is promising for mutation simulations, and
when introducing new information to the population of the GA.

Figure 5: A selection of Normal Distribution Probability Density Functions
(PDFs). Both the mean, $\mu$, and variance, $\sigma^{2}$, are varied. The red
curve on the graph represents the standard normal distribution. The key is
given on the graph (Inductive Load, 2008).

##  4 The Gaussian mutation operator

Take $x\in[a,b]$ as a real variable. The Gaussian mutation operator $M_{g}$
changes the variable $x$ to

$M_{g}(x):=min(max(N(x,\sigma),a),b)$

The operator also leverages the Gauss error function. Recall that all
solutions in the genetic algorithms are the chromosomes. These chromosomes go
through mutation (random gene alterations), crossover (swapping genetic
material), and selection (only the fittest survive). Mutation operators use a
single parent chromosome and induce mutation on some of the selected genes.
The offspring when using the Gaussian mutation operator is given by

$x^{{}^{\prime}}_{i}=\sqrt{2}\sigma(b_{i}-a_{i})erf^{-1}(u^{{}^{\prime}}i)$

, where $x^{{}^{\prime}}_{i}$ is the offspring, $\sigma$ is a fixed parameter
for all variables, and $[a_{i},b_{i}]$ represents the range of random gene
$x_{i}$. $erf^{-1}$ represents the inverse of $erf$, which is the Gauss error
function, defined by $erf(y)=\frac{2}{\sqrt{\pi}}\int_{0}^{y}e^{-t^{2}}\,dt$.
The inverse Gauss error function is applied to a calculation ui,, which is a
random value $u_{i}$ within the range of $(0,1)$, which is inside of a formula
that gives $u_{i}$, (”Mutation Algorithms for Real-Valued Parameters (GA)”,
2018). This formula can be defined using IFTTT (if this then that) notation.

With evolutionary programming, the manipulation of parameters is done through
self-adaptation. This learning format allows the GA to constantly tweak the
strategy parameters when it’s iterating through different generations to find
the optimum, essentially searching for the value of the parameters themselves
to increase their performance levels. Parametric optimization values are
subject to change; their fluctuations depend on the solved function. With
Gaussian mutation, convergence is completed much more efficiently, and domain
searching and solution tuning are much better (Hartfield, 2010), thereby
increasing the performance of the GA.

The Gaussian mutation’s shape can also be controlled by applying a parameter
to it. This mutation operator is called the q-Gaussian mutation, where the
shape of the mutation distribution is controlled via a real parameter q (Tinós
et all., 2011). The q-Gaussian, therefore, operates on a different
distribution than the traditional Gaussian mutation. This is known as the
q-analog of the normal curve, and its distribution is symmetrical at about 0.
Figure 6 shows its probability density function.

Figure 6: Probability density plots of q-Gaussian distributions. This maps the
probability density function of the q-Gaussian distribution itself
(IkamusumeFan, 2014).

q-Gaussians allow the mutation distribution shape of the GA self-adaptive,
which makes them effective in solving dynamic optimization problems (Yang et
al., 2010), where the variable values are subject to change over time.

##  5 Implementation of gaussian mutation results

The implementation of Gaussian mutation has been successful in solving a
variety of optimization problem types. A notable instance of this is in An
Improved Real-Coded Genetic Algorithm Using the Heuristical Normal
Distribution and Direction-Based Crossover, where a multi-offspring improved
real-coded genetic algorithm (MOIRCGA) using the Gaussian distribution (in
combination with direction-based crossover was proven to solve constrained
optimization problems. The MOIRCGA finished with supremacy over a real-coded
genetic algorithm (RCGA) not leveraging Gaussian distribution to find a
globally optimal solution to the problem (Wang et al., 2019), in part due to
its much faster conversion speed than the RCGA (Figure 7).

Figure 7: Optimization results for the MOIRCGA as compared to the RCGA in the
An Improved Real-Coded Genetic Algorithm Using the Heuristical Normal
Distribution and Direction-Based Crossover Experiment (Wang et al., 2019). The
MOIRCGA was superior to the RCGA and some of the other benchmarked algorithms.

Within this experiment, the algorithms were tested and compared with other
literature, and the problem was used in the parameter optimization of the
Cantilever Beam design with a discrete cross-sectional area to determine the
validity of the MOIRGCA, and sixteen test functions were evaluated as well.
The experimentation concludes with, ”optimization results show that the
function value obtained with MOIRCGA is superior to that obtained with RCGA”
(Wang et al., 2019).

Within this experiment, the algorithms were tested and compared with other
literature, and the problem was used in the parameter optimization of the
Cantilever Beam design with a discrete cross-sectional area to determine the
validity of the MOIRGCA, and sixteen test functions were evaluated as well.
The experimentation concludes with, ”optimization results show that the
function value obtained with MOIRCGA is superior to that obtained with RCGA”
(Wang et al., 2019).

###  5.1 Robotic appendages.

Genetic algorithms have been cited as one of the most useful forms of
artificial intelligence (Cheng, 2011). The algorithms have extensive
implications in robotics, especially within movement or functionality
(Davidor, 1991). The robotics inverse kinematics problem (IKP) describes
”finding a vector of joint variables which produce [the] desired end-effector
location” (DeMers et al., 1997). The problem has an infinite amount of
solutions, making it difficult to solve using classical computation. However,
genetic algorithms have been presented to create the best robotic arm (Sekaj
et al., 2014) based on trajectory control. In Optimization of Robotic Arm
Trajectory Using Genetic Algorithm, a genetic algorithm that could optimize
for energy consumption, operating time, rotation changes, and trajectory was
proposed for the industrial robot ABB IRB 6400FHD (Števo et all., 2014). In
Genetic Algorithm Based Approach for Autonomous Mobile Robot Path Planning, a
GA with an improved crossover operator was suggested to be applied to a path
planning problem for a computer vision autonomous robot for movement and
navigation.

###  5.2 Molecular geometry.

In the introductory section of this independent study, an example of
determining carbon fullerenes as the lowest energy configuration of sixty
carbon atom clusters. This example for a genetic algorithm was extracted from
Molecular Geometry Optimization with a Genetic Algorithm, in which a genetic
algorithm was used to determine the lowest energy geometry for carbon
compounds, ranging up to $C_{6}0$. The experiment, however, wasn’t limited, as
their method was capable of determining the lowest energy structure of an
atomic cluster in any arbitrary model potential. Original solution candidates
generated were used as a reference for fitness to generate a new generation of
solutions to produce lower energy solutions. The genetic model itself was
cited to be more accurate than the simulated annealing, and was stated to have
”dramatically outperformed” it, and decreased the computational time of
concluding carbon fullerenes by over 200$\%$ when compared to a human
biochemist (Deaven et al., 1995). In First-Principles Molecular Structure
Search with a Genetic Algorithm, a genetic algorithm was built using Conda’s
RDKit library in Python to accelerate computational chemistry. The algorithm
focused on identifying low-energy conformers for a given molecule using first
principle searches (Supady et al., 2015). The algorithm was successful (Figure
8) in doing so.

Figure 8: ”Evaluation of the results for the subset of the Astex Diverse Set.
(A) The relative energy of all found conformers as a function of the RMSD to
the reference ligand (blue circles). The green squares depict the reference
ligand structures after DFT optimization. (B) An overlay between the reference
ligand (green) and the best match (blue) is presented together with the
corresponding RMSD value” (Supady et al., 2015).

###  5.3 Sentiment Lexicon Optimization.

Another use of genetic algorithms is in analyzing emotions, also known as
opinion mining or sentiment analysis. In ALGA: Adaptive lexicon learning using
genetic algorithm for sentiment analysis of microblogs, the common problem in
sentiment analysis, ”improving polarity classification of sentiments in
microblogs by building adaptive sentiment lexicons” (Keshvarz et al., 2017),
was solved using a novel genetic algorithm. In the six datasets through which
it was tested, the algorithm achieved over 80$\%$ accuracy. In summary, the
algorithm was able to build the most accurate collection of words, and how
their semantics relate to their sentiment orientation (Keshvarz et al., 2017).

##  6 Conclusion

## Acknowledgements

This research was conducted independently and presented at a UPenn SEAS
functional computing paper competition. Author Okezue Bell is affiliated with
both Harvard-MIT and Moravian Academy, though this research was conducted
independently.

## References

  * [1] Keshavarz H, Abadeh MS (2017) ALGA: Adaptive lexicon learning using genetic algorithm for sentiment analysis of microblogs. Retrieved December 19, 2021, from https://doi.org/10.1016/j.knosys.2017.01.028.
  * [2] Deaven DM, Ho KM (1995) Molecular Geometry Optimization with a Genetic Algorithm. In: Physical Review Letters. Retrieved December 19, 2021, from https://doi.org/10.1103/PhysRevLett.75.288.
  * [3] Lamini C, Benhlima S, Elbekri A (2018) Genetic Algorithm Based Approach for Autonomous Mobile Robot Path Planning. In: Procedia Computer Science. Retrieved December 19, 2021, from https://doi.org/10.1016/j.procs.2018.01.113.
  * [4] Inverse Kinematics Problem. In: Inverse Kinematics Problem - an overview — ScienceDirect Topics. Retrieved December 19, 2021, from https://www.sciencedirect.com/topics/engineering/inverse-kinematics-problem.
  * [5] Števo S, Sekaj I, Dekan M (2016) Optimization of Robotic Arm Trajectory Using Genetic Algorithm. In: IFAC Proceedings Volumes. Retrieved December 19, 2021, from https://doi.org/10.3182/20140824-6-ZA-1003.01073.
  * [6] Genetic Algorithms and Robotics. In: World Scientific. Retrieved December 19, 2021, from https://doi.org/10.1142/1111.
  * [7] VII Jornadas Zaragoza-Pau de Matemática Aplicada y estadística : Jaca (Huesca), 17-18 de septiembre de 2001, 2003-01-01, ISBN 84-96214-04-4, pags. 27-36.
  * [8] CSDL: IEEE Computer Society. In: CSDL — IEEE Computer Society. Retrieved December 19, 2021, from https://doi.org/10.1109/SBRN.2010.46.
  * [9] R. Tins and S. Yang, ”Evolution Strategies with q-Gaussian Mutation for Dynamic Optimization Problems,” in Neural Networks, Brazilian Symposium on, Sao Bernardo do Campo, Sao Paulo Brazil, 2010 pp. 223-228. Retrieved December 19, 2021, from https://doi.org/10.1109/SBRN.2010.46.
  * [10] Pontes EAS (2018) A Brief Historical Overview Of the Gaussian Curve: From Abraham De Moivre to Johann Carl Friedrich Gauss. In: ijesl.org. Retrieved December 19, 2021, from http://www.ijesi.org/papers/Vol(7)i6/Version-5/D0706052834.pdf.
  * [11] Hinterding R, Michalewicz Z, Eiben AE (2005) Adaptation in Evolutionary Computation: A Survey. Retrieved December 19, 2021, from https://cs.adelaide.edu.au/ zbyszek/Papers/p34.pdf.
  * [12] Aspeel A, Dasnoy D, Jungers Raphael M, Macq B (2019) Optimal Intermittent Measurements for Tumor Tracking in X-ray Guided Radiotherapy. Retrieved December 19, 2021, from https://arxiv.org/pdf/1903.08990.pdf.


<EMAIL_ADDRESS>††thanks: NHFP Einstein fellow<EMAIL_ADDRESS><EMAIL_ADDRESS>
# The relative contribution to heavy metals production from binary neutron
star mergers and neutron star-black hole mergers

Hsin-Yu Chen  [ Salvatore Vitale  [ LIGO Laboratory, Massachusetts Institute
of Technology, 185 Albany St, Cambridge, MA 02139, USA  Department of Physics
and Kavli Institute for Astrophysics and Space Research, Massachusetts
Institute of Technology, 77 Massachusetts Ave, Cambridge, MA 02139, USA
Francois Foucart  [ Department of Physics & Astronomy, University of New
Hampshire, 9 Library Way, Durham NH 03824, USA

###### Abstract

The origin of the heavy elements in the Universe is not fully determined.
Neutron star-black hole (NSBH) and binary neutron star (BNS) mergers may both
produce heavy elements via rapid neutron-capture (r-process). We use the
recent detection of gravitational waves from NSBHs, improved measurements of
the neutron star equation-of-state, and the most modern numerical simulations
of ejected material from binary collisions to measure the relative
contribution of NSBHs and BNSs to the production of heavy elements. As the
amount of r-process ejecta depends on the mass and spin distribution of the
compact objects, as well as on the equation-of-state of the neutron stars, we
consider various models for these quantities, informed by gravitational-wave
and pulsar data. We find that in most scenarios, BNSs have produced more
r-process elements than NSBHs over the past 2.5 billion years. If black holes
have preferentially small spins, BNSs can produce at least twice of the amount
of r-process elements than NSBHs. If black hole spins are small and there is a
dearth of low mass ($<5M_{\odot}$) black holes within NSBH binaries, BNSs can
account for the near totality of the r-process elements from binaries. For
NSBH to produce large fraction of r-process elements, black holes in NSBHs
must have small masses and large aligned spins, which is disfavored by current
data.

##  1\. Introduction

Neutron star-black hole (NSBH) and binary neutron star (BNS) mergers are among
possible astrophysical formation sites of heavy elements through rapid
neutron-capture (r-process) nucleosynthesis (Lattimer & Schramm, 1974). The
tidal disruption of a neutron star (NS) by its black hole (BH) companion, and
the subsequent ejection of neutron-rich material into the interstellar medium
was originally proposed by Lattimer & Schramm (1974, 1976). Their predictions
have largely been confirmed by modern numerical simulations: as long as the BH
is of sufficiently low mass for tidal disruption to occur,
$\sim(0.01-0.1)M_{\odot}$ of neutron-rich material can be ejected during a
NSBH merger (Foucart et al., 2013; Kyutoku et al., 2015). BNS mergers
typically do not lead to as much mass ejection ($\lesssim 0.01M_{\odot}$;
Hotokezaka et al. (2013); Dietrich & Ujevic (2017)), unless their mass ratios
are highly unequal (possibly $\lesssim 0.03M_{\odot}$ for systems with NS
masses $\gtrsim 1.1M_{\odot}$; Dietrich et al. (2017); Kiuchi et al. (2019)).
On the other hand, disrupting NSBH systems and most BNS systems can both
produce compact remnants surrounded by massive accretion disks, with a
significant fraction of these disks expected to become unbound within a few
seconds of the merger (Fernández & Metzger, 2013; Siegel & Metzger, 2017;
Christie et al., 2019). Roughly $\sim 0.01M_{\odot}$ of mass can be ejected
during the post-merger evolution.

The relative importance of BNSs, NSBHs, and other potential sources of
r-process elements (e.g. collapsars; Surman et al. (2006); Siegel et al.
(2019)), whose potential as sources of r-process elements is still under
debate (Miller et al., 2020; Fujibayashi et al., 2020), and magnetorotational
core-collapse supernovae( Winteler et al. (2012); Mösta et al. (2018); Yong et
al. (2021), which however cannot easily produce the heavier r-process
elements) remains highly uncertain. To determine the relative contribution of
BNSs and NSBHs to the production of r-process elements one needs to quantify
a) how much r-process material is ejected by each system as a function of its
parameters and b) the merger rate of BNS and NSBH as a function of the
system’s parameters. Predictions might be compared to measurements of
r-process abundances on Earth (Paul et al., 2001; Wallner et al., 2021), in
the solar system (Meyer, 1993; Côté et al., 2021), and in stars other than the
Sun (Roederer & Lawler, 2012; Ji et al., 2016; Frebel, 2019; Holmbeck et al.,
2020).

The observation of a kilonova following the first gravitational wave (GW)
detection of a BNS, GW170817 (Abbott et al., 2017, 2017; Coulter et al.,
2017), is consistent with the models for an emission powered by radioactive
decays of heavy elements produced through r-process nucleosynthesis (see e.g.
Metzger (2020) for a review) with $\lesssim 0.05M_{\odot}$ of ejected matter.
Significant uncertainties exist arising from the details of the nuclear
physics processes, as well as the composition and complex 3D geometry of the
outflows. The observations of ancient dwarf galaxies also favor r-process
enrichment from rare events, such as BNSs, producing copious amounts of
r-process material (Ji et al., 2016). It thus seems very likely that some
r-process nuclei are produced in BNSs. Indeed, using the astrophysical rate of
BNS mergers inferred from the GW observations, one can estimate the BNS
contribution to the production of r-process elements (Metzger et al., 2010;
Abbott et al., 2017). On the other hand, the recent discovery of two NSBHs by
the LIGO-Virgo-Kagra (LVK) collaboration, GW200105 and GW200115 (Abbott et
al., 2021a), has provided the first direct evidence of this type of systems.
The merger rate and intrinsic properties inferred from these new sources
enable novel constraints on the relative contribution to the production of
r-process elements from BNSs and NSBHs.

##  2\. Methods

The amount of mass ejecta from binary mergers, and thus the r-process yield,
depends on the mass and spin distribution of compact objects, as well as on
the equation-of-state (EoS) of NSs. As none of these are exactly known, we
consider several possible models, which are consistent with GW and pulsar
measurements. For each model, we generate populations of BNSs and NSBHs, use
fits to numerical simulations to assess the amount of mass ejecta from each
binary, and calculate the total contribution from BNSs and NSBHs. Each step of
the analysis is described in turn here below.

###  2.1. Simulated populations

We consider several distributions of mass and spin for BNSs and NSBHs. For the
NSs (in both BNSs and NSBHs) we use the mass distribution of Farr &
Chatziioannou (2020), in which the observed Galactic NS masses were modelled
as a bimodal distribution (See also Alsing et al. (2018)). We further restrict
the NS mass to the range $1M_{\odot}<m_{\rm NS}<m_{\rm TOV}$, where $m_{\rm
TOV}$ is the EoS dependent maximum mass of a cold and non-rotating NS. Recent
work has also shown that a mixture of two Gaussians fits well the NSs detected
in LVK data (Landry & Read, 2021). For the BHs, we consider three different
mass distributions: a) uniform-in-log in the mass range $[5,40]M_{\odot}$, b)
uniform-in-log in the mass range $[m_{\rm TOV},40]M_{\odot}$, and c) the Power
Law +Peak with GW190814 mass distribution from Abbott et al. (2021c). The
models a) and b) consider a generic BH mass function with and without the
observational mass gap between NSs and BHs (Kreidberg et al., 2012). The
scenario c) adopts the primary BH mass distribution of LVK binary black hole
mergers (BBHs) assuming the secondary component mass of GW190814 is a BH
(Abbott et al., 2021c). Not only does this model excludes the existence of a
mass gap (Kreidberg et al., 2012), but it also allows for more BHs with masses
in the range $\sim 2-3M_{\odot}$, potentially generated by BNS mergers (Lu et
al., 2021; Liu & Lai, 2021).

We treat the NS as non-spinning. For the BH spins, we consider two
possibilities: a) dimensionless spin magnitudes uniformly distributed in the
range $[0,0.95]$ and aligned with the binary total angular momentum; b) spin
magnitudes and tilt angles that follow the primary BH spin distribution
reported by the LVK in Abbott et al. (2021c).

In Table 1 we list all six mass and spin models. For each model, we generated
100,000 simulated BNSs and NSBHs and calculate the ejected mass as described
below.

Label | $m_{1}$ | $|\chi_{1}|$ | Tilt | $M_{\rm{ej,NSBH}}/M_{\rm{ej,Total}}$
---|---|---|---|---
Gap+aligned spin | Uniform in log, $[5,40]M_{\odot}$ | Uniform in [0,0.95] | Aligned | 30%
Gap+BBH-like spin | Uniform in log, $[5,40]M_{\odot}$ | BBH-like | BBH-like | 1%
No gap+aligned spin | Uniform in log, $[m_{\rm TOV},40]M_{\odot}$ | Uniform in [0,0.95] | Aligned | 49%
No gap+aligned spin | Uniform in log, $[m_{\rm TOV},40]M_{\odot}$ | BBH-like | BBH-like | 11%
BBH-like mass+aligned spin | BBH-like | Uniform in [0,0.95] | Aligned | 77%
BBH-like mass+spin | BBH-like | BBH-like | BBH-like | 35%
Table 1Summary of BH mass and spin models explored in this paper. For the NSs
in both BNSs and NSBHs, we use a bimodal distribution fitted to the Galactic
NS population (Alsing et al., 2018; Farr & Chatziioannou, 2020). The label
BBH-like represents the primary BH distribution inferred from the LVK binary
black hole merger (BBH) observations (Abbott et al., 2021c), including
GW190814. We stress that this model predicts more BHs with masses in the range
$\sim 2-3M_{\odot}$ and small spin-aligned magnitudes. The last column reports
the highest possible NSBH mass ejecta fraction given the 90% upper limit of
the NSBH/BNS astrophysical rate ratio.

###  2.2. Estimation of ejected mass

BNSs and NSBHs eject mass via the tidal disruption of NSs, the disk outflows
in the post-merger remnant phase, and, for BNSs, during the collision of the
two NSs. The post-merger outflows themselves can be further subdivided into
early outflows from spiral arms in the remnant (Nedora et al., 2021),
magnetically-driven winds (Siegel & Metzger, 2017; Christie et al., 2019),
neutrino-driven winds (Just et al., 2015), and thermal outflows in the
advection-dominated disk formed late in the evolution of the remnant
(Fernández & Metzger, 2013). Each outflow component may have different
composition, temperature, and velocity, impacting both the outcome of
r-process nucleosynthesis (Wanajo et al., 2014; Lippuner & Roberts, 2015) and
the properties of the associated kilonova (Kasen et al., 2013). While a
reliable model of all of the outflow components is not currently available,
analytical fits to numerical simulations of BNSs and NSBHs exist. These
provide estimates for the total amount of mass of the dynamical ejecta (mass
ejected during the first few milliseconds following the merger) and of post-
merger disks, as well as for the fraction of the disk that will be unbound
after the merger. As we focus on the total mass of matter unbound by merger
events, these estimates will be sufficient for our purpose.

We model the total mass of ejecta from a given binary as:

$m_{\rm{ej}}=\alpha_{\rm{dyn}}m_{\rm{dyn}}+f_{\rm{loss}}m_{\rm{disk}}$ (1)

where $m_{\rm{dyn}}$ represents the mass of the dynamical ejecta,
$m_{\rm{disk}}$ the mass of the disk formed in the post-merger phase, and
$f_{\rm{loss}}$ the fraction of mass ejected from the disk. We also introduce
a scaling factor $\alpha_{\rm{dyn}}$ that will be varied to account for
uncertainties in the knowledge of $m_{\rm{dyn}}$.

In order to estimate $m_{\rm{dyn}}$ and $m_{\rm{disk}}$, we use analytical
fits to numerical simulations. For BNSs, we use Equation 6 of Krüger & Foucart
(2020) for the dynamical ejecta and Equation 4 of Krüger & Foucart (2020) for
the disk. For NSBHs, we use Equation 9 of Krüger & Foucart (2020) for the
dynamical ejecta and Equation 4 of Foucart et al. (2018) for the disk 111We
note that the original formula for the disk only takes into account the BH
spin magnitude (not its orientation) due to the spin-aligned numerical
results. Therefore, we replace the spin magnitude with the component of the
spin along the orbital angular momentum when dealing with tilted BH spins.
When the baryonic mass of NSs is needed in these formulae, we estimate it from
the gravitational mass of the NS using Equation 33 of Lattimer & Prakash
(2001).

The analytical formulae for NSBH binaries are valid for the range of EoS used
in this study, as long as the aligned component of the dimensionless BH spin
is $\chi_{\rm BH}\lesssim 0.9$. In these cases, we expect $\sim(10-20)\%$
accuracy (Krüger & Foucart, 2020; Foucart et al., 2018). On the other hand,
there are no simulations for very large mass ratios ($M_{\rm BH}/M_{\rm
NS}\gtrsim 8$) or very compact stars. However, this is not a serious
limitation since it’s known that the NS is not disrupted in these cases
(except for extreme BH spin magnitudes). The analytical formulae we use
correctly predict a lack of mass ejection for large mass ratios. For BNS
mergers, uncertainties are typically larger. Even in the regions of parameter
space where numerical simulations are available, $(30-50)\%$ errors in the
fitting formulae are to be expected (Krüger & Foucart, 2020). As a result, we
choose a range of $\alpha_{\rm dyn}$ that covers the expected uncertainty in
$m_{\rm dyn}$, i.e. $\alpha_{\rm dyn,BNS}\in[0.5,1.5]$ for BNSs and
$\alpha_{\rm dyn,NSBH}\in[0.8,1.2]$ for NSBHs. Additionally, not enough
simulations are available for asymmetric BNSs ($m_{1}/m_{2}\gtrsim 1.3$) to
reliably predict the amount of mass ejection in that regime. No numerical
simulation involving NSs of mass $\geq 1.1M_{\odot}$ has found dynamical mass
ejection $>0.03M_{\odot}$, while extrapolation of the fitting formulae towards
asymmetric binaries would lead to much more ejecta ($\gtrsim 0.06M_{\odot}$).
To avoid issues with a small number of asymmetric binaries ($<15\%$ in our
simulations) creating unphysically large amounts of r-process elements, we
impose a cap $m_{\rm cap}=0.03M_{\odot}$ on the dynamical ejecta for BNSs.

We vary $f_{\rm loss}$ to account for the uncertainty in the physics of post-
merger remnant and in the value of $m_{\rm disk}$. Recent studies of BH-disk
systems indicate that $(5-20)\%$ of the mass of the disk will be ejected in
magnetically-driven winds (Siegel et al., 2019; Christie et al., 2019), a
comparable amount will be ejected through thermal outflows on longer
timescales (Fernández et al., 2020), and additional mass ejection is expected
during the circularization of the accretion disk (Kiuchi et al., 2015). In NS-
disk systems, additional mass ejection is possible due to spiral-arm
instabilities in the disk (Nedora et al., 2021), and most of the disk may be
unbound in the presence of a long-lived NS remnant (Metzger & Fernández, 2014;
Fahlman & Fernández, 2018). Therefore we expect higher $f_{\rm loss,BNS}$ if a
long-lived NS-disk system is formed after the BNS merger instead of a promptly
collapsed BH-disk system (Metzger & Fernández, 2014). On the other hand, even
if the BNS collapses into a BH shortly after the merger, the resulting BH-disk
systems would have lower mass than the systems formed from NSBH mergers and
are expected to yield higher $f_{\rm loss}$ (Fernández et al., 2020).
Accordingly, we take $f_{\rm loss,NSBH}\in[0.15,0.60]$ for NSBH systems (which
always form BH-disk systems), $f_{\rm loss,BNS}\in[0.15,1.0]$ for BNS systems
(which may form BH-disk or NS-disk systems depending on when/if the remnant
collapses to a BH), and we assume $f_{\rm loss,BNS}\geq f_{\rm loss,NSBH}$.

###  2.3. Choice of neutron star equation-of-state

Figure 1.— Neutron star EoS provided by the Xtreme Astrophysics Group (Özel et
al., 2016; Bogdanov et al., 2016; Özel & Freire, 2016) with mass-radius
relations consistent with the combined measurement presented in Raaijmakers et
al. (2020).

The amount of ejecta is highly dependent on the compactness of NSs ($C_{\rm
NS}\equiv GM/Rc^{2}$, where $M$ and $R$ are the NS mass and radius,
respectively), which is an EoS dependent characteristic. The EoS also affects
some of our models setting the minimum mass of the BH mass function (Table 1).

EoS measurements based on pulsars (Cromartie et al., 2020; Riley et al., 2019;
Miller et al., 2019) and on GW observations of BNSs (Abbott et al., 2018a)
have been recently combined to yield joint constrains (Miller et al., 2019;
Raaijmakers et al., 2020; Miller et al., 2021; Riley et al., 2021; Raaijmakers
et al., 2021). In order to explore the effect of the choice of EoS, we choose
the EoSs yielding mass-radius relations consistent with the 95% confidence
interval of the mass-radius posterior found in Raaijmakers et al. (2020)
222Note that Raaijmakers et al. (2020) presented two different model-dependent
mass-radius posteriors. To be conservative, we pick all EoSs that are
consistent with either of the posteriors. among the EoSs provided by the
Xtreme Astrophysics Group (Özel et al., 2016; Bogdanov et al., 2016; Özel &
Freire, 2016). The selected EoSs are plotted in Figure 1.

###  2.4. Estimation of astrophysical rates

The merger rate estimates published by the LVK (Abbott et al., 2021c, a)
assume a model for the mass and spin distribution of compact objects. Since we
want to impose different astrophysical models, we don’t use the published
numbers directly, but instead perform hierarchical Bayesian inference, with
the various models detailed in Table 1. Specifically, we aim at estimating the
posterior probability density of the merger rate $\mathcal{R}$ for each
population. For a population model parametrized by a scale parameter
$\mathcal{R}$ (which determines the overall volumetric merger rate of sources)
and shape parameters $\Lambda$ (which determines the distribution of sources’
parameters), this can be written as (Loredo & Lamb, 2002; Mandel et al., 2019;
Vitale et al., 2020):

$p(\mathcal{R}|\Lambda,\\{{d}\\})\propto\pi(\mathcal{R})e^{-N^{s}_{\uparrow}}\left(\frac{N_{\uparrow}^{s}}{\alpha(\Lambda)}\right)^{N^{\mathrm{tr}}}\prod_{i=1}^{N^{\mathrm{tr}}}p\left(d\mid\Lambda,\mathcal{R}\right)$
(2)

Where $N_{\uparrow}^{s}$ is the number of sources detectable during the
observing period, and $N^{\mathrm{tr}}$ is the number of triggers that were
actually detected. Both these quantities depend on $\mathcal{R}$.
$\alpha(\Lambda)$ is the _fraction_ of events that are detectable for a given
value of shape parameters, also known as the selection function. We generate
samples from this posterior distribution using the gwpopulation (Talbot et
al., 2019) software, and a Jeffrey prior –
$\pi(\mathcal{R})\propto\mathcal{R}^{-1}$ – on the rate $\mathcal{R}$.

The LVK have not released all the data that would be required to calculate the
selection function $\alpha(\Lambda)$ in a way that is identical to what’s done
for the LVK papers 333Specifically, the LVK did not release results of
simulation campaigns that assess the efficiency of search algorithms in the
first observing run for BNS, and in any observing runs for NSBH.. Therefore,
we rely on a standard procedure in the literature and assess the detection
efficiency by using the optimal signal-to-noise ratio (SNR) $\rho$. More
specifically, we assume that a binary merger event with true parameters
$\vec{\theta}$ is detectable if its sky and orientation averaged SNR is higher
than 8 in a single LIGO detector (Dominik et al., 2013). To account for the
diverse sensitivity of the GW detectors in their first three observing runs,
we average the detection efficiency of the 3 observing runs as described in
Sec 5.2 of Vitale et al. (2020). To calculate the SNR we use the LIGO-
Livingston power spectral density at the time of GW150914 (for the first
observing run), GW170817 (for the second observing run) and GW190425 (for the
third observing run), as released by Abbott et al. (2021b). We do not include
spins when calculating the SNR, which has been shown to only impact the
detectabilty by less than a few percent (Ng et al., 2018). Finally, we assume
that sources are distributed uniformly in comoving volume, and with isotropic
sky positions and orientations. All mass models use the same distribution for
the NS in the binary, as described at the beginning of Methods section. The
mass spectrum of the BHs is described for each model in Table 1. We stress
that the shape parameters of our models are thus entirely fixed: the only
unknown parameter of each population is its overall merger rate $\mathcal{R}$.
To infer the rate of BNSs we only use GW170817 and GW190425, whereas when
calculating the rate of NSBH we only use GW200105 and GW200115.

##  3\. Results

We estimate the amount of ejecta for the 100,000 simulated BNSs and NSBHs in
each scenario listed in Table 1, add up the ejecta and scale the sum by the
rate ratio of NSBH and BNS, ($R_{\rm NSBH}/R_{\rm BNS}$), to obtain their
relative ejecta ratio ($M_{\rm{ej,NSBH}}/M_{\rm{ej,BNS}}$). We then estimate
the fraction of NSBH ejecta as $M_{\rm{ej,NSBH}}/M_{\rm{ej,Total}}\equiv
M_{\rm{ej,NSBH}}/(M_{\rm{ej,BNS}}+M_{\rm{ej,NSBH}})$. Since these estimations
are subject to different sources of uncertainty, we discuss them separately
below.

Figure 2.— The NSBH ejecta fraction for different ($\alpha_{\rm{dyn,BNS}}$,
$\alpha_{\rm{dyn,NSBH}}$, $f_{\rm{loss,BNS}}$ and $f_{\rm{loss,NSBH}}$)
[labeled respectively in the horizontal axis] when the BNS and NSBH
astrophysical rates are equivalent. We assume the BBH-like mass+spin model and
ap3 NS EoS in this example.

First we show how $\alpha_{\rm{dyn}}$ and $f_{\rm{loss}}$ (Eq. 1) affects the
results. We use the BBH-like mass+spin model and the ap3 NS EoS. Figure 2
shows that smaller $\alpha_{\rm{dyn,BNS}}$ and larger $\alpha_{\rm{dyn,NSBH}}$
lead to a larger NSBH fraction as expected. On the other hand, we find that
the impact of $f_{\rm{loss,BNS}}$ and $f_{\rm{loss,NSBH}}$ on the NSBH ejecta
fraction varies between different NS EoS and the mass and spin models.

Figure 3.— The NSBH ejecta fraction for different EoS shown in Figure 1 when
the BNS and NSBH astrophysical rates are equivalent. We assume the BBH-like
mass+spin model, $\alpha_{\rm{dyn,BNS}}=\alpha_{\rm{dyn,NSBH}}=1$,
$f_{\rm{loss,BNS}}=0.8$, and $f_{\rm{loss,NSBH}}=0.4$ in this example.

Next, we show the effect of varying the NS EoS. We use the BBH-like mass+spin
model and fix $\alpha_{\rm{dyn,BNS}}=\alpha_{\rm{dyn,NSBH}}=1$,
$f_{\rm{loss,BNS}}=0.8$, and $f_{\rm{loss,NSBH}}=0.4$. Figure 3 shows the NSBH
ejecta fraction for different EoS shown in Figure 1. In general, stiffer EoS
leads to more ejecta for both BNSs and NSBHs. Their impact on the NSBH ejecta
fraction depends on the relative ratio between the dynamical ejecta and disk,
which vary with $\alpha_{\rm{dyn}}$, $f_{\rm{loss}}$, and the mass and spin
models.

Figure 4.— The upper bound of NSBH ejecta fraction as a function of the NSBH
and BNS astrophysical rate ratio. The six mass and spin models are summarized
in Table 1. We iterate through different combinations of NS EoS,
$\alpha_{\rm{dyn,BNS}}$, $\alpha_{\rm{dyn,NSBH}}$, $f_{\rm{loss,BNS}}$, and
$f_{\rm{loss,NSBH}}$, and then use the combinations that lead to the largest
NSBH ejecta fractions to set the upper bounds in this figure. We also mark the
intersection between the upper bound and the 90% upper limit of rate ratio for
each models. This intersection represents the upper limit on the NSBH
contribution to the binary merger production of r-process elements. We
summarize the upper limit in Table 1.

In order to determine the overall uncertainty of the ejecta fraction, we
combine different sources of error. Since there is no available study on the
correlation between the sources of error we discuss above, we conservatively
assume they are independent. For each mass and spin model listed in Table 1,
we iterate through different combinations of NS EoS, $\alpha_{\rm{dyn,BNS}}$,
$\alpha_{\rm{dyn,NSBH}}$, $f_{\rm{loss,BNS}}$, and $f_{\rm{loss,NSBH}}$. We
then use the combination that leads to the largest NSBH ejecta fraction to set
the upper bound of our estimations. In Figure 4 we show the upper bound of the
fraction of NSBH ejecta as a function of the rate ratio for each scenario.

We use the LVK observations of BNSs and NSBHs to estimate the astrophysical
rate ratio between these two classes of sources. In Figure 4 we mark the 90%
upper limit of the rate ratio for each scenario. The intersection between the
upper bound of NSBH ejecta fraction and the 90% upper limit of the rate ratio
is our constraint on the fraction of NSBH ejecta. Most of the mass ejected
from BNSs and NSBHs is expected to become r-process elements, with a small
fraction of iron-peak elements that are more likely to be produced in BNSs
than NSBHs if the ejecta is less neutron-rich (electron fraction $Y_{e}\gtrsim
0.3-0.4$, see e.g. Korobkin et al. (2012); Wanajo et al. (2014); Lippuner &
Roberts (2015); Wollaeger et al. (2021)). We thus take the upper limit on the
ejected mass fraction to be a good estimate of the limit on NSBH’s
contribution to the production of r-process elements. The constraints for
different scenarios are summarized in the last column of Table 1.

NSBHs result in more mass ejected when the mass of BHs is small
($<5M_{\odot}$) or the aligned component of BH spins is large. We find that
NSBHs can account for up to 77% of the r-process-element production. The
highest fraction is obtained if there is an excess of 2-3$M_{\odot}$ BHs with
spins aligned with orbital angular momentum and uniformly distributed (i.e.
for the BBH-like mass+aligned spin model) merging in NSBHs at a high rate.
Smaller BH spin magnitudes or larger tilt angles both reduce the aligned
component of the spins, leading to less NSBH ejecta. In these cases, the NSBH
ejecta fraction decreases to 35% or less (the No gap+BBH-like spin and BBH-
like mass+spin models). These models are also more consistent with the NSBHs
observed by LVK so far (Abbott et al., 2021a). If not only BH spins are small,
but there also is a dearth of low-mass ($<5M_{\odot}$) BHs in NSBHs, then BNSs
produce virtually all of the r-process elements from compact binary mergers.
We therefore do not explore models with even larger BH mass gap or smaller BH
spins since they would lead to even less NSBH ejecta.

##  4\. Discussion

In this paper we have combined GW detections and pulsar observations to place
constraints on the relative contribution of BNSs and NSBHs to the production
of r-process elements. NSBHs contribution can be as high as 77%, if BHs in
NSBHs can have low mass (2-3$M_{\odot}$) and large spins aligned with the
orbital angular momentum. However, this low-mass and high-spin BHs scenario
seems disfavoured by GW observations of NSBHs. If most black holes have masses
in excess of $\sim 5~{}M_{\odot}$ and/or small spins, BNS will contribute the
nearly entirety of r-process elements from compact binary mergers. Different
EoS and ejecta models can lower the upper limits by up to 30%, which is
smaller than the differences arising from different mass and spin models.

The relative contribution of BNSs and NSBHs to the production of heavy
elements depends thus significantly on the distribution of black holes masses
and spins in binaries. These will be better measured by upcoming LVK
observations (Abbott et al., 2018b), leading to more precise estimates of the
ejecta ratio. In addition, current LVK observations of BNSs and NSBHs are
local ($z<0.2$), thus our inference constrains the binary merger production of
the r-process elements over the past 2.5 billion years. Future GW detections
at higher redshifts will enable measuring the evolution of BNS and NSBH
astrophysical rates and therefore their redshift-dependent contribution to the
production of heavy elements.

Although our analysis yields an estimate of the BNS ejecta fraction ($\equiv
1-M_{\rm{ej,NSBH}}/M_{\rm{ej,Total}}$), we do not present a lower limit on the
BNS ejecta since other astrophysical channels could still produce some of the
heavy elements. The constraints we have obtained are mainly limited by the
number of available numerical simulations. Whenever different options existed,
we have conservatively take those that yielded the largest error bars. A more
extensive coverage of the physical parameter space by numerical simulations of
BNS and NSBH mergers will thus reduce the uncertainty in the estimation of the
mass ejecta.

This work demonstrates the potential of combining electromagnetic and
gravitational-wave observations of compact objects. It also underlines the
need for more extensive numerical simulations of the mass ejecta from compact
objects. As advanced gravitational-wave detectors will measure with more
precision the mass and spin distribution of compact objects, whilst
gravitational-wave and pulsar observations will yield improved constraints on
the NS EoS, we will be able to significantly refine the constraints presented
here in the months and years to come.

Acknowledgments – We thank Katerina Chatziioannou and Will Farr for help using
the data release associated with Farr & Chatziioannou (2020). We thank Justin
Alsing, Emanuele Berti, Zoheyr Doctor, Anna Frebel, and Brian Metzger for
valuable comments and suggestions. HYC is supported by NASA through NASA
Hubble Fellowship grants No. HST-HF2-51452.001-A awarded by the Space
Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., for NASA, under contract
NAS5-26555. S.V. acknowledges support of the National Science Foundation and
the LIGO Laboratory. LIGO was constructed by the California Institute of
Technology and Massachusetts Institute of Technology with funding from the
National Science Foundation and operates under Cooperative Agreement No.
PHY-1764464. F.F. gratefully acknowledges support from the NSF through grant
PHY-1806278, from the DOE through grant DE-SC0020435, and from NASA through
grant 80NSSC18K0565. This research has made use of data, software and/or web
tools obtained from the Gravitational Wave Open Science Center
(https://www.gw-openscience.org/ ), a service of LIGO Laboratory, the LIGO
Scientific Collaboration and the Virgo Collaboration. LIGO Laboratory and
Advanced LIGO are funded by the United States National Science Foundation
(NSF) as well as the Science and Technology Facilities Council (STFC) of the
United Kingdom, the Max-Planck-Society (MPS), and the State of
Niedersachsen/Germany for support of the construction of Advanced LIGO and
construction and operation of the GEO600 detector. Additional support for
Advanced LIGO was provided by the Australian Research Council. Virgo is
funded, through the European Gravitational Observatory (EGO), by the French
Centre National de Recherche Scientifique (CNRS), the Italian Istituto
Nazionale di Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions
by institutions from Belgium, Germany, Greece, Hungary, Ireland, Japan,
Monaco, Poland, Portugal, Spain.

## References

  * Abbott et al. (2017) Abbott, B. P., et al. 2017, Astrophys. J. Lett., 850, L39
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 119, 161101
  * Abbott et al. (2017) Abbott, B. P., et al. 2017, Astrophys. J. Lett., 848, L12
  * Abbott et al. (2018a) —. 2018a, Phys. Rev. Lett., 121, 161101
  * Abbott et al. (2018b) —. 2018b, Living Rev. Rel., 21, 3
  * Abbott et al. (2021a) Abbott, R., et al. 2021a, Astrophys. J. Lett., 915, L5
  * Abbott et al. (2021b) —. 2021b, SoftwareX, 13, 100658
  * Abbott et al. (2021c) —. 2021c, Astrophys. J. Lett., 913, L7
  * Alsing et al. (2018) Alsing, J., Silva, H. O., & Berti, E. 2018, Monthly Notices of the Royal Astronomical Society, 478, 1377
  * Bogdanov et al. (2016) Bogdanov, S., Heinke, C. O., Özel, F., & Güver, T. 2016, 831, 184
  * Christie et al. (2019) Christie, I. M., Lalakos, A., Tchekhovskoy, A., et al. 2019, Mon. Not. Roy. Astron. Soc., 490, 4811
  * Côté et al. (2021) Côté, B., Eichler, M., Yagüe López, A., et al. 2021, Science, 371, 945
  * Coulter et al. (2017) Coulter, D. A., Foley, R. J., Kilpatrick, C. D., et al. 2017, Science, 358, 1556
  * Cromartie et al. (2020) Cromartie, H. T., Fonseca, E., Ransom, S. M., et al. 2020, Nature Astronomy, 4, 72
  * Dietrich & Ujevic (2017) Dietrich, T., & Ujevic, M. 2017, Class. Quant. Grav., 34, 105014
  * Dietrich et al. (2017) Dietrich, T., Ujevic, M., Tichy, W., Bernuzzi, S., & Bruegmann, B. 2017, Phys. Rev. D, 95, 024029
  * Dominik et al. (2013) Dominik, M., Belczynski, K., Fryer, C., et al. 2013, 779, 72
  * Fahlman & Fernández (2018) Fahlman, S., & Fernández, R. 2018, Astrophys. J. Lett., 869, L3
  * Farr & Chatziioannou (2020) Farr, W. M., & Chatziioannou, K. 2020, Research Notes of the American Astronomical Society, 4, 65
  * Fernández et al. (2020) Fernández, R., Foucart, F., & Lippuner, J. 2020, Monthly Notices of the Royal Astronomical Society, 497, 3221
  * Fernández & Metzger (2013) Fernández, R., & Metzger, B. D. 2013, Monthly Notices of the Royal Astronomical Society, 435, 502
  * Foucart et al. (2018) Foucart, F., Hinderer, T., & Nissanke, S. 2018, 98, 081501
  * Foucart et al. (2013) Foucart, F., Deaton, M. B., Duez, M. D., et al. 2013, 87, 084006
  * Frebel (2019) Frebel, A. 2019, Annals of Physics, 410, 167909
  * Fujibayashi et al. (2020) Fujibayashi, S., Shibata, M., Wanajo, S., et al. 2020, Phys. Rev. D, 102, 123014
  * Holmbeck et al. (2020) Holmbeck, E. M., Hansen, T. T., Beers, T. C., et al. 2020, Astrophysical Journal Supplemets, 249, 30
  * Hotokezaka et al. (2013) Hotokezaka, K., Kiuchi, K., Kyutoku, K., et al. 2013, 87, 024001
  * Ji et al. (2016) Ji, A. P., Frebel, A., Chiti, A., & Simon, J. D. 2016, Nature, 531, 610
  * Just et al. (2015) Just, O., Bauswein, A., Ardevol Pulpillo, R., Goriely, S., & Janka, H. T. 2015, Monthly Notices of the Royal Astronomical Society, 448, 541
  * Kasen et al. (2013) Kasen, D., Badnell, N. R., & Barnes, J. 2013, Astrophys. J., 774, 25
  * Kiuchi et al. (2019) Kiuchi, K., Kyutoku, K., Shibata, M., & Taniguchi, K. 2019, Astrophys. J. Lett., 876, L31
  * Kiuchi et al. (2015) Kiuchi, K., Sekiguchi, Y., Kyutoku, K., et al. 2015, Phys. Rev. D, 92, 064034
  * Korobkin et al. (2012) Korobkin, O., Rosswog, S., Arcones, A., & Winteler, C. 2012, Mon. Not. Roy. Astron. Soc., 426, 1940
  * Kreidberg et al. (2012) Kreidberg, L., Bailyn, C. D., Farr, W. M., & Kalogera, V. 2012, 757, 36
  * Krüger & Foucart (2020) Krüger, C. J., & Foucart, F. 2020, 101, 103002
  * Kyutoku et al. (2015) Kyutoku, K., Ioka, K., Okawa, H., Shibata, M., & Taniguchi, K. 2015, Phys.Rev.D, 92, 044028
  * Landry & Read (2021) Landry, P., & Read, J. S. 2021, arXiv e-prints, arXiv:2107.04559
  * Lattimer & Prakash (2001) Lattimer, J. M., & Prakash, M. 2001, 550, 426
  * Lattimer & Schramm (1974) Lattimer, J. M., & Schramm, D. N. 1974, ApJ, 192, L145
  * Lattimer & Schramm (1976) —. 1976, 210, 549
  * Lippuner & Roberts (2015) Lippuner, J., & Roberts, L. F. 2015, 815, 82
  * Liu & Lai (2021) Liu, B., & Lai, D. 2021, Mon. Not. Roy. Astron. Soc., 502, 2049
  * Loredo & Lamb (2002) Loredo, T. J., & Lamb, D. Q. 2002, 65, 063002
  * Lu et al. (2021) Lu, W., Beniamini, P., & Bonnerot, C. 2021, Mon. Not. Roy. Astron. Soc., 500, 1817
  * Mandel et al. (2019) Mandel, I., Farr, W. M., & Gair, J. R. 2019, Mon. Not. Roy. Astron. Soc., 486, 1086
  * Metzger (2020) Metzger, B. D. 2020, Living Rev. Rel., 23, 1
  * Metzger & Fernández (2014) Metzger, B. D., & Fernández, R. 2014, Mon. Not. Roy. Astron. Soc., 441, 3444
  * Metzger & Fernández (2014) Metzger, B. D., & Fernández, R. 2014, Mon. Not. Roy. Astron. Soc., 441, 3444
  * Metzger et al. (2010) Metzger, B. D., Martínez-Pinedo, G., Darbha, S., et al. 2010, Mon. Not. Roy. Astron. Soc., 406, 2650
  * Meyer (1993) Meyer, B. S. 1993, Meteoritics, 28, 399
  * Miller et al. (2020) Miller, J. M., Sprouse, T. M., Fryer, C. L., et al. 2020, Astrophys. J., 902, 66
  * Miller et al. (2019) Miller, M. C., Lamb, F. K., Dittmann, A. J., et al. 2019, Astrophysical Journal Letters, 887, L24
  * Miller et al. (2021) —. 2021, arXiv e-prints, arXiv:2105.06979
  * Mösta et al. (2018) Mösta, P., Roberts, L. F., Halevi, G., et al. 2018, Astrophys. J., 864, 171
  * Nedora et al. (2021) Nedora, V., Bernuzzi, S., Radice, D., et al. 2021, Astrophys. J., 906, 98
  * Ng et al. (2018) Ng, K. K. Y., Vitale, S., Zimmerman, A., et al. 2018, 98, 083007
  * Özel & Freire (2016) Özel, F., & Freire, P. 2016, Annual Reviews of Astronomy and Astrophysics, 54, 401
  * Özel et al. (2016) Özel, F., Psaltis, D., Güver, T., et al. 2016, 820, 28
  * Paul et al. (2001) Paul, M., et al. 2001, Astrophys. J. Lett., 558, L133
  * Raaijmakers et al. (2020) Raaijmakers, G., Greif, S. K., Riley, T. E., et al. 2020, Astrophysical Journal Letters, 893, L21
  * Raaijmakers et al. (2021) Raaijmakers, G., Greif, S. K., Hebeler, K., et al. 2021, arXiv e-prints, arXiv:2105.06981
  * Riley et al. (2019) Riley, T. E., Watts, A. L., Bogdanov, S., et al. 2019, Astrophysical Journal Letters, 887, L21
  * Riley et al. (2021) Riley, T. E., Watts, A. L., Ray, P. S., et al. 2021, arXiv e-prints, arXiv:2105.06980
  * Roederer & Lawler (2012) Roederer, I. U., & Lawler, J. E. 2012, Astrophysical Journal, 750, 76
  * Siegel et al. (2019) Siegel, D. M., Barnes, J., & Metzger, B. D. 2019, Nature, 569, 241
  * Siegel & Metzger (2017) Siegel, D. M., & Metzger, B. D. 2017, Phys. Rev. Lett., 119, 231102
  * Surman et al. (2006) Surman, R., McLaughlin, G. C., & Hix, W. R. 2006, Astrophys. J., 643, 1057
  * Talbot et al. (2019) Talbot, C., Smith, R., Thrane, E., & Poole, G. B. 2019, 100, 043030
  * Vitale et al. (2020) Vitale, S., Gerosa, D., Farr, W. M., & Taylor, S. R. 2020, arXiv e-prints, arXiv:2007.05579
  * Wallner et al. (2021) Wallner, A., Froehlich, M. B., Hotchkis, M. A. C., et al. 2021, Science, 372, 742
  * Wanajo et al. (2014) Wanajo, S., Sekiguchi, Y., Nishimura, N., et al. 2014, Astrophys.J.Lett., 789, L39
  * Winteler et al. (2012) Winteler, C., Käppeli, R., Perego, A., et al. 2012, Astrophysical Journal Letters, 750, L22
  * Wollaeger et al. (2021) Wollaeger, R. T., Fryer, C. L., Chase, E. A., et al. 2021, arXiv e-prints, arXiv:2105.11543
  * Yong et al. (2021) Yong, D., Kobayashi, C., Da Costa, G. S., et al. 2021, Nature, 595, 223