spacing=nonfrench

# Sliding down over a horizontally moving semi-sphere

Roberto A. Lineros1 Departamento de Física, Universidad Católica del Norte,
Avenida Angamos 0610, Casilla 1280, Antofagasta, Chile.
<EMAIL_ADDRESS>
(August 27, 2024)

###### Abstract

We studied the dynamics of an object sliding down on a semi-sphere with radius
$R$. We consider the physical setup where the semi-sphere is free to move over
a flat surface. For simplicity, we assume that all surfaces are friction-less.
We analyze the values for the last contact angle $\theta^{\star}$,
corresponding to the angle when the object and the semi-sphere detach one of
each other. We consider all possible scenarios with different combination of
mass values: $m_{A}$ and $m_{B}$, and the initial velocity of the sliding
object $A$. We found that the last contact angle only depends on the ratio
between the masses, and it is independent of the acceleration of gravity and
semi-sphere’s radius. In addition, we found that the largest possible value of
$\theta^{\star}$ is $48.19^{\circ}$ that coincides with the case of a fixed
semi-sphere. On the opposite case, the minimum value of $\theta^{\star}$ is
$0^{\circ}$ and it occurs then the object on the semi-sphere is extremely
heavy, occurring the detachment as soon as the sliding body touches the semi-
sphere. In addition, we found that if the initial kinetic energy of the
sliding object $A$ is half the value of the potential energy with respect to
the floor. The object detaches at the top of the semi-sphere.

††: Eur. J. Phys.

##  1 Introduction

In courses of Newtonian mechanics for engineers and physics students at
undergraduate level, the concepts behind Newton’s laws are key for
understanding the kinematics of objects under the effects of forces. Indeed,
in courses focuses on engineering, there is preference to solve problem using
only Newtonian mechanics instead other possible approaches like Lagrangian or
Hamiltonian mechanics. No matter the approach, it is troublesome for the
students to understand the interplay between objects in contact due to the
presence of reaction forces in cases where the contact surfaces are not flat.
At the undergraduate course level, the focus on the study of vectorial
mechanics is crucial for dealing with challenging problems like ones involving
many bodies and the use of different vectorial basis at the same time, like
cartesian, cylindrical, and spherical basis.

The problem of an object sliding down on a circular path is an academic
example to teach such concepts [1, 2, 3]. Similarly, the problem of an object
descending on a flat surface with a slope is another common example oriented
to teach relative motion in terms of relative position, velocity, and
acceleration. However, sometimes there are variations on this class of
problems by allowing the inclined plane to be affected by the reaction force
between the object and the inclined plane (see problem 15-98 of [3]). In the
literature, similar problems have been addressed with different approaches
like the use of lagrangian mechanics [4], including friction [5, 6, 7], or the
experimental demonstration [8].

In this manuscript, we consider the situation of a moving semi-sphere, which
is initially at rest, where its movement is due to the reaction force of an
object sliding down on top of it.

The manuscript is organized as follows: In section 2, we present the solution
for the standard case, where the semi-sphere remains still, in addition, we
use the result as benchmark for the moving setup. In section 3, We present the
solution and analysis for the moving semi-sphere. Finally, section 4 are the
conclusions.

##  2 System with fixed semi-sphere.

Figure 1: (Left) Physical situation on the object $A$ sliding down over a
fixed semi-sphere $B$. The semi-sphere $B$ does not present friction. Forces
acting on $A$ are shown in blue. (Right) Vector basis description.

At a first stage, we consider the situation when semi-sphere $B$ remains still
(see figure 1). This is a known problem taught in courses of Newtonian
Mechanics at the university level. The equations of motion for the object $A$
are constructed using the second Newton’s Law and correspond to:

$\displaystyle\sum\vec{F}_{A}=\vec{N}_{A}+\vec{W}_{A}=m_{A}\vec{a}_{A}\,,$ (1)

where $\vec{N}_{A}$ is the reaction force between objects and
$\vec{W}_{A}=-m_{A}g\,\hat{k}$ is the weight with $g$ the gravity’s
acceleration. While the object $A$ is in contact with the semi-sphere $B$, it
moves along the surface following a circular path. The acceleration
$\vec{a}_{A}$ is then obtained by

$\vec{a}_{A}=\vec{\alpha}\times\vec{r}+\vec{\omega}\times\left(\vec{\omega}\times\vec{r}\right)\,,$
(2)

where $\vec{\alpha}$ and $\vec{\omega}$ are the vectors of angular
acceleration and angular velocity, respectively. For simplicity, the movement
of $A$ is 2-dimensional and occurs along the plane defined by the unit vector
$\,\hat{\imath}$ and $\,\hat{k}$. Therefore, we will use a vectorial
cylindrical basis defined by:

$\displaystyle\,\hat{r}$ $\displaystyle=$
$\displaystyle\cos\theta\,\hat{k}+\sin\theta\,\hat{\imath}\,,$ (3)
$\displaystyle\,\hat{\theta}$ $\displaystyle=$
$\displaystyle-\sin\theta\,\hat{k}+\cos\theta\,\hat{\imath}\,,$ (4)

where $\theta$ is the angle between $\,\hat{k}$ and $\,\hat{r}$. Notice that
the 3-dimensional vectorial basis maintains the right-hand rule convention,
such as
$\,\hat{k}\times\,\hat{\imath}=\,\hat{r}\times\,\hat{\theta}=\,\hat{\jmath}$.
In addition, this basis is different to the standard cylindrical basis, so the
reader might be cautioned about this fact.

The acceleration $\vec{a}_{A}$ in the cylindrical basis corresponds to

$\vec{a}_{A}=\alpha R\,\hat{\theta}-\omega^{2}R\,\hat{r}\,,$ (5)

when $\vec{\alpha}=\alpha\,\hat{\jmath}$ and
$\vec{\omega}=\omega\,\hat{\jmath}$.

The equations of motions are in terms of the $\,\hat{r}$ and $\,\hat{\theta}$
components:

$\displaystyle N_{A}-m_{A}g\cos\theta$ $\displaystyle=$ $\displaystyle-
m_{A}\omega^{2}R\,,$ (6) $\displaystyle m_{A}g\sin\theta$ $\displaystyle=$
$\displaystyle m_{A}\alpha R\,.$ (7)

For the setup in figure 1, the angular velocity and acceleration are related
to the angle $\theta$ by:

$\displaystyle\omega$ $\displaystyle=$ $\displaystyle\dot{\theta}\,,$ (8)
$\displaystyle\alpha$ $\displaystyle=$
$\displaystyle\frac{d\dot{\theta}}{d\theta}\dot{\theta}\,.$ (9)

Using the latter expressions, the equations of motion are reduced to a couple
of differential equations:

$\displaystyle\frac{g}{R}\cos\theta-\frac{N_{A}}{m_{A}R}$ $\displaystyle=$
$\displaystyle\dot{\theta}^{2}\,,$ (10) $\displaystyle\frac{g}{R}\sin\theta$
$\displaystyle=$ $\displaystyle\frac{d\dot{\theta}}{d\theta}\dot{\theta}\,.$
(11)

These equations are further simplified via the substitution:
$f(\theta)=\dot{\theta}^{2}$, $f^{\prime}(\theta)=df/d\theta=2\ddot{\theta}$,
and $\kappa=g/R$; obtaining:

$\displaystyle\kappa\cos\theta-\frac{N_{A}}{m_{A}R}$ $\displaystyle=$
$\displaystyle f(\theta)\,,$ (12) $\displaystyle 2\kappa\sin\theta$
$\displaystyle=$ $\displaystyle f^{\prime}(\theta)\,.$ (13)

Since the function $f(\theta)$ corresponds to the square of the angular
velocity, the kinetic energy of $A$ corresponds to

$\displaystyle T=\frac{1}{2}m_{A}R^{2}f(\theta)\,,$ (14)

providing that the total mechanical energy is:

$\displaystyle E=\frac{1}{2}m_{A}R^{2}f(\theta)+m_{A}gR\cos\theta\,.$ (15)

Depending on the program of content in a vectorial mechanics course, the
concept of mechanical energy may not be seen. However, besides the resolution
using equation of motion, this problem can be also solve using in addition
conservation of energy.

###  2.1 Solving the equations of motion

These equations are simply solved by integrating over the $\theta$ angle in
equation 13 and after by replacing in equation 12. When considering the
initial conditions:

$\displaystyle\theta(t=0)$ $\displaystyle=$ $\displaystyle\theta_{0}=0\,,$
(16) $\displaystyle\dot{\theta}^{2}(t=0)$ $\displaystyle=$ $\displaystyle
f(\theta_{0})=2\kappa\epsilon\,,$ (17)

then the reaction force and the angular velocity squared are:

$\displaystyle N_{A}$
$\displaystyle=m_{A}g\left(3\cos\theta-2(1+\epsilon)\right)\,,$ (18)
$\displaystyle f(\theta)$
$\displaystyle=\dot{\theta}^{2}=2\kappa\left(1+\epsilon-\cos\theta\right)\,.$
(19)

We introduce the parameter $\epsilon$ to further simplify the solutions,
however, it is related with the initial kinetic energy:

$\displaystyle\epsilon=\frac{T_{0}}{m_{A}gR}\,,$ (20)

and it can be interpreted as the ration between the initial kinetic energy and
the initial potential energy. It is important to remark that the initial
position, $\theta(t=0)=0$, is an unstable equilibrium point. In order to break
the symmetrical evolution of sliding down to any side of the semi-sphere, it
is important to indicate an initial direction of movement by means of the
value of $\epsilon$.

The equations of motion are valid only for the regimen when $N_{A}\geq 0$ and
describe the object $A$ moving over the semi-sphere. The case $N_{A}=0$ sets
the last contact angle $\theta^{\star}$ that in this case corresponds to:

$\cos\theta^{\star}=\frac{2}{3}\left(1+\epsilon\right)$ (21)

This general solution allows us to set a maximum value of $\epsilon$ for which
the object $A$ detaches at the beginning of the movement
($\cos\theta^{\star}=1$):

$\epsilon_{\rm max}=\frac{1}{2}\,,$ (22)

which sets the maximum angular velocity to be:

$\dot{\theta}(t=0)=\sqrt{\kappa}\,.$ (23)

This sets the maximum initial kinetic energy to be:

$\displaystyle T_{0}^{\rm max}=\frac{1}{2}m_{A}gR\,,$ (24)

which corresponds to half of the initial potential energy.

In the limit, $\epsilon\rightarrow 0$, which means that movement starts at
rest, the last contact angle becomes independent of $R$ and $g$, and gives the
value:

$\cos\theta^{\star}=\frac{2}{3}\,\rightarrow\,\theta^{\star}\simeq
48.19^{\circ}\,.$ (25)

This value is a known result in the frictionless scenario, and the problem
corresponds to an exercise in books [1, 3, 2]. A more complete case, that
includes friction, can be seen in [5]. This result represents a benchmark
value to compare with the following case of a moving semi-sphere.

##  3 System with a moving semi-sphere

In this part, we allow the semi-sphere to freely move over a friction-less
surface. Therefore, the aim is to calculate the value of the last contact
angle $\theta^{\star}$ under these new conditions. Figure 2 shows the physical
setup.

Figure 2: (Left) Physical situation on the object $A$ sliding down over a
moving semi-sphere $B$. The semi-sphere $B$ does not present friction with
both: object $A$, and the flat surface. Forces acting on $A$ are shown in
blue. Forces acting on $B$ are shown in dark pink. Accelerations are displayed
in orange. (Right) Unit vector basis description.

###  3.1 Equations of motions

Similarly to the previously presented case, the object $A$ is affected by its
weight and the reaction force with the semi-sphere. Therefore, the equation of
motion for $A$ is:

$\sum\vec{F}_{A}=\vec{N}_{A}+\vec{W}_{A}=m_{A}\vec{a}_{A}\,.$ (26)

On the other hand, the semi-sphere $B$ is now free to move and its equation of
motion is given by:

$\sum\vec{F}_{B}=\vec{N}_{B}-\vec{N}_{A}+\vec{W}_{B}=m_{B}\vec{a}_{B}\,,$ (27)

