# The chain length of anisotropic paramagnetic particles in a rotating field

Jānis Užulis  Jānis Cīmurs<EMAIL_ADDRESS>mmml.lu.lv Laboratory of
Magnetic Soft Materials
Faculty of Physics, Mathematics and Optometry
University of Latvia
Jelgavas iela 3-530, Riga, Latvia, LV-1004

###### Abstract

In this article the maximal length of a chain of paramagnetic particles with
magnetic anisotropy in a rotating magnetic field is studied. The theory of
paramagnetic particle chains usually assumes that the particles are
magnetically isotropic and do not rotate in a rotating field. In experiments
it is seen that spherical paramagnetic particles rotate, which can be
explained by small magnetic anisotropy. In this article, the maximal chain
length is calculated for paramagnetic particles with magnetic anisotropy in a
rotating magnetic field. Results show that the maximal chain length as a
function of field frequency has the same trend for isotropic magnetic
particles and particles with magnetic anisotropy if the field frequency is
much higher or much lower than the critical frequency of an individual
particle. blue Initially randomly distributed particles will form chains that
will collide and exchange with particles till they obtain a typical chain
length.The typical chain length of a small cluster is shorter than the maximal
chain length of an isolated chain for the same field frequency. The
distribution of chain lengths in a small cluster of chains is narrower for
particles with higher magnetic anisotropy. Due to the narrower distribution of
chain lengths, particles with magnetic anisotropy can suit better for mass-
production. This article will show how magnetic anisotropy parameters of
paramagnetic particles influence chain length of chains which form in a
rotating magnetic field.

††preprint: APS/123-QED

##  I Introduction

It is shown that paramagnetic particles in a rotating magnetic field become
attractive [1, 2]. If the field frequency is small, particles form chains [1,
3], where the chain length depends on the field frequency [4, 5, 6, 7] and
fluid and particle parameters.

Magnetic particles chains can be used for cargo transport [2, 8]. The
advantage of the magnetic particles in cargo transport is that the motion can
be fully controlled from the outside using a pattern of an external magnetic
field. Functionalized magnetic particles can be linked together creating
magnetic filament [9, 10] which can be used as magnetic swimmer [11]. For
these applications, the possibility of creating a large amount of chains with
well defined length is advantageous.

Magnetic fluid, which consists of superparamagnetic particles suspended in a
liquid, also exhibits a similar formation of chains in a rotating magnetic
field [12, 13]. Despite the fact that in a phase separated magnetic fluid
surface tension and thermal motion of the particles play an important role,
the behaviour of magnetic fluid in a rotating field is similar to the
behaviour of a small cluster of magnetic particles. In slow rotating magnetic
fields, magnetic fluid forms elongated drops with narrow distribution of
lengths [12].

The theory of spherical paramagnetic particles in a rotating magnetic field
usually does not take into account magnetic anisotropy of particles, therefore
particles are assumed non-rotating. In reality it can be seen that
paramagnetic particles rotate in a rotating field [2, 14]. In [2, 14],
rotation of spherical paramagnetic particles is explained by finite relaxation
time of magnetic moment. The rotation of the paramagnetic particles can be
explained also by magnetic anisotropy of the susceptibility of the particles.

In this article, chain formation of paramagnetic particles with anisotropy in
a rotating magnetic field is investigated. The critical chain length is
compared to theory for chains of isotropic paramagnetic particles to see if
measurable difference can be observed. Also, the hypothesis that particles
with magnetic anisotropy have more narrow chain length distribution in
initially randomly distributed cluster of particles is tested. The results and
conclusions show advantages for using paramagnetic particles with anisotropy
to form paramagnetic chains.

##  II Model

In this article, we will describe paramagnetic particles with spherical form
and with uniaxial magnetic anisotropy. The magnetic anisotropy could be caused
by the crystallographic structure [15] or by the fact that particles are not
ideally spherical or a particle could be a magnetic rod in a spherical shell
[16] or due to magnetic Janus particles, where some regions have different
magnetic susceptibility [17]. When the particles are subjected to a weak
magnetic field $\vec{H}$, the magnetic dipole moment $\vec{m}_{i}$ is induced
in each particle. The dipole moment of $i$-th particle can be calculated using
a linear particle magnetization and field relation:

$\vec{m}_{i}=V_{m}\left(\chi_{\perp}\vec{H}_{i}+\Delta\chi(\vec{n}_{i}\cdot\vec{H}_{i})\vec{n}_{i}\right)\text{
,}$ (1)

where $V_{m}$ is the magnetic volume of the particle, $\vec{H}_{i}$ is the
magnetic field at the point of the particle, $\vec{n}_{i}$ is the unit vector
in the direction of the anisotropy axis of the particle,
$\Delta\chi=\chi_{\|}-\chi_{\perp}$ is difference of the magnetic
susceptibilities, where $\chi_{\|}$ and $\chi_{\perp}$ are magnetic
susceptibility of the particle in the direction of anisotropy axis and
perpendicular to it respectively. For ellipsoidal paramagnetic particles,
susceptibilities $\chi_{\perp}$ and $\chi_{\|}$ are related to material
susceptibility $\chi$, as shown by Stoner and Osborn [18, 19].
$\chi_{\|,\perp}=\chi_{i}/(1-N_{\|,\perp}\chi_{i})$, where $\chi_{i}$ is
material susceptibility in the given direction and $N_{\|,\perp}$ is a
demagnetization factor which depends on the shape of the particle. For
example, an isotropic spherical particle with large magnetic susceptibility
has $N_{\|}=N_{\perp}=1/3$ and $\chi_{\|}=\chi_{\perp}=3$. Parameter
$\frac{\chi_{\|}}{\chi_{\perp}}$ is used to characterize the anisotropy of the
particle in this article. It is assumed that magnetic relaxation time is much
shorter than the rotation period of the external magnetic field. As a result,
the magnetic moment calculated by equation (1) does not depend on the rotation
frequency. The direction of magnetization $\vec{m}$ is illustrated in Figure 1
for different external fields $\vec{H}$ and anisotropy axis $\vec{n}$
configurations.

