∎

11institutetext: Department of Physics, Faculty of Science, Okayama
University, Okayama 700-8530, Japan

#  Fidelity-mediated analysis of the transverse-field $XY$ chain with the
long-range interactions: Anisotropy-driven multi-criticality

Yoshihiro Nishiyama

(Received: date / Accepted: date)

###### Abstract

The transverse-field $XY$ chain with the long-range interactions was
investigated by means of the exact-diagonalization method. The algebraic decay
rate $\sigma$ of the long-range interaction is related to the effective
dimensionality $D(\sigma)$, which governs the criticality of the transverse-
field-driven phase transition at $H=H_{c}$. According to the large-$N$
analysis, the phase boundary $H_{c}(\eta)$ exhibits a reentrant behavior
within $2<D<3.065\dots$, as the $XY$-anisotropy $\eta$ changes. On the one
hand, as for the $D=(2+1)$ and $(1+1)$ short-range $XY$ magnets, the
singularities have been determined as $H_{c}(\eta)-H_{c}(0)\sim|\eta|$ and
$0$, respectively, and the transient behavior around $D\approx 2.5$ remains
unclear. As a preliminary survey, setting $(\sigma,\eta)=(1,0.5)$, we
investigate the phase transition by the agency of the fidelity, which seems to
detect the singularity at $H=H_{c}$ rather sensitively. Thereby, under the
setting $\sigma=4/3$ ($D=2.5$), we cast the fidelity data into the crossover-
scaling formula with the properly scaled $\eta$, aiming to determine the
multi-criticality around $\eta=0$. Our result indicates that the multi-
criticality is identical to that of the $D=(2+1)$ magnet, and $H_{c}(\eta)$’s
linearity might be retained down to $D>2$.

##  1 Introduction

The $XY$ chain with the transverse field is attracting much attention
Maziero10 ; Sun14 ; Karpat14  in the context of the quantum information theory
Luo18 ; Steane98 ; Bennett00 . A key ingredient is that the model covers both
the $XX$\- and $XY$-symmetric cases, and a variety of phase transitions occur,
as the transverse field and the $XY$-anisotropy change. Meanwhile, its
extention to the long-range-interaction case has been made Adelhardt20 ;
particularly, the limiting cases such as the transverse-field Ising and $XXZ$
chains with the long-range interactions have been investigated in depth
Defenu17 ; Dutta01 ; Koffel12 ; Campana10 ; Gong16 ; Fey16 ; Maghrebi17 ;
Frerot17 ; Roy18 ; Puebla19 . A notable point is that the effective
dimensionality Defenu17

$D=2/\sigma+1,$ (1)

varies, as the decay rate $\sigma$ of the long-range interactions obeying the
power law, $1/r^{1+\sigma}$ ($r$: distance between spins), changes. The
effective dimensionality $D$ governs the criticality of the transverse-field-
driven phase transition. In fact, the $D$ dependence on the criticality for
the classical counterpart Fisher72 ; Sak73 ; Gori17 ; Angelini14 ; Joyce66
has been studied extensively. A peculiarity of the quantum magnet is that the
long-range interaction induces the asymmetry between the real-space and
imaginary-time directions characterized by the dynamical critical exponent
$z\neq 1$ Defenu17 . Hence, it is significant to access to the ground state
(infinite imaginary-time system size) directly so as to get rid of the
influences caused by the aspect ratio between the real-space and imaginary-
time system sizes.

In this paper, we investigate the transverse-field $XY$ chain with the long-
range interactions Adelhardt20  by means of the exact-diagonalization method,
which enables us to access to the ground state directly. We devote ourselves
to the anisotropy-driven multi-criticality around the $XX$-symmetric point. So
far, the transverse-field-driven criticality with the fixed anisotropy has
been explored in detail Adelhardt20 . As for the large-$N$ magnet, the multi-
criticality has been studied Wald15 , and an intriguing reentrant behavior was
observed in low dimensions, $2<D<3.065\dots$. On the contrary, as for the $XY$
magnet. only the cases of $D=(2+1)$ Wald15 ; Henkel84 ; Jalal16 ; Nishiyama19
and $(1+1)$ Mukherjee11  have been studied, and the transient behavior in
between, $D\approx 2.5$, remains unclear. The aim of this paper is to shed
light on such a fractional-effective-dimensionality regime by adjusting the
algebraic decay rate $\sigma$ of the long-range interactions carefully with
the aid of the $\sigma\leftrightarrow D$ relation, Eq. (1).

To be specific, we present the Hamiltonian for the transverse-field $XY$ chain
with the long-range interactions

${\cal H}=-\frac{1}{{\cal N}}\sum_{i\neq
j}J_{ij}((1+\eta)S_{i}^{x}S_{j}^{x}+(1-\eta)S^{y}_{i}S^{y}_{j})-H\sum_{i=1}^{N}S^{z}_{i}.$
(2)

Here, the quantum spin-$1/2$ operator ${\mathbf{S}}_{i}$ is placed at each
lattice point, $i=1,2,\dots,N$. The parameters, $H$ and $\eta$, denote the
transverse field and $XY$ anisotropy, respectively. The $XY$ interaction
between the $i$-$j$ spins decays algebraically as

$J_{ij}=1/\sin(\pi|i-j|/N)^{1+\sigma},$ (3)

with the decay rate $\sigma$, and the periodic boundary condition is imposed,
i.e. ${\mathbf{S}}_{N+1}={\mathbf{S}}_{1}$, among the alignment of spins,
$\\{{\mathbf{S}}_{i}\\}$. The denominator ${\cal N}$ stands for the Kac factor
Homrighausen17 ; Vanderstraeten18

${\cal N}=\frac{2}{N}\sum_{i\neq j}1/\sin(\pi|i-j|/N)^{1+\sigma}.$ (4)

As shown in Eq. (1), the effective dimensionality $D=2/\sigma+1$ depends on
the decay rate $\sigma$, which thus governs the criticality of the transverse-
field-driven phase transition. Even quantitatively, the critical exponent for
the $D=3$ Ising model was pursued Goll18  by the $\sigma\leftrightarrow D$
correspondence. In Fig. 1, we present the criticality chart for the $XY$ case
($\eta\neq 1$) Defenu17 ; Dutta01 . For small $\sigma\leq 2/3$, the effective
dimensionality $D=4$ is realized, and the transverse-field-driven criticality
belongs to the mean-field type. On the contrary, for $\sigma\gtrsim 1.75$, the
renormalization-group analysis Defenu17  indicates that the long-range
interaction becomes irrelevant, and the criticality reduces to that of the
short-range magnet; namely, the $D=2[=(1+1)]$ universality class is realized
in this regime. The threshold $\sigma\approx 1.75$ depends on the internal
symmetry group, namely, either the $XX$\- ($\eta=0$) or $XY$-symmetric
($\eta\neq 0$) type Defenu17 . For the intermediate regime $2/3<\sigma\lesssim
1.75$, the effective dimensionality ranges within $2<D<4$. Rather technically,
around both upper and lower thresholds, there emerge notorious logarithmic
corrections Fey16 ; Luijten02 ; Brezin14 ; Defenu15 , and these regimes lie
off the present concern nonetheless.

Figure 1:  The criticality chart Defenu17 ; Dutta01  for the transverse-field
$XY$ chain with the long-range interactions (2) is presented. The $XY$ case
$\eta\neq 0$ is assumed. As the algebraic decay rate $\sigma$ varies, there
appear a number of universality classes as to the transverse-field-driven
phase transition successively. For the small decay rate $\sigma\leq 2/3$, the
criticality belongs to the mean-field type with the effective dimensionality
$D=4$. On the contrary, for $\sigma\gtrsim 1.75$, the ordinary short-range
$D=2[=(1+1)]$ universality comes out. In the intermediate regime
$2/3<\sigma\lesssim 1.75$, the criticality changes continuously, characterized
by the effective dimensionality ranging within $2<D<4$. The upper and lower
thresholds, $\sigma\approx 1.75$ and $2/3$, are affected by notorious
logarithmic corrections Fey16 ; Luijten02 ; Brezin14 ; Defenu15 , and these
regimes lie off the present concern nonetheless.

So far, the transverse-field-driven phase transition at $H=H_{c}(\eta)$ with
the fixed anisotropy $\eta$ has been investigated extensively by means of the
exact-diagonalization Fey16 ; Homrighausen17 , series-expansion Adelhardt20 ;
Fey16 , matrix-product-state Koffel12 ; Vanderstraeten18  and density-matrix-
renormalization-group Frerot17 ; Zhu18  methods. On the contrary, little
attention has been paid to the anisotropy-driven criticality, namely, the
multi-criticality at the $XX$-symmetric point $\eta=0$. As shown in Fig. 2,
for the large-$N$ magnet Wald15 , the phase boundary $H_{c}(\eta)$ exhibits a
reentrant behavior in low dimensions, $(1+1)<D<(2.065+1)$. Such a feature
indicates a counterintuitive picture that lower symmetry group indices the
disorder phase in the vicinity of the multi-critical point. On the one hand,
as for the short-range $XY$ magnet, only the cases of $D=(2+1)$ Wald15 ;
Henkel84 ; Jalal16 ; Nishiyama19  and $(1+1)$ Mukherjee11  have been
considered. The former shows the multi-criticality Riedel69 ; Pfeuty74
$H_{c}(\eta)-H_{c}(0)\sim|\eta|^{1/\phi}$ with the crossover exponent
$\phi=1$, whereas the latter exhibits no singularity $H_{c}(\eta)-H_{c}(0)=0$
at all. Then, there arises a problem how the crossover exponent behaves at a
transient point $D=2.5$. It might be anticipated that $H_{c}(\eta)$’s slope
grows monotonically with $\phi=1$ retained, as the dimensionality $D$
increases. However, in principle, the phase boundary can be curved convexly,
accompanied with a suppressed exponent $\phi<1$ in such a low-dimensionality
regime. The aim of this paper is to explore the multi-criticality for the
transient dimensionality by adjusting the decay rate $\sigma$ for the long-
range $XY$ magnet (2).

Figure 2:  A schematic drawing of the transverse-field-driven phase boundary
$H_{c}(\eta)$ with the $XY$ anisotropy $\eta$ for the transverse-field $XY$
model, is presented. The short-range $D=(2+1)$-dimensional $XY$ magnet Wald15
; Henkel84 ; Jalal16 ; Nishiyama19  exhibits a linear increase of
$H_{c}(\eta)-H_{c}(0)\sim|\eta|^{1/\phi}$ with the crossover exponent $\phi=1$
Riedel69 ; Pfeuty74 , whereas the $D=(1+1)$ model Mukherjee11  shows no
singularity, i.e., $H_{c}(\eta)={\rm const.}$, at all. It is not clear whether
the crossover exponent $\phi$ changes for the intermediate dimensionality such
as $D=2.5$; in principle, it can be curved convexly, $\phi<1$. Actually, as
for the large-$N$ magnet Wald15 , there appear a reentrant behavior in low
dimensions $(1+1)<D<(2.065+1)$.

As a probe to detect the phase transition Quan06 ; Zanardi06 ; HQZhou08 ; Yu09
; You11 , we resort to the fidelity Uhlmann76 ; Jozsa94 ; Peres84 ; Gorin06

$F(H,\Delta H)=|\langle H|H+\Delta H\rangle|,$ (5)

with the ground states, $|H\rangle$ and $|H+\Delta H\rangle$, for the
proximate interaction parameters, $H$ and $H+\Delta H$, respectively. The
fidelity (5) is readily accessible via the exact-diagonalization method, which
admits the ground-state vector $|H\rangle$ explicitly. Moreover, the fidelity
does not rely on any presumptions as to the order parameters involved Wang15 ,
and it is sensitive to generic types of phase transitions such as the $XX$\-
($\eta=0$) and $XY$-symmetric ($\eta\neq 0$) cases. In fairness, it has to be
mentioned that the information-theoretical quantifier, the so-called genuine
multipartite entanglement, detects the phase boundary $H_{c}(\eta)$ clearly
for rather restricted system sizes, $N\leq 20$ Roy18 . In this paper, we
treated the cluster with $N\leq 32$ spins, taking the advantage in that the
fidelity (5) is computationally less demanding,