where $\vec{N}_{B}$ is the reaction force between the floor and the semi-
sphere and $\vec{W}_{B}=-m_{B}g\,\hat{k}$ is its weight.

The movement of the semi-sphere $B$ is only horizontal, therefore, the
acceleration corresponds to $\vec{a}_{B}=a_{Bx}\,\hat{\imath}$ as well as its
velocity $\vec{v}_{B}=v_{Bx}\,\hat{\imath}$. It is also important to remark
that $\vec{a}_{B}=\vec{a}_{O}$ because the semi-sphere $B$ also non-rotating.

On the other hand, the acceleration of the object $A$ with respect to the
floor i.e. the rest frame is computed through the acceleration of $B$ and the
relative acceleration between objects:

$\displaystyle\vec{a}_{A}$ $\displaystyle=$
$\displaystyle\vec{a}_{B}+\vec{a}_{A/B}\,,$ (28) $\displaystyle\vec{a}_{A}$
$\displaystyle=$
$\displaystyle\vec{a}_{B}+\vec{\alpha}\times\vec{r}+\vec{\omega}\times\left(\vec{\omega}\times\vec{r}\right)\,.$
(29)

While the object $A$ is in contact with the semi-sphere $B$, the relative
movement of $A$ with respect to $B$ is a circular motion. Therefore, we can
write the acceleration of $A$ in term of the cylindrical components:

$\displaystyle\vec{a}_{A}$ $\displaystyle=a_{Bx}\,\hat{\imath}+\alpha
R\,\hat{\theta}-\omega^{2}R\,\hat{r}\,.$ (30)

Up to this point, the acceleration of $A$ in the cylindrical basis is:

$\vec{a}_{A}=\left(a_{Bx}\sin\theta-\omega^{2}R\right)\,\hat{r}+\left(a_{Bx}\cos\theta+\alpha
R\right)\,\hat{\theta}\,,$ (31)

which allow us to get the full set of equations of motion. 2 equations for the
object $A$:

$\displaystyle N_{A}-m_{A}g\cos\theta$ $\displaystyle=$ $\displaystyle
m_{A}\left(a_{Bx}\sin\theta-\omega^{2}R\right)\,,$ (32) $\displaystyle
m_{A}g\sin\theta$ $\displaystyle=$ $\displaystyle
m_{A}\left(a_{Bx}\cos\theta+\alpha R\right)\,,$ (33)

and 2 more for the semi-sphere $B$:

$\displaystyle-N_{A}\sin\theta$ $\displaystyle=$ $\displaystyle
m_{B}a_{Bx}\,,$ (34) $\displaystyle N_{B}-N_{A}\cos\theta-m_{B}g$
$\displaystyle=$ $\displaystyle 0\,.$ (35)

After inspection, the last 2 equations lead to:

$\displaystyle a_{Bx}$ $\displaystyle=$
$\displaystyle-\frac{N_{A}\sin\theta}{m_{B}}\,,$ (36) $\displaystyle N_{B}$
$\displaystyle=$ $\displaystyle m_{B}g+N_{A}\cos\theta\,,$ (37)

and these can be used to reduce the full set of equations of motion to only
two equations:

$\displaystyle\frac{g}{R}\cos\theta-\frac{N_{A}}{m_{A}R}\left(1+\frac{m_{A}\sin^{2}\theta}{m_{B}}\right)$
$\displaystyle=$ $\displaystyle f(\theta)\,,$ (38)
$\displaystyle\frac{2g}{R}\sin\theta+\frac{2N_{A}}{m_{B}R}\sin\theta\cos\theta$
$\displaystyle=$ $\displaystyle f^{\prime}(\theta)\,,$ (39)

where $f(\theta)=\dot{\theta}^{2}$ and $f^{\prime}(\theta)=2\ddot{\theta}$.
Notice, we use the same relations for the angular acceleration and velocity
with respect to the angle $\theta$: $\alpha=\ddot{\theta}$, and
$\omega=\dot{\theta}$.

The function $f(\theta)$ is related to the kinetic energy of $A$ only if the
semi-sphere $B$ is at rest. In a general case, the kinetic energy of $A$
depends on the velocity of the semi-sphere and the relative velocity between
objects:

$\displaystyle
T_{A}=\frac{1}{2}m_{A}v_{A}^{2}=\frac{1}{2}m_{A}\left(\vec{v}_{B}+\vec{v}_{A/B}\right)^{2}=\frac{1}{2}m_{A}\left(\vec{v}_{B}+\vec{\omega}\times\vec{r}\right)^{2}\,,$
(40)

where in terms of $f(\theta)$, it corresponds to:

$\displaystyle T_{A}=\frac{1}{2}m_{A}\left(v_{B}^{2}+2v_{B}\cos\theta
R\sqrt{f(\theta)}+R^{2}f(\theta)\right)\,,$ (41)

and with it, we can obtain the total mechanical energy of the system:

$\displaystyle
E=\frac{m_{B}}{2}v_{B}^{2}+\frac{m_{A}}{2}\left(v_{B}^{2}+2v_{B}\cos\theta
R\sqrt{f(\theta)}+R^{2}f(\theta)\right)+m_{A}gR\cos\theta\,.$ (42)

###  3.2 Solving the equations of motion

The solution of the system of equation cannot be performed as in the fixed
semi-sphere case because the reaction force $N_{A}$ is present in both
equations, and it depends on the angle. Nevertheless, the reaction force
$N_{A}$ can be isolated from the equation 38 and be written in terms of the
dynamical variables:

$N_{A}(\theta)=m_{A}R\,\frac{\kappa\cos\theta-f(\theta)}{1+\beta\sin^{2}\theta}\,,$
(43)

where $\beta=m_{A}/m_{B}$ is the ratio between the masses and $\kappa=g/R$ is
the ratio between the acceleration of gravity and the radius of the semi-
sphere. Here, the reaction force depends on the square of the angular velocity
encoded in $f(\theta)$.

After removing the explicit dependence of $N_{A}$ in the equation 39, we get
the differential equation for $f(\theta)$:

$\displaystyle\frac{1+\beta\sin^{2}\theta}{2\sin\theta}\,f^{\prime}(\theta)+\beta\cos\theta\,f(\theta)-\kappa(1+\beta)=0\,.$
(44)

This differential equation is analytically solvable and has the solution:

$\displaystyle
f(\theta)=\frac{2\kappa(1+\beta)(1-\cos\theta)+\epsilon}{1+\beta\sin^{2}\theta}\,,$
(45)

when the initial conditions: $\theta(t=0)=0$ and
$f(\theta_{0})=2\kappa\epsilon$ are included. Here, the parameter $\epsilon$
is related to the kinetic energy of the object $A$ in the same way as in the
previous section (equation 20). This is because the semi-sphere $B$ is
initially at rest.

Let’s remark that this solution is valid for, $N_{A}(\theta)\geq 0$ and it
holds for angles $0\leq\theta\leq\theta^{\star}\leq\pi/2$ where
$\theta^{\star}$ corresponds to the last contact angle.

###  3.3 Finding the last contact angle

The angle $\theta^{\star}$ can be obtained by solving the following equation:

$N_{A}(\theta^{\star})=m_{A}R\,\frac{\kappa\cos\theta^{\star}-f(\theta^{\star})}{1+\beta\sin^{2}\theta^{\star}}=0\,,$
(46)

where $f(\theta^{\star})$ is the solution of the differential equation
evaluated at $\theta^{\star}$ (equation 45). This leads to the following
equation to solve for $\theta^{\star}$:

$\displaystyle\kappa\cos\theta^{\star}=f(\theta^{\star})=\frac{2\kappa(1+\beta)(1-\cos\theta^{\star})+\epsilon}{1+\beta\sin^{2}\theta^{\star}}\,.$
(47)

This equation can be written as a depressed cubic equation for
$\xi=\cos\theta^{\star}$ as follows:

$H(\xi)=\sin^{2}{\left(\frac{\pi}{2}\tau\right)}\,\xi^{3}-3\xi+2+2\epsilon\cos^{2}{\left(\frac{\pi}{2}\tau\right)}=0\,,$
(48)

where the $H(\xi)$ function is just a reparametrization of equation 47 in
terms of adimensional parameters.

We introduce the $\tau$-parameter ranging $0\leq\tau\leq 1$ that simplifies
the expressions better than the mass ratio $m_{A}/m_{B}$ such as:

$\beta=\frac{m_{A}}{m_{B}}=\tan^{2}\left(\frac{\pi}{2}\tau\right)\,.$ (49)

At the extremes values of $\tau$, we obtain that $\tau\rightarrow 0$
corresponds to $\beta\rightarrow 0$ meaning the case of a heavy semi-sphere
which is equivalent to the fixed semi-sphere case. Similarly, when
$\tau\rightarrow 1$ corresponds to $\beta\rightarrow\infty$ meaning the
scenario where object $A$ is extreamly heavy with respect to the semi-sphere.
In figure 3, we present the function $H(\xi)$ for various values of $\tau$.

Figure 3:  $H(\xi)$ versus $\xi$ for combination of the parameters:
$\tau=0,1/2,1$ and $\epsilon=0$ (solid lines), $1/3$ (dashed lines). The green
dashed line correspond to $\epsilon=1/2$ and $\tau=0$. The intersection of
each curve $H(\xi)$ with the horizontal black line ($H(\xi)=0$) produces the
solution where $\xi=\cos\theta^{\star}$.

The analysis of the equation 48 reveals the maximum value that $\epsilon$ can
take. To find this, we evaluate the function $H(\xi)$ in $\xi=1$
($\theta^{\star}=0$):

$H(1)=(2\epsilon_{\rm max}-1)\cos^{2}{\left(\frac{\pi}{2}\tau\right)}=0\,,$
(50)

obtaining that for the range $0\leq\tau<1$, the maximum value is:

$\epsilon_{\rm max}=\frac{1}{2}\,,$ (51)

which is same the same limit obtained in the fixed semi-sphere case. This
means that the kinetic energy of $A$ can be at most half of the initial
potential energy in order to get the object $A$ sliding down the semi-sphere
$B$. This is regardless of the mass ratio $\beta$.

As stated before, the limit of fixed semi-sphere is reached when
$\tau\rightarrow 0$ and it corresponds to the mass limit: $m_{A}\rightarrow 0$
or $m_{B}\rightarrow\infty$. In this limit, equation 48 corresponds to:

$-3\xi+2+2\epsilon=0\,,$ (52)

with solution :

$\displaystyle\xi=\cos\theta^{\star}=\frac{2}{3}(1+\epsilon)\,,$ (53)

corresponding exactly to equation 21. This case agrees with the solution for
the fixed semi-sphere described in section 2, and it indicates the limit
$\beta$ is equivalent to restrict the movement of $B$ to be fixed in a point.

The other limit, $\tau\rightarrow 1$, gives the equation:

$\xi^{3}-3\xi+2=(\xi+2)(\xi-1)^{2}=0\,,$ (54)

with 3 real solutions: $\xi_{1}=1$, $\xi_{2}=1$, and $\xi_{3}=-2$. The
solution with physical meaning are $\xi_{1,2}$, both correspond to an angle
$\theta^{\star}=0$. This means when $m_{A}\rightarrow\infty$ the angle of last
contact between $A$ and $B$ is at the top of the semi-sphere and happens at
the beginning of the movement. In this case, the semi-sphere moves fast enough
to not be in contact with the object $A$.

For the intermediate cases, $0<\tau<1$ and $\epsilon=0$, the solutions for
equation 48 corresponds to:

$\displaystyle\xi_{1}$ $\displaystyle=$
$\displaystyle\frac{2}{1+2\cos{\left(\frac{\pi}{3}\tau\right)}}\,,$ (55)
$\displaystyle\xi_{2}$ $\displaystyle=$
$\displaystyle\frac{\sqrt{3}\cos{\left(\frac{\pi}{6}\tau\right)}-\sin{\left(\frac{\pi}{6}\tau\right)}}{\sin{\left(\frac{\pi}{2}\tau\right)}}\,,$
(56) $\displaystyle\xi_{3}$ $\displaystyle=$
$\displaystyle-\frac{\sqrt{3}\cos{\left(\frac{\pi}{6}\tau\right)}+\sin{\left(\frac{\pi}{6}\tau\right)}}{\sin{\left(\frac{\pi}{2}\tau\right)}}\,,$
(57)

