# A Kirchhoff-like theory for hard magnetic rods under geometrically nonlinear
deformation in three dimensions

Tomohiko G. Sano111Current affiliation: Laboratory for Structural
Instabilities in Geometry and Materials, Department of Mechanical Engineering,
Keio University, 3-14-1 Hiyoshi, Kohoku-ku, Yokohama, Kanagawa, 223-8522,
Japan Matteo Pezzulla222Current affiliation: Slender Structures Lab,
Department of Mechanical and Production Engineering, Århus University, Inge
Lehmanns Gade 10, 8000 Århus C, Denmark Pedro M. Reis<EMAIL_ADDRESS>Flexible Structures Laboratory, Institute of Mechanical Engineering, École
polytechnique fédérale de Lausanne, Switzerland.

###### Abstract

Magneto-rheological elastomers (MREs) are functional materials that can be
actuated by applying an external magnetic field. MREs comprise a composite of
hard magnetic particles dispersed into a nonmagnetic elastomeric (soft)
matrix. By applying a strong magnetic field, one can magnetize the structure
to program its deformation under the subsequent application of an external
field. There is a variety of types of MREs depending on the value of their
coercivity (_i.e._ the necessary field strength to erase the magnetization)
that can be broadly classified into soft or hard MREs. Hard MREs, whose
coercivities are large, have been receiving particular attention because the
programmed magnetization remains unchanged upon actuation. Hence, once a
structure made of a hard MRE is magnetized, it can be regarded as magnetized
permanently. Motivated by a new realm of applications, there have been
significant theoretical developments in the continuum (3D) description of hard
MREs. By reducing the 3D description into 1D or 2D via dimensional reduction,
several theories of hard magnetic slender structures such as linear beams,
elastica, and shells have been recently proposed. In this paper, we derive an
effective theory for MRE rods (slender, mono-dimensional structures) under
geometrically nonlinear 3D deformation. Our theory is based on reducing the 3D
magneto-elastic energy functional for the hard MREs into a 1D Kirchhoff-like
description (centerline-based). Restricting the theory to 2D, we reproduce
previous works on planar deformations. For further validation in the general
case of 3D deformation, we perform precision experiments with both naturally
straight and curved rods under either constant or constant-gradient magnetic
fields. Our theoretical predictions are in excellent agreement with both
discrete simulations and precision-model experiments. Finally, we discuss some
limitations of our framework, as highlighted by the experiments, where long-
range dipole-dippole interactions, which are neglected in the theory, can play
a role.

###### keywords:

Magnetorheological elastomer , Kirchhoff rod theory , Slender structures ,
Buckling

††journal:  Journal of the Mechanics and Physics of Solids

##  1 Introduction

Leveraging the coupling between internal degrees of freedom and material
strains to open up new avenues for structural instabilities has been receiving
much recent attention in the mechanics and advanced functional material
communities [1, 2]. The magneto-elastic coupling is particularly interesting
since it can be exploited to actuate soft magnetic structures without any
mechanical contact. Magneto-rheological elastomers (MREs) comprise a composite
of magnetized (metal) particles and a soft elastic matrix that can respond to
an applied magnetic field. MREs have long been studied, from both mathematical
modeling and application perspectives [2, 3]. For instance, tunable mechanical
functionalities of magnetic structures have been used in the applications such
as sensing medical devices [1], engine mounts [4], soft robotics [5],
actuators [6], and minimally invasive procedures [7].

The mechanical responses of MREs to an applied field can be categorized into
the following (i)-(v) cases, depending on their ingredients, such as the
coercivity, which is the necessary field strength to erase the magnetization.
(i) Super-para-magnetic and soft-ferromagnetic materials have large magnetic
susceptibility, enabling the temporary magnetization by an external magnetic
field [8, 9, 10, 11, 12, 13, 14]; a super-para-magnetic material loses its
magnetization, once the external field is removed. As a result, the
magnetization of these materials depends strongly on the external fields. (ii)
Hard-ferromagnetic materials have a high coercivity. An external field does
not change the magnetization within the range of practical actuation. Hence,
their magnetization is independent of external fields. (iii) Chains of hard-
ferromagnetic spheres are purely magnetic systems without elasticity [15, 16,
17]. The complex internal interaction among beads yields a variety of
equilibrium configurations. Interestingly, these chains posses an effective
bending stiffness [16]. (iv) MREs embedded with soft-ferromagnetic particles
or soft magnetic elastomer [18, 4, 19, 1, 20, 21, 22, 23, 24, 25, 26, 27]
deform upon the application of an external magnetic field due to elasto-
magnetic couplings. The magnetic particles with low coercivity dispersed in a
soft matrix align themselves with the applied field, thereby forming particle
chains. To enhance the elasto-magnetic coupling, it is common to apply the
magnetic field to MRE during the curing processes. (v) MREs embedded with
hard-ferromagnetic particles or hard magnetic elastomer have both high
coercivity and flexibility in shape-programming [28, 29, 30, 31, 32, 33, 34,
35, 7]. Once they are saturated magnetically, they retain permanent
magnetization. For the remainder of this paper, we will focus on this latter
type (v).

The high coercivity of hard MREs enables the design of highly functional
mechanical systems, such as auxetic metamaterials [29], programmable materials
[33], micro-swimmers [36, 37, 38], micro or soft-robotics [5, 39, 6] and
haptic devices [40]. Motivated by these emerging applications of hard MREs, it
is timely to formulate structural theories for elementary building blocks made
of hard MREs to serve as predictive tools in the design process and aid in
subsequent analysis.

During the past decade, there have been significant developments in the
fundamental theory of hard magnetic elastomers. In a recent pioneering work by
Zhao et al. [30], a nonlinear theory for finite deformation of 3D (bulk) hard
MRE was developed, based on the nonlinear elasticity complemented by the
magneto-mechanical constitutive laws. In this framework, the Helmholtz free
energy density function comprises elastic (neo-Hookean) and magneto-elastic
parts. The relationship between the induced magnetic flux density in the
material and the external field strength is assumed to be linear, with a
permeability close to that of vacuum. The free energy density for the magneto-
elastic part is modeled such that the magnetic moment is embedded in a soft
matrix. A simulation framework using the finite element method developed by
the same authors, based on nonlinear elasticity in 3D, was shown to be in
quantitative agreement with their experimental results [30].

Based on the framework reported in Ref. [30], theories for hard magnetic
linear beams and elastica have been derived and validated against experiments
under either constant magnetic field [31] or a field with constant gradient
[35]. Furthermore, it was shown how a hard magnetic beam (or elastica) can
buckle under an applied field, analogous to an elastica under compression [41,
28, 32, 31, 35]. Methodologically, the derivations of all of the above reduced
models for hard magnetic linear beam and elastica always boils down to the
formal procedure of dimensional reduction [42, 43], which is summarized next.
More specifically, the total energy functional of a beam is modeled as Ref.
[30], with elastic and magneto-elastic parts. The integrand is expanded with
respect to the thickness of the cross-section, only retaining the leading
order terms. Upon integration over the cross-section, a one-dimensional
(1D)-reduced energy is derived. The reduced framework is consistent with the
analogous procedure based on force and moment balance [28], which was also
implemented in the study of the inverse problem [32]. A similar strategy based
on dimensional reduction has also been employed to predict the axisymmetric
deformation of pressurized hard magnetic shells [34].

While theoretical frameworks for one-dimensional hard magnetic slender
structures are now well established for planar 2D deformations (e.g., beams
and elastica), the modeling of hard magnetic rods, with natural curvatures,
undergoing 3D deformation is still lacking. However, the extension of the
formulation in 2D toward 3D deformation is not straightforward because of the
complex interplay between elasticity and geometry [42, 43, 44, 45, 46, 47].
The precise modeling of a hard magnetic elastomeric rod is crucial to
simulate, for example, micro-magnetic-swimmers [36, 37, 38] in a 3D complex
channel, and endovascular probes [7].

In this paper, we establish a theory for the deformation of hard magnetic
elastomeric rods, combining dimensional reduction and simulations, both of
which are validated against precision experiments. For our theory, we start
from the total energy functional comprising the sum of the Hookean elastic
energy [43] and the elasto-magnetic free-energy proposed in Ref. [30]. We then
employ the centerline-based kinematics utilizing the Cosserat frame and
Darboux vector [47]. Along the (circular) cross-section, we span the (local)
polar coordinate around the centerline position to parametrize the material
points of the hard magnetic rod. Subsequently, we expand the integrand of the
elasto-magnetic energy for the rod diameter, retaining only the leading order
terms. After integrating the integrand over the cross-section, we obtain the
reduced elasto-magnetic functional. Based on this reduced elasto-magnetic
functional, we then apply the principle of virtual work to arrive at the
governing nonlinear ordinary differential equations (strong form of the
problem). The derived equilibrium equations are analogous to the Kirchhoff rod
equations, but with additional elasto-magnetic terms; hence, we call them the
magnetic Kirchhoff rod equations. The nonlinear deformations based on the
magnetic Kirchhoff rod equations are computed with the discrete simulation
method [48, 49, 50, 51, 52, 53, 54, 55, 56, 57]. Both the dimensional
reduction procedure and the simulations are then validated against precision
experiments. Since our theory is based on the assumption that the magnetic
torque is induced by the misalignment between the magnetization and the
applied field, the long-range dipole interactions are not modeled. Still, we
find excellent agreement between theory and experiments as long as the long-
range dipole interactions are negligible. We quantify when and how the long-
range interactions become important in the experiments.

Our paper is organized as follows. In Sec. 2, we define the problem. The
derivation of the magnetic Kirchhoff equations follows in Sec. 3, where we
also show that our theory reproduces the previous results of planar
deformation of hard magnetic elastica and beams. To compute 3D geometrically
nonlinear deformation, we numerically solve the derived set of equations as
detailed in Sec. 4. In Sec. 5, the experimental fabrication, apparatus, and
protocols are presented. We report the comparison between theory and
experiment in Sec. 6. We further discuss the limitation of our theory
systematically in Sec. 7. In Sec. 8, we discuss and summarize our findings.

##  2 Definition of the problem

In this section, we define the problem and identify the relevant variables. We
consider a naturally curved and twisted (inextensible) rod of circular cross-
section with diameter $d$ (area $A=\pi d^{2}/4$) and total length $L$ (Fig.
1(a)) [44, 45, 43, 46, 47]. The goal of our study was to derive a reduced-
order theoretical framework to describe the deformation of the centerline,
${\bm{r}}(s)$, of this hard magnetic rod, upon magnetic actuation. The arc-
length of the centerline is $s$ $(0\leq s\leq L)$. In particular, we seek to
derive the set of governing ordinary differential equations (ODEs); i.e., the
strong form of the problem depicting force and moment balance of the system.
This section details the kinematics and constitutive relations of our system,
in Sec. 2.1 and 2.2, respectively.

Figure 1: Definition of the problem. (a) A naturally curved and twisted rod is
magnetized permanently along ${\bm{B}}^{\rm r}$. (b) Subsequently, when a
magnetic field ${\bm{B}}^{\rm a}$ is applied, the rod deforms due to the
induced magnetic force and torque. The centerline of the rod is described
uniquely by the rotation of the Cosserat frame
$(\hat{\bm{d}}_{1},\hat{\bm{d}}_{2},\hat{\bm{d}}_{3})$ along the arc-length
$s$. (c) The cross section of the rod is assumed circular. Any material point
of the 3D rod is identified by the cylindrical coordinates $(s,\rho,\Theta)$,
where polar coordinates $(\rho,\Theta)$ are set on the
$\hat{\bm{d}}_{1}-\hat{\bm{d}}_{2}$ plane for each value of $s$.

###  2.1 Kinematics

Our rod is slender such that the diameter is significantly smaller than its
length; $d/L\ll 1$. Throughout this paper, we denote the quantities in the
reference configuration with $\accentset{\circ}{({\cdot})}$. For example, we
write the centerline position in the reference configuration as
$\mathring{\bm{r}}(s)$. The 2D orthogonal basis
$\hat{\bm{d}}_{1}$-$\hat{\bm{d}}_{2}$ is defined on the cross-section normal
to the tangent vector ${\hat{\bm{d}}}_{3}={\bm{r}}^{\prime}(s)$, where
$(\cdot)^{\prime}$ represents differentiation with respect to $s$. With the
help of the Cosserat frame basis ${\hat{\bm{d}}}_{a}~{}(a=1,2,3)$ defined
here, one can uniquely identify the configuration of the rod centerline (other
than rigid body translations and rotations). The rotation of the Cosserat
frame can be represented using the Darboux vector
${\bm{\Omega}}=\Omega_{a}\hat{\bm{d}}_{a}$ as

$\displaystyle{\hat{\bm{d}}}_{a}^{\prime}$ $\displaystyle=$
$\displaystyle{\bm{\Omega}}\times{\hat{\bm{d}}}_{a},$ (1)

where $\Omega_{a}$ are the rotation rates of the Cosserat frame around
$\hat{\bm{d}}_{a}$ [44, 45, 43, 46, 47]. We can then define the slenderness
parameters as $\mathring{\Omega}_{a}d$.

The 3D reference configuration of the rod can be described uniquely by
introducing the cylindrical coordinates $\xi_{\mu}=s,\rho,\Theta$
($0\leq\rho\leq d/2$, $0\leq\Theta\leq 2\pi$), whose origin corresponds to
${\bm{r}}(s)$ as shown in Fig. 1 (c). For example, the position vector of the
3D rod is

$\displaystyle{\bm{R}}(s,\rho,\Theta)={\bm{r}}(s)+\rho\cos\Theta~{}{\hat{\bm{d}}}_{1}(s)+\rho\sin\Theta~{}{\hat{\bm{d}}}_{2}(s).$
(2)

Once the geometry is set up, we can obtain information about the deformation
of the 3D body, which is characterized by the deformation gradient ${\bf F}$,
defined as

$\displaystyle{\bf
F}\equiv\frac{\partial{\bm{R}}}{\partial\mathring{\bm{R}}}={\bm{t}}_{\mu}\otimes\mathring{\bm{t}}_{\mu}^{\dagger}.$
(3)

The tangent vector along the curvilinear coordinate, ${\bm{t}}_{\mu}$, is
expressed as ${\bm{t}}_{\mu}=\partial{\bm{R}}/\partial\xi_{\mu}$, where
$\xi_{\mu}=s,\rho,\Theta$ for $\mu=1,2,3$, respectively. The reciprocal basis
of ${\bm{t}}_{\mu}$ denoted by, ${\bm{t}}_{\mu}^{\dagger}$, is defined as
${\bm{t}}_{\mu}^{\dagger}\cdot{\bm{t}}_{\nu}=\delta_{\mu\nu}$ with the
Kronecker delta $\delta_{\mu\nu}$, which is equivalently written as
${\bm{t}}_{\mu}^{\dagger}\equiv\partial\xi_{\mu}/\partial{\bm{R}}$. By using
the kinematic relations, we can calculate the tangent vectors. Note, for
example, that when $\mathring{\bm{t}}_{\mu}^{\dagger}$ and ${\bm{t}}_{\mu}$
are expressed in terms of the Cartesian basis $\hat{\bm{e}}_{i}(i=x,y,z)$,
each component of $\mathbf{F}$ is calculated as
$F_{ij}=({\bm{t}}_{\mu}\cdot\hat{\bm{e}}_{i})(\mathring{\bm{t}}_{\mu}^{\dagger}\cdot\hat{\bm{e}}_{j})$.
Throughout this paper, $a,b,c$ and $i,j,k$ are indices of the Cosserat frame
and the Cartesian basis, respectively, if not specified.

###  2.2 Constitutive relation and Elasto-magnetic potential

We assume that the rod is made of a Hookean material, where the internal
stress is linearly proportional to the strain [43]. The rod is permanently
magnetized in the reference configuration with a magnetic flux density
${\bm{B}}^{\rm r}$. We investigate the deformation of the rod upon the
application of an external applied magnetic flux density vector ${\bm{B}}^{\rm
a}({\bm{r}})$ (Fig. 1(b)). Note that the applied field $\bm{B}^{\rm a}$ is set
to vary in space. We will specialize our theory in the cases of either
homogeneous $\partial\bm{B}^{\rm a}/\partial r_{i}=0$ or inhomogenous fields
$\partial\bm{B}^{\rm a}/\partial r_{i}\neq 0$, making use of the precision
experimental framework introduced in section 5.2.

As a starting point for the derivation of our Kirchhoff-like theory for hard
magnetic rods, we will make use of the Helmholtz free energy for ideal hard-
magnetic soft materials that was introduced recently in Refs. [29, 30]; which
we review next. Once the rod is magnetized, all magnetized moments are
parallel with $\bm{B}^{\rm r}$ (fully saturated). Hence, it is reasonable to
assume that the magnetic permeability of the hard magnetic material is close
to that of vacuum $\mu_{0}$. The magnetized moments will not affect the
(surrounding) applied field $\bm{B}^{\rm a}$ upon actuation. Besides, the
permeability of the elastomer that we will use in the experiments (section
5.2) is nearly the same as that of the vacuum $\mu_{0}$. The magnetization
density vector per unit volume in the deformed configuration is modeled as the
“rotation” of ${\bm{B}}^{\rm r}$ as $\bm{\mathcal{M}}_{\rm 3D}={\bf
F}{\bm{B}}^{\rm r}/\mu_{0}$. Thus, the total elasto-magnetic energy for a hard
magnetic rod can be simplified as the work to align the residual magnetic
moment of the material with ${\bm{B}}^{a}$:

