# Magnetic stiffness calculation for the corresponding force between two
current-carrying circular filaments arbitrarily oriented in the space

Kirill Poletkin<EMAIL_ADDRESS>Slobodan Babic<EMAIL_ADDRESS>Innopolis University, 1, Universitetskaya Str., Innopolis, 420500, Russia  The
Institute of Microstructure Technology, Karlsruhe Institute of Technology,
Hermann-von-Helmholtz-Platz 1, 76344 Eggenstein-Leopoldshafen, Germany  53
Berlioz 101, H3E 1N2, Montréal, Québec, Canada

###### Abstract

In this article, sets of analytical formulas for calculation of nine
components of magnetic stiffness of corresponding force arising between two
current-carrying circular filaments arbitrarily oriented in the space are
derived by using Babic’s method and the method of mutual inductance
(Kalantarov-Zeitlin’s method). Formulas are presented through integral
expressions, whose kernel function is expressed in terms of the elliptic
integrals of the first and second kinds. Also, we obtained an additional set
of expressions for calculation of components of magnetic stiffness by means of
differentiation of Grover’s formula of the mutual inductance between two
circular filaments with respect to appropriate coordinates. The derived sets
of formulas were mutually validated and results of calculation of components
of magnetic stiffness agree well to each other.

###### keywords:

Magnetic stiffness, Circular filaments, Mutual inductance, Magnetic force,
Magnetic torque

††journal: Journal of Magnetism and Magnetic Materials

##  1 Introduction

Analytical and semi-analytical methods in the calculation of self- and mutual-
inductances of conducting elements of electrical circuits and magnetic force
interactions between these elements have become powerful mathematical
instruments in development of power transfer, wireless communication, sensing
and actuation and have been applied in a broad fields of science, including
electrical and electronic engineering, medicine, physics, nuclear magnetic
resonance, mechatronics and robotics, to designate the most prominent.
Although, a number of efficient numerical methodes implemented in the
commercially developed software are available, analytical and semi-analytical
methods allow to obtaining the result of calculations in the form of a final
formula with a finite number of input parameters, which when applicable may
significantly reduce computation effort. Providing the direct access to a
calculational formula for a user in such methods facilitate mathematical
analysis of obtained results of calculation and opens an opportunity for their
further optimization.

Analytical methods applied to the calculation of mutual inductance between two
circular filaments and arising magnetic force, magnetic torque and
corresponding magnetic stiffness when such the filament system carries
electric currents is a prime example. These methods have proved their
efficiency and have been successfully employed in an increasing number of
applications, including electromagnetic levitation [1, 2], superconducting
levitation [3], calculation of mutual inductance between thick coils [4],
magnetic force and torque calculation between circular coils [5, 6, 7],
wireless power transfer [8, 9, 10], electromagnetic actuation [11, 12, 13],
micro-machined contactless inductive suspensions [14, 15, 16] and hybrid
contactless suspensions [17, 18, 19, 20], biomedical applications [21, 22],
topology optimization [23], nuclear magnetic resonance [24, 25], indoor
positioning systems [26], navigation sensors [27], non-contact gap measurement
sensors [28], wireless power transfer systems [29, 30], magneto-inductive
wireless communications [31] and others.

In the present article, the set of formulas for calculation of nine components
of magnetic stiffness of corresponding force arising between two current-
carrying circular filaments arbitrarily oriented in the space are derived by
using two methods, namely, Babic’s method and the method of mutual inductance.
In the first one, the components of the magnetic field at an arbitrary point
of the secondary circular filament generated by the primary coil carrying
eletric current are calculated and then after taking the first derivatives of
these components with respect to the appropriate coordinates, the set of
analytical formulas for calculation of magnetic stiffness appeared in the
integral form whose kernel function is expressed in terms of the elliptic
integrals of the first and second kinds is derived. In the second method, the
calculation of components of magnetic stiffness of the corresponding magnetic
force is performed by means of finding the second derivatives of the function
of mutual inductance between two circular filaments recieved by using
Kalantarov-Zeitlin’s method [32] with respect to the appropriated coordinates.

The article is organized in the following way. In section 2 of the paper,
Babic’s method and its basic expressions are introduced and the set of
analytical formulas for calculation of nine components of magnetic stiffness
is derived based on this method. In section 3, the method of mutual inductance
is presented. The section includes the preliminary discussion, where the set
of coordinate frames necessary for determining the position of the secondary
coil with respect to the primary one by using Grover’s angles is given. Also,
the relationship between constants of the inclined plane equation employing to
define the angular misaliment of the secondary circular filament with respect
to the primary one in Babic’s approach and Grover’s angles is shown. In
section 4, sets of analytical formulas for calculation of nine components of
magnetic stiffness recieved by means of the introduced two approaches are
mutually verified via a number of designed examples. In section 5, conclusions
about obtained results are discussed. In the appendix, in addition to Babic’s
method and the method of mutual inductance, Grover’s method is introduced and
a set of analytical formulas for calculation of magnetic stiffness based on
Grover’s method is obtained.

##  2 Babic’s method (BM)

Figure 1: General scheme of arbitrarily positioning two current-carrying
circular filaments with respect to each other.

Let us take into consideration two current-carrying circular filaments as
showed in Fig. 1, where the center of the larger loop (primary coil) of the
radius $R_{p}$ is placed at the plane $XOY$ whose center is $O$ $(0,0,0)$. The
smaller circular loop (secondary coil) of the radius $R_{s}$ is placed in an
inclined plane whose general equation is,

${\lambda\equiv ax+by+cz+d=0,}$ (1)

where $a$, $b$, $c$ and $d$ are the components of the normal $\vec{N}$ on the
inclined plane in the center of the secondary circular segment $C$
$(x_{c},y_{c},z_{c})$.

###  2.1 Basic expressions

The segments are with the currents $I_{p}$ and $I_{s}$, respectively. For
circular filaments (see Fig. 1) we define, [33, 34, 35, 36]:

1)

Since, the primary circular filament is placed in the plane $XOY$ $(Z=0)$ with
the center at $O$ $(0,0,0)$. Hence, an arbitrary point $P$
$(x_{p},y_{p},z_{p})$ of this filament has the following parametric
coordinates (see Fig.1):

$x_{p}=R_{p}\cos\phi,\;y_{p}=R_{p}\sin\phi,\;z_{p}=0,\;\phi\in[0,2\pi].$ (2)

2)

The differential of the primary circular filament is given by

$d\vec{l}_{p}=R_{p}\\{-\sin\phi,\cos\phi,0\\}d\phi,\;\phi\in[0,2\pi].$ (3)

3)

The secondary circular filament of radius $R_{s}$ is placed in the inclined
plane (1) with the center at $C$ $(x_{c},y_{c},z_{c})$. The unit vector
$\vec{N}$ (the unit vector of the $z$-axis) at the point $C$, which is the
center of the secondary circular filament, laying in the plane $\lambda$ is
defined by

$\vec{N}=\left\\{\frac{a}{L},\frac{b}{L},\frac{c}{L}\right\\},\;L=\sqrt{a^{2}+b^{2}+c^{2}}.$
(4)

4)

The unit vector between two points C and S they are placed in the plane (1) is

$\vec{u}=\left\\{u_{x},u_{y},u_{z}\right\\}=\left\\{-\frac{ab}{Ll},\frac{l}{L},-\frac{cb}{Ll}\right\\},\;l=\sqrt{a^{2}+c^{2}}.$
(5)

5)

We define the unite vector $\vec{v}$ lying in the plane (1) and mutually
perpendicular on the unit vectors $\vec{N}$ and $\vec{u}$ as the cross-product
as follows

$\vec{v}=\vec{N}\times\vec{u}=\left\\{v_{x},v_{y},v_{z}\right\\}=\left\\{-\frac{c}{l},0,\frac{a}{l}\right\\}.$
(6)

6)

An arbitrary point $S$ $(x_{s},y_{s},z_{s})$ of the secondary circular
filament has parametric coordinates

$\begin{array}[]{l}x_{s}=x_{c}+R_{s}u_{x}\cos\vartheta+R_{s}v_{x}\sin\vartheta;\\\
y_{s}=y_{c}+R_{s}u_{y}\cos\vartheta+R_{s}v_{y}\sin\vartheta;\\\
z_{s}=z_{c}+R_{s}u_{z}\cos\vartheta+R_{s}v_{z}\sin\vartheta,\;\vartheta\in[0,2\pi].\end{array}$
(7)

This is well-known parametric equation of circle in 3D space. The filamentary
circular filaments are the part of this circle.

7)

The differential element of the secondary circular filament is given by,

$d\vec{l}_{s}=R_{s}\\{l_{xs},l_{ys},l_{zs}\\}d\vartheta,\;\;\vartheta\in[0,2\pi],$
(8)

where

$\begin{array}[]{l}l_{xs}=-u_{x}\sin\vartheta+v_{x}\cos\vartheta;\\\
l_{ys}=-u_{y}\sin\vartheta+v_{y}\cos\vartheta;\\\
l_{zs}=-u_{z}\sin\vartheta+v_{z}\cos\vartheta.\end{array}$ (9)

###  2.2 Stiffness calculation

To calculate the stiffness between two inclined circular loops as the prime
interest of this article we use the analytical formulas for calculating the
magnetic field produced by the primary current carrying with the current
$I_{p}$ at the arbitrary point $S$ $(x_{s},y_{s},z_{s})$ of the secondary
inclined current carrying loop with the current $I_{s}$ [36]. Hence, the
components of the field can be calculated as follows,

$B_{x}=\frac{\mu_{0}I_{p}z_{s}k}{16\pi
p^{2}\sqrt{R_{p}p}\left(1-k^{2}\right)}I_{x};$ (10)
$B_{y}=\frac{\mu_{0}I_{p}z_{s}k}{16\pi
p^{2}\sqrt{R_{p}p}\left(1-k^{2}\right)}I_{y};$ (11)
$B_{z}=-\frac{\mu_{0}I_{p}k}{16\pi p\sqrt{R_{p}p}\left(1-k^{2}\right)}I_{z},$
(12)

where $p=\sqrt{x_{s}^{2}+y_{s}^{2}}$, $I_{x}=x_{s}A$, $I_{y}=y_{s}A$ and
$I_{z}=D$ with

$\begin{array}[]{l}{\displaystyle
k^{2}=\frac{4R_{p}p}{(R_{p}+p)^{2}+z_{s}^{2}};}\\\
A=-2\left[(k^{2}-2)E(k)+(2-2k^{2})K(k)\right];\\\
D=-2\left[\left(k^{2}(R_{p}+p)-2p\right)E(k)+p\left(2-2k^{2}\right)K(k)\right].\end{array}$

In given expressions $K(k)$ and $E(k)$ are the complete elliptic integrals of
the first and the second kind, respectively.

Let us find the first derivatives of the components of the field with respect
to coordinates $x_{s}$, $y_{s}$ and $z_{s}$, we can write

$\begin{array}[]{l}{\displaystyle\frac{\partial B_{x}}{\partial
g}=\frac{\mu_{0}I_{p}}{16\pi\sqrt{R_{p}}}\frac{p^{-\frac{3}{2}}}{\left(1-k^{2}\right)}\frac{z_{s}}{p}\left\\{\left[-\frac{5}{2p}k\frac{\partial
p}{\partial g}+\frac{1+k^{2}}{1-k^{2}}\frac{\partial k}{\partial
g}\right]I_{x}+k\frac{\partial I_{x}}{\partial g}\right\\};}\\\
{\displaystyle\frac{\partial B_{y}}{\partial
g}=\frac{\mu_{0}I_{p}}{16\pi\sqrt{R_{p}}}\frac{p^{-\frac{3}{2}}}{\left(1-k^{2}\right)}\frac{z_{s}}{p}\left\\{\left[-\frac{5}{2p}k\frac{\partial
p}{\partial g}+\frac{1+k^{2}}{1-k^{2}}\frac{\partial k}{\partial
g}\right]I_{y}+k\frac{\partial I_{y}}{\partial g}\right\\};}\\\
{\displaystyle\frac{\partial B_{z}}{\partial
g}=\frac{\mu_{0}I_{p}}{16\pi\sqrt{R_{p}}}\frac{p^{-\frac{3}{2}}}{\left(1-k^{2}\right)}\left\\{\left[-\frac{3}{2p}k\frac{\partial
p}{\partial g}+\frac{1+k^{2}}{1-k^{2}}\frac{\partial k}{\partial
g}\right]I_{z}+k\frac{\partial I_{z}}{\partial g}\right\\};}\\\
g=x_{s},y_{s},z_{s},\end{array}$ (13)