which are obtained by analytically solving the depressed cubic equation.
However, in order to get the real values the solutions of equation 48 need to
be rephased by $e^{i\pi/3}$ when the equation is solved via the Vieta’s
substitution [9]. From the 3 roots, $\xi_{1}$ has a physical meaning:
$\xi_{1}=\cos\theta^{\star}$ while the other 2 roots give values outside the
physical range $0\leq\xi\leq 1$.

The approach to analytically solve $H(\xi)=0$ might be hard to perform by
students. An alternative approach and much easier might be by using a
numerical code, for example, written in python (See appendix A).

The dependence of the last contact angle $\theta^{\star}$ in terms of $\tau$
and $\epsilon$ is shown in figure 4. We observe that the solution of the cubic
equation include the extreme limits, like the fixed semi-sphere ($\tau=0$) for
different values of $\epsilon$.

In addition, the value of $\theta^{\star}$ indicates that the largest possible
value for any configuration of masses corresponds to the fixed semi-sphere
case and the lowest corresponds when $m_{A}\rightarrow\infty$ or
$m_{B}\rightarrow 0$ producing an extreme situation where the semi-sphere $B$
and the object $A$ get immediately detached upon the first contact.

Moreover, the dependence of $\theta^{\star}$ in terms of the value of
$\epsilon$ gives a interesting situation that implies that the object $A$ can
only remain in contact with the semi-sphere if its initial kinetic energy is,
at most, half of the initial potential energy. This can be appreciated in
figure 4 where the blue line gives the value $\theta^{\star}=0$ for all
possible values of $\tau$.

Figure 4:  Last contact angle $\theta^{\star}$ versus $\tau$. Red line shown
the functional dependence in $\tau$ for $\epsilon=0$, the brown line
corresponds to the case when $\epsilon=1/3$, and blue line for $\epsilon=1/2$.
The gray dashed line indicate $\theta^{\star}=48.19^{\circ}$ which corresponds
to the solution of the fixed-sphere case.

##  4 Conclusions

We present the analytical solution of finding the last contact angle for the
problem of an object sliding down on a semi-sphere of radius $R$ where the
semi-sphere is on a friction-less surface. The approach used to solve the
problem is in terms of Newtonian mechanics and concepts of vectorial
mechanics, which are topics familiar to undergraduate students of engineering
and bachelor in physics. The key effect to consider is the reaction forces
between the object and the semi-sphere. This force provokes the semi-sphere to
move horizontally and the object $A$ to descend, keeping contact with the
semi-sphere. If the velocity of the object $A$ is larger enough that the
reaction force between the object and the surface is null, then both objects
detach one of each other. The angle in which that occurs is the last contact
angle, and it was calculated analytically.

We found that the case of a fixed semi-sphere (equivalent to $m_{B}\gg m_{A}$)
gives the maximum possible last contact angle among any physical configuration
of masses. In addition, we found that this angle depends only on the ratio of
the masses $m_{A}$ and $m_{B}$ and it is independent of the value of the
acceleration of gravity or the radius of the semi-sphere.

In addition, we include the effect produced by the initial velocity that the
sliding object might have. This effect is parameterizes in terms of the
quantity, $\epsilon$, which is the ratio between the initial kinetic energy
and the potential energy. We found also that if the kinetic energy is larger
than 1/2 of the initial potential energy, then the sliding object detaches at
the very beginning of the movement.

The problem discussed in this manuscript present a general physical scenario
that might be worth to be used as an example for enthusiast students in
courses of Newtonian mechanics at the undergraduate or graduate level. In
addition, to verify experimentally might also be interesting.

We thank Nicolas Rojas, Fernando Guzmán, Julio Yañez, and Eduardo Peinado for
useful discussions and comments. We also thank the comments and suggestions
from the anonymous referees. RL is supported by Universidad Católica del Norte
through the Publication Incentive program No. CPIP20180343 and CPIP20200063.

##  Appendix A Python code

Here, we present an example of a code in python for solving equation 48 and
return the value of the last contact angle in degrees.

⬇

import numpy as np

import scipy as sp

import matplotlib.pyplot as plt

import scipy.optimize as spop

def H(xi,tau,eps):

PI = np.pi #value of pi

S2 = (np.sin(PI*tau/2))**2 #sin^2

C2 = (np.cos(PI*tau/2))**2 #cos^2

out = S2*xi**3 - 3*xi + 2 + 2*eps*C2 # function H(xi)

return out

# initial values for tau and epsilon

tau = 0.0 #tau

eps = 0.5 #epsilon

root = spop.fsolve(H, 0.5, args=(tau,eps))

print("The last contact angle for eps =", eps, "and tau =",tau, "is:",
np.arccos(root[0])*180/np.pi, "degrees")

## References

## References

  * [1] Beer F P, Johnston E R, Mazurek D F, Cornwell P J and Self B P 2021 Vector Mechanics for Engineers ISBN 9781259638091
  * [2] Hibbeler R C 2017 Engineering mechanics (Harlow: Pearson Education) ISBN 978-1-292-08872-3
  * [3] Riley W 1996 Engineering mechanics : statics (New York: John Wiley) ISBN 0-471-05333-3
  * [4] Balart L and Belmar-Herrera S 2019 American Journal of Physics 87 982–985 URL https://doi.org/10.1119/10.0000037
  * [5] Prior T and Mele E J 2007 American Journal of Physics 75 423–426 URL https://doi.org/10.1119/1.2410018
  * [6] González-Cataldo F, Gutiérrez G and Yáñez J M 2017 American Journal of Physics 85 108–114 URL https://doi.org/10.1119/1.4966628
  * [7] del Pino L A and Curilef S 2018 American Journal of Physics 86 470–471 URL https://doi.org/10.1119/1.5029480
  * [8] Vermillion R E and Cook G O 1988 American Journal of Physics 56 438–439 URL https://doi.org/10.1119/1.15573
  * [9] Waerden B L and van der Waerden (Bartel L) 1985 A History of Algebra ISBN 9780387136103


2022

[1]Jakub Mojsiejuk

2]Institute of Electronics, AGH University of Science and Technology, Al.
Mickiewicza 30, Kraków, 30-059, Poland

2]Faculty of Physics and Applied Computer Science, AGH University of Science
and Technology, Al. Mickiewicza 30, Kraków, 30-059, Poland

# A Comprehensive Simulation Package for Analysis of Multilayer Spintronic
Devices