$\displaystyle{E}_{\rm mag}=-\iiint{\bm{{\mathcal{M}}}}_{\rm 3D}{\bm{B}}^{\rm
a}\rho d\rho d\Theta ds=-\iiint\frac{1}{\mu_{0}}({\bf F}{\bm{B}}^{\rm
r})\cdot{\bm{B}}^{\rm a}\rho d\rho d\Theta ds\,.$ (4)

Note here that we do not consider the internal stress or long-range forces
induced by a dipole-dipole interaction [15, 58, 16, 59]. The validity of this
simplifying assumption will be evaluated thoroughly in Sec. 7. In the
following, by integrating Eq. (4) over the cross-section, we seek to obtain
the reduced expression for the elasto-magnetic energy density of the rod per a
unit length, up to the leading order in the slenderness.

##  3 Theory of a Kirchhoff-like equation for hard magnetic elastic rods

In this section, we derive the equilibrium equations (ODEs) describing the
geometrically nonlinear deformation of hard magnetic rods. The main results of
this section are Eqs. (20) and (21), that is the magnetic Kirchhoff rod
equations. To arrive at this set of ODEs, in Sec. 3.1, we will follow the
procedure of dimensional reduction. In Sec. 3.2, we review the derivation of
the non-magnetic Kirchhoff rod equations based on the principle of virtual
work. We derive the elasto-magnetic force $\bm{p}_{\rm mag}$ and torque
$\bm{q}_{\rm mag}$ through variational calculus in Sec. 3.3, thus obtaining
the magnetic Kirchhoff rod equations for hard magnetic rods. In Sec. 3.4, we
show that our framework reproduces previous works on the planar deformation of
hard magnetic beams and elastica [28, 32, 31, 35], whereas in Sec. 6 we will
employ our framework to study cases where the deformation is three-
dimensional.

###  3.1 Reduced elasto-magnetic and elastic energy

Integrating Eq. (4) over the cross-section will yield the 1D reduced energy of
the rod system, to leading order of slenderness $O(1)$
($\mathring{\Omega}_{a}d,\Omega_{a}d\ll 1$). To proceed, we first expand the
deformation gradient as
$\mathbf{F}(s,\rho,\Theta)=\mathbf{D}(s)+O(\mathring{\Omega}_{a}d,\Omega_{a}d)$,
where

$\displaystyle\mathbf{D}(s)\equiv{\hat{\bm{d}}}_{a}(s)\otimes\mathring{\hat{\bm{d}}}_{a}(s)$
(5)

is the reduced deformation gradient, mapping the undeformed tangent space to
the deformed tangent space along the centerline of the rod. The deformation
gradient is now reduced to the tensor product between deformed and undeformed
directors, as given by Eqs. (3) and (5). Given that $\mathbf{D}$ is not a
function of $\rho$ nor $\Theta$, the integral along the cross-section is
replaced by the constant cross-section area $A$. Hence, from Eq. (4), we
obtain the reduced elasto-magnetic energy functional as

$\displaystyle{E}_{\rm mag}=-\int_{0}^{L}\bm{\mathcal{M}}(s)\cdot{\bm{B}}^{\rm
a}({\bm{r}}(s))ds\,,$ (6)

where we have introduced the magnetization density vector per a unit length
$\bm{\mathcal{M}}$ (by contrast to the 3D counterpart,
$\bm{{{\mathcal{M}}}}_{\rm 3D}$, introduced above) defined as

$\displaystyle\bm{\mathcal{M}}(s)\equiv\frac{{A}}{\mu_{0}}({\bf
D}(s){\bm{B}}^{\rm r}).$ (7)

It is important to note that, when the applied field ${\bm{B}}^{\rm a}$ is not
uniform, ${\bm{B}}^{\rm a}$ at the material point $\bm{R}(s,\rho,\Theta)$ is
represented by the value at the centerline position ${\bm{r}}(s)$ as
$\bm{B}^{\rm a}(\bm{R}(s,\rho,\Theta))\simeq\bm{B}^{\rm a}(\bm{r}(s))$. The
spatial variation of ${\bm{B}}^{\rm a}$ along the cross section is at higher
order of slenderness $O(\mathring{\Omega}_{a}d,\Omega_{a}d)$ as $\bm{B}^{\rm
a}(\bm{R})=\bm{B}^{\rm a}(\bm{r})+O(\mathring{\Omega}_{a}d,\Omega_{a}d)$, and
can, therefore, be neglected.

We now want to derive the equilibrium equations, augmenting the non-magnetic
Kirchhoff rod equations. In the absence of magnetic effects, we recall that,
from the classic literature [43, 44, 45, 46, 47], the total elastic energy for
a naturally curved and twisted elastic (Hookean) rod in the absence of a
magnetic field is given by

$\displaystyle{E}_{\rm
el}=\int_{0}^{L}\frac{EI}{2}\Bigl{\\{}\Bigl{(}\Omega_{1}-\mathring{\Omega}_{1}\Bigr{)}^{2}+\Bigl{(}\Omega_{2}-\mathring{\Omega}_{2}\Bigr{)}^{2}\Bigr{\\}}+\frac{GJ}{2}\Bigl{(}\Omega_{3}-\mathring{\Omega}_{3}\Bigr{)}^{2}ds,$
(8)

where $EI$ and $GJ$ are the bending and twist moduli, respectively, with $E$
as the Young modulus, $G=E/\\{2(1+\nu)\\}$ as the shear modulus with the
Poisson ratio $\nu$, $I=\pi d^{4}/64$ as the second moment of area and $J=\pi
d^{4}/32$ as the torsional constant.

The set of equations of mechanical equilibrium is derived through the
principle of virtual work as

$\displaystyle\delta E_{\rm el}+\delta{E}_{\rm mag}=\delta W_{\rm ex},$ (9)

where $W_{\rm ex}$ is the work done by the external forces and torques,
$({\bm{F}}_{\rm ex}(0),{\bm{M}}_{\rm ex}(0))$ and $({\bm{F}}_{\rm
ex}(L),{\bm{M}}_{\rm ex}(L))$, applied at $s=0$ and $s=L$, respectively,
together with the external forces $\bm{p}(s)$ and torques $\bm{q}(s)$ per unit
length (e.g., gravity, viscous force and torque).

The magnetic Kirchhoff rod equations can be derived directly from Eq. (9),
which will be tackled in Sec. 3.3. Before describing our theory for magnetic
rods, we review the well-established and classic (non-magnetic) Kirchhoff rod
equations. Note that $\delta E_{\rm el}=\delta W_{\rm ex}$ corresponds to the
classic case in the absence of the magnetic field. Given that $E_{\rm mag}$ is
additive in the total energy, the variation of $E_{\rm mag}$ alone will yield
the expressions of the elasto-magnetic force $\bm{p}_{\rm mag}$ and torque
$\bm{q}_{\rm mag}$ per unit length, which must be added to the force and
moment balance equations as an external force and torque, respectively.

###  3.2 Review of the (non-magnetic) Kirchhoff rod equations

In the previous section, we have derived the reduced expressions for the
elasto-magnetic potential. In the absence of the magnetic field $E_{\rm
mag}=0$, the variations of $E_{\rm el}$ and $W_{\rm ex}$ follow the classic
(non-magnetic) Kirchhoff rod equations. In this section, for completeness, we
review their classic derivation. The Kirchhoff rod equations can be derived by
considering the restricted variation of the total energy functional such that
the centerline does not stretch: $\delta(ds)=0$ (i.e., the differentiation
with respect to $s$, $d/ds$, and the variation, $\delta$, commute). A detailed
account of this derivation can be found in Ref. [45]; in what follows, we
provide just a brief sketch.

The variation of the elastic energy functional $\delta E_{\rm el}$ is written
using the variation of $\Omega_{a}$, $\delta\Omega_{a}$:

$\displaystyle\delta E_{\rm el}$ $\displaystyle=$
$\displaystyle\int_{0}^{L}EI(\Omega_{1}-\mathring{\Omega}_{1})\delta\Omega_{1}+EI(\Omega_{2}-\mathring{\Omega}_{2})\delta\Omega_{2}+GJ(\Omega_{3}-\mathring{\Omega}_{3})\delta\Omega_{3}ds.$
(10)

To rewrite Eq. (10) with respect to the angular parameters and infinitesimal
angles $\bm{\epsilon}$ [45], the infinitesimal rotation vector of the Cosserat
frame ${\bm{\epsilon}}$ is introduced as
${\bm{\epsilon}}\equiv\delta\chi_{a}\hat{\bm{d}}_{a}$ whose components are
expressed as
$\delta\chi_{1}\equiv(\delta\hat{\bm{d}}_{2})\cdot\hat{\bm{d}}_{3}=-(\delta{\bm{r}})^{\prime}\cdot\hat{\bm{d}}_{2}$,
$\delta\chi_{2}\equiv(\delta\hat{\bm{d}}_{3})\cdot\hat{\bm{d}}_{1}=(\delta{\bm{r}})^{\prime}\cdot\hat{\bm{d}}_{1}$,
and
$\delta\chi_{3}\equiv(\delta\hat{\bm{d}}_{1})\cdot\hat{\bm{d}}_{2}=-(\delta\hat{\bm{d}}_{2})\cdot\hat{\bm{d}}_{1}$.
The end-forces and moments are applied at $s=0$ and $s=L$ as ${\bm{F}}_{\rm
ex}(0),{\bm{M}}_{\rm ex}(0),{\bm{F}}_{\rm ex}(L)$, and ${\bm{M}}_{\rm ex}(L)$,
respectively. The rod is subjected to the external force and torque per unit
length ${\bm{p}}(s)$ and ${\bm{q}}(s)$. The variation of the work done by the
external forces and torques (in the absence of magnetic fields) is

$\displaystyle\delta W_{\rm ex}$ $\displaystyle=$ $\displaystyle{\bm{F}}_{\rm
ex}(0)\cdot\delta{\bm{r}}(0)+{\bm{F}}_{\rm
ex}(L)\cdot\delta{\bm{r}}(L)+{\bm{M}}_{\rm
ex}(0)\cdot{\bm{\epsilon}}(0)+{\bm{M}}_{\rm ex}(L)\cdot{\bm{\epsilon}}(L)$
(11)
$\displaystyle+\int_{0}^{L}\left\\{{\bm{p}}(s)\cdot\delta{\bm{r}}(s)+{\bm{q}}(s)\cdot{\bm{\epsilon}}(s)\right\\}ds.$

The internal and external virtual works expressed with respect to
$\delta{\bm{r}}$ and ${\bm{\epsilon}}$ are then equated as $\delta E_{\rm
el}=\delta W_{\rm ex}$. Collecting the terms associated with $\delta\bm{r}(s)$
and $\bm{\epsilon}(s)$ yields the (non-magnetic) Kirchhoff rod equations

$\displaystyle{\bm{F}}^{\prime}(s)+{\bm{p}}(s)=0,$ (12)
$\displaystyle{\bm{M}}^{\prime}(s)+\hat{\bm{d}}_{3}(s)\times{\bm{F}}(s)+{\bm{q}}(s)=0,$
(13)

with the internal moment vector
${\bm{M}}(s)=M_{a}\hat{\bm{d}}_{a}=EI\\{(\Omega_{1}-\mathring{\Omega}_{1})\hat{\bm{d}}_{1}+(\Omega_{2}-\mathring{\Omega}_{2})\hat{\bm{d}}_{2}\\}+GJ(\Omega_{3}-\mathring{\Omega}_{3})\hat{\bm{d}}_{3}$,
and $\bm{F}(s)$ is the internal force acting over the cross section at $s$.
Equations (12) and (13) are the equilibrium equations of forces and moments
for the rod, respectively. The boundary conditions can be also derived from
variational calculus [45] as

$\displaystyle{\bm{M}}(0)=-{\bm{M}}_{\rm
ex}(0),~{}~{}{\bm{M}}(L)={\bm{M}}_{\rm ex}(L),~{}~{}{\bm{F}}(0)=-{\bm{F}}_{\rm
ex}(0),~{}~{}{\bm{F}}(L)={\bm{F}}_{\rm ex}(L).$ (14)

Next, we will tackle the derivation of the magnetic Kirchhoff rod equations
based on the same procedure sketched above, while the elastic terms Eqs. (12)
and (13) remain unchanged. The last task is to compute the variation of the
elasto-magnetic potential, $E_{\rm mag}$, with respect to $\delta\bm{r}(s)$
and $\bm{\epsilon}(s)$, $\delta E_{\rm mag}$, which is addressed next.

###  3.3 Derivation of the magnetic Kirchhoff rod equations using the
principle of virtual work

We now set out to derive the magnetic Kirchhoff rod equations. First, we will
compute the variation of the elasto-magnetic potential, $\delta E_{\rm mag}$.
Equating the total internal and external virtual works as Eq. (9), we will
obtain the magnetic Kirchhoff rod equations. Given that $\delta E_{\rm
el},\delta W_{\rm ex}$, and $\delta E_{\rm mag}$ are additive in Eq. (9), the
variation of the elasto-magnetic potential $\delta E_{\rm mag}$ provides the
elasto-magnetic force $\bm{p}_{\rm mag}$ and torque $\bm{q}_{\rm mag}$
explicitly, which will be readily added to Eqs. (12) and (13) to capture
magnetic effects.

Before calculating $\delta E_{\rm mag}$, we compute the variation of the
Cosserat frame basis $\hat{\bm{d}}_{a}$, $\delta\hat{\bm{d}}_{a}$. Given that
${\bm{\epsilon}}$ represents the infinitesimal rotation around
$\hat{\bm{d}}_{a}$, the variation of $\hat{\bm{d}}_{a}$ is
$\delta\hat{\bm{d}}_{a}={\bm{\epsilon}}\times\hat{\bm{d}}_{a}=-\varepsilon_{abc}\delta\chi_{b}\hat{\bm{d}}_{c}$,
where $\varepsilon_{abc}$ is the Eddington epsilon. We aim to express $\delta
E_{\rm mag}$ with respect to $\delta{\bm{r}}$ and ${\bm{\epsilon}}$. Note that
$\delta\bm{B}^{\rm r}$ and $\delta\mathring{\hat{\bm{d}}}_{a}$ are zero, since
they do not depend on the state variables. We first determine the variation of
the reduced deformation gradient ${\bf D}$ for both $\delta{\bm{r}}$ and
${\bm{\epsilon}}$ as

$\displaystyle\delta{{\bf D}}$ $\displaystyle=$
$\displaystyle\delta(\hat{\bm{d}}_{a}\otimes\mathring{\hat{\bm{d}}}_{a})=(\delta\hat{\bm{d}}_{a})\otimes\mathring{\hat{\bm{d}}}_{a}$
(15) $\displaystyle=$
$\displaystyle\\{(\delta{\bm{r}}^{\prime}\cdot\hat{\bm{d}}_{2})\\}\mathbf{\Delta}_{1}-\\{(\delta{\bm{r}}^{\prime}\cdot\hat{\bm{d}}_{1})\\}\mathbf{\Delta}_{2}-\\{({\bm{\epsilon}}\cdot\hat{\bm{d}}_{3})\\}\mathbf{\Delta}_{3}$
$\displaystyle=$
$\displaystyle-\\{({\bm{\epsilon}}\cdot\hat{\bm{d}}_{a})\\}\mathbf{\Delta}_{a}\,,$

where we have defined the new tensor $\mathbf{\Delta}_{a}~{}(a=1,2,3)$
associated with the magnetic torque as

$\displaystyle\mathbf{\Delta}_{a}(s)\equiv\varepsilon_{abc}\hat{\bm{d}}_{b}(s)\otimes\mathring{\hat{\bm{d}}}_{c}(s).$
(16)

To obtain Eq. (15), we have used the orthogonality of the basis
$\hat{\bm{d}}_{b}\times\hat{\bm{d}}_{c}=\varepsilon_{abc}\hat{\bm{d}}_{a}$. We
also used the fact that, for an arbitrary vector ${\bm{a}}={\bm{a}}(s)$,
$\delta{\bm{r}}^{\prime}\cdot{\bm{a}}=\delta(\hat{\bm{d}}_{3})\cdot{\bm{a}}=({\bm{\epsilon}}\times\hat{\bm{d}}_{3})\cdot{\bm{a}}={\bm{\epsilon}}\cdot(\hat{\bm{d}}_{3}\times{\bm{a}})$
holds, where we set ${\bm{a}}=\hat{\bm{d}}_{1}$ or
${\bm{a}}=\hat{\bm{d}}_{2}$. Given that the variation of ${\bm{B}}^{\rm
a}(\bm{r})$ acts as $\delta{\bm{B}}^{\rm a}=(\partial{\bm{B}}^{\rm a}/\partial
r_{i})\delta r_{i}\equiv(\nabla{\bm{B}}^{\rm a})\delta{\bm{r}}$, we compute
the variation of $E_{\rm mag}$, $\delta E_{\rm mag}$ as

