# Ferromagnetic–antiferromagnetic coexisting ground states and exchange bias
effects in $\bf{MnBi_{4}Te_{7}}$ and $\bf{MnBi_{6}Te_{10}}$

Xiaolong Xu,1,2,3,† Shiqi Yang,1,4,† Huan Wang,5 Roger Guzman,6,7 Yaozheng
Zhu,1 Yuxuan Peng,1 Zhihao Zang,1 Ming Xi,5 Shangjie Tian,5 Yanping Li,1
Hechang Lei,5 Zhaochu Luo,1 Jinbo Yang,1 Tianlong Xia,5,⋆ Wu Zhou,6,7,⋆ Yuan
Huang,3,⋆ and Yu Ye,1,2,8,⋆ 1State Key Laboratory for Mesoscopic Physics and
Frontiers Science Center for Nano-optoelectronics, School of Physics, Peking
University, Beijing 100871, China  2Collaborative Innovation Center of Quantum
Matter, Beijing 100871, China  3Advanced Research Institute of
Multidisciplinary Science, Beijing Institute of Technology, Beijing 100081,
China  4Academy for Advanced Interdisciplinary Studies, Peking University,
Beijing 100871, China  5Department of Physics, Renmin University of China,
Beijing 100872, China  6School of Physical Sciences, University of Chinese
Academy of Sciences, Beijing 100049, China  7CAS Centre for Excellence in
Topological Quantum Computation, University of Chinese Academy of Sciences,
Beijing 100049, China  8Yangtze Delta Institute of Optoelectronics, Peking
University, Nantong 226010 Jiangsu, China

††† These authors contribute equally.††⋆Corresponding to<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>and<EMAIL_ADDRESS>
Natural superlattice structures
$\bf{(MnBi_{2}Te_{4})(Bi_{2}Te_{3})_{\emph{n}}(\emph{n}=1,2,...)}$, in which
magnetic $\bf{MnBi_{2}Te_{4}}$ layers are separated by nonmagnetic
$\bf{Bi_{2}Te_{3}}$ layers, hold band topology, magnetism and reduced
interlayer coupling, providing a promising platform for the realization of
exotic topological quantum states. However, their magnetism in the two-
dimensional limit, which is crucial for further exploration of quantum
phenomena, remains elusive. Here, complex ferromagnetic (FM)–antiferromagnetic
(AFM) coexisting ground states that persist up to the 2-septuple layers (SLs)
limit are observed and comprehensively investigated in $\bf{MnBi_{4}Te_{7}}$
(n = 1) and $\bf{MnBi_{6}Te_{10}}$ (n = 2). The ubiquitous Mn-Bi site mixing
modifies or even changes the sign of the subtle inter-SL magnetic
interactions, yielding a spatially inhomogeneous interlayer coupling. Further,
a tunable exchange bias effect is observed in
$\bf{(MnBi_{2}Te_{4})(Bi_{2}Te_{3})_{\emph{n}}(\emph{n}=1,2)}$, arising from
the coupling between the FM and AFM components in the ground state. Our work
highlights a new approach toward the fine-tuning of magnetism and paves the
way for further study of quantum phenomena in
$\bf{(MnBi_{2}Te_{4})(Bi_{2}Te_{3})_{\emph{n}}(\emph{n}=1,2,...)}$ as well as
their magnetic applications.

The interplay between magnetism and topology inaugurates a new horizon in
exploring exotic quantum phenomena, such as the quantum anomalous Hall effect
(QAHE), axion insulators and magnetic Weyl semimetalsN40 ; N1 ; N2 ; N3 ; N4 ;
N5 . Recently, $\rm{MnBi_{2}Te_{4}}$ was discovered to be an intrinsic
stoichiometric antiferromagnetic (AFM) topological insulator (TI)N7 ; N8 ; N10
; N6 ; N9 ; N11 ; N39 . The intertwined band topology and magnetic order in
A-type AFM $\rm{MnBi_{2}Te_{4}}$ (MBT) pose great challenges to the
realization of its exotic topological physics and subsequent devices, because
either the layer number of the material needs to be strictly controlled, or a
high magnetic field ($\sim$ 6 to 8 T) is required N6 ; N9 ; N39 ; N11 .
Therefore, there is an urgent need to develop a magnetic topological insulator
of the MBT family that is magnetically insensitive to the number of layers and
has a small saturation field, in which engineering of the interlayer magnetic
interaction is the key. Now, the natural superlattice structures
$\rm{(MnBi_{2}Te_{4})(Bi_{2}Te_{3})_{\emph{n}}(\emph{n}=1,2,...)}$ provides an
opportunity to modify the interlayer exchange coupling using non-magnetic
$\rm{Bi_{2}Te_{3}}$ (BT) quintuple layer (QL) as spacer layersN37 ; N12 ; N13
; N41 ; N14 ; N15 . With the increase of n, the interlayer AFM coupling
gradually weakens, while the global spin-orbit coupling strength gradually
increases as the increase of Bi contentN12 ; N13 ; N15 ; N14 ; N41 ; N44 . For
$\rm{MnBi_{4}Te_{7}}$ (n = 1) and $\rm{MnBi_{6}Te_{10}}$ (n = 2), neutron
scattering experiments and theoretical calculations show that they possess
weak interlayer magnetic coupling, but still maintain the A-type AFM
structureN14 ; N15 ; N16 ; N17 ; N18 , that is, in each $\rm{MnBi_{2}Te_{4}}$
septuple layer (SL), the Mn magnetic moments are ferromagnetically aligned,
while the adjacent SLs are antiferromagnetically aligned. Although tremendous
efforts have been devoted to studying the magnetic properties of the
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$, most studies have only
focused on their bulk phaseN12 ; N13 ; N14 ; N41 ; N15 ; N19 ; N35 ; N42 ; N46
. However, the magnetic property and topological phase in MBT family often
exhibit intricate thickness dependence, posing a significant influence for the
realization of exotic topological physicsN6 ; N9 ; N11 ; N21 ; N39 ; N22 ; N20
. Therefore, determining the magnetism of
$\rm{(MnBi_{2}Te_{4})(Bi_{2}Te_{3})_{\emph{n}}(\emph{n}=1,2)}$ at their 2D
limit is crucial, but has so far remained elusive.

In this work, we systematically investigate the magnetic properties of
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ thin flakes down to 1-SL by
employing polar reflective magnetic circular dichroism (RMCD) spectroscopy. We
demonstrate that the odd-even-layer oscillation of compensated and
uncompensated AFM states, a common effect in ideal A-type AFM materialsN11 ;
N20 ; N23 ; N24 , vanishes in atomically thin $\rm{MnBi_{4}Te_{7}}$ and
$\rm{MnBi_{6}Te_{10}}$. In all the measured samples above 1 SL, in addition to
the AFM spin-flip transition, a significant FM hysteresis loop centered at
zero field is observed. We reveal that the peculiar hysteresis loop arises
from the complex FM-AFM coexisting ground state. Atomic-resolution electron
energy loss spectroscopy (EELS) mapping as well as single-crystal X-ray
diffraction (SC-XRD) reveal the ubiquitous Mn-Bi site mixing in the crystals.
The spins of the randomly distributed $\rm{Mn_{Bi}}$ antisite defects in each
SL are antiferromagnetically coupled to the spins of the main Mn layerN19 ;
N25 ; N18 , that further modify or even change the sign of the subtle inter-SL
magnetic interactions, yielding a spatially inhomogeneous interlayer coupling.
As a result, the energy gained through the formation of magnetic domains
compensates for the energy cost in the appearance of domain walls, yielding a
complex FM-AFM coexisting ground state. Further, a tunable exchange bias
effect is observed in $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$,
arising from the coupling between the FM and AFM components in the ground
state. The direction of this exchange bias can be easily tuned by the
historical polarization field and does not require warming and field cooling
processes. Our results reveal the nontrivial magnetic ground states in
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$, keeping the promise for
further investigation of exotic topological quantum states.

Figure 1: FM-AFM coexisting magnetic ground state in $\bf{MnBi_{6}Te_{10}}$.
a, A typical optical image of the exfoliated $\rm{MnBi_{6}Te_{10}}$ sample
with different thicknesses and terminations. b, Atomic force microscopy height
image and height line profiles measured in the area marked by the black box in
a, showing a thicknesses of $\sim$1.4 nm for $\rm{MnBi_{2}Te_{4}}$ SL and
$\sim$1.0 nm for $\rm{Bi_{2}Te_{3}}$ QL. c, RMCD signals versus $\mu_{0}H$ for
1 to 4-SL (terminated by 2-QL protection layer) $\rm{MnBi_{6}Te_{10}}$ flakes
at 2 K. d, RMCD sweeps for the 3-SL $\rm{MnBi_{6}Te_{10}}$ flake at a
temperature range that passes through its $T_{\rm{N}}$. e, 2D plot of RMCD
signal as functions of temperature and magnetic field (descending sweep) for
the 3-SL $\rm{MnBi_{6}Te_{10}}$ flake. Two AFM spin-flip fields
($H_{\rm{f-}}^{1}$ and $H_{\rm{f-}}^{2}$) and one FM spin-flip field
($H_{\rm{c-}}$) are extracted and superimposed on the plot. f,g, FM
($f_{\rm{FM}}$, left axis) and AFM ($f_{\rm{AFM}}$, right axis) component
proportions as functions of thickness (f) and temperature (g). Cross marks in
f indicate the expected values for the single-domain case. The vertical error
bars are estimated from the uncertainty of the measured signal.

FM-AFM coexisting magnetic ground states in $\bf{MnBi_{4}Te_{7}}$ and
$\bf{MnBi_{6}Te_{10}}$ samples. High-quality $\rm{MnBi_{4}Te_{7}}$ and
$\rm{MnBi_{6}Te_{10}}$ bulk crystals are grown by flux methodN33  and
confirmed by XRD results (Fig. S1). The field-cooled (FC) and zero-field-
cooled (ZFC) magnetic susceptibilities of $H\parallel c$ ($\chi^{c}$) of
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ crystals show that their
long-range AFM orders occur at Néel temperature ($T_{\rm{N}}$) of 12.1 K and
10.9 K, respectively (Fig. S2). Compared with $\rm{MnBi_{2}Te_{4}}$N16 ; N17
($T_{\rm{N}}$ of $\sim$ 24.5 K), the reduction of $T_{\rm{N}}$ in
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ manifests the weakened
interlayer coupling. Scrutiny of the $M-H$ curve at 2 K under $H\parallel c$,
we find that $\rm{MnBi_{4}Te_{7}}$ ($\rm{MnBi_{6}Te_{10}}$) undergoes a spin
transition at a very low magnetic field and quickly enters the forced FM state
with a field of about 0.25 T (0.21 T). This is in sharp contrast to
$\rm{MnBi_{2}Te_{4}}$, where the spin-flop transition occurs at about 3.5 T
and its magnetic moment eventually saturates under external magnetic fields
larger than 8 TN6 ; N9 ; N11 ; N20  (Fig. S3a). The much smaller saturation
fields again indicate a significantly reduced interlayer AFM exchange
interactions in $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$.