where

$\begin{array}[]{l}{\displaystyle\frac{\partial k}{\partial
x_{s}}=\frac{x_{s}k^{3}}{8R_{p}p^{3}}\left[R_{p}^{2}+z_{s}^{2}-p^{2}\right],\;\frac{\partial
k}{\partial
y_{s}}=\frac{y_{s}k^{3}}{8R_{p}p^{3}}\left[R_{p}^{2}+z_{s}^{2}-p^{2}\right],\;\frac{\partial
k}{\partial z_{s}}=-\frac{z_{s}k^{3}}{4R_{p}p};}\\\
{\displaystyle\frac{\partial p}{\partial
x_{s}}=\frac{x_{s}}{p},\;\frac{\partial p}{\partial
y_{s}}=\frac{y_{s}}{p},\;\frac{\partial p}{\partial z_{s}}=0;}\end{array}$
$\begin{array}[]{l}{\displaystyle\frac{\partial I_{x}}{\partial
x_{s}}=A+Cx_{s}\frac{\partial k}{\partial x_{s}},\;\frac{\partial
I_{x}}{\partial y_{s}}=Cx_{s}\frac{\partial k}{\partial
y_{s}},\;\frac{\partial I_{x}}{\partial z_{s}}=Cx_{s}\frac{\partial
k}{\partial z_{s}};}\\\ {\displaystyle\frac{\partial I_{y}}{\partial
x_{s}}=Cy_{s}\frac{\partial k}{\partial x_{s}},\;\frac{\partial
I_{y}}{\partial y_{s}}=A+Cy_{s}\frac{\partial k}{\partial
y_{s}},\;\frac{\partial I_{y}}{\partial z_{s}}=Cy_{s}\frac{\partial
k}{\partial z_{s}};}\\\ {\displaystyle\frac{\partial I_{z}}{\partial
x_{s}}=D\frac{x_{s}}{p}=T,\;\frac{\partial I_{z}}{\partial
y_{s}}=D\frac{y_{s}}{p}+T,\;\frac{\partial I_{z}}{\partial
z_{s}}=T;}\end{array}$ $\begin{array}[]{l}{\displaystyle
C=-6x_{s}k\left[E(k)-K(k)\right];}\\\ {\displaystyle T=-2k\frac{\partial
k}{\partial x_{s}}\left[3(R_{p}+p)E(k)+(R_{p}+3p)K(k)\right].}\end{array}$

Accounting for the fact that

$\frac{\partial l_{xs}}{\partial g}=\frac{\partial l_{ys}}{\partial
g}=\frac{\partial l_{zs}}{\partial g}=0,\;g=x_{s},y_{s},z_{s},$ (14)

magnetic stiffness calculation for the corresponding force between two
current-carrying circular filaments arbitrarily oriented in the space is given
by the following formulas [36]:

$S_{xx}=-\frac{\partial F_{x}}{\partial
x_{s}}=-I_{s}R_{s}\int_{0}^{2\pi}\left[l_{ys}\frac{\partial B_{z}}{\partial
x_{s}}-l_{zs}\frac{\partial B_{y}}{\partial x_{s}}\right]d\vartheta;$ (15)
$S_{xy}=-\frac{\partial F_{x}}{\partial
y_{s}}=-I_{s}R_{s}\int_{0}^{2\pi}\left[l_{ys}\frac{\partial B_{z}}{\partial
y_{s}}-l_{zs}\frac{\partial B_{y}}{\partial y_{s}}\right]d\vartheta;$ (16)
$S_{xz}=-\frac{\partial F_{x}}{\partial
z_{s}}=-I_{s}R_{s}\int_{0}^{2\pi}\left[l_{ys}\frac{\partial B_{z}}{\partial
z_{s}}-l_{zs}\frac{\partial B_{y}}{\partial z_{s}}\right]d\vartheta;$ (17)
$S_{yx}=-\frac{\partial F_{y}}{\partial
x_{s}}=I_{s}R_{s}\int_{0}^{2\pi}\left[l_{xs}\frac{\partial B_{z}}{\partial
x_{s}}-l_{zs}\frac{\partial B_{x}}{\partial x_{s}}\right]d\vartheta;$ (18)
$S_{yy}=-\frac{\partial F_{y}}{\partial
y_{s}}=I_{s}R_{s}\int_{0}^{2\pi}\left[l_{xs}\frac{\partial B_{z}}{\partial
y_{s}}-l_{zs}\frac{\partial B_{x}}{\partial y_{s}}\right]d\vartheta;$ (19)
$S_{yz}=-\frac{\partial F_{y}}{\partial
z_{s}}=I_{s}R_{s}\int_{0}^{2\pi}\left[l_{xs}\frac{\partial B_{z}}{\partial
z_{s}}-l_{ys}\frac{\partial B_{x}}{\partial z_{s}}\right]d\vartheta;$ (20)
$S_{zx}=-\frac{\partial F_{z}}{\partial
x_{s}}=-I_{s}R_{s}\int_{0}^{2\pi}\left[l_{xs}\frac{\partial B_{y}}{\partial
x_{s}}-l_{ys}\frac{\partial B_{x}}{\partial x_{s}}\right]d\vartheta;$ (21)
$S_{zy}=-\frac{\partial F_{z}}{\partial
y_{s}}=-I_{s}R_{s}\int_{0}^{2\pi}\left[l_{xs}\frac{\partial B_{y}}{\partial
y_{s}}-l_{ys}\frac{\partial B_{x}}{\partial y_{s}}\right]d\vartheta;$ (22)
$S_{zz}=-\frac{\partial F_{z}}{\partial
z_{s}}=-I_{s}R_{s}\int_{0}^{2\pi}\left[l_{xs}\frac{\partial B_{y}}{\partial
z_{s}}-l_{ys}\frac{\partial B_{x}}{\partial z_{s}}\right]d\vartheta.$ (23)

Thus, all magnetic stiffness components (15)-(23) between two inclined
current-carrying loops are given in the simple integral form, over the
complete elliptic integrals of the first and the second kind. These
expressions can be used for resolving the singular cases. It is necessary to
use expressions (15)-(23) with the following conditions $(a=c=0,l=0,L=|b|)$
for the unit vectors $\vec{u}=\left\\{-1,0,0\right\\}$ and
$\vec{v}=\left\\{0,0,1\right\\}$. The loops are perpendicular mutually.

It is clear that $S_{xy}=S_{yx}$, $S_{xz}=S_{zx}$, $S_{yz}=S_{zy}$, so that
the calculation can be simplified by finding only sixth stiffness: $S_{xx}$,
$S_{yy}$, $S_{zz}$, $S_{xy}$, $S_{xz}$, and $S_{yz}$. Doing further
investigation one can find that [37],

$S_{xx}+S_{yy}+S_{zz}=0,$ (24)

so that the problem of the stiffness calculation can be limited to find only
five components, for instance, $S_{xx}$, $S_{yy}$, $S_{xy}$, $S_{xz}$, and
$S_{yz}$.

##  3 Mutual Inductance Method (MIM)

In this section, the mutual inductance method as an alternative to Babic’s
method discussed above is presented. The essence of the method is that the
calculation of the stiffness of the corresponding magnetic force is performed
by means of finding the second derivatives of the function of mutual
inductance between two circular filaments with respect to the appropriated
coordinates.

###  3.1  Preliminary discussion

The general scheme of arbitrarily positioning of two current-carrying circular
filaments with respect to each other is considered as shown in Fig. 1. The
linear misalignment of the secondary circle with respect to the primary one is
defined by the coordinates of the centre $C$ ($x_{c},y_{c},z_{c}$). The
angular misalignment of the secondary circle can be defined by using Grover’s
angles [38, page 207]. Namely, the angle of $\theta$ and $\eta$ corresponds to
the angular rotation around an axis passing through the diameter of the
secondary circle, and then the rotation of this axis lying on the surface
$x^{\prime}By^{\prime}$ around the vertical $z^{\prime}$ axis, respectively,
as it is shown in Figure 2(a). Accounting for Eq. (4), these two angles have
the following relationship with constants of inclined plane (1):

$\theta=\arccos\left(\frac{c}{L}\right),\;\eta=\arccos\left(\frac{-b}{L\sin\theta}\right).$
(25)

(a)

Figure 2: Two manners for determining the angular position of the secondary
circle with respect to the primary one: $x^{\prime}y^{\prime}z^{\prime}$ is
the auxiliary CF the axes of which are parallel to the axes of $XYZ$,
respectively; $x^{\prime\prime}y^{\prime\prime}z^{\prime\prime}$ is the
auxiliary CF defined in such a way that the $x^{\prime}$ and
$x^{\prime\prime}$ are coincide, but the $z^{\prime\prime}$ and
$y^{\prime\prime}$ axis is rotated by the $\alpha$ angle with respect to the
$z^{\prime}$ and $y^{\prime}$ axis, respectively.

The same angular misalignment can be determined through the $\alpha$ and
$\beta$ angle, which corresponds to the angular rotation around the
$x^{\prime}$ axis and then around the $y^{\prime\prime}$ axis, respectively,
as it is shown in Figure 2(b). This additional second manner is more
convenient in a case of study dynamics and stability issues, for instance,
applying to axially symmetric inductive levitation systems [14, 16] in
compared with Grover’s manner. These two pairs of angles have the following
relationship with respect to each other such as [32]:

$\left\\{\begin{array}[]{l}\sin\beta=\sin\eta\sin\theta;\\\
\cos\beta\sin\alpha=\cos\eta\sin\theta.\end{array}\right.$ (26)

Then, the mutual inductance between two circular filaments can be calculated
by the following formulas, which were derived by using Kalantarov-Zeitlin’s
approach in work [32] for two cases. Introducing the following dimensionless
coordinates:

${\displaystyle{x}=\frac{x_{c}}{R_{s}};\;{y}=\frac{y_{c}}{R_{s}};\;{z}=\frac{z_{c}}{R_{s}};{s}=\sqrt{{x}^{2}+{y}^{2}},}\\\
$ (27)

for the first case when the $\theta$ angle is lying in an interval of
$0\leq\theta<\pi/2$, the formula can be written as

$M=\frac{\mu_{0}\sqrt{R_{p}R_{s}}}{\pi}\int_{0}^{2\pi}{r}\cdot
U\cdot\Phi(k)d\varphi,$ (28)

where

${r}={r}(\theta,\eta)=\frac{\cos\theta}{\sqrt{\sin^{2}(\varphi-\eta)+\cos^{2}\theta\cos^{2}(\varphi-\eta)}},$
(29)
$U=U({x},{y},\theta,\eta)=\frac{R}{{\rho}^{1.5}}=\frac{{r}+t_{1}\cdot\cos\varphi+t_{2}\cdot\sin\varphi}{{\rho}^{1.5}},$
(30)
$\begin{array}[]{l}t_{1}=t_{1}({x},{y},\theta,\eta)={x}+0.5{r}^{2}\tan^{2}\theta\sin(2(\varphi-\eta))\cdot{y},\\\
t_{2}=t_{2}({x},{y},\theta,\eta)={y}-0.5{r}^{2}\tan^{2}\theta\sin(2(\varphi-\eta))\cdot{x},\\\
{\rho}={\rho}({x},{y},\theta,\eta)=\sqrt{{r}^{2}+2{r}\cdot\left({x}\cos(\varphi)+{y}\sin(\varphi)\right)+{s}^{2}},\end{array}$
(31)
$\Phi(k)=\frac{1}{k}\left[\left(1-\frac{k^{2}}{2}\right)K(k)-E(k)\right],$
(32)

and

$\begin{array}[]{l}{\displaystyle
k^{2}=k^{2}({x},{y},{z},\theta,\eta)=\frac{4\nu{\rho}}{(\nu{\rho}+1)^{2}+\nu^{2}{z}_{\lambda}^{2}},}\\\
{\displaystyle\nu=R_{s}/R_{p},\;{z}_{\lambda}={z}+{r}\tan\theta\sin(\varphi-\eta)}.\end{array}$
(33)

For the second case when the $\theta$ angle is equal to $\pi/2$ and the two
filament circles are mutually perpendicular to each other, the formula becomes

$\begin{array}[]{l}{\displaystyle
M=\frac{\mu_{0}\sqrt{R_{p}R_{s}}}{\pi}\left\\{\int_{-1}^{1}U\cdot\Phi(k)d\bar{\ell}\right.}\\\
{\displaystyle\left.+\int_{1}^{-1}U\cdot\Phi(k)d\bar{\ell}\right\\}},\end{array}$
(34)

where

$U=U({x},{y},\eta)=\frac{R}{{\rho}^{1.5}}=\frac{t_{1}-t_{2}}{{\rho}^{1.5}},$
(35)
$\begin{array}[]{l}t_{1}=t_{1}({x},\eta)=\sin\eta\cdot({x}+\bar{\ell}\cos\eta),\\\
t_{2}=t_{2}({y},\eta)=\cos\eta\cdot({y}+\bar{\ell}\sin\eta),\\\
{\rho}={\rho}({x},{y},\eta)=\sqrt{{s}^{2}+2\bar{\ell}\cdot\left({x}\cos(\eta)+{y}\sin(\eta)\right)+\bar{\ell}^{2}},\end{array}$
(36)

and $\bar{\ell}=\ell/R_{s}$ is the dimensionless variable. The functions
$\Phi(k)$ and $k=k({x},{y},{z},\eta)$ in formula (34) have the same structures
as defined by Eq. (32) and (33), respectively. Besides that, in the elliptic
module $k=k({x},{y},{z},\eta)$, the ${z}_{\lambda}$ function is governed as
follows

${z}_{\lambda}={z}\pm\sqrt{1-\bar{\ell}^{2}},$ (37)

Note that integrating formula (34) between $-1$ and $1$, Eq. (37) is
calculated with the positive sign and for the other direction the negative
sign is taken.

###  3.2  Stiffness calculation

Assuming that the primary and secondary circular filaments carry the currents
of $I_{p}$ and $I_{s}$, respectively, hence, the magnetic stiffness
corresponding to the force arising between these two current-carring circular
filaments can be calculated by taking the second derivatives of the function
of the magnetic energy stored in such the system with respect to the
appropriate coordinates. Hence, all nine components of the magnetic stiffness
can be calculated by

${\displaystyle S_{gq}=-I_{p}I_{s}\frac{\partial^{2}M}{\partial g\partial
q},}$ (38)

where $g,q=x_{c}$, $y_{c}$, or $z_{c}$. Thus, to derive formulas for
calculation of the magnetic stiffness between two arbitrarily oriented
circular filaments, the second derivatives of formulas of mutual inductance,
namely, represented by Eq. (28) and (34) must be taken. Similar to the
calculation of magnetic force in such the filament system [39], finding the
second derivatives of mutual inductance is reduced to taking the second
derivatives of their kernel functions.

####  3.2.1 For the case of $0\leq\theta<\pi/2$

Formula (28) for calculation of mutual inductance is considered. Its kernel is
defined as

$\mathrm{Kr}={r}\cdot U\cdot\Phi(k).$ (39)

According to the definitions of functions $r$, $U$, $\Phi(k)$ and $k$ given in
Eq. (29), (30), (32) and (33), respectively, the second $x_{c}$-and
$y_{c}$-derivative of kernel $\mathrm{Kr}$ can be written as

$\begin{array}[]{c}{\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial
q^{2}}=\frac{\partial^{2}\mathrm{Kr}}{\partial
g^{2}}\frac{1}{R_{s}^{2}}=\frac{r}{R_{s}^{2}}\cdot\left[\frac{\partial^{2}U}{\partial
g^{2}}\cdot\Phi(k)+2\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial g}\right.}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial
k}{\partial
g}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
g^{2}}\right)\right],}\end{array}$ (40)

where $q=x_{c},y_{c}$, and $g=x,y$,

$\begin{array}[]{l}{\displaystyle\frac{\partial U}{\partial
g}=\left({\displaystyle\frac{\partial R}{\partial g}}\cdot\rho-1.5\cdot
R\cdot\frac{\partial\rho}{\partial g}\right)\bigg{/}{{\rho}^{2.5}},}\\\
{\displaystyle\frac{\partial^{2}U}{\partial
g^{2}}=\frac{{\displaystyle\left[-0.5\frac{\partial R}{\partial
g}\frac{\partial\rho}{\partial g}-1.5R\frac{\partial^{2}\rho}{\partial
g^{2}}\right]\cdot\rho-2.5\cdot\left[{\displaystyle\frac{\partial R}{\partial
g}}\rho-1.5R\frac{\partial\rho}{\partial
g}\right]\cdot\frac{\partial\rho}{\partial g}}}{{{\rho}^{3.5}}},}\\\
{\displaystyle\frac{\partial R}{\partial g}=\frac{\partial t_{1}}{\partial
g}\cdot\cos\varphi+\frac{\partial t_{2}}{\partial g}\cdot\sin\varphi},\\\
\end{array}$ (41) $\begin{array}[]{l}{\displaystyle\frac{\partial k}{\partial
g}=\frac{G}{H}\cdot\nu\frac{\partial\rho}{\partial g},}\\\
{\displaystyle\frac{\partial^{2}k}{\partial
g^{2}}=\frac{{\displaystyle\frac{\partial G}{\partial g}H-A\frac{\partial
H}{\partial g}}}{H^{2}}\cdot\nu\frac{\partial\rho}{\partial
g}+\frac{G}{H}\cdot\nu\frac{\partial^{2}\rho}{\partial g^{2}},}\\\
{\displaystyle
G=2/k-k(\nu{\rho}+1),\;H=(\nu{\rho}+1)^{2}+\nu^{2}{z}_{\lambda}^{2},}\\\
{\displaystyle\frac{\partial G}{\partial
g}=-\left[2/k^{2}+\nu{\rho}+1\right]\frac{\partial k}{\partial
g}-k\cdot\nu\frac{\partial\rho}{\partial g},\frac{\partial H}{\partial
g}=2(\nu{\rho}+1)\nu\frac{\partial\rho}{\partial g},}\end{array}$ (42)
$\begin{array}[]{l}{\displaystyle\frac{d\Phi(k)}{dk}=\frac{1}{k^{2}}\left[\frac{2-k^{2}}{2(1-k^{2})}E(k)-K(k)\right]},\\\
{\displaystyle\frac{d^{2}\Phi(k)}{dk^{2}}=-\frac{\left(4-7k^{2}+k^{4}\right)E(k)+\left(-4+9k^{2}-5k^{4}\right)K(k)}{2k^{3}(k^{2}-1)^{2}}},\end{array}$
(43)

and when $g=x$ we have

$\begin{array}[]{l}{\displaystyle\frac{\partial t_{1}}{\partial
x}=1,\;\frac{\partial t_{2}}{\partial
x}=-0.5{r}^{2}\tan^{2}\theta\sin(2(\varphi-\eta))},\\\
{\displaystyle\frac{\partial\rho}{\partial
x}=\left(r\cdot\cos\varphi+x\right)\big{/}{{\rho}},\;\frac{\partial^{2}\rho}{\partial
x^{2}}=\left(\rho-(r\cdot\cos\varphi+x)\frac{\partial\rho}{\partial
x}\right)\big{/}{{\rho}^{2}},}\end{array}$ (44)

and when $g=y$ we have

$\begin{array}[]{l}{\displaystyle\frac{\partial t_{1}}{\partial
y}=0.5{r}^{2}\tan^{2}\theta\sin(2(\varphi-\eta)),\;\frac{\partial
t_{2}}{\partial y}=1},\\\ {\displaystyle\frac{\partial\rho}{\partial
y}=\left(r\cdot\sin\varphi+y\right)\big{/}{{\rho}},\;\frac{\partial^{2}\rho}{\partial
y^{2}}=\left(\rho-(r\cdot\sin\varphi+y)\frac{\partial\rho}{\partial
y}\right)\big{/}{{\rho}^{2}}.}\end{array}$ (45)

The second $z_{c}$-derivative of kernel $\mathrm{Kr}$ is

${\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial
z_{c}^{2}}=\frac{\partial^{2}\mathrm{Kr}}{\partial
z^{2}}\frac{1}{R_{s}^{2}}=\frac{r}{R_{s}^{2}}\cdot
U\left[\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial k}{\partial
z}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
z^{2}}\right],}$ (46)

where

$\begin{array}[]{l}{\displaystyle\frac{\partial k}{\partial
z}=-\sqrt{4\nu\rho}\cdot\frac{\nu^{2}{z}_{\lambda}}{\left((\nu{\rho}+1)^{2}+\nu^{2}{z}_{\lambda}^{2}\right)^{3/2}}},\\\
{\displaystyle\frac{\partial^{2}k}{\partial
z^{2}}=\nu^{2}\sqrt{4\nu\rho}\cdot\frac{2\nu^{2}{z}_{\lambda}^{2}-(\nu{\rho}+1)^{2}}{\left((\nu{\rho}+1)^{2}+\nu^{2}{z}_{\lambda}^{2}\right)^{5/2}}}.\end{array}$
(47)

Note that in Eq. (46), the first and second derivative of function $\Phi$ with
respect to $k$ are the same as in Eqs (43), respectively.

For further differentiation of the kernel, it is recognized that

${\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial q\partial
z_{c}}=\frac{\partial^{2}\mathrm{Kr}}{\partial z_{c}\partial
q},\;\frac{\partial^{2}\mathrm{Kr}}{\partial x_{c}\partial
y_{c}}=\frac{\partial^{2}\mathrm{Kr}}{\partial y_{c}\partial x_{c}},}$ (48)

where $q=x_{c},y_{c}$. Using properties (48) the derivation of the second
derivatives, where the variable $z_{c}$ is involved, can be simplified by
taking the first derivative of the kernel with respect to $z_{c}$. Following
this, we have

$\begin{array}[]{c}{\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial
q\partial z_{c}}=\frac{\partial^{2}\mathrm{Kr}}{\partial g\partial
z}\frac{1}{R_{s}^{2}}=\frac{r}{R_{s}^{2}}\cdot\left[\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial z}\right.}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial
k}{\partial g}\cdot\frac{\partial k}{\partial
z}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial g\partial
z}\right)\right],}\end{array}$ (49)

where $q=x_{c},y_{c}$ and $g=x,y$,

$\begin{array}[]{l}{\displaystyle\frac{\partial k}{\partial
z}=-{2}\cdot\nu^{2.5}{z}_{\lambda}\cdot\frac{\sqrt{\rho}}{H^{3/2}}},\\\
{\displaystyle\frac{\partial^{2}k}{\partial g\partial
z}=\frac{-\nu^{2.5}{z}_{\lambda}}{\sqrt{\rho}}\cdot\frac{{\displaystyle\frac{\partial\rho}{\partial
g}\cdot H-3{\rho}\cdot\frac{\partial H}{\partial g}}}{H^{5/2}},}\\\
{\displaystyle H=(\nu{\rho}+1)^{2}+\nu^{2}{z}_{\lambda}^{2},\;\frac{\partial
H}{\partial g}=2(\nu{\rho}+1)\nu\frac{\partial\rho}{\partial g}.}\end{array}$
(50)

The second derivative with respect to variables $x_{c}$ and $y_{c}$ is

$\begin{array}[]{c}{\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial
x_{c}\partial y_{c}}=\frac{\partial^{2}\mathrm{Kr}}{\partial x\partial
y}\frac{1}{R_{s}^{2}}=\frac{r}{R_{s}^{2}}\cdot\left[\frac{\partial^{2}U}{\partial
x\partial y}\Phi(k)+\frac{\partial U}{\partial
x}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial y}\right.}\\\
{\displaystyle\left.+\frac{\partial U}{\partial
y}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial
x}+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial k}{\partial
x}\cdot\frac{\partial k}{\partial
y}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial x\partial
y}\right)\right],}\end{array}$ (51)

where