The rest of this paper is organized as follows. In Sec. 2, the numerical
results are presented. In the last section, we address the summary and
discussions.

##  2 Numerical results

In this section, we present the numerical results for the transverse-field
$XY$ chain with the long-range interactions (2). We employed the exact-
diagonalization method for the cluster with $N\leq 32$ spins. Because the
exact-diagonalization method yields the ground-state vector $|H\rangle$
explicitly, one is able to calculate the fidelity, namely, the overlap between
the proximate parameters, $F=|\langle H|H+\Delta H\rangle|$ (5),
straightforwardly. Thereby, we evaluated the fidelity susceptibility Quan06 ;
Zanardi06 ; HQZhou08 ; Yu09 ; You11 ,

$\chi_{F}=-\frac{1}{N}\partial_{\Delta H}^{2}F(H,H+\Delta H)|_{\Delta H=0},$
(6)

in order to detect the signature for the criticality. The fidelity
susceptibility yields rather reliable estimates for the criticality, even
though the available system size is restricted Yu09 . According to Ref.
Albuquerque10 , the fidelity susceptibility obeys the scaling formula

$\chi_{F}=N^{\alpha_{F}/\nu}f\left((H-H_{c})N^{1/\nu}\right),$ (7)

with a certain scaling function $f$. Here, the indices $\alpha_{F}$ and $\nu$
denote the fidelity-susceptibility and correlation-length ($\xi$) critical
exponents, respectively, and these exponents describe the power-law
singularities of the respective quantities such as
$\chi_{F}\sim|H-H_{c}|^{-\alpha_{F}}$ and $\xi\sim|H-H_{c}|^{-\nu}$.

In fairness, it has to be mentioned that the fidelity susceptibility (6) was
utilized successfully for the analysis of the multi-criticality in $D=(1+1)$
dimensions Mukherjee11 . In this elaborated work, the authors took a direct
route toward the multi-critical point; here, the signal of the fidelity
susceptibility splits into sequential subpeaks, reflecting the intermittent
level crossings along $\eta=0$. Here, we took an indirect route to the multi-
critical point with the properly scaled $\eta$ for each $N$ through resorting
to the crossover-scaling theory Riedel69 ; Pfeuty74 . Before commencing
detailed crossover-scaling analyses of the multi-criticality, we demonstrate
the performance of the $\chi_{F}$-mediated simulation scheme with the fixed
anisotropy $\eta$.

###  2.1  Fidelity-susceptibility analysis of the critical point $H_{c}$ for
the fixed $(\sigma,\eta)=(1,0.5)$

As a preliminary survey, by the agency of the fidelity susceptibility (6), we
analyze the critical point $H_{c}$ with the fixed interaction parameters to
$\sigma=1$ and $\eta=0.5$, for which an elaborated series-expansion result
Adelhardt20  is available.

In Fig. 3, we present the fidelity susceptibility $\chi_{F}$ for various
values of the transverse field $H$ and the system sizes, ($+$) $N=28$,
($\times$) $30$ and ($*$) $32$, with the fixed interaction parameters,
$(\sigma,\eta)=(1,0.5)$. The fidelity susceptibility exhibits a pronounced
peak around $H\approx 0.55$, indicating that the transverse-field-driven phase
transition separating the order ($H<H_{c}$) and disorder ($H>H_{c}$) phases
takes place.

Figure 3:  The fidelity susceptibility $\chi_{F}$ (6) is plotted for various
values of the transverse field $H$ and the system sizes, ($+$) $N=28$,
($\times$) $30$, and ($*$) $32$, with the algebraic decay rate $\sigma=1$ and
the $XY$ anisotropy $\eta=0.5$. The fidelity-susceptibility peak around
$H\approx 0.55$ indicates the onset of the transverse-field-driven phase
transition.

Aiming to estimate the critical point precisely, in Fig. 4, we present the
approximate critical point $H_{c}^{*}(N)$ for $1/N^{1/\nu}$ Albuquerque10
with the fixed $(\sigma,\eta)=(1,0.5)$. Here, the approximate critical point
$H_{c}^{*}(N)$ denotes the location of the fidelity-susceptibility peak

$\partial_{H}\chi_{F}(N)|_{H=H_{c}^{*}(N)}=0,$ (8)

for each $N$, and as mentioned above, the index $\nu$ describes correlation-
length’s singularity, i.e., $\xi\sim|H-H_{c}|^{-\nu}$. The abscissa scale
$1/N^{1/\nu}$ Albuquerque10  comes from this formula; actually, the relation
$H_{c}^{*}(N)-H_{c}\sim N^{1/\nu}$ follows immediately, because the
correlation length has the same scaling dimension as that of the system size
Albuquerque10 , $\xi\sim N$. The inverse-correlation-length critical exponent
is expressed as

$1/\nu=\sigma/2+1/3,$ (9)

from the scaling formula $\nu=2\nu_{SR}(D)/\sigma$ ($\nu_{SR}(D)$:
correlation-length exponent for the short-range $D$-dimensional counterpart)
Defenu17 , the $\epsilon$-expansion result $1/\nu_{SR}(D)=2-(4-D)/3$ Amit05 ,
and Eq. (1). This approximate expression $1/\nu$ (9) is not used in the
crossover-scaling analysis in Sec. 2.3, which is the main concern of this
paper; note that the multi-critical behavior is not identical to that of the
aforementioned transverse-field-driven one.

Figure 4:  The approximate critical point $H_{c}^{*}(N)$ (8) is plotted for
$1/N^{1/\nu}$ with the fixed interaction parameters $(\sigma,\eta)=(1,0.5)$,
and the inverse correlation-length exponent $1/\nu$ (9). The least-squares fit
to these data yields an estimate $H_{c}=0.62665(5)$ in the thermodynamic limit
$N\to\infty$. A possible systematic error is considered in the text. This
result agrees with the series-expansion result (11).

The least-squares fit to the data in Fig. 4 yields an estimate
$H_{c}=0.62665(5)$ in the thermodynamic limit $N\to\infty$. Clearly, the
estimate may be affected by the extrapolation errors. In order to appreciate a
possible systematic error, we carried out the same extrapolation scheme with
an alternative $1/\nu[=\sigma/(2-0.037)/0.63]\approx 0.81$ obtained via Eq.
(32) of Ref. Defenu17  and the 3D-Ising result Deng03 . Replacing the abscissa
scale with this value $1/\nu=0.81$, we arrived at an alternative one
$H_{c}=0.62877$. The deviation from the aforementioned estimate
$H_{c}=0.62665$ appears to be $\approx 2\cdot 10^{-3}$, which seems to
dominate the aforementioned least-squares-fitting error $5\cdot 10^{-5}$;
namely, the deviation between these independent extrapolation schemes
indicates an appreciable systematic error. Hence, considering the former
$2\cdot 10^{-3}$ as the error margin, we estimate the critical point as

$H_{c}=0.6267(20).$ (10)

This result appears to agree with the series-expansion result Adelhardt20

$H_{c}\approx 0.627,$ (11)

for $\sigma=1$ and $\eta=0.5$. This value (11) was read off by the present
author from Fig. 3 of Ref Adelhardt20 ; here, the Kac factor ${\cal
N}=4\zeta(2)$ (4) has to be multiplied so as to remedy the energy-scale
difference.

A few remarks are in order. First, the agreement between the sophisticated-
series-expansion estimate $H_{c}\approx 0.627$ [Eq. (11)], and ours
$H_{c}=0.6267(20)$ [Eq. (10)] confirms the validity of the fidelity-
susceptibility-mediated simulation scheme. Second, as seen from Fig. 3, the
back-ground contributions (non-singular part) as to the fidelity-
susceptibility peak appear to be rather suppressed. Actually, as argued in the
next section, the critical exponent of the fidelity susceptibility is
substantially larger than that of the specific heat, and the fidelity
susceptibility exhibits a pronounced signature for the criticality. Such a
character was reported by the exact-diagonalization analysis for the two-
dimensional $XXZ$ magnet with the restricted $N\leq 20$ Yu09 . Last, we
explain the reason why only the large system sizes $N\approx 30$ were treated
in the extrapolation analysis as in Fig. 4. Because we are considering the
crossover critical phenomenon, the series of finite-size data exhibit two
types of scaling behaviors, such as the Ising- and $XX$-type singularities, as
$N$ increases. In our preliminary survey, the large system sizes $N\approx 30$
were found to capture the desired scaling behavior coherently, at least, for
the critical domain undertaken in this simulation study. Because of this
reason, the extrapolation scheme in Fig. 4 cannot be straightforwardly
replaced with the sophisticated extrapolation sequences.

###  2.2  Scaling analysis of the fidelity susceptibility for
$(\sigma,\eta)=(1,0.5)$

In this section, we investigate the critical behavior of the fidelity
susceptibility, following the analyses in the preceding section. For that
purpose, we rely on the scaling theory (7) for the fidelity susceptibility
developed in Ref. Albuquerque10 ; afterward, this scaling theory is extended
in order to investigate the multi-criticality Riedel69 ; Pfeuty74  around
$\eta=0$. According to the scaling argument Albuquerque10  the scaling
dimension $\alpha_{F}/\nu$ of the fidelity susceptibility satisfies the
relation

$\alpha_{F}/\nu=\alpha/\nu+z,$ (12)

with the dynamical critical exponent Defenu17

$z=\sigma/2,$ (13)

and the specific-heat critical exponent $\alpha$; namely, the specific heat
exhibits the singularity as $C\sim|H-H_{c}|^{-\alpha}$. Notably enough, the
scaling dimension $\alpha_{F}/\nu$ of the fidelity susceptibility is larger
than that of the specific-heat $\alpha/\nu$, indicating that the former should
exhibit a pronounced singularity, as compared to the latter. As a consequence,
we arrive at the expression

$\alpha_{F}/\nu=\sigma-1/3,$ (14)

from the hyper-scaling relation $\alpha=2-(1+z)\nu$ Albuquerque10 , and Eq.
(9), (12), and (13). The critical indices associated with the above scaling
formula (7) are all fixed, and now, we are able to carry out the scaling
analysis of the fidelity susceptibility without any adjustable parameters.

In Fig. 5, we present the scaling plot,
$(H-H_{c})N^{1/\nu}$-$N^{-\alpha_{F}/\nu}\chi_{F}$, for various system sizes,
($+$) $N=28$, ($\times$) $30$, and ($*$) $32$, with the fixed $\sigma=1$ and
$\eta=0.5$. Here, the scaling parameters, $H_{c}$, $1/\nu$ and
$\alpha_{F}/\nu$, are given by Eq. (10), (9), and (14), respectively. We
observe that the scaled data collapse into a scaling curve satisfactorily,
confirming the validity of the analysis in Sec. 2.1 as well as the scaling
argument Defenu17  introduced above.

Figure 5:  The scaling plot of the fidelity susceptibility,
$(H-H_{c})N^{1/\nu}$-$N^{-\alpha_{F}/\nu}\chi_{F}$, is presented for various
system sizes, ($+$) $N=28$, ($\times$) $30$, and ($*$) $32$, with
$(\sigma,\eta)=(1,0.5)$. Here, the scaling parameters, critical point, inverse
correlation-length exponent and $\chi_{F}$’s scaling dimension, are set to
$H_{c}=0.6267$ (10), $1/\nu=5/6$ (9), and $\alpha_{F}/\nu=2/3$ (14),
respectively. The fidelity-susceptibility data obey the scaling theory (7)
rather satisfactorily.

A number of remarks are in order. First, no adjustable parameter was
incorporated in the scaling analysis of Fig. 5. Actually, the scaling
parameters, $H_{c}$, $\nu$, and $\alpha_{F}/\nu$, were fixed in prior to
undertaking the scaling analysis. Last, the scaling data, Fig. 5, appear to be
less influenced by the finite-size artifact, indicating that the simulation
data already enter into the scaling regime. It is a benefit of the fidelity
susceptibility that it is not influenced by corrections to scaling very
severely Yu09 . Encouraged by this finding, we explore the multi-critical
behavior via $\chi_{F}$ in the next section.

###  2.3  Crossover-scaling analysis of the fidelity susceptibility for
$\sigma=4/3$ ($D=2.5$) around the $XX$-symmetric point