Figure 2: STEM characterizations and an intuitive model of inter-SL FM and AFM
orders in $\bf{MnBi_{4}Te_{7}}$ and $\bf{MnBi_{6}Te_{10}}$. a, Atomic-
resolution HAADF image of the cross-section of the $\rm{MnBi_{6}Te_{10}}$
crystal along the [100] direction. Experimental (shaded area) and simulated
(grey curve) integral HAADF intensity profiles along the $c$-axis show clear
discrepancy, showing $\rm{Bi_{Mn}}$ antisite defects in SL (red arrow) and
$\rm{Mn_{Te}}$ in QL (grey arrow). b, Atomic structure and HAADF image of the
$\rm{MnBi_{6}Te_{10}}$ crystal with the corresponding EEELs mapping of the Mn
element ($\rm{L_{2,3}}$ edge). Clear Mn signals are present in the Bi layers
of SL and QL. The orange dashed lines indicate the van der Waals gaps. c,d,
Illustrations of the inter-SL FM and AFM regions in $\rm{MnBi_{4}Te_{7}}$ and
$\rm{MnBi_{6}Te_{10}}$, in which the dominant magnetic exchange interactions
are labelled.

To explore how the weakened interlayer exchange interaction affects the
magnetic ground order, we investigated their layer-number-dependent magnetism
of $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ using the RMCD
spectroscopy. Atomically thin flakes down to 1 SL are mechanically exfoliated
from their bulk crystals onto gold substrates using the standard Scotch tape
methodN27 , and then protected with a layer of polymethyl methacrylate (PMMA)
(Fig. 1a). The number of layers is confirmed by atomic force microscopy
characterizations (see Fig. S4 for $\rm{MnBi_{4}Te_{7}}$ samples and Fig. S5
for $\rm{MnBi_{6}Te_{10}}$ samples). Typical height line profiles of the
stepped $\rm{MnBi_{6}Te_{10}}$ flakes indicate SL and QL thicknesses of
approximately 1.4 and 1.0 nm, respectively (Fig. 1b), consistent with previous
reports N22 ; N28 . The layer-number-dependent magnetic behavior in
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ is very similar, so here we
show 1-4 SLs $\rm{MnBi_{6}Te_{10}}$ as examples, and provide the full set of
other samples in Figs. S6 to S9. 1-SL $\rm{MnBi_{6}Te_{10}}$ shows a distinct
RMCD signal at zero field and a clear hysteresis loop (Fig. 1c), confirming
its FM nature. With increasing temperature, the hysteresis loop shrinks and
disappears at 10 K, indicating an FM to paramagnetic phase transition (Fig.
S8a). Compared with the intrinsic 1-SL $\rm{MnBi_{2}Te_{4}}$ ($T_{\rm{C}}$ =
14.5 K)N20 , the decreased $T_{\rm{C}}$ may be due to the increased
$\rm{Bi_{Mn}}$ sites, as the intralayer exchange coupling decreases with the
increase of the average distance between the Mn sites. Surprisingly, with
increasing thickness, no odd-even layer-number oscillation occurs, which is a
typical feature of ideal A-type AFM materials. Instead, both FM and AFM
signals are present in all the measured samples with layer number $N$>1 (Fig.
1c). In the 2-SL samples, one AFM and one FM spin-flip transitions occur
during the magnetic reversal process, signifying the coexistence of ↑↑ and ↑↓
domains. With increasing thickness, the magnetic reversal curve evolves into
two AFM and one FM spin-flip transitions with different magnetic moment
amplitudes.

As the temperature increases, the value of the spin-flip field
$H_{\rm{f-}}^{1}$ in the 3-SL sample changes from negative to positive, and
then slowly decreases to approach zero (Fig. 1d, e). A positive value of
$H_{\rm{f-}}^{1}$ signifies the AFM nature of this spin-flip process, since
interlayer AFM coupling prefers to make the magnetic moment in adjacent layers
be antiparallel, while the Zeeman energy tends to keep the magnetic moment
parallel to the magnetic field. Since the AFM interlayer coupling is stronger
in the $\rm{MnBi_{4}Te_{7}}$ samples, $H_{\rm{f-}}^{1}$ usually occurs at
larger positive values (see Figs. S6 and S9). Different coupling types and
coupling strengths between AFM and FM components will result in distinct
hysteresis behaviors. The FM spin-flip field $H_{\rm{c-}}$ of the
$\rm{MnBi_{6}Te_{10}}$ samples increase monotonically from a negative value
and approach zero with increasing temperature (Figs. 1d, e and Fig. S8b).
However, in some $\rm{MnBi_{4}Te_{7}}$ samples (Fig. S9d), $H_{\rm{c-}}$ first
increases from negative to positive and then jumps back to negative with
increasing temperature, signifying the correlated coupling between the FM and
AFM components. Specifically, we build a five-layer macrospin model to
interpret the hysteresis loops in the $\rm{MnBi_{4}Te_{7}}$ and
$\rm{MnBi_{6}Te_{10}}$ samples (Fig. S10), which captures the main features.

In real materials, the FM and AFM distributions and proportions can be very
complex due to the spatial inhomogeneity. We denote the proportions of FM and
AFM components by $f_{\rm{FM}}=\frac{\text{RMCD(FM flip
height)}}{\text{RMCD(up saturation)}-\text{RMCD(down saturation)}}$ and
$f_{\rm{AFM}}=1-f_{\rm{FM}}$, respectively. As expected, $f_{\rm{AFM}}\sim 0$
($f_{\rm{FM}}\sim 1$) in the 1-SL sample (Fig. 1f). The $f_{\rm{AFM}}$
($f_{\rm{FM}}$) of $\rm{MnBi_{6}Te_{10}}$ samples increases (decreases) as the
number of layers increases. In most samples, the values of $f_{\rm{AFM}}$ do
not coincide with those expected for single-domain case (cross marks in Fig.
1f), indicating a multi-domain configuration within the laser spot size. The
values of $f_{\rm{FM}}$ and $f_{\rm{AFM}}$ are nearly temperature independent
over the entire measurement temperature range (Fig. 1g).

After confirming the FM-AFM coexistence ground states, the puzzling magnetism
of bulk $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ can also be well
explained. It is worth noting that at low temperatures, the
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ crystals exhibit non-zero
magnetization at zero field (Fig. S3), and the magnetic reversal completes
through three sluggish spin-flip transitions (marked by the arrows in Fig.
S3)N7 ; N12 ; N13 ; N14 ; N16 ; N22 ; N41 . The non-zero magnetization at the
zero field suggests that there may be some FM components in the AFM sate. The
FM-AFM coexisting magnetic order is also confirmed by the bifurcations of the
ZFC and FC curves at temperatures slightly below the Néel temperature (Fig.
S2). The large difference in the values of $f_{\rm{AFM}}$ between thin flakes
and bulk crystals suggests that the AFM-FM coexisting magnetic order is
possibly a surface-related effect. The complex multi-domain structure in the
bulk samples smears the distinct signatures of the FM-AFM coexisting ground
state, resulting in sluggish changes in the magnetic moment, which also masks
the potential applications of this unique magnetic order.

Figure 3: Spatial RMCD mappings of $\bf{(MBT+2BT)_{3}}$ sample. a, RMCD signal
versus external magnetic field of a 3-SL $\rm{MnBi_{6}Te_{10}}$ flake at 4 K.
Four magnetic fields (labeled as 1 to 4) corresponding to four different spin
configurations are selected for RMCD maping. b, Optical image of the stepped
$\rm{MnBi_{6}Te_{10}}$ sample. The selected area for RMCD mapping is labeled
by the dashed red box. c-e, Spatial distributions of the RMCD signal changes
in two AFM spin-flip transitions (1 to 2 and 3 to 4) and one FM spin-flip
transition (2 to 3). The homogeneous signals indicate the uniform FM-AFM
coexistence under the experimental spatial resolution. The bare substrate and
etched hole are marked by black dashed lines.

Origin of FM-AFM coexisting magnetic ground state. Due to the weak interlayer
coupling, the energy difference between the AFM and FM spin configurations is
very small N14 ; N15 , providing a breeding ground for the magnetism tuningN43
; N45 . Cross-sectional atomic-resolution high-angle annular dark-field
scanning transmission electron microscopy (HAADF-STEM) image along the [100]
direction shows that two $\rm{Bi_{2}Te_{3}}$ layers are inserted between
$\rm{MnBi_{2}Te_{4}}$ layers in $\rm{MnBi_{6}Te_{10}}$ (Figs. 2a and b).
However, the experimental and simulated integrated HAADF intensity profiles
along the $c$-axis show clear discrepancy, showing $\rm{Bi_{Mn}}$ antisite
defects in SL (red arrow in Fig. 2a) and $\rm{Mn_{Te}}$ in QL (grey arrow in
Fig. 2a). In addition, atomic electron energy loss spectroscopy (EELS) mapping
shows clear Mn signals in the Bi layers in SL and also QL (Fig. 2b),
demonstrating the existence of the $\rm{Mn_{Bi}}$ antisite defects (see Fig.
S11 for $\rm{MnBi_{4}Te_{7}}$).