Figure 1:  Magnetic moment $\vec{m}$ direction depending on the magnetic field
$\vec{H}$ direction. The magnetization $|\vec{m}|$ is normalized to the field
strength $|\vec{H}|$ for each anisotropy case. In upper two examples $\vec{n}$
is an easy magnetization axis ($\frac{\chi_{\|}}{\chi_{\perp}}=2$,
$\frac{\chi_{\|}}{\chi_{\perp}}=1.5$) and in the lower one example $\vec{n}$
is a hard magnetization axis ($\frac{\chi_{\|}}{\chi_{\perp}}=0.5$). The
direction of the magnetization $\vec{m}$ and the external field $\vec{H}$
direction coincide when the direction of $\vec{H}$ is along the anisotropy
axis $\vec{n}$ or perpendicular to it.

Figure 1 shows that the magnetization $\vec{m}$ is closer to the anisotropy
axis if $\chi_{\|}>\chi_{\perp}$ and the anisotropy axis is an easy axis. If
$\chi_{\perp}>\chi_{\|}$ the anisotropy axis is a hard axis and magnetization
of the particle is closer to the equator of the particle.

In the model these particles are dispersed in liquid with viscosity $\eta$.
Each particle rotates due to torque
$\vec{\tau}_{i}=\mu_{0}\vec{m}_{i}\times\vec{H}_{i}$, which is compensated by
viscous drag. The rotation of the particle can be calculated as [20]:

$\dfrac{\differential\vec{n}}{\differential
t}=\dfrac{\mu_{0}V_{m}\Delta\chi}{\xi_{r}}\left(\vec{n}\cdot\vec{H}_{i}\right)\left[\vec{H}_{i}-\left(\vec{n}\cdot\vec{H}_{i}\right)\vec{n}\right]\text{
,}$ (2)

where $\vec{n}$ is the magnetic anisotropy axis, which is fixed to the
particle, $\mu_{0}$ is the magnetic constant (vacuum permeability),
$\xi_{r}=8\pi\eta R^{3}$ is the rotational drag coefficient, where $R$ is the
hydrodynamic radius of the particle and $\eta$ is viscosity of the liquid. It
is assumed that viscosity dominates over inertia and Reynolds number is below
unity and inertial terms can be neglected.

The magnetic field at the point of the particle $\vec{H}_{i}$ is a sum of the
external field $\vec{H}$ and the field produced by all other particles. $j$-th
particle at the point of $i$-th particle produces field:

$\displaystyle\vec{H}_{ij}(\vec{r}_{ij})=\frac{1}{4\pi}\left({\frac{3\vec{r}_{ij}(\vec{m_{j}}\cdot\vec{r}_{ij})}{r_{ij}^{5}}}-{\frac{\vec{m_{j}}}{r_{ij}^{3}}}\right)\text{
,}$ (3)

where $\vec{r}_{ij}=\vec{r}_{i}-\vec{r}_{j}$ is distance vector from the point
of $j$-th particle to the point of $i$-th particle. The magnetic field at the
point of the $i$-th particle is
$\vec{H}_{i}=\vec{H}+\sum\limits_{j}\vec{H}_{ij}$

Since the field produced by particles is not homogeneous, particles also move.
The force exerted by $j$-th particle on the $i$-th particle can be calculated
as

$\begin{split}\displaystyle\vec{F}_{ij,M}&={\dfrac{3\mu_{0}}{4\pi
r_{ij}^{5}}}\bigg{[}(\vec{m}_{i}\cdot\vec{r}_{ij})\vec{m}_{j}+(\vec{m}_{j}\cdot\vec{r}_{ij})\vec{m}_{i}\\\
&+(\vec{m}_{i}\cdot\vec{m}_{j})\vec{r}_{ij}-{\dfrac{5(\vec{m}_{i}\cdot\vec{r}_{ij})(\vec{m}_{j}\cdot\vec{r}_{ij})}{r^{2}_{ij}}}\vec{r}_{ij}\bigg{]}\text{
.}\end{split}$ (4)

The non-overlap of the particles is fulfilled using the repulsive part of
Lenard-Jones potential. The repulsive force is:

$\vec{F}_{ij,R}=\begin{cases}-G\left(\frac{1}{r_{ij}^{13}}-\frac{r_{ij}-2R}{r_{max}-2R}\frac{1}{r_{max}^{13}}\right)\frac{\vec{r}_{ij}}{r_{ij}}\text{
if }r_{ij}\leq r_{max}\\\ 0\text{ if }r_{ij}>r_{max}\\\ \end{cases}\text{.}$
(5)

The repulsive force is proportional to the particle distance in power of $-13$
as obtained from the Lenard-Jones potential repulsive part. Constant $G$ is
chosen so that the maximal attractive magnetic force $\vec{F}_{ij,M}$ between
two particles at the distance of $|r_{ij}|=2R$ between particle centers is
compensated by repulsive force $\vec{F}_{ij,R}$ and
$\vec{F}_{ij,M,max}+\vec{F}_{ij,R}=0$. The maximal attractive magnetic force
is obtained in a situation where $\vec{n}_{i}\|\vec{H}_{i}$ and
$\chi_{\|}>\chi_{\perp}$ and $\vec{r}_{ij}\|\vec{m}_{i}\|\vec{m}_{j}$. The
maximal magnetic force is
$\vec{F}_{ij,M,max}=\frac{3V_{m}^{2}\chi_{\|}^{2}H_{i}^{2}\mu_{0}}{32\pi
R^{4}}$ and consequently $G=F_{ij,M,max}(2R)^{13}$. The parameter $r_{max}$ is
a bit larger than $2R$ and determines the softness of the particle. The
equation part $\frac{r_{ij}-2R}{r_{max}-2R}\frac{1}{r_{max}^{13}}$ makes the
repulsive force 0 at the particle distance $r_{ij}=r_{max}$ and makes the
repulsion function without steps, to avoid numerical instabilities.

The sum of all forces, the magnetic force (4) and the repulsive force (5),
will cause the particle to move. In the non-inertial limit the driving force
is compensated by drag force for spherical particle $\vec{F}_{D}=-6\pi\eta
R\vec{v}$, where $\vec{v}$ is velocity of the particle. Velocity of the $i$-th
particle is calculated as:

$\vec{v}_{i}=\frac{1}{6\pi\eta R}\sum\limits_{j}\vec{F}_{ij},$ (6)

where $\vec{F}_{ij}=\vec{F}_{ij,M}+\vec{F}_{ij,R}$ is a sum of all the forces.

Magnetic particles are subjected by an external rotating magnetic field:
$\vec{H}=H(\cos(\omega_{H}t),\sin(\omega_{H}t),0)$, where $\omega_{H}$ is a
field rotation frequency.

