# Euler scheme for approximation of solution of nonlinear ODEs under inexact
information

Natalia Czyżewska  ,  Paweł M. Morkisz  and  Paweł Przybyłowicz  AGH
University of Science and Technology, Faculty of Applied Mathematics, al.
Mickiewicza 30, 30-059 Kraków, Poland<EMAIL_ADDRESS>corresponding
author AGH University of Science and Technology, Faculty of Applied
Mathematics, al. Mickiewicza 30, 30-059 Kraków, Poland<EMAIL_ADDRESS>AGH University of Science and Technology, Faculty of Applied Mathematics, al.
Mickiewicza 30, 30-059 Kraków, Poland<EMAIL_ADDRESS>
###### Abstract.

We investigate error of the Euler scheme in the case when the right-hand side
function of the underlying ODE satisfies nonstandard assumptions such as local
one-side Lipschitz condition and local Hölder continuity. Moreover, we assume
two cases in regards to information availability: exact and noisy with respect
to the right-hand side function. Optimality analysis of the Euler scheme is
also provided. Lastly, we present the results of some numerical experiments.

Mathematics Subject Classification: 65L05, 65L70

###### Key words and phrases:

noisy information, ODEs, nonstandard assumptions, local conditions, one-side
Lipschitz condition, Hölder continuity, Euler scheme

###### Contents

  1. 1 Introduction
  2. 2 Preliminaries
  3. 3 Euler scheme under exact information
  4. 4 Euler scheme under noisy information
  5. 5 Lower bounds and optimality of the Euler scheme in the worst-case setting
  6. 6 Numerical results
    1. 6.1 Additive example
    2. 6.2 Multiplicative example
  7. 7 Appendix

##  1\. Introduction

In this paper we investigate error of the Euler scheme for the following
initial-value problems