Below ${T\rm{{}_{N}}}$, the spins of Mn atoms are ferromagnetically coupled
($-J_{1}$) within the SL, but antiferromagnetically coupled ($J_{4}$) with
those of Mn atoms in the adjacent SL. The spins of $\rm{Mn_{Bi}}$ in SL are
ferromagnetically coupled ($-J_{2}$) within the sublayer, but
antiferromagnetically coupled ($J_{3}$) to the main Mn layerN25 , forming a
ferrimagnetic SL configuration (Fig. 2c). Flipping the spins of $\rm{Mn_{Bi}}$
in SL to align with those in the Mn layer requires a magnetic field above 50
TN25 , so all spins in SL flip together under a small magnetic field. This
stable ferrimagnetic SL configuration results in a reduced measured saturation
magnetic moment than the theoretically valueN14  (Fig. S3), and determines
that the magnetic moment change at low magnetic fields (from $-$0.2 T to 0.2
T) comes only from that of the main Mn layer. Considering that Mn atoms doped
into the $\rm{Bi_{2}Te_{3}}$ always tend to be ferromagnetically coupledN29 ;
N30 ; N31 , the spins of $\rm{Mn_{Bi}}$ are expected to be ferromagnetically
coupled ($-J_{5}$) in the adjacent SLs and within QLs. The magnetic moment in
the QL to flip to align with the magnetic filed requires a magnetic field of
about 8 T and 6 T in $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$N19 ,
respectively (marked by the arrows in Fig. S3). For simplification, we do not
additionally consider the coupling in QL and accommodate it in $J_{5}$, which
effectively strengthens the FM coupling between the $\rm{Mn_{Bi}}$ antisites
in the two neighboring SLs. Then, the Hamiltonian of the 2-SL system can be
expressed as:

$\begin{split}H&=-J_{1}\sum_{<ij>}\vec{M}_{i}^{1}\cdot\vec{M}_{j}^{1}-J_{1}\sum_{<ij>}\vec{M}_{i}^{2}\cdot\vec{M}_{j}^{2}\\\
&-J_{2}\sum_{<ij>}\vec{m}_{i}^{1}\cdot\vec{m}_{j}^{1}-J_{2}\sum_{<ij>}\vec{m}_{i}^{2}\cdot\vec{m}_{j}^{2}\\\
&+J_{4}\sum_{<ij>}\vec{M}_{i}^{1}\cdot\vec{M}_{j}^{2}-J_{5}\sum_{<ij>}\vec{m}_{i}^{1}\cdot\vec{m}_{j}^{2}\\\
&+J_{3}\sum_{<ij>}\vec{M}_{i}^{1}\cdot\vec{m}_{j}^{1}+J_{3}\sum_{<ij>}\vec{M}_{i}^{2}\cdot\vec{m}_{j}^{2}\end{split}$
(1)

where $m_{i}$ and $M_{i}$ are the magnetic moments at each $\rm{Mn_{Bi}}$ and
Mn lattice sites, superscripts 1 and 2 represent the SL layer indices, and the
sum is over all the nearest-neighboring lattice sites. In the presence of
$\rm{Mn_{Bi}}$ antisite defects, the energy difference per unit cell between
the AFM and FM configurations (Figs. 2c and d) can be expressed as:

$\displaystyle
U_{\rm{AFM}}-U_{\rm{FM}}=2(J_{5}\sum_{<ij>}m_{i}^{1}m_{j}^{2}-J_{4}\sum_{<ij>}M_{i}^{1}M_{j}^{2})$
(2)

This energy difference increases with the amount of Mn-Bi site mixing (the
numbers of $M_{i}$ sites decreases and the number of $m_{i}$ sites increases),
eventually leading to a change in the sign of the interlayer magnetic coupling
(from AFM to FM). The randomly distributed $\rm{Mn_{Bi}}$ antisites lead to a
spatially inhomogeneous interlayer coupling. As a result, the FM-AFM
coexisting ground states are expected to emerge when the energy gain from
forming FM domains exceeds the energy cost from forming domain walls. In
general, the lattice defect concentration of the MBT family is difficult to
precisely control due to the narrow growth temperature window and the large
size difference between Mn and Bi ions. Then, in the isostructural
$\rm{MnSb_{2}Te_{4}}$ with a larger growth temperature windowN32 , we observe
the evolution of A-type AFM to FM-AFM coexistence and finally to the FM ground
state with varying the Mn-Sb site-mixing concentration (Fig. S12).

Domain size and distribution characterizations. To evaluate the domain sizes
of the FM and AFM components, we characterize the magnetic spatial homogeneity
by RMCD mapping. In the typical RMCD-$\mu_{0}H$ curve of a 3-SL
$\rm{MnBi_{6}Te_{10}}$ sample (Fig. 3a), we map the RMCD signals in a selected
area (Fig. 3b) under four selected magnetic fields (0.1 T, 0 T, $-$0.05 T, and
$-$0.1 T, respectively) corresponding to four different spin configurations
(see Fig. S13 for $\rm{MnBi_{4}Te_{7}}$ sample). Changes in the RMCD signal at
the two AFM spin-flip transitions (Fig. 3c and Fig. 3e) and the FM spin-flip
transition (Fig. 3d) are uniform across the scanned sample area, indicating
homogeneous AFM-FM coexistence at a spatial resolution limited by the laser
spot size of $\sim$ 2 $\mu$m. The small domain size of the FM and AFM
components is consistent with the large number of Mn-Bi site mixing
characterized by single-crystal XRD (see Tables. I and II).

Figure 4:  Exchange bias effects of the FM component in $\bf{MnBi_{4}Te_{7}}$
and $\bf{MnBi_{6}Te_{10}}$ samples. a,b, Exchange bias effects of the FM
component in the $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ measured at
4 K. The direction of the exchange bias can be tuned by the historical
polarization field. The exchange bias effect observed in
$\rm{MnBi_{6}Te_{10}}$ is opposite to that of $\rm{MnBi_{4}Te_{7}}$. Colored
dashed lines represent $H_{\text{E}}$. c,d, Temperature-dependent exchange
bias effects of $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$, with the
full hysteresis loops being presented as the references. The red and green
arrows represent the $H_{\rm{c-}}$ and $H_{\rm{c+}}$ for the biased hysteresis
loop, respectively, and the blue dashed lines represent $H_{\text{E}}$. e,f,
Evolution of $H_{\rm{c}}$ and $H_{\text{E}}$ with temperature in
$\rm{MnBi_{4}Te_{7}}$ (blue squares) and $\rm{MnBi_{6}Te_{10}}$ (green
squares).

Tunable exchange bias from the FM-AFM coexisting ground state. From the
temperature-dependent hysteresis loops, we conclude that the coexisting FM and
AFM components are correlated, which provides us with a unique system to study
exchange bias effects. The large-field full hysteresis loops at 4 K (grey data
in Figs. 4a-d) are plotted as references for the minor hysteresis loops of the
FM components. Historically polarized by a large positive saturation magnetic
field, the minor hysteresis loop of the FM component shifts to the right side
in $\rm{MnBi_{4}Te_{7}}$ sample (blue data in Fig. 4a). The hysteresis
behavior of the magnetic field descending sweeping is time reversal to that of
the magnetic field ascending sweeping. Historically polarized by a large
negative saturation magnetic field, then the minor hysteresis loop of the FM
component shifts to the left, indicating a negative exchange bias (orange data
in Figs. 4b). The direction of this exchange bias can be easily tuned by the
historical polarization field and does not require warming and field cooling
processes. Moreover, this exchange bias is very stable under multiple back-
and-forth magnetic field sweeps, with no training effect (Fig. S14). The
exchange bias effect observed in $\rm{MnBi_{6}Te_{10}}$ is opposite to that of
$\rm{MnBi_{4}Te_{7}}$, possibly due to the different magnetic structures at
the FM/AFM interfaces because of the different interlayer coupling strengths
(see discussion following Fig. S10).

As the temperature increases, the coercive field $H_{\text{c}}$ for the biased
hysteresis loop, defined as $1/2(H_{\text{c+}}-H_{\text{c-}})$, gradually
shrinks and the exchange bias field $H_{\text{E}}$, defined as
$1/2(H_{\text{c+}}+H_{\text{c-}})$, slightly increases (Figs. 4c-f). The
$H_{\text{c}}$ shows an abnormal increase at 10 K for $\rm{MnBi_{4}Te_{7}}$
and 8 K for $\rm{MnBi_{6}Te_{10}}$, respectively, accompanied by the
vanishment of the exchange bias effect. In typical AFM and FM heterostructure
systems, the exchange bias effect occurs only when
$(K_{\text{AFM}}t_{\text{AFM}})/{J_{\text{int}}}\geq 1$, where
$K_{\text{AFM}}$ is the anisotropy energy of the AFM component,
$t_{\text{AFM}}$ is the thickness of the AFM component, and $J_{\text{int}}$
is the exchange coupling between the AFM and FM components. As the temperature
increase, the $K_{\text{AFM}}$ in the MBT system decreases, and eventually,
the AFM pinning layer flips collectively with the FM spins, resulting in the
vanishment of the exchange bias and the abnormal increase in the coercive
field. At high temperatures, the two AFM spin-flip transitions lead to
identical RMCD signal changes, also confirming the collective flipping of the
AFM and FM components during the FM spin-flip transition.

In summary, we systematically study the magnetism of atomically thin
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ flakes in the parameter space
of layer number, temperature and applied magnetic field using RMCD
spectroscopy. The complex FM-AFM coexisting ground state is observed and
verified. The weakened interlayer coupling and inhomogeneously distributed
ubiquitous Mn-Bi site mixing have been attributed to result in such a unique
magnetic ground state. A tunable exchange bias effect is observed in
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$, arising from the coupling
between the FM and AFM components. The demonstrated tunable exchange bias by
introducing the coexistence of interlayer AFM and FM coupling provides design
principles and materials for spintronic devices. Using sophisticated
techniques, synthetic antiferromagnets composed of two or more ferromagnetic
layers separated by non-magnetic spacers can be precisely prepared. Due to the
weak interlayer exchange coupling in synthetic antiferromagnets, introducing
spatial inhomogeneity in thickness of the spacer layer or disorders can
therefore tune the interaction and form an FM-AFM coexisting ground state,
allowing for precise manipulation of the exchange bias effect in this
systemN34 . By unraveling the puzzling magnetic states in
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ flakes, our findings pave the
way for further investigation of quantum phenomena intertwined with the
magnetic orders.

##  I Methods