In this article thermal fluctuations are omitted. Thermal fluctuations are
irrelevant if the displacement due to thermal fluctuations is much smaller
than the displacement due to magnetic forces and torques. We can rewrite this
statement as $\langle\vec{r}^{2}\rangle\gg D_{T}t$ and
$\langle\vec{\theta}^{2}\rangle\gg D_{R}t$, where $\langle\vec{r}^{2}\rangle$
is displacement and $\langle\vec{\theta}^{2}\rangle$ is rotation of the
particle in characteristic time $t$, $D_{T}$ and $D_{R}$ are translation and
rotation diffusion constants respectively. If we take $1/\omega_{H}$ as the
characteristic time and $R$ as the characteristic distance. we can obtain that
the radius of the particle should fulfill
$R^{3}\gg(k_{B}T)/(4\pi\eta\omega_{H})$. The typical values of
$\sqrt[3]{(k_{B}T)/(4\pi\eta\omega_{H})}$ are in the range $[100nm,1\mu m]$.
So the particle size should be at least 10 times larger than the value of
$\sqrt[3]{(k_{B}T)/(4\pi\eta\omega_{H})}$ which would inherit that $R^{3}$ is
1000 times larger than $(k_{B}T)/(4\pi\eta\omega_{H})$.

##  III Theory

The theory by Melle et al. [4] shows that isotropic particles in a slow
rotating magnetic field form chains with length proportional to
$\frac{1}{\sqrt{\omega_{H}}}$. In this article we will show that if a particle
has an anisotropy of magnetic susceptibility, then the chain length in the
rotating magnetic field deviates from the trend $\frac{1}{\sqrt{\omega_{H}}}$
if the frequency $\omega_{H}$ is close to the critical frequency of a single
particle $\omega_{C,1}$, where

$\omega_{C,1}=\frac{\mu_{0}V_{m}\Delta\chi
H^{2}}{2\xi_{r}}=\frac{\mu_{0}V_{m}\Delta\chi H^{2}}{16\pi\eta R^{3}}$

and $H$ is the field strength of the external magnetic field.

We call $\omega_{C,1}$ the critical frequency of a single particle because an
isolated paramagnetic particle rotates synchronously with the field if
$\omega_{H}<\omega_{C,1}$ and asynchronously if $\omega_{H}>\omega_{C,1}$.[21]
If the particles magnetic anisotropy axis is a hard axis
($\chi_{\|}<\chi_{\perp}$), then the anisotropy axis lies perpendicularly to
the plane of the rotation of the magnetic field and the directions of
$\vec{m}$ and $\vec{H}$ coincide. The rotation of the particle is no longer
happening and the chain formation is the same as for isotropic particles where
chain length, as shown by Melle et al [4], is proportional to
$\frac{1}{\sqrt{\omega_{H}}}$.

For example, a spheroidal particle whose symmetry axis is 1% longer than the
other two axis and which is made of material with magnetic susceptibility
$\xi=7$ (in SI units), in water, in magnetic field with strength $H=1.6kA/m$,
will have the critical frequency $\omega_{C,1}=11rad/s$. The critical
frequency does not depend on the size of the particle but on the shape of the
magnetic core. The critical frequency is highly effected by the field strength
and can also be manipulated by viscosity of the fluid

Further in this article we will use the term critical frequency $\omega_{C,N}$
to emphasize difference between the current rotation frequency $\omega_{H}$
and the maximal possible field frequency for the given chain $\omega_{C,N}$.
If the external rotation frequency $\omega_{H}<\omega_{C,N}$, isolated chain
of length $N$ rotates with the field, but if $\omega_{H}>\omega_{C,N}$, the
isolated chain breaks. As a result, the theory of Melle [4] becomes
$N\propto\frac{1}{\sqrt{\omega_{C,N}}}$

###  III.1 Slow field

For particles with a magnetic anisotropy easy axis ($\chi_{\|}>\chi_{\perp}$),
the anisotropy axis is in the plane of the rotating field and the magnetic
moment of the particle depends on the angle between the external field
$\vec{H}$ and the anisotropy axis of the particle $\vec{n}$. In a slow
rotating magnetic field $\omega_{H}<\omega_{C,1}$ the particle follows the
field and the magnetic moment of the particle is constant:

$\absolutevalue{\vec{m}}=V_{m}\sqrt{\chi_{\perp}^{2}\vec{H}^{2}+2\chi_{\perp}\Delta\chi\left(\vec{n}\cdot\vec{H}\right)^{2}+\Delta\chi^{2}\left(\vec{n}\cdot\vec{H}\right)^{2}}\text{
,}$

where [20]

$\left(\vec{n}\cdot\vec{H}\right)^{2}=\frac{H^{2}}{2}\left(1+\sqrt{1-\left(\frac{\omega_{H}}{\omega_{C,1}}\right)^{2}}\right)\text{
.}$

From the theory of Melle et al. [4] chain will break if the angle between the
chain and the external magnetic field exceeds the critical value
$\frac{\pi}{4}$. In other words, the chain will break if hydrodynamic drag
$6\pi\eta R^{2}\omega_{C,N}N^{2}$ overcomes magnetic force
$\frac{3\mu_{0}}{4\pi(2R)^{4}}|\vec{m}|^{2}$. Putting all together, it gives a
formula connecting the chain length $N$ with the critical frequency
$\omega_{C,N}$ in a slow rotating field:

$N=\sqrt{\frac{V_{m}\omega_{C,1}}{\Delta\chi 16\pi
R^{3}\omega_{C,N}}}\sqrt{\chi_{\perp}^{2}+\chi_{\|}^{2}+(\chi_{\|}^{2}-\chi_{\perp}^{2})\sqrt{1-\frac{\omega_{C,N}^{2}}{\omega_{C,1}^{2}}}}$
(7)

It can be shown that the chain length $N$ deviates from the trend
$N\propto\frac{1}{\sqrt{\omega_{C,N}}}$ only for values of
$\frac{\omega_{C,N}}{\omega_{C,1}}$ close to unity. If the rotation frequency
is small ($\omega_{C,N}\ll\omega_{C,1}$), the equation (7) becomes

$N=\sqrt{\frac{V_{m}\chi_{\|}^{2}}{8\pi\Delta\chi
R^{3}}}\sqrt{\frac{\omega_{C,1}}{\omega_{C,N}}}=\sqrt{\frac{\mu_{0}V_{m}^{2}H^{2}\chi_{\|}^{2}}{16\pi\xi_{r}R^{3}}}\sqrt{\frac{1}{\omega_{C,N}}}\text{
,}$ (8)