$\displaystyle\left\\{\begin{array}[]{ll}z^{\prime}(t)=g(t,z(t)),&t\in[a,b],\\\
z(a)=\xi,\end{array}\right.$ (1.3)

where $-\infty<a<b<+\infty$, $\xi\in\mathbb{R}^{d}$ and a right-hand side
function $g:[a,b]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$, $d\in\mathbb{Z}_{+}$.
Comparing to the classical assumptions such us global Lipschitz condition we
assume that the right-hand side function $f$ is continuous and of at most
linear growth but it satisfies the one-sided Lipschitz condition together with
the Hölder condition only locally.

Numerical approximation of solutions of ordinary differential equations (ODEs)
has drawn the attention of numerous scientists [4, 5, 8, 10, 11, 13].
Furthermore, solving ODE problems with noisy information was studied in [2, 3,
14, 16, 17, 25, 26, 29, 31, 32]. Also, approximation of the solutions of
partial differential equations or stochastic differential equations with both
exact and inexact information was a subject of extensive study, [1, 21, 22,
33, 35, 36] and [20, 24, 19, 27, 28] respectively. Nevertheless, when it comes
to the nonstandard assumptions about the right-hand side function of the
initial-value problem, only the existence and uniqueness of the solutions were
examined, see [9]. There is a lack of research in the literature about
numerical methods for such cases under nonstandard assumptions.

If we take a closer look at the computation on a digital computer, we quickly
conclude that the computer uses only a finite set of numbers because of the
finite memory. So all computer calculations are burdened with errors, e.g.,
from rounding. Also, the novel CPUs and GPUs have an option to speed up the
computations when there is lower precision of the processed numbers. From the
mathematical point of view, lowering precision can be interpreted as adding
noisy information with the control for the relative error of the noise. Hence,
all the applied numerical methods should have a strict error analysis with
noisy information setting.

In this paper, we deal with nonstandard assumptions about the right-hand side
function $g$, and we allow both exact and inexact information about the
initial value $\xi$ and values of the right-hand side function $g$. We use
slightly different assumptions in comparison to the auxiliary ordinary
differential equation in the [7], i.e., there we assume locally Hölder
continuity and satisfying global one-side Lipschitz conditions and now we
consider locally Hölder continuity and satisfying one-side Lipschitz condition
only in the certain neighborhood. This weakening of the assumptions was caused
by the necessity that came from a phase change model of metallic materials [6,
7].

According to our best knowledge, there is a lack of literature on numerical
solving of initial-value problems under nonstandard assumptions with noisy
information. We focus on the most basic explicit Euler scheme. We believe that
this paper might be helpful for further research on implicit numerical
methods, higher-order methods, or even numerical methods for delay
differential equations with noisy information under nonstandard assumptions on
the right-hand side function.

To do so, we start with an error analysis for the classical explicit Euler
scheme under nonstandard assumptions for the right-hand side function $f$
mentioned before, i.e., local Hölder continuity and one-side Lipschitz
condition (Theorem 3.1). Next, we introduce a corruption for the right-hand
side function and then we do an error analysis for the Euler scheme (Theorem
4.1). In both cases, we show the dependence of the error on Hölder exponents
of the right-hand side function. Furthermore, we establish lower bounds for
all algorithms based on exact and inexact information in the certain class of
right-hand side functions (Proposition 5.1 and Proposition 5.2 respectively).
Consequently, we can provide a result concerning the optimality of the Euler
scheme in a particular case (Theorem 5.3). The final contribution is
performing numerical experiments confirming our theoretical results.

The paper is organized as follows. The problem posing with the nonstandard
assumptions is given in Section 2. Section 3 on exact information includes a
definition of the Euler scheme and a result about the upper bound of the error
of this scheme under mentioned nonstandard assumptions. The meaning of the
corrupting noise is explained in Section 4, along with the definition of the
Euler scheme in case of inexact information and the upper bound of the Euler
scheme’s error with only noisy available standard information about the right-
hand side function. Finally, Section 5 contains the definition of the setting
in information-based complexity spirit (see [34]), the lower bounds of the
Euler scheme, and the optimality of this scheme as the main result of this
paper. Section 6 is dedicated to numerical experiments performed using our
implementation of the Euler scheme. Last but not least, we put all the
auxiliary results (including, i.e. a proof of the existence and uniqueness of
the solutions of (1.3)) in the Appendix.

##  2\. Preliminaries

For $x,y\in\mathbb{R}^{d}$ we take $\displaystyle{\langle
x,y\rangle=\sum\limits_{k=1}^{d}x_{k}y_{k}}$ and $\|x\|=\langle
x,x\rangle^{1/2}$. For $R\in[0,+\infty]$ and $x\in\mathbb{R}^{d}$, we denote
by $B_{x}(R)=\\{y\in\mathbb{R}^{d}\ |\ \|y-x\|\leq R\\}$. Of course
$B_{x}(+\infty)=\mathbb{R}^{d}$. Moreover, for $H\in\mathbb{R}$ we denote by
$H^{+}=\max\\{0,H\\}$.

We assume that the right-hand side function
$g:[a,b]\times\mathbb{R}^{d}\to\mathbb{R}^{d}$ satisfies the following
conditions:

  * (G1)

$g\in C([a,b]\times\mathbb{R}^{d};\mathbb{R}^{d})$.

  * (G2)

There exists $K\in(0,+\infty)$ such that for all
$(t,y)\in[a,b]\times\mathbb{R}^{d}$

$\|g(t,y)\|\leq K(1+\|y\|).$

  * (G3)

For every $U\in\textrm{comp}(\mathbb{R}^{d})$ there exists
$H_{U}\in\mathbb{R}$ such that for all $t\in[a,b]$, $y_{1},y_{2}\in U$

$\langle y_{1}-y_{2},g(t,y_{1})-g(t,y_{2})\rangle\leq
H_{U}\|y_{1}-y_{2}\|^{2}.$

  * (G4)

There exist $\alpha,\beta\in(0,1]$ such that for every
$U\in\textrm{comp}(\mathbb{R}^{d})$ there exists $L_{U}\in(0,\infty)$ such
that for all $t_{1},t_{2}\in[a,b]$, $y_{1},y_{2}\in U$

$\|g(t_{1},y_{1})-g(t_{2},y_{2})\|\leq
L_{U}\Bigl{(}|t_{1}-t_{2}|^{\alpha}+\|y_{1}-y_{2}\|^{\beta}\Bigr{)}.$

The assumption (G2) is the classical global linear growth condition. However,
instead of a global Lipschitz condition, we assume that the right-hand side
function $g$ satisfies locally one-side Lipschiz assumption (G3). It follows
from Lemma 7.2 that under the assumptions (G1)-(G3) the equation (1.3) has a
unique solution $z\in C^{1}([a,b];\mathbb{R}^{d})$. The additional assumption
(G4), called the $(\alpha,\beta)$-local Hölder condition, is needed for
obtaining an upper error bound for the Euler scheme that depends on
$(\alpha,\beta)$. We stress that under the assumptions (G1)-(G4) the right-
hand side function might be highly nonlinear, see, for example, Section 6.

In the next section we investigate error of the classical Euler scheme under
such nonstandard assumptions. Moreover, we consider the cases when exact and
inexact information about the right-hand side function $g$ is available.

##  3\. Euler scheme under exact information

Let $n\in\mathbb{Z}_{+}$, $t_{k}=a+kh$ for $k=0,1,\dots,n$ where
$h=\frac{b-a}{n}$. The Euler scheme is defined as follows

$\displaystyle\left\\{\begin{array}[]{ll}y_{0}=\xi,\\\ y_{k+1}=y_{k}+h\cdot
g(t_{k},y_{k})\quad\hbox{for}&k=0,1,\dots,n-1,\end{array}\right.$ (3.3)

The approximation of $z=z(t)$ in $[a,b]$ is given by a piecewise linear
function $l_{n}:[a,b]\to\mathbb{R}^{d}$, defined in the following way

$l_{n}(t)=l_{k,n}(t)\ \hbox{for}\ t\in[t_{k},t_{k+1}],\ k=0,1,\ldots,n-1,$
(3.4)

where

$l_{k,n}(t)=y_{k}+(t-t_{k})\cdot g(t_{k},y_{k}).$ (3.5)

The aim of this section is to prove the following result.

######  Theorem 3.1.

Let $\xi\in\mathbb{R}^{d}$ and let us assume that $g$ satisfies (G1)-(G4).
Then there exists $\bar{C}_{1}\in[0,\infty)$ such that for all
$n\in\mathbb{N}$

$\sup\limits_{a\leq t\leq
b}\|z(t)-l_{n}(t)\|\leq\bar{C}_{1}(1+\|\xi\|)(h^{\alpha}+h^{\beta}).$ (3.6)

###### Proof.

By the assumption (G2) we have that for all $k=0,1,\ldots,n-1$

$\|y_{k+1}\|\leq\|y_{k}\|+h\|g(t_{k},y_{k})\|\leq(1+hK)\|y_{k}\|+hK.$

Hence, by the discrete version of Gronwall’s lemma we get that for all
$k=0,1,\ldots,n$ that

$\|y_{k}\|\leq(1+Kh)^{k}\|\xi\|+(1+Kh)^{k}-1\leq e^{K(b-a)}(1+\|\xi\|),$ (3.7)

and therefore

$\max_{n\in\mathbb{Z}_{+}}\max_{0\leq k\leq n}\|y_{k}\|\leq C_{1}(1+\|\xi\|),$
(3.8)

where $C_{1}=C_{1}(a,b,K)=e^{K(b-a)}$. Note that $C_{1}\leq C_{0}$, where
$C_{0}$ is defined in Lemma 7.2. Now for $k=0,1,\ldots,n-1$ we consider the
following local ODE

$z_{k}^{\prime}(t)=g(t,z_{k}(t)),\quad t\in[t_{k},t_{k+1}],\quad
z_{k}(t_{k})=y_{k}.$ (3.9)

By Lemma 7.2 there exists the unique solution
$z_{k}:[t_{k},t_{k+1}]\to\mathbb{R}^{d}$ of (3.9). From the assumption (G2)
and by (3.8) we get for all $t\in[t_{k},t_{k+1}]$, $k=0,1,\ldots,n-1$ that

$\|z_{k}(t)\|\leq
e^{K(b-a)}(1+\|\xi\|)+K(b-a)+K\int\limits_{t_{k}}^{t}\|z_{k}(s)\|ds,$ (3.10)

and, by the Gronwall’s lemma,

$\max_{n\in\mathbb{Z}_{+}}\max_{k=0,1,\dots,n-1}\sup_{t_{k}\leq t\leq
t_{k+1}}\|z_{k}(t)\|\leq K_{1},$ (3.11)

where $K_{1}=K_{1}(a,b,K,\xi)=C_{2}(\|\xi\|+1)$, and
$C_{2}=C_{2}(a,b,K)=e^{K(b-a)}[e^{K(b-a)}+K(b-a)]$. We have that $C_{2}\geq
C_{0}$. Moreover, by (G2) and (3.11) we have for $t\in[t_{k},t_{k+1}]$,
$k=0,1,\ldots,n-1$

$\|z_{k}(t)-y_{k}\|\leq hK\Bigl{(}1+\sup\limits_{t_{k}\leq t\leq
t_{k+1}}\|z_{k}(t)\|\Bigr{)}\leq hC_{3}(1+\|\xi\|),$ (3.12)

where $C_{3}=C_{3}(a,b,K)=K(1+C_{2})=K(1+e^{K(b-a)}[e^{K(b-a)}+K(b-a)])$.
Therefore, by (3.8), (3.11), (7.1) and the fact that $z_{k}(t_{k})=y_{k}$ we
have

$y_{k},z_{k}(t),z(t)\in B_{0}(C_{2}(1+\|\xi\|))$ (3.13)

for all $t\in[t_{k},t_{k+1}]$, $k=0,1,\ldots,n-1$, $n\in\mathbb{N}$. This and
(G4) imply that for $U=B_{0}(C_{2}(1+\|\xi\|))$ that there exists
$L_{1}\in(0,+\infty)$ such that for all $k=0,1,\ldots,n-1$, $n\in\mathbb{N}$

$\displaystyle\|z_{k}(t_{k+1})-y_{k+1}\|\leq\int_{t_{k}}^{t_{k+1}}\|g(s,z_{k}(s))-g(t_{k},y_{k})\|ds$
$\displaystyle\leq
L_{1}\int_{t_{k}}^{t_{k+1}}(|s-t_{k}|^{\alpha}+\|z_{k}(s)-y_{k}\|^{\beta})ds\leq
L_{1}C_{4}(1+\|\xi\|)h(h^{\alpha}+h^{\beta}),$ (3.14)

where $C_{4}=C_{4}(a,b,K,\beta)=\max\\{1,2C_{3}^{\beta}\\}$. Furthermore, for
$k=0,1,\dots,n-1$ we have

$\|z(t_{k+1})-y_{k+1}\|\leq\|z(t_{k+1})-z_{k}(t_{k+1})\|+\|z_{k}(t_{k+1})-y_{k+1}\|,$
(3.15)

and for $t\in[t_{k},t_{k+1}]$, $k=0,1,\ldots,n-1$ it holds that

$\frac{d}{dt}\|z(t)-z_{k}(t)\|^{2}=2\langle
z(t)-z_{k}(t),g(t,z(t))-g(t,z_{k}(t))\rangle.$ (3.16)

Therefore, by (G3) applied to $U=B_{1}$ we receive that there exists
$H_{1}\in\mathbb{R}$ such that

$\frac{d}{dt}\|z(t)-z_{k}(t)\|^{2}\leq 2H_{1}\|z(t)-z_{k}(t)\|^{2}\\\ \leq
2H_{1}^{+}\|z(t)-z_{k}(t)\|^{2}$

for all $t\in[t_{k},t_{k+1}]$, $k=0,1,\ldots,n-1$, $n\in\mathbb{N}$. Since
$z_{k}(t_{k})=y_{k}$, we get

$\|z(t)-z_{k}(t)\|^{2}\leq\|z(t_{k})-y_{k}\|^{2}+2H^{+}_{1}\int\limits_{t_{k}}^{t}\|z(s)-z_{k}(s)\|^{2}ds,$
(3.17)

for $t\in[t_{k},t_{k+1}]$ and, by the Gronwall’s lemma,

$\|z(t)-z_{k}(t)\|\leq e^{H_{1}^{+}(t-t_{k})}\|z(t_{k})-y_{k}\|.$ (3.18)

Combining (3), (3.15), (3.18) we receive for $k=0,1,\ldots,n-1$

$\|z(t_{k+1})-y_{k+1}\|\leq
e^{H_{1}^{+}h}\|z(t_{k})-y_{k}\|+L_{1}C_{4}(1+\|\xi\|)h(h^{\alpha}+h^{\beta}).$
(3.19)

Denoting by $e_{k}=z(t_{k})-y_{k}$ we arrive at

$\|e_{k+1}\|\leq
e^{H_{1}^{+}h}\|e_{k}\|+L_{1}C_{4}(1+\|\xi\|)h(h^{\alpha}+h^{\beta}),$

with $e_{0}=z(a)-y_{0}=0$. If $H_{1}^{+}>0$ then by Gronwall’s lemma we obtain

$\|e_{k}\|\leq\frac{e^{H_{1}^{+}(b-a)}-1}{H_{1}^{+}}L_{1}C_{4}(1+\|\xi\|)(h^{\alpha}+h^{\beta}),$
(3.20)

for $k=0,1,\ldots,n$, where $C_{4}=C_{4}(a,b,K,\beta)$. In the particular
case, when $H_{1}^{+}=0$ we simply have that for $k=0,1,\ldots,n$

$\|e_{k}\|\leq
kL_{1}C_{4}(1+\|\xi\|)(h^{\alpha}+h^{\beta})\leq(b-a)L_{1}C_{4}(1+\|\xi\|)(h^{\alpha}+h^{\beta}).$
(3.21)

For $t\in[t_{k},t_{k+1}]$, $k=0,1,\ldots,n$ we have by (3.18) that

$\displaystyle\|z(t)-l_{n}(t)\|=\|z(t)-l_{k,n}(t)\|\leq\|z(t)-z_{k}(t)\|+\|z_{k}(t)-l_{k,n}(t)\|$
$\displaystyle\leq e^{H_{1}^{+}h}\|e_{k}\|+\|z_{k}(t)-l_{k,n}(t)\|,$ (3.22)

where by (3)

$\|z_{k}(t)-l_{k,n}(t)\|\leq\int\limits_{t_{k}}^{t_{k+1}}\|g(s,z_{k}(s))-g(t_{k},y_{k})\|ds\leq
L_{1}C_{4}(1+\|\xi\|)h(h^{\alpha}+h^{\beta}).$ (3.23)

Combining (3), (3.23), (3.20), (3.21) we get the thesis. ∎

##  4\. Euler scheme under noisy information

In this section we show the upper error bound for the error of the Euler
scheme in the case of inexact information. Namely, we consider the corrupted
right-hand side functions of the following form

$\tilde{g}(t,y)=g(t,y)+\delta_{g}(t,y),\quad(t,y)\in[a,b]\times\mathbb{R}^{d},$
(4.1)

where $\delta_{g}:[a,b]\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ is a
corrupting function that is Borel measurable and

$\|\delta_{g}(t,y)\|\leq\delta(1+\|y\|)$ (4.2)

for all $(t,y)\in[a,b]\times\mathbb{R}^{d}$ and $\delta\in[0,1]$. We refer to
$\delta$ as to precision parameter.

In the case of inexact evaluations of $g$, the Euler scheme has the following
form

$\displaystyle\left\\{\begin{array}[]{ll}\tilde{y}_{0}=\tilde{\xi},\\\
\tilde{y}_{k+1}=\tilde{y}_{k}+h\cdot\tilde{g}(t_{k},\tilde{y}_{k}),&k=0,1,\dots,n-1,\end{array}\right.$
(4.5)

where $\tilde{\xi}\in B_{\xi}(\delta)$. The approximation of $z=z(t)$ in
$[a,b]$ is given by a piecewise linear function
$\tilde{l}_{n}:[a,b]\to\mathbb{R}^{d}$, defined as

$\tilde{l}_{n}(t)=\tilde{l}_{k,n}(t)\ \hbox{for}\ t\in[t_{k},t_{k+1}],\
k=0,1,\ldots,n-1,$ (4.6)

where

$\tilde{l}_{k,n}(t)=\tilde{y}_{k}+(t-t_{k})\cdot\tilde{g}(t_{k},\tilde{y}_{k}).$
(4.7)

Of course when the information is exact (i.e., $\delta=0$) then
$\tilde{y}_{k}=y_{k}$ for $k=0,1,\ldots,n$, and $l_{n}\equiv\tilde{l}_{n}$.
Moreover, the Euler algorithm $\tilde{l}_{n}$ uses $n$ (noisy) evaluations of
$\tilde{g}$.

In the theorem below we show the upper error bound for the Euler scheme in the
case of inexact information about the function $g$. We stress that we have
different bounds in the local Lipschitz case ($\beta=1$) and in the local
Hölder case ($\beta\in(0,1)$). This is due to the fact that locally Lipschitz
continuous function (i.e. $\beta=1$ in (G4)) satisfies also locally one-sided
Lipschitz condition (G3). This is not the case when $g$ is only locally Hölder
continuous. We stress that the condition (G3) is crucial for obtaining the
upper error bound in the case when $\beta\in(0,1)$, see also Remark 4.2.

######  Theorem 4.1.

Let $\xi\in\mathbb{R}^{d}$, let $g$ satisfy (G1)-(G4), and let $\delta_{g}$ be
a Borel measurable corrupting function satisfying (4.2). Then there exists
$\bar{C}_{2}\in[0,\infty)$ such that for all $n\in\mathbb{N}$, $n\geq\lceil
2(b-a)\rceil+1$, $\delta\in[0,1]$, $\tilde{\xi}\in B_{\xi}(\delta)$ the
following holds:

  * •

if $\beta=1$ then

$\sup\limits_{a\leq t\leq
b}\|z(t)-\tilde{l}_{n}(t)\|\leq\bar{C}_{2}(1+\|\xi\|)(h^{\alpha}+\delta),$
(4.8)

  * •

if $\beta\in(0,1)$ then

$\sup\limits_{a\leq t\leq
b}\|z(t)-\tilde{l}_{n}(t)\|\leq\bar{C}_{2}(1+\|\xi\|)\Bigl{(}h^{\alpha}+h^{\beta}+\mathbf{1}_{(0,1]}(\delta)(h^{1/2}+\delta)\Bigr{)}.$
(4.9)

###### Proof.

By (G2) we have for all $k=0,1,\dots,n-1$ that

$\|\tilde{y}_{k+1}\|\leq\|\tilde{y}_{k}\|+h(K+1)(1+\|\tilde{y}_{k}\|)\\\
\leq(1+h(K+1))\|\tilde{y}_{k}\|+h(K+1).$

Therefore, by a discrete version of Gronwall’s lemma we get for
$k=0,1,\ldots,n-1$, $n\in\mathbb{N}$, $\delta\in[0,1]$ that

$\displaystyle\|\tilde{y}_{k}\|\leq(1+h(K+1))^{k}\|\tilde{y}_{0}\|+(1+h(K+1))^{k}-1$
$\displaystyle\leq e^{(b-a)(K+1)}(\|\tilde{\xi}\|+1)\leq K_{2}(1+\|\xi\|)$
(4.10)

where $K_{2}=K_{2}(a,b,K)=2e^{(b-a)(K+1)}$. Since $K_{2}\geq C_{1}$, by
(3.8),(4) we have

$y_{k},\tilde{y}_{k}\in B_{0}(K_{2}(1+\|\xi\|))$ (4.11)

for all $k=0,1,\ldots,n-1$, $n\in\mathbb{N}$, $\delta\in[0,1]$. Let us denote
$e_{k}=\tilde{y}_{k}-y_{k}$ for $k=0,1,\dots,n$. Then we have that
$\|e_{0}\|=\|\tilde{y}_{0}-y_{0}\|=\|\tilde{\xi}-\xi\|\leq\delta$ and for
$k=0,1,\dots,n-1$

$e_{k+1}=e_{k}+h\left[g(t_{k},\tilde{y}_{k})-g(t_{k},y_{k})\right]+h\delta_{g}(t_{k},\tilde{y}_{k}).$
(4.12)

We consider two cases, when $\beta=1$ and $\beta\in(0,1)$.

If $\underline{\beta=1}$ then (G4) implies that for
$U=B_{0}(K_{2}(1+\|\xi\|))$ there exist $L_{2}\in(0,+\infty)$ such that for
$k=0,1,\ldots,n-1$, $n\in\mathbb{N}$, $\delta\in[0,1]$

$\|e_{k+1}\|\leq(1+hL_{2})\|e_{k}\|+h\delta(1+\|\xi\|)(1+K_{2}).$ (4.13)

By the discrete version of the Gronwall’s lemma we get for all
$n\in\mathbb{N}$, $\delta\in[0,1]$

$\max\limits_{0\leq k\leq
n}\|y_{k}-\tilde{y}_{k}\|\leq(1+\|\xi\|)\Bigl{(}e^{L_{2}(b-a)}+\frac{e^{L_{2}(b-a)}-1}{L_{2}}(1+K_{2})\Bigr{)}\delta.$
(4.14)

Hence, by Theorem 3.1 we get that there exists $\tilde{C}_{1}\in[0,+\infty)$
such that for all $t\in[t_{k},t_{k+1}]$, $k=0,1,\ldots,n$, $n\in\mathbb{N}$,
$\delta\in[0,1]$

$\|z(t)-\tilde{l}_{n}(t)\|\leq\bar{C}_{1}(1+\|\xi\|)h^{\alpha}+\|\tilde{l}_{k,n}(t)-l_{k,n}(t)\|,$
(4.15)

and, by (4.14),

$\displaystyle\|\tilde{l}_{k,n}(t)-l_{k,n}(t)\|\leq(1+hL_{2})\|y_{k}-\tilde{y}_{k}\|+h\delta(1+\|\tilde{y}_{k}\|)$
$\displaystyle\leq\tilde{C}_{1}(1+\|\xi\|)\delta+h\delta(1+\|\xi\|)(1+K_{2}).$
(4.16)

From (4.15), (4.16) we get (4.8).

If $\underline{\beta\in(0,1)}$ then we have to proceed in a different way, see
also Remark 4.2. From (4.12) we have for $k=0,1,\ldots,n-1$ that

$\|e_{k+1}-h\delta_{g}(t_{k},\tilde{y}_{k})\|^{2}=\|e_{k}+h\left[g(t_{k},\tilde{y}_{k})-g(t_{k},y_{k})\right]\|^{2}.$
(4.17)

We have

$\displaystyle\|e_{k+1}-h\delta_{g}(t_{k},\tilde{y}_{k})\|^{2}=\|e_{k+1}\|^{2}-2h\langle
e_{k+1},\delta_{g}(t_{k},\tilde{y}_{k})\rangle+h^{2}\|\delta_{g}(t_{k},\tilde{y}_{k})\|^{2}$
$\displaystyle\geq\|e_{k+1}\|^{2}-2h\langle
e_{k+1},\delta_{g}(t_{k},\tilde{y}_{k})\rangle.$ (4.18)

On the other hand, (4.11) and (G3), (G4) imply that for
$U=B_{0}(K_{2}(1+\|\xi\|))$ there exist $L_{2}\in(0,+\infty)$,
$H_{2}\in\mathbb{R}$ such that for $k=0,1,\ldots,n-1$, $n\in\mathbb{N}$,
$\delta\in[0,1]$

$\displaystyle\|e_{k}+h\left[g(t_{k},\tilde{y}_{k})-g(t_{k},y_{k})\right]\|^{2}$
$\displaystyle=\|e_{k}\|^{2}+2h\left\langle
e_{k},g(t_{k},\tilde{y}_{k})-g(t_{k},y_{k})\right\rangle+h^{2}\|g(t_{k},\tilde{y}_{k})-g(t_{k},y_{k})\|^{2}$
$\displaystyle\leq\|e_{k}\|^{2}+2hH_{2}^{+}\|e_{k}\|^{2}+h^{2}L_{2}^{2}\|e_{k}\|^{2\beta}\leq(1+2hH_{2}^{+})\|e_{k}\|^{2}+2h^{2}L_{2}^{2}(1+\|e_{k}\|^{2})$
$\displaystyle\leq(1+K_{3}h)\|e_{k}\|^{2}+2h^{2}L^{2},$ (4.19)

where $K_{3}=2(H_{2}^{+}+L_{2}^{2})$. (We have also used the fact that
$x^{2\beta}=(x^{2})^{\beta}\leq 1+x^{2}$ for all $x\in\mathbb{R}$,
$\beta\in(0,1]$.) Combining (4.17), (4.18) and (4.19) we get

$\|e_{k+1}\|^{2}\leq(1+K_{3}h)\|e_{k}\|^{2}+2h^{2}L_{2}^{2}+2h\langle
e_{k+1},\delta_{g}(t_{k},\tilde{y}_{k})\rangle.$ (4.20)

Moreover,

$\langle
e_{k+1},\delta_{g}(t_{k},\tilde{y}_{k})\rangle\leq\|e_{k+1}\|\cdot\|\delta_{g}(t_{k},\tilde{y}_{k})\|\leq\frac{1}{2}\|e_{k+1}\|^{2}+\delta^{2}(1+K_{2}^{2})(1+\|\xi\|)^{2}.$
(4.21)

Combining (4.20) with (4.21) we end up with

$\|e_{k+1}\|^{2}\leq(1+K_{3}h)\|e_{k}\|^{2}+2h^{2}L_{2}^{2}+h\|e_{k+1}\|^{2}+2h\delta^{2}(1+K_{2}^{2})(1+\|\xi\|)^{2}.$
(4.22)

Since $h\in(0,\frac{1}{2})$, we have $0<\frac{1}{1-h}\leq 1+2h\leq 2$ and for
all $k=0,1,\dots,n-1$

$\displaystyle\|e_{k+1}\|^{2}\leq(1+2h)(1+K_{3}h)\|e_{k}\|^{2}+K_{4}(1+\|\xi\|)^{2}h(h+\delta^{2})$
$\displaystyle\leq(1+K_{5}h)\|e_{k}\|^{2}+K_{4}(1+\|\xi\|)^{2}h(h+\delta^{2}),$
(4.23)

where $K_{4}=4\max\\{L_{2}^{2},1+K_{2}^{2}\\}$, $K_{5}=2+3K_{3}>0$. By the
discrete version of the Gronwall’s lemma we obtain for all $k=0,1,\dots,n$

$\displaystyle\|e_{k}\|^{2}\leq(1+K_{5}h)^{k}\|e_{0}\|^{2}+\frac{(1+K_{5}h)^{k}-1}{K_{5}}K_{4}(1+\|\xi\|)^{2}(h+\delta^{2})$
$\displaystyle\leq
e^{K_{5}(b-a)}\delta^{2}+\frac{e^{K_{5}(b-a)}-1}{K_{5}}K_{4}(1+\|\xi\|)^{2}(h+\delta^{2})$
$\displaystyle\leq K_{6}(1+\|\xi\|)^{2}(h+\delta^{2}),$ (4.24)

where $K_{6}=2\max\\{e^{K_{5}(b-a)},\frac{e^{K_{5}(b-a)}-1}{K_{5}}K_{4}\\}$.
Therefore, for all $n\in\mathbb{N}$, $\delta\in[0,1]$

$\max_{k=0,1,\dots,n}\|y_{k}-\tilde{y}_{k}\|\leq
K_{7}(1+\|\xi\|)(h^{\frac{1}{2}}+\delta),$ (4.25)

with $K_{7}=K_{6}^{1/2}$. Again, by Theorem 3.1 we get that there exists
$\tilde{C}_{2},\tilde{C}_{3}\in[0,+\infty)$ such that for all
$t\in[t_{k},t_{k+1}]$, $k=0,1,\ldots,n$, $n\in\mathbb{N}$, $\delta\in[0,1]$

$\|z(t)-\tilde{l}_{n}(t)\|\leq\bar{C}_{1}(1+\|\xi\|)(h^{\alpha}+h^{\beta})+\|\tilde{l}_{k,n}(t)-l_{k,n}(t)\|,$
(4.26)

and, by (4.25),

$\displaystyle\|\tilde{l}_{k,n}(t)-l_{k,n}(t)\|\leq\|y_{k}-\tilde{y}_{k}\|+hL_{2}\|y_{k}-\tilde{y}_{k}\|^{\beta}+h\delta(1+\|\tilde{y}_{k}\|)$
$\displaystyle\leq\tilde{C}_{2}(1+\|\xi\|)(h^{1/2}+\delta+h^{1+\frac{\beta}{2}}+h\delta^{\beta})+h\delta(1+\|\xi\|)(1+K_{2})$
$\displaystyle\leq\tilde{C}_{3}(1+\|\xi\|)(h^{1/2}+\delta).$ (4.27)

Finally, from (4.26), (4) we get (4.9). ∎

######  Remark 4.2.

Note that if $\beta\in(0,1)$ then we can only get that

$\|e_{k+1}\|\leq\|e_{k}\|+hL_{2}\|e_{k}\|^{\beta}+h\delta(1+\|\xi\|)(1+K_{2})$
(4.28)

for $k=0,1,\ldots,n-1$. Since $\|e_{k}\|^{\beta}\leq 1+\|e_{k}\|$, it holds

$\|e_{k+1}\|\leq(1+hL_{2})\|e_{k}\|+hL_{2}+h\delta(1+\|\xi\|)(1+K_{2}),$
(4.29)

which in turns implies that

$\max\limits_{0\leq k\leq n}\|e_{k}\|\leq
e^{L_{2}(b-a)}\delta+\frac{e^{L_{2}(b-a)}-1}{L_{2}}(L_{2}+\delta)(1+\|\xi\|)(1+K_{2})=O(1),$
(4.30)

as $\delta\to 0^{+}$. This estimate is obviously too rough.

##  5\. Lower bounds and optimality of the Euler scheme in the worst-case
setting

In this section we provide optimality analysis of the Euler algorithm in the
case of inexact information about the right-hand side function. We settle our
problem in the worst-case model of error (see [34]), which, in particular,
allows us to investigate lower error bounds.

Let $\alpha,\beta\in(0,1]$, $K,L\in(0,+\infty)$, $H\in\mathbb{R}$,
$R\in[0,+\infty]$. We consider the class
$F_{R}^{\alpha,\beta}=F^{\alpha,\beta}_{R}(a,b,d,K,H,L)$ of pairs $(\xi,g)$
that satisfy the following conditions:

  * (A0)

$\|\xi\|\leq K$,

  * (A1)

$g\in C([a,b]\times\mathbb{R}^{d};\mathbb{R}^{d})$,

  * (A2)

$\|g(t,y)\|\leq K(1+\|y\|)$ for all $(t,y)\in[a,b]\times\mathbb{R}^{d}$,

  * (A3)

$\langle y_{1}-y_{2},g(t,y_{1})-g(t,y_{2})\rangle\leq H\|y_{1}-y_{2}\|^{2}$
for all $t\in[a,b]$, $y_{1},y_{2}\in B_{0}(R)$,

  * (A4)

$\|g(t_{1},y_{1})-g(t_{2},y_{2})\|\leq
L\Bigl{(}|t_{1}-t_{2}|^{\alpha}+\|y_{1}-y_{2}\|^{\beta}\Bigr{)}$ for all
$t_{1},t_{2}\in[a,b]$, $y_{1},y_{2}\in B_{0}(R)$.

We refer to $\alpha,\beta\in(0,1]$, $-\infty<a<b<+\infty$,
$d\in\mathbb{Z}_{+}$, $K,L\in[0,\infty)$, $H\in\mathbb{R}$, $R\in[0,+\infty]$
as to the parameters of the class. Except for $a,b,d$ the parameters are not
known and will not be used by the defined algorithms.

The class $F^{\alpha,\beta}_{\infty}$ consists of functions $g$ that satisfy
conditions (A3), (A4) globally on the whole domain
$[a,b]\times\mathbb{R}^{d}$, while functions from $F^{\alpha,\beta}_{0}$ might
not satisfy the one-side Lipschitz condition and the $(\alpha,\beta)$-Hölder
condition even locally. We have that $F_{\infty}^{\alpha,\beta}\subset
F_{R_{1}}^{\alpha,\beta}\subset F_{R_{2}}^{\alpha,\beta}\subset
F_{0}^{\alpha,\beta}$ for all $0\leq R_{2}\leq R_{1}\leq+\infty$.

We assume that an algorithm may use exact or noisy evaluations of the function
$g$. Consequently, we can have access to the function $g$ through its possibly
perturbed values $\tilde{g}$ at given points $(t,y)$. To be more precise we
consider the following set of corrupting functions

$\displaystyle\mathcal{K}(\delta)=\\{\tilde{\delta}:[a,b]\times\mathbb{R}^{d}\to\mathbb{R}^{d}\
:\ \tilde{\delta}-\hbox{Borel measurable},\
\|\tilde{\delta}(t,y)\|\leq\delta(1+\|y\|)$
$\displaystyle\quad\quad\quad\quad\quad\hbox{for all}\
t\in[a,b],y\in\mathbb{R}^{d}\\}.$ (5.1)

Of course $\mathcal{K}(\delta_{1})\subset\mathcal{K}(\delta_{2})$ for all
$0\leq\delta_{1}\leq\delta_{2}\leq 1$. We also set

$V_{(\xi,g)}(\delta)=B_{\xi}(\delta)\times V_{g}(\delta),$ (5.2)

where

$V_{g}(\delta)=\bigcup_{\tilde{\delta}\in\mathcal{K}(\delta)}\\{g+\tilde{\delta}\\}.$
(5.3)

For $0\leq\delta_{1}\leq\delta_{2}\leq 1$ we have
$V_{(\xi,g)}(\delta_{1})\subset V_{(\xi,g)}(\delta_{2})$ and
$V_{(\xi,g)}(0)=\\{(\xi,g)\\}$.

Let $(\xi,g)\in F_{0}^{\alpha,\beta}$ and $(\tilde{\xi},\tilde{g})\in
V_{(\xi,g)}(\delta)$. A vector of noisy information about $(\xi,g)$ is of the
following form

$\mathcal{N}(\tilde{\xi},\tilde{g})=[\tilde{g}(t_{0},y_{0}),\tilde{g}(t_{1},y_{1}),\dots,\tilde{g}(t_{i-1},y_{i-1}),\tilde{\xi}]$

where $i\in\mathbb{N}$ is a total number of noisy evaluations of $g$. The
points evaluation points $(t_{k},y_{k})$ can be chosen in adaptive way. This
means that there exist functions $\psi_{j}$ such that

$(t_{0},y_{0})=\psi_{0}(\tilde{\xi}),$ (5.4)

and for $j=1,2,\ldots,i-1$

$(t_{j},y_{j})=\psi_{j}\Bigl{(}\tilde{g}(t_{0},y_{0}),\ldots,\tilde{g}(t_{j-1},y_{j-1}),\tilde{\xi}\Bigr{)}.$
(5.5)

An algorithm that approximates $z$ and uses the vector of (noisy) information
$\mathcal{N}$ is defined as

$\mathcal{A}(\tilde{\xi},\tilde{g})=\varphi(\mathcal{N}(\tilde{\xi},\tilde{g})),$
(5.6)

for some function

$\varphi:\mathbb{R}^{(i+1)d}\to B([a,b];\mathbb{R}^{d}),$ (5.7)

where $B([a,b];\mathbb{R}^{d})$ is the linear space of all functions
$l:[a,b]\to\mathbb{R}^{d}$ bounded in the supremum norm. For a given
$n\in\mathbb{N}$ we define $\Phi_{n}$ as a class of all algorithms of the form
(5.6) for which the total number of evaluations $i$ is at most $n$. For
example, $\tilde{l}_{n}\in\Phi_{n}$.

For a fixed $(\xi,g)\in F^{\alpha,\beta}_{0}$ the error of
$\mathcal{A}\in\Phi_{n}$ is defined as follows

$e(\mathcal{A},\xi,g,\delta)=\sup\limits_{(\tilde{\xi},\tilde{g})\in
V_{(\xi,g)}(\delta)}\sup\limits_{a\leq t\leq
b}\|z(\xi,g)(t)-\mathcal{A}(\tilde{\xi},\tilde{g})(t)\|,$ (5.8)

while the worst-case error of the algorithm $\mathcal{A}$ in a subclass
$\mathcal{G}\subset F^{\alpha,\beta}_{0}$

$e(\mathcal{A},\mathcal{G},\delta)=\sup\limits_{(\xi,g)\in\mathcal{G}}e(\mathcal{A},\xi,g,\delta).$
(5.9)

We investigate behavior of the $n$th minimal error defined as below

$e_{n}(\mathcal{G},\delta)=\inf_{\mathcal{A}\in\Phi_{n}}e(\mathcal{A},\mathcal{G},\delta).$
(5.10)

The following worst-case results in the class $F^{\alpha,\beta}_{R}$ (with
suitable $R$’s) for the Euler algorithm $\tilde{l}_{n}$, follow from the
proofs of Theorems 3.1, 4.1. Note that the error bounds in the Theorems 3.1,
4.1 are shown for a fixed $(\xi,g)$, while now we are considering the worst-
case error in the class $F^{\alpha,\beta}_{R}$. For the clarity we present
upper bounds on the error of the Euler algorithm in the case of exact and
inexact information separately.

######  Proposition 5.1.

(The case of exact information) Let
$R_{1}=e^{K(b-a)}(1+K)\left(e^{K(b-a)}+K(b-a)\right)$. Then there exists
$C\in(0,+\infty)$, depending only on the parameters of the class
$F^{\alpha,\beta}_{R_{1}}$, such for all $(\xi,g)\in
F_{R_{1}}^{\alpha,\beta}$, $n\in\mathbb{N}$ it holds

$\sup\limits_{a\leq t\leq b}\|z(t)-l_{n}(t)\|\leq C(h^{\alpha}+h^{\beta}).$
(5.11)

######  Proposition 5.2.

(The case of inexact information) Let
$R_{2}=\max\\{R_{1},K_{2}(1+K)\\}=e^{K(b-a)}(1+K)\max\\{e^{K(b-a)}+K(b-a),2\\}$.
Then there exists $C\in(0,+\infty)$, depending only on the parameters of the
class $F^{\alpha,\beta}_{R_{2}}$, such that for all $n\geq\lceil
2(b-a)\rceil+1$, $\delta\in[0,1]$, $(\xi,g)\in F^{\alpha,\beta}_{R_{2}}$,
$(\tilde{\xi},\tilde{g})\in V_{(\xi,g)}(\delta)$ the following holds

  * •

if $\beta=1$ then

$\sup\limits_{a\leq t\leq b}\|z(t)-\tilde{l}_{n}(t)\|\leq
C(h^{\alpha}+\delta).$ (5.12)

  * •

if $\beta\in(0,1)$ then

$\sup\limits_{a\leq t\leq b}\|z(t)-\tilde{l}_{n}(t)\|\leq
C\Bigl{(}h^{\alpha}+h^{\beta}+\mathbf{1}_{(0,1]}(\delta)\cdot(h^{1/2}+\delta)\Bigr{)}.$
(5.13)

In the special case when $\beta=1$ we get the following sharp bounds on the
$n$th minimal error in the class $F^{\alpha,1}_{R_{2}}$.

######  Theorem 5.3.

Let $\alpha\in(0,1]$, $\beta=1$ and $R_{2}$ defined as in Proposition 5.2.
Then it holds

$e_{n}(F^{\alpha,1}_{R_{2}},\delta)=\Theta(n^{-\alpha}+\delta),$ (5.14)

as $n\to+\infty$, $\delta\to 0^{+}$. The Euler algorithm $\tilde{l}_{n}$ is
the optimal one.

###### Proof.

The upper bound follows from (5.12) in Proposition 5.2 and the fact that
$\tilde{l}_{n}\in\Phi_{n}$.

Now we turn to the lower bound. We consider the following class

$\mathcal{G}^{\alpha,1}_{R_{2}}=\\{(0,g)\in F^{\alpha,1}_{R_{2}}\ |\
g(t,y)=g(t,0)\textrm{ for all }(t,y)\in[a,b]\times\mathbb{R}^{d}\\},$

which is a subclass of $F^{\alpha,1}_{R_{2}}$. Moreover, for
$(\xi,g)\in\mathcal{G}^{\alpha}_{R_{2}}$ we have that
$\displaystyle{z(\xi,g)(t)=\int\limits_{a}^{t}g(s,0)ds}$, $t\in[a,b]$. Hence,

$e_{n}(F^{\alpha,1}_{R_{2}},\delta)\geq
e_{n}(\mathcal{G}^{\alpha,1}_{R_{2}},0)\geq\inf_{\mathcal{A}\in\Phi_{n}}\sup\limits_{(0,g)\in\mathcal{G}^{\alpha,1}_{R_{2}}}\|z(0,g)(b)-\mathcal{A}(0,g)(b)\|=\Omega(n^{-\alpha}),$
(5.15)

which follows from the result on lower bound for a Riemann integration problem
of Hölder continuous functions under exact information, see Proposition 1.3.9,
page 34 in [29]. Moreover, using analogous argumentation as in the proof of
Theorem 3 in [2] we get, in the case of inexact information, that

$e_{n}(F^{\alpha,1}_{R_{2}},\delta)\geq(b-a)\delta.$ (5.16)

From (5.15) and (5.16) we get the lower bound in (5.14). ∎

######  Remark 5.4.

By considering the autonomous problem $z^{\prime}(t)=g(z(t))$, $z(a)=\xi$,
$t\in[a,b]$, with globally Lipschitz continuous function $g$ we know, from the
results of [14], that we can only get the lower bound $\Omega(h)$. Hence,
under our assumptions there is a gap between the upper bound $O(h^{\beta})$
and the lower bound $\Omega(h)$ \- we conjecture that the optimal error bound
is $\Theta(h^{\beta})$.

##  6\. Numerical results

We have conducted various numerical tests to further prove (5.11). The exact
solutions for the analyzed problems are unknown, so we proceeded as follows to
estimate the empirical rate of convergence. In parallel, we computed the
approximated solution and a numerical reference solution with a smaller step
size $h$ (1000 times smaller than in an approximated solution). Then we
computed an error as a difference between those two solutions.

The theoretical convergence rate is known. We presented the test results in
the log-log plots to test a numerical convergence rate. On the y axis there is
$-\log_{10}\textrm{(err)}$ and on the x axis the $\log_{10}n$, where err and
$n$ are given. Then, the linear regression slope corresponds to the empirical
convergence rate and is presented for the case with exact information.

We computed only one approximated solution to get the numerical convergence
rate for the case without noise. To gain the numerical convergence rate for a
noisy case, we conducted 200 Monte Carlo tries, and then we took the biggest
value of the error to present in the plot to mimic the worst-case setting
behavior.

###  6.1. Additive example

For all $(t,x)\in[a,b]\times\mathbb{R}$ let us consider the function

$g(t,x)=A|t|^{\alpha}x\sin(x^{2}+1)+B_{1}\textrm{sgn}(x)|x|^{\varrho_{1}}+B_{2}\textrm{sgn}(x)|x|^{\varrho_{2}},$
(6.1)

where $A\in\mathbb{R}$, $B_{1},B_{2}\in(-\infty,0]$,
$\alpha,\varrho_{1},\varrho_{2}\in(0,1]$. Note that the function (6.1)
satisfies assumptions (G1)-(G4), see Fact 7.4.

In the first experiment we take the following values of the parameters $a=0$,
$b=10$, $A=2$, $B_{1}=-2$, $B_{2}=-0.5$, $\alpha=1$ and $\rho_{1}=1$. We
manipulate with $\rho_{2}$ and noise level $\delta$.

In Figure 1 we have several reference solutions for
$\rho_{2}\in\\{1.0,0.75,0.5,0.25\\}$ and $\delta=0.0$. In Figure 2 we have
test results for $\delta=0$ (exact information) and
$\delta\in\\{0.1,0.01,0.001\\}$ (noisy information).

Figure 1. An approximated solution computed on the dense mesh, treated as the
reference solution of the model (6.1) with: $\rho_{2}=1$, $\rho_{2}=0.75$,
$\rho_{2}=0.5$, and $\rho_{2}=0.25$.

i $-\log_{10}\textrm{(err)}$ vs. $\log_{10}n$ for $\rho_{2}=1$.

ii $-\log_{10}\textrm{(err)}$ vs. $\log_{10}n$ for $\rho_{2}=0.75$.

iii $-\log_{10}\textrm{(err)}$ vs. $\log_{10}n$ for $\rho_{2}=0.5$.

iv $-\log_{10}\textrm{(err)}$ vs. $\log_{10}n$ for $\rho_{2}=0.25$.

Figure 2. $-\log_{10}\textrm{(err)}$ vs. $\log_{10}n$ for the model (6.1).

In the next experiment we take the following values of the parameters: $a=0$,
$b=10$, $A=2$, $B_{1}=-2$, $B_{2}=-0.5$, $\rho_{1}=1$, $\rho_{2}=0.75$ but
$\alpha=0.4$. In Figure 3 we compare reference solutions for
$\alpha\in\\{0.4,1.0\\}$ with $\delta=0$. In Figure 4 we have test results for
$\delta=0$ and $\delta\in\\{0.1,0.01,0.001\\}$.

Figure 3. An approximated solution computed on the dense mesh, treated as the
reference solution of the model (6.1) with: $\alpha=0.4$, $\alpha=1.0$.
Figure 4. $-\log_{10}\textrm{(err)}$ vs. $\log_{10}n$ for the model (6.1) with
$\alpha=0.4$.

###  6.2. Multiplicative example

Let us consider a family of the following functions

$g(t,x)=\gamma(t)h(x)f(x),\quad(t,x)\in[a,b]\times\mathbb{R}$ (6.2)

where $h(x)=-\textrm{sgn}(x)|x|^{\varrho}$, $\varrho\in(0,1]$, and for
$\gamma:[a,b]\to\mathbb{R}$, $f:\mathbb{R}\to\mathbb{R}$ the following
conditions are satisfied:

  * •

there exists $L_{\gamma}\geq 0$, $\alpha\in(0,1]$ such that for all
$t,s\in[a,b]$

$|\gamma(t)-\gamma(s)|\leq L_{\gamma}|t-s|^{\alpha},$ (6.3)

  * •

for all $t\in[a,b]$

$\gamma(t)\geq 0,$ (6.4)

  * •

there exists $D\geq 0$ such that for all $x\in\mathbb{R}$

$0\leq f(x)\leq D,$ (6.5)

  * •

there exists $L_{f}\geq 0$ such that for all $x,y\in\mathbb{R}$

$|f(x)-f(y)|\leq L_{f}|x-y|.$ (6.6)

It turns out that the function (6.2) satisfies the assumptions (G1)-(G4), see
Fact 7.5.

For the testing purpose we take the following particular function

$g(t,x)=-0.3\cdot|\sin(\pi\cdot
t+1)|^{2/3}\cdot\textrm{sgn}(x)|x|^{2/3}\cdot|\cos(x)|,$ (6.7)

with $\varrho=\frac{2}{3}$, $a=-3.5$, $b=3.5$, and with the initial-condition
$\xi=3$.

In Figure 5 we show a reference solution for exact information ($\delta=0$).
In Figure 6 we present the test results for $\delta=0$ and
$\delta\in\\{0.1,0.01,0.001\\}$.

Figure 5. A reference solution of the model (6.2).
Figure 6. $-\log_{10}\textrm{(err)}$ vs. $\log_{10}n$ for the model (6.2).

##  7\. Appendix

In this section we give some auxiliary facts about solution $z$ of the ODE
(1.3) and properties of the right-hand side function (6.1). We use the
following version of Peano’s theorem.

######  Theorem 7.1 (Peano’s theorem; [9], page 292., Theorem 70.4).

If the function $g$ satisfies the assumptions (G1) and (G2), then the initial
value problem (1.3) has at least one solution that belongs to
$C^{1}([a,b]\times\mathbb{R}^{d};\mathbb{R}^{d})$.

######  Lemma 7.2.

Let $\xi\in\mathbb{R}^{d}$ and let us assume that $g$ satisfies the
assumptions (G1), (G2) and (G3). Then there exists a unique solution
$z=z(\xi,g)\in C^{1}([a,b];\mathbb{R}^{d})$ of (1.3) and

$\sup_{a\leq t\leq b}\|z(t)\|\leq K_{0},$ (7.1)

where $K_{0}=K_{0}(a,b,K,\xi)=(\|\xi\|+K(b-a))e^{K(b-a)}\leq C_{0}(1+\|\xi\|)$
and $C_{0}=C_{0}(a,b,K)=\max\\{1,K(b-a)\\}e^{K(b-a)}$.

###### Proof.

By (G1), (G2) and from Theorem 7.1 the equation (1.3) has at least one
solution $z=z(\xi,g)\in C^{1}([a,b];\mathbb{R}^{d})$. Hence, by (G2) we get
for all $t\in[a,b]$

$\|z(t)\|\leq\|\xi\|+K(b-a)+K\int\limits_{a}^{t}\|z(s)\|ds,$ (7.2)

and from the Gronwall’s lemma we obtain that

$\sup_{a\leq t\leq b}\|z(t)\|\leq K_{0},$ (7.3)

where $K_{0}=K_{0}(a,b,K,\xi)=(\|\xi\|+K(b-a))e^{K(b-a)}$. Hence, $z(t)\in
B_{0}(K_{0})$ for all $t\in[a,b]$.

Let us consider two solutions $z=z(\xi,g)$ and $\tilde{z}=\tilde{z}(\xi,g)$ of
(1.3). We have that $z(t),\tilde{z}(t)\in B_{0}(K_{0})$ for all $t\in[a,b]$.
By (G3) applied to $U=B_{0}(K_{0})$ we have that there exists
$H_{U}\in\mathbb{R}$ such that for all $t\in[a,b]$

$\displaystyle\frac{d}{dt}\|z(t)-\tilde{z}(t)\|^{2}$ $\displaystyle=$
$\displaystyle\frac{d}{dt}\left\|\int_{a}^{t}(g(s,z(s))-g(s,\tilde{z}(s)))ds\right\|^{2}$
(7.4) $\displaystyle\leq$ $\displaystyle 2\langle
z(t)-\tilde{z}(t),g(t,z(t))-g(t,\tilde{z}(t))\rangle$ $\displaystyle\leq$
$\displaystyle 2H_{U}\|z(t)-\tilde{z}(t)\|^{2}\leq
2H_{U}^{+}\|z(t)-\tilde{z}(t)\|^{2}.$

Hence, the $C^{1}$-function $[a,b]\ni
t\to\varphi(t)=\|z(t)-\tilde{z}(t)\|^{2}\in[0,+\infty)$, where $\varphi(a)=0$,
satisfies the following integral inequality

$\varphi(t)\leq 2H_{U}^{+}\int_{a}^{t}\varphi(s)ds.$

From the Gronwall’s lemma we therefore get that $\varphi(t)=0$ for all
$t\in[a,b]$, which in turn implies that $z(t)=\tilde{z}(t)$ for all
$t\in[a,b]$. ∎

######  Remark 7.3.

An existence and a uniqueness of the solution of (1.3) assuming (G1)-(G3)
follow also form Theorem 3.27. page 182 in [30].

######  Fact 7.4.

The function (6.1) satisfies the assumptions (G1)-(G4) with $\alpha=\alpha$
and $\beta=\min\\{\varrho_{1},\varrho_{2}\\}$.

###### Proof.

We have that $g(t,x)=\sum\limits_{i=1}^{3}h_{i}(t,x)$ where
$h_{1}(t,x)=A|t|^{\alpha}x\sin(x^{2}+1)$ and
$h_{i}(t,x)=B_{i}\textrm{sgn}(x)|x|^{\varrho_{i}}$ for $i=2,3$,
$(t,x)\in[a,b]\times\mathbb{R}$. Of course $h_{i}\in
C([a,b]\times\mathbb{R};\mathbb{R})$ for $i=1,2,3$ and therefore $g$ satisfies
(G1). Moreover, for all $(t,x)\in[a,b]\times\mathbb{R}$

$|g(t,x)|\leq|A|\max\limits_{a\leq t\leq
b}|t|^{\alpha}|x|+(|B_{1}|+|B_{2}|)(1+|x|)\leq\Bigl{(}|A|\max\limits_{a\leq
t\leq b}|t|^{\alpha}+|B_{1}|+|B_{2}|\Bigr{)}(1+|x|),$ (7.5)

where $D_{1}=\max\limits_{a\leq t\leq b}|t|^{\alpha}$. Therefore, the function
$g$ satisfies (G2). Next, for all $t_{1},t_{2}\in[a,b]$,
$x_{1},x_{2}\in\mathbb{R}$ we have that

$|h_{1}(t_{1},x_{1})-h_{1}(t_{2},x_{2})|\leq|h_{1}(t_{1},x_{1})-h_{1}(t_{1},x_{2})|+|h_{1}(t_{1},x_{2})-h_{1}(t_{2},x_{2})|,$
(7.6)

where

$|h_{1}(t_{1},x_{1})-h_{1}(t_{1},x_{2})|\leq\frac{3}{2}|A|\cdot
D_{1}\cdot(1+x_{1}^{2}+x_{2}^{2})\cdot|x_{1}-x_{2}|,$ (7.7)

and, by Lemma A1.3 in [7],

$|h_{1}(t_{1},x_{2})-h_{1}(t_{2},x_{2})|\leq|A|\cdot|x_{2}|\cdot|t_{1}-t_{2}|^{\alpha}.$
(7.8)

This implies, in particular, that for all $t\in[a,b]$,
$x_{1},x_{2}\in\mathbb{R}$

$\displaystyle(x_{1}-x_{2})(h_{1}(t,x_{1})-h_{1}(t,x_{2}))\leq|(x_{1}-x_{2})(h_{1}(t,x_{1})-h_{1}(t,x_{2}))|$
$\displaystyle=|x_{1}-x_{2}|\cdot|h_{1}(t,x_{1})-h_{1}(t,x_{2})|\leq\frac{3}{2}|A|\cdot
D_{1}\cdot(1+x_{1}^{2}+x_{2}^{2})\cdot|x_{1}-x_{2}|^{2}.$ (7.9)

Since $B_{1},B_{2}\leq 0$, from the proof of Fact A.4 in [7] we have that for
all $t\in[a,b]$, $x_{1},x_{2}\in\mathbb{R}$, $i=2,3$

$(x_{1}-x_{2})(h_{i}(t,x_{1})-h_{i}(t,x_{2}))\leq 0,$ (7.10)

and for all $t_{1},t_{2}\in[a,b]$, $x_{1},x_{2}\in\mathbb{R}$, $i=2,3$

$|h_{i}(t_{1},x_{1})-h_{i}(t_{2},x_{2})|\leq
2|B_{i}|\cdot|x_{1}-x_{2}|^{\varrho_{i}}.$ (7.11)

Hence, by (7), (7.10) we get for all $t\in[a,b]$, $x_{1},x_{2}\in\mathbb{R}$

$\displaystyle(x_{1}-x_{2})(g(t,x_{1})-g(t,x_{2}))=(x_{1}-x_{2})\cdot\sum\limits_{i=1}^{3}(h_{i}(t,x_{1})-h_{i}(t,x_{2})$
$\displaystyle\leq(x_{1}-x_{2})(h_{1}(t,x_{1})-h_{1}(t,x_{2}))\leq\frac{3}{2}|A|\cdot
D_{1}\cdot(1+x_{1}^{2}+x_{2}^{2})\cdot|x_{1}-x_{2}|^{2}.$ (7.12)

Since for every $U\in\textrm{comp}(\mathbb{R})$ there exists a ball
$B_{0}(R)$, with the radius $R\in(0,+\infty)$, such that $U\subset B_{0}(R)$,
for all $t\in[a,b]$, $x_{1},x_{2}\in U$ we have

$(x_{1}-x_{2})(g(t,x_{1})-g(t,x_{2}))\leq\frac{3}{2}|A|\cdot
D_{1}\cdot(1+2R^{2})\cdot|x_{1}-x_{2}|^{2},$ (7.13)

and, by (7.7),(7.8), (7.11), for all $t_{1},t_{2}\in[a,b]$,
$x_{1},x_{2}\in\mathbb{R}$

$\displaystyle|g(t_{1},x_{1})-g(t_{2},x_{2})|\leq\sum\limits_{i=1}^{3}|h_{i}(t_{1},x_{1})-h_{i}(t_{2},x_{2})|$
$\displaystyle\leq|A|\cdot
R\cdot|t_{1}-t_{2}|^{\alpha}+\max\Bigl{\\{}\frac{3}{2}|A|D_{1}(1+2R^{2}),2|B_{1}|,2|B_{2}\Bigr{\\}}\cdot
3(1+2R)\cdot|x_{1}-x_{2}|^{\min\\{\varrho_{1},\varrho_{2}\\}}.$

Hence, $g$ satisfies also (G3), (G4), and this ends the proof. ∎

######  Fact 7.5.

The function (6.2) satisfies the assumptions (G1)-(G4) with $\alpha=\alpha$
and $\beta=\varrho$.

###### Proof.

Of course, $g\in C([a,b]\times\mathbb{R};\mathbb{R})$ and therefore $g$
satisfies (G1). Furthermore, for all $(t,x)\in[a,b]\times\mathbb{R}$

$\displaystyle|g(t,x)|$ $\displaystyle\leq
D\sup_{t\in[a,b]}|\gamma(t)|(1+|x|),$

hence $g$ satisfies (G2) with $K=D\sup\limits_{t\in[a,b]}|\gamma(t)|$.
Moreover, $h$ is monotonically decreasing, thus for all $t\in[a,b]$ and
$x,y\in\mathbb{R}$

$\displaystyle(x-y)\Bigl{(}g(t,x)-g(t,y)\Bigr{)}$
$\displaystyle=\gamma(t)(x-y)\Bigl{(}h(x)f(x)-h(y)f(y)\Bigr{)}$
$\displaystyle=\gamma(t)f(x)(x-y)(h(x)-h(y))+\gamma(t)(x-y)h(y)(f(x)-f(y))$
$\displaystyle\leq\gamma(t)(x-y)h(y)(f(x)-f(y))\leq|\gamma(t)(x-y)h(y)(f(x)-f(y))|$
$\displaystyle\leq
L_{f}\sup_{t\in[a,b]}|\gamma(t)|\cdot(1+|y|)\cdot|x-y|^{2}.$

For all $U\in\textrm{comp}(\mathbb{R})$ there exists a ball $B_{0}(R)$, $R>0$,
such that $U\subset B_{0}(R)$ and for all $t\in[a,b]$, $x,y\in U$ we have

$(x-y)\Bigl{(}g(t,x)-g(t,y)\Bigr{)}\leq
L_{f}\sup_{t\in[a,b]}|\gamma(t)|\cdot(1+R)\cdot|x-y|^{2}.$

Hence, $g$ satisfies the assumption (G3). Furthermore, for all
$t_{1},t_{2}\in[a,b]$ and $y_{1},y_{2}\in\mathbb{R}$ we have

$|g(t_{1},y_{1})-g(t_{2},y_{2})|\leq|g(t_{1},y_{1})-g(t_{2},y_{1})|+|g(t_{2},y_{1})-g(t_{2},y_{2})|,$

where

$|g(t_{1},y_{1})-g(t_{2},y_{1})|\leq
L_{\gamma}D\cdot(1+|y_{1}|)\cdot|t_{1}-t_{2}|^{\alpha},$ (7.14)

and

$\displaystyle|g(t_{2},y_{1})-g(t_{2},y_{2})|$ $\displaystyle\leq
L_{f}\sup_{t\in[a,b]}|\gamma(t)|\cdot(1+|y_{1}|)\cdot|y_{1}-y_{2}|+2D\sup_{t\in[a,b]}|\gamma(t)|\cdot|y_{1}-y_{2}|^{\varrho}$
$\displaystyle\leq\sup_{t\in[a,b]}|\gamma(t)|\cdot\max\\{L_{f},D\\}\cdot\Bigl{[}(1+|y_{1}|)\cdot(1+|y_{1}|+|y_{2}|)+2\Bigr{]}\cdot|y_{1}-y_{2}|^{\varrho}.$

Thus, since for all $U\in\textrm{comp}(\mathbb{R})$ there exists a ball
$B_{0}(R)$, $R>0$, such that $U\subset B_{0}(R)$, for all
$t_{1},t_{2}\in[a,b]$, $y_{1},y_{2}\in U$ we have

$\displaystyle|g(t_{2},y_{1})$ $\displaystyle-$ $\displaystyle
g(t_{2},y_{2})|\leq L_{\gamma}D(1+R)|t_{1}-t_{2}|^{\alpha}$
$\displaystyle+\sup_{t\in[a,b]}|\gamma(t)|\cdot\max\\{L_{f},D\\}\Bigl{[}(1+R)(1+2R)+2\Bigr{]}|y_{1}-y_{2}|^{\varrho}.$

This implies that $g$ satisfies (G4) with $\alpha=\alpha$ and $\beta=\varrho$,
and the proof is completed. ∎

## References

  * [1] W. F. Ames, Numerical methods for partial differential equations. 3rd ed., Academic Press, San Diego, 1992.
  * [2] T. Bochacik, P. Przybyłowicz, On the randomized Euler schemes for ODEs under inexact information, to appear in Numer. Algor., https://doi.org/10.1007/s11075-022-01299-7
  * [3] T. Bochacik, M. Goćwin, P. M. Morkisz, P. Przybyłowicz, Randomized Runge–Kutta method—Stability and convergence under inexact information, J. Complexity 65 (2021), 101554.
  * [4] J. C. Butcher, Numerical methods for ordinary differential equations. 3rd ed., Wiley, Chichester, 2016.
  * [5] E. W. Cheney, D. R. Kincaid, Numerical Mathematics and Computing, 7th ed., Brooks/Cole Cengage Learning, Boston, 2012.
  * [6] N. Czyżewska, J. Kusiak, P. Morkisz, P. Oprocha, M. Pietrzyk, P. Przybyłowicz, Ł. Rauch, D. Szeliga, On mathematical aspects of evolution of dislocation density in metallic materials, IEEE Access 10 (2022), 86793–86812.
  * [7] N. Czyżewska, P. M. Morkisz, P. Przybyłowicz, Approximation of solutions of DDEs under nonstandard assumptions via Euler scheme, to appear in Numer. Algor., DOI:10.1007/s11075-022-01324-9.
  * [8] R. D. Driver, Ordinary and delay differential equations. Springer-Verlag, New York, 1977.
  * [9] L. Górniewicz, R. S. Ingarden, Mathematical Analysis for Physicists (in Polish). Wydawnictwo Naukowe UMK, 2012.
  * [10] Griffiths D. F., Higham D. J., Numerical Methods for Ordinary Differential Equations. Initial Value Problems. Springer-Verlag, London, 2010.
  * [11] E. Hairer and G. Wanner, Solving Ordinary Differential Equations II. Stiff and Differential-Algebraic Problems, 2nd rev. ed., Springer-Verlag, Heidelberg, New York, 2010.
  * [12] S. Heinrich, Complexity of initial value problems in Banach spaces, Zh. Mat. Fiz. Anal. Geom. 9 (2013), 73–101.
  * [13] Z. Jackiewicz, General linear methods for ordinary differential equations. Wiley, Hoboken, New Jersey, 2009.
  * [14] B. Kacewicz, Optimal solution of ordinary differential equations, J. Complex. 3 (1987), 451–465.
  * [15] B. Kacewicz, M. Milanese, A. Vicino, Conditionally optimal algorithms and estimation of reduced order models, J. Complex. 4 (1988), 73–85.
  * [16] B. Kacewicz, L. Plaskota, On the minimal cost of approximating linear problems based on information with deterministic noise, Numer. Funct. Anal. and Optimiz. 11 (1990), 511-528.
  * [17] B. Kacewicz, P. Przybyłowicz, On the optimal robust solution of IVPs with noisy information, Numer. Algor. 71 (2016), 505–518.
  * [18] B. Kacewicz, P. Przybyłowicz, Efficient finite-dimensional solution of initial value problems in infinite-dimensional Banach spaces, J. Math. Anal. Appl. 471 (2019), 322-341.
  * [19] A. Kałuża, P. M. Morkisz, P. Przybyłowicz, Optimal approximation of stochastic integrals in analytic noise model, Appl. Math. and Comput. 356 (2019), 74–91.
  * [20] P. E. Kloeden, E. Platen, Numerical Solution of Stochastic Differential Equations. Springer-Verlag, Berlin, 3rd ed., 1999.
  * [21] S. Larsson, V. Thomée, Partial differential equations with numerical methods. Springer-Verlag, Berlin, 2003.
  * [22] S. Mazumder, Numerical methods for partial differential equations. Finite difference and finite volume methods. Academmic Press, London, 2016.
  * [23] M. Milanese, A. Vicino, Optimal estimation theory for dynamic systems with set membership uncertainty: an overview, Automatica 27 (1991), 997–1009.
  * [24] G. N. Milstein, M. V. Tretyakov, Stochastic Numerics for Mathematical Physics. Springer-Verlag, Berlin, 2004.
  * [25] P. M. Morkisz, L. Plaskota, Approximation of piecewise Hölder functions from inexact information, J. Complex. 32 (2016), 122–136.
  * [26] P. M. Morkisz, L. Plaskota, Complexity of approximating Hölder classes from information with varying Gaussian noise, to appear in J. Complex. 60 (2020), ID 101497.
  * [27] P. M. Morkisz, P. Przybyłowicz, Optimal pointwise approximation of SDE’s from inexact information, J. Comput. Appl. Math. 324 (2017), 85–100.
  * [28] P. M. Morkisz, P. Przybyłowicz, Randomized derivative-free Milstein algorithm for efficient approximation of solutions of SDEs under noisy information, J. Comput. Appl. Math. 383 (2021), ID 113112.
  * [29] E. Novak, Deterministic and Stochastic Error Bounds in Numerical Analysis, Lecture Notes in Mathematics, vol. 1349, New York, Springer–Verlag, 1988.
  * [30] E. Pardoux, A. Rascanu, Stochastic Differential Equations, Backward SDEs, Partial Differential Equations, Stochastic Modelling and Applied Probability, Springer International Publishing Switzerland, 2014.
  * [31] L. Plaskota, Noisy Information and Computational Complexity, Cambridge Univ. Press, Cambridge, 1996.
  * [32] L. Plaskota, Noisy information: optimality, complexity, tractability, in Monte Carlo and quasi-Monte Carlo Methods 2012, J. Dick, F.Y. Kuo, G.W. Peters, I.H. Sloan (Eds.), Springer 2013, 173–209.
  * [33] A.A. Samarskii, The Theory of Difference Schemes. Marcel Dekker, New York, 2001.
  * [34] J.F. Traub, G.W. Wasilkowski, H. Woźniakowski, Information-Based Complexity, Academic Press, New York, 1988.
  * [35] A.G. Werschulz, The complexity of definite elliptic problems with noisy data. J. Complex. 12 (1996), 440-473.
  * [36] A.G. Werschulz, The complexity of indefinite elliptic problems with noisy data. J. Complex. 13 (1997), 457-479.


# Low-frequency signature of magnetization nutation in nanomagnets

M. P. Adams1<EMAIL_ADDRESS>R. Bastardis2<EMAIL_ADDRESS>A. Michels1<EMAIL_ADDRESS>H. Kachkachi2<EMAIL_ADDRESS>1Department of Physics and Materials Science, University of Luxembourg, 162A
avenue de la Faiencerie, L-1511 Luxembourg, Grand Duchy of Luxembourg
2Université de Perpignan via Domitia, Lab. PROMES CNRS UPR8521, Rambla de la
Thermodynamique, Tecnosud, 66100 Perpignan, France

###### Abstract

In this work, we show that surface anisotropy in nanomagnets induces a
nutational motion of their magnetization at various frequencies, the lowest of
which can be described by the macrospin model whose dynamics is governed by an
effective energy potential. We derive analytical expressions for the
precession and nutation frequencies and amplitudes as functions of the size of
the nanomagnet and its atomistic parameters, such as the exchange coupling and
the onsite anisotropy. Our analytical model predicts a reduction of the
precession frequency with increased surface anisotropy. We also simulate the
dynamics of the corresponding atomistic many-spin system and compare the
results with the effective model. We thereby show that the first nutation mode
induced by the finite size and surface anisotropy occurs at a frequency that
is four times larger than the precession frequency, thus lending itself to a
relatively easy detection by standard experiments of magnetic resonance.

##  I introduction

Nutation is the well known motion of a gyroscope described by classical
mechanics; it develops whenever an external force tilts the rotation axis away
from the direction of the gravity field, since then the rotation axis no
longer coincides with the direction of the angular momentum.

In the case of a magnetic material, magnetization nutation occurs whenever a
time-dependent magnetic field (rf or microwave field) is present in the
system. Accordingly, it has been demonstrated that there appears the
fundamental effect of transient nutations in NMR [1], EPR [2, 3, 4] and
optical resonance [5]. Spin nutation was first predicted in Josephson
junctions [6, 7, 8, 9, 10] and was later developed using various approaches
based on relativistic quantum mechanics, first principles [11, 12], and
electronic structure theory [13, 14, 15, 16, 17]. A more recent macroscopic
approach [18, 19, 20, 21] deals with the magnetization nutation by accounting
for magnetic inertia through the addition of a second-order time derivative of
the magnetization in the Landau-Lifshitz-Gilbert equation [22, 23]. On the
experimental side, evidence of the effects of nutation has been reported for
thin films [24, 25, 26, 27] in the THz regime. Indeed, magnetization nutation
meets with a strong and expanding interest within the magnetism community
owing to the real possibility of addressing and controlling the ultra-fast
magnetization dynamics, in view of potential practical applications in high-
speed data processing and storage technologies based on magnetic materials.
Considering the strong potential of achieving high storage densities, this
issue acquires still more enthusiasm for nanoscaled magnetic systems.

The dynamic behavior of nanomagnets requires a good understanding of the
intricate interplay between different factors influencing the spin dynamics,
such as finite-size, boundary, and surface effects. In particular, surface
anisotropy emerges as a crucial parameter that can significantly influence the
dynamics of nanomagnets as it affects the potential energy and thereby the
relaxation rates [28, 29]. Recently, it was shown [30, 31] that spin
misalignment, induced by surface effects, triggers nutational motion of the
magnetization of a nanomagnet, with frequencies ranging from tens of GHz to
THz, depending on various physical parameters, the predominant among them
being the spin-spin exchange coupling and the on-site surface anisotropy. For
example, in Refs. [32, 31] it was shown that surface anisotropy in nanocubes
triggers absorption modes of higher frequency than the ordinary uniform
precessional motion with frequencies $f_{\mathrm{p}}$$\sim$$10\,\mathrm{GHz}$
(conventional ferromagnetic resonance). More precisely, the authors of Ref.
[31] observed a resonance at a frequency $f_{\mathrm{c}}$ that is 4 times
$f_{\mathrm{p}}$ and another frequency $f_{\mathrm{n}}$ in the THz regime.
While $f_{\mathrm{n}}$ is related to atomic spin fluctuations and their
exchange coupling, the frequency $f_{\mathrm{c}}$ emerges from the cubic
anisotropy induced by spin inhomogeneities at the surface of the nanomagnet
[33, 34, 35]. As such, $f_{\mathrm{c}}$ corresponds to the nutational mode
with the lowest frequency and should be more easily detected in resonance
experiments on, _e.g._ arrays of nanomagnets with acute surface effects, such
as platelets, pillars, spheres, or cubes.

In the present work, we focus on the frequency $f_{\mathrm{c}}$ and derive its
analytical expression using the effective macrospin model whose potential
energy is a polynomial in the components of the net magnetic moment, with
coefficients that are functions of the spin-spin exchange coupling and on-site
anisotropy constants [33, 34, 35, 36]. The results obtained using the
effective model are compared with those of numerical simulations for the
corresponding many-spin nanomagnet. The effective model, thus validated,
allows us to derive analytical expressions for the precession and nutation
frequencies and the nutation amplitude and to study them in terms of the
nanomagnet size and surface anisotropy.

Plan of the paper: We first define the many-spin and macrospin approaches and
establish the connection between them. Then, we discuss the numerical results
from both models for the time trajectories of the net magnetic moment. In
particular, we discuss the results regarding the magnetization nutation and
its dependence on the particle size and surface anisotropy. This study allows
us to validate the effective model which is then used to derive the analytical
expressions for the precession and nutation frequencies and amplitudes. From
these expressions and further simulations, we draw the main conclusions of the
present work.

##  II Model for the magnetization dynamics

###  II.1 Many-spin approach (MSP)

In the many-spin problem (MSP) or many-spin approach, or still atomistic
approach, the nanomagnet is regarded as a crystallite of a given shape and
size, containing $\mathcal{N}$ atomic spins
$\boldsymbol{S}_{i}=\mu_{\mathrm{a}}\boldsymbol{s}_{i}$ located at the nodes
of a given crystalline structure; $\mu_{\mathrm{a}}$ denotes the atomic
magnetic moment and $\boldsymbol{s}_{i}$ is a unit vector in the direction of
$\boldsymbol{S}_{i}$. The magnetic state of such a system is studied with the
help of the atomistic Hamiltonian [37, 38, 39, 40, 41, 42, 43, 44, 45, 46]:

$\displaystyle\mathcal{H}$
$\displaystyle=\mathcal{H_{\mathrm{exch}}}+\mathcal{H}_{\mathrm{\mathrm{Z}}}+\mathcal{H}_{\mathrm{A}}$
$\displaystyle=-\frac{1}{2}J\sum_{\left\langle
i,j\right\rangle}\boldsymbol{s}_{i}\cdot\boldsymbol{s}_{j}-\mu_{\mathrm{a}}\mathbf{B}_{0}\cdot\sum_{i=1}^{\mathcal{N}}\boldsymbol{s}_{i}+\sum_{i=1}^{\mathcal{N}}\mathcal{H}_{\mathrm{A},i},$
(1)

where $\mathcal{H_{\mathrm{exch}}}$ is the nearest-neighbor (n.n.) exchange
energy with $J>0$, $\mathcal{H}_{\mathrm{\mathrm{Z}}}$ denotes the Zeeman
energy with $\mathbf{B}_{0}=\mu_{0}\mathbf{H}_{0}$ being the homogeneous
externally applied magnetic field, and $\mathcal{H}_{\mathrm{A}}$ the magnetic
anisotropy energy. For the core spins, we assume the anisotropy to be of
uniaxial symmetry with constant $K_{\mathrm{c}}$, while for surface spins we
adopt the model proposed by Néel [47] with constant $K_{\mathrm{s}}$.
$\mathcal{H}_{\mathrm{A},i}$ is then given by

$\mathcal{H}_{\mathrm{A},i}=\begin{cases}-K_{\mathrm{c}}\left(\boldsymbol{s}_{i}\cdot\mathbf{e}_{z}\right)^{2},&i\in{\mathrm{core}}\\\
+\frac{1}{2}K_{\mathrm{s}}{\sum_{j\in{\mathrm{n.n.}}}}\left(\boldsymbol{s}_{i}\cdot\mathbf{u}_{ij}\right)^{2},&i\in\mathrm{surface},\end{cases}$
(2)

where the unit vectors
$\mathbf{u}_{ij}=(\mathbf{r}_{i}-\mathbf{r}_{j})/\|\mathbf{r}_{i}-\mathbf{r}_{j}\|$
connect the nearest-neighbor spins “$i$” and “$j$”. All physical parameters
are measured in units of energy per atom.

The (undamped) dynamics of the many-spin system is governed by the following
Larmor equation:

$\frac{d\boldsymbol{s}_{i}}{d\tau}=\boldsymbol{s}_{i}\times\bm{b}_{{\rm
eff},i}$ (3)

with the (normalized) local effective field $\bm{b}_{{\rm eff},i}$ acting on
$\bm{s}_{i}$ being defined by $\bm{b}_{{\rm
eff},i}=-\delta\mathcal{H}/\delta\bm{s}_{i}$; $\tau$ is the reduced time
defined by:

$\tau\equiv\frac{t}{\tau_{{\rm s}}},$ (4)

where $\tau_{{\rm s}}=\mu_{\mathrm{a}}/\left(\gamma J\right)$ is a
characteristic time of the system’s dynamics. By way of example, for cobalt
$J=8{\rm\ meV}$, leading to $\tau_{{\rm s}}=70\,{\rm fs}$. Henceforth, we will
only use the dimensionless time $\tau$. Accordingly, the frequency
$\omega=2\pi f=2\pi/T$ is measured in units of $\tau_{{\rm s}}^{-1}$ and, as
such, in the sequel $\nu\equiv\tau_{{\rm s}}\omega$. Note that in these units,
$b_{{\rm eff}}$ is equal to the effective field (in Tesla) multiplied by
$\mu_{\mathrm{a}}/J$ and is thus dimensionless.

To study the dynamics of this MSP, for arbitrary values of all physical
parameters, one resorts to numerical methods for solving Eq. (3) using, for
instance, the iterative routine based on the $4^{\mathrm{th}}$-order Runge-
Kutta scheme in combination with the projection step
$\boldsymbol{s}_{i}^{\nu+1}=\boldsymbol{s}_{i}^{\mathrm{RK4}}/\|\boldsymbol{s}_{i}^{\mathrm{RK4}}\|$
to preserve the constraint $\|\boldsymbol{s}_{i}\|=1$. The net magnetic moment
$\boldsymbol{m}$ (unit vector) is then computed using

$\boldsymbol{m}=\frac{\sum_{i}\boldsymbol{s}_{i}}{\|\sum_{i}\boldsymbol{s}_{i}\|}.$
(5)

For the initial state of the system, we choose a coherent spin configuration
along the orientation $\left(\theta_{0},\phi_{0}\right)$, _i.e._
$\bm{s}_{i}=\left(\sin\theta_{0}\cos\phi_{0},\sin\theta_{0}\sin\phi_{0},\cos\theta_{0}\right)$,
for all $i=1,\ldots,\mathcal{N}$.

Nutational motion was investigated in Ref. [31] using this MSP model and the
results were compared with those of the macrospin model studied with the help
of the augmented Landau-Lifshitz-Gilbert equation. In the present work, our
aim is to derive analytical expressions for the precession and nutation
frequencies in terms of the atomistic parameters
$J,K_{\mathrm{c}},K_{\mathrm{s}}$ and the size of the nanomagnet. For this
purpose, we resort to the effective one-spin problem derived in Refs. [33, 34,
35, 36] and compare the corresponding results with those obtained for the MSP
in Ref. [31] in regard with the nutational motion.

###  II.2 Effective macrospin approach (EOSP)

The effective macrospin approach, or effective one-spin problem (EOSP), is
obtained from the many-spin problem under certain conditions regarding the
surface anisotropy, which should not be too strong with respect to the
exchange coupling, and the particle size, which should not be too small (see
Refs. [33, 34, 36, 35] for details). The EOSP consists of the net magnetic
moment Eq. (5) evolving in an effective energy potential given by

$\mathcal{H}_{\mathrm{eff}}\simeq-
k_{2}m_{z}^{2}+k_{4}\sum_{\alpha=x,y,z}m_{\alpha}^{4}.$ (6)

The values and signs of the (effective) coefficients $k_{2}$ and $k_{4}$ are
functions of the atomistic parameters $J,K_{\mathrm{c}},K_{\mathrm{s}}$, in
addition to the size and shape of the nanomagnet and its underlying crystal
lattice [35]. In the following, we will use the dimensionless constants
defined by $k_{\mathrm{c}}=K_{\mathrm{c}}/J$ and
$k_{\mathrm{s}}=K_{\mathrm{s}}/J$, so that $k_{2}$ and $k_{4}$ are
dimensionless.

For a nanomagnet cut out of a simple cubic lattice, we have [33, 48]

$k_{2}=k_{\mathrm{c}}\frac{N_{\mathrm{c}}}{\mathcal{N}},\quad
k_{4}\simeq\begin{cases}\kappa\frac{k_{\mathrm{s}}^{2}}{z},&\mathrm{sphere},\\\
\\\
\left(1-0.7/\mathcal{N}^{1/3}\right)^{4}\frac{k_{\mathrm{s}}^{2}}{z},&\mathrm{cube},\end{cases}$
(7)

where $z=6$ is the coordination number, $\kappa$ represents a (dimensionless)
surface integral, and $N_{\mathrm{c}}$ is the number of atoms in the core of
the nanomagnet (with full coordination), while $\mathcal{N}$ is the total
number of atoms (including both the core and surface). Note that $\kappa$ was
derived in Ref. 33 in the absence of core anisotropy. As shown in Ref. [49],
in the presence of the latter, the spin misalignment caused by surface
anisotropy does not propagate to the center of the nanomagnet; it is “fended
off” by the uniaxial anisotropy in the core which tends to align all spins
together. The result of this competition is that $k_{4}$ (or $\kappa$) is
multiplied by the factor $N_{\mathrm{c}}/\mathcal{N}$. Likewise, for both the
cube and sphere, $k_{2}$ may be approximated [34] by the first expression in
Eq. (7).

From the effective Hamiltonian in Eq. (6), we derive the following effective
field that drives the dynamics of the magnetic moment $\boldsymbol{m}$

$\boldsymbol{b}_{\mathrm{eff}}=2k_{2}m_{z}\boldsymbol{e}_{z}-4k_{4}\left(m_{x}^{3}\boldsymbol{e}_{x}+m_{y}^{3}\boldsymbol{e}_{y}+m_{z}^{3}\boldsymbol{e}_{z}\right)$
(8)

through the Larmor equation of motion for $\boldsymbol{m}$ as follows:

$\frac{d\boldsymbol{m}}{d\tau}=\boldsymbol{m}\times\boldsymbol{b}_{\mathrm{eff}}.$
(9)

We would like to emphasize that the EOSP described here is not some variant of
a single-domain macrospin model. Indeed, it is a macrospin model, but with
coefficients that are functions of the atomistic physical parameters. In other
terms, it has the advantage to be a "simple" macroscopic model which inherits
(to some extent) the nanomagnet atomistic features, especially spin
inhomogeneities.

##  III magnetization nutation

###  III.1 MSP versus EOSP

Figure 1: Components of the net magnetic moment for both MSP (continued lines) and EOSP (dashed lines and symbols) for the spherical nanomagnet with a sc lattice and size $\mathcal{N}=1088$. The constants $k_{\mathrm{c}}$ and $k_{\mathrm{s}}$ are indicated in the legend and the effective coefficients are given in Table 1. On the right, $T_{\mathrm{p}},T_{\mathrm{c}},T_{\mathrm{n}}$ denote the periods corresponding to the respective frequencies $f_{\mathrm{p}},f_{\mathrm{c}},f_{\mathrm{n}}$ mentioned in the text. $\mathcal{N}=1088$ ($N_{\mathrm{c}}=706)$ | $k_{2}$ | $k_{4}$ |  $\omega_{\mathrm{c}}$(GHz)
---|---|---|---
$k_{\mathrm{c}}=0.001,k_{\mathrm{s}}=0.01$ | $6.5\times 10^{-4}$ | $5.8\times 10^{-6}$ | 51.5
$k_{\mathrm{c}}=0.01,k_{\mathrm{s}}=0.1$ | $6.5\times 10^{-3}$ | $5.81\times 10^{-4}$ | 515
Table 1: Effective parameters used for the simulations on the nanospheres
[from Eq. (7)].

Figure 2: Components of the net magnetic moment for both MSP (continued lines) and EOSP (dashed lines and symbols) for a cube-shaped nanomagnet with a sc lattice and sizes $\mathcal{N}=15^{3}$ (left) and $\mathcal{N}=9^{3}$ (right). The constants $k_{\mathrm{c}}$ and $k_{\mathrm{s}}$ are indicated in the legend and the effective coefficients are given in Table 2.  $\mathcal{N}$ | $N_{\mathrm{c}}$ | $k_{2}$ | $k_{4}$ |  $\omega_{\mathrm{c}}$(GHz)
---|---|---|---|---
$729=9^{3}$ | 343 | $4.7\times 10^{-4}$ | $12\times 10^{-6}$ | 355
$3375=15^{3}$ | 2197 | $6.5\times 10^{-4}$ | $13.8\times 10^{-6}$ | 52
Table 2: Effective parameters used for the simulations on the nanocubes [from
Eq. (7)].

Figure 3: Components of the net magnetic moment for both MSP (continued lines) and EOSP (dashed lines and symbols) for the truncated octahedron with fcc lattice and size $\mathcal{N}=2075$. The constants $k_{\mathrm{c}}$ and $k_{\mathrm{s}}$ are indicated in the legend and the effective coefficients are given in Table 3. $\mathcal{N}=2075$ | $k_{2}$ | $k_{4}$ |  $\omega_{\mathrm{c}}$(GHz)
---|---|---|---
$k_{\mathrm{c}}=0.001,k_{\mathrm{s}}=0.01$ | $6.7\times 10^{-4}$ | $-4.5\times 10^{-6}$ | 54.6
$k_{\mathrm{c}}=0.01,k_{\mathrm{s}}=0.1$ | $6.7\times 10^{-3}$ | $-2.9\times 10^{-4}$ | 546
Table 3: Effective parameters used for the simulations on the truncated
octahedrons obtained by a fitting procedure. Note that $k_{4}$ is negative for
the fcc lattice [35].

In this section, we compare the dynamics of the net magnetic moment of
different many-spin nanomagnets with that of the corresponding effective
models for spheres [50], cubes [51, 52, 53, 54, 55, 56, 32] and truncated
octahedrons [57, 58, 35]. Before discussing the results in more detail, let us
summarize the procedure followed for the comparison between the MSP and EOSP.
For a given nanomagnet with a given size, shape, underlying lattice, and
energy parameters ($J,K_{\mathrm{c}},K_{\mathrm{s}}$), we solve the undamped
Larmor equation (3) for all atomic spins $\bm{s}_{i}$ starting from the
initial state $\left(\theta_{0},\phi_{0}\right)$ taken here to be
$\left(\pi/4,\pi/4\right)$. Then, we plot the components $m_{\alpha}$ (with
$\alpha=x,y,z$) of the net magnetic moment (of the MSP) as functions of
time111Note that in all plots we denote by $s_{\alpha}$, $\alpha=x,y,z$, the
Cartesian components of the net magnetic moment of the MSP and by
$m_{\alpha}$, $\alpha=x,y,z$, those of the magnetic moment of the EOSP.. Next,
for a cube or a sphere with a simple cubic lattice, we use the expressions for
$k_{2}$ and $k_{4}$ in Eqs. (7) to solve the Larmor Eq. (9) for the EOSP and
obtain the components $m_{\alpha}$ as functions of time. For other shapes
and/or lattice structures, such as the truncated octahedron or cubes with
_e.g._ body-centered-cubic lattice [see Fig. 3], $k_{4}$ is obtained by
fitting the MSP data. Note that in the case of an underlying bcc or fcc
lattice, $k_{4}$ becomes negative [35].

The results of such a comparison confirm the following facts : the net
magnetic moment of the many-spin nanomagnet exhibits three oscillation modes
with different frequencies [see Fig. 1 right]: the precession frequency
$f_{\mathrm{p}}$ for the $x,y$ components and the nutation frequencies
$f_{\mathrm{c}}$ and $f_{\mathrm{n}}$ for the $z$ component. The data in Figs.
1, 2, 3 are plotted as continued lines, in black and red for $s_{x,y}$ and in
blue for $s_{z}$. The corresponding EOSP magnetic moment exhibits only two
oscillation frequencies: the precession frequency $f_{\mathrm{p}}$ for
$m_{x,y}$ (in red and black symbols) and the (smaller) nutation frequency
$f_{\mathrm{c}}$ for $m_{z}$ (in magenta). This result obtains in all cases of
size, shape and other physical parameters, as confirmed in Figs. 1, 2 and 3.
We see that as far as $f_{\mathrm{p}}$ and $f_{\mathrm{c}}$ are concerned, and
this corresponds to the validity domain of the EOSP (i.e._,_ relatively weak
surface anisotropy), the dynamics of EOSP perfectly matches that of the MSP.
However, the latter exhibits (in blue) an extra wiggling motion in the
component $s_{z}$ (with frequency $f_{\mathrm{n}}$) on top of a signal with
frequency $f_{\mathrm{c}}$. As mentioned in the introduction and discussed in
Ref. [31], the appearance of the frequency $f_{\mathrm{c}}$ is due to the
quartic term in the Hamiltonian (6), which is a consequence of spin disorder
induced by surface anisotropy in the EOSP regime; the coefficient of this term
($k_{4}$) is a function of $K_{\mathrm{s}}$, as can be seen from Eq. (7). On
the other hand, $f_{\mathrm{n}}$ is due to the fluctuations of the individual
atomic spins with frequencies on the order of the exchange coupling $J$.

Further analysis of Figs. 1, 2, 3 reveals further useful information. First,
the two frequencies $f_{\mathrm{p}}$ and $f_{\mathrm{c}}$ change with size,
shape, and energy parameters. As the size decreases, the ratio
$N_{\mathrm{c}}/\mathcal{N}$ decreases, and the surface contribution to the
overall energy increases [49]. Second, as the surface anisotropy constant
$K_{\mathrm{s}}$ increases, the nutation frequency $f_{\mathrm{n}}$ is more
clearly identified with an increasing amplitude. This effect is more clearly
seen in Figs. 1 and 3 (plots on the right). However, the value of the
frequency itself changes with the exchange coupling. Likewise, when going from
the cube to the more rounded geometry of the sphere, through the truncated
octahedron, we see that the nutational motion is progressively enhanced and
its amplitude increases.

###  III.2 Nutation frequency and amplitude: analytical approach

From the previous section we conclude that the effective model (EOSP)
perfectly describes the dynamics of the many-spin nanomagnet as far as the
validity conditions are met (_i.e._ small surface anisotropy). On the other
hand, we have seen that, in addition to the nutation frequency
$f_{\mathrm{c}}$, the MSP exhibits another nutation mode with a much higher
frequency $f_{\mathrm{n}}$; the amplitude of this mode increases with the
enhancement of surface effects, either through an increase of the anisotropy
constant or that of the number of surface spins.

Now, with the macroscopic EOSP model at hand, we can derive analytical
expressions for the frequency and amplitude of the precession and nutation
modes and study their behavior in terms of surface anisotropy and size.
Instead of the Cartesian components $m_{x},m_{y},m_{z}$ of $\boldsymbol{m}$,
it is more convenient to parametrize the vector $\boldsymbol{m}$ using the
spherical coordinates $\boldsymbol{q}=\left(\theta,\phi\right)$, with $\theta$
being the polar angle and $\phi$ the azimuth angle, _i.e._ ,
$\boldsymbol{m}\left(\boldsymbol{q}\right)=\left(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta\right)$.
Then, in this coordinate system, the Larmor Eq. (9) becomes

$\frac{d\boldsymbol{q}}{d\tau}=\boldsymbol{b}_{\mathrm{eff}}^{\prime}\left(\boldsymbol{q}\right),$
(10)

where we have defined
$\boldsymbol{b}_{\mathrm{eff}}^{\prime}=\left(h_{\theta},h_{\phi}\right)$ as
the new effective driving field. More explicitly, the components of the
spherical driving field $\boldsymbol{b}_{\mathrm{eff}}^{\prime}$ are given by

$\displaystyle b_{\mathrm{eff},\theta}^{\prime}\left(\boldsymbol{q}\right)$
$\displaystyle=k_{4}\sin\left(4\phi\right)\sin^{3}\theta,$ (11) $\displaystyle
b_{\mathrm{eff},\phi}^{\prime}\left(\boldsymbol{q}\right)$
$\displaystyle=2k_{2}\cos\theta-
k_{4}\cos\theta[7\cos^{2}\theta-\sin^{2}\theta\cos\left(4\phi\right)-3].$

In Eq. (11), the appearance of trigonometric functions with the argument
$4\phi$ already hints to the fact that we should observe a mode with a
frequency four times that of the “fundamental” mode, namely the precession
mode.

To analytically solve the Larmor equation (10), we make the following Ansatz
for the angular motion:

$\boldsymbol{q}_{p}\left(\tau\right)=\left(\begin{array}[]{c}\theta_{p}\left(\tau\right)\\\
\phi_{p}\left(\tau\right)\end{array}\right)=\left(\begin{array}[]{c}\theta_{0}\\\
\phi_{0}+\nu_{\mathrm{p}}\tau\end{array}\right),$ (12)

where $\omega_{\mathrm{p}}$ is the (circular) precession frequency. Note that
this Ansatz is only valid for $k_{4}=0$. For finite but small $k_{4}$, we
assume that the dynamics is approximately governed by the equation

$\frac{d\boldsymbol{q}}{d\tau}\simeq\boldsymbol{b}_{\mathrm{eff}}^{\prime}\left[\boldsymbol{q}_{p}\left(\tau\right)\right].$
(13)

More explicitly, this yields:

$\displaystyle\frac{d\theta}{d\tau}$ $\displaystyle\simeq-
k_{4}\sin^{3}\theta_{0}\sin\left(4\nu_{\mathrm{p}}\tau\right),$ (14)
$\displaystyle\frac{d\phi}{d\tau}$
$\displaystyle\simeq-\cos\theta_{0}[2k_{2}+3k_{4}-7k_{4}\cos^{2}\theta_{0}]-k_{4}\cos\theta_{0}\sin^{2}\theta_{0}\cos\left(4\nu_{\mathrm{p}}\tau\right).$

From these approximate solutions, we infer the following equations:

$\displaystyle\nu_{\mathrm{p}}$
$\displaystyle=\cos\theta_{0}\left[2k_{2}+\left(3-7\cos^{2}\theta_{0}\right)k_{4}\right],$
(15) $\displaystyle\nu_{\mathrm{c}}$ $\displaystyle=4\nu_{\mathrm{p}},$
$\displaystyle a_{\phi}$
$\displaystyle=\frac{1}{4}\,\frac{k_{4}\sin^{2}\theta_{0}}{2k_{2}+\left(3-7\cos^{2}\theta_{0}\right)k_{4}},\,a_{\theta}=a_{\phi}\tan\theta_{0},$

where $\nu_{\mathrm{p}}$ and $\nu_{\mathrm{c}}$ are the (dimensionless)
precession and nutation frequencies, and $a_{\theta},a_{\phi}$ the amplitudes.
Note that in situations where $k_{4}$ becomes negative, the latter should be
replaced by $\left|k_{4}\right|$ in the expression of the amplitude
$a_{\phi}$.

Integrating Eqs. (14) yields the time trajectories of the net magnetic moment:

$\displaystyle\boldsymbol{q}\left(\tau\right)$
$\displaystyle\simeq\left(\begin{array}[]{c}\theta_{0}\\\
\phi_{0}+\nu_{\mathrm{p}}\tau\end{array}\right)-\left(\begin{array}[]{c}a_{\theta}\left[1-\cos\left(\nu_{\mathrm{c}}\tau\right)\right]\\\
a_{\phi}\sin\left(\nu_{\mathrm{c}}\tau\right)\end{array}\right).$ (20)

The set of Eqs. (15) and (20) constitutes the original result of this work.
Before discussing the conclusions we can draw from them, we summarize how they
can be used in practice. Given a (many-spin) nanomagnet of a certain size and
shape and a chosen set of energy parameters $J,K_{\mathrm{c}},K_{\mathrm{s}}$,
so that the spin configuration is not too much disordered, we can characterize
the (modified) precession and the first nutation modes by their frequencies
and amplitudes given by the corresponding effective macrospin model. If the
nanomagnet is a nanocube or a nanosphere with the underlying simple-cubic
lattice, the effective coefficients $k_{2}$ and $k_{4}$ are given by Eqs. (7).
In the general case, one may fit the MSP curves $m_{\alpha}\left(\tau\right)$
to Eq. (20) and infer the quantities in Eqs. (15). In fact, once the constant
$k_{4}$ is obtained, all the other quantities
($\nu_{\mathrm{p}},\nu_{\mathrm{c}},a_{\theta},a_{\phi}$) can easily be
derived. By way of illustration, we have applied the latter procedure to the
much studied iron nanocubes [59, 60, 51, 52, 61] with an underlying bcc
lattice. The results are shown in Fig. 4, where $k_{4},\omega_{\mathrm{c}}$
and $a_{\theta}$ are plotted against the (linear) cube size $N$; the fit for
$k_{4}$ yields $k_{4}\simeq-12.5k_{\mathrm{s}}^{2}/8\left(1+N\right)^{5/3}$.

From Eq. (15) and the numerical results in Fig. 4, we make the following
observations:

  * •

The frequencies $\nu_{\mathrm{p}}$ and $\nu_{\mathrm{c}}$ increase when the
size increases but decrease with surface anisotropy, _i.e._ when
$K_{\mathrm{s}}$ ($k_{4}\propto k_{\mathrm{s}}^{2}$) increases, for a small
initial polar angle $\theta_{0}$ (precession angle). For a relatively large
angle $\theta_{0}$ ($\gtrsim 49^{\circ}$), the frequency becomes an increasing
function of $K_{\mathrm{s}}$ (plot not shown).

  * •

The nutation frequency $\nu_{\mathrm{c}}$ is four times the precession
frequency $\nu_{\mathrm{p}}$, independently of the nanomagnet’s size.

Regarding the amplitude, from Eqs. (15) we see that it increases with the
surface anisotropy constant $K_{\mathrm{s}}$, as already concluded from the
results of the MSP simulations in Figs. 1, 2, 3. On the other hand, its
behavior as a function of the size depends on the shape. In Fig. 4, we see
that it is a decreasing function of $N$.

For a cube with simple-cubic lattice,
$N_{\mathrm{c}}/\mathcal{N}=\left(1-2/N\right)^{3}$ and thereby
$a_{\phi},a_{\theta}\sim 1/N$, which means that when the size increases, the
nutation mode dies out. This has a very important implication. In general, any
bulk or thin magnetic film presents an anisotropy energy that is an expansion
in the components of its net magnetic moment, similar to the Hamiltonian in
Eq. (6). However, the dynamics of a bulk magnetic material with such a
Hamiltonian does not necessarily exhibit magnetic nutation. Here, we see that
nanoscaled magnetic systems, such as nanocubes, do develop such a nutational
motion at a frequency that should be observable in resonance experiments on
arrays of such objects. For a sphere, using $k_{2}$ and $k_{4}$ from Eq. (7),
Eq. (15) implies that the amplitudes $a_{\phi},a_{\theta}$ do not depend on
the size, as has also been confirmed by the MSP numerical simulations. For an
arbitrary shape and/or lattice structure, it is not easy to derive simple
formulae for the coefficients $k_{2},k_{4}$. Nevertheless, as was discussed
above, they may be obtained by fitting the dynamics of the net magnetic moment
of the many-spin nanomagnet to Eqs. (20). This is what has been done for the
bcc cubic nanomagnets, studied by experimental groups [51, 52, 53, 54, 55, 56,
32], and the results are shown in Fig. 4.

Figure 4:  EOSP parameter $k_{4}$, nutation frequency $\omega_{\mathrm{c}}$,
and amplitude $a_{\theta}$ as a function of the linear size $N$ of a bcc iron
cube of size $N\times N\times N$ ($k_{\mathrm{c}}=0.001$ and
$k_{\mathrm{s}}=0.01$). The parameters were obtained from a fit of the many-
spin nanomagnet results to Eqs. (15) and (20). The full curves are guides to
the eye.

##  IV Conclusion

We have shown that the magnetization of a nanoscaled magnetic material may
exhibit nutational motion with frequencies in the GHz to THz range, the lowest
of which is four times the precession frequency. We have provided analytical
expressions for precession and nutation frequencies and the amplitude of the
nutational oscillations as functions of the size of the nanomagnet and
atomistic energy parameters, such as the spin-spin exchange coupling and the
on-site surface and core magneto-crystalline anisotropies.

This has been possible owing to the correspondence we have firmly established
between the many-spin approach and the effective macrospin approach, a
correspondence that has been validated through numerical simulations of the
dynamics of various nanomagnets (cubes, spheres, and truncated octahedrons)
with different energy parameters and lattice structures. The results of these
simulations have confirmed that the effective model recovers very well the
dynamics of the many-spin nanomagnet in the low-frequency regime, when surface
anisotropy is not too strong. In this regime, we observe the precession
frequency $f_{\mathrm{p}}$ and the lowest nutation frequency $f_{\mathrm{c}}$
that is four times $f_{\mathrm{p}}$ and whose existence is related to the
effective cubic anisotropy induced by spin inhomogeneities. The many-spin
nanomagnet also exhibits a much higher frequency $f_{\mathrm{n}}$ that can be
related to the spin-spin exchange coupling. Even the precession frequency
$f_{\mathrm{p}}$ turns out to be affected by surface anisotropy, among other
parameters. Within the effective macrospin model, we have demonstrated that
the amplitude of the nutational motion is enhanced by surface anisotropy, but
its behavior with the size depends on the shape of the nanomagnet.

All in all, in this work, we want to emphasize the fact that surface-induced
magnetization nutation in nanoscaled magnets, such as nanocubes and
nanospheres, already occurs at the lowest frequency that is four times the
precession frequency, which itself is altered by spin misalignments. This low-
frequency nutational oscillation lends itself to a detection and
characterization using standard experimental techniques of magnetic resonance,
such as the network analyzer with varying field and frequency, at low
temperature. An array of well separated platelets or nanocubes should provide
the adequate conditions for such measurements. Inelastic neutron spectroscopy
[62] might also be able to resolve the predicted nutational spin dynamics in
nanomagnets.

###### Acknowledgements.

M.P. Adams and A. Michels thank the National Research Fund of Luxembourg for
financial support (AFR Grant No. 15639149). H. Kachkachi thanks A. Michels and
the Department of Physics (Univ. Luxembourg) for the hospitality extended to
him during his stay.

## References

  * Torrey [1949] H. C. Torrey, Phys. Rev. 76, 1059 (1949).
  * Verma and Fessenden [1973] N. C. Verma and R. W. Fessenden, J. Chem. Phys. 58, 2501 (1973).
  * P. W. Atkins, A. J. Dobbs, and K. A. McLauchlan [1974] P. W. Atkins, A. J. Dobbs, and K. A. McLauchlan, Chem. Phys. Lett. 25, 105 (1974).
  * Fedoruk [2002] G. G. Fedoruk, J. Appl. Spectroscopy 69, 161 (2002).
  * Hocker and Tang [1968] G. B. Hocker and C. L. Tang, Phys. Rev. Lett. 21, 592 (1968).
  * Zhu and Fransson [2006] J.-X. Zhu and J. Fransson, J. Phys.: Condens. Matter 18, 9929 (2006).
  * Fransson and Zhu [2008] J. Fransson and J.-X. Zhu, New Journal of Physics 10, 013017 (2008).
  * Fransson [2008] J. Fransson, Nanotechnology 19, 285714 (2008).
  * Nussinov et al. [2005] Z. Nussinov, A. Shnirman, D. P. Arovas, A. V. Balatsky, and J. X. Zhu, Phys. Rev. B 71, 214520 (2005).
  * Zhu et al. [2004] J.-X. Zhu, Z. Nussinov, A. Shnirman, and A. V. Balatsky, Phys. Rev. Lett. 92 (2004).
  * Mondal et al. [2017] R. Mondal, M. Berritta, A. K. Nandy, and P. M. Oppeneer, Phys. Rev. B 96, 024425 (2017).
  * Mondal and Oppeneer [2020] R. Mondal and P. M. Oppeneer, J. Phys.: Condens. MatterB 32, 455802 (2020).
  * Bhattacharjee et al. [2012] S. Bhattacharjee, L. Nordström, and J. Fransson, Phys. Rev. Lett. 108, 057204 (2012).
  * F"ahnle et al. [2011] M. Fähnle, D. Steiauf, and C. Illg, Phys. Rev. B 84, 172403 (2011), URL https://link.aps.org/doi/10.1103/PhysRevB.84.172403.
  * Kikuchi and Tatara [2015] T. Kikuchi and G. Tatara, Phys. Rev. B 92, 184410 (2015).
  * Thonig D., Eriksson O., Pereiro M. [2017] Thonig D., Eriksson O., Pereiro M., Scientific Reports 7, 931 (2017).
  * R. Cheng, X. Wu, and D. Xiao [2017] R. Cheng, X. Wu, and D. Xiao, Phys. Rev. B 96, 054409 (2017).
  * M.-C. Ciornei, J. M. Rubi, and J. E. Wegrowe [2011] M.-C. Ciornei, J. M. Rubi, and J. E. Wegrowe, Phys. Rev. B 83, 020410 (2011).
  * E. Olive, Y. Lansac, and J. E. Wegrowe [2012] E. Olive, Y. Lansac, and J. E. Wegrowe, Appl. Phys. Lett. 100, 192407 (2012).
  * Olive et al. [2015] E. Olive, Y. Lansac, M. Meyer, M. Hayoun, and J.-E. Wegrowe, Journal of Applied Physics 117, 213904 (2015), eprint https://doi.org/10.1063/1.4921908, URL https://doi.org/10.1063/1.4921908.
  * E. Olive and J. E. Wegrowe [2016] E. Olive and J. E. Wegrowe, J. Phys.: Condens. Mat. 28, 106001 (2016).
  * L. D. Landau and E. M. Lifshitz [1935] L. D. Landau and E. M. Lifshitz, Phys. Z. Sowjetunion 8, 153 (1935).
  * Gilbert [1956] T. L. Gilbert, Ph.D. thesis, Illinois Institute of Technology, Chicago (1956).
  * Li et al. [2015] Y. Li, A.-L. Barra, S. Auffret, U. Ebels, and W. E. Bailey, Phys. Rev. B 92, 140413 (2015).
  * K. Neeraj et al. [2021] K. Neeraj et al., Nature Physics 17, 245 (2021).
  * U.Unikandanunni et al. [2022] U.Unikandanunni et al., Phys.Rev.Lett. 129, 237201 (2022).
  * A. De et al. [2024] A. De, J. Schlegel, A. Lentfert, L. Scheuer, B. Stadtmüller1, Ph. Pirro, G. von Freymann, U. Nowak, and M Aeschlimann, arXiv:2405.01334 , (2024).
  * P.-M. Déjardin, H. Kachkachi, Yu. Kalmykov [2008] P.-M. Déjardin, H. Kachkachi, Yu. Kalmykov, J. Phys. D 41, 134004 (2008).
  * Vernay et al. [2014] F. Vernay, Z. Sabsabi, and H. Kachkachi, Phys. Rev. B 90, 094416 (2014).
  * B"ottcher and Henk [2012] D. Böttcher and J. Henk, Phys. Rev. B 86, 020404 (2012).
  * Bastardis, R. and Vernay, F. and Kachkachi, H. [2018] Bastardis, R. and Vernay, F. and Kachkachi, H., Phys. Rev. B 98, 165444 (2018).
  * Bastardis et al. [2017] R. Bastardis, F. Vernay, D. A. Garanin, and H. Kachkachi, J. Phys. C: Condens. Matter 29, 025801 (2017).
  * Garanin and Kachkachi [2003] D. A. Garanin and H. Kachkachi, Phys. Rev. Lett. 90, 065504 (2003), URL http://link.aps.org/doi/10.1103/PhysRevLett.90.065504.
  * H. Kachkachi and E. Bonet, Surface-induced cubic anisotropy in nanomagnets [2006] H. Kachkachi and E. Bonet, Surface-induced cubic anisotropy in nanomagnets, Phys. Rev. B 73, 224402 (2006).
  * R. Yanes, O. Fesenko-Chubykalo, H. Kachkachi, D.A. Garanin, R. Evans, R. W. Chantrell [2007] R. Yanes, O. Fesenko-Chubykalo, H. Kachkachi, D.A. Garanin, R. Evans, R. W. Chantrell, Phys. Rev. B 76, 064416 (2007).
  * H. Kachkachi, Effects of spin non-collinearities in magnetic nanoparticles [2007] H. Kachkachi, Effects of spin non-collinearities in magnetic nanoparticles, J. Magn. Magn. Mater. 316, 248 (2007).
  * D.A. Dimitrov and Wysin, Effects of surface anisotropy on hysteresis in fine magnetic particles [1994] D.A. Dimitrov and Wysin, Effects of surface anisotropy on hysteresis in fine magnetic particles, Phys. Rev. B 50, 3077 (1994).
  * R.H. Kodama and A.E. Berkovitz, Atomic-scale magnetic modeling of oxide nanoparticles [1999] R.H. Kodama and A.E. Berkovitz, Atomic-scale magnetic modeling of oxide nanoparticles, Phys. Rev. B 59, 6321 (1999).
  * Kachkachi and Garanin [2001a] H. Kachkachi and D. A. Garanin, Physica A Statistical Mechanics and its Applications 300, 487 (2001a), eprint cond-mat/0001278.
  * Kachkachi and Garanin [2001b] H. Kachkachi and D. A. Garanin, European Physical Journal B 22, 291 (2001b), eprint cond-mat/0012017.
  * O. Iglesias and A. Labarta, Finite-size and surface effects in maghemite nanoparticles: Monte Carlo simulations [2001] O. Iglesias and A. Labarta, Finite-size and surface effects in maghemite nanoparticles: Monte Carlo simulations, Phys. Rev. B 63, 184416 (2001).
  * H. Kachkachi and M. Dimian, Hysteretic properties of a magnetic particle with strong surface anisotropy [2002] H. Kachkachi and M. Dimian, Hysteretic properties of a magnetic particle with strong surface anisotropy, Phys. Rev. B 66, 174419 (2002).
  * H. Kachkachi and D.A. Garanin [2005] H. Kachkachi and D.A. Garanin, in _Surface effects in magnetic nanoparticles_ , edited by D. Fiorani (Springer, Berlin, 2005), p. 75.
  * N. Kazantseva, D. Hinzke, U. Nowak, R. W. Chantrell, U. Atxitia, and O. Chubykalo-Fesenko, Towards multiscale modeling of magnetic materials: Simulations of FePt [2008] N. Kazantseva, D. Hinzke, U. Nowak, R. W. Chantrell, U. Atxitia, and O. Chubykalo-Fesenko, Towards multiscale modeling of magnetic materials: Simulations of FePt, Phys. Rev. B 77, 184428 (2008).
  * M. P. Adams, A. Michels and H. Kachkachi [2022] M. P. Adams, A. Michels and H. Kachkachi, J. Appl. Cryst. 55, 1475 (2022).
  * Adams et al. [2024] M. P. Adams, E. P. Sinaga, H. Kachkachi, and A. Michels, Phys. Rev. B 109, 024429 (2024), URL https://link.aps.org/doi/10.1103/PhysRevB.109.024429.
  * L. Néel, Anisotropie magnétique superficielle et surstructures d’orientation [1954] L. Néel, Anisotropie magnétique superficielle et surstructures d’orientation, J. Phys. Radium 15, 225 (1954).
  * D. A. Garanin [2023] D. A. Garanin, Phys. Rev. B 98, 054427 (2023).
  * M. P. Adams, A. Michels and H. Kachkachi [2023] M. P. Adams, A. Michels and H. Kachkachi, Phys. Scr. 98, 105512 (2023).
  * Salzemann et al. [2023] C. Salzemann, V. Russier, M. Pancaldi, P. Vavassori, A. Berger, and I. Lisiecki, Colloids and Surfaces A: Physicochemical and Engineering Aspects 678, 132473 (2023).
  * Snoeck et al. [2008] E. Snoeck, C. Gatel, L. M. Lacroix, T. Blon, S. Lachaize, J. Carrey, M. Respaud, and B. Chaudret, Nano Letters 8, 4293 (2008).
  * Trunova et al. [2008] A. V. Trunova, R. Meckenstock, I. Barsukov, C. Hassel, O. Margeat, M. Spasova, J. Lindner, and M. Farle, Journal of Applied Physics 104, 093904 (2008).
  * Jiang et al. [2010] F. Jiang, C. Wang, Y. Fu, and R. Liu, Journal of Alloys and Compounds 503, L31 (2010), ISSN 0925-8388.
  * Mehdaoui et al. [2010a] B. Mehdaoui, A. Meffre, L.-M. Lacroix, J. Carrey, S. Lachaize, M. Gougeon, M. Respaud, and B. Chaudret, J. Magn. Magn. Mater. 322, L49 (2010a).
  * Kronast et al. [2011] F. Kronast, N. Friedenberger, K. Ollefs, S. Gliga, L. Tati-Bismaths, R. Thies, A. Ney, R. Weber, C. Hassel, F. M. Römer, et al., Nano Letters 11, 1710 (2011).
  * O’Kelly et al. [2012] C. O’Kelly, S. J. Jung, A. P. Bell, and J. J. Boland, Nanotechnology 23, 435604 (2012).
  * M. Jamet, W. Wernsdorfer, C. Thirion, D. Mailly, V. Dupuis, P. Mélinon, and A. Pérez [2001] M. Jamet, W. Wernsdorfer, C. Thirion, D. Mailly, V. Dupuis, P. Mélinon, and A. Pérez, Phys. Rev. Lett. 86, 4676 (2001).
  * M. Jamet, W. Wernsdorfer, C. Thirion, V. Dupuis, P. Mélinon, A. Pérez, and D. Mailly [2004] M. Jamet, W. Wernsdorfer, C. Thirion, V. Dupuis, P. Mélinon, A. Pérez, and D. Mailly, Phys. Rev. B 69, 24401 (2004).
  * P. Tartaj, P. Morales, S. Veintemillas-Verdaguer, T. Gonzalez-Carreno, and C. Serna [2003] P. Tartaj, P. Morales, S. Veintemillas-Verdaguer,T. Gonzalez-Carreno, and C. Serna, J. Phys. D: Appl. Phys. 36, R182 (2003).
  * Tran [2006] M. Tran, Master’s thesis, Institut National des Sciences Appliquees de Toulouse, Toulouse (2006).
  * Mehdaoui et al. [2010b] B. Mehdaoui, A. Meffre, L.-M. Lacroix, J. Carrey, S. Lachaize, M. Gougeon, M. Respaud, and B. Chaudret, Journal of Magnetism and Magnetic Materials 322, L49 (2010b), eprint 0907.4063.
  * Ch. Franz. S. S"aubert, A. Wendl, F. X. Haslbeck, O. Soltwedel, J. K. Jochum, L. Spitz, J. Kindervater, A. Bauer, P. B"oni, Ch. Pfleiderer [2019] Ch. Franz. S. Säubert, A. Wendl, F. X. Haslbeck, O. Soltwedel, J. K. Jochum, L. Spitz, J. Kindervater, A. Bauer, P. Böni, Ch. Pfleiderer, J. Phys. Soc. Japan 88, 081002 (2019).