Crystal growth and magnetic characterizations.
$\rm{(MnBi_{2}Te_{4})(Bi_{2}Te_{3})_{\emph{n}}(\emph{n}=1,2,...)}$ single
crystals were grown by flux method N33 . Mn powder, Bi lump and Te lump were
weighed with the ratio of Mn: Bi: Te = 1: 8: 13 (MnTe: $\rm{Bi_{2}Te_{3}}$ =
1: 4). The mixture was loaded into a corundum crucible sealed into a quartz
tube. The tube was then placed into a furnace and heated to 1100 °C for 20 h
to allow sufficient homogenization. After a rapid cooling to 600 °C at 5 °C/h,
the mixture was cooled slowly to 585 °C (581 °C) at 0.5 °C/h for
$\rm{MnBi_{4}Te_{7}}$ ($\rm{MnBi_{6}Te_{10}}$) and kept at this temperature
for 2 days. Finally, the single crystals were obtained after centrifuging. The
centimeter-scale plate-like $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$
single crystals can be easily exfoliated. Magnetic measurements of
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ single crystals were measured
by the vibrating sample magnetometer (VSM) option in a Quantum Design Physical
Properties Measurement System (PPMS-9 T). The temperature-dependent
magnetization measurements are described in detail in Supplementary Materials.

Preparation of the ultra-thin samples. The $\rm{MnBi_{4}Te_{7}}$ and
$\rm{MnBi_{6}Te_{10}}$ flakes with different thicknesses were first
mechanically exfoliated on a polydimethylsiloxane (PDMS) substrate by the
Scotch tape method. The exfoliated samples on PDMS substrates were then dry-
transferred onto 285 nm $\rm{SiO_{2}}$/Si substrates with evaporated gold
films. Then, a layer of PMMA was spin-coated on the thin flakes for
protection.

AFM characterization. The thickness of the ultra-thin samples wwas verified by
the atomic force microscopy characterization using the Oxford Cypher S system
in tapping mode. According to the height line profiles, the
$\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ were confirmed to possess an
alternated lattice structure of BT (1 nm) + MBT (1.4 nm) and BT (1 nm) + BT (1
nm) + MBT (1.4 nm). See more details in Supplementary Materials.

RMCD measurements. The RMCD measurements were performed based on the Attocube
closed-cycle cryostat (attoDRY2100) down to 1.6 K and up to 9 T in the out-of-
plane direction. The linearly polarized light of 633 nm HeNe laser was
modulated between left and right circular polarization by a photoelastic
modulator (PEM) and focused on the sample through a high numerical aperture
(0.82) objective. The reflected light was detected by a photomultiplier tube
(THORLABS PMT1001/M). The magnetic reversal under external magnetic field was
detected by the RMCD signal determined by the ratio of the a.c. component of
PEM at 50.052 kHz and the a.c. component of chopper at 779 Hz (dealt by a two-
channel lock-in amplifier Zurich HF2LI). The errors of ratio of FM and AFM
components are determined by the instability of the data acquired during RMCD
measurements.

STEM characterization. Samples for cross-sectional investigations were
prepared by standard lift-out procedures using an FEI Helios NanoLab G3 CX
focused ion beam system. To minimize sidewall damage and make the samples
sufficiently thin to be electronically transparent, final milling was carried
out at a voltage of 5 kV and a fine milling at 2 kV. Aberration-corrected STEM
imaging was performed using a Nion HERMES-100 operating at an acceleration
voltage of 60 kV and a probe forming semi-angle of 32 mrad. HAADF images were
acquired using an annular detector with a collection semi-angle of 75-210
mrad. EELS measurements were performed using a collection semi-angle of 75
mrad, an energy dispersion of 0.3 eV per channel, and a probe current of
$\sim$20 pA. The Mn-$L$ (640 eV) and Te-$M$ (572 eV) absorption edges were
integrated for elemental mapping after background subtraction. The original
spectrum images were processed to reduce random noise using a principal
component analysis (PCA) tool. HAADF image simulations were computed using the
STEM_CELL software simulation package matching the microscope experimental
settings described above and using a supercell with a thickness $\sim$20 nm.

##  II DATA AVAILABILITY

The data that support the findings of this study will be available at an open-
access repository with a doi link, when accepted for publishing.

##  III ACKNOWLEDGEMENT

This work was supported by the National Key R&D Program of China (Grants No.
2018YFA0306900, No. 2017YFA0206301, 2019YFA0308602, No. 2019YFA0308000, and
2018YFA0305800), the National Natural Science Foundation of China (Grants No.
62022089, No. 12074425, and No. 11874422), Strategic Priority Research Program
(B) of the Chinese Academy of Sciences (Grant No. XDB33000000), Beijing
Natural Science Foundation (Grant No. JQ21018 and BJJWZYJH01201914430039), and
the fundamental Research Funds for the Central Universities (E1E40209).

##  IV Author contributions

Y.Y., S.Y., and X.X. conceived the project, designed the experiments, analyzed
the results and wrote the manuscript. S.Y. and X.X. conducted the RMCD
measurements. H.W. and T.X. grew the $\rm{MnBi_{4}Te_{7}}$ and
$\rm{MnBi_{6}Te_{10}}$ bulk crystals. M.X., S.T., and H.L. grew the
$\rm{MnSb_{2}Te_{4}}$ bulk crystal. Y.H. prepared the few-layer samples. Y.P.
and J.Y. performed the magnetic characterizations of the bulk crystals. R.G.
performed the STEM characteristics under the supervision of W.Z. Y.Z. and Z.L.
helped with results analysis. All authors discussed the results and
contributed to the manuscript.

##  V ADDITIONAL INFORMATION

Competing interests: The authors declare no competing financial interests.

## References

  * (1) Tokura, Y., Yasuda, K. & Tsukazaki, A.  Magnetic topological insulators.  _Nature Reviews Physics_ 1, 126–143 (2019).
  * (2) Nomura, K. & Nagaosa, N.  Surface-quantized anomalous hall current and the magnetoelectric effect in magnetically disordered topological insulators.  _Physical review letters_ 106, 166802 (2011).
  * (3) Chang, C.-Z. _et al._ Experimental observation of the quantum anomalous hall effect in a magnetic topological insulator.  _Science_ 340, 167–170 (2013).
  * (4) Li, R., Wang, J., Qi, X.-L. & Zhang, S.-C.  Dynamical axion field in topological magnetic insulators.  _Nature Physics_ 6, 284–288 (2010).
  * (5) Liu, E. _et al._ Giant anomalous hall effect in a ferromagnetic kagome-lattice semimetal.  _Nature physics_ 14, 1125–1131 (2018).
  * (6) Liu, D. _et al._ Magnetic weyl semimetal phase in a kagomé crystal.  _Science_ 365, 1282–1285 (2019).
  * (7) Li, J. _et al._ Intrinsic magnetic topological insulators in van der waals layered $\rm{MnBi_{2}Te_{4}}$-family materials.  _Sci. Adv._ 5, eaaw5685 (2019).
  * (8) Otrokov, M. M. _et al._ Prediction and observation of an antiferromagnetic topological insulator.  _Nature_ 576, 416–422 (2019).
  * (9) Deng, H. _et al._ High-temperature quantum anomalous hall regime in a $\rm{MnBi_{2}Te_{4}}$/$\rm{Bi_{2}Te_{3}}$ superlattice.  _Nature Physics_ 17, 36–42 (2021).
  * (10) Deng, Y. _et al._ Quantum anomalous hall effect in intrinsic magnetic topological insulator $\rm{MnBi_{2}Te_{4}}$.  _Science_ 367, 895–900 (2020).
  * (11) Liu, C. _et al._ Robust axion insulator and chern insulator phases in a two-dimensional antiferromagnetic topological insulator.  _Nat. Mater._ 19, 522–527 (2020).
  * (12) Ovchinnikov, D. _et al._ Intertwined topological and magnetic orders in atomically thin chern insulator $\rm{MnBi_{2}Te_{4}}$.  _Nano Lett._ 21, 2544–2550 (2021).
  * (13) Ge, J. _et al._ High-chern-number and high-temperature quantum hall effect without landau levels.  _National science review_ 7, 1280–1287 (2020).
  * (14) Wu, J. _et al._ Toward $\rm{2D}$ magnets in the ($\rm{MnBi_{2}Te_{4}}$) $\rm{({Bi_{2}Te_{3}})_{n}}$ bulk crystal.  _Advanced Materials_ 32, 2001815 (2020).
  * (15) Hu, C. _et al._ Realization of an intrinsic ferromagnetic topological state in $\rm{MnBi_{8}Te_{13}}$.  _Science advances_ 6, eaba4275 (2020).
  * (16) Wu, J. _et al._ Natural van der waals heterostructural single crystals with both magnetic and topological properties.  _Science advances_ 5, eaax9989 (2019).
  * (17) Roychowdhury, S. _et al._ Giant topological hall effect in the noncollinear phase of two-dimensional antiferromagnetic topological insulator $\rm{MnBi_{4}Te_{7}}$.  _Chemistry of Materials_ 33, 8343–8350 (2021).
  * (18) Hu, C. _et al._ A van der waals antiferromagnetic topological insulator with weak interlayer magnetic coupling.  _Nature communications_ 11, 1–8 (2020).
  * (19) Xie, H. _et al._ The mechanism exploration for zero-field ferromagnetism in intrinsic topological insulator $\rm{MnBi_{2}Te_{4}}$ by $\rm{Bi_{2}Te_{3}}$ intercalations.  _Applied Physics Letters_ 116, 221902 (2020).
  * (20) Vidal, R. _et al._ Orbital complexity in intrinsic magnetic topological insulators $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$.  _Physical Review Letters_ 126, 176403 (2021).
  * (21) Ding, L. _et al._ Crystal and magnetic structures of magnetic topological insulators $\rm{MnBi_{2}Te_{4}}$ and $\rm{MnBi_{4}Te_{7}}$.  _Physical Review B_ 101, 020412 (2020).
  * (22) Yan, J.-Q. _et al._ Crystal growth and magnetic structure of $\rm{MnBi_{2}Te_{4}}$.  _Phys. Rev. Mater._ 3, 064202 (2019).
  * (23) Yan, J.-Q. _et al._ A-type antiferromagnetic order in $\rm{MnBi_{4}Te_{7}}$ and $\rm{MnBi_{6}Te_{10}}$ single crystals.  _Physical Review Materials_ 4, 054202 (2020).
  * (24) Hu, C. _et al._ Tuning magnetism and band topology through antisite defects in $\rm{Sb}$-doped $\rm{MnBi_{4}Te_{7}}$.  _Physical Review B_ 104, 054422 (2021).
  * (25) Shi, M. _et al._ Magnetic and transport properties in the magnetic topological insulators $\rm{MnBi_{2}Te_{4}}$ $\rm{({Bi_{2}Te_{3}})_{n}}$(n= 1, 2).  _Physical Review B_ 100, 155144 (2019).
  * (26) Chen, B. _et al._ Coexistence of ferromagnetism and topology by charge carrier engineering in the intrinsic magnetic topological insulator $\rm{MnBi_{4}Te_{7}}$.  _Physical Review B_ 104, 075134 (2021).
  * (27) Tan, A. _et al._ Metamagnetism of weakly coupled antiferromagnetic topological insulators.  _Physical Review Letters_ 124, 197201 (2020).
  * (28) Vidal, R. C. _et al._ Topological electronic structure and intrinsic magnetization in $\rm{MnBi_{4}Te_{7}}$: a $\rm{Bi_{2}Te_{3}}$ derivative with a periodic $\rm{Mn}$ sublattice.  _Physical Review X_ 9, 041065 (2019).
  * (29) Xu, L. _et al._ Persistent surface states with diminishing gap in $\rm{MnBi_{2}Te_{4}}$/$\rm{Bi_{2}Te_{3}}$ superlattice antiferromagnetic topological insulator.  _Science Bulletin_ 65, 2086–2093 (2020).
  * (30) Yang, S. _et al._ Odd-even layer-number effect and layer-dependent magnetic phase diagrams in $\rm{MnBi_{2}Te_{4}}$.  _Phys. Rev. X_ 11, 011003 (2021).
  * (31) Huang, B. _et al._ Layer-dependent ferromagnetism in a van der waals crystal down to the monolayer limit.  _Nature_ 546, 270–273 (2017).
  * (32) Zang, Z. _et al._ Layer-number-dependent antiferromagnetic and ferromagnetic behavior in $\rm{MnSb_{2}Te_{4}}$.  _Physical Review Letters_ 128, 017201 (2022).
  * (33) Lai, Y., Ke, L., Yan, J., McDonald, R. D. & McQueeney, R. J.  Defect-driven ferrimagnetism and hidden magnetization in $\rm{MnBi_{2}Te_{4}}$.  _Physical Review B_ 103, 184429 (2021).
  * (34) Aliev, Z. S. _et al._ Novel ternary layered manganese bismuth tellurides of the $\rm{MnTe-{Bi_{2}Te_{3}}}$ system: Synthesis and crystal structure.  _Journal of Alloys and Compounds_ 789, 443–450 (2019).
  * (35) Huang, Y. _et al._ Universal mechanical exfoliation of large-area $\rm{2D}$ crystals.  _Nature communications_ 11, 1–9 (2020).
  * (36) Wu, X. _et al._ Distinct topological surface states on the two terminations of $\rm{MnBi_{4}Te_{7}}$.  _Physical Review X_ 10, 031013 (2020).
  * (37) Shao, J. _et al._ Pressure-tuned intralayer exchange in superlattice-like $\rm{MnBi_{2}Te_{4}}$/ $\rm{({Bi_{2}Te_{3}})_{n}}$ topological insulators.  _Nano Letters_ 21, 5874–5880 (2021).
  * (38) Xie, H. _et al._ Charge carrier mediation and ferromagnetism induced in $\rm{MnBi_{6}Te_{10}}$ magnetic topological insulators by antimony doping.  _Journal of Physics D: Applied Physics_ 55, 104002 (2021).
  * (39) Hor, Y. _et al._ Development of ferromagnetism in the doped topological insulator $\rm{Bi_{2-x}Mn_{x}Te_{3}}$.  _Physical Review B_ 81, 195203 (2010).
  * (40) Henk, J. _et al._ Topological character and magnetism of the dirac state in $\rm{Mn}$-doped $\rm{Bi_{2}Te_{3}}$.  _Physical Review Letters_ 109, 076801 (2012).
  * (41) Zimmermann, S. _et al._ Spin dynamics and magnetic interactions of $\rm{Mn}$ dopants in the topological insulator $\rm{Bi_{2}Te_{3}}$.  _Physical Review B_ 94, 125205 (2016).
  * (42) Liu, Y. _et al._ Site mixing for engineering magnetic topological insulators.  _Physical Review X_ 11, 021033 (2021).
  * (43) Chen, B. _et al._ All-oxide–based synthetic antiferromagnets exhibiting layer-resolved magnetization reversal.  _Science_ 357, 191–194 (2017).