<EMAIL_ADDRESS>Sławomir Ziętek  Krzysztof Grochot  Witold Skowroński
Tomasz Stobiecki  [  [

###### Abstract

We present cmtj – a comprehensive simulation package that allows large-scale
macrospin simulations for a variety of multilayer spintronics devices. Apart
from conventional static simulations, such as magnetoresistance and
magnetisation hysteresis loops, cmtj implements a mathematical model of
dynamic experimental techniques commonly used for spintronics devices
characterisation, for instance: spin diode ferromagnetic resonance, pulse-
induced microwave magnetometry, or harmonic Hall voltage measurements. We
demonstrate the accuracy of the macrospin simulations on a variety of
examples, accompanied by some experimental results.

###### keywords:

macrospin, Landau–Lifshitz-Gilbert-Slonczewski, spintronics, numerical

##  1 Introduction

Modern development of spintronic devices dieny_opportunities_2020  requires
careful design and a time-consuming fabrication process, for which preparation
is often carried out with the help of simulations. The use of magnetic
materials and multilayer structures with certain material parameters enables
the design of a spintronic device with optimal functionality, depending on the
application: memory (high endurance, low energy writing)
ikegawa_magnetoresistive_2020 , sensor (high sensitivity, low noise) xu_macro-
spin_2017  or high frequency component (high quality factor, low energy
consumption, and low noise) hirohata_review_2020 .

As we climb from the macromagnetic models, steadily increasing the resolution
of the phenomena, and finally reaching the ab initio simulations, we do so
with increasing computational cost. This computational cost is closely
correlated with the number and complexity of the simulation parameters, which
hampers the speed of prototyping and the design of the experiment. As a
result, we are choosing between a slower, but accurate approach and a faster,
but not as precise method. In this paper, we do not undertake to solve this
dilemma, but introduce tools that bring flexibility into the process of making
that decision.

Micromagnetic packages, such as MuMax3 vansteenkiste_design_2014  and OOMMF
m._j._donahue_oommf_1999 , offer a lot of plasticity in modelling magnetic
interactions of complex structures while maintaining an acceptable
computational cost. On the other side of the spectrum, VASP kresse_ab_1993
and Quantum Espresso giannozzi_quantum_2009 ; giannozzi_advanced_2017  lead
the way as acknowledged self-consistent solvers. In the macromagentic regime,
a notable research direction was devoted to creating magnetic tunnel junction
(MTJ) behavioural models in Verilog, for example, in the work of R. Garg et
al. garg_behavioural_2012  or T. Chen et al. chen_comprehensive_2015 ;
tingsu_chen_comprehensive_2015 . Such models attempt to capture the electronic
nature of spintronics devices, describing them as discrete elements or
electronic elements in a circuit, making it easy to prototype devices composed
of discrete components. However, to our best knowledge, so far there have not
been any wide-spread, easy-to-use, computationally efficient, open-sourced
library that allows for quick testing and verification of ideas using
macromagnetic simulations for magnetic multilayer systems including, but not
limited to, static and dynamic characteristics of spin valves and magnetic
tunnel junctions (MTJs). Thus, our main objective is to fill a gap in the
magnetic simulation hierarchy by providing a standardised Python package,
cmtj, for rapid prototyping and simulation. cmtj enables simulation of
current-induced magnetisation dynamics calculations originating from spin
transfer torque (STT) slonczewski_current-driven_1996 ; ogrodnik_study_2021
such as spin torque oscillator (STO) kiselev_microwave_2003 , STT-induced
magnetisation switching myers_current-induced_1999  or voltage-controlled
magnetic anisotropy (VCMA) phenomena nozaki_electric-field-induced_2012 . Some
of those motivating examples are shown in Fig.1. In this work, we demonstrate
a couple of standard use cases and workflows where cmtj provides valuable
information, showing its ability to verify experimental setups. As an
additional asset, we discuss PyMag, a graphical user interface (GUI) that is
beneficial also for educational purposes, profiting from a user-friendly
navigable front-end, but using a fast C++ back-end provided by cmtj.

First, we discuss the main theoretical components of the macromagnetic
simulation. Then, we describe the structure of the simulation in the cmtj and
critical software functionalities, such as software and driver systems.
Finally, we discuss a couple of more advanced simulations that may be
performed with cmtj. In the Appendix, we develop a short walk-through example
for PyMag.

Figure 1: Example trajectories extracted by simulating a single ferromagnetic
layer with different parameters and excitations. (a) depicts a stable STO
(STT-excited system) trajectory under a constant current density. (b)
demonstrates a trajectory obtained with exciting the magnetisation with VCMA.
In (c) we see a trajectory under a pulse excitation of the Oersted field.
Finally, (d) is a thin, bilayer ferromagnetic system with a low energy barrier
that changes states from parallel to anti-parallel under thermal noise. All
trajectories may be reproduced from the example codes available in
documentation to cmtj.

##  2 cmtj backend

The cmtj core is written in C++ and provides a simple header-only library
interface. This approach weighs on accessibility and compatibility with most
platforms, as the user does not have to go through a complicated installation
process. If the user wishes to benefit from the provided python interface, the
setup involves a standard pip installation process. The scope of the python
interface covers all basic functionalities of the C++ core library, exposing
the functions using the open source solution Pybind11
noauthor_pybindpybind11_2022 . In addition, the python package offers utility
functions that complement frequently used operations such as unit conversions,
parallelism, parameter sweep, filtering, energy and resistance calculations,
or template procedures for spin diode ferromagnetic resonance (SD-FMR) or
pulse-induced microwave magnetometry (PIMM). The core library is composed of
couple classes – mainly the Junction and Layer classes that define a basic
magnetic component in the MTJ simulation and the Driver class that contains
definitions of various excitations that influence the magnetisation dynamics
of the system. More detailed descriptions of the cmtj functionality, as well
as implementation details, can be found in the library documentation. We
describe the Layer, Junction, and Driver interfaces in the following sections;
details of additional functionality can be found in the documentation.

###  2.1 Magnetic contributions

The pivotal equation of the magnetic macrospin simulations is called Landau-
Lifshitz-Gilbert-Slonczewski (LLGS) and takes the following form
gilbert_classics_2004 ; ralph_spin_2008 ; slonczewski_current-driven_1996 ;
nguyen_spinorbit_2021 ; ament_solving_2017 :

$\frac{\textrm{d}\textbf{m}}{\textrm{dt}}=-\gamma_{0}\textbf{m}\times\textbf{H}_{\mathrm{eff}}+\alpha_{\textrm{G}}\textbf{m}\times\frac{\textrm{d}\textbf{m}}{\textrm{dt}}\\\
-\gamma_{0}H_{\textrm{FL}}(\textbf{m}\times\textbf{p})-\gamma_{0}H_{\textrm{DL}}(\textbf{m}\times\textbf{m}\times\textbf{p})$
(1)

where  $\ignorespaces\mathbf{m}=\frac{\mathbf{M}}{M_{s}}$ is a normalised
magnetisation vector, with  $\ignorespaces M_{\textrm{s}}$ as the
magnetisation saturation, $\alpha_{\textrm{G}}$ as the dimensionless Gilbert
damping parameter, $\mathbf{p}$ is the polarisation vector, and
$\ignorespaces\gamma_{0}$ is the gyromagnetic factor. The factors
$H_{\textrm{DL}}$ and $H_{\textrm{FL}}$ are so-called damping- and field-like
torque amplitudes, respectively, and generally take the form of STT or SOT,
depending on the experiment. Each of those forms has a separate set of
parameters and differs fundamentally in interpretation; for details, consult
the documentation.

The effective field vector $\mathbf{H}_{\mathrm{eff}}$ is usually composed of
various field contributions that, depending on the context of the simulation,
may be added or disabled. The simulation package in question already provides
a range of such contributions:

$\mathbf{H}_{\textrm{eff}}=\mathbf{H}_{\textrm{ext}}+\mathbf{H}_{\textrm{IEC}}+\mathbf{H}_{\textrm{Oe}}+\mathbf{H}_{\textrm{K}}+\mathbf{H}_{\textrm{demag}}+\mathbf{H}_{\textrm{dipole}}+\mathbf{H}_{\textrm{th}}^{*}+\mathbf{H}_{\textrm{1/f}}^{*}$
(2)

where each component corresponds, respectively, to the applied external field,
interlayer exchange coupling (IEC), Oersted field, anisotropy field,
demagnetising, dipole, thermal, and $1/f$ noise field perez_extending_nodate ;
voss_1f_1978 . Contributions marked with ∗, if enabled in the simulation,
require a stochastic solver, which is described in Sect.2.3. Each of the
contributions that constitute $\mathbf{H}_{\textrm{eff}}$ may be varied over
time using a driver system, as laid out in Sect.2.5.

###  2.2 Simulation design

The simulations in cmtj are conceptually divided into three levels of
abstraction: single-layer simulations, multilayer simulations, and stacked
device simulations. The first of those levels is designed to provide a basic
interface of the ferromagnetic (FM) layer which, at low level, is described by
a time-dependent magnetisation vector.

Multilayer simulation applications revolve most commonly around simple bilayer
structures consisting of heavy metal (HM) and FM layers, or MTJs composed of
FM layers, each separated by a thin tunnel barrier (TB). The non-ferromagnetic
layers such as HM are naturally not simulated in the package, but their
inclusion in the experiment has direct consequences in the simulation by its
SOT property. For example, between two FM layers separated with TB there is
IEC that varies with thickness katayama_interlayer_2006 ; koziol-
rachwal_interlayer_2017 . Similarly, layers can be coupled in the multilayer
system via a longer-range dipole interaction. HMs are typically simulated
indirectly with the addition of STT or SOT terms liu_spin-torque_2012 ;
miron_perpendicular_2011 . It is also possible to save computing time with the
inclusion of pinned layers, experimentally obtained by the exchange-bias
structure enforced by the presence of an antiferromagnet. In such a case, the
solver is not run for a pinned/reference layer, but the effects of SOT or STT
are still modelled via a constant reference vector. Finally, stack
simulations, described in more detail in Sect.2.6, refer to the composites of
multilayer devices, which exhibit various electric coupling.

###  2.3 Solver methods

The core solver for most of the systems in cmtj is the Runge-Kutta 4
algorithm, as it balances decent convergence with speed. However, if the user
specifies a temperature component for the system, cmtj switches to the
stochastic solver, using the Euler-Heun method, and solves the Stratonovich
form of the LLG equation bertotti_chapter_2009 :

$\mathrm{d}\mathbf{m}(t)=\mathbf{f}(\mathbf{m}(t),t)dt+\mathbf{g}(\mathbf{m}(t),t)\circ\mathrm{d}\mathbf{W_{t}}$
(3)

where the non-stochastic part $\mathbf{f}(\mathbf{m}(t),t)$ is the LL form of
the LLG equation, while the stochastic part, $\mathbf{g}(\mathbf{m}(t),t)$
contains thermal and other stochastic contributions:

$\mathbf{g}(\mathbf{m}_{t},t)\circ\mathrm{d}\mathbf{W_{t}}=-\frac{\sigma\gamma}{1+\alpha^{2}}[\mathbf{m}\times\mathrm{d}\mathbf{W_{t}}+\alpha\mathbf{m}\times(\mathbf{m}\times\mathrm{d}\mathbf{W_{t}})]$
(4)

The $\mathrm{d}\mathbf{W_{t}}$ is the random unit vector whose components are
sampled from the normal distribution with zero mean, $\mathcal{N}(0,1)$.

###  2.4 Junction class system

In cmtj each FM layer is by default independently solved, unless there is some
coupling involved. To facilitate multilayer experiments, multiple FM layers
may be combined into a Junction with global control of certain magnetic and
dynamic parameters. For instance, one can apply the same magnetic field
drivers to all or selected layers of the multilayer system using the Junction
methods. This approach simplifies the code and provides a consistent API for
both individual Layer objects and whole Junction.

###  2.5 Driver class system

The dynamic pathway to excite any system in cmtj takes place through the
driver system. The user can define any excitation in the form of a driver, an
adequate Scalar driver, or a vector (Axial) driver. The latter are just
compositions of the former along each dimension (the x, y, and z axes
independently). The Scalar driver is calculated at each time step, and
influences the selected field contributions. For instance, one may define a
sinusoidal driver and use it as a driving excitation of the anisotropy,
leading to simulation of the VCMA effect. Likewise, the externally applied
magnetic field can be added through the AxialDriver class, by specifying field
contributions along each axis independently.

###  2.6 Stack class system

As mentioned above, devices composed of layers can be further combined into
Stack objects. Each Stack may be of type Parallel or Series, which determines
the resistance calculation. The main consequence of composing a stack device
type is the electrical coupling by a current passing through. The model for
this kind of coupling follows the work of Taniguchi et al.
taniguchi_mutual_2018 . For a parallel or series stack of two junctions, the
current depends on the free magnetizations $\mathbf{m}$ and the pinning layers
$\mathbf{p}$ of the junctions $i_{th}$ and $i+1_{th}$ as follows:

$\displaystyle I(t)=I_{0}(t)+\chi
I_{0}(t)(\mathbf{m_{i}}\cdot\mathbf{p}_{i}+\mathbf{m}_{i+1}\cdot\mathbf{p}_{i+1})\quad\textrm{(series)}$
(5) $\displaystyle I(t)=I_{0}(t)+\chi
I_{0}(t)(\mathbf{m_{i}}\cdot\mathbf{p}_{i}-\mathbf{m}_{i+1}\cdot\mathbf{p}_{i+1})\quad\textrm{(parallel)}$
(6)

where $I_{0}(t)$ represents the value of the uncoupled current. The coupling
strength $\chi$ may be positive or negative, most of the time strictly much
less than 1 in absolute value.

##  3 Examples

In this section, we focus on reproducing selected interesting experimental
data and, where possible, comparing the simulation result with the
experimental data. The examples are arranged by the level of complexity,
starting from the simplest towards more advanced. We have chosen the following
techniques to showcase in this section: standard M(H) and R(H) loops in
Sect.3.1, magnetoresistance (MR) based SD-FMR and PIMM in Sect.3.2, harmonic
Hall voltage detection, angular variant in Sect.3.3, current-induced
magnetisation switching (CIMS) in Sect.3.4, electrical coupling of two MTJs in
Sect.3.5.

All the examples presented in this section can be reproduced by running the
Jupyter notebooks from the examples section on the GitHub repository of the
cmtj package.

###  3.1 M(H) and R(H) loops

The M(H) and R(H) are the basic magnetic characterisation methods that enable
the determination of various material parameters of the multilayer system,
such as magnetisation saturation, magnetic anisotropy, or magnetoresistance
ratio. The M(H) loop simulates magnetometry measurements (such as a vibrating
sample magnetometer (VSM) or magnetooptical Kerr effect (MOKE)), whereas the
R(H) loops are indispensable in adjusting the magnetoresistance magnitude.
Furthermore, by performing angular sweep simulations in different planes, one
can also determine a type of magnetoresistance, because anisotropic MR (AMR),
giant MR (GMR)/tunnelling MR (TMR) and spin-Hall MR (SMR) are characterised by
different angular dependencies kim_spin_2016 ; cho_large_2015 ;
rzeszut_biaxial_2018 . Other magnetoresistance flavours, such as those
presented in, for example, the works by avci_unidirectional_2015 ;
velez_hanle_2016 , may be easily added after the magnetisation dynamics has
been computed by cmtj.

Figure 2: M(H) (a-c) and R(H) (d-f) loop of a Co/Ru/Co trilayer system
described in the example. Parameters taken from Tab.3. The magnetic moment is
normalised to unity, whereas the SMR, AMR and GMR magnitudes are set to -0.24
$\Omega$, -0.045 $\Omega$ and 2 $\Omega$, respectively (please note that the
$m_{y}$ and $m_{z}$ components show very little variation with field H applied
along y-axis).

Magnetic field-dependent simulations of the Co(4 nm)/Ru(0.53 nm)/Co(4 nm)
system are presented in Fig.2(a-f). In this example, two FMs are strongly
antiferromagnetically coupled through a Ru spacer mckinnon_spacer_2021 ,
resulting in a zero-net magnetic moment without an external magnetic field.
Increasing H along the y direction leads to a scissor-shaped magnetisation
vector alignment, which is saturated at around 600 kA/m. This is made evident
by the loop M(H) for the $m_{y}$ component illustrated in Fig.2b), with the
remaining, $m_{x}$ and $m_{z}$ having very small amplitudes. Furthermore, we
show that, due to antiferromagnetic coupling (AFM), the two layers oscillate
in antiphase, as shown in Fig.3.

Figure 3: Selected magnetisation trajectory components (a-c) for a
$H_{\textrm{ext}}=-537$\mathrm{k}\mathrm{A}\mathrm{/}\mathrm{m}$$ for the
Co(4nm)/Ru(0.53nm)/Co(4nm) system, marked with the dashed white line in (d).
Colour denotes a given layer, blue is top and red is the bottom layer. (c)
shows that the $m_{z}$ components of the two layers oscillate in antiphase. In
Fig.3(d) we display Kittel dispersion relation obtained from SD-FMR
measurements (described in Sect.4), with experimental points given by blue
dots, and phase difference $\Delta\phi_{z}$ between $m_{z}$ components of two
layers (solid white line) across the external magnetic field range. The region
between approximately -600 $\mathrm{k}\mathrm{A}\mathrm{/}\mathrm{m}$ and 600
$\mathrm{k}\mathrm{A}\mathrm{/}\mathrm{m}$, where the optical branch is
visible, exhibits antiphase oscillation of $m_{z}$ components. Parameters
taken from Tab.3.

###  3.2 GMR system

A more complex example section begins with a study of a
CoFe(2.1nm)/Cu(1.9–2.37nm)/CoFe(1nm)/NiFe(5nm) GMR sample with variable
thickness of the Cu spacer layer. This system is characterised by oscillatory
coupling with the spacer thickness well described in terms of the Rudermann-
Kittel-Kasyua-Yosida (RKKY) interaction between the magnetisations of the
reference and the free layers zietek_influence_2015 . In the case of GMR/TMR,
the magnetoresistance $R(\theta)$ is calculated according to the resistances
in parallel ($R_{P}$) and antiparallel states ($R_{AP}$):

$R(\theta)=R_{P}+\frac{1}{2}(R_{AP}-R_{P})(1-\mathbf{m}_{free}\cdot\mathbf{m}_{reference})$
(7)

where $\theta$ is the angle between the magnetisation vectors of the free and
reference layer, $m_{free}$ is the magnetisation of the free layer and
$m_{reference}$ is the magnetisation of the reference layer. Using the
parameters for the system from the papers zietek_rectification_2015 ;
zietek_influence_2015 , also summarised in Tab.4, we perform spin SD-FMR and
PIMM simulations.

####  3.2.1 SD-FMR

Figure 4:  Kittel dispersion relations (a-c) for 3 different values of Cu
thickness in CoFe(2.1nm)/Cu(1.9–2.37nm)/CoFe(1nm)/NiFe(5nm) sample, with the
blue dots marking the experimental data. Resistance loops at
$\phi_{H}=90^{\circ}$ (d) and at $\phi_{H}=0^{\circ}$ (e) applied external
magnetic field with $\phi$ denoting the polar angle. Parameters of the system
may be found in Tab.4.

In the simulation setup discussed, a representative of SD-FMR is modelled
using the 2D map, $V_{\textrm{mix}}$ vs. $H$ and $f$, $V_{mix}(f,H)$, where
$V_{mix}$ is the mixing voltage of $R(\theta)$ and $I_{\textrm{AC}}$ of a
frequency $f$, under an applied external field $H$. The results are presented
in Fig.4 with open blue circles as the experimental points for comparison.
Afterwards, we turn to the analysis of magnetisation dynamics in the GMR
system by investigating its response using the SD-FMR method. In this
experiment an AC current, $I_{\textrm{AC}}$, is supplied in a given frequency
range, typically between 2 and 40 GHz, while performing a sweep of an external
magnetic field. As a result, the DC mixing voltage $V_{\textrm{mix}}$ is
obtained as a function of both the frequency and the magnitude of the external
magnetic field tulapurkar_spin-torque_2005 . This method allows for the
calculation of the ferromagnetic resonance frequency for a given set of system
parameters in the form of a Kittel dispersion relation. Furthermore, analysis
of the shape of the simulated signal can be used to determine the spin torque
values, in the current perpendicular to the plane STT sankey_measurement_2008
and current in the plane SOT geometries liu_spin-torque_2012 .

Figure 5: PIMM of the CoFe(2.1nm)/Cu(1.9–2.37nm)/CoFe(1nm)/NiFe(5nm) system
described in Sec.3.2. (a) PIMM simulations for a range of IEC values in
$(-1\,$\mathrm{m}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$,1\,$\mathrm{m}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$)$
combined into a single map. Here, the colour indicates the total frequency
amplitude summed over all PIMM simulations. (b) illustrates the same PIMM, but
with marked maximum amplitude resonance line for each $J$ value mapped
accordingly to the colour scale. The colour scale attached represents
different IEC values, from
$-1\,$\mathrm{m}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$$ to
$1\,$\mathrm{m}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$$. We observe the shift of
the resonance curve towards the centre for smaller absolute IEC values. Note
that at $H=0$ the system still exhibits non zero oscillation. Parameters used
for simulation can be found in Tab. 4. Field applied along at $45^{\circ}$
angle between x and y axis.

####  3.2.2 PIMM

Apart from the SD-FMR method described above, magnetisation dynamics can also
be investigated by analysing the magnetisation response to picosecond magnetic
field pulses. In the experiment, this time domain method is performed
indirectly using PIMM silva_inductive_1999 ; serrano-guisan_inductive_2011 ;
banasik_magnetic_2015 . The damped oscillatory response of the multilayer
system to the short magnetic field pulse contains information on the
characteristic resonance. From the simulation perspective, PIMM is simulated
as an analysis of the free oscillations contrary to the dynamics induced by
the external alternating signal, as is the case for the FMR. If the FFT
computation of the pulse response for each value of the applied magnetic field
from the sweep is repeated, one can obtain a spectrogram representing the
dispersion relation where each pixel represents the FFT magnitude for a given
frequency and magnetic field. For the GMR system, multiple PIMM simulations
are conducted for a range of different IEC values, resulting in Fig.5(c).
Clearly, as the magnitude of the IEC coupling increases, the resonance curves
shift away from the zero field. This example shows the ease of predicting
characteristic resonance frequencies in multilayer systems.

###  3.3 Angular harmonic Hall voltage detection

Figure 6: Second Harmonic Hall voltage detection in angular variation, dots
represent experiment, lines demark the simulation data. The y-axis is
normalised. Measurement data collected from Pt(4.3nm)/FeCoB(2nm) device.

There are several experimental methods that allow the determination of SOT
components, such as the SOT-FMR line shape liu_spin-torque_2011  and the line
width analysis pai_spin_2012 , SOT switching current analysis hao_giant_2015 ,
or field-dependent harmonic Hall voltage techniques kim_layer_2013 . Another
method, less susceptible to various artefacts, such as the anomalous Nernst
effect, is called the angular-dependent harmonic Hall voltage
nguyen_spinorbit_2021 . Our model has been thoroughly verified using the
field-dependent method zietek_numerical_2022 , here we present a simplified
angular variation of the harmonic Hall voltage detection. We show the
simulation result, along with the experimental results for the corresponding
Pt(4.3nm)/FeCoB(2nm) sample skowronski_angular_2021 , in Fig.6. Each line is
produced by sweeping with the azimuth angle
$\phi\in[-180^{\circ},180^{\circ}]$ at the frequency excitation of the AC
current $f$ and measuring the second harmonic output in the mixing voltage of
the input current (for details, see skowronski_angular_2021 ). Very good
agreement is found between the experimental data and the simulation results.
The parameters for this system can be found in Tab.5.

###  3.4 CIMS

Next, an example of the SOT-induced magnetisation switching in the HM/FM
bilayer is discussed. The experimental data were obtained in the multilayer
system: Pt(4nm)/Co(1nm)/MgO (A1) from the work by Grochot et al.
grochot_current-induced_2021 . In this case, we reproduce the theoretical
switching behaviour of the system using the field-like and damping-like SOT,
and magnetic parameters obtained from the experiment. The result is shown in
Fig.7. We approximated the critical current analytically using a formula from
Lee et al. lee_threshold_2013 :

$j_{\textrm{sw}}\approx\frac{2e\mu_{0}M_{\textrm{s}}t_{\textrm{FM}}}{\hbar\theta_{SH}}(\frac{H_{\textrm{K,eff}}}{2}-\frac{H_{x}}{\sqrt{2}})$
(8)

where $e$ is the electron charge, $\mu_{0}$ is the magnetic permeability in
vacuum, $M_{\textrm{s}}$ is the magnetisation saturation, $\hbar$ is the
reduced Planck constant, $\theta_{SH}$ is the effective spin Hall angle,
$H_{\textrm{K,eff}}$ is the effective perpendicular anisotropy value, and
$H_{x}$ is the applied field along the $x$ direction.

Figure 7: Analytical and simulated critical current density vs. external field
$H_{\textrm{x}}$ for Pt(4nm)/Co(1nm)/MgO SOT device.

For simulations, we used the trapezoidal impulse shape, with a rising and
falling edge of 1 $\mathrm{n}\mathrm{s}$ and a flat edge of 3
$\mathrm{n}\mathrm{s}$. We also normalise the damping- and field-like torque
magnitudes obtained from the experiment by the current density with which they
have been measured. During the current density sweep, they will scale
proportionally, giving the correct values of $H_{\textrm{DL}}$ and
$H_{\textrm{FL}}$ at each step.

###  3.5 Electrical coupling

Figure 8: The electric synchronisation of two example MTJs, with electric
coupling constant set to $\chi=-0.12$ (a) and $\chi=0.1$ (b). Solid blue and
red lines in (a) indicate the primary resonance modes of individual MTJs from
the stack. The Fig.8(c) shows the same device with no coupling is present.

Finally, this example illustrates how electric coupling can be simulated using
cmtj. First, we created two MTJs with slightly different magnetic and
resistance parameters, emulating experimental dispersion. Then, using the
Stack class interface (see Sect.2.6), we couple them in a parallel connection,
setting the coupling value $\chi$. Finally, we sweep the external magnetic
field and measure the frequency response of the stack to constant current. The
parameters of two systems have been collected in Tab.6. In our case, we
perform a sweep across a range of $H_{\textrm{ext}}$ values, applied at the
azimuthal angle of $5^{\circ}$. Then, we measure the frequency response of the
coupled system to the constant current density excitation. For larger values
of the applied external magnetic field, and negative coupling constant $\chi$,
we observe how two main resonance lines, each corresponding to a separate MTJ,
converge towards a common resonance mode. Ultimately, around 250
$\mathrm{k}\mathrm{A}\mathrm{/}\mathrm{m}$, the two MTJs desynchronise and
their main resonance lines separate again. The result of that electric
synchronisation of two MTJs is depicted in Fig.8(a), while Fig.8(b-c)
illustrate a situation with positive coupling coefficient and coupling
disabled, respectively. This example shows that the software presented can
also be used in more complex systems, for example, for the analysis of neural
computing platforms romera_binding_2022 .

###  3.6 Benchmarks

The cmtj is meant to run also on personal computers, thus in Tab.1 we present
example execution times that were measured on 2020 MacBook Pro with Apple
Silicon M1 and 16 GB RAM. The cmtj library includes additional helper
functions that allow for easy parallelism if the simulation setup allows it.
Running certain simulations in parallel without considering initial position
effects may produce drastically different results. For instance, while
sweeping with an external magnetic field, we can treat each next step as a
slight perturbation relative to the previous one, and thus we can use the
final step magnetisation from $H_{i}$ as the starting position for $H_{i+1}$
simulation. To compensate for this, for example in the parallel setup, one can
extend the simulation with an initial relaxation period to stabilise the
magnetisation first.

Benchmark | Steps [$10^{6}$] | Distributed | Time [s] | Steps/s [$10^{3}$/s]
---|---|---|---|---
M(H), R(H), PIMM | 50 | NO | 42 | 1190
CIMS | 187.5 | YES | 355 | 528
GMR VSD | 1872 | YES | 367 | 5101
GMR PIMM | 360 | YES | 315 | 1143
Stack synchronisation | 120 | NO | 62 | 1935
Table 1: The execution times and the number of steps per experiment are given
in the examples section (Sect.3). The numbers summarise all sweeps over
parameters, fields, and frequencies, and thus represent the total measured
time required to reproduce the experiment. The timing was done on a MacBook
Pro 2020 with an Apple M1 processor and 16 GB of RAM. Selected examples
benefitted from parallel computation to show the ease of use of the included
package parallelism method based on the Python multiprocessing library. Some
experiments, such as, for instance, CIMS (15 ns), involve more parameter
scanning than individual simulation steps, so they are less efficient than
simulations with longer in-simulation time, such as GMR VSD (60 ns).

##  4 Conclusions

In this article, we presented cmtj, fast modelling software for macrospin
simulation of magnetic heterostructures. Its core is grounded in the LLGS
equation, and the package is capable of calculating both the static and
dynamic characteristics of spintronic devices. In the spirit of the modern
modular development approach, suggested, for example, in locatelli_spin-
torque_2014 ; camsari_modular_2015 , cmtj expands its usability to connect
multiple spintronic structures using the Junction or Stack system.
Furthermore, we show that the results obtained from cmtj agree well with the
experiment and that the calculated simulation derivatives are shown to be
consistent in both the static and dynamic regimes.

##  5 Acknowledgements

We thank P. Ogrodnik for a fruitful and insightful discussion. The research
project is partly supported by the Excellence initiative – research university
programme for the AGH University of Science and Technology. T.S. and K.G. were
supported by the National Science Centre, Poland, Grant No.
UMO-2016/23/B/ST3/01430 (SPINORBITRONICS).

##  Appendix A PyMag example

Figure 9: View of PyMag UI (a) Control panel, (b) sample layer and stimulus
settings panel, (c) R-H plots, (d) M-H plots, (e) map of PIMM results, (f) map
of SD-FMR results, (g) simulations management panel, (h) experimental or
external simulation management panel.

PyMag is a supplementary GUI programme that incorporates certain experimental
procedures using the cmtj interface, such as M(H), R(H), PIMM or SD-FMR in a
user-friendly interface. In this section, we discuss step by step how to
reproduce a simple simulation with PyMag software. PyMag allows for the
definition of the material parameters and stimuli of the spintronics device,
the simulation process, the management of the simulation list, previewing and
analysis of the simulation results, the management of loaded experimental
data, the export of selected results, and the saving and loading of the entire
simulation state. PyMag interface was based on open source libraries: PyQt5
and pyqtgraph Python. In addition, PyMag allows the user to compare their
experimental data with the simulated results using the Experiment management
tab.

The structure parameters of each layer of the device must be defined in the
table; see Fig.9(f). Each row corresponds to a single ferromagnetic layer and
contains a variety of parameters. For convenience we put them here, placing
the units expected by PyMag in brackets: magnetisation saturation
$\mu_{0}M_{\mathrm{s}}($\mathrm{T}$)$, uniaxial anisotropy
$K_{\mathrm{u}}($\mathrm{J}\mathrm{/}{m^{3}}$$), interlayer exchange coupling,
first order $J($\mathrm{J}\mathrm{/}{m^{2}}$)$ and second order
$J_{2}($\mathrm{J}\mathrm{/}{m^{2}}$)$, geometrical properties (thickness,
width , length $d($\mathrm{m}$)$, $w($\mu\mathrm{m}$)$, $l($\mu\mathrm{m}$)$,
respectively), anisotropic magnetoresistance ($\mathrm{AMR}$), spin Hall
magnetoresistance ($\mathrm{SMR}$), anomalous Hall effect ($\mathrm{AHE}$)
(all in $\Omega$), and the following dynamic parameters: damping factor
$\alpha$, field-like torque $H_{\mathrm{FL}}$
($\mathrm{A}\mathrm{/}\mathrm{m}$), damping-like torque $H_{\mathrm{DL}}$
($\mathrm{A}\mathrm{/}\mathrm{m}$). The Oersted field excitation
$H_{\textrm{Oe}}$ is given in the table in $\mathrm{1}\mathrm{/}\mathrm{m}$,
as it is multiplied by the current value entered in the stimulus section. The
following parameters are entered as vector quantities (each being a vector of
dimension 3): $\mathbf{K}_{\mathrm{udir}}$, the direction of uniaxial
anisotropy, $\mathbf{N}_{\textrm{demag}}$, the diagonal of the demagnetisation
tensor, $\mathbf{p}$, the polarisation vector in the SOT torque, and
$\mathbf{H}_{\textrm{Oedir}}$, the direction of the Oersted field.

For a basic, demonstrative example we have chosen a trilayer structure, with
two ferromagnetic layers separated with a spacer layer, Co/Ru/Co, which
exhibits a quadratic interlayer exchange coupling, producing optical and
acoustic lines in the FMR spectrum. The multilayer parameters used in the
simulations are grouped in Tab.3.

Then we move on to the definition of a stimulus, which can be defined in the
Stimulus tab. In PyMag we follow a physical spherical coordinate system that
designates $\theta$ as the polar angle (measured from the positive z-axis),
and $\phi$ as the azimuth angle. In our software, the external field
excitation may be swept either by magnitude or by any of those two angles.
Furthermore, there is a DC ($I_{\mathrm{DC}}$) and AC ($I_{\mathrm{AC}}$)
current that can be supplied in a selected frequency range described by three
parameters: $f_{\mathrm{min}}$, $f_{\mathrm{max}}$, $f_{\mathrm{steps}}$). For
the current example, we keep $I_{\textrm{DC}}$ fixed and sweep the magnetic
field along the y-axis between -600 and 600 kA/m.