which is proportionality $N\propto\frac{1}{\sqrt{\omega_{C,N}}}$. For
particles which are almost spherical, the critical rotation frequency
$\omega_{C,1}$ can become extremely small and frequencies $\omega_{C,N}$
corresponding to $N>1$ are much higher than the frequency $\omega_{C,1}$ and
the formula (7) is not suitable.

###  III.2 Fast field

Particles, which are almost isotropic, have low critical field frequency
$\omega_{C,1}$ and do not follow the field. This regime is called asynchronous
rotation because the angle between the external field $\vec{H}$ and the
particle direction $\vec{n}$ changes. Numerical simulations show that, if the
field frequency is large enough, particles synchronise with each other.

The force between two neighbouring synchronised particles would be:

$\vec{F}_{ij}=\dfrac{3\mu_{0}}{4\pi(2R)^{4}}\left[2(\vec{m}\cdot\hat{r})\vec{m}+\vec{m}^{2}\hat{r}-5(\vec{m}\cdot\hat{r})^{2}\hat{r}\right]\text{
,}$ (9)

where attractive radial force is

$F_{r}=\dfrac{3\mu_{0}}{4\pi(2R)^{4}}\left[\vec{m}^{2}-3(\vec{m}\cdot\hat{r})^{2}\right]$

and the force which makes particles to rotate along common center is

$\vec{F}_{\theta}=\dfrac{3\mu_{0}}{2\pi(2R)^{4}}(\vec{m}\cdot\hat{r})(\vec{m}\cdot\hat{\theta})\text{
,}$

where $\hat{r}$ is a unit vector in the direction to the other particle and
$\hat{\theta}$ is a unit vector perpendicular to $\hat{r}$ and in the
direction of the rotation. Using Eq. (1) and introducing angle $\theta$
between $\vec{H}$ and $\vec{r}$ and angle $\varphi$ between $\vec{n}$ and
$\hat{r}$ gives relations for magnetic moment $\vec{m}$:

$\hat{r}\cdot\vec{m}=V_{m}H\Big{(}\chi_{\perp}\cos\theta+\Delta\chi\cos\varphi\cos(\theta-\varphi)\Big{)}$
$\vec{m}^{2}=V_{m}^{2}H^{2}\Big{(}\chi_{\perp}^{2}+(\chi_{\|}^{2}-\chi_{\perp}^{2})\cos^{2}(\theta-\varphi)\Big{)}$
$\vec{m}\cdot\hat{\theta}=V_{m}H\Big{(}-\chi_{\perp}\sin\theta+\Delta\chi\sin\varphi\cos(\theta-\varphi)\Big{)}$

It will be assumed that in high frequency $\theta$ changes slowly and
$\varphi$ increases linearly with time. This provides a possibility to average
trigonometric functions over linearly increasing $\varphi$. The average inter-
particle force will read

$\langle
F_{r}\rangle=-\dfrac{3\mu_{0}V_{m}^{2}H^{2}}{32\pi(2R)^{4}}\left(2(\chi_{\perp}^{2}+\chi_{\|}^{2})+3(\chi_{\perp}+\chi_{\|})^{2}\cos(2\vartheta)\right)$
$\langle
F_{\theta}\rangle=-\dfrac{3\mu_{0}V_{m}^{2}H^{2}}{16\pi(2R)^{4}}(\chi_{\perp}+\chi_{\|})^{2}\sin(2\vartheta)$

For reasonable $\frac{\chi_{\|}}{\chi_{\perp}}$, the average force is close to
an isotropic situation, where $\chi_{\|}=\chi_{\perp}$.

The conclusion in high frequency range is similar to the conclusion for
isotropic particles [4]. In an anisotropic situation the angle $\theta$, at
which $\langle F_{r}\rangle$ becomes repulsive, is larger than the magic
angle, for any $\frac{\chi_{\|}}{\chi_{\perp}}$ value. Therefore, breakage of
the chain will happen when friction force will overcome maximal $\langle
F_{\theta}\rangle$. As a result, the maximal chain length in a high frequency
range is:

$N=\sqrt{\frac{\mu_{0}V_{m}^{2}H^{2}(\chi_{\|}+\chi_{\perp})^{2}}{64\pi
R^{3}\xi_{r}}}\sqrt{\frac{1}{\omega_{C,N}}}$ (10)

As can be seen, in high frequencies, the chain length is also proportional to
$\frac{1}{\sqrt{\omega_{C},N}}$. The proportionality constant in the fast
field is different from the proportionality constant found in the slow field.
The proportionality constant is higher for the slow field than for the fast
field if $\chi_{\|}>\chi_{\perp}$.

When the field frequency gets close to the critical frequency, the given
equation of chain length Eq. 10 will not work. To calculate a more precise
relation, the back and forth rotation of each particle should be considered in
the derivation, which gives field frequency dependent mean values for
trigonometric functions of $\varphi$. For example, mean value
$\langle\sin(2\varphi)\rangle=\frac{\omega_{H}-\sqrt{\omega_{H}^{2}-\omega_{C,1}^{2}}}{\omega_{C,1}}$
[21]. In addition, in asynchronous rotation of the particles, the assumption
that the angle $\theta$ between the chain and the field does not change does
not hold. The fact is that the angle $\theta$ changes do not allow to obtain
maximal chain length in the closed form in medium frequencies in the
asynchronous regime.

###  III.3 Long chain correction

The results obtained so far were based on the assumption that farther
particles do not influence the magnetic field of middle particles. This is not
true even for two particles. For longer chains, the influence of all particles
should be considered. Usually in theory an infinitely long chain is
considered. In an infinitely long straight chain the force and the torque on
each particle can be calculated, and they involve two angles ($\varphi$ and
$\theta$) in different trigonometric combinations. This does not allow to
obtain a solution in the closed form even for low field frequencies. In
simulations, it can be seen that the chain is not straight but S-shape. As a
result, further analysis of chain lengths is based on simulations results.

##  IV Results

The model described in section II was implemented in a C++ computer program,
to do numerical simulations. Translation motion Eq.(6) and rotation motion
Eq.(2) for each particle is being calculated by solving first order ordinary
differential equation system of $6N$ equations, where $N$ is the number of
particles. GNU Scientific Library [22], Runge-Kutta 4th order method with
adaptive time step were used to solve the differential equation system.

###  IV.1 Soft VS hard particles

By choosing parameter $r_{max}$ value in Eq.(5), steepness of the repulsion
force function can be regulated. The smaller $r_{max}$ value we choose, the
steeper the growth of the repulsive force, as can be seen in Fig. 2.