# The Impact of Evolutionary Computation on Robotic Design:
A Case Study with an Underactuated Hand Exoskeleton

Baris Akbas1, Huseyin Taner Yuksel1, Aleyna Soylemez1, Mazhar Eid Zyada2,
Mine Sarac2, , and Fabio Stroppa1 This work is funded by TUBİTAK project
number 123M690 and partially funded by TUBİTAK project number 121C145 and
121C147.1Computer Engineering, Kadir Has University, İstanbul, Turkey. E-mail:
<EMAIL_ADDRESS>Engineering, Kadir Has University,
İstanbul, Turkey.

###### Abstract

Robotic exoskeletons can enhance human strength and aid people with physical
disabilities. However, designing them to ensure safety and optimal performance
presents significant challenges. Developing exoskeletons should incorporate
specific optimization algorithms to find the best design. This study
investigates the potential of Evolutionary Computation (EC) methods in robotic
design optimization, with an underactuated hand exoskeleton (U-HEx) used as a
case study. We propose improving the performance and usability of the U-HEx
design, which was initially optimized using a naive brute-force approach, by
integrating EC techniques such as Genetic Algorithm and Big Bang-Big Crunch
Algorithm. Comparative analysis revealed that EC methods consistently yield
more precise and optimal solutions than brute force in a significantly shorter
time. This allowed us to improve the optimization by increasing the number of
variables in the design, which was impossible with naive methods. The results
show significant improvements in terms of the torque magnitude the device
transfers to the user, enhancing its efficiency. These findings underline the
importance of performing proper optimization while designing exoskeletons, as
well as providing a significant improvement to this specific robotic design.

##  I Introduction

Exoskeleton robotic devices are often used to augment users’ strength and
endurance during physically demanding tasks [1, 2, 3], to allow users control
a secondary robotic device during teleoperation scenarios [4], or to aid
limited movements for patients with neurological and physical disabilities
[5]. Depending on the application, such exoskeletons can be designed for the
whole body [6] or for specific body locations such as arms [7], legs [8],
wrists [9], or hands [10]. Regardless of the application or the body location,
exoskeletons are very challenging to be designed, implemented, and controlled
[11, 12]. Safety is the primary and most important issue: exoskeleton joints
(i) must align perfectly with anatomical joints to avoid potential harm, (ii)
should work effectively within the workspace of human anatomical joints, and
(iii) should allow the exoskeleton to follow users’ behavior without creating
discomfort. Finally, these devices –- especially assistive ones –- should be
as compact and lightweight as possible to enhance wearability.

Figure 1: U-HEx: (a) A user wearing the first prototype of U-HEx from the
state-of-the-art, which is bulky and cumbersome [10]. Note that the picture
only shows the device with a single finger. (b) Schema of the comparison
between the approaches in terms of run time and optimality of the solution.

These challenges can be overcome by designing exoskeletons that can reach high
output forces and feature effective power transmission despite using small,
lightweight actuators. Thus, the design process of such exoskeletons should be
integrated with various optimization algorithms; that is, the search for the
best element within a set of alternatives based on specific criteria.
Optimization is a common tool for solving engineering problems [13, 14, 15].
While the most conventional strategies focus on numerical and calculus-based
methods [16], they might not be the best solution for engineering designs due
to their properties such as non-discrete domains, non-differentiability,
multi-modality, discontinuity, reliability, and robustness. Alternatively, the
nature-inspired methods of Evolutionary Computation (EC) appear to be a common
and effective way to deal with engineering optimization problems [17] – and
often with exoskeleton design [12]. Unfortunately, the integration between
design and its optimization is not always straightforward. Roboticists, who
often have mechanical/mechatronics backgrounds, might not know the latest
trends or trade-offs in optimization. This is even more exacerbated by the
lack of systematic studies in the literature on the impact of different
optimization techniques for robotic design.

In this work, we attempt to fill the gap in the literature by raising
awareness in the robotic community while displaying the impact of EC on
robotic designs. We provide a systematical analysis of its performance on a
poorly optimized design from the literature [10]. Compared to a naive-
deterministic search approach (i.e., iteratively exploring each possible
combination of design parameters – “brute-force”), we hypothesize that EC
techniques (i) might converge to a more precise and better solution (H1) and
(ii) will reduce the convergence time despite exploring border search spaces
(H2). In addition, using various EC algorithms might offer different solutions
at different convergence times (H3).

Specifically, we applied the brute-force optimization approach and two EC
algorithms on the design of an Underactuated Hand Exoskeleton (U-HEx [10],
shown in Fig. 1 (a)), firstly on the same (limited) search space already
explored in previous works, and then on a wider search space – which could not
be achieved with brute force. We specifically chose U-HEx due to its complex
mechanism. Unlike conventional serial-link robots, its interconnected
kinematics model makes it challenging to predict the contribution of each link
length to the achieved range of motion. We believe that this complexity
highlights the differences between optimization methods and helps us formulate
a clearer and better discussion. The comparison of the final designs will
analyze both the optimality of the solution (i.e., effective force
transmission) and the run time, as shown in Fig. 1 (b).

##  II Background

###  II-A The Underactuated Hand Exoskeleton (U-HEx)

Fig. 1 (a) shows U-HEx – a wearable robotic device for the hand to
rehabilitate stroke survivors through physical therapy [10]. U-HEx is designed
with a single actuator to control two finger joints through underactuation
[18]. With no external forces on the phalanges, the actuator opens and closes
the finger naturally. As the user interacts with physical objects, the
underactuated mechanism modifies the transmitted forces to each finger joint.
Ultimately, U-HEx can automatically adjust its behavior based on the
interaction forces – allowing users to grasp objects with different shapes and
sizes using a single actuator with no prior mechanical or control adjustments
[19].

Furthermore, U-HEx promotes enhanced safety in multiple aspects of human-robot
interaction. Firstly, the finger phalanges are considered as a part of the
kinematic chain– such that the device is self-adaptable to a predefined range
of hand sizes. Secondly, the underactuated kinematics inherently decouples the
mechanical joints from the anatomical ones, so there is no need for
calibration. The forces acting on the finger joints are transferred through
the complicated design of mechanical links of U-HEx from the single actuator.
Therefore, deciding the mechanical link lengths is of the utmost importance to
ensure effective force transmission despite its complex kinematics chain.

###  II-B Optimization and Evolutionary Computation