In this section, we investigate the crossover-scaling (multi-critical)
behavior around the $XX$-symmetric point $\eta=0$ via the fidelity
susceptibility. Here, we set the decay rate to $\sigma=4/3$, which corresponds
to $D=2.5$ according to Eq. (1). As mentioned in Introduction, the $D=2$
Mukherjee11  and $3$ Wald15 ; Henkel84 ; Jalal16 ; Nishiyama19  cases have
been studied, and the transient behavior in between remains unclear.

For that purpose, we incorporate a yet another parameter $\eta$ accompanied
with the crossover exponent $\phi$. Then, the aforementioned expression (7) is
extended to the crossover-scaling formula Riedel69 ; Pfeuty74

$\chi_{F}=N^{\dot{\alpha}_{F}/\dot{\nu}}g\left((H-H_{c}(\eta))N^{1/\dot{\nu}},\eta
N^{\phi/\dot{\nu}}\right),$ (15)

with a certain scaling function $g$. Here, the symbol $H_{c}(\eta)$ denotes
the critical point for each $\eta$, and the indices $\dot{\nu}$ and
$\dot{\alpha}_{F}$ denote the correlation-length and fidelity-susceptibility
critical exponents, respectively, right at the multi-critical point $\eta=0$;
namely, respective singularities are given by
$\xi\sim|H-H_{c}(0)|^{-\dot{\nu}}$ and
$\chi_{F}\sim|H-H_{c}(0)|^{-\dot{\alpha}_{F}}$. The former relation together
with $N\sim\xi$ immediately yields Riedel69 ; Pfeuty74

$H_{c}(\eta)-H_{c}(0)\sim|\eta|^{1/\phi},$ (16)

because the second argument of the crossover-scaling formula (15), $\eta
N^{\phi/\dot{\nu}}$, should be dimensionless (scale-invariant). Hence, the
crossover exponent $\phi$ governs the power-law singularity of the phase
boundary. As in Eq. (12), these critical indices satisfy the scaling relation
Albuquerque10

$\dot{\alpha}_{F}/\dot{\nu}=\dot{\alpha}/\dot{\nu}+\dot{z},$ (17)

with the specific-heat and dynamical critical exponents, $\dot{\alpha}$ and
$\dot{z}$, respectively, at the multi-critical point.

Before commencing the crossover-scaling analyses of $\chi_{F}$, we fix the set
of the multi-critical indices appearing in Eq. (15). At the $XX$-symmetric
point, the critical indices for the transverse-field-driven phase transition
were determined Zapf14  as $\dot{\alpha}=1/2$, $\dot{z}=2z$ and

$1/\dot{\nu}=\sigma.$ (18)

This index (18) is taken from Eq. (19) of Ref. Defenu17 , as it means the
mean-field value Zapf14 . Hence, from Eq. (13) and (17), we arrive at

$\dot{\alpha}_{F}/\dot{\nu}=3\sigma/2.$ (19)

The above argument completes the prerequisite for the crossover-scaling
analysis. The index $\phi$ has to be determined so as to attain a good data
collapse of the crossover-scaling plot, based on the formula (15).

In Fig. 6, we present the crossover-scaling plot,
$(H-H_{c}(\eta))N^{1/\dot{\nu}}$-$N^{-\dot{\alpha}_{F}/\dot{\nu}}\chi_{F}$,
for various system sizes, ($+$) $N=28$, ($\times$) $30$, and and ($*$) $32$,
with $\sigma=4/3$, $1/\dot{\nu}=\sigma$ (18), and
$\dot{\alpha}_{F}/\dot{\nu}=3\sigma/2$ (19). Here, the second argument of the
crossover-scaling formula (15) is fixed to $\eta
N^{\phi/\dot{\nu}}=15.2(\approx 0.15\cdot 32^{\phi/\dot{\nu}})$ under the
optimal crossover exponent $\phi=1$, and the critical point $H_{c}$ was
determined via the same scheme as that of Sec. 2.1. From Fig. 6, we see that
the crossover-scaled data fall into a scaling curve satisfactorily, indicating
that the choice $\phi=1$ should be a feasible one.

Figure 6:  The crossover-scaling plot of the fidelity susceptibility,
$(H-H_{c}(\eta))N^{1/\dot{\nu}}$-$N^{-\dot{\alpha}_{F}/\dot{\nu}}\chi_{F}$, is
presented for various system sizes, ($+$) $N=28$, ($\times$) $30$, and ($*$)
$32$, with the fixed decay rate $\sigma=4/3$ corresponding to $D=2.5$ [Eq.
(1)]. Here, the second argument of the crossover-scaling formula (15) is fixed
to $\eta N^{\phi/\dot{\nu}}=15.2$ with the optimal crossover exponent
$\phi=1$, and the other multi-critical indices are set to $1/\dot{\nu}=\sigma$
(18) and $\dot{\alpha}_{F}/\dot{\nu}=3\sigma/2$ (19). The crossover-scaled
data collapse into the scaling curve rather satisfactorily under the optimal
setting $\phi=1$.

Setting the crossover exponent to a slightly large value $\phi=1.15$, in Fig.
7, we present the crossover-scaling plot,
$(H-H_{c}(\eta))N^{1/\dot{\nu}}$-$N^{-\dot{\alpha}_{F}/\dot{\nu}}\chi_{F}$,
for various system sizes $N=28,30,32$; the symbols and the critical indices,
$1/\dot{\nu}$ and $\dot{\alpha}_{F}/\dot{\nu}$, are the same as those of Fig.
6. Here, the second argument of the crossover-scaling formula (15) is set to
$\eta N^{\phi/\dot{\nu}}=30.5$ with $\phi=1.15$. For such a large vale of
$\phi=1.15$, the hilltop data get scattered, as compared to those of Fig. 6.
Likewise, in Fig. 8, for a small value of $\phi=0.85$, we present the
crossover-scaling plot,
$(H-H_{c}(\eta))N^{1/\dot{\nu}}$-$N^{-\dot{\alpha}_{F}/\dot{\nu}}\chi_{F}$,
for various system sizes $N=28,30,32$; the symbols and the critical indices,
$1/\dot{\nu}$ and $\dot{\alpha}_{F}/\dot{\nu}$, are the same as those of Fig.
6. Here, the second argument of the crossover-scaling formula (15) is set to
$\eta N^{\phi/\dot{\nu}}=7.62$ with $\phi=0.85$. For such a small value of
$\phi=0.85$, the right-side-slope data seem to split off. Considering that the
cases, Fig. 7 and 8, set the upper and lower bounds, respectively, for $\phi$,
we conclude that the crossover exponent lies within

$\phi=1.00(15).$ (20)

Figure 7:  The crossover-scaling plot of the fidelity susceptibility,
$(H-H_{c}(\eta))N^{1/\dot{\nu}}$-$N^{-\dot{\alpha}_{F}/\dot{\nu}}\chi_{F}$, is
presented for various system sizes, ($+$) $N=28$, ($\times$) $30$, and ($*$)
$32$, with the fixed decay rate $\sigma=4/3$. Here, the second argument of the
crossover-scaling formula (15) is fixed to $\eta N^{\phi/\dot{\nu}}=30.5$ with
a slightly large crossover exponent $\phi=1.15$; the other multi-critical
indices, $\dot{\nu}$ and $\dot{\alpha}_{F}/\dot{\nu}$, are the same as those
of Fig. 6. For such a large value of $\phi=1.15$, the hilltop data get
scattered.  Figure 8:  The crossover-scaling plot of the fidelity
susceptibility,
$(H-H_{c}(\eta))N^{1/\dot{\nu}}$-$N^{-\dot{\alpha}_{F}/\dot{\nu}}\chi_{F}$, is
presented for various system sizes, ($+$) $N=28$, ($\times$) $30$, and ($*$)
$32$, with the fixed decay rate $\sigma=4/3$. Here, the second argument of the
crossover-scaling formula (15) is fixed to $\eta N^{\phi/\dot{\nu}}=7.62$ with
a slightly small crossover exponent $\phi=0.85$; the other multi-critical
indices, $\dot{\nu}$ and $\dot{\alpha}_{F}/\dot{\nu}$, are the same as those
of Fig. 6. For such a small value of $\phi=0.85$, the right-side-slope data
start to split off.

This result (20) indicates that the phase boundary for $\sigma=4/3$ ($D=2.5$)
rises up linearly, $H_{c}(\eta)\sim|\eta|$, around the multi-critical point
$\eta=0$; see Fig. 1. Namely, the multi-criticality is identical to that of
$D=3$ Wald15 ; Henkel84 ; Jalal16 ; Nishiyama19 . Hence, it is suggested that
the linearity is robust against the dimensionality $D$, and merely, the slope
grows up monotonically, as the dimensionality $D$ increases from $D=2$. In
other words, no exotic character such as the reentrant behavior predicted by
the large-$N$ theory occurs, at least, for the $XY$ magnet.

A number of remarks are in order. First, underlying physics behind the
crossover-scaling plot, Fig. 6, differs significantly from that of the
fixed-$\eta$ scaling plot, Fig. 5. Actually, the former scaling dimension,
$\dot{\alpha}_{F}/\dot{\nu}=3\sigma/2$ (19), is larger than that of the
latter, $\alpha_{F}/\nu=\sigma-1/3$ (14). Hence, the data collapse in Fig. 6
is by no means accidental, and accordingly, the crossover exponent $\phi$ has
to be adjusted rather carefully. Second, we stress that the critical indices
other than $\phi$ were fixed in prior to performing the crossover-scaling
analyses. Last, in the $\chi_{F}$-mediated analysis, no presumptions as to the
order parameters are made. Because in our study, the crossover between the
$XX$\- and $XY$-symmetric cases is concerned, it is significant that the
quantifier is sensitive to both order parameters in a systematic manner.

##  3 Summary and discussions

We investigated the transverse-field $XY$ chain with the long-range
interactions (2) by means of the exact-diagonalization method. Because the
method allows us to access directly to the ground state, one does not have to
care about the anisotropy between the real-space and imaginary-time directions
rendered by the dynamical critical exponent $z(\neq 1)$ (13). As a preliminary
survey, setting the decay rate and the $XY$ anisotropy to $\sigma=1$ and
$\eta=0.5$, respectively, we analyzed the transverse-field-driven phase
transition by the agency of the fidelity susceptibility (6). Our result
$H_{c}=0.6267(20)$ [Eq. (10)] agrees with the preceeding series-expansion
result $H_{c}\approx 0.627$ [Eq. (11)], confirming the validity of our
simulation scheme. Actually, as the scaling relation (12) indicates, the
scaling dimension of the fidelity susceptibility, $\alpha_{F}/\nu$, is larger
than that of the specific heat, $\alpha/\nu$, and the former exhibits a
pronounced peak, as shown in Fig. 3. We then turn to the analysis of the
multi-criticality around the $XX$-symmetric point $\eta=0$ under the setting
$D=2.5$ ($\sigma=4/3$). Scaling the the $XY$-anisotropy parameter $\eta$
properly, we cast the fidelity-susceptibility data into the crossover-scaling
formula (15). The crossover-scaled data fall into the scaling curve
satisfactorily under the setting, $\phi=1.00(15)$ [Eq. (20)]. The result
indicates that the phase boundary $H_{c}(\eta)$ rises up linearly,
$\sim|\eta|$, around the multi-critical point for $D=2.5$. Such a character is
identical to that of the $D=(2+1)$ case, suggesting that $H_{c}(\eta)$’s
linearity is retained in low dimensions, at least, for the $XY$ magnet.

Then, there arises a problem whether the crossover exponent is influenced by
the extention of the internal symmetry group. Actually, as for the large-$N$
magnet, the phase boundary should exhibit a reentrant behavior within
$2<D<3.065\dots$. We conjecture that for sufficiently large internal-symmetry
group, the phase boundary gets curved convexly, $\phi\approx 0.5$, for
$D\approx 2.5$, and above this threshold, the reentrant behavior eventually
sets in. The SU$(4)$ magnet Itoi00  subjected to the transverse field would be
a promising candidate to examine this scenario, and this problem is left for
the future study.

## Acknowledgment

This work was supported by a Grant-in-Aid for Scientific Research (C) from
Japan Society for the Promotion of Science (Grant No. 20K03767).