$\begin{array}[]{l}{\displaystyle\frac{\partial^{2}U}{\partial y\partial
x}=\left(\frac{\partial G_{x}}{\partial
y}\rho-2.5G_{x}\frac{\partial\rho}{\partial
y}\right)\bigg{/}{{\rho}^{3.5}},}\\\ {\displaystyle
G_{x}={\displaystyle\frac{\partial R}{\partial x}}\cdot\rho-1.5\cdot
R\cdot\frac{\partial\rho}{\partial x},}\\\ {\displaystyle\frac{\partial
G_{x}}{\partial y}={\displaystyle\frac{\partial R}{\partial
x}}\cdot\frac{\partial\rho}{\partial y}-1.5\cdot{\displaystyle\frac{\partial
R}{\partial y}}\cdot\frac{\partial\rho}{\partial
x}-1.5R\frac{\partial^{2}\rho}{\partial y\partial x},}\\\
{\displaystyle\frac{\partial^{2}\rho}{\partial y\partial
x}=-\frac{\partial\rho}{\partial x}\cdot\frac{\partial\rho}{\partial
y}\cdot\frac{1}{\rho}},\end{array}$ (52)
$\begin{array}[]{l}{\displaystyle\frac{\partial^{2}k}{\partial y\partial
x}=\frac{{\displaystyle\frac{\partial G}{\partial y}H-G\frac{\partial
H}{\partial y}}}{H^{2}}\cdot\nu\frac{\partial\rho}{\partial
x}+\frac{G}{H}\cdot\nu\frac{\partial^{2}\rho}{\partial y\partial x},}\\\
{\displaystyle
G=2/k-k(\nu{\rho}+1),\;H=(\nu{\rho}+1)^{2}+\nu^{2}{z}_{\lambda}^{2},}\\\
{\displaystyle\frac{\partial G}{\partial
y}=-\left[2/k^{2}+\nu{\rho}+1\right]\frac{\partial k}{\partial
y}-k\cdot\nu\frac{\partial\rho}{\partial y},\frac{\partial H}{\partial
y}=2(\nu{\rho}+1)\nu\frac{\partial\rho}{\partial y}.}\end{array}$ (53)

Other derivatives of functions $\Phi$, $k$, $U$ in Eqs (49) and (51) are
defined in the same way as it has been shown above.

Hence, for this particular case when $0\leq\theta<\pi/2$, according to Eqs
(40), (46), (49) and (51) all nine components of magnetic stiffness can be
calculated as follows:

$\begin{array}[]{l}{\displaystyle
S_{qq}=-\frac{\mu_{0}I_{p}I_{s}\sqrt{R_{p}}}{\pi
R_{s}^{3/2}}\int_{0}^{2\pi}{r}\cdot\left[\frac{\partial^{2}U}{\partial
g^{2}}\cdot\Phi(k)+2\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial g}\right.}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial
k}{\partial
g}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
g^{2}}\right)\right]d\varphi;}\end{array}$ (54)
$\begin{array}[]{l}{\displaystyle
S_{z_{c}z_{c}}=-\frac{\mu_{0}I_{p}I_{s}\sqrt{R_{p}}}{\pi
R_{s}^{3/2}}\int_{0}^{2\pi}{r}\cdot
U\left[\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial k}{\partial
z}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
z^{2}}\right]d\varphi;}\end{array}$ (55) $\begin{array}[]{l}{\displaystyle
S_{z_{c}q}=S_{qz_{c}}=-\frac{\mu_{0}I_{p}I_{s}\sqrt{R_{p}}}{\pi
R_{s}^{3/2}}\int_{0}^{2\pi}{r}\cdot\left[\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial z}\right.}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial
k}{\partial g}\cdot\frac{\partial k}{\partial
z}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial g\partial
z}\right)\right]d\varphi;}\end{array}$ (56) $\begin{array}[]{l}{\displaystyle
S_{x_{c}y_{c}}=S_{y_{c}x_{c}}=-\frac{\mu_{0}I_{p}I_{s}\sqrt{R_{p}}}{\pi
R_{s}^{3/2}}\int_{0}^{2\pi}r\cdot\left[\frac{\partial^{2}U}{\partial x\partial
y}\Phi(k)+\frac{\partial U}{\partial
x}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial y}\right.}\\\
{\displaystyle\left.+\frac{\partial U}{\partial
y}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial
x}+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial k}{\partial
x}\cdot\frac{\partial k}{\partial
y}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial x\partial
y}\right)\right]d\varphi,}\end{array}$ (57)

where $q=x_{c}$ or $y_{c}$, $g=x$ or $y$, respectively.

####  3.2.2 The second case: $\theta=\pi/2$

For this case, formula (34) for calculation of mutual inductance is used. Its
kernel is defined as

$\mathrm{Kr}=U\cdot\Phi(k).$ (58)

Accounting for that in this case the function $U$ is defined by Eq. (35), then
the second $x_{c}$-and $y_{c}$-derivative of kernel $\mathrm{Kr}$ can be
written as

$\begin{array}[]{c}{\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial
q^{2}}=\frac{\partial^{2}\mathrm{Kr}}{\partial
g^{2}}\frac{1}{R_{s}^{2}}=\frac{1}{R_{s}^{2}}\cdot\left[\frac{\partial^{2}U}{\partial
g^{2}}\cdot\Phi(k)+2\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial g}\right.}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial
k}{\partial
g}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
g^{2}}\right)\right],}\end{array}$ (59)

where $q=x_{c},y_{c}$, and $g=x,y$. The first and second derivatives of
function $U$ with respect to $g$ are taken analogically as shown in Eqs (41).
Similar to Eqs (42), the derivatives of function $k$ are defined. When $g=x$
we have

$\begin{array}[]{l}{\displaystyle\frac{\partial R}{\partial x}=\frac{\partial
t_{1}}{\partial x}=\sin\eta,\frac{\partial\rho}{\partial
x}=\frac{x+\bar{\ell}\cos\eta}{\rho}}\\\
{\displaystyle\frac{\partial^{2}\rho}{\partial
x^{2}}=\frac{\rho-(x+\bar{\ell}\cos\eta){\displaystyle\frac{\partial\rho}{\partial
x}}}{\rho^{2}}},\end{array}$ (60)

and when $g=y$ we can write

$\begin{array}[]{l}{\displaystyle\frac{\partial R}{\partial y}=\frac{\partial
t_{2}}{\partial y}=-\cos\eta,\frac{\partial\rho}{\partial
y}=\frac{y+\bar{\ell}\sin\eta}{\rho}}\\\
{\displaystyle\frac{\partial^{2}\rho}{\partial
y^{2}}=\frac{\rho-(y+\bar{\ell}\sin\eta){\displaystyle\frac{\partial\rho}{\partial
y}}}{\rho^{2}}}.\end{array}$ (61)

The second $z_{c}$-derivative of the kernel $\mathrm{Kr}$ is defined similar
to Eq. (46) as

${\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial
z_{c}^{2}}=\frac{\partial^{2}\mathrm{Kr}}{\partial
z^{2}}\frac{1}{R_{s}^{2}}=\frac{1}{R_{s}^{2}}\cdot
U\left[\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial k}{\partial
z}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
z^{2}}\right],}$ (62)

where the derivatives of function $k$ have the same meaning as in Eqs. (47).
Using properties (48) the derivation of the second derivatives with respect to
the variables $z_{c}$ and $x_{c}$, and $z_{c}$ and $y_{c}$ can be written
similar to Eq. (49) as follows

$\begin{array}[]{c}{\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial
q\partial z_{c}}=\frac{\partial^{2}\mathrm{Kr}}{\partial g\partial
z}\frac{1}{R_{s}^{2}}=\frac{1}{R_{s}^{2}}\cdot\left[\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial z}\right.}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial
k}{\partial g}\cdot\frac{\partial k}{\partial
z}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial g\partial
z}\right)\right],}\end{array}$ (63)

where $q=x_{c},y_{c}$ and $g=x,y$, taking into account Eqs. (60) and (61) the
first and the second derivatives of $U$, $k$ are determined by Eqs. (41) and
(50), respectively.

The second derivative of the kernel with respect to variables $x_{c}$ and
$y_{c}$ is defined in a similar way to Eq. (51) and is as follows

$\begin{array}[]{c}{\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial
x_{c}\partial y_{c}}=\frac{\partial^{2}\mathrm{Kr}}{\partial x\partial
y}\frac{1}{R_{s}^{2}}=\frac{1}{R_{s}^{2}}\cdot\left[\frac{\partial^{2}U}{\partial
x\partial y}\Phi(k)+\frac{\partial U}{\partial
x}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial y}\right.}\\\
{\displaystyle\left.+\frac{\partial U}{\partial
y}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial
x}+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial k}{\partial
x}\cdot\frac{\partial k}{\partial
y}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial x\partial
y}\right)\right].}\end{array}$ (64)

In Eq. (64) all derivatives of functions $U$ and $k$ are determined in Eq.
(51), but the definitions of function $R$ and $\rho$ and their respective
derivatives given in Eqs. (60) and (61) must be taken into account.

Hence, for this case when $\theta=\pi/2$, according to Eqs (59), (63) and (64)
all nine components of magnetic stiffness can be calculated as follows:

$\begin{array}[]{l}{\displaystyle
S_{qq}=-\frac{\mu_{0}I_{p}I_{s}\sqrt{R_{p}}}{\pi
R_{s}^{3/2}}\left[\int_{-1}^{1}\frac{\partial^{2}U}{\partial
g^{2}}\cdot\Phi(k)+2\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial g}\right.}\\\
{\displaystyle+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial
k}{\partial
g}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
g^{2}}\right)d\bar{\ell}+\int_{1}^{-1}\frac{\partial^{2}U}{\partial
g^{2}}\cdot\Phi(k)}\\\ {\displaystyle\left.+2\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial
g}+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial k}{\partial
g}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
g^{2}}\right)d\bar{\ell}\right]};\end{array}$ (65)
$\begin{array}[]{l}{\displaystyle
S_{z_{c}z_{c}}=-\frac{\mu_{0}I_{p}I_{s}\sqrt{R_{p}}}{\pi
R_{s}^{3/2}}\left\\{\int_{-1}^{1}U\left[\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial
k}{\partial
z}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
z^{2}}\right]d\bar{\ell}\right.}\\\
{\displaystyle\left.+\int_{1}^{-1}U\left[\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\left(\frac{\partial
k}{\partial
z}\right)^{2}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial
z^{2}}\right]d\bar{\ell}\right\\}};\end{array}$ (66)
$\begin{array}[]{l}{\displaystyle
S_{z_{c}q}=S_{qz_{c}}=-\frac{\mu_{0}I_{p}I_{s}\sqrt{R_{p}}}{\pi
R_{s}^{3/2}}\left[\int_{-1}^{1}\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial z}\right.}\\\
{\displaystyle+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial
k}{\partial g}\cdot\frac{\partial k}{\partial
z}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial g\partial
z}\right)d\bar{\ell}+\int_{1}^{-1}\frac{\partial U}{\partial
g}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial z}}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial
k}{\partial g}\cdot\frac{\partial k}{\partial
z}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial g\partial
z}\right)d\bar{\ell}\right];}\end{array}$ (67)
$\begin{array}[]{l}{\displaystyle
S_{x_{c}y_{c}}=S_{y_{c}x_{c}}=-\frac{\mu_{0}I_{p}I_{s}\sqrt{R_{p}}}{\pi
R_{s}^{3/2}}\left[\int_{-1}^{1}\frac{\partial^{2}U}{\partial x\partial
y}\Phi(k)+\frac{\partial U}{\partial
x}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial y}\right.}\\\
{\displaystyle+\frac{\partial U}{\partial
y}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial
x}+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial k}{\partial
x}\cdot\frac{\partial k}{\partial
y}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial x\partial
y}\right)d\bar{\ell}}\\\
{\displaystyle+\int_{1}^{-1}\frac{\partial^{2}U}{\partial x\partial
y}\Phi(k)+\frac{\partial U}{\partial
x}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial y}+\frac{\partial
U}{\partial y}\cdot\frac{d\Phi(k)}{dk}\cdot\frac{\partial k}{\partial x}}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Phi(k)}{dk^{2}}\cdot\frac{\partial
k}{\partial x}\cdot\frac{\partial k}{\partial
y}+\frac{d\Phi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial x\partial
y}\right)d\bar{\ell}\right],}\end{array}$ (68)

where $q=x_{c}$ or $y_{c}$, $g=x$ or $y$, respectively.

Thus, the set of formulas (54)-(57) and (65)-(68) for calculation of all nine
components of magnetic stiffness of the corresponding force arising between
two current-carrying circular filaments arbitrarily oriented in the space are
derived by using the mutual inductance method. The derived formulas are
intuitively understandable for application, they can be easily programmed. For
this purpose, the Matlab language was used. The Matlab files with the
implemented formulas are available as supplementary materials to this article.

##  4 Numerical verification of derived formulas