Optimization is the mathematical process of searching for a set of decision
variables to minimize or maximize one or more specific objective functions
while satisfying certain constraints. Optimization problems can be solved
using methods that systematically and efficiently create and compare solutions
to find the best outcome – namely, Local Search techniques. These methods can
be exact, providing the precise optimal solution based on direct or gradient-
based approaches [16], or they can be approximate, yielding sub-optimal
solutions that are acceptable approximations of the global optimum. In
engineering design problems, approximate methods are often favored due to the
complexity of objective functions [20], the preference for robust and reliable
solutions over global ones [21], and the presence of uncertainty in the search
space [22].

Evolutionary Computation (EC), a sub-field of soft computing, offers popular
approximate methods for engineering [12]. Inspired by natural selection, EC
techniques generate a population of potential solutions and evolve them toward
the optimal solution using different metaheuristics (i.e., based on genetic
recombination, a storm of birds foraging for food, a community of bugs
building a colony, etc.). These population-based techniques allow for
parallelized search and the retrieval of multiple optimal solutions,
particularly in multi-modal or conflicting multi-objective problems. Engineers
find EC methods beneficial because they can effectively deal with objective
functions defined implicitly through simulations and not rely on a specific
mathematical model. This flexibility allows the same algorithmic
implementation to be applied to different problems.

##  III Exoskeleton Design Optimization

###  III-A U-HEx Kinematics

Figure 2: Kinematic model of U-HEx depicting all link lengths, passive joints
(rotational and linear), and the active linear actuator. The colored links
represent the decision variables of the optimization problem (i.e., the links
to be optimized): in yellow, the ones that were used in its original design
[10]; in green, the additional three that were included in the current work
thanks to evolutionary computation. The rest of the links are set to fixed
values to ensure the kinematic chain is closed. Only joint $O$ is actuated.

Fig. 2 shows the kinematic model for U-HEx with important points depicted with
letters as detailed in previous works [10]. Each gray dot represents a passive
rotational joint (or anatomical finger joints) while empty circles ($H$, $E$,
and $C$) are fixed points along rigid links with $90^{o}$ fixed angle. The
system has only one linear actuator between the points $O$ and $A$. The
exoskeleton is fixed on the hand from points $K$ and $O$ with variable lengths
in the $x$ and $y$ directions. Finally, the exoskeleton is attached to the
finger phalanges from points $I$ and $J$ – with passive linear joint sliders
represented as $c_{1}$ and $c_{2}$, respectively.

Defining the closed-loop kinematic chains with the letters depicted in Fig. 2,
inverse kinematics are computed via numerical methods to compute $8$ unknown
variables ($l_{OA}$, $q_{O}$, $q_{A}$, $q_{B}$, $q_{G}$, $q_{D}$, $c_{1}$,
$c_{2}$) for given finger pose ($q_{MCP},q_{PIP}$ – from fully open to fully
closed) and the set of link lengths. Once the kinematics are computed, the
distribution from the unit actuator forces (1 N) to the torques around the
finger joints ($\tau_{MCP},\tau_{PIP}$) can be obtained either through static
equations or the Jacobian of the system.

###  III-B Link-Length Optimization Problem

The mechanical design of U-HEx should be optimized by searching for the set of
link lengths (decision variables) that maximizes its force transmission
(objective function) – i.e., the amount of force transferred to the finger
joints ($\tau_{i}$), as formulated in Eqn. (1). The problem is also subjected
to the following physical constraints: (i) U-HEx is connected to the user’s
fingers through passive linear sliders ($c1$ and $c2$ in Fig. 2), whose
movements are limited by the user’s finger size and (ii) the ratio between the
torques exerted on two finger joints must be within a reasonable range at
different orientations of the finger (between $0.05$ and $20$).

maximize $\displaystyle\sqrt{\tau_{MCP}+\tau_{PIP}}$ (1)

While the hand exoskeleton has many link lengths that needed to be set, the
previous study only optimized the important set of link lengths that were
found through sensitivity analysis [10]. These important link lengths are
depicted with yellow-lined links in Fig. 2 (i.e., $\overline{BC}$,
$\overline{CD}$, $\overline{DE}$, $\overline{EF}$, $\overline{FG}$, and
$\overline{GH}$). Green-lined links indicate the additional decision variables
that were included in the current study thanks to EC ($\overline{BK}$,
$\overline{CI}$, and $\overline{EJ}$). Gray-lined links indicate the link
lengths that are always kept constant.

All optimization methods are performed using MATLAB script. For a decided set
of link lengths to be tried, a Simulink model is executed with a fixed-step
solver to compute (i) the inverse kinematics through numerical methods and
(ii) the statics through analytical methods as the finger joints are iterated
from fully open (0 deg each) to fully closed (80 deg for MCP and 90 deg for
PIP). Once the Simulink file is terminated, we first check the constraints on
the passive sliders ($c_{1}$ and $c_{2}$) and then the ratio between the
transmitted torques ($\frac{\tau_{MCP}}{\tau_{PIP}}$). If the given set of
link lengths satisfies these constraints, the objective function is computed
for the finger pose fully closed.

###  III-C Optimization Methods

####  III-C1 Brute Force (BF)

The optimization method originally implemented to design U-HEx aimed at
evaluating every solution in the search space [10]. This is a brute-force
method to solve an optimization problem; therefore, it is highly inefficient
regarding both time and computational resources. Since the decision variables
are lengths (i.e., measured in millimeters), their domain is considered almost
continuous, making the search space too wide (theoretically infinite) for
being treated with brute force – even when the variables are bounded to
specific lower and upper limits. Due to the enormity of this search space, the
designers were forced to introduce the following limitations:

  * •

discretize the continuous domain of the decision variables by sampling with a
fixed interval/step;

  * •

increase the step between two contiguous discretized values for each decision
variable (i.e., reduce the precision in millimeters); and

  * •

reduce the number of decision variables, fixing the values of specific link
lengths – which were selected through sensitivity analysis to identify the
ones that do not significantly affect the output performance.

The time complexity of such an algorithm is $\mathcal{O}(n^{d})$, where $d$ is
the number of decision variables, and $n$ is the cardinality of their domain
(assumed to be the same for each variable or equal to the variable with the
highest cardinality). When U-HEx was originally designed, this execution took
approximately three days, during which the system crashed several times due to
excessive processing (tested on a 2014 machine with a 2.40 GHz CPU and 16 GB
RAM).

####  III-C2 Genetic Algorithms (GAs)

GAs are the most popular methods of EC techniques as they directly implement
the process of natural selection and survival of the fittest [23, 24]. About
Algorithm 1, they (i) generate a population of random solutions $P$ within the
search space of the problem, (ii) assign a fitness value to each solution
based on the objective function, and (iii) generate new solutions $Q$ by
mixing the values of the ones in the current population (i.e., a process named
crossover).

By allowing only the most-fitting solutions $M$ to perform crossover and be
preserved in the next generations, GAs evolve their population to converge to
the optimum of the problem. Like in biology, the operation of crossover
exploits the features of good solutions (parents) to produce similar new
solutions (offspring) and speed up convergence to an optimum; however, there
is no guarantee for the optimum to be global rather than local. Therefore, GAs
implement an additional operator inspired by genetic mutation, randomly
modifying values of a newly generated solution to favor search-space
exploration and escape local optima.

Besides the common advantages of population-based methods (see Sec. II-B), GAs
are easy to implement and very efficient to converge. On the other hand, their
many genetic operators come with many parameters and different types, and
their fine-tuning might be non-trivial and primarily based on trial and error.
Furthermore, since GAs are iterative stochastic methods, they might be
inefficient for real-time problem solving – which is not the case for this
study.

input : Population size n, number of generations g

output : The most fitting solution $\textnormal{{P}}(1)$

1

2begin

3

4
$\textnormal{{P}}\leftarrow\textnormal{{randomInitialization}}(\textnormal{{n}})$;

5 $\textnormal{{P}}\leftarrow\textnormal{{evaluation}}(\textnormal{{P}})$;

6 for  _$i\in[1,\textnormal{{g}}]$_ do

7 $\textnormal{{M}}\leftarrow\textnormal{{selection}}(\textnormal{{P}})$;

8 $\textnormal{{Q}}\leftarrow\textnormal{{variation}}(\textnormal{{M}})$;

9 $\textnormal{{Q}}\leftarrow\textnormal{{evaluation}}(\textnormal{{Q}})$;

10
$\textnormal{{P}}\leftarrow\textnormal{{survival}}(\textnormal{{P}},\textnormal{{Q}})$;

11

12

13 return _P_ ;

14

Algorithm 1 Genetic Algorithm

input : Population size n, number of generations g

output : The most fitting solution $\textnormal{{P}}(1)$

1

2begin

3

4
$\textnormal{{P}}\leftarrow\textnormal{{randomInitialization}}(\textnormal{{n}})$;

5 for  _$i\in[1,\textnormal{{g}}]$_ do

6 if  _$i\neq 1$_ then

7 $\textnormal{{P}}\leftarrow\textnormal{{bang}}(\textnormal{{cm}},i)$;

8

9 $\textnormal{{P}}\leftarrow\textnormal{{evaluation}}(\textnormal{{P}})$;

10 $\textnormal{{cm}}\leftarrow\textnormal{{crunch}}(\textnormal{{P}})$;

11

12

13 return _P_ ;

14

Algorithm 2 Big Bang-Big Crunch Algorithm

####  III-C3 Big Bang-Big Crunch Algorithm (BB-BC)

BB-BC [25] is inspired by the evolution of the universe through two phases of
explosion and implosion: (i) energy dissipation producing disorder and
randomness, and (ii) randomness drawn back into a (different) order. With
reference to Algorithm 2, BB-BC creates an initial random population $P$
uniformly spread throughout the search space (the explosion, or big bang),
evaluates them, and collects them into their center of mass $cm$ (the
implosion, or big crunch). These two phases are repeated throughout the
execution, spreading new solutions closer to the center of mass as the number
of iterations increases. Re-iterating this procedure leads the center of mass
to converge to the optimal solution of the problem. BB-BC is known to
outperform GAs in terms of convergence speed; thus, we considered applying it
to our problem due to the high run time of the objective function.

##  IV Results of the Comparisons

We conducted two main experiments on a computer with 16-core 5.4 GHz CPU and
64 GB RAM using (i) three optimization methods (BF, GA, and BB-BC) following
the original optimization settings for U-HEx design [10] within the original
search space and (ii) two evolutionary optimization methods (GA and BB-BC) in
a wider search space. These values have been chosen empirically through
observations and sensitivity analysis on the feasibility of the solution
retrieved (i.e., the solutions are infeasible outside those bounds). Their
results will be compared regarding the optimality of the solution and run
time.