Author contribution statement

The presented idea was conceived by Y.N. He also performed the computer
simulations, analyzed the data, and wrote the manuscript.

## References

  * (1) J. Maziero, H. C. Guzman, L. C. Céleri, M. S. Sarandy, and R. M. Serra, Phys. Rev. A 82 (2010) 012106.
  * (2) Z.-Y. Sun, Y.-Y. Wu, J. Xu, H.-L. Huang, B.-F. Zhan, B. Wang, and C.-B. Duanpra, Phys. Rev. A 89 (2014) 022101.
  * (3) G. Karpat, B. Çakmak, and F. F. Fanchini, Phys. Rev. B 90 (2014) 104431.
  * (4) Q. Luo, J. Zhao, and X. Wang, Phys. Rev. E 98 (2018) 022106.
  * (5) A. Steane, Rep. Prog. Phys. 61 (1998) 117.
  * (6) C.H. Bennett and D.P. DiVincenzo, Nature 404 (2000) 247.
  * (7) P. Adelhardt, J. A. Koziol, A. Schellenberger, and K. P. Schmidt, Phys. Rev. B 102 (2020) 174424.
  * (8) N. Defenu, A. Trombettoni, and S. Ruffo Phys. Rev. B 96 (2017) 104432.
  * (9) A. Dutta and J. K. Bhattacharjee, Phys. Rev. B 64 (2001) 184106.
  * (10) T. Koffel, M. Lewenstein, and L. Tagliacozzo, Phys. Rev. Lett. 109 (2012) 267203.
  * (11) L. S. Campana, L. De Cesare, U. Esposito, M. T. Mercaldo, and I. Rabuffo, Phys. Rev. B 82 (2010) 024409.
  * (12) Z.-X. Gong, M. F. Maghrebi, A. Hu, M. Foss-Feig, P. Richerme, C. Monroe, and A. V. Gorshkov, Phys. Rev. B 93 (2016) 205115.
  * (13) S. Fey and K. P. Schmidt, Phys. Rev. B 94 (2016) 075156.
  * (14) M. F. Maghrebi, Z.-X. Gong, and A. V. Gorshkov, Phys. Rev. Lett. 119 (2017) 023001.
  * (15) I. Frérot, P. Naldest, and T. Roscilde, Phys. Rev. B 95 (2017) 245111.
  * (16) S. S. Roy and H. S. Dhar, Phys. Rev. A 99 (2019) 062318.
  * (17) R. Puebla, O. Marty, and M. B. Plenio, Phys. Rev A 100 (2019) 032115.
  * (18) M. E. Fisher, S.-k. Ma, and B. G. Nickel, Phys. Rev. Lett. 29 (1972) 917.
  * (19) J. Sak, Phys. Rev. B 8 (1973) 281.
  * (20) G. Gori, M. Michelangeli, N. Defenu, and A. Trombettoni Phys. Rev. E 96 (2017) 012108.
  * (21) M. C. Angelini, G. Parisi, and F. Ricchi-Tersenghi, Phys. Rev. E 89 (2014) 062120.
  * (22) J. S. Joyce, Phys. Rev. 146 (1966) 349.
  * (23) S. Wald and M. Henkel, J. Stat. Mech.: Theory and Experiment (2015) P07006.
  * (24) M. Henkel, J. Phys. A: Mathematical and Theoretical 17 (1984) L795.
  * (25) S. Jalal, R. Khare, and S. Lal, arXiv:1610.09845.
  * (26) Y. Nishiyama, Eur. Phys. J. B 92 (2019) 167.
  * (27) V. Mukherjee, A. Polkovnikov, and A. Dutta, Phys. Rev. B 83 (2011) 075118.
  * (28) I. Homrighausen, N. O. Abeling, V. Zauner-Stauber, J. C. Halimeh, Phys. Rev. B 96 (2017) 104436.
  * (29) L. Vanderstraeten, M. Van Damme, H. P. Büchler, F. Verstraete, Phys. Rev. Lett. 121 (2018) 090603.
  * (30) R. Goll and P. Kopietz, Phys. Rev. E 98 (2018) 022135.
  * (31) E. Luijten and H. W. J. Blöte, Phys. Rev. Lett. 89 (2002) 025703.
  * (32) E. Brezin, G. Parisi, and F. Ricci-Tersenghi, J. Stat. Phys. 157 (2014) 855.
  * (33) N. Defenu, A. Trombettoni, and A. Codello, Phys Rev. E 92 (2015) 052113.
  * (34) Z. Zhu, G. Sun, W.-L. You, and D.-N. Shi, Phys. Rev. A 98 (2018) 023607.
  * (35) E.K. Riedel and F. Wegner, Z. Phys. 225 (1969) 195.
  * (36) P. Pfeuty, D. Jasnow, and M. E. Fisher, Phys. Rev. B 10 (1974) 2088.
  * (37) H. T. Quan, Z. Song, X. F. Liu, P. Zanardi, and C. P. Sun, Phys. Rev. Lett. 96 (2006) 140604.
  * (38) P. Zanardi and N. Paunković, Phys. Rev. E 74 (2006) 031123.
  * (39) H.-Q. Zhou, and J. P. Barjaktarevic̃, J. Phys. A: Math. Theor. 41 (2008) 412001.
  * (40) W.-C. Yu, H.-M. Kwok, J. Cao, and S.-J. Gu, Phys. Rev. E 80 (2009) 021108.
  * (41) W.-L. You and Y.-L. Dong, Phys. Rev. B 84 (2011) 174426.
  * (42) A. Uhlmann, Rep. Math. Phys. 9 (1976) 273.
  * (43) R. Jozsa, J. Mod. Opt. 41 (1994) 2315.
  * (44) A. Peres, Phys. Rev. A 30 (1984) 1610.
  * (45) T. Gorin, T. Prosen, T. H. Seligman, and M. Žnidarič, Phys. Rep. 435 (2006) 33.
  * (46) L. Wang, Y.-H. Liu, J. Imriška, P. N. Ma, and M. Troyer, Phys. Rev. X 5 (2015) 031007.
  * (47) A. F. Albuquerque, F. Alet, C. Sire, and S. Capponi, Phys. Rev. B 81 (2010) 064418.
  * (48) D.J. Amit and V. Martín-Mayor, Field Theory, the renormalization group, and critical phenomena, World Scientific 2005
  * (49) Y. Deng and H.W.J. Blöte, Phys. Rev. E 68 (2003) 036125.
  * (50) V. Zapf, M. Jaime, and C. D. Batista, Rev. Mod. Phys. 86 (2014) 563.
  * (51) C. Itoi, S. Qin, and I. Affleck, Phys. Rev. B 61 (2000) 6747.


# BEATLE—Self-Reconfigurable Aerial Robot: Design, Control and Experimental
Validation

Junichiro Sugihara, Moju Zhao, Takuzumi Nishio, Kei Okada, and Masayuki Inaba
Junichiro Sugihara (corresponding author), Takuzumi Nishio, Kei Okada, and
Masayuki Inaba are with the Department of Mechano-Infomatics, The University
of Tokyo, Bunkyo-ku, Tokyo 113-8656, Japan
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>inaba@jsk.imi.i.u-tokyo.ac.jp).Moju Zhao is
with the Department of Mechanical-Engineering, The University of Tokyo,
Bunkyo-ku, Tokyo 113-8656, Japan (e-mail:chou@jsk.imi.i.u-tokyo.ac.jp).

###### Abstract

Modular self-reconfigurable robots (MSRRs) offer enhanced task flexibility by
constructing various structures suitable for each task. However, conventional
terrestrial MSRRs equipped with wheels face critical challenges, including
limitations in the size of constructible structures and system robustness due
to elevated wrench loads applied to each module. In this work, we introduce a
Aerial MSRR (A-MSRR) system named BEATLE, capable of merging and separating
in-flight. BEATLE can merge without applying wrench loads to adjacent modules,
thereby expanding the scalability and robustness of conventional terrestrial
MSRRs. In this article, we propose a system configuration for BEATLE,
including mechanical design, a control framework for multi-connected flight,
and a motion planner for reconfiguration motion. The design of a docking
mechanism and housing structure aims to balance the durability of the
constructed structure with ease of separation. Furthermore, the proposed
flight control framework achieves stable multi-connected flight based on
contact wrench control. Moreover, the proposed motion planner based on a
finite state machine (FSM) achieves precise and robust reconfiguration motion.
We also introduce the actual implementation of the prototype and validate the
robustness and scalability of the proposed system design through experiments
and simulation studies.

###### Index Terms:

Aerial robots, reconfigurable robots, distributed systems, contact wrench
control.

##  I Introduction

Figure 1:  The proposed robot platform BEATLE: BEnding-moment-free self-
reconfigurable Aerial roboT moduLE. This chronophotography visualizes
reconfiguration motion of BEATLE.

In recent years, significant progress has been made in the field of modular
self-configurable robots (MSRRs)[1, 2, 3, 4, 5], which possess the ability to
autonomously merge to new create structures. The principal objective of MSRRs
is to address various tasks utilizing a single type of robot, achieved through
the construction of structures tailored to each task. For example, mobility
can be enhanced by constructing wheel or leg structures[6, 7, 8], while
manipulation capabilities can be augmented through the creation of arm
structures[9, 10, 11]. Furthermore, the integration of multiple bodies can
yield the formation of bridge-like structures[12, 11], thereby facilitating
the movement of other robots.

However, these terrestrial MSRRs encounter the drawback of being constrained
in the scale of structures they can construct. This limitation is particularly
evident in cantilever configurations such as arm and bridge structures. In the
scenario of a cantilever arrangement, a module imposes a contact wrench on
adjacent modules due to its own weight. Therefore, considering the module on
the edge as the $1st$ module, the $i-th$ module is subject to a wrench load
induced by the weight of $(i-1)$ bodies. Consequently, as the scale of the
intended structure increases, the impact of contact wrench becomes more
pronounced. These contact wrenches can lead to overloading of actuators and
mechanical components, thereby reducing system responsiveness and structural
robustness. Hence, to address this problem, it is necessary for each module to
compensate for its own weight and control its orientation independently.

In this study, we introduce an aerial Modular Self-Reconfigurable Robot
(A-MSRR) platform named BEATLEFig. 1, designed for autonomous merging and
separation in mid-air. Similar to terrestrial MSRRs, BEATLE enhances available
force and torque by rigidly interconnecting with one another, thereby
augmenting task capabilities. Furthermore, the interlinking of aerial modules
enables the autonomous assembly of structures such as bridges and ladders at
elevated altitudes. Although terrestrial MSRRs are equipped with rotational
joints[4] at the docking points, aerial robots can execute diverse tasks
without joints[13, 14] since their bodies possess inherent three-dimensional
mobility. Additionally, BEATLE possesses the capability to exert force in an
arbitrary direction through its thrust vectoring rotor, enabling it to
compensate its own weight in any orientation without imposing torque loads on
neighboring bodies, during a multi-connected state. Hence, BEATLE facilitates
amalgamation with an arbitrary number of moudules, substantially enhancing the
scalability of MSRR.

###  I-A Related works

####  I-A1 Docking mechanism

Docking mechanism is one of the most crucial hardware components in MSRRs, and
various coupling methods have been proposed in prior research. For example,
the A-MSRR presented in [15] features a magnetic docking mechanism. Magnetic
coupling proves suitable for A-MSRRs as magnetic force can compensate for
misalignment between two docking parts, simplyfing aerial docking. However,
[15] lacks the capability for undocking motion due to the absence of a
mechanism to release magnetic connections. Although the A-MSRR in [16],
equipped with a much smaller magnetic connector, can perform aerial undocking
motion by applying a momentary bending moment to the docking part with rotor
thrusts, this smaller connector conversely complicates the docking motion. In
contrast, TRADY, an A-MSRR proposed in our previous work [17], achieves both
docking and undocking motion by employing a pole-switching mechanism with a
servo motor.