$\displaystyle\delta E_{\rm mag}$ $\displaystyle=$
$\displaystyle-\int_{0}^{L}[\delta{\bm{r}}\cdot{\bm{p}}_{\rm
mag}(s)+{\bm{\epsilon}}\cdot{\bm{q}}_{\rm mag}(s)]ds,$ (17)
$\displaystyle{\bm{p}}_{\rm mag}(s)$ $\displaystyle\equiv$
$\displaystyle\frac{A}{\mu_{0}}({\bf D}(s){\bm{B}}^{\rm
r})\cdot(\nabla{\bm{B}}^{\rm a})=\bm{\mathcal{M}}\cdot(\nabla{\bm{B}}^{\rm
a}),$ (18) $\displaystyle{\bm{q}}_{\rm mag}(s)$ $\displaystyle\equiv$
$\displaystyle-\frac{A}{\mu_{0}}\left\\{(\mathbf{\Delta}_{a}(s){\bm{B}}^{\rm
r})\cdot{\bm{B}}^{\rm
a}\right\\}\hat{\bm{d}}_{a}(s)=\bm{\mathcal{M}}\times{\bm{B}}^{\rm a}.$ (19)

We have introduced ${\bm{p}}_{\rm mag}$ and ${\bm{q}}_{\rm mag}$ as the
elasto-magnetic force and torque, respectively. When the applied field is
homogeneous $\partial\bm{B}^{\rm a}/\partial r_{i}=0$, the magnetic torque
$\bm{q}_{\rm mag}\neq 0$ drives the deformation (i.e., $\bm{p}_{\rm mag}=0$),
while both $\bm{p}_{\rm mag}$ and $\bm{q}_{\rm mag}$ contribute to the
deformation if the applied field is inhomogeneous.

Adding the elasto-magnetic force and torque terms from Eqs. (18) and (19) into
the classic (non-magnetic) Kirchhoff rod equations (12) and (13), yields the
force and moment balance equations:

$\displaystyle{\bm{F}}^{\prime}(s)+{\bm{p}}(s)+{\bm{p}}_{\rm mag}(s)=0,$ (20)
$\displaystyle{\bm{M}}^{\prime}(s)+\hat{\bm{d}}_{3}(s)\times{\bm{F}}(s)+{\bm{q}}(s)+{\bm{q}}_{\rm
mag}(s)=0,$ (21)

where the magnetic force ${\bm{p}}_{\rm mag}$ and torque ${\bm{q}}_{\rm mag}$
were defined in Eqs. (18) and (19), respectively. Note that the boundary
conditions specified in Eq. (14) remain unchanged in the elasto-magnetic case
because they do not appear in Eq. (17). We refer to the set of Eqs. (20) and
(21) as the magnetic Kirchhoff rod equations, where the coupling between
elasticity and magnetism is captured by ${\bm{p}}_{\rm mag}$ and
${\bm{q}}_{\rm mag}$.

###  3.4 Reduction of the magnetic Kirchhoff rod equations to planar
deformations

In this subsection, for verification purposes, we show that the magnetic
Kirchhoff equations derived above can reproduce existing results in the
literature for the hard magnetic beams and elastica [31, 35], loaded under a
constant or constant gradient magnetic fields. In Ref. [31], a theoretical
model under constant field was derived within the continuum framework and
validated against experiments. The governing equations for a more general
description that can tackle either under constant or constant gradient fields
were developed through dimensional reduction and validated experimentally in
Ref. [35].

Consider a naturally straight rod clamped at $s=0$ (${\bm{r}}(0)={\bm{0}}$),
while the other end at $s=L$ is set to be force and momentum free. The rod is
free from any external forces and torques other than magnetic fields, that is
${\bm{p}}={\bm{q}}=0$. We will derive the governing equation for the bending
angle $\theta(s)$ in the $x$-$y$ plane such that $\theta(s)=0$ holds in the
reference configuration. The constitutive law is
${\bm{M}}(s)=EI\theta^{\prime}(s)\hat{\bm{e}}_{z}$. We will study two cases: a
naturally straight hard magnetic rod loaded under a constant or constant-
gradient fields, which have been studied in Ref. [31] and Ref. [35],
respectively.

First, we consider the case of a rod under a constant magnetic field
${\bm{B}}^{\rm a}=B^{\rm a}\hat{\bm{e}}_{y}$; see Fig. 2(a). The rod is
magnetized as ${\bm{B}}^{r}=-B^{r}\hat{\bm{e}}_{y}$ along the center-line (
$\mathring{\hat{\bm{d}}}_{1}=\hat{\bm{e}}_{x}$,
$\mathring{\hat{\bm{d}}}_{2}=\hat{\bm{e}}_{z}$,
$\mathring{\hat{\bm{d}}}_{3}=-\hat{\bm{e}}_{y}$) with $B^{r}>0$ and $B^{\rm
a}>0$. In the absence of the applied field, the rod is naturally straight
along $\hat{\bm{e}}_{y}$ as $\mathring{\hat{d}}_{3}=-\hat{\bm{e}}_{y}$. The
Cosserat frame in the reference configuration is given by
$\hat{\bm{d}}_{1}=(\cos\theta,\sin\theta,0)$,
$\hat{\bm{d}}_{2}=\hat{\bm{e}}_{z}$ and
$\hat{\bm{d}}_{3}=(\sin\theta,-\cos\theta,0)$. Since the applied field is
homogeneous, we have ${\bm{p}}_{\rm mag}=0$ and $\bm{F}(s)=0$. To calculate
${\bm{q}}_{\rm mag}$ in Eq. (19), we first need to express $\bm{\mathcal{M}}$
with respect to $\theta$. From Eq. (7), $\bm{\mathcal{M}}$ is simplified as

$\displaystyle\bm{\mathcal{M}}=\frac{AB^{\rm r}}{\mu_{0}}\hat{\bm{d}}_{3}.$
(22)

Substituting ${\bm{F}}(s)=0$ and Eq. (22) into the moment balance equation
(Eq. (21)), we can reproduce the equation for the hard magnetic elastica [31]
as

$\displaystyle EI\theta^{\prime\prime}(s)+\frac{AB^{\rm r}B^{\rm
a}}{\mu_{0}}\sin\theta(s)=0,$ (23)

with the boundary conditions $\theta(0)=0$ and $\theta^{\prime}(L)=0$. As
discussed in Ref. [31], Eq. (23) is mathematically equivalent to the governing
equation of the clamped-free elastica. In this configuration, the magnetic rod
is known to buckle at the critical applied field [31]

$\displaystyle B^{{\rm a}*}_{\rm
bend}=\left(\frac{\pi}{2}\right)^{2}\frac{\mu_{0}EI}{AB^{r}L^{2}}\,,$ (24)

a result that will be used in Sec. 6.2 to connect our 3D analytical results
with the existing 2D results in literature [31].

Figure 2: Reduction of magnetic Kirchhoff rod equations to planar (2D) cases.
(a) Schematics of a hard magnetic rod under constant field that deforms in 2D.
The rod is magnetized along the tangent of the rod $\hat{\bm{d}}_{3}$. (b)
Schematics of a hard magnetic rod under constant gradient field. The rod is
magnetized along $\hat{\bm{d}}_{1}$.

Second, we consider the case of a rod in a constant gradient magnetic field:
${\bm{B}}^{\rm a}({\bm{r}})=by\hat{\bm{e}}_{y}$, where $b$ is a constant. The
rod is clamped along $\hat{\bm{e}}_{x}$ ($y=0$) and magnetized as
${\bm{B}}^{\rm r}=B^{\rm r}\hat{\bm{e}}_{y}$ (Fig. 2 (b)). The Cosserat frame
in the reference configuration is given by
$\mathring{\hat{\bm{d}}}_{1}=\hat{\bm{e}}_{y}$,
$\mathring{\hat{\bm{d}}}_{2}=\hat{\bm{e}}_{z}$, and
$\mathring{\hat{\bm{d}}}_{3}=\hat{\bm{e}}_{x}$. We can write the Cosserat
frame in the deformed configuration as
${\hat{\bm{d}}}_{1}=(-\sin\theta,\cos\theta,0)$,
$\hat{\bm{d}}_{2}=\hat{\bm{e}}_{z}$, and
${\hat{\bm{d}}}_{3}=(\cos\theta,\sin\theta,0)$. Hence, the magnetization
vector is simplified from Eq. (7) as

$\displaystyle\bm{\mathcal{M}}=\frac{AB^{\rm r}}{\mu_{0}}\hat{\bm{d}}_{1}.$
(25)

Given that the rod is in a constant gradient field, $\bm{p}_{\rm mag}$ is non-
zero. From Eqs. (18) and (25), we find $\bm{p}_{\rm mag}(s)=(AB^{\rm
r}b/\mu_{0})\cos\theta(s)~{}\hat{\bm{e}}_{y}$. The force balance in Eq. (20)
can now be integrated as

$\displaystyle{\bm{F}}(s)$ $\displaystyle=$
$\displaystyle{\bm{F}}(L)+\int_{s}^{L}{\bm{p}}_{\rm
mag}(s^{\prime})ds^{\prime}$ (26) $\displaystyle=$
$\displaystyle\frac{AB^{r}b}{\mu_{0}}\hat{\bm{e}}_{y}\int_{s}^{L}\cos\theta(s^{\prime})ds^{\prime}.$

Plugging this result into the moment balance in Eq. (21), we obtain the
following integro-differential equation

$\displaystyle
EI\frac{d^{2}{\theta}}{d{s}^{2}}-\frac{AB^{r}b}{\mu_{0}}\left\\{{y}({s})\sin{\theta}({s})-\cos{\theta}({s})\int_{{s}}^{1}\cos{\theta}(s^{\prime})ds^{\prime}\right\\}=0,$
(27)

whose second and third terms originate from magnetic torque and force,
respectively. Equation (27) has been recently derived and validated
experimentally in Refs. [31, 35], verifying our reduced framework presented in
this section.

We note that the inhomogeneous applied field can be interpreted as a
distributed load on the magnetic rod. Indeed, in the limit of small bending
angle $\theta\ll 1$, Eq. (27) is identical to that of an elastica under
gravity [41]. By analogy with the gravitational case, we can introduce the
magneto-bending length $\ell_{m}$ as

$\displaystyle\ell_{m}\equiv\left(\frac{\mu_{0}EI}{AB^{\rm r}b}\right)^{1/3},$
(28)

which can be regarded as a (persistence) length quantifying the relative
importance of the magnetic and elastic effects. When $\lambda_{m}\equiv
L/\ell_{m}\ll 1$, the gradient of the magnetic field is negligible, while in
the case of $\lambda_{m}\gg 1$, the magnetic gradient plays an important role
in the bending of the magnetic rod.

We have shown that the magnetic Kirchhoff rod equations derived in Sec. 3.3
are able to reproduce existing results reported in Refs. [31, 35] for the
planar (2D) deformation of the magnetic rod. The specialization to the 2D
cases verifies our theory, even if only partially. In the following sections,
we will combine discrete simulation and precision experiments to study non-
planar (3D) deformations of the rod, towards the validation of our theoretical
framework.

##  4 Discrete simulations of magneto-elastic rods

Although our magnetic Kirchhoff rod model reproduces the governing equation of
planar hard magnetic beams and elastica, computing the 3D nonlinear geometric
deformation requires, in general, a numerical approach. Specifically, we
discretize the rod-centerline as a set of connected particles via straight
rigid segments [48, 49, 50, 51, 52, 53, 54, 55, 56, 57]. This discrete method
is widely used to simulate the large deformation of rod-like structures such
as nano-springs [49], DNA [48], helical bacteria [50], flagella [52, 51],
human hairs [53], and gift-wrapping ribbons [54, 55, 56, 57]. Within this
computational framework, the rod centerline is regarded as the chain of
straight segments of natural length $\ell_{0}$ and an equal mass $m=\rho
AL/N$. At each node $i=1,2,...,N$, we assigned a discrete Cosserat frame basis
vectors $(\hat{\bm{d}}_{1,i},\hat{\bm{d}}_{2,i},\hat{\bm{d}}_{3,i})$,
corresponding to $\hat{\bm{d}}_{a}$ in the continuum model. The elastic
bending and twist deformations are represented by the Euler angles between
adjacent discretized Cosserat frame basis vectors [48, 49, 50, 52, 51, 53, 54,
55, 56, 57].

In the present section, we detail the simulation method, with a particular
focus on the discrete versions for the elasto-magnetic forces,
$\bm{\mathcal{P}}$, and torques, $\mathcal{Q}$, in Sec. 4.1. The equations of
motion for the discrete segments follow in Sec. 4.2.

###  4.1 Elasto-magnetic force and torque in the discrete simulation

To compute the nonlinear geometric deformation of our rod-like magnetic
structures, we adopt dynamic simulations, instead of solving the static
equilibrium equations. Static analysis packages for ODEs (e.g., AUTO [60])
might have been able to compute the deformation of our MRE rods. However, it
is not trivial to implement the elasto-magnetic interactions into such
packages due to the geometric non-linearity of our MRE rods. To facilitate
numerical convergence to the mechanical equilibrium, we perform the dynamic
simulation by including the inertial terms in Eqs. (20) and (21).
Specifically, we will derive the equations of motion for ${\bm{r}}$ and the
twist angle $\chi_{3}$ (see Sec. 3.2 for the definition of $\chi_{3}$). We
recall that the expressions of the elasto-magnetic force ${\bm{p}}_{\rm mag}$
and torque ${\bm{q}}_{\rm mag}$, provided in Eqs. (17)-(19), have been derived
from the continuum framework, through the variation $\delta E_{\rm mag}$ with
respect to the set of $(\delta{\bm{r}},{\bm{\epsilon}})$. The set of
infinitesimal bending and twist angles were represented by
${\bm{\epsilon}}=\delta\chi_{a}\hat{\bm{d}}_{a}$ as in Sec. 3.2. In the
discrete simulation framework [48, 49, 50, 51, 52, 53, 54, 55, 56, 57], it is
sufficient to derive the equation of motion of the centerline position and the
twist angle $(\bm{r},\chi_{3})$ alone. In other words, the remaining bending
angles $\chi_{1}$ and $\chi_{2}$ are determined from the centerline positions
$\bm{r}$, as we explain below [45]. To include the elasto-magnetic force and
torque in the equations of motion for $(\bm{r},\chi_{3})$, we will introduce
the discrete version of $({\bm{p}}_{\rm mag},{\bm{q}}_{\rm mag})$, as
$(\bm{{\mathcal{P}}},{{\mathcal{Q}}})$, by rewriting $\delta E_{\rm mag}$ (Eq.
(17)) with respect to the set of $(\bm{r},\chi_{3})$. The discrete version of
the elasto-magnetic torque ${{\mathcal{Q}}}$ (a scalar) will correspond to the
elasto-magnetic torque around $\hat{\bm{d}}_{3}$, while the remaining
components of $\bm{q}_{\rm mag}$ (i.e., ${q}_{\rm mag1}$ and ${q}_{\rm mag2}$)
will be included in the equation of motion for $\bm{r}$, thereby defining
$\bm{{\mathcal{P}}}$, as we detail below.

Expressing $\delta\chi_{1}$ and $\delta\chi_{2}$ from Eq. (17) with the aid of
$\delta\chi_{1}=-(\delta{\bm{r}})^{\prime}\cdot\hat{\bm{d}}_{2}$ and
$\delta\chi_{2}=(\delta{\bm{r}})^{\prime}\cdot\hat{\bm{d}}_{1}$, we find the
following equation after the partial integration of
$(\delta{\bm{r}})^{\prime}$:

$\displaystyle\delta E_{\rm mag}$ $\displaystyle=$
$\displaystyle-\int_{0}^{L}[\delta{\bm{r}}(s)\cdot{\bm{p}}_{\rm
mag}(s)+{\bm{\epsilon}}(s)\cdot{\bm{q}}_{\rm mag}(s)]ds$ (29) $\displaystyle=$
$\displaystyle-\int_{0}^{L}[\delta{\bm{r}}\cdot{\bm{p}}_{\rm
mag}+\delta\chi_{1}q_{\rm mag1}+\delta\chi_{2}q_{\rm
mag2}+\delta\chi_{3}q_{\rm mag3}]ds$ $\displaystyle=$
$\displaystyle-\int_{0}^{L}[\delta{\bm{r}}\cdot{\bm{p}}_{\rm
mag}+\\{-(\delta{\bm{r}})^{\prime}\cdot\hat{\bm{d}}_{2}\\}q_{\rm
mag1}+\\{(\delta{\bm{r}})^{\prime}\cdot\hat{\bm{d}}_{1}\\}q_{\rm
mag2}+\delta\chi_{3}q_{\rm mag3}]ds$ $\displaystyle=$
$\displaystyle-\int_{0}^{L}\left[\delta{\bm{r}}\cdot\bm{{\mathcal{P}}}+\delta\chi_{3}{\mathcal{Q}}\right]ds+\left[\delta{\bm{r}}\cdot\left\\{{q}_{{\rm
mag}1}\hat{\bm{d}}_{2}-{q}_{{\rm
mag}2}\hat{\bm{d}}_{1}\right\\}\right]_{0}^{L},$

where we introduce the $a(=1,2,3)$-th component of the magnetic torque
$q_{{\rm mag}a}={\bm{q}}_{\rm mag}\cdot\hat{\bm{d}}_{a}$. Here, we define the
discrete version of the elasto magnetic force and torque as