Developed sets of formulas for calculation of nine components of magnetic
stiffness of corresponding force between two current-carrying circular
filaments derived by means of Babic’s method Eqs (15)-(23) and the method of
mutual inductance (Kalantarov-Zeitlin’s method) Eqs (54)-(57) and Eqs
(65)-(68) are mutually verified to each other through applying developed
formulas to a number of examples designed in this section. In all examples
bellow, it is assumed that the carrying currents in both coils are equal to
one ampere ($I_{p}=I_{s}=$1\text{\,}\mathrm{A}$$). In addition to the
calculation of components of magnetic stiffness of the considered filament
system is supported by the set of expressions (83)-(85) derived from Grover’s
formula for calculation of mutual inductance [38, page 207, Eq. (179)]. The
derivations of these expressions are shown in A. All calculations for
considered cases proved the robustness and efficiency of developed formulas.

Figure 3: Geometrical scheme of circular filaments with parallel axes denoted
via Grover’s notation: $\rho$ is the distance between axes; $d$ is the
distance between the coils’ planes $d=z_{c}$.

###  4.1 Magnetic stiffness between circular filaments with parallel axes

The scheme for calculation of magnetic stiffness between circular filaments
with parallel axes is shown in Fig. 3. The linear misalignment in the Grover
notation can be defined by the geometrical parameter, $d$, which is the
distance between the planes of circles and the parameter, $\rho$, is the
distance between their axes. These parameters have the following relationship
to the notation defined in this article, namely, $z_{c}=d$ and
$\rho=\sqrt{x_{c}^{2}+y_{c}^{2}}$. Fig. 3 shows the particular case, when
$\rho=y_{c}$.

####  Example 1: (Example 16, page 74 in Babič’s work [36])

Two coaxial circular filaments for which the primary coil has a radius of
$R_{p}=2$ $\mathrm{m}$, and the secondary one $R_{s}=1$ $\mathrm{m}$. The
axial distance between filaments is $d=z_{c}=1$$\mathrm{m}$. The results of
calculation of diagonal and non-diagonal components of magnetic stiffness are
as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (54), (55) |  Eqs (83), (84)
$S_{xx}$ | $1.032010586220236\times 10^{-7}$ | $1.032010586220216\times 10^{-7}$ | $1.032010586220211\times 10^{-7}$
$S_{yy}$ | $1.032010586220238\times 10^{-7}$ | $1.032010586220216\times 10^{-7}$ | $1.032010586220211\times 10^{-7}$
$S_{zz}$ | $-2.064021172440475\times 10^{-7}$ | $-2.064021172440485\times 10^{-7}$ | $-2.064021172440499\times 10^{-7}$

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (56), (57) |  Eqs (85)
$S_{xy}$ | $-1.346233158563486\times 10^{-142}$ | $8.298183971090421\times 10^{-23}$ | Not Applicable (NA)
$S_{yz}$ | $2.184515917637124\times 10^{-142}$ | $-6.010510193398827\times 10^{-23}$ | $1.739957438966064\times 10^{-23}$
$S_{xz}$ | $-2.533487747768601\times 10^{-142}$ | $-5.887846720064157\times 10^{-23}$ | $1.739957438966064\times 10^{-23}$

Analysis of the results of calculation shows that they are agree well to each
other. The difference arises after thirteenth digit in a resulting number. In
the calculation of non-diagonal components, the order of magnitude
corresponding to minus twenty three can be considered as approximately equal
zero. In Babic’s method alternatively to Eqs (16), (20), (17) the other set of
formulas, namely, Eqs (18), (22), (21) can be used for calculation of non-
diagonal components. The results are the same and equal to zeros. Also, note
that the sum of diagonal components in all methods is almost zero and the
condition (24) is held with an accuracy of minus twenty one order of
magnitude. Worth noting that in some cases this condition helps us also to
restore the missed component of stiffness in the orthogonal direction in
Grover’s method.

####  Example 2

Let us consider the coils having the same radii as in the previous example 1,
but the center of the secondary coil is located at the point
$x_{c}=0$$\mathrm{m}$, $y_{c}=0.5$$\mathrm{m}$ and $z_{c}=1$$\mathrm{m}$,
which corresponds to the following Grover parameters: $\rho=0.5$$\mathrm{m}$
and $d=1$$\mathrm{m}$. The results of calculation of diagonal and non-diagonal
components of magnetic stiffness are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (54), (55) |  Eqs (83), (84)
$S_{xx}$ | $1.402100143236235\times 10^{-7}$ | $1.402100143236235\times 10^{-7}$ | $1.402100143236236\times 10^{-7}$
$S_{yy}$ | $2.118309158188127\times 10^{-7}$ | $2.118309158188122\times 10^{-7}$ | $2.118309158188126\times 10^{-7}$
$S_{zz}$ | $-3.520409301424362\times 10^{-7}$ | $-3.520409301424361\times 10^{-7}$ | $-3.520409301424361\times 10^{-7}$

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (56), (57) |  Eqs (85)
$S_{xy}$ | $3.873036680241143\times 10^{-141}$ | $-2.629904868295323\times 10^{-22}$ | NA
$S_{yz}$ | $9.040026778347652\times 10^{-8}$ | $9.040026778347634\times 10^{-8}$ | $9.040026778347627\times 10^{-8}$
$S_{xz}$ | $-1.38528525192781\times 10^{-140}$ | $9.028031637431708\times 10^{-23}$ | NA

The component $S_{xx}$ in GM was restored form the condition (24).

####  Example 3

The coils having the same radii as in the previous examples and the center of
the secondary coil is located at the point $x_{c}=0$$\mathrm{m}$,
$y_{c}=0$$\mathrm{m}$ and $z_{c}=0$$\mathrm{m}$. Grover parameters are zeros.
The results of calculation of diagonal and non-diagonal components of magnetic
stiffness are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (54), (55) |  Eqs (83), (84)
$S_{xx}$ | $-6.367128613342259\times 10^{-7}$ | $-6.367128613342232\times 10^{-7}$ | $-6.367128613342277\times 10^{-7}$
$S_{yy}$ | $-6.367128613342259\times 10^{-7}$ | $-6.36712861334223\times 10^{-7}$ | $-6.367128613342218\times 10^{-7}$
$S_{zz}$ | $1.273425722668452\times 10^{-6}$ | $1.273425722668452\times 10^{-6}$ | $1.273425722668449\times 10^{-6}$

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (56), (57) |  Eqs (85)
$S_{xy}$ | $-1.752675194540019\times 10^{-141}$ | $-2.271727526158087\times 10^{-22}$ | NA
$S_{yz}$ | 0 | $0$ | 0
$S_{xz}$ | 0 | $0$ | 0

Figure 4: Geometrical scheme of mutually perpendicular current-carrying
circular filaments (the second case: $\theta=\pi/2$).

###  4.2 Magnetic stiffness between mutually perpendicular circular filaments

In this section, magnetic stiffness between mutually perpendicular current-
carrying filaments is calculated. The general scheme is shown in Fig. 4. For
the mutual inductance method, it corresponds to the second case, when
$\theta=\pi/2$ and the angular misalignment is only characterized by the angle
$\eta$, and the set of formulas (65)-(68) is used.

####  Example 4

The two mutually perpendicular coils have the following radii, namely,
$R_{p}=0.2$ $\mathrm{m}$, and the secondary one $R_{s}=0.1$ $\mathrm{m}$. The
center of the secondary coil is located at the origin as shown in Fig. 5(a).
The angle $\eta$ is zero $\mathrm{rad}$ for MIM. For BM, the angular
misalignment is characterized by the following components of plane equation
(1): $a=c=0$ and $b=1$. The results of calculation of diagonal and non-
diagonal components of magnetic stiffness are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (65), (66) |  Eqs (83), (84)
$S_{xx}$ | 0 | 0 | Not a Number (NaN)
$S_{yy}$ | 0 | 0 | NaN
$S_{zz}$ | 0 | 0 | NaN

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (67), (68) |  Eqs (85)
$S_{xy}$ | 0 | 0 | NA
$S_{yz}$ | $2.706560599934499\times 10^{-6}$ | $2.706560599933974\times 10^{-6}$ | NaN
$S_{xz}$ | 0 | $0$ | NA

Note that for Grover’s method it is the singular case.

(a) Example 4

(b) Example 5

Figure 5: Schemes of configuration of mutually perpendicular current-carrying
coils.

####  Example 5

The same arrangement of coils as in example 4 is used, but the secondary coil
is turned on the angle $\eta=\pi/2$ $\mathrm{rad}$ as shown in Fig. 5(b). For
BM, the angular misalignment is characterized by the following components of
plane equation (1): $b=c=0$. The results of calculation are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (65), (66) |  Eqs (83), (84)
$S_{xx}$ | 0 | 0 | Not a Number (NaN)
$S_{yy}$ | 0 | 0 | NaN
$S_{zz}$ | 0 | 0 | NaN

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (67), (68) |  Eqs (85)
$S_{xy}$ | 0 | 0 | NA
$S_{yz}$ | 0 | $1.657290387703741\times 10^{-22}$ | NaN
$S_{xz}$ | $-2.706560599933975\times 10^{-6}$ | $-2.706560599933974\times 10^{-6}$ | NA

Note that for Grover’s method it is the singular case.

####  Example 6

The two mutually perpendicular coils have the following radii, namely,
$R_{p}=1.0$ $\mathrm{m}$, and the secondary one $R_{s}=0.5$ $\mathrm{m}$. The
center of the secondary coil is located at the point $x_{c}=0$$\mathrm{m}$,
$y_{c}=2$$\mathrm{m}$ and $z_{c}=0$$\mathrm{m}$ . The angle $\eta$ is zero
$\mathrm{rad}$ for MIM. For BM, the angular misalignment is characterized by
the following components of plane equation (1): $a=c=0$ and $b=1$. The results
of calculation of diagonal and non-diagonal components of magnetic stiffness
are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (65), (66) |  Eqs (83), (84)
$S_{xx}$ | 0 | 0 | $-3.68315244876568\times 10^{-24}$
$S_{yy}$ | $1.19859223988025\times 10^{-143}$ | 0 | $8.191673408378001\times 10^{-24}$
$S_{zz}$ | $4.794368959521\times 10^{-142}$ | 0 | $-4.508520959612321\times 10^{-24}$

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (67), (68) |  Eqs (85)
$S_{xy}$ | 0 | 0 | NA
$S_{yz}$ | $-1.368742764885786\times 10^{-7}$ | $-1.368742764885786\times 10^{-7}$ | $-1.368742764885786\times 10^{-7}$
$S_{xz}$ | 0 | $5.151865880056137\times 10^{-24}$ | NA

####  Example 7

The two mutually perpendicular coils have the following radii, namely,
$R_{p}=1.0$ $\mathrm{m}$, and the secondary one $R_{s}=0.5$ $\mathrm{m}$. The
center of the secondary coil is located at the point $x_{c}=0$$\mathrm{m}$,
$y_{c}=2$$\mathrm{m}$ and $z_{c}=3$$\mathrm{m}$ . The angle $\eta$ is zero
$\mathrm{rad}$ for MIM. For BM, the angular misalignment is characterized by
the following components of plane equation (1): $a=c=0$ and $b=1$. The results
of calculation of diagonal and non-diagonal components of magnetic stiffness
are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (65), (66) |  Eqs (83), (84)
$S_{xx}$ | $-2.262682905005172\times 10^{-9}$ | $-2.262682905005171\times 10^{-9}$ | $-2.262682905005174\times 10^{-9}$
$S_{yy}$ | $-2.710919377082796\times 10^{-9}$ | $-2.710919377082797\times 10^{-9}$ | $-2.710919377082796\times 10^{-9}$
$S_{zz}$ | $4.973602282087952\times 10^{-9}$ | $4.973602282087968\times 10^{-9}$ | $4.97360228208797\times 10^{-9}$

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (67), (68) |  Eqs (85)
$S_{xy}$ | $-2.364186955151571\times 10^{-143}$ | $-2.347855544296417\times 10^{-25}$ | NA
$S_{yz}$ | $3.101402573517489\times 10^{-9}$ | $3.101402573517489\times 10^{-9}$ | $3.101402573517489\times 10^{-9}$
$S_{xz}$ | $7.984780962755199\times 10^{-143}$ | $8.433113791758558\times 10^{-26}$ | NA

Figure 6: Geometrical scheme of mutually perpendicular current-carrying
circular filaments for Example 8.

####  Example 8