However, a remaining issue with magnetic connections is the lack of resistance
against bending moments. Particularly, to enable high-torque activities and
the assembly of large structures, a more robust mechanical coupling is
necessary. Regarding mechanical connections, certain mechanisms for
terrestrial MSRRs adopt this approach. For instance, in [18], a docking
mechanism comprised of locking “jaws” capable of hooking onto the
corresponding “jaws” on the other side is proposed. Additionally, a mechanism
consisting of movable pins that can be inserted into the holes of other
modules is also suggested in [19]. However, both these “jaws” and pins are
relatively small compared to the overall body size, which may limit their
practical durability. Therefore, this article introduces a hybrid docking
mechanism composed of both magnetic and mechanical connectors. The mechanical
connection is achieved by inserting sticks which are supported by durable
structure. Additionally, to improve the rigidity of the assembled structure
and ensure safer reconfiguration motion, a housing design for aerial modules
is proposed.

####  I-A2 Control method

Figure 2: The proposed control framework for BEATLE: (A) Fundamental Control
Flow and (B) Contact Wrench Compensation. Blocks highlighted in bold font
denote processes essential for multi-connected flight. Green arrows signify
reference inputs from upstream sources, while red arrows indicate information
exchange with other modules. Each blue block represents an individual module.
Figure 3: The overview of the proposed design for BEATLE. Top: visualization
of BEATLE module and its fundamental components. Bottom: visualization of
docking mechanisms and housing plates where (a) represents the male side and
(b) represents the female side.

The key points regarding the stability and scalability of the A-MSRR system
include robustness during multi-connected flight and stability of aerial
reconfiguration motion. These attributes are characterized by the
controllability of each module. These attributes depend on the control degree
of freedom (CDoF) of the assembled structure and each module. Concerning the
controllability of the assembled structure, A-MSRRs proposed in [15, 16, 20,
21],are underactuated systems during both single flight and multi-connected
flight. According to the general control theory for underactuated aerial robot
[22], translational control is achieved through the rotational motion of the
airframe. However, as the scale of the assembled structure increases, the
responsiveness of rotational motion decreases. Hence, when more than a certain
number of modules are connected, flight control becomes critically unstable.
On the other hand, in the case of TRADY [17], while each module is
underactuated, the merged structure is over-actuated. Consequently, TRADY can
address the instability caused by the low responsiveness of rotational motion.

On the other hand, concerning the controllability of each module, all
previously proposed A-MSRRs are underactuated. These underactuated A-MSRRs
cannot explicitly control the inter-module contact forces that arise during
multi-connected flight because underactuated modules do not have extra CDoF.
Therefore, the flight control of A-MSRR in the merged state is achieved by
employing a dynamics model that considers the entire connected structure as a
single aircraft. However, a significant issue with this approach is that the
dynamics model required for control needs updating during reconfiguration.
This updating results in discontinuous changes in the target thrust during
reconfiguration. These discrete changes lead to instability in the aircraft,
such as sudden ascents or descents [15, 17].

Another challenge arises from the accumulation of unexpected contact wrench
between adjacent modules. This is because the control process is
decentralized, with each module independently carrying out its control
procedures. Consequently, variations in target control outputs among modules
arise from deviations in self-state estimation and subtle disparities in
aircraft characteristics. The interference between modules resulting from this
discrepancy in target outputs leads to the accumulation of contact wrenches.
As the number of modules increases, these contact wrenches become more
pronounced, destabilizing flight control and reconfiguration operations, thus
significantly impacting the scalability of the system.

Therefore, BEATLE proposed in this study features a thrust configuration,
enabling each module to perform overactuated control independently. This
aircraft configuration, along with the proposed control framework, enables
explicit control of inter-module contact wrenches during multi-connected
flight. Thus, the problems of dynamics model switching and contact wrench
accumulation that were problems with previous A-MSRRs can be solved, greatly
improving system stability and scalability. This control framework is
summarized in Fig. 2.

###  I-B Contributions of this work

In summary, this paper presents an extended A-MSRR system, offering the
following novel contributions:

  1. 1.

We propose a mechanical design for A-MSRR aimed at constructing structures
with heightened rigidity.

  2. 2.

We introduce a flight control framework facilitating scalable multi-connected
flight and ensuring stable self-reconfiguration based on contact wrench
control.

  3. 3.

We validate the efficacy of the proposed A-MSRR system through a comprehensive
evaluation encompassing real-machine experiments and simulation studies.

To the best of our knowledge, this is the first work that introduces an A-MSRR
capable of conducting multi-connected flight utilizing a contact wrench
controller.

###  I-C Notation

From this section, nonbold lowercase symbols (e.g., $m$) represent scalars,
and bold symbols (e.g., $\bm{u}$) represent vectors. Superscripts (e.g,
${}^{\\{{CoG}\\}}\bm{p}$) represent the frame in which the vector or matrics
is expressed, and subscripts represent the target frame or an axis,
e.g.,${}^{\\{W\\}}\bm{r}_{\\{CoG\\}}$ represents a vector point from
${\\{W\\}}$ to ${\\{CoG\\}}$ w.r.t. ${\\{W\\}}$, whereas $u_{x}$ denotes the
$x$ component of the vector $\bm{u}$. Furthermore, superscripts (e.g.,
${}^{t}\bm{f}$, $\bm{\Lambda}^{T}$) represent transposed vectors.

###  I-D Organization

Section. II introduces the module configuration method, covering thrust
configuration and mechanical design. Subsequently, Section. III introduces the
design of basic flight controllers for the modules. Following this, Section.
IV describes the motion planning method for reconfiguration operations, while
Section. V proposes a method to convert target state and a contact wrench
control framework for multi-connected flight. Then, Section. VI evaluates the
proposed methods through simulations, and Section. VII elaborates on the
mechanical implementation of the actual system and the experimental results
obtained using it before concluding in Section. VIII.

##  II Design

###  II-A Rotor Configuration

As highlighted in Section. I, achieving overactuation for BEATLE necessitates
the use of vectoring rotors. While the minimum requirement is three rotors,
this study opts for four rotors to enhance redundancy and maintain symmetry
across the entire module structure. As depicted in Fig. 3, each rotor can
rotate its thrust direction about an axis parallel to its rotor arm within the
range of $[-\frac{\pi}{2},\frac{\pi}{2}]$. To simplify the rotor structure, we
employ the direct drive method for rotor tilting, foregoing the use of gears
or pulleys in the transmission mechanism. In summary, with four rotors and
four servos, the BEATLE module has a total control input dimension of eight.

###  II-B Docking Mechanism

Figure 4: The diagram illustrating the operation of the coupling mechanism.
Some parts are depicted in cross-sectional view for ease of viewing.
$l_{1}$-$l_{3}$ are key parameters of design.

The entire coupling system design is illustrated in Fig. 3, with the male side
depicted in a) and the female side in b). Furthermore, the operation of the
docking system is illustrated in Fig. 4. This system features two key aspects.

Initially, the female side mechanism is equipped with a switchable permanent
magnet module as proposed in [17], which is driven by a servo motor. The
tensile strength of this magnetic connection is $220\text{\,}\mathrm{N}$.
Additionally, the female side mechanism features a cone-shaped structure known
as a drogue, which compensates for alignment errors between modules. $l_{1}$
in Fig. 4 corresponds to the tolerable range of error that can be compensated
by the drogue. In this study, $l_{1}$ is set to
$70\text{\,}\mathrm{m}\mathrm{m}$ based on the position control accuracy
obtained through experiments.

Secondly, the male side incorporates a stick insertion system utilizing a
slider-crank mechanism, which is also driven by a servo motor. These sticks
are composed of CFRP pipe and PLA core, and they are supported by sockets
firmly fixed to BEATLE’s skeleton. This structure efficiently converts most of
the applied bending load on the sticks into tensile force. Thus, the torque
resistance of the proposed mechanism primarily depends on the tensile strength
of the sticks. In the case of CFRP, its tensile strength is
$2200\text{\,}\mathrm{M}\mathrm{P}\mathrm{a}$ and is durable enough for the
expected load. $l_{2}$ and $l_{3}$ are parameters representing the length of
the inserted stick and the height of the coupling mechanism, respectively,
with $l_{2}$ dependent on $l_{3}$. In this study, $l_{3}$ is set to
$170\text{\,}\mathrm{m}\mathrm{m}$ as the minimum value to ensure clearance
for the propeller, and consequently, $l_{2}$ was set to
$110\text{\,}\mathrm{m}\mathrm{m}$.

###  II-C Housing

BEATLE’s housing shown in Fig. 3 serves three main purposes: firstly, to
enhance the rigidity of the docking connection by expanding the contact area;
secondly, to prevent propeller collision during reconfiguration motion; and
lastly, to increase the utility of structures constructed with BEATLE, such as
bridges. However, as the housing surrounds the rotor, it potentially obstructs
air intake and affects thrust efficiency. To address this issue, appropriate
lightening holes are made in the housing to ensure airflow. In this study,
approximately $65\text{\,}\%$ of the housing’s footprint is cut out. Table. I
compares the thrust during hovering with and without the rotor covered by the
housing. Since the housing weighs $0.4\text{\,}\mathrm{k}\mathrm{g}$, the
thrust loss of each rotor caused by aerodynamic interference is calculated as
$0.7\text{\,}\mathrm{N}$. This value is $6.2\text{\,}\%$ of the hovering
thrust, which is not considered to be a significant problem.

TABLE I: Effect of housing on thrust. | With housing | Without housing
---|---|---
Hovering Thrust | $10.5\text{\,}\mathrm{N}$ | $12.2\text{\,}\mathrm{N}$

##  III Fundamental Flight Control

In this section, we present the design of the fundamental flight controller
for each module. The process introduced in this section corresponds to
“Fundamental Control Flow” in Fig. 2.

###  III-A Single Module Dynamics Model

Figure 5: Coordinate model for BEATLE module. TABLE II: Definition of quantities for dynamics model. Symbol | Definition
---|---
$m$ | Mass of body
I | Moment of inertia of body
$\bm{r}$ | Position of CoG
$\bm{\alpha}=[\begin{array}[]{ccc}\theta&\phi&\psi\\\ \end{array}]$ | Euler angles (roll, pitch, yaw)
$\bm{p_{i}}$ | Positiopn of i-th rotor
$\sigma_{i}\in\left[-1,\;1\right]$ | Counter torque of each rotor
$c$ | Yaw-thrust ratio
$\bm{f}$ | Resultant force
$\bm{\tau}$ | Resultant torque
$\bm{\omega}$ | Angular velociy

The coordinates model for BEATLE is depicted in Fig. 5, with each quantity
defined as shown in Table. II.

Originally, the dynamics model including rotor thrust vectoring are nonlinear
system. However, [23] proposed a method to linearize the dynamics model by
assuming that rotor thrust vectoring occurs momentarily and the intermediate
rotational motion can be disregarded. In this study, we also adopt this
approximation.

First of all, as illustrated in Fig. 5, the rotor thrust of the $i-th$ rotor
can be expressed as following:

${}^{\\{Rotor_{i}\\}}\bm{\Lambda}_{i}=\left[\begin{array}[]{ccc}\lambda_{i,x}&0&\lambda_{i,z}\end{array}\right]^{T},$
(1)

where $\bm{\Lambda}_{i}$ is the thrust force vector of $i-th$ rotor, while
$\lambda_{i,x}$ and $\lambda_{i,z}$ are rotor thrust elements aligned with the
$x$ and $z$ axis of $\\{Rotor_{i}\\}$, respectively.

Here, the coordinates conversion of thrust vector from $i-th$ rotor arm
coordinates $\\{Rotor_{i}\\}$ to center-of-gravity coordinates $\\{CoG\\}$ can
be written as follows:

$\displaystyle{}^{\\{CoG\\}}\bm{\Lambda}_{i}$
$\displaystyle={}^{\\{CoG\\}}R_{\\{Rotor_{i}\\}}{}^{\\{Rotor_{i}\\}}\bm{\Lambda}_{i}$
(2)

where ${}^{\\{CoG\\}}R_{\\{Rotor_{i}\\}}$ is the transformation matrix from
$\\{Rotor_{i}\\}$ to $\\{CoG\\}$.

Based on (12), the wrench-force ${}^{\\{{CoG}\\}}\bm{f}$ and torque
${}^{\\{{CoG}\\}}\bm{\tau}$ ) can be written as