Tab. I lists each algorithm’s parameters and respective values. Most of them
did not require any pre-evaluation and were empirically set to a value or a
type. For example, the number of generations/iterations was set to $50$ after
observing that both algorithms converge earlier (around 35th iteration for GA,
and 20th for BB-BC, after which the most fitting solution does not improve
more than $0.5$ Nm). The probability of performing mutation in GA and the
population size require further investigation through preliminary tests to be
determined to increase the efficacy of each EC method.

TABLE I: Parameters of GA and BB-BC

PARAMETER | PRE. TEST | EXPERIMENT VALUE
---|---|---
Max Num. Generations (GA, BB-BC) |  | $50$
Population Size (GA, BB-BC) |  $150$, $300$ | $300$
Selection Type (GA) |  | Binary Tournament [26]
Crossover Type (GA) |  | blx-$\alpha$ ($\alpha=0.5$) [27]
Crossover Probability (GA) |  | $1.0$
Mutation Type (GA) |  | Polynomial [28]
Mutation Probability (GA) |  $0.2$, $0.4$, $0.6$ | $0.2$
Survival Type (GA) |  | Elitist ($\mu+\lambda$ scheme) [29]
Crunch Method (BB-BC) |  | Best Fit [30]
Constraint Handling (GA, BB-BC) |  | Deb’s method [31]

###  IV-A Preliminary Tests on Optimization Methods Parameters

There is no strict rule for setting the probability of performing mutation of
the GA ($PoM$), even though the literature suggests keeping this value low to
favor the exploitation of a good solution [32]. We executed the GA ten times
for each different value ($3_{PoM}=0.2$, $0.4$, $0.6$) on a population of 150
solutions and observed the respective designs. The results of a one-way
Analysis of Variation (ANOVA) indicate that there are no statistically
significant differences among them ($F(2,8)=0.913$, $p=0.419$,
$\eta^{2}=0.092$).

Regarding the population size ($P$), larger values correspond to higher
search-space exploration but negatively affect the run time as the number of
evaluations increases linearly. We executed the GA and the BB-BC independently
ten times (each) for different population sizes ($2_{P}=150$, $300$) and
observed their effect on the solutions. For GA, our t-test results indicate
that population 300 obtains statistically significantly higher forces than 150
($p=0.012$) against a higher run time ($p=0.004$). In contrast, for BB-BC, our
t-test results indicate that population 300 obtains a statistically
significant difference from 150 in terms of run time ($p=0.005$) but not the
obtained forces ($p=0.162$).

Based on these results, we performed the main experiments with a fixed value
for ${PoM=0.2}$ and the population size ${P=300}$, as also summarized in Table
I.

TABLE II: Link Lengths Bounds (Experiment 1) $\overline{BC}$ | $\overline{CD}$ | $\overline{DE}$ | $\overline{EF}$ | $\overline{FG}$ | $\overline{GH}$
---|---|---|---|---|---
$38\rightarrow 60$ | $10\rightarrow 30$ | $15\rightarrow 51$ | $15\rightarrow 51$ | $20\rightarrow 56$ | $64\rightarrow 100$
TABLE III: Link Lengths Bounds (Experiment 2) $\overline{BK}$ | $\overline{CI}$ | $\overline{EJ}$
---|---|---
$20\rightarrow 50$ | $10\rightarrow 17$ | $20\rightarrow 50$
Figure 3: Comparison between results of Experiments 1 and 2. Plots include
median, interquartile range, and outliers. Details are zoomed in to appreciate
variance amongst different conditions.

###  IV-B Experiment 1: Comparison with Previous Work

We evaluated the impact of EC on U-HEx design by comparing the optimality of
their retrieved solutions against the original BF design [10]. For a valid
comparison, the search space of the chosen 6 decision variables is kept
constant as in our previous work, summarized in Table II. To compare the run
time as well, we re-executed BF. Since BF is expected to provide the same
output at every run, we executed it only once, whereas we executed both GA and
BB-BC twenty times each. At each execution, we recorded the set of optimized
link lengths (i.e., values of decision variable for the most fitting solution)
and compared the optimization methods in terms of (i) the optimality measure
(i.e., torque magnitude) and (ii) the run time until convergence (i.e., when
the improvement in the most fitting solution’s value was smaller than $0.5$
Nm).

TABLE IV: Numerical results of Experiment 1 (6 DVs) | Optimality (Nm) | Run Time (s)
---|---|---
Brute Force | $31.53\pm 0.00$ | $94830.49\pm 0.00$
GA | $31.63\pm 0.01$ | $860.51\pm 168.25$
BB-BC | $31.65\pm 0.04$ | $312.50\pm 204.63$

##### Optimality

Fig. 3 up-left and Table IV show the optimized joint torques for all
optimization methods. The findings of one-way ANOVA indicate these methods to
be statistically significantly different than each other ($F(2,8)=280.515$,
$p<0.001$, $\eta^{2}=0.986$). We then performed a post-hoc analysis with the
Bonferroni test: the optimal solution obtained with BF is found to be
significantly worse than GA ($p<0.001$) and than BB-BC ($p<0.001$) but not
between GA and BB-BC ($p=0.269$).

##### Run Time

Fig. 3 down-left and Table IV show the run time until convergence for all
optimization methods. The findings of one-way ANOVA indicate these methods to
be statistically significantly different than each other
($F(2,18)=17882238.147,p<0.001,\eta^{2}=1.00$). Further post-hoc (Bonferroni)
analysis shows that the run time with BF is significantly higher than GA
($p<0.001$) and than BB-BC ($p<0.001$), while the run time until convergence
is significantly higher with GA than BB-BC ($p<0.001$).

###  IV-C Experiment 2: Inclusion of More Decision Variables

The previous experiment shows that EC methods improve the optimum performance
measures and the run time compared to a naive BF. With the run time decreasing
more than 10 times, we can now include other decision variables, which were
constant in the previous experiment. We included three additional decision
variables (DVs) to the algorithms: BK, EJ, and CI, with the domains reported
in Table II. To emphasize the impact of this comparison, we performed a two-
way ANOVA with factors defined as EC ($2_{EC}$ = GA, BB-BC) and DVs ($2_{DVs}$
= 6 DVs, 9 DVs).

TABLE V: Numerical results of Experiment 2 (9 DVs) | Optimality (Nm) | Run Time (s)
---|---|---
GA | $56.20\pm 0.41$ | $1208.59\pm 78.31$
BB-BC | $55.99\pm 0.54$ | $802.50\pm 298.03$
TABLE VI: Best Design retrieved with Brute Force vs Evolutionary Algorithms | Link Lengths (mm) | Magnitude (Nm)
---|---|---
| $\overline{BC}$ | $\overline{CD}$ | $\overline{DE}$ | $\overline{EF}$ | $\overline{FG}$ | $\overline{GH}$ | $\overline{BK}$ | $\overline{CI}$ | $\overline{EJ}$
Brute Force | 58.00 | 10.00 | 15.00 | 51.00 | 56.00 | 100.00 | 35.00 | 16.00 | 37.00 | 31.53
Evolutionary | 60.00 | 10.00 | 15.00 | 51.00 | 56.00 | 91.37 | 48.50 | 10.98 | 36.54 | 56.86

##### Optimality

Fig. 3 up-right and Tab. V show the results of the experiment with 9 DVs in
terms of torque magnitudes. We found statistical significance between
different DVs ($F(1,36)=46655.129$, $p<0.001$, $\eta^{2}=0.999$), but not
between EC ($F(1,36)=0.739$, $p=0.396$, $\eta^{2}=0.02$) or interactions
($F(1,36)=1.131$, $p=0.295$, $\eta^{2}=0.03$).

##### Run Time

Fig. 3 up, right and Table V show the results of the second experiment with 9
DVs in terms of the convergence run time. We found statistical significance
for the main factors of DVs ($F(1,36)=38.281$, $p<0.001$, $\eta^{2}=0.515$)
and for EC ($F(1,36)=49.615$, $p<0.001$, $\eta^{2}=0.58$), but not between
interactions ($F(1,36)=1.098$, $p=0.302$, $\eta^{2}=0.30$).

##  V Discussions

Our main motivation was to systematically compare EC algorithms to naive
optimization methods like BF from the perspective of an engineering design
problem – e.g., the mechanical design of a robotic device (U-HEx). With
Experiment 1, we compared a previously implemented method with two EC methods
and observed that EC methods are statistically significantly better than BF –
with the given domain and restrictions regarding the optimality of the
obtained solution [10]. Similarly, the run time recorded with BF is at least
ten times higher than EC methods. Therefore, our first two hypotheses (H1) EC
provides better and more optimal solutions than practical BF and (H2) EC
methods provide an optimal solution significantly faster than BF hold true.

The most relevant practical limitation of BF is its run time. Although running
BF was computationally possible (with 26 hours run time), three main
restrictions were originally made to permit it: discretizing the continuous
domain, increasing the step between contiguous values, and reducing the number
of decision variables. Based on our results, the first two restrictions
hindered more fitting solutions. It is evident that, with the same search
space, BF would retrieve the optimal solution, outperforming any optimization
method in terms of optimality. However, in practice, even with a discrete
space and a step of 1 between solutions
($905\mathpunct{\raisebox{0.0pt}{,}}219\mathpunct{\raisebox{0.0pt}{,}}763$
combinations), we estimate a run time of $4\mathpunct{\raisebox{0.0pt}{,}}190$
days on the same machine – with no guarantee of retrieving the same
(sub)optimal EC solution reported in Table VI and shown in Fig. 4, which
featured a (pseudo)-continuous search space. In other words, it is possible
that the BF could yield solutions that are not statistically different than
EC, but EC’s superiority in run time would hold true regardless.

Figure 4: CAD models of U-HEx with the lengths retrieved from (a) BF, similar
to the original prototype [10] and (b) EC – specifically BB-BC.

The third restriction was to reduce the number of decision variables from the
search space and keep their values constant. Particularly, the sensitivity
analysis performed in the original work [10] led designers to remove some
decision variables from the problem. Thanks to EC’s faster run time, we were
able to include three further decision variables in Experiment 2 (Sec. IV-C).
Our results show that adding more decision variables slightly increases the
convergence run time to favor the optimality of the retrieved solution.