Figure 2: (a) Repulsive force Eq. (5) $F_{R}$ depending on particle distance
$r$, with different parameter $r_{max}$, normalized to maximal magnetic force
of two particles at distance $2R$. Steepness of the repulsive force function
$F_{R}$ (the softness of the particle) can be regulated by choosing parameter
$r_{max}$. (b) A two particle chain and a situation where the third particle
is at a distance where the force is repulsive (particle and forces in black)
and a situation where the third particle is at a distance where the force is
attractive to the chain (particle and forces in gray). Sizes of the force
vectors are drawn in a scale where $F_{i,M,max}=1$.

If there is a chain of two particles (a dimer) and the third particle comes
closer to the dimer, like illustrated in Fig.2, the magnetic force acting on
the third particle is repulsive if its position $x>2R$. If the third particle
position $x<2R$, then the magnetic force on the third particle is attractive
to the dimer and they form a trimer. The distance between the center of the
third particle and the line of centers of dimer particles, where the force
changes direction, is $\sqrt{5}R$. If the value of $r_{max}$ is larger than
$\sqrt{5}R$, then in a cluster of particles the probability for particles to
form a trimer is small and these three particles will form a chain structure.
The increase of $r_{max}$ not only increases the distance of the repulsive
force but also makes particles softer. The particles close to each other
amplify the field and the repulsive force at the distance 2R does not
compensate the attractive magnetic force completely, and the force equilibrium
distance is at a distance a bit less than 2R, which makes the force acting on
the third particle flip direction at a distance a bit less than $\sqrt{5}R$.
If larger $r_{max}$ value is used, the particles are softer than with small
$r_{max}$ value.

In the article [23] the order of the chain in the external field is being
regulated by salt concentration in the liquid. The higher the salt
concentration, the less disordered the chain gets. In our model similar
structural effects are observed by changing $r_{max}$ value. Higher $r_{max}$
value corresponds to higher salt concentration. Despite the fact that changing
salt concentration gives qualitatively similar results as changing $r_{max}$
value, the dynamics of such charged colloids could differ from the results
obtained in this article. Movement of the charged particle would lead to a
shifted ionic layer around the particle that creates electrostatic
interactions not included in the model discussed in this article.

###  IV.2 The critical length of a single chain

At the beginning, $N$ particles are placed in a straight chain along the
direction of a stationary field. At this phase particles obtain equilibrium
distance between each other. Then the field starts to rotate with a certain
frequency $\omega_{H}$. If the frequency is higher than the critical frequency
$\omega_{C,N}$, the chain breaks after some periods, but, if the frequency is
below the critical - the chain does not break. There is a very well observable
boundary between the frequency where the chain breaks in a few periods and the
frequency where the chain does not break in even a couple of hundred periods.

The angle between the line of the end particles of the chain and the external
field (Fig.3) is being monitored to detect when the chain breaks. If the angle
gets larger than $\frac{\pi}{2}$, then the chain is definitely broken. By
monitoring this angle even breaking of a two particle chain can be detected.

Figure 3:  S-shaped chain at the critical field rotation frequency
$\omega_{C,10}$. The angle between the line of the end particles of the chain
and the external field $\vec{H}$ are being monitored to detect breaking of the
chain. The field rotates clockwise. The easy magnetization axis is shown as a
black line in a particle.

If the field rotation frequency $\omega_{H}$ is slightly smaller than the
critical rotation frequency $\omega_{C,N}$, for a given chain, 4 structurally
different rotation regimes can be observed. If $\omega_{H}<\omega_{C,1}$, then
each particle rotates synchronously with the field (mode 1. Fig.4). If
$\omega_{H}$ is above $\omega_{C,1}$, then two scenarios can be observed. In
the first scenario, the end particles of the chain start to lag behind the
field and their rotation becomes asynchronous with the field, but the middle
particles rotate synchronously with the field (mode 2.a Fig. 4). In the second
scenario, the middle particle of the chain starts to lag the field (mode 2.b
Fig. 4). If the frequency is increased further, then more than two particles
rotate asynchronously with the field, and the rotation of the middle particles
is faster than the rotation of the end particles (mode 3. Fig. 4). If the
field frequency is increased even more, then the rotation frequency for all
particles becomes so small that it can be assumed that the particles
synchronize (mode 4. Fig. 4). Visual representation of all modes can be found
in the supplementary video.

Figure 4:  Chain rotation mode classification depending on the field rotation
frequency $\omega_{H}$ if the field rotation frequency is close to the
critical frequency $\omega_{C,N}$. 1. Synchronous mode (each particle rotates
synchronously with the field). 2.a Both end particles rotate asynchronously
with the field. 2.b Middle particle of the chain rotates asynchronously with
the field. 3. More than two particles rotate asynchronously with the field and
their rotation frequencies do not match. 4. All particles rotate slower than
the field and they rotate synchronously with each other. Examples of all modes
in video format can be found in supplementary materials. The positions of the
boundaries between the modes depend on the particle and field properties.

Fig.5 shows the critical chain length $N$ dependence on the field frequency
for particles with different magnetic anisotropy properties. For particles
with $\frac{\chi_{\|}}{\chi_{\perp}}=1.5$ and $\chi_{\|}=2$ there is a
transition between mode 1. and mode 2.a for a chain length of 13 and 12
particles (a 13 particles chain at its critical frequency $\omega_{C,13}=1.49$
rotates in mode 1 and a 12 particle chain at its critical frequency
$\omega_{C,12}=1.6$ rotates in mode 2.a). This mode change is not particularly
pronounced in Fig. 5 if the chain is long, but, for shorter chains, like in
example $\frac{\chi_{\|}}{\chi_{\perp}}=1.5$ and $\chi_{\|}=0.25$, the mode
transition is visible as a step in a plot between chain length of $5$ and $4$
particles.

For particles $\frac{\chi_{\|}}{\chi_{\perp}}=1.25$ and $\chi_{\|}=0.5$ there
is a transition from mode 1. to mode 2.b at chain length of 10 to 9, and an 8
particle chain rotates in mode 3. The transition between mode 2.a and mode 3.
is visible in example $\frac{\chi_{\|}}{\chi_{\perp}}=1.5$ and $\chi_{\|}=2$
for a chain length of 10 and 9 particles. Steepness of the plot changes as the
mode changes.