${}^{\\{{CoG}\\}}\bm{f}=\sum\limits_{i=1}^{4}{}^{\\{CoG\\}}\bm{\Lambda}_{i},$
(3)
${}^{\\{{CoG}\\}}\bm{\tau}=\sum\limits_{i=1}^{4}\left(\left[{}^{\\{{CoG}\\}}\bm{p}_{i}\times\right]-c\>\sigma_{i}{}^{\\{{CoG}\\}}\right){}^{\\{CoG\\}}\bm{\Lambda}_{i}.$
(4)

From (3) and (4), the translational and rotational dynamics of a multirotor
unit are given by the Newton-Euler equation as followings:

$m{}^{\\{W\\}}\ddot{\bm{r}}_{\\{CoG\\}}={}^{\\{W\\}}R_{\\{CoG\\}}{}^{\\{{CoG}\\}}\bm{f}-m\bm{g},$
(5)
$\begin{split}{}^{\\{{CoG}\\}}I{}^{\\{{CoG}\\}}\dot{\bm{\omega}}&={}^{\\{{CoG}\\}}\bm{\tau}+\bm{\sigma}{}^{\\{{CoG}\\}}\bm{f}\\\
&\quad-{}^{\\{{CoG}\\}}\bm{\omega}\times{}^{\\{{CoG}\\}}I{}^{\\{{CoG}\\}}\bm{\omega},\end{split}$
(6)

where ${\\{W\\}}$ frame represent the world coordinate system and $\bm{g}$ is
gravity acceleration vector.

###  III-B Pose PID Control

Considering equations (3) and (4), a common PID controller is applied to
obtain the desired wrench on the $\\{CoG\\}$ as follows:

$\begin{split}\bm{f}^{des}=m\;{}^{\\{W\\}}&\bm{R}_{\\{CoG\\}}^{-1}\left(K_{f,p}\bm{e_{r}}+K_{f,i}\int\bm{e_{r}}dt+K_{f,d}\dot{\bm{e_{r}}}\right)\\\
&+\bm{f}^{ff},\end{split}$ (7)
$\begin{split}\bm{\tau}^{des}=\;{}^{\\{{CoG}\\}}\bm{I}\left(K_{\tau,p}\bm{e_{\alpha}}+K_{\tau,i}\int\bm{e_{\alpha}}dt+\;K_{\tau,d}\dot{\bm{e_{\alpha}}}\right)\\\
+{}^{\\{{CoG}\\}}\bm{\omega}\times\;{}^{\\{{CoG}\\}}\bm{I}\;{}^{\\{{CoG}\\}}\bm{\omega}+\bm{\tau}^{ff},\end{split}$
(8)

where $\bm{e_{r}}$ is a position error defined as
$\bm{e_{r}}={}^{\\{W\\}}\bm{r}^{des}_{\\{CoG\\}}-{}^{\\{W\\}}\bm{r}_{\\{CoG\\}}$,
and $\bm{K}_{.,p}$, $\bm{K}_{.,i}$, $\bm{K}_{.,d}$ are gains for controller.
Additionally, $f^{ff}$ and $\tau^{ff}$ represent the feedforward terms
generated in the contact wrench compensation flow, as depicted in Fig. 2 B).

###  III-C Control Allocation

Using (3) and (4), allocation from the thrust force $\bm{\lambda}$ to the
resultant wrench ${}^{\\{{CoG}\\}}\bm{W}$ can be given by following:

${}^{\\{{CoG}\\}}\bm{W}=\left(\begin{array}[]{c}{}^{\\{{CoG}\\}}\bm{f}\\\
{}^{\\{{CoG}\\}}\bm{\tau}\end{array}\right)=Q{}^{\\{CoG\\}}\bm{\Lambda}_{i},$
(9)

where

$Q=\left(\begin{array}[]{c}E\\\
\left[{}^{\\{{CoG}\\}}\bm{p}_{i}\times\right]-c\>\sigma_{i}{}^{\\{{CoG}\\}}\end{array}\right),$
(10)
$\bm{\Lambda}=\left[\begin{array}[]{cccc}{}^{\\{CoG\\}}\bm{\Lambda}_{1}^{T}&{}^{\\{CoG\\}}\bm{\Lambda}_{2}^{T}&{}^{\\{CoG\\}}\bm{\Lambda}_{3}^{T}&{}^{\\{CoG\\}}\bm{\Lambda}_{4}^{T}\\\
\end{array}\right]^{T}.$ (11)

Consequently, the desired thrust allocation can be obtained as follows:

${}^{\\{{CoG}\\}}\bm{\Lambda}^{des}=Q^{\\#}{}^{\\{{CoG}\\}}\bm{W},$ (12)

where $Q^{\\#}$ is the pseudo-inverse matrix of $Q$. Therefore, using ,
desired $i-th$ thrust vector can be obtained as follows:

${}^{\\{{Rotor_{i}}\\}}\bm{\Lambda}^{des}={}^{\\{Rotor_{i}\\}}R_{\\{CoG\\}}{}^{\\{{CoG}\\}}\bm{\Lambda}^{des},$
(13)

Lastly, the final control outputs including target thrust $\Lambda^{out}_{i}$
and tilt angle $\beta_{i}$ for each rotor can be calculate as following:

$\displaystyle\Lambda^{out}_{i}=\|\bm{\Lambda}_{i}\|,$
$\displaystyle\beta_{i}=\text{atan}(\lambda_{i,z},\lambda_{i,x}).$ (14)

##  IV Reconfiguration Motion Planning

Figure 6: The motion planning framework for aerial reconfiguration is
depicted. In A), the comprehensive structure of the Finite State Machine (FSM)
is delineated. Each block delineated with a standard line represents an
individual state, whereas those with a bold line signify a Sub-FSM. These Sub-
FSMs encompass states requisite for unitary reconfiguration motion. In B), the
motion of robots throughout the processes within Sub-FSM $i$ is visually
represented.

In this section, we propose motion planning strategies aimed at achieving
stable in-flight reconfiguration. As mentioned in our previous work [17], pose
control becomes unstable during aerial reconfiguration motion. Hence, Finite
State Machine (FSM) style motion planning approach is appropriate since it can
easily incorporate recovery actions.

As depicted in Fig. 6 A), in this study, the planning process consists of
multiple sub-FSMs to support reconfiguration motion with arbitrary number of
modules. Each sub-FSM corresponds to a specific assembly motion, with the
process transitioning to the next sub-FSM upon completion of the current one.

In a sub-FMS, the assembly motion is segmented into three states, as
visualized in Fig. 6 B). Initially, in the “Approach” state, the followr
module is guided to a position and orientation that satisfies the following
relationship with the leader structure:

${}^{\\{GCoS\\}}\bm{r}_{\\{CoG_{i}\\}}=\left[\begin{array}[]{ccc}d^{st}&0&0\end{array}\right]^{T},$
(16) ${}^{\\{GCoS\\}}R_{\\{CoG_{i}\\}}=\bm{0}$ (17)

where $d^{st}$ is the desired distance between leader and follower.

If this relationship is achieved within an acceptable margin of error, the
procedure progresses to the “Approach” state, where the follower moves towards
the leader and endeavors to establish contact. During this motion, if the pose
error exceeds an acceptable threshold, the reconfiguration process reverts to
the “Alignment” state. Through this mechanism, it becomes possible to promptly
rectify a hazardous pose situation between the two agents and resume the
reconfiguration motion without interruption.

Finally, in the “Assembly” state, the leader and follower are fixed by
activating the docking mechanism, followed by the updating of parameters such
as the quantity of interconnected modules and their arrangement.

In the event of a disassembly motion, this can be achieved by transitioning
these states in reverse order.

##  V Multi-Connected Flight

In this section, we introduce the methods to achieve stable multi-connected
flight and aerial reconfiguration. Firstly, we present the method of target
state conversion, followed by the framework for contact wrench compensation.

###  V-A Target State Conversion

During multi-connected flight, the target position and attitude of the
geometric center of the entire structure (GCoS) is transmitted from upstream.
Hence, for the initial step, the target state of the GCoS should be
transformed into the CoG of each module. This process corresponds to “Target
State Conversion” block depicted in Fig. 2.

To begin, the position of the GCoS is defined as the average of the positions
of CoGs of connected modules, as follows:

$^{\\{W\\}}\bm{r}_{\\{GCoS\\}}=\frac{1}{n}\sum^{n}_{i=1}{}^{\\{W\\}}\bm{r}_{\\{CoG\\},i},$
(18)

where $n$ represents the number of connected modules. Fig. 7 visualizes the
relationship between $\\{CoG_{i}\\}$ in the case of four modules.
Subsequently, the conversion is derived using a homogeneous transformation
matrix as follows:

$\displaystyle\left[\begin{array}[]{c}^{\\{W\\}}\bm{r}^{des}_{\\{CoG\\},i}\\\
1\\\
\end{array}\right]=H_{i}\left[\begin{array}[]{c}^{\\{W\\}}\bm{r}^{des}_{\\{GCoS\\}}\\\
1\\\ \end{array}\right],$ (23) $\displaystyle
H_{i}=\left(\begin{array}[]{cc}{}^{\\{W\\}}\bm{R}^{des}_{\\{GCoS\\}}&\bm{d}^{3\times
1}_{i}\\\ \bm{0}^{1\times 3}&1\\\ \end{array}\right),$ (26)
$\displaystyle\bm{d}_{i}={}^{\\{W\\}}\bm{r}_{\\{CoG\\},i}-^{\\{W\\}}\bm{r}_{\\{GCoS\\}},$

where, ${}^{\\{W\\}}\bm{R}^{des}_{\\{GCoS\\}}$ represents the rotation matrix
equivalent to the target Euler angle of the entire structure.

###  V-B Contact Wrench Compensation

Figure 7: Dynamics and coordinates model during four-connected flight.

The conversion process outlined in (23) theoretically generates target
positions and orientations for each module that do not interfere with each
other. However, in practice, unexpected contact wrenches accumulate between
adjacent modules due to deviations in target poses among modules. This
deviation arises from slight offsets in the sensors and subtle differences in
aircraft characteristics between modules. The contact wrenches not only affect
the stability of normal flight but also significantly destabilize the assembly
and disassembly motion by abruptly appearing or disappearing during the
process. Hence, we require a framework to dynamically and explicitly control
this contact wrench as illustrated in B) block in Fig. 2. Following this, we
elaborate on the design of the contact wrench filter and outline the method
for calculating the feed-forward term.

####  V-B1 Dynamics in Multi-connected State

To simplify the dynamics model, we assume that the connection between modules
is sufficiently strong and that no relative motion is generated between them.
Therefore, according to the action-reaction law, the external wrench
$\bm{W}^{ext}_{i}$ applied to the $i-th$ module can be expressed as follows:

$\bm{W}^{ext}_{i}=\bm{W}_{i}^{cont}-\bm{W}_{i+1}^{cont}+\bm{W}^{dist},$ (27)

where
$\bm{W}_{i}^{cont}=\left[{}^{t}\bm{f}^{cont}_{i},\;{}^{t}\bm{\tau}^{cont}_{i}\right]^{T}$
is the contact wrench applied to the $i-th$ module by the
$\left(i-1\right)-th$ module as illustrated in Fig. 7, and $\bm{W}^{dist}$
represents the disturbance wrench like aerodynamic effects. Note that in this
study, we assume the disturbance wrench affects all modules equally.

####  V-B2 External Wrench Estimation

Generally, two methods are used to estimate external wrenches like
$\bm{W}^{ext}$ in (27): momentum-based, utilizing velocity data, and
acceleration-based. Due to substantial acceleration noise in the IMU sensors
on Beatle, we employ the full momentum-based observer, as also utilized in
[24]. By rewriting (3) and (4) to incorporate external wrenches, we derive the
following relationship:

$\displaystyle\dot{\bm{p}}$
$\displaystyle=M\dot{\bm{v}}=J^{t}\bm{W}+\bm{W}^{ext}-N,$ (28) $\displaystyle
M$ $\displaystyle=\left(\begin{array}[]{cc}mE_{3\times 3}&0_{3\times 3}\\\
0_{3\times 3}&I\\\
\end{array}\right),J=\left(\begin{array}[]{cc}R^{T}&0_{3\times 3}\\\
0_{3\times 3}&E_{3\times 3}\\\ \end{array}\right),$ (33) $\displaystyle N$
$\displaystyle=\left(\begin{array}[]{c}m\bm{g}\\\
\bm{\omega}\times\;\bm{I}\;\bm{\omega}\\\ \end{array}\right),$ (36)