The two mutually perpendicular coils have the following radii, namely,
$R_{p}=1.0$ $\mathrm{m}$, and the secondary one $R_{s}=0.5$ $\mathrm{m}$. The
center of the secondary coil is located at the point $x_{c}=1$$\mathrm{m}$,
$y_{c}=2$$\mathrm{m}$ and $z_{c}=3$$\mathrm{m}$. The angle $\eta$ is $\pi/2$
$\mathrm{rad}$ for MIM as shown in Fig. 6. For BM, the angular misalignment is
characterized by the following components of plane equation (1): $b=c=0$ and
$a=1$. The results of calculation are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (65), (66) |  Eqs (83), (84)
$S_{xx}$ | $2.444411106760408\times 10^{-9}$ | $2.444411106760407\times 10^{-9}$ | NA
$S_{yy}$ | $-6.546286516635751\times 10^{-10}$ | $-6.546286516635743\times 10^{-10}$ | NA
$S_{zz}$ | $-1.789782455096833\times 10^{-9}$ | $-1.789782455096833\times 10^{-9}$ | NA

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (67), (68) |  Eqs (85)
$S_{xy}$ | $1.042889962848133\times 10^{-9}$ | $1.042889962848133\times 10^{-9}$ | NA
$S_{yz}$ | $-2.190346410345056\times 10^{-9}$ | $-2.190346410345057\times 10^{-9}$ | NA
$S_{xz}$ | $1.067688019471112\times 10^{-9}$ | $1.067688019471113\times 10^{-9}$ | NA

###  4.3 Magnetic stiffness between circular filaments arbitrarily positioned
in the space

In this section, using the equation of inclined plane (1) for BM and its
relationship with the Grover’s angles (25) for MIM to define different angular
misalignments of the secondary coil with respect to the primary one, a number
of examples with different arrangements of coils for calculation of magnetic
stiffness are designed and considered below. The calculation is accompanied by
the evaluation of stiffness by means of Grover’s formulas Eqs (83), (84) and
(85) when they are applicable.

Figure 7: Geometrical scheme for Example 9.

####  Example 9

The two coils have the following radii, namely, the primery one $R_{p}=0.4$
$\mathrm{m}$, and the secondary one $R_{s}=0.05$ $\mathrm{m}$. The center of
the secondary coil is located at the point $x_{c}=0.1$$\mathrm{m}$,
$y_{c}=0.15$$\mathrm{m}$ and $z_{c}=0$$\mathrm{m}$. For BM, the angular
misalignment is characterized by the following plane equation, namely,
$3x+2y+z=0.6$. According to Eqs. (25), it corresponds to the angle
$\eta=2.15879893034246$ $\mathrm{rad}$
($123.69006752598$$\mathrm{\SIUnitSymbolDegree}$) and
$\theta=1.30024656381632$ $\mathrm{rad}$
($74.498640433063$$\mathrm{\SIUnitSymbolDegree}$) in notations of MIM. The
coils’ arrangement is shown in Fig. 7. The results of calculation are as
follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (54), (55) |  Eqs (83), (84)
$S_{xx}$ | $-5.327433787498592\times 10^{-8}$ | $-5.32743378749859\times 10^{-8}$ | NA
$S_{yy}$ | $-7.00453721025121\times 10^{-8}$ | $-7.004537210251179\times 10^{-8}$ | NA
$S_{zz}$ | $1.23319709977498\times 10^{-7}$ | $1.23319709977498\times 10^{-7}$ | NA

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (56), (57) |  Eq (85)
$S_{xy}$ | $-1.500101622143787\times 10^{-8}$ | $-1.500101622143775\times 10^{-8}$ | NA
$S_{yz}$ | $-2.121823074384375\times 10^{-7}$ | $-2.121823074384376\times 10^{-7}$ | NA
$S_{xz}$ | $-2.167771909796546\times 10^{-7}$ | $-2.167771909796545\times 10^{-7}$ | NA

Figure 8: Geometrical scheme for Example 10.

####  Example 10

The primery coil has a radius of $R_{p}=4$ $\mathrm{m}$, and the secondary one
has a radius of $R_{s}=2$ $\mathrm{m}$. The center of the secondary coil is
located at the point $x_{c}=1$$\mathrm{m}$, $y_{c}=1$$\mathrm{m}$ and
$z_{c}=-1$$\mathrm{m}$. For BM, the angular misalignment is characterized by
the following plane equation, namely, $x+2y+3z=0$. According to Eqs. (25), the
angles $\eta$ and $\theta$ are $2.67794504458899$ $\mathrm{rad}$
($153.434948822922$$\mathrm{\SIUnitSymbolDegree}$) and $0.640522312679424$
$\mathrm{rad}$ ($36.6992252004899$$\mathrm{\SIUnitSymbolDegree}$),
respectively, in notations of MIM. The coils’ arrangement is shown in Fig. 8.
The results of calculation are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (54), (55) |  Eqs (83), (84)
$S_{xx}$ | $2.181662870952764\times 10^{-8}$ | $2.181662870952587\times 10^{-8}$ | NA
$S_{yy}$ | $3.44815074198134\times 10^{-8}$ | $3.448150741981526\times 10^{-8}$ | NA
$S_{zz}$ | $-5.629813612934105\times 10^{-8}$ | $-5.629813612934103\times 10^{-8}$ | NA

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (56), (57) |  Eqs (85)
$S_{xy}$ | $5.713524254486024\times 10^{-8}$ | $5.713524254486125\times 10^{-8}$ | NA
$S_{yz}$ | $-1.070006627660674\times 10^{-7}$ | $-1.07000662766067\times 10^{-7}$ | NA
$S_{xz}$ | $-1.105341440732599\times 10^{-7}$ | $-1.1053414407326\times 10^{-7}$ | NA

Figure 9: Geometrical scheme for Example 11.

####  Example 11

The radii of the coils are the same as in the previous example 10, but the
center of the secondary coil is located at the point $x_{c}=-1$$\mathrm{m}$,
$y_{c}=1$$\mathrm{m}$ and $z_{c}=1$$\mathrm{m}$. The angular misalignment for
BM is characterized by the following plane equation, namely, $3x+2y+1z=0$.
According to Eqs. (25), the angles $\eta$ and $\theta$ are $2.15879893034246$
$\mathrm{rad}$ ($123.69006752598$$\mathrm{\SIUnitSymbolDegree}$) and
$1.30024656381632$ $\mathrm{rad}$
($74.498640433063$$\mathrm{\SIUnitSymbolDegree}$), respectively, in notations
of MIM. The coils’ arrangement is shown in Fig. 9. The results of calculation
are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (54), (55) |  Eqs (83), (84)
$S_{xx}$ | $2.944440845967626\times 10^{-8}$ | $2.944440845966999\times 10^{-8}$ | NA
$S_{yy}$ | $-5.401491525386883\times 10^{-8}$ | $-5.401491525387853\times 10^{-8}$ | NA
$S_{zz}$ | $2.457050679419257\times 10^{-8}$ | $2.457050679419274\times 10^{-8}$ | NA

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (56), (57) |  Eqs (85)
$S_{xy}$ | $-7.826596371018449\times 10^{-9}$ | $-7.826596371018387\times 10^{-9}$ | NA
$S_{yz}$ | $-2.426184704880834\times 10^{-8}$ | $-2.426184704880832\times 10^{-8}$ | NA
$S_{xz}$ | $-1.147221843342089\times 10^{-7}$ | $-1.14722184334209\times 10^{-7}$ | NA

Figure 10: Geometrical scheme for Example 12: the linked coils.

####  Example 12

The radii of each coil are the same and equal to $1$ $\mathrm{m}$. The center
of the secondary coil is located at the point $x_{c}=-1$$\mathrm{m}$,
$y_{c}=1$$\mathrm{m}$ and $z_{c}=0$$\mathrm{m}$. The angular misalignment for
BM is characterized by the following plane equation, namely, $x+y+z=0$.
According to Eqs. (25), the angles $\eta$ and $\theta$ are $2.35619449019234$
$\mathrm{rad}$ ($135$$\mathrm{\SIUnitSymbolDegree}$) and $0.955316618124509$
$\mathrm{rad}$ ($54.7356103172453$$\mathrm{\SIUnitSymbolDegree}$),
respectively, in notations of MIM. The arrangement of the linked coils is
shown in Fig. 10. The results of calculation are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (54), (55) |  Eqs (83), (84)
$S_{xx}$ | $1.046477792966786\times 10^{-7}$ | $1.046477792966789\times 10^{-7}$ | NaN
$S_{yy}$ | $1.046477792966786\times 10^{-7}$ | $1.046477792966785\times 10^{-7}$ | NaN
$S_{zz}$ | $-2.092955585933573\times 10^{-8}$ | $-2.092955585933573\times 10^{-7}$ | NaN

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (56), (57) |  Eqs (85)
$S_{xy}$ | $2.014629707091004\times 10^{-7}$ | $2.014629707091005\times 10^{-7}$ | NA
$S_{yz}$ | $-2.577031542995681\times 10^{-7}$ | $-2.577031542995682\times 10^{-7}$ | NaN
$S_{xz}$ | $-2.577031542995681\times 10^{-7}$ | $-2.577031542995681\times 10^{-7}$ | NA

####  Example 13

The radii of coils and their angular orienation with respect to each other are
the same as in Example 12. The center of the secondary coil is located at the
point $x_{c}=0$$\mathrm{m}$, $y_{c}=1$$\mathrm{m}$ and $z_{c}=-1$$\mathrm{m}$.
The results of calculation are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (54), (55) |  Eqs (83), (84)
$S_{xx}$ | $2.986090969216481\times 10^{-7}$ | $2.986090969216481\times 10^{-7}$ | $2.986090969216485\times 10^{-7}$
$S_{yy}$ | $-6.600031210636617\times 10^{-8}$ | $-6.600031210636671\times 10^{-8}$ | $-6.600031210636627\times 10^{-8}$
$S_{zz}$ | $-2.326087848152819\times 10^{-8}$ | $-2.326087848152818\times 10^{-7}$ | $-2.326087848152822\times 10^{-7}$

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (56), (57) |  Eqs (85)
$S_{xy}$ | $-3.199237389569022\times 10^{-8}$ | $-3.199237389569028\times 10^{-8}$ | NA
$S_{yz}$ | $8.115657627384119\times 10^{-7}$ | $8.115657627384113\times 10^{-7}$ | $8.11565762738412\times 10^{-8}$
$S_{xz}$ | $2.749069587302089\times 10^{-7}$ | $2.749069587302088\times 10^{-7}$ | NA

Figure 11: Geometrical scheme for Example 14.

####  Example 14

The primary circular filament has a radius of $R_{p}=0.2$ $\mathrm{m}$, while
the secondary one has $R_{s}=0.1$ $\mathrm{m}$. The centre of secondary coil
is located at the point $C$ having the following coordinates
$x_{c}=0.1$$\mathrm{m}$, $y_{c}=0.1$$\mathrm{m}$ and $z_{c}=0.1$$\mathrm{m}$.
The angular misalignment for BM is defined by the following plane equation,
namely, $x+y+z=0.3$. According to Eqs. (25), the angles $\eta$ and $\theta$
are $2.35619449019234$ $\mathrm{rad}$ ($135$$\mathrm{\SIUnitSymbolDegree}$)
and $0.955316618124509$ $\mathrm{rad}$
($54.7356103172453$$\mathrm{\SIUnitSymbolDegree}$), respectively, in notations
of MIM. The arrangement of the coils is shown in Fig. 11. The results of
calculation are as follows

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (15), (19), (23) |  Eqs (54), (55) |  Eqs (83), (84)
$S_{xx}$ | $2.868431152918931\times 10^{-5}$ | $2.868431152918925\times 10^{-5}$ | NA
$S_{yy}$ | $2.868431152918931\times 10^{-5}$ | $2.868431152918923\times 10^{-5}$ | NA
$S_{zz}$ | $-5.736862305837862\times 10^{-5}$ | $-5.736862305837814\times 10^{-5}$ | NA

|  BM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  MIM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$ |  GM, $\mathrm{N}\text{\,}{\mathrm{m}}^{-1}$
---|---|---|---
|  Eqs (16), (20), (17) |  Eqs (56), (57) |  Eqs (85)
$S_{xy}$ | $2.397140500000452\times 10^{-5}$ | $2.397140500000446\times 10^{-5}$ | NA
$S_{yz}$ | $1.724024033513611\times 10^{-6}$ | $1.724024033513529\times 10^{-6}$ | NA
$S_{xz}$ | $1.724024033513611\times 10^{-6}$ | $1.724024033513543\times 10^{-6}$ | NA