Figure 5:  Chain length $N$ and its critical field rotation frequency
$\frac{\omega_{C,N}}{\omega_{C,1}}$ for chains with different magnetic
anisotropy properties. Transitions between chain rotation modes are visible
like step or change of steepness of the graph.

For particles with low $\omega_{C,1}$, which corresponds to almost isotropic
particles, where $\frac{\omega_{C,N}}{\omega_{C,1}}>4$, mode 4. are observed.
The critical chain length in mode 4. can be well approximated as a function
proportional to $\frac{1}{\sqrt{\omega_{C.N}}}$. In Fig.6 linear fit is shown
for the critical chain length for chains where the critical chain length is in
mode 4. In Fig.6 for susceptibility values $\chi_{\|}/\chi_{\perp}=1.5$ and
$\chi_{\|}=2$ also linear fit in the region $\omega_{C,N}/\omega_{C,1}>1$
where the critical chain is in mode 1. can be observed. The linear fits have
different slopes in mode 1. and mode 4. for the same particles.

Figure 6: The critical chain length $N_{critical}$ as a function of
$\sqrt{\frac{\omega_{C,1}}{\omega_{C,N}}}$ for chains with different magnetic
anisotropy properties. Linear fit to $N\propto\frac{1}{\sqrt{\omega_{C,N}}}$
in mode 4. for small values (equation (10)) and mode 1. for large values
(equation (8)) is shown as dotted lines. (a) Fast rotation (large rotation
frequency $\omega_{C,N}$). (b) Slow rotation

Similarly the critical chain length is proportional to
$\frac{1}{\sqrt{\omega_{C.N}}}$ in mode 1., where particles have large
magnetic anisotropy and $\frac{\omega_{C,N}}{\omega_{C,1}}<1$. Linear fit in
mode 1. is shown in Fig. 6. These findings agree with the theory from section
III.

In modes 2. and 3. the critical chain length deviates from proportionality to
$\frac{1}{\sqrt{\omega_{C.N}}}$ as can be seen in Fig. 6 and Fig. 6 for
$\frac{\omega_{C,N}}{\omega_{C,1}}$ values close to unity. In modes 2. and 3.
there can be observed a switch from the slope of mode 1. seen in Fig. 6 and in
Fig. 6 to the right to the slope of mode 4. seen in Fig. 6. Since particles
rotate asynchronously between each other, the critical chain length is a bit
shorter than in neighbouring modes 1. and 4.

###  IV.3 Distribution of chain lengths in a small cluster

If at the beginning the particles are randomly distributed in a limited area
and a rotating field is being applied, then particles get organized in planes
where distribution of chain lengths can be observed. If $r_{max}$ is smaller
than $\sqrt{5}R$, particles create chains with the width larger than 1 and
other large structures. The region where $r_{max}<\sqrt{5}R$ is not studied.
If $r_{max}$ value is larger than $\sqrt{5}R$, in simulations $r_{max}=2.3R$,
particles create one particle thick chains. Chains rotate next to each other,
collide and exchange with particles. As a result, chains of different length
can be observed.

The length of each chain is being measured and chain length distribution is
obtained. If the distance between two particles is less than $2.1R$, then they
are being considered neighbors. If a particle has one neighbor, then it is an
end particle of a chain. If a particle has two neighbors, then it is one of
the middle particles of the chain. If a particle has more than two neighbors,
then the structure is not a 1 particle thin chain and all the neighboring
particles are excluded from the chain distribution results. This gives an
algorithm to obtain chain length distribution.

The simulations were done in an infinite box with no boundaries. In most
cluster simulations 61 or 91 particles were used. To verify the results, some
simulations were done also with 127 and 169 particles (numbers 61,91,127,169
are centred hexagonal numbers). Initially particle centers were randomly
distributed in a cylindrical volume of a height $2R$ and a diameter that
corresponds to a hexagonal crystal with a lattice constant $8R$, i.e. the
diameter of the initial cluster is $80R$ for 61 particles and the diameter of
the initial cluster is $96R$ for 91 particles. The tests showed that larger or
smaller initial density of particles increased the time to reach a stable
histogram.

The development of the distribution of chain lengths can be seen in Fig. 7 for
particle parameters $\frac{\chi_{\|}}{\chi_{\perp}}=1.5$ and $\chi_{\|}=0.25$
at field frequency $\frac{\omega_{H}}{\omega_{C,1}}=0.1$. The histogram of
chain lengths is obtained by counting chains over 400 subsequent periods and
100 frames in each period. As can be seen in the histogram, in the first 2000
periods a typical chain length of 12 particles developed.

Figure 7:  (a) The chain length $N$ distribution development over the time in
the cluster of 91 particles. (b) The chain length $N$ distribution dependence
on the number of particles in the cluster. Simulation gives typical chain
length $N_{ensemble}=12\pm 0.5$. Each histogram contains the sum of chain
length distribution histograms time interval $t$ given in the number of
periods from the beginning of the calculation and 100 frames per period.
$\frac{\chi_{\|}}{\chi_{\perp}}=1.5$ and $\chi_{\|}=0.25$ at field frequency
$\frac{\omega_{H}}{\omega_{C,1}}=0.1$.

To simulate an ensemble of particles an infinite number of particles or
periodic boundary conditions should be used in simulations. Our simulations
show that 2-3 times more particles than the critical chain length is enough to
get the distribution of chain lengths in the cluster. If more particles are
being used or the initial layout of particles, or particles density are
changed then only the time for reaching the stable configuration changes, but
the final distribution of chains lengths is not significantly changed. The
distribution of chain lengths in a cluster dependence on the initial number of
particles is shown in Fig. 7. The histogram of 61 particles is similar to the
histograms with higher particle count could be because 61 particles is enough
to create all possible chain length configurations in different time steps and
averaging over time is similar to averaging over space. 37 particles create
three rotating chains which poorly collide with each other.

In mode 1. there is a certain chain length that dominates in the distribution
of chains lengths, but in all other modes the width of the dominating peak is
larger. The error bars are being calculated as a half-width at half of
maximum.

The ratio between the critical chain length of a single chain and the typical
chain length in a cluster of chains has been calculated. Numerical simulations
show that the critical and the typical chain length ratio is constant for
certain type of particle as can be seen in Fig. 8 graphs A1, A2, B1 and B2.
Ratio mean values coincide within its standard deviations in cases A1 and B1
($1.65\pm 0.06$ for A1 and $1.68\pm 0.07$ for B1). For larger $\chi_{\|}$ the
ratio between the critical and the typical chain length decrease, as can be
seen in Fig. 8 graph A2 (mean value $1.57\pm 0.07$) and B2 (mean value
$1.48\pm 0.15$). Graph B2 includes particle motion not only in mode 1., but
also modes 2.b and 3. In this case the ratio does not change much as the modes
change.