where $\bm{p}$ and $\bm{v}$ represent the generalized momentum and velocity,
respectively, while $\bm{W}$ and $\bm{W}^{ext}$ denote the wrench generated
from rotor thrust and the external wrench, respectively. Subscripts and
superscripts consistent with (3) and (4) are omitted. According to [25], the
estimated external wrench can be expressed as the following residual term:

$\hat{\bm{W}}^{ext}=K_{I}\left[\bm{p}\left(t\right)-\bm{p}\left(t_{0}\right)-\int_{t_{0}}^{t}\left(J^{T}\bm{W}+\hat{\bm{W}}^{ext}-N\right)dt\right],$
(37)

where, $t$ and $t_{0}$ represent the current and initial time instances
respectively, and $K_{I}\in R^{6\times 6}$ is a positive-definite gain matrix.
Taking the time derivative of (37), we obtain the following dynamics for the
residual vector:

$\dot{\hat{\bm{W}}}^{ext}+K_{I}\hat{\bm{W}}^{ext}=K_{I}\bm{W}^{ext}.$ (38)

(38) reveals a first-order low-pass filter that converges to $\bm{W}^{ext}$ as
$t\to\infty$ approaches infinity, with the convergence speed dependent on the
matrix $K_{I}$.

Note that this method only requires the estimation of the generalized velocity
$\bm{v}$, which is less susceptible to sensor offsets than pose estimation.

####  V-B3 Contact Wrench Filtering

The quantity required for $i-th$ module to compensate for the contact wrench
is $\bm{W}^{cont}_{i}$. However, $\hat{\bm{W}}^{ext}_{i}$ obtained in (37)
includes not only $\bm{W}^{cont}_{i}$ but also $\bm{W}^{cont}_{i+1}$ and
$\bm{W}^{dist}$ as implied by (28). Therefore, a filter to remove these values
from $\hat{\bm{W}}^{ext}_{i}$ is necessary. To decompose
$\hat{\bm{W}}^{ext}_{i}$ into each element, the proposed filter utilizes the
values of $\hat{\bm{W}}^{ext}_{k}$ from all other connected modules. According
to (28), $\bm{W}^{dist}$ can be derived by summing up these estimated external
wrench values as follows:

$\bm{W}^{dist}=\frac{1}{n}\sum^{n}_{k=1}\hat{\bm{W}}^{ext}_{k}.$ (39)

Note that each contact wrench cancels out, leaving only the $\bm{W}^{dist}$.
Now, as $1st$ module is positioned at the left end, we have
$\bm{W}^{cont}_{1}=0$. Therefore, from (28), $\bm{W}^{cont}_{2}$ is obtained
as follows.

$\bm{W}^{cont}_{2}=\hat{\bm{W}}^{ext}_{1}-\bm{W}^{dist}$ (40)

In a same manner,
$\bm{W}^{cont}_{3},\;\bm{W}^{cont}_{4},\;\dots\bm{W}^{cont}_{n}$ are derived
sequentially.

Note that this process necessitate communication among modules, hence the
control frequency of the feedforward term relies on the bandwidth of this
communication. However, it is crucial to emphasize that this process operates
independently of normal flight control, thus any time delay does not adversely
affect the frequency of flight control execution. Moreover, since we assume
that the contact wrench is approximately steady in the case of BEATLE, any
incurred time delay is not considered significant in ensuring stability.

####  V-B4 Feedforward Term for Wrench Compensation

To control the contact wrench, an appropriate feedforward term
$\bm{W}^{ff}=\left[{}^{t}\bm{f}^{ff}\;\;{}^{t}\bm{\tau}^{ff}\right]^{T}$ is
incorporated into the fundamental pose control output of each module as shown
in (7) and (8).

One important consideration is that this feedforward term has a cascading
effect on other aircraft as well. First of all, the relationship between the
$\bm{W}^{cont}_{i}$ at discrete times $t$ and $t+1$ is expressed by the
following recurrence relation:

$\bm{W}^{cont}_{i}(t+1)=\bm{W}^{cont}_{i}(t)-\bm{W}^{ff}_{i}(t+1)+\bm{W}^{ff}_{i-1}(t+1),$
(41)

where $\bm{W}^{ff}_{i}(t)$ is a feedforward wrench add to $i-th$ module at
time $t$. Note that (41) includes $\bm{W}^{ff}_{i-1}$ which is the feedforward
term added to $(i-1)-th$ module. Here, (41) is subject to the following
constraints:

$\displaystyle\bm{W}^{cont}_{i}(t+1)$
$\displaystyle=^{des}\bm{W}^{cont}_{i}(t+1),$ (42)
${}^{des}\bm{W}^{cont}_{1}(t+1)$
$\displaystyle=\bm{W}^{cont}_{1}(t)=\bm{W}^{ff}_{0}(t+1)=\bm{0},$ (43)

where ${}^{des}\bm{W}^{cont}_{i}(t+1)$ is the desired contact wrench occurring
between module $(i-1)$ and module $i$, which is set to $\bm{0}$ in this study.

Under the conditions of (42) and (43), solving the recurrence equation (41)
yields the target feedforward term at time $t+1$ as follows:

$\bm{W}^{ff}_{i}(t+1)=\sum^{i}_{k=1}\\{\bm{W}^{cont}_{k}(t)-^{des}\bm{W}^{cont}_{k}(t+1)\\}.$
(44)

##  VI Simulation Studies

In this study, we conduct simulation experiments that are challenging to carry
out in the real world due to safety and resource constraints.

###  VI-A Simulation Setup

Regarding the simulation conditions, we utilize Gazebo, ensuring that each
physical parameters matches those of the real machine. Furthermore, we
introduce thfe Gaussian noise to the pose information obtained from Gazebo to
replicate sensor noise in the real machine. The mean and standard deviation of
the noise are $0.05$ and $0.01$, respectively.

###  VI-B Validation of Contact Wrench Compensation

Figure 8: The Results of multi-connected hovering and in-flight separation
simulation without the proposed wrench compensation method: (a) and (b) show
the positional and rotational errors of the left module. Note that only the
errors in the left module are shown for simplicity; (c) and (d) represent the
contact force and torque occurring between two modules. (i) through (iii) in
the graph correspond to (i) through (iii) in the snapshots.

Firstly, we validate our hypothesis regarding contact wrench accumulation and
the effectiveness of the proposed method in eliminating this contact wrench.
Note that, similar to real world, the self-position estimation of each module
is slightly deviated due to sensor noise in the simulation. Fig. 8 displays
the outcomes of multi-connected flight and disassembly motion without the
proposed contact wrench compensation. According to Fig. 8 (c) and (d), it is
evident that the contact wrench gradually accumulates at the coupling between
two modules. As Fig. 8 (a) and (b) reveal, when the coupling is released, the
accumulated wrench is released, causing the divergence of pose control of each
module.

###  VI-C Scaling Feasibility Verification

Figure 9: The results of takeoff and hovering simulation with different numbers of connected modules: (a), (b), and (c) show the results for one, four, and eight moudules, respectively. Regarding the positional and rotational errors, the values for the module closest to the center of the structure are shown as representative values. TABLE III: RMSE during hovering in each numer of connected modules n | RMSE
---|---
$e_{x}$[m] |  $e_{y}$[m] |  $e_{z}$[m] |  $e_{\theta}$[rad] |  $e_{\phi}$[rad] |  $e_{\psi}$[rad]
1 | 2.0e-3 | 1.9e-3 | 3.8e-3 | 1.2e-3 | 1.5e-3 | 1.3e-3
4 | 4.0e-3 | 3.1e-3 | 4.0e-3 | 1.4e-3 | 3.1e-3 | 3.5e-3
8 | 1.4e-3 | 3.9e-3 | 3.0e-3 | 1.0e-3 | 1.1e-3 | 1.5e-3

Next, we verify the scalability of our system by comparing the stability of
flight with different numbers of connected modules. In this simulation,
takeoff and hovering motions with proposed contact wrench compensation method
are observed. The results are shown in Fig. 9 and Table. III. From the results
in Table. III, it is evident that increasing the number of connected aircraft
does not affect the accuracy of pose control. This implies that the proposed
system is scalable to an arbitrary number of modules.

##  VII Experiment

###  VII-A Implementation

Figure 10: The snapshot of the prototype of BEATLE module. (a): Exterior. (b):
Detailed components. (c): Male-side docking mechanism. (d): Female-side
docking mechanism.

####  VII-A1 Hardware Configuration

Based on the design proposed in Section. II, we implemente the prototype of
BEATLE as depicted in Fig. 10.

The skeleton and housing primarily consist of CFRP, while other components are
constructed from PLA. The CFRP plate has a thickness of
$2\text{\,}\mathrm{m}\mathrm{m}$, and the CFRP cylinder has a thickness of
$1\text{\,}\mathrm{m}\mathrm{m}$.

The rotor module consists of a 9-inch propeller (GEMFAN 9045-3), a brushless
motor (T-Motor AT2814), and a servo motor (KRS3304R2). Each rotor is driven by
an ESC (T-Motor AIR 40A). With this configuration, each rotor generates thrust
ranging from approximately $1\text{\,}\mathrm{N}$ to $20\text{\,}\mathrm{N}$
at a voltage of $25.2\text{\,}\mathrm{V}$ ($6\text{\,}\mathrm{S}$).

Regarding processors, our original flight control unit, named “Spinal,” is
installed, equipped with an onboard IMU and an STM32 MUC. The “Fundamental
Control Flow” in Fig. 2 runs at a frequency of
$500\text{\,}\mathrm{H}\mathrm{z}$ on the Spinal. Additionally, as a high-
level processor, an onboard PC (VIM 4) is equipped, where “Trajectory
Conversion” and “Contact Wrench Compensation” in Fig. 2 operate at
$40\text{\,}\mathrm{H}\mathrm{z}$.

The primary characteristics of the developed BEATLE prototype are presented in
Table. IV.

####  VII-A2 External Components

In the real-machine experiment, we utilize a motion capture system for
position estimation and a laptop to transmit upstream commands. Communication
between these external components and the robot occurs via a Wi-Fi network.

TABLE IV: main parameters of robot Parameter | Value
---|---
Mass of main body | $4.2\text{\,}\mathrm{k}\mathrm{g}$
Size of main body |  $0.52\text{\,}\mathrm{m}$ $\times$$0.52\text{\,}\mathrm{m}$
Battery | $6000\text{\,}\mathrm{m}\mathrm{A}\mathrm{h}$
Continuous flight time | $6\text{\,}\mathrm{m}\mathrm{i}\mathrm{n}$

###  VII-B Trajectory Tracking

TABLE V: RMSE during trajectory tracking in each numer of connected modules n | RMSE
---|---
$e_{\theta}$[rad] |  $e_{\phi}$[rad] |  $e_{\psi}$[rad]
1 | 1.5e-2 | 1.9e-2 | 4.3e-2
2 | 1.7e-2 | 2.8e-2 | 4.8e-2
3 | 2.1e-2 | 4.8e-2 | 4.3e-2
Figure 11: The results of trajectory tracking experiments with different
numbers of connected modules. (a), (b), and (c) show the results for one, two,
and three modules, respectively. Left: visualized trajectory tracking; right:
in-flight attitude control error.

We begin by introducing the trajectory tracking experiments. In this
experiment, various numbers of connected modules track a lemniscate trajectory
defined as $r^{2}=\sqrt{\cos\left(2\theta\right)}$. Here, $\theta$ transitions
$0\to 2\pi$ over $100\text{\,}\mathrm{s}$ at a constant rate. Consequently,
the target translational acceleration changes nonlinearly.

In this experiment, we utilize 1 to 3 modules. The results for each number of
connected modules are depicted in Fig. 11 and Table. V.

From these results, it can be observed that even with an increase in the
number of connected modules, the modules can achieve various accelerations
while maintaining flight stability. Although the rotational error slightly
increases as the number of modules increases, this is primarily due to the
increase in aerodynamic interference caused by the surrounding walls and is
considered specific to this experimental environment.

###  VII-C In-flight assembly and disassembly