Finally, the precision of the simulation is governed by the total simulation
time in nanoseconds: $\textrm{LLG}_{\mathrm{time}}$ and the number of
integration steps ($\textrm{LLG}_{\mathrm{steps}}$) in the RK4 method. We find
that, for most cases, the integration step should not be greater than
$10^{-12}$ to ensure the correctness and proper convergence of the simulation.
However, this depends on multiple factors. For instance, large IEC coupling,
as in the case of the example, may lead to very complex trajectories that
require much shorter integration times to converge properly. We propose
setting 500000 steps per $5$\mathrm{n}\mathrm{s}$$, which corresponds to the
integration step of $10^{-14}$\mathrm{s}$$. When inspecting the Convergence
tab, the user can conveniently observe the regions of good and poor
convergence, defined as a running difference in magnitude of
$\frac{d\mathbf{m}}{dt}$.

The control panel (Fig.9(g-h)) is used to monitor the simulation status. The
user can load or save selected simulations, displaying the experimental data
directly on the simulated runs. The progress bar informs the user about the
stage of the run, and the stop and start buttons control the simulation flow.
To run a simulation, we click on Run simulation button in the control panel.
The graphs in the visualisation tabs are populated online as the simulation
progresses. In PyMag, the results of this aggregation can be found in the PIMM
tab (Fig.9(e))