In figures A3 and B3 Fig. 8 typical chain length in a cluster becomes close to
the critical chain length and the ratio becomes one. The ratio becomes one if
the critical chain length is in mode 2.a. The mode 2.a can be observed in A3
and B3 Fig. 8 for frequencies up to $\frac{\omega_{H}}{\omega_{C,1}}=1.79$ in
graph A3 and up to $\frac{\omega_{H}}{\omega_{C,1}}=1.72$ in graph B3. If the
frequency is increased, the mode 2.a gradually changes to mode 3., and the
ratio becomes larger. The explanation for the fact that in mode 2.a typical
chain length becomes equal to typical chain length in a cluster could be
following: In a cluster, when two identical chains meet with each other, for a
moment they form one twice as long chain, in mode 2.a the chain will split in
half because only the end particles of these two chains were disorientated.
Unlike in all other modes where all the particles have similar disordered
state and have more equal probability become the breaking point.

Rotation frequencies in Fig.8 were chosen to coincide with the critical
frequencies which for most of the parameters can be found in Fig.5. There are
two reasons why error bars in Fig.8 are larger for larger frequencies.
Firstly, the chain distribution is wider in the asynchronous regime. Secondly,
the chains become shorter for large frequencies. The width of the distribution
can not be smaller than 1 which gives error 0.5 for $N_{critical}=2$.

In addition, simulations using isotropic particles were done. Since
$\omega_{C,1}=0$ for an isotropic particle, field frequency $\omega_{H}$ was
characterized by Mason number:

$Mn^{\prime}=\dfrac{\omega_{H}\zeta_{r}}{\mu_{0}V_{m}H^{2}}=\frac{\omega_{H}8\pi\eta
R^{3}}{\mu_{0}V_{m}H^{2}}$ (11)

We use apostrophe here ($Mn^{\prime}$) because different Mason number
definitions have been seen in literature. Our definition disagrees these
definitions: $Mn^{\prime}=4Mn$ if $Mn$ from [4, 24] is used;
$Mn^{\prime}=Mn/8$ if $Mn$ from [25] is used.

Summary of the results of a cluster of isotropic particles can be seen in
graph C Fig. 8 for the critical chain lengths from 8 (for larger
$Mn^{\prime}$) till 18 ( for smaller $Mn^{\prime}$). As can be seen, the ratio
$N_{critical}/N_{ensemble}=1.5$ also for isotropic particles. The width of the
chain length distribution (error-bars in graph C Fig. 8) for isotropic
particles is similar to that obtained in synchronous regime of anisotropic
particles as can be seen in low frequencies in graphs A1, B1 and A2 Fig, 8.

Figure 8:  Ratio of the critical length of one chain and the typical length of
particles $\frac{N_{critical}}{N_{ensemble}}$ for particles with different
magnetic properties. Weighted mean value of ratio is illustrated as a solid
line and weighted standard deviation is illustrated with a dashed line. Error
bars represent the width of the distribution of $N$ in a cluster of particles.
In graph C, shown ratio $\frac{N_{critical}}{N_{ensemble}}$ for isotropic
particles. In graph C, Mason number (11) is used as x-axis. Ratios in graphs
A1, A2, B1, B2 and C are constant within standard deviations. Ratios in graphs
A1 and B1 coincide within standard deviation. The larger magnetic
susceptibility the smaller the ratio is.

##  V Conclusions and discussion

Numerical simulations show that the critical chain length dependence on field
frequency for paramagnetic particles with magnetic anisotropy obeys trend
observed for isotropic particles if the field frequency is not close to the
critical frequency of individual particle $\omega_{C,1}$. For both particles,
magnetically isotropic and anisotropic, the critical chain length in rotating
field can be approximated as a function inverse proportional to square root of
rotation frequency $N\propto\frac{1}{\sqrt{\omega_{H}}}$. The deviation from
this trend for paramagnetic particles with magnetic anisotropy is the highest
if rotation frequency is slightly above the critical frequency of individual
particle $\omega_{C,1}$. There is a small difference in the proportionality
factor in $N\propto\frac{1}{\sqrt{\omega_{H}}}$ above and below the critical
frequency $\omega_{C,1}$. In conclusion, the observation of the critical chain
length does not allow us to identify that the particles have magnetic
anisotropy if the field frequency is not close to the critical frequency of an
individual particle $\omega_{C,1}$.

Simulations show that the typical length in a cluster of chains is smaller
than the maximal (the critical) length of the single chain, for the same
frequency. Typical chain length is $1.5$ times smaller than the critical chain
length. The exception is the situation where only two end particles of the
single chain rotate asynchronously with the field, when the field frequency is
close to the critical. In this situation, typical chain length is almost equal
to the critical chain length. This can be observed if the rotation frequency
is slightly higher than the critical frequency of an individual particle
$\omega_{C,1}$ and particles have high magnetic susceptibility. As a result,
the fact that the chain length in a cluster is smaller than the critical chain
length should be considered, when comparing a theory with an experiment.

In low frequency regime, where field rotation frequency $\omega_{H}$ is
smaller than the critical frequency of an individual particle $\omega_{C,1}$,
and particles rotate synchronously with the field, distribution of chain
lengths in a cluster is narrower than in high frequency range, where
$\omega_{H}>\omega_{C,1}$. All particles, which have anisotropy, will rotate
slowly and the chain length distribution will narrow down. Therefore, to
create large amounts of paramagnetic chains with equal length it is better to
use paramagnetic particles with magnetic anisotropy then magnetically almost
isotropic particles. On contrary, fully isotropic particles will have similar
width of distribution of chain lengths as highly anisotropic particles.

In this article a small cluster in an infinite medium is used. This allowed
the particles to settle in the optimal configuration where attraction of the
particles is compensated by the chain breakage. In contrast, if the large
ensemble of particles would have been used, the obtained results would depend
on the concentration of the particles. For higher particle concentration, as
used in magnetorheological fluids, there would not be enough space for the
critical chain length to exist and the typical chain length would be smaller
than obtained in simulation in a small cluster. In an ensemble also long range
interactions play a significant role, which are negligible in a small cluster.
The role of particle concentration, long range interactions and temperature in
the chain formation and stability in an ensemble of anisotropic particles are
left open for further research.