##  5 Conclusion

In this article, sets of analytical formulas for calculation of nine
components of magnetic stiffness of corresponding force arising between two
current-carrying circular filaments arbitrarily oriented in the space have
been derived by using Babic’s method and the method of mutual inductance
(Kalantarov-Zeitlin’s method). Formulas are presented through integral
expressions, whose kernel function is expressed in terms of the elliptic
integrals of the first and second kinds. Also, the additional set of
expressions for calculation of components of magnetic stiffness by means of
differentiation of Grover’s formula with respect to appropriate coordinates
has been obatined. Grover’s method provides the most simplest approach for
calculation of magnetic stiffness, however the calculation is constrained by
four components only, namely, two diaganal and two non-diaganal components
instead of nine ones. Also, the GM suffers from singular cases shown, for
instance, in Examples 4, 5 and 12, which limit the applicability of the
method. In opposite to the GM, the set of formulas (15)-(23), (54)-(57) and
(65)-(68) is deduced by BM and MIM, respectively, is universally applicable
for calculation of the magnetic stiffness and covers all possible arrangements
between two current-carrying circular filaments . The derived sets of formulas
were mutually validated and results of calculation of components of magnetic
stiffness agree well to each other.

The set of formulas (54)-(57) and (65)-(68) obtained by means of MIM and the
set of expressions (83)-(85) obtained by means of GM were programmed by using
the Matlab language. The Matlab files with the implemented formulas are
available as supplementary materials to this article.

## Acknowledgment

Kirill Poletkin acknowledges with thanks the support from German Research
Foundation (Grant KO 1883/37-1) under the priority programme SPP 2206.

##  Appendix A Stiffness calculation. Grover’s method (GM) [38, page 207, Eq.
(179)]

Figure 12: Geometrical scheme of current-carrying circular filaments
arbitrarily oriented in the space: the Grover notations.

According to Grover’s notations, the linear misalignment of the centre of the
secondary circle is characterised by two parameters, namely, $d=z_{c}$ and
$\rho=\sqrt{x_{c}^{2}+y_{c}^{2}}$ as shown in Figure 12. Besides that the
angular misalignment is defined in accordance with the first manner as shown
in Fig. 2, but keeping the original Grover’s notation the angle, $\eta$, is
replaced by $\psi$. In absence of the angular misalignment, the CF $xyz$
assigned to the secondary circle is oriented in the following way. The
$z$-axis is directed upward along the $d$-line, while the $y$-axis is parallel
to the $\rho$-line and directed in continuation of the $\rho$-line. Then
adopting the above considered notations, Grover’s formula for calculation of
mutual inductance between two circular filaments can be written as

$M=\frac{\mu_{0}\sqrt{R_{p}R_{s}}}{2\pi}\int_{0}^{2\pi}U\cdot\Psi(k)d\varphi,$
(69)

where

$U=U({\gamma},\theta,\psi)=\frac{R({\gamma},\theta,\psi)}{{V}^{1.5}}=\frac{\cos\theta-\gamma(\cos\psi\cos\varphi-\sin\psi\cos\theta\sin\varphi)}{{V}^{1.5}},$
(70)
$V=V({\gamma},\theta,\psi)=\sqrt{1-{\cos(\varphi)}^{2}{\sin(\theta)}^{2}+2\gamma(\sin\psi\sin\varphi-\cos\varphi\cos\psi\cos\theta)+\gamma^{2}},$
(71)
$\Psi(k)=\frac{2}{k}\left[\left(1-\frac{k^{2}}{2}\right)K(k)-E(k)\right],$
(72) $\begin{array}[]{l}{\displaystyle
k^{2}=k^{2}({\gamma},\Delta,\theta,\psi)=\frac{4\alpha{V}}{(\alpha{V}+1)^{2}+z^{2}},}\\\
{\displaystyle\alpha=R_{s}/R_{p},\;\Delta=d/R_{p},\;\gamma=\rho/R_{s}},\;z={\Delta}-{\alpha}\sin\theta\cos\varphi.\end{array}$
(73)

The kernel of formula (69) is

$\mathrm{Kr}=U\cdot\Psi(k).$ (74)

Accounting for (70), (71), (72) and (73), the second $\rho$-derivative of the
kernel has a similar structure to Eq. (59) and becomes as follows

$\begin{array}[]{c}{\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial\rho^{2}}=\frac{\partial^{2}\mathrm{Kr}}{\partial\gamma^{2}}\frac{1}{R_{s}^{2}}=\frac{1}{R_{s}^{2}}\cdot\left[\frac{\partial^{2}U}{\partial\gamma^{2}}\cdot\Psi(k)+2\frac{\partial
U}{\partial\gamma}\cdot\frac{d\Psi(k)}{dk}\cdot\frac{\partial
k}{\partial\gamma}\right.}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Psi(k)}{dk^{2}}\cdot\left(\frac{\partial
k}{\partial\gamma}\right)^{2}+\frac{d\Psi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial\gamma^{2}}\right)\right],}\end{array}$
(75)

where

$\begin{array}[]{l}{\displaystyle\frac{\partial
U}{\partial\gamma}=\frac{J}{V^{2.5}}=\left({\displaystyle\frac{\partial
R}{\partial\gamma}}\cdot V-1.5\cdot R\cdot\frac{\partial
V}{\partial\gamma}\right)\bigg{/}{{V}^{2.5}},}\\\
{\displaystyle\frac{\partial^{2}U}{\partial\gamma^{2}}=\frac{{\displaystyle\frac{\partial
J}{\partial\gamma}\cdot V-2.5\cdot J\cdot\frac{\partial
V}{\partial\gamma}}}{{{V}^{3.5}}},}\\\ {\displaystyle\frac{\partial
J}{\partial\gamma}={\displaystyle-0.5\frac{\partial
R}{\partial\gamma}\cdot\frac{\partial V}{\partial\gamma}}-1.5\cdot
R\cdot\frac{\partial^{2}V}{\partial\gamma^{2}}},\\\
{\displaystyle\frac{\partial
R}{\partial\gamma}=-\cos\psi\cos\varphi+\sin\psi\cos\theta\sin\varphi},\\\
{\displaystyle\frac{\partial
V}{\partial\gamma}=\frac{\sin\psi\sin\varphi-\cos\varphi\cos\psi\cos\theta+\gamma}{V}},\\\
{\displaystyle\frac{\partial^{2}V}{\partial\gamma^{2}}=\frac{V-\left(\sin\psi\sin\varphi-\cos\varphi\cos\psi\cos\theta+\gamma\right){\displaystyle\frac{\partial
V}{\partial\gamma}}}{V^{2}}},\end{array}$ (76)
$\begin{array}[]{l}{\displaystyle\frac{\partial
k}{\partial\gamma}=\frac{G}{H}\cdot\alpha\frac{\partial
V}{\partial\gamma},}\\\
{\displaystyle\frac{\partial^{2}k}{\partial\gamma^{2}}=\frac{{\displaystyle\frac{\partial
G}{\partial\gamma}H-A\frac{\partial
H}{\partial\gamma}}}{H^{2}}\cdot\alpha\frac{\partial
V}{\partial\gamma}+\frac{G}{H}\cdot\alpha\frac{\partial^{2}V}{\partial\gamma^{2}},}\\\
{\displaystyle G=2/k-k(\alpha{V}+1),\;H=(\alpha{V}+1)^{2}+{z}^{2},}\\\
{\displaystyle\frac{\partial
G}{\partial\gamma}=-\left[2/k^{2}+\alpha{V}+1\right]\frac{\partial
k}{\partial\gamma}-k\cdot\alpha\frac{\partial
V}{\partial\gamma},\frac{\partial H}{\partial\gamma}=2(\alpha
V+1)\alpha\frac{\partial V}{\partial\gamma}.}\end{array}$ (77)

Note that the first and second derivatives of function $\Psi$ with respect to
$k$ are equal to ${\displaystyle 2\frac{d\Phi}{dk}}$ and ${\displaystyle
2\frac{d^{2}\Phi}{dk^{2}}}$, respectively, where the first and second
derivatives of $\Phi$ with respect to $k$ are defined by Eqs. (43).

The second $d$-derivative of the kernel is

${\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial
d^{2}}=\frac{\partial^{2}\mathrm{Kr}}{\partial\Delta^{2}}\frac{1}{R_{p}^{2}}=\frac{1}{R_{p}^{2}}\cdot
U\left[\frac{d^{2}\Psi(k)}{dk^{2}}\cdot\left(\frac{\partial
k}{\partial\Delta}\right)^{2}+\frac{d\Psi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial\Delta^{2}}\right],}$
(78)

where

$\begin{array}[]{l}{\displaystyle\frac{\partial
k}{\partial\Delta}=\frac{G}{H},\;\frac{\partial^{2}k}{\partial\Delta^{2}}=\frac{{\displaystyle\frac{\partial
G}{\partial\Delta}H-A\frac{\partial H}{\partial\Delta}}}{H^{2}},}\\\
{\displaystyle G=-k\cdot z,\;H=(\alpha{V}+1)^{2}+{z}^{2},}\\\
{\displaystyle\frac{\partial G}{\partial\Delta}=-\frac{\partial
k}{\partial\Delta}\cdot z-k,\frac{\partial H}{\partial\Delta}=2z.}\end{array}$
(79)

Accounting for the property

${\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial\rho\partial
d}=\frac{\partial^{2}\mathrm{Kr}}{\partial d\partial\rho},}$ (80)

the second derivative of the kernel with respect of variables $d$ and $\rho$
can be written as

$\begin{array}[]{c}{\displaystyle\frac{\partial^{2}\mathrm{Kr}}{\partial\rho\partial
d}=\frac{\partial^{2}\mathrm{Kr}}{\partial\gamma\partial\Delta}\frac{1}{R_{s}R_{p}}=\frac{1}{R_{s}R_{p}}\cdot\left[\frac{\partial
U}{\partial\gamma}\cdot\frac{d\Psi(k)}{dk}\cdot\frac{\partial
k}{\partial\Delta}\right.}\\\
{\displaystyle\left.+U\left(\frac{d^{2}\Psi(k)}{dk^{2}}\cdot\frac{\partial
k}{\partial\gamma}\cdot\frac{\partial
k}{\partial\Delta}+\frac{d\Psi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial\gamma\partial\Delta}\right)\right],}\end{array}$
(81)

where

$\begin{array}[]{l}{\displaystyle\frac{\partial^{2}k}{\partial\gamma\partial\Delta}=\frac{{\displaystyle\frac{\partial
G}{\partial\gamma}H-G\frac{\partial H}{\partial\gamma}}}{H^{2}},}\\\
{\displaystyle G=-k\cdot z,\;H=(\alpha{V}+1)^{2}+{z}^{2},}\\\
{\displaystyle\frac{\partial G}{\partial\gamma}=-\frac{\partial
k}{\partial\gamma}\cdot z,\frac{\partial H}{\partial\gamma}=2(\alpha
V+1)\alpha\frac{\partial V}{\partial\gamma}.}\end{array}$ (82)

The first derivatives of functions $U$ and $k$ with respect to variables
$\gamma$ and $\Delta$ in the above equation are the same as determined in Eqs
(76), (77) and (79).

Using the derivatives of the kernel obtained above, we can gain the second
derivatives of Grover’s formula of mutual inductance with respect to the
appropriate coordinates and write corresponding formulas for the calculation
of magnetic stiffness. Taking into account Eqs (75), (78) and (81), we can
write