##  Appendix B Layer parameters

### Parameters explained

Parameter | Meaning
---|---
$\mu_{0}M_{\textrm{s}}$ | Magnetisation saturation
$J$ | IEC value
$J_{\textrm{quad}}$ | IEC value (quadratic term)
$K_{u}$ | Anisotropy value
$\mathbf{K}_{dir}$ | Anisotropy axis
$\alpha_{\textrm{G}}$ | Gilbert damping parameter
$\mathbf{N}_{\textrm{demag}}$ | diagonal of demagnetisation tensor
$t_{\textrm{FM}}$ | thickness of a FM layer
$w$ | width of a sample
$l$ | length of the sample
$H_{\textrm{DL}}$ | magnitude of the damping-like torque
$H_{\textrm{FL}}$ | magnitude of the field-like torque
$\mathbf{p}$ | polarisation vector
$\chi$ | electric coupling constant
$\lambda$ | Slonczewski spacer layer parameter
$\eta$ | spin polarisation efficiency
Table 2: Explanation of the key parameters used in the simulations.

### Example 1 – Co(4nm)/Ru(0.53nm)/Co(4nm)

Parameter | Value | Unit
---|---|---
Magnetic parameters
$\mu_{0}M_{\textrm{s}}^{1,2}$ | 1.65 | $\mathrm{T}$
$K_{u}^{1,2}$ | 1050 | $\mathrm{J}\mathrm{/}\mathrm{m}^{3}$
$\mathbf{K}_{dir}^{1,2}$ | [1, 0, 0] | -
$J^{1}$ | -1.78 | $\mathrm{m}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$
$J_{quad}^{1}$ | -0.169 | $\mathrm{m}\mathrm{J}\mathrm{/}\mathrm{m}^{2}$
$\alpha_{\textrm{G}}^{1,2}$ | 0.005 | -
$\mathbf{N}_{\textrm{demag}}^{1,2}$ | [0.001, 0.000, 0.998] | -
Dimensions
$t_{\textrm{FM}}^{1}$ | 4 | $\mathrm{n}\mathrm{m}$
$t_{\textrm{FM}}^{2(*)}$ | 3.99 | $\mathrm{n}\mathrm{m}$
$w^{1,2}$ | 20 | $\mathrm{u}\mathrm{m}$
$l^{1,2}$ | 30 | $\mathrm{u}\mathrm{m}$
Resistance parameters
$\textrm{AMR}^{1,2}$ | -0.045 | $\mathrm{O}\mathrm{h}\mathrm{m}$
$\textrm{SMR}^{1,2}$ | -0.24 | $\mathrm{O}\mathrm{h}\mathrm{m}$
$\textrm{AHE}^{1,2}$ | -2.7 | $\mathrm{O}\mathrm{h}\mathrm{m}$
$\textrm{GMR}^{1,2}$ | 2 | $\mathrm{O}\mathrm{h}\mathrm{m}$
$\textrm{R}_{\textrm{xx0}}^{1,2}$ | 100 | $\mathrm{O}\mathrm{h}\mathrm{m}$
$\textrm{R}_{\textrm{xy0}}^{1,2}$ | 1 | $\mathrm{O}\mathrm{h}\mathrm{m}$
Table 3: Parameter table for the example trilayer Co(4nm)/Ru(0.53nm)/Co(4nm)
system. The superscript denotes the layer id. If there is no parameter for the
corresponding layer id, the value is 0 by default. $(^{*})$ We distinguish
between the thicknesses of the upper and lower Co layers due to the effect of
the interfacial magnetic dead layer. The parameter explanation may be found in
the previous section, Tab.2

### Example 2 – CoFe(2.1nm)/Cu(1.9–2.37nm)/CoFe(1nm)/NiFe(5nm)

Parameter | Value | Unit
---|---|---
Magnetic parameters
$\mu_{0}M_{\textrm{s}}^{1}$ | 1.03 | $\mathrm{T}$
$\mu_{0}M_{\textrm{s}}^{2}$ | 1.65 | $\mathrm{T}$
$K_{u}^{1}$ | 0.8 | $\mathrm{k}\mathrm{J}\mathrm{/}\mathrm{m}^{3}$
$K_{u}^{2}$ | $1$ | $\mathrm{T}\mathrm{J}\mathrm{/}\mathrm{m}^{3}$
$\mathbf{K}_{dir}^{1,2}$ | [1, 0, 0] | -
$\alpha_{\textrm{G}}^{1,2}$ | 0.024 | -
$\mathbf{N}_{\textrm{demag}}^{1,2}$ | [0, 0, 1] | -
Dimensions
$t_{\textrm{FM}}^{1}$ | 2.1 | $\mathrm{n}\mathrm{m}$
$t_{\textrm{FM}}^{2}$ | 6 | $\mathrm{n}\mathrm{m}$
Resistance parameters
$\textrm{R}_{\textrm{P}}$ | 163.5 | $\mathrm{O}\mathrm{h}\mathrm{m}$
$\textrm{R}_{\textrm{AP}}$ | 176 | $\mathrm{O}\mathrm{h}\mathrm{m}$
Table 4: Parameter table for the example GMR system. The superscript denotes
the layer id.

### Example 3 – Angular harmonics for Pt(4.3nm)/FeCoB(2nm)

Parameter | Value | Unit
---|---|---
Magnetic parameters
$\mu_{0}M_{\textrm{s}}$ | 1.2 | $\mathrm{T}$
$K_{u}$ | 1.49 | $\mathrm{k}\mathrm{J}\mathrm{/}\mathrm{m}^{3}$
$\mathbf{K}_{dir}$ | [0, 0, 1] | -
$\alpha_{\textrm{G}}$ | 0.003 | -
$\mathbf{N}_{\textrm{demag}}$ | [0, 0, 1] | -
Dimensions
$t_{\textrm{FM}}$ | 2 | $\mathrm{n}\mathrm{m}$
$w$ | 10 | $\mu\mathrm{m}$
$l$ | 30 | $\mu\mathrm{m}$
Resistance parameters
AHE | 1.15 | $\mathrm{O}\mathrm{h}\mathrm{m}$
SMR | -0.125 | $\mathrm{O}\mathrm{h}\mathrm{m}$
AMR | $-10^{-4}$ | $\mathrm{O}\mathrm{h}\mathrm{m}$
Torque parameters
$H_{\textrm{DL}}$ | 1054 | $\mathrm{A}\mathrm{/}\mathrm{m}$
$H_{\textrm{FL}}$ | -164 | $\mathrm{A}\mathrm{/}\mathrm{m}$
$\mathbf{p}$ | [0, 1, 0] | -
Table 5: Parameter table for the angular harmonics example, Pt/CoFeB sample.