##  VI Acknowledgments

This work was supported by PostDoc Latvia
[Project no. 1.1.1.2/VIAA/1/16/060].

## References

  * [1] B. Yigit, Y. Alapan, M. Sitti, Cohesive self-organization of mobile microrobotic swarms, Soft Matter 16 (2020) 1996–2004.  doi:10.1039/C9SM01284B.
  * [2] H. Massana-Cid, F. Meng, D. Matsunaga, R. Golestanian, P. Tierno, Tunable self-healing of magnetically propelling colloidal carpets, Nature Communications 10 (1) (2019) 2444.  doi:10.1038/s41467-019-10255-4.
  * [3] K. Han, G. Kokot, S. Das, R. G. Winkler, G. Gompper, A. Snezhko, Reconfigurable structure and tunable transport in synchronized active spinner materials, Science Advances 6 (12) (2020).  doi:10.1126/sciadv.aaz8535.
  * [4] S. Melle, J. E. Martin, Chain model of a magnetorheological suspension in a rotating field, The Journal of Chemical Physics 118 (21) (2003) 9875–9881.  doi:10.1063/1.1570817.
  * [5] Y. Gao, M. A. Hulsen, T. G. Kang, J. M. J. den Toonder, Numerical and experimental study of a rotating magnetic particle chain in a viscous fluid, Phys. Rev. E 86 (2012) 041503.  doi:10.1103/PhysRevE.86.041503.
  * [6] A. Vázquez-Quesada, T. Franke, M. Ellero, Theory and simulation of the dynamics, deformation, and breakup of a chain of superparamagnetic beads under a rotating magnetic field, Physics of Fluids 29 (3) (2017) 032006.  doi:10.1063/1.4978630.
  * [7] M. Devi, P. P. Dutta, D. Mohanta, Analytical calculation of chain length in ferrofluids, Bulletin of Materials Science 38 (1) (2015) 221–226.  doi:10.1007/s12034-014-0812-9.
  * [8] W. Lo, C. Lu, C. Lin, C. Chen, Trajectory of a non-magnetic particle transported by a rotating magnetic particle chain, IEEE Transactions on Magnetics 55 (2) (2019) 1–4.  doi:10.1109/TMAG.2018.2865854.
  * [9] A. Cēbers, I. Javaitis, Dynamics of a flexible magnetic chain in a rotating magnetic field, Phys. Rev. E 69 (2004) 021404.  doi:10.1103/PhysRevE.69.021404.
  * [10] S. L. Biswal, A. P. Gast, Rotational dynamics of semiflexible paramagnetic particle chains, Physical review. E, Statistical, nonlinear, and soft matter physics 69 (4 Pt 1) (2004) 041406.  doi:10.1103/physreve.69.041406.
  * [11] P. Tierno, Recent advances in anisotropic magnetic colloids: realization, assembly and applications, Phys. Chem. Chem. Phys. 16 (2014) 23515–23528.  doi:10.1039/C4CP03099K.
  * [12] A. Stikuts, R. Perzynski, A. Cēbers, Spontaneous order in ensembles of rotating magnetic droplets, Journal of Magnetism and Magnetic Materials 500 (2020) 166304.  doi:10.1016/j.jmmm.2019.166304.
  * [13] C.-Y. Chen, H.-C. Hsueh, S.-Y. Wang, Y.-H. Li, Self-assembly and novel planetary motion of ferrofluid drops in a rotational magnetic field, Microfluidics and Nanofluidics 18 (5) (2015) 795–806.  doi:10.1007/s10404-014-1472-1.
  * [14] F. Martinez-Pedrero, A. Ortiz-Ambriz, I. Pagonabarraga, P. Tierno, Colloidal microworms propelling via a cooperative hydrodynamic conveyor belt, Phys. Rev. Lett. 115 (2015) 138301.  doi:10.1103/PhysRevLett.115.138301.
  * [15] A. Aharoni, Introduction to the Theory of Ferromagnetism, International Series of Monographs on Physics, Clarendon Press, 2000.
  * [16] A. Mukhtar, K. Wu, X. Cao, liyuan Gu, Magnetic nanowires in biomedical applications, Nanotechnology 31 (43) (2020) 433001.  doi:10.1088/1361-6528/aba1ba.
  * [17] E. Yammine, E. Souaid, S. Youssef, M. Abboud, S. Mornet, M. Nakhl, E. Duguet, Particles with magnetic patches: Synthesis, morphology control, and assembly, Particle & Particle Systems Characterization 37 (7) (2020) 2000111.  doi:10.1002/ppsc.202000111.
  * [18] E. C. Stoner, Xcvii. the demagnetizing factors for ellipsoids, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 36 (263) (1945) 803–821.  doi:10.1080/14786444508521510.
  * [19] J. A. Osborn, Demagnetizing factors of the general ellipsoid, Phys. Rev. 67 (1945) 351–357.  doi:10.1103/PhysRev.67.351.
  * [20] J. Cīmurs, A. Brasovs, K. Ērglis, Stability analysis of a paramagnetic spheroid in a precessing field, Journal of Magnetism and Magnetic Materials 491 (2019) 165630.  doi:10.1016/j.jmmm.2019.165630.
  * [21] B. Frka-Petesic, K. Erglis, J. Berret, A. Cebers, V. Dupuis, J. Fresnais, O. Sandre, R. Perzynski, Dynamics of paramagnetic nanostructured rods under rotating field, Journal of Magnetism and Magnetic Materials 323 (10) (2011) 1309 – 1313, proceedings of 12th International Conference on Magnetic Fluid.  doi:10.1016/j.jmmm.2010.11.036.
  * [22] B. Gough, GNU Scientific Library Reference Manual - Third Edition, 3rd Edition, Network Theory Ltd., 2009.
URL https://www.gnu.org/software/gsl/

  * [23] Z. He, H. D. Robinson, Assembly and disorder dissipation in superparamagnetic nanoparticle chains in a rotating magnetic field, arXiv (2020).  arXiv:arXiv:2003.03008.
  * [24] D. Klingenberg, J. Ulicny, M. Golden, Mason numbers for magnetorheology, Journal of Rheology 51 (5) (2007) 883–893.  doi:10.1122/1.2764089.
  * [25] S. G. Sherman, A. C. Becnel, N. M. Wereley, Relating mason number to bingham number in magnetorheological fluids, Journal of Magnetism and Magnetic Materials 380 (2015) 98–104.  doi:10.1016/j.jmmm.2014.11.010.