TABLE VI: RMSE during reconfiguration motion $e_{x}$[m] | $e_{y}[m]$ | $e_{z}[m]$ | $e_{\theta}[rad]$ | $e_{\phi}[rad]$ | $e_{\psi}[rad]$
---|---|---|---|---|---
2.3e-2 | 2.2e-2 | 2.5e-2 | 1.6e-2 | 5.2e-2 | 4.8e-2

Next, we demonstrate an in-flight reconfiguration experiment. In this
experiment, three modules merge sequentially, and subsequently, they separate
sequentially.

The result of the experiment is depicted in Fig. 12, Fig. 13, and Table. VI.
According to Fig. 12 (c)-(f), it can be observed that the contact wrench
between modules instantaneously occurs when a module docks, but this wrench is
instantly eliminated, allowing the merged modules to continue stable flight.
This result indicates that the proposed wrench compensation method functioned
properly. Fig. 12 (a) and (b) also demonstrate that pose instability, similar
to those observed in previous conventional A-MSRRs [15, 16, 17], does not
occur during reconfiguration motions.

On the other hand, slight instabilities in position and orientation can be
observed just before assembly (around $20\text{\,}\mathrm{s}$) and just after
separation (around $70\text{\,}\mathrm{s}$). These instabilities are
presumably caused by aerodynamic interference between modules that are very
close together.

In this experiment, hovering was sustained for more than
$30\text{\,}\mathrm{s}$ after the three bodies merged. However, no
accumulation of contact wrenches was observed and the position error and
attitude error were both within $0.025\text{\,}\mathrm{m}$ and
$0.048\text{\,}\mathrm{r}\mathrm{a}\mathrm{d}$, respectively.

Figure 12: The Results of the in-flight reconfiguration experiment: (a) and
(b) show the positional and rotational errors. Note that the values for the
module closest to the center of the structure are shown as representative
values; (c) and (d) represent the contact force and torque occurring between
module 1 and module 2, while (e) and (f) represent those between module 2 and
module 3. (ii) through (vii) in the graph correspond to (ii) through (vii) in
the snapshots.  Figure 13: The snapshots of the in-flight reconfiguration
experiment. From left to right, the modules are designated as module1,
module2, and module3. (i)-(iv) represent assembly motion, while (v)-(viii)
depict disassembly motion.

##  VIII Conclusion

In this paper, we present the design and control framework for a Aerial
Modular Self-Reconfigurable Robot (A-MSRR) named BEATLE, capable of connecting
and separating in-flight. BEATLE is equipped with a coupling mechanism
consisting of both magnetic and mechanical locking mechanisms, as well as a
robust housing, enabling the construction of a highly rigid structure. By
proposing a thrust configuration with vectoring rotors that enables over-
actuation and a control framework for multi-connected flight incorporating
contact wrench controller, we have significantly enhanced the scalability and
stability of conventional methods. The feasibility and scalability of the
proposed framework have been validated through simulation studies and
experiments with implemented prototype robots, which encompass hovering,
trajectory tracking, and in-flight self-reconfiguration.

Although we solely focus on A-MSRR with rigid and one-direction connections in
this work, it is expected that the proposed method can be extended to A-MSRR
systems with joint and multi-direction connections in future work. With the
extended A-MSRR system, various structures such as bridges with complicated
truss structures, multi-link arms connected to ceilings, and multi-legged
robots capable of terrestrial movement can be realized.

## References

  * [1] M. Jorgensen, E. Ostergaard, and H. Lund, “Modular atron: modules for a self-reconfigurable robot,” vol. 2. IEEE, 2004, pp. 2068–2073. [Online]. Available: http://ieeexplore.ieee.org/document/1389702/
  * [2] M. Yim, W. min Shen, B. Salemi, D. Rus, M. Moll, H. Lipson, E. Klavins, and G. Chirikjian, “Modular self-reconfigurable robot systems [grand challenges of robotics],” _IEEE Robotics and Automation Magazine_ , vol. 14, pp. 43–52, 3 2007. [Online]. Available: http://ieeexplore.ieee.org/document/4141032/
  * [3] V. Zykov, E. Mytilinaios, M. Desnoyer, and H. Lipson, “Evolved and designed self-reproducing modular robotics,” _IEEE Transactions on Robotics_ , vol. 23, pp. 308–319, 4 2007. [Online]. Available: http://ieeexplore.ieee.org/document/4154840/
  * [4] S. S. R. Chennareddy, A. Agrawal, and A. Karuppiah, “Modular self-reconfigurable robotic systems: A survey on hardware architectures,” _Journal of Robotics_ , vol. 2017, pp. 1–19, 2017. [Online]. Available: https://www.hindawi.com/journals/jr/2017/5013532/
  * [5] A. Brunete, A. Ranganath, S. Segovia, J. P. de Frutos, M. Hernando, and E. Gambao, “Current trends in reconfigurable modular robots design,” _International Journal of Advanced Robotic Systems_ , vol. 14, p. 172988141771045, 5 2017. [Online]. Available: http://journals.sagepub.com/doi/10.1177/1729881417710457
  * [6] H. Wei, Y. Chen, J. Tan, and T. Wang, “Sambot: A self-assembly modular robot system,” _IEEE/ASME Transactions on Mechatronics_ , vol. 16, pp. 745–757, 8 2011. [Online]. Available: http://ieeexplore.ieee.org/document/5634125/
  * [7] J. Davey, N. Kwok, and M. Yim, “Emulating self-reconfigurable robots - design of the smores system.” IEEE, 10 2012, pp. 4464–4469. [Online]. Available: http://ieeexplore.ieee.org/document/6385845/
  * [8] C. Liu, Q. Lin, H. Kim, and M. Yim, “Smores-ep, a modular robot with parallel self-assembly,” _Autonomous Robots_ , vol. 47, pp. 211–228, 2 2023. [Online]. Available: https://link.springer.com/10.1007/s10514-022-10078-1
  * [9] G. Liang, H. Luo, M. Li, H. Qian, and T. L. Lam, “Freebot: A freeform modular self-reconfigurable robot with arbitrary connection point - design and implementation.” IEEE, 10 2020, pp. 6506–6513.
  * [10] P. Swissler and M. Rubenstein, “Fireant3d: a 3d self-climbing robot towards non-latticed robotic self-assembly.” IEEE, 10 2020, pp. 3340–3347. [Online]. Available: https://ieeexplore.ieee.org/document/9341116/
  * [11] Y. Tu, G. Liang, and T. L. Lam, “Freesn: A freeform strut-node structured modular self-reconfigurable robot - design and implementation.” IEEE, 5 2022, pp. 4239–4245. [Online]. Available: https://ieeexplore.ieee.org/document/9811583/
  * [12] T. Tosun, J. Daudelin, G. Jing, H. Kress-Gazit, M. Campbell, and M. Yim, “Perception-informed autonomous environment augmentation with modular robots.” IEEE, 5 2018, pp. 6818–6824. [Online]. Available: https://ieeexplore.ieee.org/document/8463155/
  * [13] R. Watson, M. Kamel, D. Zhang, G. Dobie, C. MacLeod, S. G. Pierce, and J. Nieto, “Dry coupled ultrasonic non-destructive evaluation using an over-actuated unmanned aerial vehicle,” _IEEE Transactions on Automation Science and Engineering_ , vol. 19, pp. 2874–2889, 10 2022.
  * [14] R. R. Martinez, H. Paul, and K. Shimonomura, “Aerial torsional work utilizing a multirotor uav with add-on thrust vectoring device,” _Drones_ , vol. 7, p. 551, 8 2023.
  * [15] D. Saldana, B. Gabrich, G. Li, M. Yim, and V. Kumar, “Modquad: The flying modular structure that self-assembles in midair.” Institute of Electrical and Electronics Engineers Inc., 9 2018, pp. 691–698.
  * [16] D. Saldana, P. M. Gupta, and V. Kumar, “Design and control of aerial modules for inflight self-disassembly,” _IEEE Robotics and Automation Letters_ , vol. 4, pp. 3402–3409, 10 2019.
  * [17] J. Sugihara, T. Nishio, K. Nagato, M. Nakao, and M. Zhao, “Design, control, and motion strategy of trady: Tilted-rotor-equipped aerial robot with autonomous in-flight assembly and disassembly ability,” _Advanced Intelligent Systems_ , vol. 5, no. 10, p. 2300191, 2023. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/aisy.202300191
  * [18] B. Salemi, M. Moll, and W. min Shen, “Superbot: A deployable, multi-functional, and modular self-reconfigurable robotic system.” IEEE, 10 2006, pp. 3636–3641. [Online]. Available: http://ieeexplore.ieee.org/document/4058969/
  * [19] K. C. Wolfe, M. S. Moses, M. D. Kutzer, and G. S. Chirikjian, “M¡sup¿3¡/sup¿express: A low-cost independently-mobile reconfigurable modular robot.” IEEE, 5 2012, pp. 2704–2710. [Online]. Available: http://ieeexplore.ieee.org/document/6224971/
  * [20] R. Oung, A. Ramezani, and R. D’Andrea, “Feasibility of a distributed flight array,” in _Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference_ , 2009, pp. 3038–3044.
  * [21] R. Oung, F. Bourgault, M. Donovan, and R. D’Andrea, “The distributed flight array,” in _2010 IEEE International Conference on Robotics and Automation_ , 2010, pp. 601–607.
  * [22] “Opportunities and challenges with autonomous micro aerial vehicles,” vol. 31, 9 2012, pp. 1279–1291.
  * [23] “The voliro omniorientational hexacopter: An agile and maneuverable tiltable-rotor aerial vehicle,” _IEEE Robotics and Automation Magazine_ , vol. 25, pp. 34–44, 12 2018.
  * [24] M. Ryll, G. Muscio, F. Pierri, E. Cataldi, G. Antonelli, F. Caccavale, and A. Franchi, “6d physical interaction with a fully actuated aerial robot.” IEEE, 5 2017, pp. 5190–5195. [Online]. Available: http://ieeexplore.ieee.org/document/7989608/
  * [25] A. de Luca and R. Mattone, “Sensorless robot collision detection and hybrid force/motion control.” IEEE, 2005, pp. 999–1004. [Online]. Available: http://ieeexplore.ieee.org/document/1570247/

|  Junichiro Sugihara is a PhD student in the Department of Mechano-
Informatics, School of Information Science and Technology, University of
Tokyo. He received BE degree from the Department of Mechanical Engineering,
University of Tokyo, in 2023. His research interests are mechanical design,
modeling and control of aerial robots, and self-reconfigurable modular robot.
---|---
|  Moju Zhao is currently an Associate Professor at The University of Tokyo.
He received Doctor Degree from the Department of Mechano-Informatics, The
University of Tokyo, 2018. His research interests are mechanical design,
modelling and control, motion planning, and vision based recognition in aerial
robotics. His main achievement is the articulated aerial robots which have
received several awards in conference and journal, including the Best Paper
Award in IEEE ICRA 2018.
---|---
|  Takuzumi Nishio is an Assistant Professor in the Department of Mechano-
Informatics, School of Information Science and Technology, University of
Tokyo. He received PhD degree from the Department of Interdisciplinary-
Informatics, University of Tokyo, in 2022. His research interests are
mechanical design, modeling and control of aerial robots, and vision-based
recognition and motion planning of manipulators in field robotics.
---|---
|  Kei Okada received BE in Computer Science from Kyoto University in 1997. He
received MS and PhD in Information Engineering from The University of Tokyo in
1999 and 2002, respectively. From 2002 to 2006, he joined the Professional
Programme for Strategic Software Project in The University Tokyo. He was
appointed as a lecturer in the Creative Informatics at the University of Tokyo
in 2006 and an associate professor in the Department of Mechano-Informatics in
2009. His research interests include humanoids robots, real-time 3D computer
vision, and recognitionaction integrated system.
---|---
|  Masayuki Inaba is a professor in the Department of Mechano-Informatics,
School of Information Science and Technology, University of Tokyo. He received
PhD degree from the Graduate School of Information Engineering, University of
Tokyo, in 1986. He is currently a professor in the Department of Mechano-
Informatics, University of Tokyo, in 2000. His research interests include
robotic systems and software architectures to advance robotics research.
---|---