$\displaystyle\bm{{\mathcal{P}}}\equiv{\bm{p}}_{\rm mag}+\left({q}_{{\rm
mag}1}\hat{\bm{d}}_{2}-{q}_{{\rm mag}2}\hat{\bm{d}}_{1}\right)^{\prime},$ (30)
$\displaystyle{\mathcal{Q}}\equiv{q}_{{\rm mag}3}\,,$ (31)

which will be used to formulate the discrete version of the equations of
motion. Note that $q_{\rm mag1}$ and $q_{\rm mag2}$ are now included in
$\bm{{\mathcal{P}}}$, appropriately. In the next subsection, we incorporate
the elasto-magnetic force and torque through $\bm{{\mathcal{P}}}$ and
${{\mathcal{Q}}}$ into the equations of motion for the centerline position and
the twist angle $(\bm{r},\chi_{3})$.

###  4.2 Equations of motion used in the discrete simulations

In this subsection, we derive the discrete version of the equation of motion
for the MRE rod by adding $(\bm{{\mathcal{P}}},{\mathcal{Q}})$ into the
equation of motion for the centerline position and the twist angle
$(\bm{r},\chi_{3})$. We start from the magnetic Kirchhoff equations that we
derived earlier, with the sole addition of the inertial terms

$\displaystyle\rho A\frac{\partial{\bm{v}}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{\partial{\bm{F}}}{\partial
s}+{\bm{p}}+\bm{{\mathcal{P}}},$ (32)
$\displaystyle\frac{\partial{\bm{L}}}{\partial t}$ $\displaystyle=$
$\displaystyle\frac{\partial{\bm{M}}}{\partial
s}+\hat{\bm{d}}_{3}\times{\bm{F}}+{\bm{q}}+{\mathcal{Q}}\hat{\bm{d}}_{3},$
(33)

where $({\bm{p}}_{\rm mag},{\bm{q}}_{\rm mag})$ are replaced by
$(\bm{{\mathcal{P}}},{\mathcal{Q}}\hat{\bm{d}}_{3})$, as derived in the
previous section. The left hand sides of Eq. (33) represent the inertia terms.
Here, we define the velocity of the centerline
$\bm{v}(s,t)={\partial{\bm{r}}}/{\partial t}$, the angular momentum per unit
length ${\bm{L}}={I}_{a}\omega_{a}\hat{\bm{d}}_{a}$, and the angular velocity
vector ${\bm{\omega}}(s,t)=\omega_{a}\hat{\bm{d}}_{a}$, with the mass density
per unit volume $\rho$ and the principal moments of inertia of the cross
section $I_{1}=I_{2}=I$ and $I_{3}=I_{1}+I_{2}=2I$ [61]. The angular velocity
vector $\bm{\omega}$ describes the rotation rates of $\hat{\bm{d}}_{a}$ in
time as

$\displaystyle\frac{\partial\hat{\bm{d}}_{a}}{\partial
t}={\bm{\omega}}\times\hat{\bm{d}}_{a}\,.$ (34)

To ensure that the rod relaxes to the mechanical equilibrium, we include
numerical drag forces and torque acting on the rod-centerline per unit length
through ${\bm{p}}=-\gamma_{t}{\bm{v}}$ and
${\bm{q}}=-\gamma_{r}\omega_{3}\hat{\bm{d}}_{3}$. We readily derive the
dynamic equations of twist around $\hat{\bm{d}}_{3}$ by taking the inner
product between the moment balance equation and $\hat{\bm{d}}_{3}$. Then,
introducing the internal elastic forces ${\bm{f}}$ and the axial torque ${T}$
per unit length as

$\displaystyle{\bm{f}}\equiv\frac{\partial{\bm{F}}}{\partial
s},~{}~{}{T}\equiv\frac{\partial{M_{3}}}{\partial
s}+M_{2}\Omega_{1}-M_{1}\Omega_{2},$ (35)

we rewrite Eqs. (32) and (33) as

$\displaystyle\rho A\frac{\partial^{2}{\bm{r}}}{\partial t^{2}}$
$\displaystyle=$ $\displaystyle{\bm{f}}+{\bm{p}}+\bm{{\mathcal{P}}},$ (36)
$\displaystyle I_{3}\frac{\partial^{2}\chi_{3}}{\partial t^{2}}$
$\displaystyle=$ $\displaystyle{T}-\gamma_{r}\frac{\partial\chi_{3}}{\partial
t}+{{\mathcal{Q}}},$ (37)

respectively. We solve Eq. (37) for the discrete versions of the position
vector ${\bm{r}}_{i}$ and the twist angle $\chi_{3,i}$. Choosing the units of
mass, length, and force as $m\ell_{0}$, $\ell_{0}$ and $EI/\ell_{0}^{2}$,
respectively, we non-dimensionalize Eq. (37) as

$\displaystyle\frac{\partial^{2}\tilde{{\bm{r}}}_{i}}{\partial\tilde{t}^{2}}$
$\displaystyle=$
$\displaystyle\tilde{\bm{f}}_{i}+\tilde{\bm{p}}_{i}+\tilde{\bm{{\mathcal{P}}}}_{i},$
(38)
$\displaystyle\tilde{I}_{3}\frac{\partial^{2}\chi_{3,i}}{\partial\tilde{t}^{2}}$
$\displaystyle=$
$\displaystyle\tilde{T}_{i}-\tilde{\gamma}_{r}\frac{\partial\chi_{3,i}}{\partial\tilde{t}}+\tilde{{\mathcal{Q}}}_{i},$
(39)

where quantities with $(\tilde{\cdot})$ represent the corresponding
dimensionless variables.

The numerical scheme for calculating the elastic force ${\bm{f}}_{i}$ and
axial torque ${T}_{i}$ at each node $i$ follows that of Refs. [48, 49, 50, 51,
52]. In short, we employ the Euler angle representation to describe the
configuration of the discrete director frames at each point. The Darboux
vector components in the discretized model, $\Omega_{a,i}~{}(a=1,2,3)$, and
the corresponding discretized elastic energy can be expressed in terms of the
three Euler angles. The variation of the elastic energy $E_{\rm el}$ is
related to the variations of ${\bm{r}}_{i}$ and $\chi_{3,i}$ through the
variations of the corresponding Euler angles, from which we find the bending
and twisting forces acting on each node of the chain [48, 49, 50, 51, 52]. All
values of the dimensionless parameters are taken from the experiments, which
will be presented in Sec. 5, such that our simulations do not contain any free
parameters. We performed the simulation for a sufficiently long time such that
the rod relaxes to the mechanical equilibrium (typically $10^{7}-10^{9}$
discrete time steps), where the number of nodes lies within the range $40\leq
N\leq 120$. We use the two-step Adams-Bashforth method [62] to numerically
integrate the re-scaled dynamical equations in Eq. (39) with non-dimensional
time steps, typically ranging between $10^{-3}-10^{-1}$, to ensure sufficient
numerical accuracy. The director frames at each node are also updated at each
time step, and the corresponding Euler angles are calculated for a new
configuration.

##  5 Experimental fabrication, apparatus, and protocols

We perform three different sets of experiments to validate our theory
presented in Sec. 3: (i) a naturally straight rod under a constant field, (ii)
a helical rod under a constant field, and (iii) a helical rod under a constant
gradient field. We chose these three configurations because previous works
[41, 32, 31, 35] were limited to the case of straight beams, without twist
deformations. For case (i), we still consider a straight rod but also include
twist deformation. Through case (ii), we can individually validate the elasto-
magnetic torque $\bm{q}_{\rm mag}$. In case (iii), both $\bm{q}_{\rm mag}$ and
$\bm{p}_{\rm mag}$ can be validated. In this section, we present the details
of our experiments, whose results will be provided in Sec. 6. To perform the
experiments, in Sec. 5.1, we detail the procedure that we developed to
manufacture our MRE rods. In Sec. 5.2, we present the design of the magnetic
coils used to generate the magnetic field. The detailed protocols for the
cases (i)-(iii) follow in Sec. 5.3.

###  5.1 Fabrication of hard magnetic rods

Our rods were made of a magnetorheological elastomer (MRE), a composite of
hard-magnetic NdPrFeB particles and vinylpolysiloxane (VPS) polymer (Young’s
modulus $E_{\rm vps}=1.16\pm 0.03~{}{\rm MPa}$ and mass density $\rho_{\rm
vps}=1.17\pm 0.26~{}{\rm g/cm^{3}}$). The fabrication of the MRE rods involved
the following steps. First, the non-magnetized NdPrFeB particles (average size
of 5$\mu$m, mass density $\rho_{\rm mag}=7.61~{}{\rm g/cm^{3}}$
MQFP-15-7-20065-089, Magnequench) of weight $m_{\rm mag}$ were mixed with the
VPS base liquid solution (Elite Double 32, Zhermack) of weight $m_{\rm base}$
using a centrifugal mixer (ARE-250, Thinky Corporation). Secondly, the VPS
catalyst (weight $m_{\rm cat}=m_{\rm base}$) was added to the mixed solution
with a ratio of 1:1 in weight to the VPS base. The solution of total weight
$m_{\rm tot}=m_{\rm mag}+m_{\rm base}+m_{\rm cat}$ was then mixed using the
centrifugal mixer for 40 s at 2000 rpm (clockwise), and another 20 s at 2200
rpm (counterclockwise). We further degassed the solution in a vacuum chamber
(absolute pressure below 8 mbar) to remove any air bubbles that could
otherwise compromise the homogeneity of the MRE. The final solution was
injected into the molds (detailed in the next paragraph) using a syringe to
cast either straight or helical rods. We varied the mass concentration ratio
for NdPrFeB particles in the range $c=m_{\rm mag}/m_{\rm tot}=10-30\%$.

Figure 3: Fabrication procedure to manufacture straight and helical rods made
out of a MRE. (a) Schematic diagram of the mold for a straight rod. A straight
acrylic tube (green) is threaded into two parallel supporting plates. (b)
Photograph of the mold used to produce a straight rod. (c) Schematic diagram
of the mold for a helical rod. A rigid cylinder with a helical groove (inner
mold) is 3D printed. (d) Photograph of the mold for a helical rods. (e) The
mold with the cured elastomer is placed into an impulse magnetizer such that
${\bm{B}}^{\rm r}$ in (a) and (c) is aligned with the axis of the coil.
Schematics of a magnetized (f) straight and (g) helical rod. The local
magnetization vectors are depicted by the embedded “magnets,” where the red
and blue ends correspond to the north and south poles, respectively.

The schematic and a photograph of a mold used to fabricate a straight rod are
shown in Fig. 3 (a) and (b), respectively. We injected the mixed solution into
a straight acrylic tube (Plexiglas XT tubes incolore 0A070, Röhm, Switzerland,
inner diameter $2$ or $4$ mm) threaded into two coaxial holes in acrylic
plates. One of the inlets of the tube was designed to have a convex shape
pointing to the direction of the residual flux density. The convex inlet is
utilized to align the magnetization vector ${\bm{B}}^{\rm r}$ precisely.

To fabricate helical rods, we prepared a rigid cylinder of diameter $20$ mm
with a helical half-piped groove (diameter $4.5$ mm) using a 3D printer
(Printer: Form 2, Formlabs, Material: Clear Resin (RS-F2-GPCL-04)). The
schematic and the photograph of the mold with a helical groove are shown in
Fig. 3 (c) and (d), respectively. A flexible tube (diameter $4$ mm, Misumi) is
placed along the groove, injecting the mixed solution into it. The inner mold
and the tube were then positioned inside a 3D printed (outer) cylinder (inner
diameter $12$ mm and outer diameter $16$ mm). The central axis of the helical
groove and that of the outer cylinder were set co-axially. The 3D printed
molds allowed us to control the pitch angle $\psi$ of a helical MRE rod up to
$\psi\lesssim 1.49$ rad. For larger values of $\psi$, we prepared the acrylic
cylinder (without a groove) and tightly spooled a flexible tube around it with
a pitch angle of $\psi=1.51$ rad.

After injection of the mixed solution into either the straight or helical
molds, the curing of the polymer mixture occurred in approximately $20$ min,
at room temperature. The (projected) radius of the helical centerline was
chosen to be $R=10~{}{\rm mm}$ to minimize the deformation of the tube,
allowing us to have a cross-section of the helix close to circular, which
could otherwise deviate from the circular shape for tighter helices due to the
Brazier instability [63].

Before demolding, the rods were magnetized permanently by saturating the
NdPrFeB particles in the MRE using an applied uniaxial magnetic field (4.4 or
2.5 T) generated by an impulse magnetizer (IM-K-010020-A Magnet-Physik Dr.
Steingroever GmbH), as shown in Fig. 3 (e). The directions of ${\bm{B}}^{\rm
r}$ for either a straight or helical rods were set normal to the cross section
or parallel to the central axis of a helix, respectively. Schematic diagrams
for the magnetized straight and helical rods are shown in Figs. 3 (f) and (g),
respectively. The fabricated MRE rods possessed the residual flux densities of
$|{\bm{B}}^{\rm r}|=0.90c_{v}~{}{\rm T}$ (Secs. 6.1 and 6.3) and
$|{\bm{B}}^{\rm r}|=0.86c_{v}~{}{\rm T}$ (Sec. 6.2). The Young’s modulus was
$E=E_{\rm vps}/(1-c_{v}^{1/2})$ with the volume concentration of the particles
equal to $c_{v}=(1+(\rho_{\rm mag}/\rho_{vps})(m_{\rm tot}/m_{\rm
mag}-1))^{-1}$ [64]. After magnetization, the rods were demolded and cut to
the desired length $L$ (within an experimental uncertainty of $\pm 1$ mm).

###  5.2 Experimental apparatus

During the experimental tests, the samples are loaded magnetically by placing
them in between a set of two coaxial coils, which induce a steady axial
symmetric magnetic flux ${\bm{B}}^{\rm a}(r,z)$ [35]. Each coil is
manufactured by winding an aluminum circular spool with an insulated magnet
wire (Repelec Moteurs S.A.). The magnet wire (enameled wire, Isomet AG) has a
circular cross-section of diameter $1.32$ mm for the copper core and thickness
$0.102$ mm for the outer insulation layer. The final dimensions of the coil
are $86$ mm in inner diameter, $152$ mm in outer diameter, and $43$ mm ($33$
mm for wires and $10$ mm for the spool) in height; the mean radius of the coil
${\mathcal{R}}$ is ${\mathcal{R}}=59.5$ mm. The coils are powered by a DC
power supply providing a maximum current to power ratio of $25$ A/ 1.5kW (EA-
PSI 9200-25T, EA-Elektro-Automatik GmbH).

The two identical coils are set co-axially along $\hat{e}_{z}$ and their
current can flow in either the same or opposite directions. To realize a
constant (homogeneous) field, the center-to-center axial distance is set to be
$59.5$ mm$(={\mathcal{R}})$, and the current in the coils is set to flow in
the same direction; this is known as the Helmholtz coil configuration. In the
central region between the coils, the field generated by the Helmholtz coil is

$\displaystyle{\bm{B}}^{\rm a}=B^{\rm a}\hat{e}_{z}.$ (40)

In contrast, to realize the constant gradient field, the center-to-center
axial distance is set to be $103$ mm $(=\sqrt{3}{\mathcal{R}})$, and the
current in the coils is set to flow in opposite directions. This configuration
is known as the Maxwell coil, where the magnetic field induced near the coil
is

$\displaystyle{\bm{B}}^{\rm a}(r,z)=-\frac{b}{2}r\hat{e}_{r}+bz\hat{e}_{z}.$
(41)

In this Maxwell coil configuration, the gradient of the magnetic field is
$b=(\partial_{z}{\bm{B}}^{\rm a})\cdot\hat{e}_{z}$, i.e. constant. From the
Gauss law of magnetism, $\nabla\cdot\bm{B}^{\rm a}=0$, there are gradient
components in both the radial, $\bm{e}_{r}$, and the $\bm{e}_{z}$ directions.

We characterized the relationship between the flux density and the current in
the coils, $I_{c}$, by using a Teslameter (FH 55, Magnet-Physik Dr.
Steingroever GmbH) and the results are plotted in Fig. 4. The fields described
by Eqs. (40) and (41) are realized in the region of $z,r\lesssim
0.6{\mathcal{R}}$. For the Helmholtz configuration, the constant field is
induced near the center, for several values of $I_{c}$. In Fig. 4 (a), we plot
the profiles of the measured field $\bm{B}^{\rm a}$ divided by $I_{c}$, which
collapse onto a single curve. This collapse confirms that, $\bm{B}^{\rm a}$ is
proportional to $I_{a}$. By fitting the constant $B^{\rm a}$ near the center,
we find $\bm{B}^{\rm a}$ and $I_{a}$ as $B_{a}/I_{c}=8.9~{}{\rm mT/A}$. For
the Maxwell configuration, a constant gradient field is generated near the
center of the coils (Fig. 4 (b) and (c)). We measured the profiles of the
induced field both along with $\bm{e}_{z}$ and $\bm{e}_{r}$ for different
values of $I_{c}$. The data again collapse on a single curve as $\bm{B}^{\rm
a}/I_{c}$, which implies that the gradient of the field $b$ is linearly
proportional to $I_{c}$. By fitting the field along $\bm{e}_{z}$, we find
$b/I_{c}=0.133~{}{\rm T/(A\cdot m)}$. From Eq. (41), the gradient along
$\bm{e}_{r}$ is half of that along $\bm{e}_{z}$. In Fig. 4 (c), we plot the
slope of $b/2$ as the dashed line, which is in an excellent agreement with the
measurement. This experimental characterization confirms that our coils
accurately produce the fields according to Eqs. (40) and (41) when set in
either the Helmholtz or the Maxwell configurations, respectively.

Figure 4:  Experimental characterization of the field generated by the
magnetic coils in the (a) Helmholtz (constant field) and (b, c) Maxwell
(constant gradient field) configurations. (a) Profile of the (rescaled)
applied field along $x=y=0$. Profile of the applied field along (b) $x=y=0$
and (c) $y=z=0$. The dashed lines in the plot correspond to the (a) constant
and (b) (c) linear fits within the shaded regions where the experiments were
performed.

###  5.3 Experimental protocols

Having described the experimental apparatus above, we proceed by presenting
the protocols that we followed during the experiments. We placed the samples
(prepared through the procedure detailed in Sec. 5.1) between the coils (whose
field was characterized in the previous section). First, we will study the
deformation of a straight or curved rod under the constant field (Helmholtz
coil) followed by the study of the deformation of a helix under the constant
gradient field (Maxwell coil), as shown in Fig. 5 (a)-(d). A photograph of our
apparatus is shown in Fig. 5 (a). All the experimental results for the
corresponding setups will be presented in Sec. 6. The protocols for each
experiment are detailed next.

Figure 5:  (a) Photograph of the full experimental apparatus. A single power
supply generates the electric current provided to the two coils. A digital
camera set underneath the coils captures the deformation of the rod. (b)
Experimental setup used for the twist instability of a straight rod (Sec.
6.1). The straight rod magnetized along $-\hat{e}_{z}$ is clamped along
$\hat{e}_{x}$. The applied field ${\bm{B}}^{\rm a}$ is set as shown in Eq.
(40). (c) The schematic of the experimental setup for a helix under constant
field (right panel, See Sec. 6.2), with the uniformly helical magnetic rod
(left panel). The helix magnetized along $\hat{e}_{z}$ is clamped such that
its central axis is parallel to that of the coils. The magnetic helix is
immersed in a glycerol bath. (d) Photograph of a magnetic helical rod (in the
glycerol bath) under constant gradient field (Sec. 6.3).

In Sec. 6.1, we will present results on the deformation of a straight rod
under constant field. For these experiments, our MRE rods were magnetized
along $-\hat{\bm{e}}_{z}$, and a constant magnetic field was applied along
$\hat{\bm{e}}_{z}$; see Fig. 5 (b). An acrylic clamp, engraved to have the
convex shape at one of the ends of the rod, hanged the rod from above. This
clamp was also designed such that its width matched the gap between the coils
($13.5~{}{\rm mm}$). We varied the total length of the rod in the ranges
$L=40-57$ mm, the mass concentration ratio of NdPrFeB particles is
$c=10,~{}20,~{}30\%$, and the rod diameter was $d=3,5$ mm. We measured the
twist deformation of the rod (detailed in the next section) using a digital
camera, placed underneath the setup and directed upwards. The measured
twisting angles were averaged over 4 or 5 samples.

In Sec. 6.2, we will study the deformation of a helical rod under the constant
field. For these experiments, the helix was clamped vertically by a 3D-printed
rigid cylinder such that its central axis was parallel to the coil axis
($\hat{\bm{e}}_{z}$). The constant magnetic field (Eq. (40)) was applied
(anti-parallel to $\bm{B}^{\rm r}$) to the helical rod magnetized along
$\hat{\bm{e}}_{z}$, as shown in Fig. 5 (c). The experiments were performed in
a glycerol bath ($85$%, $1.23$ ${\rm g}/$cm3, Sigma-Aldrich), which density-
matched the rods, thereby minimizing the effect of gravity. We performed
experiments for the MRE rod of diameter $2$ mm, pitch angle $\psi=1.28$ rad,
and radius of curvature $R=10$ mm. The total length $L$ and the concentration
ratio $c$ were varied in the ranges $L=40-120$ mm and $c=10,20\%$,
respectively.

For the results in Sec. 6.3 and 7, we applied a constant gradient field
(according to Eq. (41)) to a helical rod, whose central axis in the reference
configuration was chosen to be aligned with the axis of the coils ($z$-axis,
perpendicular to gravity). As shown in Fig. 5 (d), we clamped the helix and
perform the experiment in a glycerol bath. During the experiments, we imaged
the deformed shape of the helix using a digital camera. In Sec. 6.3, we
measured the displacement of the free-end as a function of the applied field.
The total length $L$ and the concentration ratio $c$ were varied as
$L=74,81,103$ mm and $c=10,20\%$, respectively, while the diameter $d=2$ mm,
the radius of natural curvature $R=10$ mm and the pitch angle $\psi=1.51$ rad
were fixed. In Sec. 7, the deformation of the magnetic helix were studied,
varying the pitch angle $\psi$ and total length $L$ in the range of
$\psi=1.26-1.49$ rad and $L=65-140$ mm, while $d=2$ mm, $R=10$ mm, and
$c=20\%$ were kept fixed throughout Sec. 7.

##  6 Validation of the theory of hard magnetic rods vs experiments

In this section, we will use precision experiments on specific test
configurations to perform a detailed validation of the theoretical framework
introduced in Sec. 3.3 for hard MRE rods. In Sec. 6.1 and 6.2, the twist
instabilities for a straight rod and a helix are investigated under a constant
field. Given that $\bm{p}_{\rm mag}=0$ under the constant field, we can show
that the elasto-magnetic torque $\bm{q}_{\rm mag}$ in both naturally straight
and curved configurations captures the corresponding experimental results in
excellent agreement. Then, in Sec. 6.3, we study the deformation of a helix
under a constant gradient field to test the validity of both $\bm{p}_{\rm
mag}$ and $\bm{q}_{\rm mag}$.

###  6.1 Twist instability of a straight magnetic rod under magnetic loading

Consider a straight rod aligned with $\hat{\bm{e}}_{x}$ (clamped at $s=0$),
that is magnetized normal to the tangent, i.e., ${\bm{B}}^{\rm
r}=B^{r}\hat{\bm{e}}_{z}$. We apply the constant flux density as
${\bm{B}}^{\rm a}=-B^{a}\hat{\bm{e}}_{z}$ (see Fig. 6 (a)-(d)). Figures 6 (a)
and (c) are photographs from the experiments of the rod cross-sections at
their free-ends, parallel with $y$-$z$ planes. Three-dimensional
representations of the straight MRE rods obtained from the simulations are
shown in Figs. 6 (b) and (d). The rod remains straight when the magnitude of
the applied field $B^{\rm a}$ is sufficiently small enough (Fig. 6 (a) and
(b)), whereas the rod twists, with the centerline still straight, above the
critical flux density $B^{a*}_{\rm twist}$ (Fig. 6 (c) and (d)). In the
experiments, we measured the twist angle $\phi_{L}$ at $s=L$ by tracking the
rotation angle of a nitinol rod embedded near the free-end perpendicular to
the centerline (as shown by $\hat{n}(L)$ in Figs. 6 (a) and (c)). In the
simulation, we also measured $\phi_{L}$ by computing the rotation angle for
one of the basis vectors of the Cosserat frame $\hat{\bm{d}}_{1}$. In Fig. 6
(e), we superpose the experimental and numerical results of $\phi_{L}$ as a
function of $B^{\rm a}/B^{\rm a*}_{\rm twist}$. The error bars in the
experimental data correspond to the uncertainty in the setting of the angle of
the clamp. Experimental and numerical results are in excellent agreement with
each other across the full parameter range in $B^{\rm a}/B^{\rm a*}_{\rm
twist}$.

Figure 6: Twist instability of straight MRE rods. (a) Experimental photograph
of the cross section at $s=L$ with $B^{\rm a}<B^{\rm a*}_{\rm twist}$. The
symbols $\hat{n}(L)$ and $\mathring{\hat{n}}(L)$ represent the direction of
the nitinol rod in the deformed and reference configuration, respectively. The
corresponding simulation snapshot is shown in (b). (c) Experimental photograph
of the cross section at $s=L$ with $B^{\rm a}>B^{\rm a*}_{\rm twist}$. The
corresponding simulation snapshot is shown in (d). (e) $\phi_{L}$ as a
function of $B^{\rm a}/B^{\rm a*}_{\rm twist}$. The data points are
experimental results and the solid line is the simulation result. (f)
Logarithmic plot of $\phi_{L}$ as a function of $B^{\rm a}/B^{\rm a*}_{\rm
twist}-1$. The legend is the same as (e). The dashed line is the fitting of
the power law with the exponent $1/2$, Eq. (48), obtained through scaling
arguments.

We rationalize the experimental and numerical results using the magnetic
Kirchhoff equations that we derived in Sec. 3.3 by the perturbation against
the twist instability. Following Ref. [65], we expand the Cosserat frame with
respect to the straight configuration
$(\hat{\bm{d}}_{1}^{(0)},\hat{\bm{d}}_{2}^{(0)},\hat{\bm{d}}_{3}^{(0)})=(\hat{\bm{e}}_{y},\hat{\bm{e}}_{z},\hat{\bm{e}}_{x})$
as

$\displaystyle\hat{\bm{d}}_{a}=(\delta_{ab}+\varepsilon_{abc}\alpha_{c})\hat{\bm{d}}_{b}^{(0)},$
(42)

where $\varepsilon_{abc}\alpha_{c}$ is the asymmetric tensor and
$\alpha_{a}=\alpha_{a}(s)~{}(a=1,2,3)$ is the small (angular) perturbation
parameter $\alpha_{a}\ll 1$. The coefficient $\varepsilon_{abc}\alpha_{c}$
enforces the orthogonality of $\hat{\bm{d}}_{a}~{}(a=1,2,3)$ even up to the
1st order of the perturbation. Plugging the perturbed form of
$\hat{\bm{d}}_{a}$, Eq. (42), into the kinematic equation, Eq. (1), (see A for
details), we find

$\displaystyle\alpha_{a}^{\prime}=\Omega_{a}\,,$ (43)

from which we conclude that $\alpha_{a}$ represents the infinitesimal rotation
angle around $\hat{\bm{d}}_{a}$. Given that the rod is clamped at $s=0$ and
free at $s=L$, the corresponding boundary conditions are $\alpha_{a}(0)=0$ and
$\alpha_{a}^{\prime}(L)=\Omega(L)=0$. Furthermore, ${\bm{p}}_{\rm mag}=0$ and
thus ${\bm{F}}(s)=0$. Substituting Eq. (42) into Eq. (21) and linearizing the
result for $\alpha_{a}$ yields

$\displaystyle EI\alpha_{1}^{\prime\prime}=-\frac{AB^{\rm r}B^{\rm
a}}{\mu_{0}}\alpha_{1},$ (44) $\displaystyle EI\alpha_{2}^{\prime\prime}=0,$
(45) $\displaystyle GJ\alpha_{3}^{\prime\prime}=-\frac{AB^{\rm r}B^{\rm
a}}{\mu_{0}}\alpha_{3}.$ (46)

Equations (44)-(46) are the governing equations for the infinitesimal rotation
angles $\alpha_{a}$ along the rod. To determine the conditions for the
existence of nontrivial solutions $(\alpha_{a}(s)\neq 0)$ for Eqs. (44)-(46),
we will substitute $\alpha_{a}\propto\sin(\pi s/2L)~{}(a=1,2,3)$, which
satisfies the boundary conditions $\alpha_{a}(0)=\alpha_{a}^{\prime}(L)=0$,
into Eqs. (44)-(46). Equation (45) gives us the null amplitude
$\alpha_{2}(s)=0$, while Equations (44) and (46) are satisfied if ${AB^{\rm
r}B^{\rm a}}/{\mu_{0}}=EI(\pi/2L)^{2}$ or ${AB^{\rm r}B^{\rm
a}}/{\mu_{0}}=GJ(\pi/2L)^{2}$ hold, respectively. The smaller value of $B^{\rm
a}$ derived here corresponds to the critical value of $B^{\rm a}$ for the
twist instability, $B^{\rm a*}_{\rm twist}$. Since $GJ/EI=1/(1+\nu)<1$, the
eigenvalue of Equation (46) corresponds to the critical flux density, above
which the rod twists: Hence, we find

$\displaystyle B^{\rm a*}_{\rm
twist}=\left(\frac{\pi}{2}\right)^{2}\frac{\mu_{0}GJ}{AB^{\rm r}L^{2}}.$ (47)

Notice that $EI$ in Eq. (24) is now replaced by $GJ$ in Eq. (47). Given that
$\alpha_{3}=\alpha_{3}(s)$ is the profile of the (infinitesimal) twisting
angle (from Eq. (43)), $\alpha_{3}(L)$ corresponds to the twist angle of the
free-end measured in our experiment and $\phi_{L}=\alpha_{3}(L)$ in the limit
of $\phi_{L}\ll 1$. By analogy with the equation of clamped-free elastica [43]
(replacing $\alpha_{3}(s)$ with the bending angle in elastica), Equation (46)
is now mathematically equivalent to the amplitude equation for the pitch-fork
bifurcation. In the neighborhood of $B^{\rm a}/B^{\rm a*}_{\rm twist}\sim 1$,
we thus anticipate that $\phi_{L}$ obeys the following scaling law, where the
twist angle at the tip $\phi_{L}$ evolves as

$\displaystyle\phi_{L}\propto\sqrt{\frac{B^{\rm a}}{B^{\rm a*}_{\rm
twist}}-1}.$ (48)

We find an excellent agreement between experiments and the theoretical
description developed above $\phi_{L}$, although the transition is not sharp
due to the inevitable imperfection of the system below the critical field
$B^{\rm a}/B^{\rm a*}_{\rm twist}<1$ [66]. Furthermore, from the logarithmic
plot in Fig. 6 (f), we confirm that $|\phi_{L}|$ is consistent with the
predicted scaling law Eq. (48), in the region of $|B^{\rm a}/B^{\rm a*}_{\rm
twist}-1|\ll 1$. The above results demonstrate that our magnetic Kirchhoff
equations correctly predict the twist instability of straight rods. In
particular, we validated the elasto-magnetic torque $\bm{q}_{\rm mag}$ in the
case of straight rods. Combining the fact that our framework reproduces the
previous results on the planar deformation (Sec. 3.4) [28, 32, 31, 35] and the
results presented in this subsection, the magnetic Kirchhoff equation (20) and
(21) correctly predict the large deformation of naturally straight hard MRE
rods. Next, we focus on the deformation of hard MRE rods that are naturally
curved .

###  6.2 Buckling of a magnetic helix under constant field

We proceed by considering a naturally curved and twisted rod, i.e., a helix,
which undergoes structural (buckling) instability of a helix at a critical
field $B^{\rm a*}$, but with an inevitable twist-bend coupling. We will show
that our simulations and theoretical analyses rationalize our experimental
results with excellent quantitative agreement.

We fabricate a hard magnetic helix with the pitch angle $\psi=1.28~{}{\rm
rad}$, which is clamped such that the central axis is aligned along the coil
axis. We apply the magnetic field ${\bm{B}}^{\rm a}=-B^{\rm
a}\hat{\bm{e}}_{z}$ on a helical rod with lengths in the range of $40\leq
L~{}[{\rm mm}]\leq 120$. These helices are magnetized along their central axis
as ${\bm{B}}^{\rm r}=B^{\rm r}\hat{\bm{e}}_{z}$. Below the critical field
$B^{\rm a}<B^{\rm a*}$, the rod remains uniformly helical, as shown in the
filled data points in Fig. 7 (a), buckling occurs when $B^{\rm a}>B^{\rm a*}$,
as depicted by the open symbols. We regarded the rod as helical if the tip
moves by less than the diameter $d$ in the experiments.

We now specialize the magnetic Kirchhoff equations Eqs. (20) and (21) for a
magnetic helix to rationalize the experimental observations. The critical flux
density for buckling of the magnetic helix $B^{\rm a*}$ can be derived via
rigorous analysis of Eqs. (20) and (21). Here, to simplify the analysis and
emphasize the importance of a twist-bend coupling in the system, we will
compute the critical field strength $B^{a*}$ in the limit of small curvature
$KL\ll 1$ ($K$ is total curvature for a helix defined later) but with
$\psi<\pi/2$, which indeed provides a good prediction for our observations as
we detail later.

Figure 7: Buckling of a magnetic helix. (a) Phase diagram of the buckling
instability of a magnetic helix in experiments. Filled and empty data points
are uniform or buckled helix states, respectively. Buckling of a magnetic
helix under constant field. In the previous section. The dashed line is the
analytic prediction in the limit of $L/R\ll 1$. The simulation snapshots for
$L/R=12$ correspond to $B^{a}/B^{a*}=0.25$ (b) and $B^{a}/B^{a*}=2.5$ (c). (d)
Rescaled phase diagram of our experimental results. The simulation results and
asymptotic solutions are superposed with $\psi=1.28~{}{\rm rad}$.

Similarly to Sec. 6.1, we expand the Cosserat frame with respect to the
uniform helix as we did in Eqs. (42) and (43). The Darboux vector in the
uniform helix is
${\bm{\Omega}}=\kappa{\hat{\bm{d}}}_{2}+\tau{\hat{\bm{d}}}_{3}$ with the
natural curvature $\kappa$ and twist $\tau$. By integrating the kinematic
relation Eq. (1),
$\hat{\bm{d}}_{a}^{\prime}={\bm{\Omega}}\times\hat{\bm{d}}_{a}$, the Cosserat
frame $\hat{\bm{d}}_{a}~{}(a=1,2,3)$ for a helix in the reference
configuration is computed as

$\displaystyle\hat{\bm{d}}_{1}^{(0)}(s)$ $\displaystyle=$
$\displaystyle\cos(Ks)\hat{\bm{e}}_{x}-\sin(Ks)\hat{\bm{e}}_{y},$ (49)
$\displaystyle\hat{\bm{d}}_{2}^{(0)}(s)$ $\displaystyle=$
$\displaystyle-\cos\psi\left(\sin(Ks)\hat{\bm{e}}_{x}+\cos(Ks)\hat{\bm{e}}_{y}\right)-\sin\psi\hat{\bm{e}}_{z},$
(50) $\displaystyle\hat{\bm{d}}_{3}^{(0)}(s)$ $\displaystyle=$
$\displaystyle\sin\psi\left(\sin(Ks)\hat{\bm{e}}_{x}+\cos(Ks)\hat{\bm{e}}_{y}\right)-\cos\psi\hat{\bm{e}}_{z}.$
(51)

We choose the central axis to be parallel to $\hat{z}$ and define the pitch
angle $\psi$ as $\kappa=K\sin\psi$, $\tau=\kappa\cos\psi$, and
$K\equiv\sqrt{\kappa^{2}+\tau^{2}}$. The radius of curvature on the $x-y$
plane $R$ and $\psi$ are related as $\kappa=\sin^{2}\psi/R$ and $\tau=\sin
2\psi/2R$ (see B for the detailed derivation).

Plugging Eqs. (49)-(51) into the magnetic Kirchhoff equations (Eqs. (20) and
(21)) and using Eqs. (42) and (43), we obtain the set of equilibrium equations
for the linear perturbation of a helix in the limit of $KL\ll 1$:

$\displaystyle EI\alpha_{1}^{\prime\prime}$ $\displaystyle=$
$\displaystyle-\frac{AB^{\rm r}B^{\rm a}}{\mu_{0}}\alpha_{1}$ (52)
$\displaystyle EI\alpha_{2}^{\prime\prime}$ $\displaystyle=$
$\displaystyle\frac{AB^{\rm r}B^{\rm
a}}{\mu_{0}}\cos\psi\left(\alpha_{2}\cos\psi-\alpha_{3}\sin\psi\right)$ (53)
$\displaystyle GJ\alpha_{3}^{\prime\prime}$ $\displaystyle=$
$\displaystyle-\frac{AB^{\rm r}B^{\rm
a}}{\mu_{0}}\sin\psi\left(\alpha_{2}\cos\psi-\alpha_{3}\sin\psi\right),$ (54)

with boundary conditions that are the same as those given in Sec. 6.1:
$\alpha_{a}(0)=\alpha_{a}^{\prime}(L)=0$. Equations (52)-(54) are the
linearized moment balance equations of Eq. (21).

Assuming the functional form $\alpha_{a}\propto\sin(ks)$ with $k=\pi/2L$ (that
satisfies the boundary conditions $\alpha_{a}(0)=\alpha_{a}^{\prime}(L)=0$),
the critical magnetic field $B^{a*}$ satisfying $\alpha_{a}\neq 0$ can be
obtained from Eqs. (53) and (54) as

$\displaystyle B^{a*}=\left(\frac{\pi}{2L}\right)^{2}\frac{\mu_{0}}{AB^{\rm
r}}\frac{EI}{\cos^{2}\psi+(EI/GJ)\sin^{2}\psi}.$ (55)

Equation (55) captures the previous results on planar bending in Eq. (24) in
Sec. 3.4, as well as the twist instabilities in Eq. (47), which we
investigated in Sec. 6.1. Indeed, by taking the limit $\psi\to 0$, (i.e., rod
clamped vertically and magnetized along its tangent), we get $B^{a*}\to
B^{a*}_{\rm bend}$. By contrast, in the limit of $\psi\to\pi/2$, (i.e., rod
clamped horizontally but magnetized as in Sec. 6.1), we recover $B^{a*}\to
B^{a*}_{\rm twist}$.

We performed numerical simulations corresponding to the buckling of a hard
magnetic helix under the constant magnetic field. When the applied field is
sufficiently small $B^{\rm a}<B^{a*}$, the rod remains helical, but the helix
buckles when $B^{\rm a}>B^{a*}$. In Figs. 7 (b) and (c), we show
representative snapshots of the numerical simulations for $B^{\rm a}<B^{a*}$
and $B^{\rm a}>B^{a*}$, respectively. Based on Eq. (55), we summarize
experimental and numerical results in the rescaled phase diagram shown in Fig.
7 (d). In the simulations, we remark that the rod remains helical if the
displacement of the free-end is less than $10^{-3}\ell_{0}$, where $\ell_{0}$
is the natural length of the spring in the simulation (unit length of the
simulation). The phase boundary between the undeformed (helical) and buckled
states is in an excellent agreement between experiments and simulations.
Furthermore, the analytic prediction (Eq. (55)) correctly predicts both
experimental and numerical results. It should be noted that not only Eq. (55)
predicts the phase boundary for $L/R\ll 1$, but also that for $L/R=O(1)$,
without any adjustable parameters. A more rigorous analysis will be able to
derive the full analytic prediction of the phase boundary, wheres it is beyond
the scope of the current paper.

In this subsection, we have confirmed that the elasto-magnetic torque
$\bm{q}_{\rm mag}$ Eq. (19) derived in Sec. 3.1 correctly predicts the
experimental results for a naturally curved and twisted rod. Combining the
results from Sec. 6.1 and those in this subsection, we validated the magnetic
Kirchhoff equations (Eqs. (20) and (21)) in 3D deformation under the constant
field. In the next subsection, we consider the case of loading in a constant
gradient magnetic field to validate the elasto-magnetic force $\bm{p}_{\rm
mag}$ (Eq. (18)) to fully validate our theoretical framework.

###  6.3 Deformation of a magnetic helix under constant gradient field

In this subsection, we will validate the framework for magnetic helices under
the constant gradient field (${\bm{p}}_{\rm mag},{\bm{q}}_{\rm mag}\neq 0$).
For this third and final test configuration, the magnetic helix is clamped at
$s=0$, while setting the other end ($s=L$) to be free. As we increase the
gradient of the field $b$ (or equivalently, the rescaled gradient
$\lambda_{m}=L/\ell_{m}$, see Eq. (28)), the magnetic helix stretches. We will
present results on the displacement of the free-end as a function of
$\lambda_{m}$ and find excellent agreement between experiments and
simulations. In this subsection, we study the magnetic helix of a pitch
shorter than that in Sec. 6.2. Hence, the helix interacts with itself non-
locally via long-range interaction. Although long-range interactions are not
considered in our theory, we can predict the stretching of the helix as long
as the pitch becomes larger upon increasing $\lambda_{m}$.

As an illustrative configuration, the magnetic helix is placed with its
central axis is aligned with the coil axis, and the clamped-end is located on
the $z=0$ plane (see Fig. 8 for a photograph of the setup), within the error
of $\pm 1$ mm. The magnetic helix has geometries of $(L,R)=(10,103)$ ${\rm
mm}$ and pitch angle $\psi=1.51~{}{\rm rad}$. We suspend the magnetic helix in
a glycerol bath, horizontally. The gradient of the magnetic field $b$ is
increased from zero up to $b=1.06~{}{\rm T/m}$ (maximum value in our setup),
and then $b$ is decreased back to zero.

In the absence of the applied field, $\lambda_{m}=0$ (i.e. $b=0$), the helix
is fully contracted due to the long-range dipole-dipole attractive
interaction. As discussed in Refs. [30, 35] and Sec. 3, the applied field
gradient is analogous to a unidirectional body force of constant magnitude
such as gravity. Hence, the helix extends (or contracts) when $\lambda_{m}$ is
increased or decreased, respectively. We present the photograph of the
extended helix with $\lambda_{m}=3.0$ in Fig. 8 (a). The rescaled gradient is
then decreased to $\lambda_{m}=2.5$ and $\lambda=2.0$, as shown in Fig. 8 (b)
and (c), respectively. We performed the numerical simulations of the magnetic
helix under the constant gradient field to compare its deformed configurations
with the experimental ones. Figures 8 (d), (e) and (f) are the snapshots of
our discrete simulation with $\lambda_{m}=3.0,2.5$ and $2.0$, respectively.
Similar to our experimental observations, the helix contracts as $\lambda_{m}$
is decreased. The predicted shapes in the numerical simulation agree with the
configurations in experiments qualitatively.

Figure 8: Deformation of a hard magnetic helix. Snapshots of (a)-(c)
experiments and (d)-(f) simulations with ${B}^{\rm r}=13~{}{\rm mT}$,
$R=10~{}{\rm mm}$ and $L=103~{}{\rm mm}$. The corresponding values of rescaled
gradient are (a) (d) $\lambda_{m}=3.0$, (b) (e) $\lambda_{m}=2.5$, (c) (f)
$\lambda_{m}=2.0$.

To quantify the deformation of the helix in the present configuration, we
measured the end-to-end distance along $z$ between the clamp and the tip of
the helix: $z_{L}\equiv z(L)-z(0)$ defined in Fig. 8 (a). In Fig. 9 (a)-(c),
we plot $z_{L}$ as a function of the rescaled applied field gradient
$\lambda_{m}=L/\ell_{m}$ for different values of $L$ or $B^{\rm r}$, where
$\ell_{m}$ is the magneto-bending length defined in Eq. (28). The data points
correspond to the experimental results. As an example, we describe the
experimental results presented in Fig. 9 (a). In the absence of the applied
field $\lambda_{m}=0$, the helix is in self-contact due to long-range
attractive interactions. Upon increasing $\lambda_{m}$, the helix loses self-
contact when $\lambda_{m}=\lambda_{m}^{+}$. Past this point, $z_{L}$ increases
as the elasto-magnetic force stretches the helix. When $\lambda_{m}$ is
decreased, the helix contracts until self-contact at
$\lambda_{m}=\lambda_{m}^{-}$. For the sake of convenience, we refer to the
branch where the helix is in self-contact self-contact branch, and to the
branch without self-contact as the elasto-magnetic branch.

We find that the values of $\lambda_{m}^{\pm}$ depend on $B^{\rm r}$ and $L$.
Indeed, in Figs. 9 (a) and (c), $\lambda_{m}^{-}\neq 0$, while, in Fig. 9 (b),
$\lambda_{m}^{-}=0$. The fact that $\lambda_{m}^{-}=0$ implies that the
helical shape (without contact) is stable but it is not a unique stable state
due to the long-range interactions. The value of $\lambda_{m}^{\pm}$ can be
used to classify the hysteric behavior, depending on the magnitude of residual
flux density $|B^{\rm r}|$, the pitch angle $\psi$, and the rescaled total
length $L/R$. We will discuss the hysteric behavior in $z_{L}$ more in detail,
in Sec. 7.

Complementing the experimental data, we performed the discrete simulations
corresponding to our experimental settings to measure the end-to-end distance
$z_{L}$ as a function of $\lambda_{m}$. Given that our theory does not include
the self long-range interaction, $z_{L}$ is not hysteretic (the self-contact
branch does not appear). Still the predictions from our model are in excellent
agreement with the experimental results in the elasto-magnetic branches. In
Figs. 9 (a)-(c), we plot the predictions from the simulations as solid lines.
The width of the shaded region in Figs. 9 (a)-(c) correspond to the
propagation of errors due to the $\pm 1$ mm, uncertainty of the total length
$L$ measurement. From the excellent agreement between the experimental data
and numerical results, we conclude that the magnetic Kirchhoff equations (Eqs.
(20) and (21)) provide an accurate quantitative prediction for the deformation
of the magnetic helix under a constant gradient field, when the self long-
range interaction is not dominant.

Figure 9: The $z$-position of the free-end as a function of the rescaled
gradient $\lambda_{m}=L/\ell_{m}$ for different $B^{\rm r}$ and length $L$.
The data points and solid lines are experimental and simulation results,
respectively. The shaded areas are errors due to the clamped position computed
from the simulation. (a) $B^{\rm r}=13~{}{\rm mT}$, $L=103~{}{\rm mm}$ (b)
$B^{\rm r}=13~{}{\rm mT}$, $L=74~{}{\rm mm}$ (c) $B^{\rm r}=26~{}{\rm mT}$,
$L=81~{}{\rm mm}$.

In this section, we combined the discrete simulations with experiments, to
validate the magnetic Kirchhoff equations (Eqs. (20) and (21)) derived in Sec.
3.3, for specific cases where the long-range interaction can be neglected
compared with the elasto-magnetic forces. We showed that the numerical
simulations can predict the behavior of the elasto-magnetic branch observed
experimentally. We have validated the elasto-magnetic torque $\bm{q}_{\rm
mag}$ and force $\bm{p}_{\rm mag}$ in a naturally straight and curved rod,
thereby allowing us to validate Eqs. (20) and (21) in 3D geometrically
nonlinear deformations, at least when the long-range interaction is
negligible. Next, in Sec. 7, we perform a more systematic experimental study
of when and how the long-range interaction affects the deformation of a
magnetic helix.

##  7 Limitations of the theory

In this section, we systematically quantify the limitation of the theory
through precision experiments alone, because the self long-range interactions
is not included in Eqs. (20) and (21); doing so is a challenging endeavour
that goes beyond the scope of the present study. We classify the hysteretic
behavior reported in Fig. 9 based on the value of $\lambda_{m}^{\pm}$; the
critical field strength at which $z_{L}$ jumps discontinuously. We will
discuss the classification by considering the structural instability of the
magnetic helix in the absence of the applied field, that is $\lambda_{m}=0$.

The long-range interactions in the magnetic helix depend on the combined
effects of the magnitude of the magnetization
$|\bm{\mathcal{M}}|\propto|\bm{B}^{\rm r}|$ and the geometry of the helix
(e.g., the pitch angle $\psi$ and the rescaled total length $L/(2\pi R)$). The
magnitude of the magnetization $|\bm{\mathcal{M}}|$ is programmed at the
fabrication stage by the concentration ratio of the NdPrFeB particles $c$.
When $c$ is small, the long-range interaction will be suppressed (and vice
versa). On the other hand, the dependence of $\psi$ and $L/(2\pi R)$ is highly
nontrivial, as we discuss below.

The constant gradient field is applied to the magnetic helix as in Sec. 6.3.
The helix is extended and contracted under the constant gradient magnetic
field, for several pitch angles (in the natural configuration). We introduce
the normalized pitch angle as $\tilde{\psi}\equiv\psi/(\pi/2)$ such that
$\tilde{\psi}=1$ corresponds to the planar circle with zero twist
(experimentally non-realizable due to the finite diameter of the rod) and
$\tilde{\psi}=0$ corresponds to a straight rod. As $\tilde{\psi}$ decreases,
the pitch of the helix increases. Practically, $\tilde{\psi}_{\rm max}=0.96$
is the largest value of $\tilde{\psi}$ realizable in our experimental
conditions with $(d,R)=(2,10)$ ${\rm mm}$. Our rods are fabricated
systematically through the procedure in Sec. 5.1 with different pitch angles
in the range $0.80\leq\tilde{\psi}\leq 0.95$. We clamp the magnetized helix
with $B^{\rm r}=26~{}{\rm mT}$ at $z(0)=0$, while $s=L$ is set to be free.
Through the cyclic protocols consisting of increasing and decreasing the flux
density, the values of $\lambda_{m}$ characterizing the hysteresis
$\lambda_{m}^{\pm}$ are identified (see Fig. 9 for the definition of
$\lambda_{m}^{\pm}$).

Figure 10: Classification of the hysteretic behavior in the stretching and
contraction of a magnetic helix. The phase diagram of hysteresis is obtained
from experiments with $B^{\rm r}=26{\rm mT}$ and $R=10{\rm mm}$; $\bigcirc$:
No hysteresis ($\lambda_{m}^{+},\lambda_{m}^{-}=0$), $\triangle$: Meta-stable
helix ($\lambda_{m}^{+}\neq 0,~{}\lambda_{m}^{-}=0$) and $\times$: Fully
contracted helix ($\lambda_{m}^{+},\lambda_{m}^{-}\neq 0$).

In Fig. 10, we present the experimental phase diagram in the
$\tilde{\psi}$-$L/2\pi R$ parameter space to classify the response against
$\lambda_{m}$, finding that there are three different possible responses
(i)-(iii), depending on the geometry of the magnetic helix. (i) When $L/2\pi
R\simeq 1$ (length of the helix is nearly a single turn) or $\tilde{\psi}$ is
small, (i.e. when the pitch $p=2\pi R/\tan\psi$ is long enough), no hysteresis
is observed (circles in Fig. 10). When we plot $z_{L}$ as a function of
$\lambda_{m}$, $z_{L}$ changes continuously and reversibly for any
$\lambda_{m}$. The tip of the helix follows the same curve upon increasing or
decreasing $\lambda_{m}$. In this case, we have
$\lambda_{m}^{+}=\lambda_{m}^{-}=0$ (Fig. 10). For higher values of $\psi$, we
observe hysteresis, which is divided into two further categories;
$\lambda_{m}^{+}\neq 0,~{}\lambda_{m}^{-}=0$ or
$\lambda_{m}^{+},\lambda_{m}^{-}\neq 0$, which are shown as triangles and
cross marks in Fig. 10, respectively. (ii) In the case of $\lambda_{m}^{+}\neq
0,~{}\lambda_{m}^{-}=0$, $z_{L}$ changes discontinuously from the contracted
helix to the stretched helix at $\lambda_{m}=\lambda_{m}^{+}$, while $z_{L}$
decreases smoothly as $\lambda_{m}\to 0$, i.e. $\lambda_{m}^{-}=0$ (see Fig. 9
(b) for the corresponding behavior). (iii) By contrast, in the case of
$\lambda_{m}^{+},\lambda_{m}^{-}\neq 0$, $z_{L}$ jumps at both
$\lambda_{m}^{+}$ and $\lambda_{m}^{-}$. Thus, if we plot $z_{L}$ as a
function of $\lambda_{m}$, we obtain similar curves as those in Figs. 9 (a)
and (c).

The classification of the hysteretic responses in the $z_{L}$-$\lambda_{m}$
plots (i)-(iii) can be understood by considering the configuration at
$\lambda_{m}=0$. In the following, we discuss the stability of the uniform
helix at $\lambda_{m}=0$ against mechanical perturbations. (i) When
$\lambda_{m}^{+}=\lambda_{m}^{-}=0$ (no hysteresis), the uniform helical shape
corresponds to a unique energy minimum. The reference configuration is mono-
stable upon any mechanical perturbations, because self-interactions have
minimal effects. (ii) For $\lambda_{m}^{+}\neq 0$ and $\lambda_{m}^{-}=0$, the
uniform helix is stable but it is not a unique minimum. In this case, the
helix can contracts to the (stable) state with the normalized pitch angle of
$\tilde{\psi}=\tilde{\psi}_{\rm max}$ if perturbed; we call this state meta-
stable helix. Lastly, (iii) for $\lambda_{m}^{+},\lambda_{m}^{-}\neq 0$, the
uniform helix is no longer stable. Indeed, the contracted state is the unique
energy minimum. In this third case, the helix remains fully-contracted against
any mechanical perturbation.

The nature of the hysteretic behavior in magnetic helices depends on the
concentration ratio of NdPrFeB particles $c$ and the geometry of the helix as
described by $\psi$ and $L/(2\pi R)$. We classified the deformation due to
long-range self-interactions based on the value of $\lambda_{m}$ at which
$z_{L}(\lambda_{m})$ changes discontinuously. The hysteretic behavior in the
magnetic helix originates from the bi-stability due to the long-range
interactions. In the previous sections, we had validated the magnetic
Kirchhoff rod equations when the long-range interaction can be neglected. At
the same time, we clarified the limitation of our theory systematically. The
hysteretic behavior, which our theory cannot capture, is observed only when
the pitch angle $\tilde{\psi}$ is large and when the rod is long $L/(2\pi
R)\gtrsim 1$. We can conclude that the magnetic Kirchhoff rod equations (Eqs.
(20) and (21)) are valid as long as $\psi$, $L/(2\pi R)$, or $c$ remains
small.

##  8 Conclusion

The goal of the present study was to develop a Kirchhoff-like theory for hard
magnetic rods based on dimensional reduction. Our theoretical framework,
rooted in 3D elasticity, yielded a reduced (centerline-based) description of
the magnetic rod (Sec. 3). The set of governing equations contain elasto-
magnetic forces and torques (Eqs. (18) and (19)), in addition to the purely
elastic ones for a (non-magnetic) rod. We validated the magnetic Kirchhoff rod
equations (Eqs. (20) and (21)) through a set of precision experiments. Our
theoretical results supplemented by the discrete simulation are in excellent
agreement with the experimental results. Moreover, the magnetic Kirchhoff rod
equations reproduce the previous results on the planar deformation of hard
magnetic beams or elastica [28, 30, 35, 32]. We confirmed that the force and
moment balance equations in our theory for hard magnetic rods reduce to those
of Refs. [28, 30, 35, 32] for beams and elastica.

To validate our theory for the 3D deformation of magnetic rods, we performed
three different sets of experiments; (i) a straight magnetic rod under
constant external magnetic field (Sec. 6.1), (ii) a helical rod under constant
field (Sec. 6.2), and (iii) a helical rod under the constant gradient field
(Sec. 6.3), with the clamped-free boundary conditions. (i) In the case of a
straight rod under the constant field, we studied the twist instability. When
the rod is magnetized along $\hat{\bm{d}}_{1}$ (perpendicular to the tangent
$\hat{\bm{d}}_{3}$) and we apply an anti-parallel field, the rod twists, while
its centerline remains straight. We derived an analytical prediction for the
critical applied field for this twist instability $B^{\rm a*}_{\rm twist}$,
which is in excellent agreement with experimental results. (ii) We studied the
buckling instability of a helical rod under constant external magnetic field.
The helix was magnetized along its central axis. When we apply the field to
the magnetic helix, it buckles above the critical applied field $B^{\rm a*}$.
We showed that this instability is triggered by a balance between elastic and
elasto-magnetic torques. The critical applied fields were predicted
analytically from a linear stability analysis. The simulation results
correctly capture the experimental observations. (iii) Under a constant
gradient field, the helix magnetized along the central axis stretches (or
contracts) as we increase (or decrease) the applied field strength. Due to the
long-range self-interaction between magnetic dipoles within the rod, the
deformation is hysteretic upon extension and contraction.

Our framework correctly predicts the deformation of the magnetic helix when
its pitch is large enough or when the long-range self-interactions are
negligible. The hysteretic behavior was studied systematically through
experiments by controlling the pitch angle $\psi$ of the helix. We revealed
that there are three types of behavior; no hysteresis state, meta-stable, and
fully contracted helices, depending on the total length $L/R$ and the pitch
angle $\psi$.

Although our theory is in excellent agreement with the experimental results
when the elasto-magnetic force and torque are dominant, the non-local
interactions between the magnetization vectors $\bm{\mathcal{M}}$ are
necessary to describe the hysteretic behavior. In Refs. [15, 16, 17], the
continuum mechanics of magnetic chains has been studied by considering the
self-interactions between magnetic beads. In the future, combining the
continuum mechanics of the magnetic dipole interactions and our magnetic
Kirchhoff rod equations, the hysteretic behavior of the magnetic helix could
be studied in detail theoretically. Extending the magnetic Kirchhoff equations
to include the dipole-dipole interactions is one of the interesting directions
of future works. Also, extending our theoretical framework toward dynamics of
hard magnetic rods is also another exciting opportunity for future work. The
dynamics could be simulated by considering the (external) viscous force and
torque of a surrounding medium in Eqs. (20) and (21).

The framework established in this paper would be valuable to simulate the
large deformation of hard magnetic rods used in micro or soft-robotics [36,
39, 5, 32] or haptic devices [40]. In the future, novel functional devices of
complex geometries could be developed by combining the simple building blocks
studied in this paper.

## Acknowledgement

T.G.S. was supported by Grants-in-Aid for Japan Society for the Promotion of
Science Overseas Research Fellowship 2019-60059 and MEXT KAKENHI 18K13519.

##  Appendix A Derivation of Equation (43)

In this appendix, we provide the detailed derivation of Eq. (43);

$\displaystyle\alpha_{a}^{\prime}=\Omega_{a}\,.$ (56)

We recall that our rod is naturally curved and twisted:
${\bm{\Omega}}^{(0)}=\Omega_{a}^{(0)}\hat{\bm{d}}_{a}^{(0)}$. We expand the
Cosserat frame basis in Eq. (42) around the reference configuration
$\hat{\bm{d}}_{a}^{(0)}$. We will determine the relationship between the
first-order perturbation of the components of the Darboux vector
$\Omega_{a}^{(1)}$ defined as

$\displaystyle{\Omega}_{a}={\Omega}^{(0)}_{a}+{\Omega}^{(1)}_{a}+\cdots\,,$
(57)

and the first order perturbation of $\hat{\bm{d}}_{a}$, $\alpha_{a}$. Note
that $\Omega_{a}^{(1)}$ and $\alpha_{a}$ are infinitesimal quantities of the
same order $\Omega_{a}^{(1)}\sim\alpha_{a}$. We substitute Eqs. (42)

$\displaystyle\hat{\bm{d}}_{a}=(\delta_{ab}+\varepsilon_{abc}\alpha_{c})\hat{\bm{d}}_{b}^{(0)}\,,$
(58)

and (57) into Eq. (1);
$\hat{\bm{d}}_{a}^{\prime}=\bm{\Omega}\times\hat{\bm{d}}_{a}$, to derive Eq.
(43). The final result is consistent with that in Ref. [65] (which is based on
matrix representation), while our derivation is based on a tensor
representation.

First, the kinematic equation Eq. (1) is rewritten using the Eddington epsilon
as

$\displaystyle\hat{\bm{d}}_{a}^{\prime}=-\varepsilon_{abc}\Omega_{b}\hat{\bm{d}}_{c},$
(59)

which also holds for the base solution as
$\hat{\bm{d}}_{a}^{(0)^{\prime}}=-\varepsilon_{abc}\Omega_{b}^{(0)}\hat{\bm{d}}_{c}^{(0)}$.
The left-hand side of Eq. (59) is computed as

$\displaystyle\hat{\bm{d}}_{a}^{\prime}$ $\displaystyle=$
$\displaystyle\\{(\delta_{ab}+\varepsilon_{abc}\alpha_{c})\hat{\bm{d}}_{b}^{(0)}\\}^{\prime}$
(60) $\displaystyle=$
$\displaystyle-\varepsilon_{abc}\Omega_{b}^{(0)}\hat{\bm{d}}_{c}^{(0)}+\varepsilon_{abc}\alpha_{c}^{\prime}\hat{\bm{d}}_{b}^{(0)}-\varepsilon_{abc}\varepsilon_{bde}\alpha_{c}\Omega_{d}^{(0)}\hat{\bm{d}}_{e}^{(0)}$
$\displaystyle=$
$\displaystyle-\varepsilon_{abc}\Omega_{b}^{(0)}\hat{\bm{d}}_{c}^{(0)}+\varepsilon_{abc}\alpha_{c}^{\prime}\hat{\bm{d}}_{b}^{(0)}-(\delta_{cd}\delta_{ae}-\delta_{ce}\delta_{ad})\alpha_{c}\Omega_{d}^{(0)}\hat{\bm{d}}_{e}^{(0)}$
$\displaystyle=$
$\displaystyle-\varepsilon_{abc}\Omega_{b}^{(0)}\hat{\bm{d}}_{c}^{(0)}+\varepsilon_{abc}\alpha_{c}^{\prime}\hat{\bm{d}}_{b}^{(0)}-(\alpha_{c}\Omega_{c}^{(0)})\hat{\bm{d}}_{a}^{(0)}+\Omega_{a}^{(0)}(\alpha_{c}\hat{\bm{d}}_{c}^{(0)}).$

We expand the right-hand side of Eq. (59) as

$\displaystyle-\varepsilon_{abc}\Omega_{b}\hat{\bm{d}}_{c}$ $\displaystyle=$
$\displaystyle-\varepsilon_{abc}({\Omega}^{(0)}_{b}+{\Omega}^{(1)}_{b})(\delta_{cd}+\varepsilon_{cde}\alpha_{e})\hat{\bm{d}}_{d}^{(0)}$
(61) $\displaystyle=$
$\displaystyle-\varepsilon_{abc}\Omega_{b}^{(0)}\hat{\bm{d}}_{c}^{(0)}-\varepsilon_{abc}\Omega_{b}^{(1)}\hat{\bm{d}}_{c}^{(0)}-\varepsilon_{abc}\varepsilon_{cde}\alpha_{e}{\Omega}^{(0)}_{b}\hat{\bm{d}}_{d}^{(0)}$
$\displaystyle=$
$\displaystyle-\varepsilon_{abc}\Omega_{b}^{(0)}\hat{\bm{d}}_{c}^{(0)}-\varepsilon_{abc}\Omega_{b}^{(1)}\hat{\bm{d}}_{c}^{(0)}-(\delta_{ad}\delta_{be}-\delta_{ae}\delta_{bd})\alpha_{e}{\Omega}^{(0)}_{b}\hat{\bm{d}}_{d}^{(0)}$
$\displaystyle=$
$\displaystyle-\varepsilon_{abc}\Omega_{b}^{(0)}\hat{\bm{d}}_{c}^{(0)}-\varepsilon_{abc}\Omega_{b}^{(1)}\hat{\bm{d}}_{c}^{(0)}-(\alpha_{c}\Omega_{c}^{(0)})\hat{\bm{d}}_{a}^{(0)}+\alpha_{a}(\Omega_{c}^{(0)}\hat{\bm{d}}_{c}^{(0)}).$

Note that, in Eqs. (60) and (61), we dropped higher order terms (e.g. the
terms proportional to $\alpha_{a}\Omega_{b}^{(1)}$).

Substituting Eqs. (60) and (61) into Eq. (59), we obtain the tensorial
relation

$\displaystyle\varepsilon_{abc}\alpha_{c}^{\prime}\hat{\bm{d}}_{b}^{(0)}=-\varepsilon_{abc}\Omega_{b}^{(1)}\hat{\bm{d}}_{c}^{(0)}+\alpha_{a}(\Omega_{c}^{(0)}\hat{\bm{d}}_{c}^{(0)})-\Omega_{a}^{(0)}(\alpha_{c}\hat{\bm{d}}_{c}^{(0)}).$
(62)

Next, by taking the inner product between Eq. (62) and
$\hat{\bm{d}}_{d}^{(0)}$, and by arranging the indices, we arrive at

$\displaystyle\varepsilon_{abc}\alpha_{c}^{\prime}=\varepsilon_{abc}\Omega_{c}^{(1)}+\left(\alpha_{a}\Omega_{b}^{(0)}-\Omega_{a}^{(0)}\alpha_{b}\right).$
(63)

The presence of the term inside the parenthesis in the right-hand side of Eq.
(63) encodes the geometric non-linearity intrinsic to a naturally curved rod.
When the rod is naturally straight, $\Omega_{a}^{(0)}=0$, we recover Eq. (43).

##  Appendix B Derivation of Equations (49)-(51)

In this appendix, we derive the Cosserat frame basis for the uniform helix
Eqs. (49)-(51). To do so, we will integrate the kinematic relation Eq. (1)
with ${\bm{\Omega}}=\kappa\hat{\bm{d}}_{2}+\tau\hat{\bm{d}}_{3}$. The goal is
to express $\bm{r}(s)$ and $\hat{\bm{d}}_{a}~{}(a=1,2,3)$ in terms of the
Cartesian basis $\hat{\bm{e}}_{i}$. From the kinematic relation, we get

$\displaystyle\hat{\bm{d}}_{1}^{\prime}$ $\displaystyle=$
$\displaystyle-\kappa\hat{\bm{d}}_{3}+\tau\hat{\bm{d}}_{2},$ (64)
$\displaystyle\hat{\bm{d}}_{2}^{\prime}$ $\displaystyle=$
$\displaystyle-\tau\hat{\bm{d}}_{1},$ (65)
$\displaystyle\hat{\bm{d}}_{3}^{\prime}$ $\displaystyle=$
$\displaystyle\kappa\hat{\bm{d}}_{1},$ (66)

which are the same as the Frenet-Serret equations [43]. Differentiating Eq.
(64) with respect to $s$, and with the help of Eqs. (65) and (66), we find the
ODE only for $\hat{\bm{d}}_{1}$;

$\displaystyle\hat{\bm{d}}_{1}^{\prime\prime}=-K^{2}\hat{\bm{d}}_{1}\,.$ (67)

Here, we introduced the normalized curvature $K$ as
$K=\sqrt{\kappa^{2}+\tau^{2}}$, or equivalently

$\displaystyle\kappa=K\sin\psi,~{}~{}\tau=K\cos\psi.$ (68)

After integrating Eq.( 67) twice, we compute $\hat{\bm{d}}_{1}$ as

$\displaystyle\hat{\bm{d}}_{1}(s)=\cos(Ks)\hat{\bm{d}}_{1}(0)+\frac{\tau}{K}\sin(Ks)\hat{\bm{d}}_{2}(0)-\frac{\kappa}{K}\sin(Ks)\hat{\bm{d}}_{3}(0).$
(69)

Following a similar procedure for $\hat{\bm{d}}_{2}$ and $\hat{\bm{d}}_{3}$,
we obtain

$\displaystyle\hat{\bm{d}}_{2}(s)-\hat{\bm{d}}_{2}(0)$ $\displaystyle=$
$\displaystyle-\frac{\tau}{K}\sin(Ks)\hat{\bm{d}}_{1}(0)-\frac{\tau}{K^{2}}(1-\cos(Ks))(\tau\hat{\bm{d}}_{2}(0)-\kappa\hat{\bm{d}}_{3}(0)),$
(70) $\displaystyle\hat{\bm{d}}_{3}(s)-\hat{\bm{d}}_{3}(0)$ $\displaystyle=$
$\displaystyle\frac{\kappa}{K}\sin(Ks)\hat{\bm{d}}_{1}(0)+\frac{\kappa}{K^{2}}(1-\cos(Ks))(\tau\hat{\bm{d}}_{2}(0)-\kappa\hat{\bm{d}}_{3}(0)).$
(71)

By integrating $\boldsymbol{r}^{\prime}=\hat{\bm{d}}_{3}$, the centerline
position is obtained as

$\displaystyle\boldsymbol{r}(s)-\boldsymbol{r}(0)=\frac{\kappa}{K^{2}}(1-\cos(Ks))\hat{\bm{d}}_{1}(0)-\frac{\kappa}{K^{3}}\sin(Ks)(\tau\hat{\bm{d}}_{2}(0)-\kappa\hat{\bm{d}}_{3}(0))+\frac{\tau
s}{K^{2}}(\kappa\hat{\bm{d}}_{2}(0)+\tau\hat{\bm{d}}_{3}(0)),$ (72)

where the unit vector of the central axis of the helix $\hat{N}$ is

$\displaystyle\hat{N}\equiv\frac{\kappa}{K}\hat{\bm{d}}_{2}(0)+\frac{\tau}{K}\hat{\bm{d}}_{3}(0).$
(73)

Finally, let us rewrite Eq. (72) in a more compact form by choosing
$\hat{\bm{d}}_{a}(0)$ and $\boldsymbol{r}(0)$ such that the central axis lies
along $-\hat{\bm{e}}_{z}$; i.e., we select $\hat{N}=-\hat{\bm{e}}_{z}$,
$\hat{\bm{d}}_{1}(0)=\hat{\bm{e}}_{x}$, and
$\hat{\bm{e}}_{y}=\hat{\bm{e}}_{z}\times\hat{\bm{e}}_{x}=-(\tau\hat{\bm{d}}_{2}(0)-\kappa\hat{\bm{d}}_{3}(0))/K$.
Then, using Eq. (68), the centerline position is calculated as

$\displaystyle\boldsymbol{r}(s)=-\frac{\sin\psi}{K}\cos(Ks)\hat{\bm{e}}_{x}+\frac{\sin\psi}{K}\sin(Ks)\hat{\bm{e}}_{y}-s\cos\psi\hat{\bm{e}}_{z},$
(74)

which implies

$\displaystyle R=\frac{\sin\psi}{K}.$ (75)

The normalized curvature $K$ is now rewritten as a function of $\kappa$ and
$\psi$ only; $K=\kappa\sqrt{1+\tan^{-2}\psi}=\kappa/\sin\psi$, i.e.
$\kappa=\sin^{2}\psi/R$ and $\tau=\sin\psi\cos\psi/R=\sin 2\psi/2R$. The
choice of $\hat{N}$ allows us to obtain the Cosserat frame basis vectors as

$\displaystyle\hat{\bm{d}}_{1}(s)$ $\displaystyle=$
$\displaystyle\cos(Ks)\hat{\bm{e}}_{x}-\sin(Ks)\hat{\bm{e}}_{y},$ (76)
$\displaystyle\hat{\bm{d}}_{2}(s)$ $\displaystyle=$
$\displaystyle-\cos\psi\left(\sin(Ks)\hat{\bm{e}}_{x}+\cos(Ks)\hat{\bm{e}}_{y}\right)-\sin\psi\hat{\bm{e}}_{z},$
(77) $\displaystyle\hat{\bm{d}}_{3}(s)$ $\displaystyle=$
$\displaystyle\sin\psi\left(\sin(Ks)\hat{\bm{e}}_{x}+\cos(Ks)\hat{\bm{e}}_{y}\right)-\cos\psi\hat{\bm{e}}_{z},$
(78)

thereby reproducing Eqs. (49)-(51) that we set out to derive.

## References

  * Li et al. [2009] W. Li, K. Kostidis, Xianzhou Zhang, Yang Zhou,  Development of a force sensor working with mr elastomers,  in: 2009 IEEE/ASME International Conference on Advanced Intelligent Mechatronics, pp. 233–238.
  * Dorfmann and Ogden [2014] L. Dorfmann, R. W. Ogden, Nonlinear Theory of Electroelastic and Magnetoelastic Interactions, Springer, New York, 2014.
  * Menzel [2019] A. M. Menzel,  Mesoscopic characterization of magnetoelastic hybrid materials: magnetic gels and elastomers, their particle-scale description, and scale-bridging links,  Archive of Applied Mechanics 89 (2019) 17–45.
  * Ginder et al. [1999] J. M. Ginder, M. E. Nichols, L. D. Elie, J. L. Tardiff,  Magnetorheological elastomers: properties and applications,  in: M. R. Wuttig (Ed.), Smart Structures and Materials 1999: Smart Materials Technologies, International Society for Optics and Photonics, 1999, pp. 131–138.
  * Huang et al. [2016] H.-W. Huang, M. S. Sakar, K. Riederer, N. Shamsudhin, A. Petruska, S. Pané, B. J. Nelson,  Magnetic Microrobots with Addressable Shape Control ,  IEEE International Conference on Robotics and Automation ICRA Stockholm, Sweden (2016).
  * Hu et al. [2018] W. Hu, G. Z. Lum, M. Mastrangeli, M. Sitti,  Small-scale soft-bodied robot with multimodal locomotion,  Nature 554 (2018) 81–85.
  * Pancaldi et al. [2020] L. Pancaldi, P. Dirix, A. Fanelli, A. M. Lima, N. Stergiopulos, P. J. Mosimann, D. Ghezzi, M. S. Sakar,  Flow driven robotic navigation of microengineered endovascular probes,  Nature Communications 11 (2020) 1–14.
  * Moon and Pao [1968] F. C. Moon, Y.-H. Pao,  Magnetoelastic Buckling of a Thin Plate,  Journal of Applied Mechanics 35 (1968) 53–58.
  * Cebers [2003] A. Cebers,  Dynamics of a chain of magnetic particles connected with elastic linkers,  Journal of Physics: Condensed Matter 15 (2003) S1335–S1344.
  * Cēbers and Javaitis [2004] A. Cēbers, I. Javaitis,  Bending of flexible magnetic rods,  Physical Review E 70 (2004) 021404.
  * Dreyfus et al. [2005] R. Dreyfus, J. Baudry, M. L. Roper, M. Fermigier, H. A. Stone, J. Bibette,  Microscopic artificial swimmers,  Nature 437 (2005) 862–865.
  * Roper et al. [2006] M. Roper, R. Dreyfus, J. Baudry, M. Fermigier, J. Bibette, H. A. Stone,  On the dynamics of magnetically driven elastic filaments,  Journal of Fluid Mechanics 554 (2006) 167–190.
  * Cēbers and Cīrulis [2007] A. Cēbers, T. Cīrulis,  Magnetic elastica,  Physical Review E 76 (2007) 031504.
  * Gerbal et al. [2015] F. Gerbal, Y. Wang, F. Lyonnet, J.-C. Bacri, T. Hocquet, M. Devaud,  A refined theory of magnetoelastic buckling matches experiments with ferromagnetic and superparamagnetic rods,  Proceedings of the National Academy of Sciences 112 (2015) 7135–7140.
  * Hall et al. [2013] C. L. Hall, D. Vella, A. Goriely,  The Mechanics of a Chain or Ring of Spherical Magnets,  SIAM Journal on Applied Mathematics 73 (2013) 2029–2054.
  * Vella et al. [2014] D. Vella, E. du Pontavice, C. L. Hall, A. Goriely,  The magneto-elastica: from self-buckling to self-assembly,  Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 470 (2014) 20130609.
  * Schönke and Fried [2017] J. Schönke, E. Fried,  Stability of vertical magnetic chains,  Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 473 (2017) 20160703.
  * Rigbi and Jilkén [1983] Z. Rigbi, L. Jilkén,  The response of an elastomer filled with soft ferrite to mechanical and magnetic influences,  Journal of Magnetism and Magnetic Materials 37 (1983) 267–276.
  * Dorfmann and Ogden [2003] A. Dorfmann, R. W. Ogden,  Magnetoelastic modelling of elastomers,  European Journal of Mechanics - A/Solids 22 (2003) 497–507.
  * Danas et al. [2012] K. Danas, S. V. Kankanala, N. Triantafyllidis,  Experiments and modeling of iron-particle-filled magnetorheological elastomers,  Journal of the Mechanics and Physics of Solids 60 (2012) 120–138.
  * Loukaides et al. [2015] E. G. Loukaides, S. K. Smoukov, K. A. Seffen,  Magnetic actuation and transition shapes of a bistable spherical cap,  International Journal of Smart and Nano Materials 5 (2015) 270–282.
  * Seffen and Vidoli [2016] K. A. Seffen, S. Vidoli,  Eversion of bistable shells under magnetic actuation: a model of nonlinear shapes,  Smart Materials and Structures 25 (2016) 065010.
  * Schmauch et al. [2017] M. M. Schmauch, S. R. Mishra, B. A. Evans, O. D. Velev, J. B. Tracy,  Chained Iron Microparticles for Directionally Controlled Actuation of Soft Robots,  ACS Applied Materials & Interfaces 9 (2017) 11895–11901.
  * Ciambella et al. [2018] J. Ciambella, A. Favata, G. Tomassetti,  A nonlinear theory for fibre-reinforced magneto-elastic rods,  Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 474 (2018) 20170703.
  * Singh and Onck [2018] R. P. Singh, P. R. Onck,  Magnetic field induced deformation and buckling of slender bodies,  International Journal of Solids and Structures 143 (2018) 29–58.
  * Psarra et al. [2019] E. Psarra, L. Bodelot, K. Danas,  Wrinkling to crinkling transitions and curvature localization in a magnetoelastic film bonded to a non-magnetic substrate,  Journal of the Mechanics and Physics of Solids 133 (2019) 103734.
  * Alapan et al. [2020] Y. Alapan, A. C. Karacakol, S. N. Guzelhan, I. Isik, M. Sitti,  Reprogrammable shape morphing of magnetic soft machines,  Science Advances 6 (2020) eabc6414.
  * Lum et al. [2016] G. Z. Lum, Z. Ye, X. Dong, H. Marvi, O. Erin, W. Hu, M. Sitti,  Shape-programmable magnetic soft matter,  Proceedings of the National Academy of Sciences 113 (2016) E6007–E6015.
  * Kim et al. [2018] Y. Kim, H. Yuk, R. Zhao, S. A. Chester, X. Zhao,  Printing ferromagnetic domains for untethered fast-transforming soft materials,  Nature 558 (2018) 274–279.
  * Zhao et al. [2019] R. Zhao, Y. Kim, S. A. Chester, P. Sharma, X. Zhao,  Mechanics of hard-magnetic soft materials,  Journal of the Mechanics and Physics of Solids 124 (2019) 244–263.
  * Wang et al. [2020] L. Wang, Y. Kim, C. F. Guo, X. Zhao,  Hard-magnetic elastica,  Journal of the Mechanics and Physics of Solids 142 (2020) 104045.
  * Ciambella and Tomassetti [2020] J. Ciambella, G. Tomassetti,  A form-finding strategy for magneto-elastic actuators,  International Journal of Non-Linear Mechanics 119 (2020) 103297.
  * Chen et al. [2021] T. Chen, M. Pauly, P. M. Reis,  A reprogrammable mechanical metamaterial with stable memory,  Nature 589 (2021) 386–390.
  * Yan et al. [2021a] D. Yan, M. Pezzulla, L. Cruveiller, A. Abbasi, P. M. Reis,  Magneto-active elastic shells with tunable buckling strength,  Nature Communications 12 (2021a) 1–9.
  * Yan et al. [2021b] D. Yan, A. Abbasi, P. M. Reis,  A comprehensive framework for hard-magnetic beams: reduced-order theory, 3D simulations, and experiments,  arXiv:2106.14878[cond-mat.soft] (2021b).
  * Diller et al. [2014] E. Diller, J. Zhuang, G. Z. Lum, M. R. Edwards, M. Sitti,  Continuously distributed magnetization profile for millimeter-scale elastomeric undulatory swimming,  Applied Physics Letters 104 (2014) 174101.
  * Hu et al. [2018] N. Hu, L. Wang, W. Zhai, M. Sun, H. Xie, Z. Wu, Q. He,  Magnetically Actuated Rolling of Star-Shaped Hydrogel Microswimmer,  Macromolecular Chemistry and Physics 219 (2018) 1700540.
  * Zhang and Diller [2018] J. Zhang, E. Diller,  Untethered Miniature Soft Robots: Modeling and Design of a Millimeter-Scale Swimming Magnetic Sheet,  Soft Robotics 5 (2018) 761–776.
  * Tsumori et al. [2015] F. Tsumori, A. Saijou, T. Osada, H. Miura,  Development of actuation system for artificial cilia with magnetic elastomer,  Japanese Journal of Applied Physics 54 (2015) 06FP12.
  * Pece et al. [2017] F. Pece, J. J. Zarate, V. Vechev, N. Besse, O. Gudozhnik, H. Shea, O. Hilliges,  MagTics,  in: UIST ’17: The 30th Annual ACM Symposium on User Interface Software and Technology, ACM, New York, NY, USA, 2017, pp. 143–154.
  * Wang [1986] C. Y. Wang,  A critical review of the heavy elastica,  International Journal of Mechanical Sciences 28 (1986) 549–559.
  * Dill [1992] E. H. Dill,  Kirchhoff’s theory of rods,  Archive for History of Exact Sciences 44 (1992) 1–23.
  * Audoly and Pomeau [2010] B. Audoly, Y. Pomeau, Elasticity and Geometry, Oxford Univ. Press, 2010\.
  * Landau and Lifshitz [1980] L. D. Landau, E. M. Lifshitz, Theory of Elasticity, Pergamon Press, 1980.
  * Powers [2010] T. R. Powers,  Dynamics of filaments and membranes in a viscous fluid,  Reviews of Modern Physics 82 (2010) 1607–1631.
  * Bigoni [2015] D. Bigoni, Extremely Deformable Structures (CISM Lecture Notes No. 562), Springer-Verlag, 2015.
  * Goriely [2017] A. Goriely, The Mathematics and Mechanics of Biological Growth, Springer-Verlag, 2017\.
  * Chirico and Langowski [1994] G. Chirico, J. Langowski,  Kinetics of DNA supercoiling studied by Brownian dynamics simulation,  Biopolymers 34 (1994) 415–433.
  * Wada and Netz [2007a] H. Wada, R. R. Netz,  Stretching helical nano-springs at finite temperature,  Europhysics Letters 77 (2007a) 68001.
  * Wada and Netz [2007b] H. Wada, R. R. Netz,  Model for Self-Propulsive Helical Filaments: Kink-Pair Propagation,  Physical Review Letters 99 (2007b) 108102.
  * Vogel and Stark [2010] R. Vogel, H. Stark,  Force-extension curves of bacterial flagella,  The European Physical Journal E 33 (2010) 259–271.
  * Reichert [2006] M. Reichert, Hydrodynamic Interactions in Colloidal and Biological Systems, Ph.D. thesis, University Konstanz, 2006.
  * Bergou et al. [2008] M. Bergou, M. Wardetzky, S. Robinson, B. Audoly, E. Grinspun,  Discrete elastic rods,  in: ACM SIGGRAPH 2008 papers, ACM Press, New York, New York, USA, 2008, p. 1.
  * Morigaki et al. [2016] Y. Morigaki, H. Wada, Y. Tanaka,  Stretching an Elastic Loop: Crease, Helicoid, and Pop Out,  Physical Review Letters 117 (2016) 198003.
  * Sano et al. [2017] T. G. Sano, T. Yamaguchi, H. Wada,  Slip Morphology of Elastic Strips on Frictional Rigid Substrates,  Physical Review Letters 118 (2017) 178001.
  * Sano and Wada [2018] T. G. Sano, H. Wada,  Snap-buckling in asymmetrically constrained elastic strips,  Physical Review E 97 (2018) 013002.
  * Sano and Wada [2019] T. G. Sano, H. Wada,  Twist-Induced Snapping in a Bent Elastic Rod and Ribbon,  Physical Review Letters 122 (2019) 114301.
  * Jackson [1975] J. D. Jackson, Classical electrodynamics; 2nd ed., Wiley, New York, NY, 1975\.
  * Garcia-Gonzalez and Hossain [2021] D. Garcia-Gonzalez, M. Hossain,  Microstructural modelling of hard-magnetic soft materials: Dipoleâdipole interactions versus zeeman effect,  Extreme Mechanics Letters 48 (2021) 101382.
  * Doedel [1981] E. Doedel,  AUTO: A program for the automatic bifurcation analysis of autonomous systems,  Congressus Numerantium 30 (1981).
  * Goldstein et al. [2002] H. Goldstein, C. Poole Jr., J. Safko, Classical mechanics 3rd. Edition, Addison-Wesley, 2002.
  * Hairer et al. [1993] E. Hairer, S. P. Nørsett, G. Wanner, Solving ordinary differential equations I: Nonstiff problems (2nd ed.), Springer Verlag, Berlin, 1993.
  * Brazier [1927] L. G. Brazier,  On the flexure of thin cylindrical shells and other ”thin” sections,  Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 116 (1927) 104–114.
  * Counto [1964] U. J. Counto,  The effect of the elastic modulus of the aggregate on the elastic modulus, creep and creep recovery of concrete,  Magazine of Concrete Research 16 (1964) 129–138.
  * Goriely and Tabor [1997] A. Goriely, M. Tabor,  Nonlinear dynamics of filaments. III. Instabilities of helical rods,  Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 453 (1997) 2583–2601.
  * Strogatz [2018] S. H. Strogatz, Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering, 2nd Ed., CRC Press, Boca Raton, FL, 2018\.