$\begin{array}[]{l}{\displaystyle
S_{\rho\rho}=-\frac{\mu_{0}I_{p}I_{s}}{2\pi}\frac{\sqrt{R_{p}}}{R_{s}^{1.5}}\int_{0}^{2\pi}\frac{\partial^{2}U}{\partial\gamma^{2}}\cdot\Psi(k)+2\frac{\partial
U}{\partial\gamma}\cdot\frac{d\Psi(k)}{dk}\cdot\frac{\partial
k}{\partial\gamma}}\\\
{\displaystyle+U\left(\frac{d^{2}\Psi(k)}{dk^{2}}\cdot\left(\frac{\partial
k}{\partial\gamma}\right)^{2}+\frac{d\Psi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial\gamma^{2}}\right)\;\;d\varphi,}\end{array}$
(83) $\begin{array}[]{l}{\displaystyle
S_{dd}=-\frac{\mu_{0}I_{p}I_{s}}{2\pi}\frac{\sqrt{R_{s}}}{R_{p}^{1.5}}\int_{0}^{2\pi}U\left[\frac{d^{2}\Psi(k)}{dk^{2}}\left(\frac{\partial
k}{\partial\Delta}\right)^{2}+\frac{d\Psi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial\Delta^{2}}\right]\;d\varphi,}\end{array}$
(84) $\begin{array}[]{l}{\displaystyle S_{\rho
d}=S_{d\rho}=-\frac{\mu_{0}}{2\pi}\frac{I_{p}I_{s}}{\sqrt{R_{p}R_{s}}}\int_{0}^{2\pi}\frac{\partial
U}{\partial\gamma}\cdot\frac{d\Psi(k)}{dk}\cdot\frac{\partial
k}{\partial\Delta}}\\\
{\displaystyle+U\left(\frac{d^{2}\Psi(k)}{dk^{2}}\cdot\frac{\partial
k}{\partial\gamma}\cdot\frac{\partial
k}{\partial\Delta}+\frac{d\Psi(k)}{dk}\cdot\frac{\partial^{2}k}{\partial\gamma\partial\Delta}\right)\;\;d\varphi.}\end{array}$
(85)

## References

  * [1] E. Okress, D. Wroughton, G. Comenetz, P. Brace, J. Kelly, Electromagnetic levitation of solid and molten metals, Journal of Applied Physics 23 (5) (1952) 545–552.
  * [2] B. K. Narukullapati, T. K. Bhattacharya, Determination of a stable lateral region of a floating disc–Mathematical analysis and FEM simulation, Alexandria Engineering Journal 60 (3) (2021) 3107–3118.
  * [3] E. Paredes, G. B. Galán, J. A. Veira, J. Mosqueira, L. Romaní, G. Domarco, Measurements of the currents induced in coaxially confronted superconducting rings, Engineering Research Express 3 (1) (2021) 015010.
  * [4] R. Ravaud, G. Lemarquand, V. Lemarquand, S. Babic, C. Akyel, Mutual inductance and force exerted between thick coils, Progress In Electromagnetics Research 102 (2010) 367–380.
  * [5] S. I. Babic, C. Akyel, Magnetic force calculation between thin coaxial circular coils in air, IEEE Transactions on Magnetics 44 (4) (2008) 445–452.
  * [6] S. I. Babic, C. Akyel, Torque calculation between circular coils with inclined axes in air, International Journal of Numerical Modelling: Electronic Networks, Devices and Fields 24 (3) (2011) 230–243.
  * [7] S. I. Babic, C. Akyel, Y. Ren, W. Chen, Magnetic force calculation between circular coils of rectangular cross section with parallel axes for superconducting magnet, Progress In Electromagnetics Research 37 (2012) 275–288.  doi:10.2528/PIERB11110508.
URL http://dx.doi.org/10.2528/PIERB11110508

  * [8] U.-M. Jow, M. Ghovanloo, Design and optimization of printed spiral coils for efficient transcutaneous inductive power transmission, IEEE Transactions on biomedical circuits and systems 1 (3) (2007) 193–202.
  * [9] Y. P. Su, X. Liu, S. Y. R. Hui, Mutual inductance calculation of movable planar coils on parallel surfaces, IEEE Transactions on Power Electronics 24 (4) (2009) 1115–1123.  doi:10.1109/TPEL.2008.2009757.
  * [10] S. Y. Chu, A. T. Avestruz, Transfer-power measurement: A non-contact method for fair and accurate metering of wireless power transfer in electric vehicles, in: 2017 IEEE 18th Workshop on Control and Modeling for Power Electronics (COMPEL), 2017, pp. 1–8.  doi:10.1109/COMPEL.2017.8013344.
  * [11] A. Shiri, A. Shoulaie, A new methodology for magnetic force calculations between planar spiral coils, Progress In Electromagnetics Research 95 (2009) 39–57.
  * [12] R. Ravaud, G. Lemarquand, V. Lemarquand, Force and stiffness of passive magnetic bearings using permanent magnets. part 1: Axial magnetization, IEEE transactions on magnetics 45 (7) (2009) 2996.
  * [13] S. Obata, A muscle motion solenoid actuator, Electrical Engineering in Japan 184 (2) (2013) 10–19.
  * [14] K. Poletkin, A. Chernomorsky, C. Shearwood, U. Wallrabe, A qualitative analysis of designs of micromachined electromagnetic inductive contactless suspension, International Journal of Mechanical Sciences 82 (2014) 110–121.  doi:10.1016/j.ijmecsci.2014.03.013.
URL http://authors.elsevier.com/sd/article/S0020740314000897

  * [15] Z. Lu, K. Poletkin, B. den Hartogh, U. Wallrabe, V. Badilita, 3D micro-machined inductive contactless suspension: Testing and modeling, Sensors and Actuators A Physical 220 (2014) 134–143.  doi:10.1016/j.sna.2014.09.017.
URL http://dx.doi.org/10.1016/j.sna.2014.09.017

  * [16] K. Poletkin, Z. Lu, U. Wallrabe, J. Korvink, V. Badilita, Stable dynamics of micro-machined inductive contactless suspensions, International Journal of Mechanical Sciences 131-132 (2017) 753 – 766.  doi:https://doi.org/10.1016/j.ijmecsci.2017.08.016.
URL http://www.sciencedirect.com/science/article/pii/S0020740316306555

  * [17] K. V. Poletkin, A. I. Chernomorsky, C. Shearwood, A proposal for micromachined accelerometer, base on a contactless suspension with zero spring constant, IEEE Sensors J. 12 (07) (2012) 2407–2413.  doi:10.1109/JSEN.2012.2188831.
  * [18] K. V. Poletkin, J. G. Korvink, Modeling a pull-in instability in micro-machined hybrid contactless suspension 7 (1) (2018) 11.
  * [19] K. Poletkin, Static pull-in behavior of hybrid levitation micro-actuators: simulation, modelling and experimental study, IEEE/ASME Transactions on Mechatronics (2020) 753–764doi:10.1109/TMECH.2020.2999516.
URL https://doi.org/10.1109/TMECH.2020.2999516

  * [20] K. Poletkin, On the static pull-in of tilting actuation in electromagnetically levitating hybrid micro-actuator: Theory and experiment, Actuators 10 (10) (2021).  doi:10.3390/act10100256.
URL https://www.mdpi.com/2076-0825/10/10/256

  * [21] T. Theodoulidis, R. J. Ditchburn, Mutual impedance of cylindrical coils at an arbitrary position and orientation above a planar conductor, IEEE Transactions on Magnetics 43 (8) (2007) 3368–3370.  doi:10.1109/TMAG.2007.894559.
  * [22] M. Sawan, S. Hashemi, M. Sehil, F. Awwad, M. Hajj-Hassan, A. Khouas, Multicoils-based inductive links dedicated to power up implantable medical devices: modeling, design and experimental results, Biomedical Microdevices 11 (5) (2009) 1059.  doi:10.1007/s10544-009-9323-7.
URL https://doi.org/10.1007/s10544-009-9323-7

  * [23] S. Kuznetsov, J. K. Guest, Topology optimization of magnetic source distributions for diamagnetic and superconducting levitation, Journal of Magnetism and Magnetic Materials 438 (2017) 60 – 69.  doi:https://doi.org/10.1016/j.jmmm.2017.04.052.
URL http://www.sciencedirect.com/science/article/pii/S0304885316319515

  * [24] D. Hoult, B. Tomanek, Use of mutually inductive coupling in probe design, Concepts in Magnetic Resonance 15 (4) (2002) 262–285.  arXiv:https://onlinelibrary.wiley.com/doi/pdf/10.1002/cmr.10047, doi:10.1002/cmr.10047.
URL https://onlinelibrary.wiley.com/doi/abs/10.1002/cmr.10047

  * [25] N. Spengler, P. T. While, M. V. Meissner, U. Wallrabe, J. G. Korvink, Magnetic lenz lenses improve the limit-of-detection in nuclear magnetic resonance, PLOS ONE 12 (8) (2017) 1–17.  doi:10.1371/journal.pone.0182779.
URL https://doi.org/10.1371/journal.pone.0182779

  * [26] G. D. Angelis, V. Pasku, A. D. Angelis, M. Dionigi, M. Mongiardo, A. Moschitta, P. Carbone, An indoor ac magnetic positioning system, IEEE Transactions on Instrumentation and Measurement 64 (5) (2015) 1267–1275.  doi:10.1109/TIM.2014.2381353.
  * [27] F. Wu, J. Jeon, S. K. Moon, H. J. Choi, H. Son, Voice coil navigation sensor for flexible silicone intubation, IEEE/ASME Transactions on Mechatronics 21 (2) (2016) 851–859.  doi:10.1109/TMECH.2015.2476836.
  * [28] D. Jiao, L. Ni, X. Yu, J. Zhe, A non-contact gap measurement method for narrow and irregular targets using inductive planar sensors with simple calibration, IET Science, Measurement & Technology 15 (9) (2021) 710–718.  arXiv:https://ietresearch.onlinelibrary.wiley.com/doi/pdf/10.1049/smt2.12071, doi:https://doi.org/10.1049/smt2.12071.
URL https://ietresearch.onlinelibrary.wiley.com/doi/abs/10.1049/smt2.12071

  * [29] X. Zhang, C. Quan, Z. Li, Mutual inductance calculation of circular coils for an arbitrary position with electromagnetic shielding in wireless power transfer systems, IEEE Transactions on Transportation Electrification (2021) 1–1doi:10.1109/TTE.2021.3054762.
  * [30] S. Y. Chu, X. Cui, X. Zan, A.-T. Avestruz, Transfer-power measurement using a non-contact method for fair and accurate metering of wireless power transfer in electric vehicles, IEEE Transactions on Power Electronics 37 (2) (2021) 1244–1271.
  * [31] B. Gulbahar, A communication theoretical analysis of multiple-access channel capacity in magneto-inductive wireless networks, IEEE Transactions on Communications 65 (6) (2017) 2594–2607.  doi:10.1109/TCOMM.2017.2669995.
  * [32] K. V. Poletkin, J. G. Korvink, Efficient calculation of the mutual inductance of arbitrarily oriented circular filaments via a generalisation of the Kalantarov-Zeitlin method, Journal of Magnetism and Magnetic Materials 483 (2019) 10–20.  doi:https://doi.org/10.1016/j.jmmm.2019.03.078.
URL https://www.sciencedirect.com/science/article/pii/S0304885318337703

  * [33] S. Babic, F. Sirois, C. Akyel, C. Girardi, Mutual inductance calculation between circular filaments arbitrarily positioned in space: Alternative to Grover’s formula, IEEE Transactions on Magnetics 46 (9) (2010) 3591–3600.  doi:10.1109/TMAG.2010.2047651.
URL https://doi.org/10.1109/TMAG.2010.2047651

  * [34] S. Babic, C. Akyel, Magnetic force between inclined circular filaments placed in any desired position, IEEE transactions on magnetics 48 (1) (2012) 69–80.
  * [35] S. Babic, C. Akyel, Magnetic force between inclined circular loops (Lorentz Approach), Progress In Electromagnetics Research B 38 (2012) 333–349.
  * [36] S. Babic, Vector potential, magnetic field, mutual inductance, magnetic force, torque and stiffness calculation between current-carrying arc segments with inclined axes in air, Physics 3 (4) (2021) 1054–1087.  doi:10.3390/physics3040067.
URL https://www.mdpi.com/2624-8174/3/4/67

  * [37] S. Paul, J. Wright, J. Z. Bird, 3-D steady-state eddy-current damping and stiffness terms for a finite thickness conductive plate, IEEE Transactions on Magnetics 50 (11) (2014) 1–4.  doi:10.1109/TMAG.2014.2327097.
  * [38] F. W. Grover, Inductance calculations : working formulas and tables, special ed. prepared for instrument society of america Edition, Research Triangle Park, N.C. : Instrument Society of America, 1981, reprint. Originally published: New York : Van Nostrand, 1946. With publisher’s comment.
  * [39] K. Poletkin, Calculation of force and torque between two arbitrarily oriented circular filaments using kalantarov-zeitlin’s method (2021).  arXiv:2106.09496.
URL https://arxiv.org/abs/2106.09496