Previously, designers of U-HEx used sensitivity analysis to identify the most
impactful link lengths as decision variables to make the execution more
effective and not waste computation time [10]. Interestingly, Tables VI and II
show that four of these six link lengths ($\overline{CD}$, $\overline{DE}$,
$\overline{EF}$, and $\overline{FG}$) remain the same at the mechanical limits
imposed by the kinematic chain. Thus, pre-preparing the optimization through
search spaces and limited link lengths might require extra attention. Even
though with high dimensional spaces, sensitivity analysis can still be useful
with EC, it was unnecessary while using EC methods in this study.

Our third hypothesis (H3) various EC algorithms might offer different
solutions at different convergence times was also confirmed. BB-BC converged
significantly faster than GA – which is in line with the claims of the
literature [25] even though we observed no statistically significant
differences between GA and BB-BC in optimality. We also observed that
increasing the population size of GA significantly improves the optimality of
the retrieved solution – indicating that larger population sizes might lead to
better designs than the one reported in Table VI and shown in Fig. 4. This can
also be true for BB-BC, even though we did not find a significant difference
in changing the population size (possibly, the observation value was not large
enough). However, increasing the population size also increases the overall
run time.

Lastly, we would like to emphasize that we do not claim to have proven the
superiority of optimization methods against simple and brute enumeration (BF)
– this is a well-known advantage of numerical optimization. Instead, our main
motivation is to compare different EC methods and provide systematic evidence
for the impact of alternative approaches through a hands-on design case of a
robotic exoskeleton device (U-HEx). With the enlightenment from our findings,
we invite young roboticists, researchers, and designers to be mindful of such
alternative methods and choose efficient and effective methods to reach
optimality for their design. After all, due to this negligence, U-HEx existed
in a non-optimal shape for more than a decade.

##  VI Conclusions

In this work, we presented a comparative study to highlight the impact of EC
on robotic design. Specifically, we re-optimized the design parameters for
U-HEx, an underactuated hand exoskeleton with a numerical, complex kinematic
structure – previously optimized with a naive brute-force method. We showed
that EC allowed the device to be further optimized by adding further decision
variables. Ultimately, increasing the optimality of the device might actually
improve the usability and the efficacy of U-HEx during human interactions
during physical rehabilitation therapy for patients with hand disabilities.
The code is available for on EVO Lab’s MathWorks File Exchange
repository111www.mathworks.com/matlabcentral/fileexchange/157446.

In the future, we will compare our results with more optimization methods, and
expand the study with further objectives to optimize U-HEx (e.g., balancing
the forces on each joint to minimize the possibility of hurting the user and
reducing the size of the exoskeleton to promote comfort and portability). This
investigation will require the implementation of specific multi-objective
optimization algorithms with EC and cannot be achieved by BF. We will also
investigate the implications of the proposed improvements during the real
human-robot interaction by manufacturing both designs and studying the user
experience and interaction forces.

## References

  * [1] T. Poliero, M. Lazzaroni, S. Toxiri, C. Di Natali, D. G. Caldwell, and J. Ortiz, “Applicability of an active back-support exoskeleton to carrying activities,” _Frontiers in Robotics and AI_ , vol. 7, p. 579963, 2020.
  * [2] A. Mauri, J. Lettori, G. Fusi, D. Fausti, M. Mor, F. Braghin, G. Legnani, and L. Roveda, “Mechanical and control design of an industrial exoskeleton for advanced human empowering in heavy parts manipulation tasks,” _Robotics_ , vol. 8, no. 3, p. 65, 2019.
  * [3] T. Bützer, O. Lambercy, J. Arata, and R. Gassert, “Fully wearable actuated soft exoskeleton for grasping assistance in everyday activities,” _Soft robotics_ , vol. 8, no. 2, pp. 128–143, 2021.
  * [4] T. Koyama, I. Yamano, K. Takemura, and T. Maeno, “Multi-fingered exoskeleton haptic device using passive force feedback for dexterous teleoperation,” in _IEEE/RSJ International Conference on Intelligent Robots and Systems_ , vol. 3, 2002, pp. 2905–2910.
  * [5] F. Stroppa, C. Loconsole, S. Marcheschi, and A. Frisoli, “A robot-assisted neuro-rehabilitation system for post-stroke patients’ motor skill evaluation with alex exoskeleton,” in _Proceedings of the International Conference on NeuroRehabilitation (ICNR)_ , 2017, pp. 501–505.
  * [6] S. Marcheschi, F. Salsedo, M. Fontana, and M. Bergamasco, “Body extender: Whole body exoskeleton for human power augmentation,” in _IEEE International Conference on Robotics and Automation (ICRA)_ , 2011, pp. 611–616.
  * [7] D. Gijbels, I. Lamers, L. Kerkhofs, G. Alders, E. Knippenberg, and P. Feys, “The armeo spring as training tool to improve upper limb functionality in multiple sclerosis: A pilot study,” _Journal of Neuroengineering and Rehabilitation_ , vol. 8, pp. 1–8, 2011.
  * [8] A. B. Zoss, H. Kazerooni, and A. Chu, “Biomechanical design of the berkeley lower extremity exoskeleton (bleex),” _IEEE/ASME Transactions on Mechatronics_ , vol. 11, no. 2, pp. 128–138, 2006.
  * [9] D. Buongiorno, E. Sotgiu, D. Leonardis, S. Marcheschi, M. Solazzi, and A. Frisoli, “Wres: A novel 3 DoF WRist ExoSkeleton with tendon-driven differential transmission for neuro-rehabilitation and teleoperation,” _IEEE Robotics and Automation Letters (RA-L)_ , vol. 3, no. 3, pp. 2152–2159, 2018.
  * [10] M. Sarac, M. Solazzi, E. Sotgiu, M. Bergamasco, and A. Frisoli, “Design and kinematic optimization of a novel underactuated robotic hand exoskeleton,” _Meccanica_ , vol. 52, pp. 749–761, 2017.
  * [11] M. Sarac, M. Solazzi, and A. Frisoli, “Design requirements of generic hand exoskeletons and survey of hand exoskeletons for rehabilitation, assistive, or haptic use,” _IEEE Transactions on Haptics (ToH)_ , vol. 12, no. 4, pp. 400–413, 2019.
  * [12] F. Stroppa, A. Soylemez, H. T. Yuksel, B. Akbas, and M. Sarac, “Optimizing exoskeleton design with evolutionary computation: An intensive survey,” _Robotics_ , vol. 12, no. 4, p. 106, 2023.
  * [13] R. Sioshansi and A. J. Conejo, _Optimization in Engineering_. Cham: Springer International Publishing, 2017, vol. 120.
  * [14] R. B. Statnikov and J. B. Matusov, _Multicriteria Optimization and Engineering_. Springer Science and Business Media, 2012.
  * [15] J. Andersson, “A survey of multiobjective optimization in engineering design,” _Department of Mechanical Engineering, Linktjping University. Sweden_ , 2000.
  * [16] J.-F. Bonnans, J. C. Gilbert, C. Lemaréchal, and C. A. Sagastizábal, _Numerical Optimization: Theoretical and Practical Aspects_. Springer Science and Business Media, 2006.
  * [17] D. Dumitrescu, B. Lazzerini, L. C. Jain, and A. Dumitrescu, _Evolutionary Computation_. CRC press, 2000.
  * [18] T. Laliberte, L. Birglen, and C. Gosselin, “Underactuation in robotic grasping hands,” _Machine Intelligence and Robotic Control_ , vol. 4, no. 3, pp. 1–11, 2002.
  * [19] A. Buryanov and V. Kotiuk, “Proportions of hand segments,” _Int. J. Morphol_ , pp. 755–758, 2010.
  * [20] H. Norde, F. Patrone, and S. Tijs, “Characterizing properties of approximate solutions for optimization problems,” _Mathematical Social Sciences_ , vol. 40, no. 3, pp. 297–311, 2000.
  * [21] Y. Nomaguchi, K. Kawakami, K. Fujita, Y. Kishita, K. Hara, and M. Uwasu, “Robust design of system of systems using uncertainty assessment based on lattice point approach: Case study of distributed generation system design in a japanese dormitory town,” _International Journal of Automation Technology_ , vol. 10, no. 5, pp. 678–689, 2016.
  * [22] B. Dizangian and M. Ghasemi, “Reliability-based design optimization of complex functions using self-adaptive particle swarm optimization method,” _International Journal of Optimization in Civil Engineering_ , vol. 5, no. 2, pp. 151–165, 2015.
  * [23] D. E. Goldberg, _Genetic Algorithms in Search, Optimization, and Machine Learning_. Addison-Wesley, 1989.
  * [24] ——, _Real-Coded Genetic Algorithms, Virtual Alphabets and Blocking_ , 1991, vol. 5, no. 2.
  * [25] O. K. Erol and I. Eksin, “A new optimization method: Big Bang–Big Crunch,” _Advances in Engineering Software_ , vol. 37, no. 2, pp. 106–111, 2006.
  * [26] D. E. Goldberg and K. Deb, “A comparative analysis of selection schemes used in genetic algorithms,” in _Foundations of genetic algorithms_. Elsevier, 1991, vol. 1, pp. 69–93.
  * [27] L. J. Eshelman and J. D. Schaffer, “Real-coded genetic algorithms and interval-schemata,” in _Foundations of genetic algorithms_. Elsevier, 1993, vol. 2, pp. 187–202.
  * [28] K. Deb, M. Goyal _et al._ , “A combined genetic adaptive search (geneas) for engineering design,” _Computer Science and informatics_ , vol. 26, pp. 30–45, 1996.
  * [29] H.-G. Beyer and H.-P. Schwefel, “Evolution strategies–a comprehensive introduction,” _Natural computing_ , vol. 1, pp. 3–52, 2002.
  * [30] H. M. Genç, I. Eksin, and O. K. Erol, “Big bang-big crunch optimization algorithm hybridized with local directional moves and application to target motion analysis problem,” in _2010 IEEE International Conference on Systems, Man and Cybernetics_. IEEE, 2010, pp. 881–887.
  * [31] C. A. C. Coello, “Theoretical and numerical constraint-handling techniques used with evolutionary algorithms: A survey of the state of the art,” _Computer Methods in Applied Mechanics and Engineering_ , vol. 191, no. 11-12, pp. 1245–1287, 2002.
  * [32] D. E. Goldberg, J. Richardson _et al._ , “Genetic algorithms with sharing for multimodal function optimization,” in _Genetic Algorithms and Their Applications: Proceedings of the Second International Conference on Genetic Algorithms_ , vol. 4149, 1987.