### Example 4 – Electric coupling

Parameter | Value | Unit
---|---|---
Magnetic parameters
$\mu_{0}M_{\textrm{s}}^{1}$ | 1.6 | $\mathrm{T}$
$\mu_{0}M_{\textrm{s}}^{2}$ | 1.76 | $\mathrm{T}$
$K_{u}^{1}$ | 70 | $\mathrm{k}\mathrm{J}\mathrm{/}\mathrm{m}^{3}$
$K_{u}^{2}$ | 100 | $\mathrm{k}\mathrm{J}\mathrm{/}\mathrm{m}^{3}$
$\mathbf{K}_{dir}^{1,2}$ | [0, 0, 1] | -
$\alpha_{\textrm{G}}$ | 0.005 | -
$\mathbf{N}_{\textrm{demag}}$ | [0, 0, 1] | -
Resistance parameters
$\textrm{R}_{\textrm{P}}^{1}$ | 100 | $\mathrm{O}\mathrm{h}\mathrm{m}$
$\textrm{R}_{\textrm{P}}^{2}$ | 110 | $\mathrm{O}\mathrm{h}\mathrm{m}$
$\textrm{R}_{\textrm{AP}}^{1}$ | 200 | $\mathrm{O}\mathrm{h}\mathrm{m}$
$\textrm{R}_{\textrm{AP}}^{2}$ | 220 | $\mathrm{O}\mathrm{h}\mathrm{m}$
Torque parameters
$\lambda$ | 0.69 | -
$\eta^{1,2}$ | 1 | -
$\mathbf{p}^{1,2}$ | [1, 0, 0] | -
Stack parameters
$\chi$ | -0.12 | -
$j$ | 60 | $\mathrm{M}\mathrm{A}\mathrm{/}\mathrm{m}^{2}$
Table 6: Parameter table for the electric coupling system of two MTJs
connected in parallel. The current density was set to 60
$\mathrm{M}\mathrm{A}\mathrm{/}\mathrm{m}^{2}$. The parameters partially taken
from skowronski_microwave_2019 , with the MTJ composition of CoFeB(2.5nm)/MgO
(1nm)/CoFeB(1.8nm).

## References

  *   * (1) Dieny, B., Prejbeanu, I.L., Garello, K., Gambardella, P., Freitas, P., Lehndorff, R., Raberg, W., Ebels, U., Demokritov, S.O., Akerman, J., Deac, A., Pirro, P., Adelmann, C., Anane, A., Chumak, A.V., Hirohata, A., Mangin, S., Valenzuela, S.O., Onbaşlı, M.C., d’Aquino, M., Prenat, G., Finocchio, G., Lopez-Diaz, L., Chantrell, R., Chubykalo-Fesenko, O., Bortolotti, P.: Opportunities and challenges for spintronics in the microelectronics industry. Nature Electronics 3(8), 446–459 (2020). https://doi.org/10.1038/s41928-020-0461-5. Accessed 2022-06-07
  * (2) Ikegawa, S., Mancoff, F.B., Janesky, J., Aggarwal, S.: Magnetoresistive Random Access Memory: Present and Future. IEEE Transactions on Electron Devices 67(4), 1407–1419 (2020). https://doi.org/10.1109/TED.2020.2965403. Accessed 2022-06-07
  * (3) Xu, Y., Yang, Y., Luo, Z., Xu, B., Wu, Y.: Macro-spin modeling and experimental study of spin-orbit torque biased magnetic sensors. Journal of Applied Physics 122(19), 193904 (2017). https://doi.org/10.1063/1.4994109. Accessed 2022-06-07
  * (4) Hirohata, A., Yamada, K., Nakatani, Y., Prejbeanu, I.-L., Diény, B., Pirro, P., Hillebrands, B.: Review on spintronics: Principles and device applications. Journal of Magnetism and Magnetic Materials 509, 166711 (2020). https://doi.org/10.1016/j.jmmm.2020.166711. Accessed 2022-06-07
  * (5) Vansteenkiste, A., Leliaert, J., Dvornik, M., Helsen, M., Garcia-Sanchez, F., Van Waeyenberge, B.: The design and verification of MuMax3. AIP Advances 4(10), 107133 (2014). https://doi.org/10.1063/1.4899186. Accessed 2018-05-10
  * (6) Donahue, M.J., Porter, D.G.: OOMMF User’s Guide, Version 1.0. Interagency Report NISTIR 6376 National Institute of Standards and Technology, Gaithersburg (1999). http://math.nist.gov/oommf/
  * (7) Kresse, G., Hafner, J.: Ab initio molecular dynamics for liquid metals. Physical Review B 47(1), 558–561 (1993). https://doi.org/10.1103/PhysRevB.47.558. Accessed 2020-01-09
  * (8) Giannozzi, P., Baroni, S., Bonini, N., Calandra, M., Car, R., Cavazzoni, C., Ceresoli, D., Chiarotti, G.L., Cococcioni, M., Dabo, I., Corso, A.D., Gironcoli, S.d., Fabris, S., Fratesi, G., Gebauer, R., Gerstmann, U., Gougoussis, C., Kokalj, A., Lazzeri, M., Martin-Samos, L., Marzari, N., Mauri, F., Mazzarello, R., Paolini, S., Pasquarello, A., Paulatto, L., Sbraccia, C., Scandolo, S., Sclauzero, G., Seitsonen, A.P., Smogunov, A., Umari, P., Wentzcovitch, R.M.: QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials. Journal of Physics: Condensed Matter 21(39), 395502 (2009). https://doi.org/10.1088/0953-8984/21/39/395502. Publisher: IOP Publishing. Accessed 2020-04-04
  * (9) Giannozzi, P., Andreussi, O., Brumme, T., Bunau, O., Nardelli, M.B., Calandra, M., Car, R., Cavazzoni, C., Ceresoli, D., Cococcioni, M., Colonna, N., Carnimeo, I., Corso, A.D., Gironcoli, S.d., Delugas, P., DiStasio, R.A., Ferretti, A., Floris, A., Fratesi, G., Fugallo, G., Gebauer, R., Gerstmann, U., Giustino, F., Gorni, T., Jia, J., Kawamura, M., Ko, H.-Y., Kokalj, A., Küçükbenli, E., Lazzeri, M., Marsili, M., Marzari, N., Mauri, F., Nguyen, N.L., Nguyen, H.-V., Otero-de-la-Roza, A., Paulatto, L., Poncé, S., Rocca, D., Sabatini, R., Santra, B., Schlipf, M., Seitsonen, A.P., Smogunov, A., Timrov, I., Thonhauser, T., Umari, P., Vast, N., Wu, X., Baroni, S.: Advanced capabilities for materials modelling with Quantum ESPRESSO. Journal of Physics: Condensed Matter 29(46), 465901 (2017). https://doi.org/10.1088/1361-648X/aa8f79. Publisher: IOP Publishing. Accessed 2020-04-04
  * (10) Garg, R., Kumar, D., Jindal, N., Negi, N., Ahuja, C.: Behavioural model of Spin Torque Transfer Magnetic Tunnel Junction, Using Verilog-A 1, 7 (2012)
  * (11) Chen, T., Eklund, A., Iacocca, E., Rodriguez, S., Malm, B.G., Åkerman, J., Rusu, A.: Comprehensive and Macrospin-Based Magnetic Tunnel Junction Spin Torque Oscillator Model- Part II: Verilog-A Model Implementation. IEEE Transactions on Electron Devices 62(3), 1045–1051 (2015). https://doi.org/10.1109/TED.2015.2390676. Conference Name: IEEE Transactions on Electron Devices
  * (12) Tingsu Chen, Eklund, A., Iacocca, E., Rodriguez, S., Malm, B.G., Akerman, J., Rusu, A.: Comprehensive and Macrospin-Based Magnetic Tunnel Junction Spin Torque Oscillator Model-Part I: Analytical Model of the MTJ STO. IEEE Transactions on Electron Devices 62(3), 1037–1044 (2015). https://doi.org/10.1109/TED.2015.2390411. Accessed 2021-01-25
  * (13) Slonczewski, J.C.: Current-driven excitation of magnetic multilayers. Journal of Magnetism and Magnetic Materials 159(1-2), 1–7 (1996). https://doi.org/10.1016/0304-8853(96)00062-5. Accessed 2018-06-05
  * (14) Ogrodnik, P., Grochot, K., Karwacki, Ł., Kanak, J., Prokop, M., Chęciński, J., Skowroński, W., Ziętek, S., Stobiecki, T.: Study of Spin–Orbit Interactions and Interlayer Ferromagnetic Coupling in Co/Pt/Co Trilayers in a Wide Range of Heavy-Metal Thickness. ACS Applied Materials & Interfaces 13(39), 47019–47032 (2021). https://doi.org/10.1021/acsami.1c11675. Accessed 2022-06-09
  * (15) Kiselev, S.I., Sankey, J.C., Krivorotov, I.N., Emley, N.C., Schoelkopf, R.J., Buhrman, R.A., Ralph, D.C.: Microwave oscillations of a nanomagnet driven by a spin-polarized current. Nature 425(6956), 380–383 (2003). https://doi.org/10.1038/nature01967. Accessed 2022-06-01
  * (16) Myers, E.B., Ralph, D.C., Katine, J.A., Louie, R.N., Buhrman, R.A.: Current-Induced Switching of Domains in Magnetic Multilayer Devices. Science 285(5429), 867–870 (1999). https://doi.org/10.1126/science.285.5429.867. Accessed 2020-01-08
  * (17) Nozaki, T., Shiota, Y., Miwa, S., Murakami, S., Bonell, F., Ishibashi, S., Kubota, H., Yakushiji, K., Saruya, T., Fukushima, A., Yuasa, S., Shinjo, T., Suzuki, Y.: Electric-field-induced ferromagnetic resonance excitation in an ultrathin ferromagnetic metal layer. Nature Physics 8, 491 (2012)
  * (18) pybind/pybind11. pybind. original-date: 2015-07-05T19:46:48Z (2022). https://github.com/pybind/pybind11 Accessed 2022-06-01
  * (19) Gilbert, T.L.: Classics in Magnetics A Phenomenological Theory of Damping in Ferromagnetic Materials. IEEE Transactions on Magnetics 40(6), 3443–3449 (2004). https://doi.org/10.1109/TMAG.2004.836740. Accessed 2018-05-10
  * (20) Ralph, D.C., Stiles, M.D.: Spin transfer torques. Journal of Magnetism and Magnetic Materials 320(7), 1190–1216 (2008). https://doi.org/10.1016/j.jmmm.2007.12.019. Accessed 2018-08-27
  * (21) Nguyen, M.-H., Pai, C.-F.: Spin–orbit torque characterization in a nutshell. APL Materials 9(3), 030902 (2021). https://doi.org/10.1063/5.0041123. Accessed 2021-10-03
  * (22) Ament, S., Rangarajan, N., Parthasarathy, A., Rakheja, S.: Solving the stochastic Landau-Lifshitz-Gilbert-Slonczewski equation for monodomain nanomagnets : A survey and analysis of numerical techniques. arXiv:1607.04596 [cs] (2017). arXiv: 1607.04596. Accessed 2021-11-13
  * (23) Perez, G.: Extending the Capacity of 1 / f Noise Generation, 10
  * (24) Voss, R.F.: ”1/f noise” in music: Music from 1/f noise. The Journal of the Acoustical Society of America 63(1), 258 (1978). https://doi.org/10.1121/1.381721. Accessed 2022-07-23
  * (25) Katayama, T., Yuasa, S., Velev, J., Zhuravlev, M.Y., Jaswal, S.S., Tsymbal, E.Y.: Interlayer exchange coupling in Fe/MgO/Fe magnetic tunnel junctions. Applied Physics Letters 89(11), 112503 (2006). https://doi.org/10.1063/1.2349321. Accessed 2020-04-06
  * (26) Kozioł-Rachwał, A., Skowroński, W., Frankowski, M., Chęciński, J., Ziętek, S., Rzeszut, P., Ślęzak, M., Matlak, K., Ślęzak, T., Stobiecki, T., Korecki, J.: Interlayer exchange coupling, dipolar coupling and magnetoresistance in Fe/MgO/Fe trilayers with a subnanometer MgO barrier. Journal of Magnetism and Magnetic Materials 424, 189–193 (2017). https://doi.org/10.1016/j.jmmm.2016.09.131. Accessed 2018-05-10
  * (27) Liu, L., Pai, C.-F., Li, Y., Tseng, H.W., Ralph, D.C., Buhrman, R.A.: Spin-Torque Switching with the Giant Spin Hall Effect of Tantalum. Science 336(6081), 555–558 (2012). https://doi.org/10.1126/science.1218197. Accessed 2022-06-07
  * (28) Miron, I.M., Garello, K., Gaudin, G., Zermatten, P.-J., Costache, M.V., Auffret, S., Bandiera, S., Rodmacq, B., Schuhl, A., Gambardella, P.: Perpendicular switching of a single ferromagnetic layer induced by in-plane current injection. Nature 476(7359), 189–193 (2011). https://doi.org/10.1038/nature10309. Accessed 2021-12-09
  * (29) Bertotti, G., Mayergoyz, I.D., Serpico, C.: Chapter 10 - Stochastic Magnetization Dynamics. In: Bertotti, G., Mayergoyz, I.D., Serpico, C. (eds.) Nonlinear Magnetization Dynamics in Nanosystems. Elsevier Series in Electromagnetism, pp. 271–357. Elsevier, Oxford (2009). https://doi.org/10.1016/B978-0-08-044316-4.00012-8. https://www.sciencedirect.com/science/article/pii/B9780080443164000128 Accessed 2022-06-22
  * (30) Taniguchi, T., Tsunegi, S., Kubota, H.: Mutual synchronization of spin-torque oscillators consisting of perpendicularly magnetized free layers and in-plane magnetized pinned layers. Applied Physics Express 11(1), 013005 (2018). https://doi.org/10.7567/APEX.11.013005. Accessed 2021-09-09
  * (31) Kim, J., Sheng, P., Takahashi, S., Mitani, S., Hayashi, M.: Spin Hall Magnetoresistance in Metallic Bilayers. Physical Review Letters 116(9), 097201 (2016). https://doi.org/10.1103/PhysRevLett.116.097201. Accessed 2021-11-13
  * (32) Cho, S., Baek, S.-h.C., Lee, K.-D., Jo, Y., Park, B.-G.: Large spin Hall magnetoresistance and its correlation to the spin-orbit torque in W/CoFeB/MgO structures. Scientific Reports 5(1), 14668 (2015). https://doi.org/10.1038/srep14668. Accessed 2022-06-09
  * (33) Rzeszut, P., Skowroński, W., Ziętek, S., Ogrodnik, P., Stobiecki, T.: Biaxial Magnetic-Field Setup for Angular-Dependent Measurements of Magnetic Thin Films and Spintronic Nanodevices. IEEE Transactions on Magnetics 54(8), 1–7 (2018). https://doi.org/10.1109/TMAG.2018.2845871. Accessed 2022-07-08
  * (34) Avci, C.O., Garello, K., Ghosh, A., Gabureac, M., Alvarado, S.F., Gambardella, P.: Unidirectional spin Hall magnetoresistance in ferromagnet/normal metal bilayers. Nature Physics 11(7), 570–575 (2015). https://doi.org/10.1038/nphys3356. Accessed 2022-06-13
  * (35) Vélez, S., Golovach, V.N., Bedoya-Pinto, A., Isasa, M., Sagasta, E., Abadia, M., Rogero, C., Hueso, L.E., Bergeret, F.S., Casanova, F.: Hanle Magnetoresistance in Thin Metal Films with Strong Spin-Orbit Coupling. Physical Review Letters 116(1), 016603 (2016). https://doi.org/10.1103/PhysRevLett.116.016603. Accessed 2022-06-13
  * (36) McKinnon, T., Heinrich, B., Girt, E.: Spacer layer thickness and temperature dependence of interlayer exchange coupling in Co/Ru/Co trilayer structures. Physical Review B 104(2), 024422 (2021). https://doi.org/10.1103/PhysRevB.104.024422. Accessed 2022-07-13
  * (37) Ziętek, S., Ogrodnik, P., Skowroński, W., Wiśniowski, P., Czapkiewicz, M., Stobiecki, T., Barnaś, J.: The influence of interlayer exchange coupling in giant-magnetoresistive devices on spin diode effect in wide frequency range. Applied Physics Letters 107(12), 122410 (2015). https://doi.org/10.1063/1.4931771. Accessed 2022-05-11
  * (38) Ziętek, S., Ogrodnik, P., Frankowski, M., Chęciński, J., Wiśniowski, P., Skowroński, W., Wrona, J., Stobiecki, T., Żywczak, A., Barnaś, J.: Rectification of radio frequency current in giant magnetoresistance spin valve. Physical Review B 91(1), 014430 (2015). https://doi.org/10.1103/PhysRevB.91.014430. Accessed 2022-04-16
  * (39) Tulapurkar, A.A., Suzuki, Y., Fukushima, A., Kubota, H., Maehara, H., Tsunekawa, K., Djayaprawira, D.D., Watanabe, N., Yuasa, S.: Spin-torque diode effect in magnetic tunnel junctions. Nature 438(7066), 339–342 (2005). https://doi.org/10.1038/nature04207. Accessed 2022-07-14
  * (40) Sankey, J.C., Cui, Y.-T., Sun, J.Z., Slonczewski, J.C., Buhrman, R.A., Ralph, D.C.: Measurement of the spin-transfer-torque vector in magnetic tunnel junctions. Nature Physics 4(1), 67–71 (2008). https://doi.org/10.1038/nphys783. Accessed 2022-06-07
  * (41) Silva, T.J., Lee, C.S., Crawford, T.M., Rogers, C.T.: Inductive measurement of ultrafast magnetization dynamics in thin-film Permalloy. Journal of Applied Physics 85(11), 7849–7862 (1999). https://doi.org/10.1063/1.370596. Accessed 2021-12-09
  * (42) Serrano-Guisan, S., Skowroński, W., Wrona, J., Liebing, N., Czapkiewicz, M., Stobiecki, T., Reiss, G., Schumacher, H.W.: Inductive determination of the optimum tunnel barrier thickness in magnetic tunneling junction stacks for spin torque memory applications. Journal of Applied Physics 110(2), 023906 (2011). https://doi.org/10.1063/1.3610948. Accessed 2022-07-08
  * (43) Banasik, M., Wrona, J., Kanak, J., Ziętek, S., Skowroński, W., Żywczak, A., Czapkiewicz, M., Stobiecki, T.: Magnetic Properties and Magnetization Dynamics of Magnetic Tunnel Junction Bottom Electrode With Different Buffer Layers. IEEE Transactions on Magnetics 51(11), 1–4 (2015). https://doi.org/10.1109/TMAG.2015.2440561. Accessed 2022-07-08
  * (44) Liu, L., Moriyama, T., Ralph, D.C., Buhrman, R.A.: Spin-Torque Ferromagnetic Resonance Induced by the Spin Hall Effect. Physical Review Letters 106(3), 036601 (2011). https://doi.org/10.1103/PhysRevLett.106.036601. Accessed 2021-12-09
  * (45) Pai, C.-F., Liu, L., Li, Y., Tseng, H.W., Ralph, D.C., Buhrman, R.A.: Spin transfer torque devices utilizing the giant spin Hall effect of tungsten. Applied Physics Letters 101(12), 122404 (2012). https://doi.org/10.1063/1.4753947. Accessed 2022-06-28
  * (46) Hao, Q., Xiao, G.: Giant Spin Hall Effect and Switching Induced by Spin-Transfer Torque in a W / Co 40 Fe 40 B 20 / MgO Structure with Perpendicular Magnetic Anisotropy. Physical Review Applied 3(3), 034009 (2015). https://doi.org/10.1103/PhysRevApplied.3.034009. Accessed 2022-06-28
  * (47) Kim, J., Sinha, J., Hayashi, M., Yamanouchi, M., Fukami, S., Suzuki, T., Mitani, S., Ohno, H.: Layer thickness dependence of the current-induced effective field vector in Ta|CoFeB|MgO. Nature Materials 12(3), 240–245 (2013). https://doi.org/10.1038/nmat3522. Accessed 2022-06-28
  * (48) Ziętek, S., Mojsiejuk, J., Grochot, K., Łazarski, S., Skowroński, W., Stobiecki, T.: Numerical model of harmonic Hall voltage detection for spintronic devices. Physical Review B 106(2), 024403 (2022). https://doi.org/10.1103/PhysRevB.106.024403. Accessed 2022-07-13
  * (49) Skowroński, W., Grochot, K., Rzeszut, P., Łazarski, S., Gajoch, G., Worek, C., Kanak, J., Stobiecki, T., Langer, J., Ocker, B., Vafaee, M.: Angular Harmonic Hall Voltage and Magnetoresistance Measurements of Pt/FeCoB and Pt-Ti/FeCoB Bilayers for Spin Hall Conductivity Determination. IEEE Transactions on Electron Devices 68(12), 6379–6385 (2021). https://doi.org/%****␣sn-article.bbl␣Line␣1075␣****10.1109/TED.2021.3122999. Accessed 2022-06-28
  * (50) Grochot, K., Karwacki, Ł., Łazarski, S., Skowroński, W., Kanak, J., Powroźnik, W., Kuświk, P., Kowacz, M., Stobiecki, F., Stobiecki, T.: Current-Induced Magnetization Switching of Exchange-Biased NiO Heterostructures Characterized by Spin-Orbit Torque. Physical Review Applied 15(1), 014017 (2021). https://doi.org/10.1103/PhysRevApplied.15.014017. Accessed 2022-03-23
  * (51) Lee, K.-S., Lee, S.-W., Min, B.-C., Lee, K.-J.: Threshold current for switching of a perpendicular magnetic layer induced by spin Hall effect. Applied Physics Letters 102(11), 112410 (2013). https://doi.org/10.1063/1.4798288. Accessed 2022-03-23
  * (52) Romera, M., Talatchian, P., Tsunegi, S., Yakushiji, K., Fukushima, A., Kubota, H., Yuasa, S., Cros, V., Bortolotti, P., Ernoult, M., Querlioz, D., Grollier, J.: Binding events through the mutual synchronization of spintronic nano-neurons. Nature Communications 13(1), 883 (2022). https://doi.org/10.1038/s41467-022-28159-1. Accessed 2022-06-28
  * (53) Locatelli, N., Cros, V., Grollier, J.: Spin-torque building blocks. Nature Materials 13(1), 11–20 (2014). https://doi.org/10.1038/nmat3823. Accessed 2021-11-15
  * (54) Camsari, K.Y., Ganguly, S., Datta, S.: Modular Approach to Spintronics. Scientific Reports 5(1), 10571 (2015). https://doi.org/%****␣sn-article.bbl␣Line␣1175␣****10.1038/srep10571. Accessed 2022-06-07
  * (55) Skowroński, W., Chęciński, J., Ziętek, S., Yakushiji, K., Yuasa, S.: Microwave magnetic field modulation of spin torque oscillator based on perpendicular magnetic tunnel junctions. Scientific Reports 9(1), 19091 (2019). https://doi.org/10.1038/s41598-019-55220-9. Accessed 2021-09-09