# Large nearest neighbour balls
in hyperbolic stochastic geometry

Moritz Otto<EMAIL_ADDRESS>Department of Mathematics, Aarhus University,
Denmark Christoph Thäle<EMAIL_ADDRESS>Faculty of Mathematics, Ruhr
University Bochum, Germany

###### Abstract

Consider a stationary Poisson process in a $d$-dimensional hyperbolic space.
For $R>0$ define the point process $\xi_{R}^{(k)}$ of exceedance heights over
a suitable threshold of the hyperbolic volumes of $k$th nearest neighbour
balls centred around the points of the Poisson process within a hyperbolic
ball of radius $R$ centred at a fixed point. The point process $\xi_{R}^{(k)}$
is compared to an inhomogeneous Poisson process on the real line with
intensity function $e^{-u}$ and point process convergence in the Kantorovich-
Rubinstein distance is shown. From this, a quantitative limit theorem for the
hyperbolic maximum $k$th nearest neighbour ball with a limiting Gumbel
distribution is derived.

Keywords. Geometric extreme value theory, hyperbolic stochastic geometry,
nearest neighbour balls, Poisson process approximation.
MSC. 52A55, 60D05, 60G55.

##  1 Introduction

The study of extreme values, or more generally processes of exceedance heights
and associated order statistics, is a classical topic in probability theory. A
systematic study of extreme values for random geometric systems is more recent
and we refer, for example, to [6, 7, 10, 12] for particular results on the
Poisson-Voronoi, -Delaunay or -line tessellation, to [22, 29] for
distinguished results on random interpoint distances, and to [4, 13, 26, 27,
30, 31] for general approaches leading to various other stochastic-geometric
applications. A systematic study of random processes of exceedance heights in
stochastic geometry is the content of [4, 13, 24].

In this paper we are interested in quantitative limit theorems for so-called
large $k$th nearest neighbour balls, another classical stochastic geometry
model whose investigation goes back to [18] and which has recently been
studied in [4, 11, 17]. In particular, nearest neighbour balls (that is, $k$th
nearest neighbour balls with $k=1$) can be regarded as spatial analogues of
the concept of spacings in dimension one. In its simplest form the model can
be described as follows: Take a sequence $(X_{i})_{i\in{\mathbb{N}}}$ of
independent random points which are uniformly distributed on the
$d$-dimensional unit cube $[0,1]^{d}\subset{\mathbb{R}}^{d}$. For $n\geq 1$,
$k\in\\{1,\ldots,n\\}$ and $i\in\\{1,\ldots,n\\}$ let $r(i,n)$ be the distance
of $X_{i}$ to its $k$th nearest neighbour among the points
$X_{1},\ldots,X_{n}$, where the distance is understood in the Euclidean sense.
Then for $t\in{\mathbb{R}}$ define the random variable

$C_{n}:=\sum_{i=1}^{n}\mathds{1}\Big{\\{}{\mathcal{H}}_{e}^{d}(B_{e}(X_{i},r(i,n)))\geq{t+\log
n\over n}\Big{\\}},$

where ${\mathcal{H}}_{e}^{d}$ stands for the $d$-dimensional Hausdorff measure
and $B_{e}(X_{i},r(i,n))$ for the $d$-dimensional ball of radius $r(i,n)$
centred at $X_{i}$ with respect to the Euclidean structure on
${\mathbb{R}}^{d}$. In other words, $C_{n}$ counts the number of exceedances
of volumes of $k$th nearest neighbour balls that are larger than the threshold
$(t+\log n)/n$. It follows from the results in [4, 11, 17] that, after
suitable normalization, $C_{n}$ converges in distribution, as $n\to\infty$, to
a Poisson random variable with mean $e^{-t}$. Moreover, the rate of
convergence, measured in the total variation distance, is of order $(\log\log
n)/\log n$. This result immediately leads to a limit theorem for the maximum
volume $M_{n}:=\max\\{{\mathcal{H}}_{e}^{d}(B_{e}(X_{i},r(i,n))):1\leq i\leq
n\\}$ of the $k$th nearest neighbour balls, which says that $nM_{n}$
converges, after suitable centring, to a Gumbel distribution, as $n\to\infty$.
A similar result holds if the sample of $n$ independent random points is
replaced by a homogeneous Poisson process in $[0,1]^{d}$ with intensity $n$.

While the results and references just mentioned deal with $k$th nearest
neighbour balls in a $d$-dimensional Euclidean space, we follow another line
of current research in stochastic geometry and introduce and study a similar
model in a $d$-dimensional hyperbolic space ${\mathbb{H}}^{d}$ of constant
negative curvature $-1$. Random geometric systems in such a non-Euclidean set-
up have so far been studied in the context of random polytopes [2, 3, 16],
random graphs [5, 15, 14, 25] and tessellations [16, 21]. However, the study
of extreme values in hyperbolic stochastic geometry has so far left no trace
in the existing literature. The present paper can be understood as a first
attempt in this direction. Moreover, since the results for $k$th nearest
neighbour balls in Euclidean space have found applications in goodness-of-fit
testing for point processes [18, 19], our contributions can also be of
interest for similar studies in hyperbolic space.

In principle it is possible to rephrase the Euclidean model of large $k$th
nearest neighbour balls in a hyperbolic space, where the cube (which does not
exist in hyperbolic geometry) is replaced by a hyperbolic ball of radius one,
say. However, in this case, one can localize the problem and work with
approximations in the corresponding tangent spaces. Within these tangent
spaces the model is Euclidean again and we get back a result similar to that
in [4, 11, 17]. For this reason, we modify the set-up as follows: We start
with a stationary Poisson process in ${\mathbb{H}}^{d}$ and look at large
$k$th nearest neighbour balls associated with points in a family of hyperbolic
balls of radius $R\to\infty$. Up to a rescaling, in a Euclidean space this
set-up is the same as fixing the radius of the ball and increasing the
intensity of the Poisson process (or equivalently the number of points).
However, this is no more the case in a hyperbolic space. Even more, since the
problem in this form cannot locally be approximated by Euclidean models in
tangent spaces, we will arrive at results which are of a different nature and
‘feel’ the negative curvature of the underlying space.

In the next section we formally describe the framework we work with and
present our results.

##  2 Set-up and results

Fix a dimension parameter $d\in{\mathbb{N}}$ and consider a $d$-dimensional
hyperbolic space ${\mathbb{H}}^{d}$ together with the intrinsic (Riemannian)
metric $d_{h}$ and the corresponding $d$-dimensional Hausdorff measure
${\mathcal{H}}^{d}$. Although all our results are independent of a concrete
model for ${\mathbb{H}}^{d}$ (as can be seen from the fact that non of our
arguments or computations rely on a specific model), for concreteness one may
consider the Beltrami-Klein model in which ${\mathbb{H}}^{d}$ is identified
with the open Euclidean unit ball $\mathbb{B}^{d}$ and the Riemannian metric
is given by

${\rm d}s^{2}={{\rm d}x_{1}^{2}+\ldots+{\rm d}x_{d}^{2}\over
1-x_{1}^{2}-\ldots-x_{d}^{2}}+{(x_{1}{\rm d}x_{1}+\ldots+x_{d}{\rm
d}x_{d})^{2}\over(1-x_{1}^{2}-\ldots-x_{d}^{2})^{2}},$

see [8, 28] for details and further models for ${\mathbb{H}}^{d}$. While in
this model hyperbolic hyperplanes are non-empty intersections of Euclidean
hyperplanes with $\mathbb{B}^{d}$, hyperbolic balls are represented by
Euclidean ellipsoids, see Figure 1. For $z\in\mathbb{H}^{d}$ and $r>0$ let
$B(z,r):=\\{x\in\mathbb{H}^{d}:\,d_{h}(x,z)\leq r\\}$ denote the closed
hyperbolic ball of radius $r$ centred at $z$. We abbreviate $B_{r}:=B(p,r)$,
where $p\in\mathbb{H}^{d}$ is some arbitrary fixed point, referred to as the
origin of ${\mathbb{H}}^{d}$.

Figure 1: Construction of nearest neighbour balls ($k=1$) in the Beltrami-
Klein model for the hyperbolic plane. The ball $B_{R}$ is shown in black with
a white centre, the black points are points from $\eta$ together with their
nearest neighbour balls. The area of the two blue balls exceed the value
$v_{1}(R)$.

Let $\eta$ be a Poisson process in $\mathbb{H}^{d}$, $d\geq 2$, with intensity
measure $\mathcal{H}^{d}$. We note that $\eta$ is stationary in the sense that
its distribution is invariant under all isometries of the hyperbolic space.
For $R>0$ and $k\in{\mathbb{N}}$ let

$\displaystyle
v_{k}(R):=R\,(d-1)+{(k-1)}\log(R\,(d-1))-\log\left(\frac{{(k-1)}!2^{d-1}(d-1)}{\omega_{d}}\right),$
(2.1)

where $\omega_{d}={2\pi^{d/2}\over\Gamma({d\over 2})}$ is the surface area of
the $(d-1)$-dimensional Euclidean unit sphere. For $x\in\mathbb{H}^{d}$,
$k\in{\mathbb{N}}$ and a general, locally finite and simple counting measure
$\mu$ on $\mathbb{H}^{d}$ let $r(x,k,\mu)$ denote the hyperbolic distance to
the $k$th nearest neighbour of $x$ in $\mu$. In the focus of our results is
the point process

$\displaystyle\xi_{R}^{(k)}:=\sum_{x\in\eta\cap
B_{R}}\delta_{\mathcal{H}^{d}(B_{r(x,k,\eta-\delta_{x})})-v_{k}(R)},\qquad
R>0,$ (2.2)

on the real line ${\mathbb{R}}$. This process describes the exceedance heights
over the threshold $v_{k}(R)$ of the volumes of $k$th nearest neighbour balls
centred around the points of $\eta$ within a ball $B_{R}$ of radius $R$, see
Figure 1. We give a brief heuristic argument which explains at least for $k=1$
that $v_{k}(R)$ given in (2.1) is the correct threshold to expect Poisson
approximation for $\xi_{R}^{(k)}$. Given some point of the Poisson process
$\eta$, the probability that the hyperbolic ball of volume $v_{1}(R)+c$ around
this point does not contain any further points of $\eta$ is
$\exp(-v_{1}(R)+c)$. Assuming that all such balls around all points of
$\eta\cap B_{R}$ behave asymptotically independently (which is, of course, not
true and requires justification), the probability $\mathbb{P}[\mathcal{E}(R)]$
of the event $\mathcal{E}(R)$ that in $\eta\cap B_{R}$ there is no 1-nearest
neighbour ball of radius larger than $v_{1}(R)+c$ should be approximately
equal to $\exp(-\mathcal{H}^{d}(B_{R})\exp(-v_{1}(R)-c))$, where
$\mathcal{H}^{d}(B_{R})$ is for large $R$ approximately the number of points
of $\eta\cap B_{R}$. Using that $\mathcal{H}^{d}(B_{R})\exp(-v_{1}(R))\to 1$
as $R\to\infty$ (see (3.3)) we obtain from the classical Poisson limit theorem
that $\mathbb{P}[\mathcal{E}(R)]$ converges to the Gumbel limit
$\exp(-e^{-c})$ as $R\to\infty$.

Our main results quantify on the positive real half-axis $\mathbb{R}_{+}$ the
approximation of $\xi_{R}^{(k)}$ by a suitable Poisson process, where we
distinguish the cases $k=1$ and $k\geq 2$. The distance is thereby measured by
means of the so-called Kantorovich-Rubinstein distance which for two finite
simple counting measures $\mu_{1}$ and $\mu_{2}$ on ${\mathbb{H}}^{d}$ is
given by

$\mathbf{d_{KR}}(\mu_{1},\mu_{2}):=\sup_{h}\big{(}{\mathbb{E}}[h(\mu_{1})]-{\mathbb{E}}[h(\mu_{2})]\big{)},$

where the supremum is taken over all measurable $1$-Lipschitz functions with
respect to the total variation distance on the space of finite simple counting
measures on ${\mathbb{H}}^{d}$.

######  Theorem 1.

Let $\zeta$ be an inhomogeneous Poisson process on ${\mathbb{R}}$ with
intensity measure $\mathbb{E}\zeta$ given by
$\mathbb{E}\zeta((u,\infty))=e^{-u},\,{u\in\mathbb{R}}$. Let
${c\in\mathbb{R}}$.

  * (i)

Suppose that $k=1$. Then there are constants $C_{1,d},R_{1,d}>0$ only
depending on $d$ and $c$ such that for all $R\geq R_{1,d}$,

$\displaystyle\mathbf{d_{KR}}(\xi_{R}^{(1)}\cap(c,\infty),\zeta\cap(c,\infty))\leq\begin{cases}C_{1,d}\,Re^{-R(d-1)/2}&:d\leq
5\\\ C_{1,d}\,e^{-2R}&:d\geq 6.\end{cases}$

  * (ii)

Suppose that $k\geq 2$. Then there are constants $C_{k,d},R_{k,d}>0$ only
depending on $d$, $k$ and on $c$ such that for all $R\geq R_{k,d}$,

$\displaystyle\mathbf{d_{KR}}(\xi_{R}^{(k)}\cap(c,\infty),\zeta\cap(c,\infty))\leq{C_{k,d}{\log
R\over R}}.$

In particular, for any $k\in{\mathbb{N}}$ the point process $\xi_{R}^{(k)}$
converges in distribution to the Poisson process $\zeta$, as $R\to\infty$.

######  Remark 2.

  * (i)

As a generalization of Theorem 1, one can prove that the marked point
processes

$\displaystyle\sum_{x\in\eta\cap
B_{R}}\delta_{(x,\mathcal{H}^{d}(B_{r(x,k,\eta-\delta_{x})})-v_{k}(R))},\qquad
R>0,$

restricted to some interval $(c,\infty)$ converge to a Poisson process on the
product space $\mathbb{H}^{d}\times\mathbb{R}$, as $R\to\infty$.

  * (ii)

The distinction between $k=1$ and $k\geq 2$ and the qualitatively different
results in these cases reflect, in a sense, the growth of $v_{k}(R)-R(d-1)$.
While $v_{1}(R)-R(d-1)$ is constant in $R>0$, we find that $v_{k}(R)-R(d-1)$
grows logarithmically in $R$ for $k\geq 2$.

  * (iii)

We leave it as an open problem to decide whether (or not) the bounds in
Theorem 1 are optimal. However, we remark at this point that the bound in
Theorem 1 for $k\geq 2$ are in accordance with the bounds of the analogous
problem in a Euclidean space (see [4, Theorem 6.4] and [11, Theorem 1.2]).
Indeed, note that by Lemma 4 the statement of Theorem 1(ii) is equivalent to

$\displaystyle\mathbf{d_{KR}}(\xi_{R}^{(k)}\cap(c,\infty),\zeta\cap(c,\infty))$
$\displaystyle\leq\tilde{C}_{k,d}\,{\log\log\mathcal{H}^{d}(B_{R})\over\log\mathcal{H}^{d}(B_{R})},\qquad
R>R_{k,d},$

with some positive constants $\tilde{C}_{k,d}$ depending on $c$, $k$ and $d$.

  * (iv)

If we replace $v_{1}(R)$ in (2.1) by
$\tilde{v}_{1}(R):=\log\mathcal{H}^{d}(B_{R})$, then the bound in Theorem 1(i)
can be improved to

$\mathbf{d_{KR}}(\xi_{R}^{(1)}\cap(c,\infty),\zeta\cap(c,\infty))\leq\tilde{C}_{1,d}\,Re^{-R(d-1)/2},\qquad
R>\tilde{R}_{1,d},$

or equivalently,

$\displaystyle\mathbf{d_{KR}}(\xi_{R}^{(1)}\cap(c,\infty),\zeta\cap(c,\infty))\leq\hat{C}_{1,d}\,{\log\mathcal{H}^{d}(B_{R})\over\mathcal{H}^{d}(B_{R})^{1/2}},\qquad
R>\tilde{R}_{1,d},$

for all $d\geq 2$, where
$\tilde{C}_{1,d},\hat{C}_{1,d},\tilde{R}_{1,d}\in(0,\infty)$ are constants
only depending on $c$ and $d$. In this form, the bound is in accordance with
[4] if in the proof of Theorem 6.4 therein one systematically exploits that
$k=1$111In fact, in [4, Equation (6.12)] the integral term vanishes and the
two remaining terms decay exponentially in $n$. The bounds $E_{1}\leq
c_{1}n^{-1}$ and $E_{2}\leq c_{2}n^{-1}\log n$ remain unchanged, while
$E_{3,1}=0$ and $E_{3,2}\leq c_{3}n^{-1/2}\log n$ in the terminology of [4].
This eventually yields a bound of order $n^{-1/2}\log n$ for the Kantorovich-
Rubinstein distance.. We note that the numerically much more tractable
expression $v_{1}(R)$ in (2.1) can be regarded as a first-order approximation
of $\tilde{v}_{1}(R)$. It is the error in this approximation, which makes the
case distinction in Theorem 1(i) unavoidable. We decided to work with
$v_{k}(R)$ as given by (2.1) in order to ease comparison with the Euclidean
case.

The following extreme value statement for the distribution of the maximum
volume of $k$th nearest neighbour balls in $\eta\cap B_{R}$ is a direct
consequence of Theorem 1. It can be understood as the hyperbolic analogue to
results for the asymptotic distribution of maximum $k$th nearest neighbour
balls in Euclidean space, see [4], Section 6.2, and [11] for general
$k\in\mathbb{N}$ as well as [17] for an elementary proof in the special case
$k=1$.

######  Corollary 3.

Let $k\in{\mathbb{N}}$, $c\in\mathbb{R}$ and denote by $C_{k,d},R_{k,d}>0$ the
constants appearing in Theorem 1.

  * (i)

Suppose that $k=1$. Then for all $R\geq R_{1}$ we have that

$\displaystyle\Bigl{|}{\mathbb{P}}\Big{(}\max_{x\in\eta\cap
B_{R}}\mathcal{H}^{d}(B_{r(x,1,\eta-\delta_{x})})-v_{1}(R)\leq
c\Big{)}-\exp({-}e^{-c})\Bigl{|}\leq\begin{cases}C_{1,d}\,Re^{-R(d-1)/2}&:d\leq
5\\\ C_{1,d}\,e^{-2R}&:d\geq 6.\end{cases}$

  * (ii)

Suppose that $k\geq 2$. Then for all $R\geq R_{k}$ we have that

$\displaystyle\Bigl{|}{\mathbb{P}}\Big{(}\max_{x\in\eta\cap
B_{R}}\mathcal{H}^{d}(B_{r(x,k,\eta-\delta_{x})})-v_{k}(R)\leq
c\Big{)}-\exp({-}e^{-c})\Bigl{|}\leq{C_{k,d}{\log R\over R}}.$

In particular, for any $k\in{\mathbb{N}}$ the random variable
$\max_{x\in\eta\cap
B_{R}}\mathcal{H}^{d}(B_{r(x,k,\eta-\delta_{x})})-v_{k}(R)$ converges in
distribution to a Gumbel distribution, as $R\to\infty$.

###### Proof.

In view of the inequality

$\displaystyle|{\mathbb{P}}(\xi_{R}^{(k)}\cap(c,\infty)=\emptyset)-{\mathbb{P}}(\zeta\cap(c,\infty)=\emptyset)|\leq\mathbf{d_{KR}}(\xi_{R}^{(k)}\cap(c,\infty),\zeta\cap(c,\infty)),$

the claim follows directly from Theorem 1. ∎

The remaining parts of this paper are structured as follows. In Section 3 we
provide some necessary background material from hyperbolic geometry and prove
some auxiliary geometric estimates. We also rephrase there a general bound for
quantitative Poisson approximation from [4] on which the proof of Theorem 1 is
based. The latter is provided in Section 4.

##  3 Preliminaries

###  3.1 Hyperbolic geometry

In this section we collect some preliminary materials from hyperbolic
geometry, which are relevant in our context, and refer to the monograph [28]
and the survey article [8] for further information. We recall that $B(z,r)$
stands for a $d$-dimensional geodesic ball of radius $r>0$ centred at
$z\in{\mathbb{H}}^{d}$. If $z=p$ we simply write $B_{r}$ for $B(p,r)$. The
volume of $B_{r}$ is given by

$\displaystyle{\mathcal{H}}^{d}(B_{r})=\omega_{d}\int_{0}^{r}\sinh^{d-1}(u)\,\mathrm{d}u,$
(3.1)

where we recall that $\omega_{d}$ is the surface area of the
$(d-1)$-dimensional Euclidean unit sphere, see [28, Equation (3.26)]. Identity
(3.1) is a consequence of the polar integration formula in hyperbolic geometry
[9, pp. 123-125], which says that

$\displaystyle\int_{{\mathbb{H}}^{d}}f(x)\,{\mathcal{H}}^{d}({\rm
d}x)=\omega_{d}\int_{\mathbb{S}_{p}^{d-1}}\int_{0}^{\infty}\sinh^{d-1}(u)\,f(\exp_{p}(uv))\,{\rm
d}u\sigma_{p}({\rm d}v),$ (3.2)

where $\mathbb{S}_{p}^{d-1}$ is the $(d-1)$-dimensional unit sphere in the
tangent space $T_{p}$ at $p$, $\sigma_{p}$ the normalized spherical Lebesgue
measure on $\mathbb{S}_{p}^{d-1}\subset T_{p}$ and $\exp_{p}(uv)$ stands for
the point in ${\mathbb{H}}^{d}$ arising by applying the exponential map
$\exp_{p}:T_{p}\to{\mathbb{H}}^{d}$ to the point $uv\in T_{p}$. In particular,
we have the following bounds for ${\mathcal{H}}^{d}(B_{r})$.

######  Lemma 4.

Let $d\geq 2$. Then there are constants $\Gamma_{d},\gamma_{d}>0$ only
depending on $d$ such that
$\gamma_{d}e^{r(d-1)}\leq{\mathcal{H}}^{d}(B_{r})\leq\Gamma_{d}e^{r(d-1)}$ for
all $r\geq 2$.

###### Proof.

It is elementary to check that $\sinh(u)\geq u$ for any $u\geq 0$ and that
$\sinh(u)\geq e^{u}/3$ for $u\geq 1$. Since $d\geq 2$ we conclude that

$\displaystyle{\mathcal{H}}^{d}(B_{r})$
$\displaystyle\geq\omega_{d}\int_{0}^{1}u^{d-1}\,\mathrm{d}u+{\omega_{d}\over
3^{d-1}}\int_{1}^{r}e^{(d-1)u}\,\mathrm{d}u$ $\displaystyle={\omega_{d}\over
d}+{\omega_{d}\over(d-1)3^{d-1}}(e^{r(d-1)}-e^{d-1})$
$\displaystyle\geq{\omega_{d}\over(d-1)3^{d-1}}(e^{r(d-1)}-e^{d-1})$
$\displaystyle\geq\gamma_{d}e^{r(d-1)}$

with the choice $\gamma_{d}={\omega_{d}\over 2(d-1)3^{d-1}}$, where for the
last inequality we used that $r\geq 2$. On the other hand, $\sinh(u)\leq
e^{u}/2$ for all $u\geq 0$ and we obtain

$\displaystyle{\mathcal{H}}^{d}(B_{r})\leq{\omega_{d}\over
2^{d-1}}\int_{0}^{r}e^{(d-1)u}\,\mathrm{d}u\leq\Gamma_{d}e^{r(d-1)}$

with $\Gamma_{d}={\omega_{d}\over(d-1)2^{d-1}}$. ∎

The following lemma is an essential ingredient of the proof of Theorem 1. It
provides a bound for the volume of the difference of two nearby hyperbolic
balls with the same radius.

Figure 2: Illustration in the Beltrami-Klein model for the hyperbolic plane of
the argument used in the proof of Lemma 5. Shown are the two balls $B(x,r)$
and $B(z,r)$ with $d_{h}(x,z)<r$, as well as two perpendiculars through points
of $B_{r}(z,x)$.

######  Lemma 5.

Let $x,z\in\mathbb{H}^{d}$ and $0<s:=d_{h}(x,z)\leq r$ with $r-s/2\geq 2$.

  * (i)

It holds that

$\displaystyle\alpha_{1}se^{(d-{1})(r-s/2)}\leq\mathcal{H}^{d}(B(z,r)\setminus
B(x,r))\leq\alpha_{2}se^{(d-{1})r},$

where $\alpha_{1},\alpha_{2}\in(0,\infty)$ are constants only depending on
$d$.

  * (ii)

For $s>0$ we have

$\displaystyle\omega_{d}\left[\frac{e^{s(d-1)}}{(d-1)2^{d-1}}-\frac{(d-1)e^{s(d-3)}}{(d-3)2^{d-1}}\right]\leq\mathcal{H}^{d}(B_{s})\leq\omega_{d}\frac{e^{s(d-1)}}{(d-1)2^{d-1}}.$
(3.3)

###### Proof.

We start by observing that the boundary of the intersection $B(z,r)\cap
B(x,r)$ is a $(d-2)$-dimensional sphere of radius
$\overline{r}:=\operatorname{arcosh}\left(\frac{\cosh(r)}{\cosh(s/2)}\right)$
according to [28, Theorem 3.5.3]. Let $B_{r}(z,x)$ be the corresponding
$(d-1)$-dimensional ball. For each $y\in B_{r}(z,x)$ let $L(y)$ be the
hyperbolic line through $y$ which is orthogonal to the hyperbolic hyperplane
containing $B_{r}(z,x)$. By construction, the set $(B(z,r)\setminus
B(x,r))\cap L(y)$ is a hyperbolic segment of length $r-(r-s)=s$ and it follows
that

$\displaystyle\mathcal{H}^{d}(B(z,r)\setminus B(x,r))$
$\displaystyle\geq\int_{B_{r}(z,x)}{\mathcal{H}}^{1}((B(z,r)\setminus
B(x,r))\cap L(y))\,\cosh(d_{h}(L(y),m))\,{\mathcal{H}}^{d-1}(\mathrm{d}y)$
$\displaystyle=s\int_{B_{r}(z,x)}\cosh(d_{h}(L(y),m))\,{\mathcal{H}}^{d-1}(\mathrm{d}y),$

see Figure 2, where $m$ stands for the centre of $B_{r}(z,x)$ and
$d_{h}(L(y),m)$ for the hyperbolic distance from $L(y)$ to $m$.

To derive a lower bound for the integral we use that $\overline{r}\geq r-s/2$
according to [20, Lemma 6] and our assumption. Using the polar integration
formula (3.2) within the hyperbolic hyperplane containing $B_{r}(x,z)$ this
leads to

$\displaystyle\int_{B_{r}(z,x)}\cosh(d_{h}(L(y),m))\,{\mathcal{H}}^{d-1}(\mathrm{d}y)={\omega_{d-1}}\int_{0}^{r-s/2}\cosh(t)\,\sinh^{d-2}(t)\,{\rm
d}t,$

which implies the asserted lower bound for $d=2$, since $\cosh(t)$ is the
derivative of $\sinh(t)$ and $\sinh(t)\geq e^{t-3}$ for $t\geq 2$. For $d\geq
3$ we use now additionally that $\cosh(t)\geq e^{t}/2$ for all $t\geq 0$ and
that $\sinh(t)\geq t$ for all $0\leq t\leq 2$. This gives

$\displaystyle\mathcal{H}^{d}(B(z,r)\setminus B(x,r))$ $\displaystyle\geq
s\omega_{d-1}\Big{(}\int_{2}^{r-s/2}\cosh(t)\,\sinh^{d-2}(t)\,{\rm
d}t+\int_{0}^{2}\cosh(t)\,\sinh^{d-2}(t)\,{\rm d}t\Big{)}$
$\displaystyle\geq\frac{s\omega_{d-1}}{2}\Big{(}e^{-3(d-2)}\int_{2}^{r-s/2}e^{(d-1)t}\,{\rm
d}t+\int_{0}^{2}t^{d-2}\,{\rm d}t\Big{)}$
$\displaystyle=\frac{s\omega_{d-1}}{2(d-1)}\Big{(}\frac{e^{(d-1)(r-s/2)}}{e^{3(d-2)}}-e^{-d+4}+2^{d-1}\Big{)}$
$\displaystyle\geq\alpha_{1}se^{(d-1)(r-s/2)},$

where we used the fact that $2^{d-1}-e^{-d+4}>0$ for $d\geq 3$.

To obtain an upper bound, let $\tilde{B}_{r}(z,x)$ be the $(d-1)$-dimensional
ball with radius $r$ that is contained in the hyperbolic hyperplane through
$B_{r}(z,x)$. We use the polar integration formula (3.2) in the hyperbolic
hyperplane containing $B_{r}(z,x)$ and obtain

$\displaystyle\mathcal{H}^{d}(B(z,r)\setminus B(x,r))$ $\displaystyle\leq
s\int_{\tilde{B}_{r}(z,x)}\cosh(d_{h}(L(y),m))\,{\mathcal{H}}^{d-1}(\mathrm{d}y)$
$\displaystyle=s{\omega_{d-1}}\int_{0}^{r}\cosh(t)\,\sinh^{d-2}(t)\,{\rm d}t$
$\displaystyle\leq{s\omega_{d-1}\over 2^{d-2}}\int_{0}^{r}e^{(d-1)t}\,{\rm
d}t$ $\displaystyle={s\omega_{d-1}\over(d-1)2^{d-2}}(e^{(d-1)r}-1)$
$\displaystyle\leq\alpha_{2}se^{(d-1)r},$

where we additionally applied the inequalities $\cosh(t)\leq e^{t}$ and
$\sinh(t)\leq e^{t}/2$ for all $t\geq 0$. We thus conclude the proof of part
(i). The inequalities in (ii) follow from (3.1) and the fact that
$\sinh(t)=\frac{e^{t}-e^{-t}}{2}$ for $t\in{\mathbb{R}}$. ∎

###  3.2 Poisson approximation

In this section we rephrase a special case of [4, Theorem 4.1] which we will
use to prove Theorem 1. We work with a Polish space $\mathbb{X}$ and denote by
${\bf N}_{\mathbb{X}}$ the space of locally finite, simple counting measures
on $\mathbb{X}$. As usual, we identify each element in ${\bf N}_{\mathbb{X}}$
with its support. Let $f:\mathbb{X}\times{\bf N}_{\mathbb{X}}\to{\mathbb{R}}$
and $g:\mathbb{X}\times{\bf N}_{\mathbb{X}}\to\\{0,1\\}$ be a measurable
functions and define for $\mu\in{\bf N}_{\mathbb{X}}$ the random point process

$\xi[\mu]:=\sum_{x\in\mu}g(x,\mu)\delta_{f(x,\mu)},$ (3.4)

whose intensity measure is denoted by ${\mathbb{E}}\xi[\mu](\,\cdot\,)$. Let
$\mathcal{F}$ denote the system of closed sets in $\mathbb{X}$ endowed with
the Fell topology, see [23, p. 256]. To each point $x\in\mathbb{X}$ we
associate in a measurable way a stopping set
$\mathcal{S}(x,\,\cdot\,):\mathbf{N}_{\mathbb{X}}\to\mathcal{F}$ as well as a
closed set $S_{x}\subset\mathbb{X}$ with $x\in S_{x}$, where we recall that
the stopping property of $S(x,\,\cdot\,)$ means that
$\\{\mu\in\mathbf{N}_{\mathbb{X}}:S(x,\mu)\subseteq
K\\}=\\{\mu\in\mathbf{N}_{\mathbb{X}}:S(x,\mu\cap K)\subseteq K\\}$ for all
compact $K\subseteq\mathbb{X}$. It is assumed that $f$ and $g$ are localized
in the sense that for $\mathcal{S}(x,\mu)\subset S_{x}$,

$\displaystyle g(x,\mu)$ $\displaystyle=g(x,\mu\cap S_{x})$ $\displaystyle
f(x,\mu)$ $\displaystyle=f(x,\mu\cap S_{x})\quad\text{if}\quad g(x,\mu)=1.$

To rephrase the result from [4] we need the total variation distance
$\mathbf{d_{TV}}(\nu_{1},\nu_{2})$ between two measures $\nu_{1},\nu_{2}$ on
$\mathbb{X}$, which is defined as

$\mathbf{d_{TV}}(\nu_{1},\nu_{2}):=\sup_{B}|\nu_{1}(B)-\nu_{2}(B)|,$

where the supremum runs over all Borel subsets $B$ of $\mathbb{X}$ satisfying
$\nu_{1}(B),\nu_{2}(B)<\infty$. Now, let $\eta$ be a Poisson process on
$\mathbb{X}$ with intensity measure ${\mathbb{E}}\eta(\,\cdot\,)$, and define
the quantities

$\displaystyle E_{1}$
$\displaystyle:=\int_{\mathbb{X}}{\mathbb{E}}[g(x,\eta+\delta_{x})\mathds{1}\\{\mathcal{S}(x,\eta)\not\subset
S_{x}\\}]\,{\mathbb{E}}\eta({\rm d}x),$ $\displaystyle E_{2}$
$\displaystyle:=\int_{\mathbb{X}}\int_{\mathbb{X}}\mathds{1}\\{S_{x}\cap
S_{z}\neq\varnothing\\}{\mathbb{E}}[g(x,\eta+\delta_{x})]{\mathbb{E}}[g(z,\eta+\delta_{z})]\,{\mathbb{E}}\eta({\rm
d}z){\mathbb{E}}\eta({\rm d}x),$ $\displaystyle E_{3}$
$\displaystyle:=\int_{\mathbb{X}}\int_{\mathbb{X}}\mathds{1}\\{S_{x}\cap
S_{z}\neq\varnothing\\}{\mathbb{E}}[g(x,\eta+\delta_{x}+\delta_{z})g(z,\eta+\delta_{x}+\delta_{z})]\,{\mathbb{E}}\eta({\rm
d}z){\mathbb{E}}\eta({\rm d}x),$

where $\delta_{(\cdot)}$ denotes the Dirac measure. Assuming finally that
${\mathbb{E}}\xi[\eta]({\mathbb{R}})<\infty$, [4, Theorem 4.1] says that the
Kantorovich-Rubinstein distance $\mathbf{d_{KR}}(\xi[\eta],\zeta)$ between
$\xi[\eta]$ and a Poisson process $\zeta$ on ${\mathbb{R}}$ with finite
intensity measure ${\mathbb{E}}\eta(\,\cdot\,)$ can be estimated from above by

$\displaystyle\mathbf{d_{KR}}(\xi[\eta],\zeta)\leq\mathbf{d_{TV}}({\mathbb{E}}\xi[\eta],{\mathbb{E}}\zeta)+2(E_{1}+E_{2}+E_{3}).$
(3.5)

We refer to [4, 13] for a a dual formulation as well as more details on the
Kantorovich-Rubinstein distance between random points measures. Let us also
remark that if in (3.5) the Kantorovich-Rubinstein distance is replaced by the
total variation distance, the result without the factor $2$ on the right-hand
side can be found in [1].

##  4 Proof of Theorem 1

Our goal is to apply the Poisson approximation bound (3.5). To this end, we
need to specify $\mathbb{X}$, $\eta$, the functions $f$ and $g$ as well as the
sets $S(x,\,\cdot\,)$ and $S_{x}$. For the space $\mathbb{X}$ we take the
$d$-dimensional hyperbolic space ${\mathbb{H}}^{d}$ and for $\eta$ a Poisson
process on ${\mathbb{H}}^{d}$ with intensity measure ${\mathcal{H}}^{d}$. To
ensure finiteness of the involved intensity measures we fix some arbitrary
$c\in\mathbb{R}$ and define the two functions $g,f:{\mathbb{H}}^{d}\times{\bf
N}_{{\mathbb{H}}^{d}}\to{\mathbb{R}}$ by

$\displaystyle g(x,\mu):=\mathds{1}\\{x\in
B_{R}\\}\mathds{1}\\{\mathcal{H}^{d}(B_{r(x,k,{\mu})})-v_{k}(R)>c\\},$
$\displaystyle f(x,\mu):=\mathcal{H}^{d}(B_{r(x,k,{\mu})})-v_{k}(R),$

where we recall that $r(x,k,\mu)$ denotes the distance to the $k$th nearest
neighbour of a point $x$ in the support of a simple counting measure $\mu$.
Then the point process of exceedances (2.2) restricted to the interval
$(c,\infty)$ has the same distribution as $\xi[\eta]$ in (3.4) using the
functions $f$ and $g$ as just defined. Next, for $x\in\mathbb{H}^{d}$ let
$\mathcal{S}(x,\mu):=B(x,r(x,k,\mu))$ and $S_{x}:=B(x,r_{c^{\prime}})$ be the
closed ball with centre at $x$ and $\mathcal{H}^{d}$-measure
$c^{\prime}+v_{k}(R)>0$ for some $c^{\prime}>\max(c,0)$ to be specified below.
We emphasize that this choice implicitly determines the radius
$r_{c^{\prime}}$ via (3.1). By construction, the functions $f$ and $g$ are
localized to the sets $\mathcal{S}(x,\mu)$.

Taking $s:=r_{c^{\prime}}$ in Lemma 5(ii) we obtain

$\begin{split}\log(c^{\prime}+v_{k}(R))&\leq
r_{c^{\prime}}(d-1)+\log\left(\frac{\omega_{d}}{(d-1)2^{d-1}}\right)\\\
&\leq\log\left(\frac{(d-3)e^{2r_{c}^{\prime}}}{(d-3)e^{2r_{c}^{\prime}}-(d-1)^{2}}\right)+\log(c^{\prime}+v_{k}(R)).\end{split}$
(4.1)

Finally, we let $\zeta$ be an inhomogeneous Poisson process on ${\mathbb{R}}$
with intensity measure ${\mathbb{E}}\zeta((u,\infty))=e^{-u}$, $u>c$, as in
the statement of Theorem 1. We can now apply (3.5) to conclude that for all
$R>0$,

$\displaystyle\mathbf{d_{KR}}(\xi_{R}^{(k)}\cap(c,\infty),\zeta\cap(c,\infty))\leq\mathbf{d_{TV}}(\mathbb{E}\xi_{R}^{(k)}\cap(c,\infty),\mathbb{E}\zeta\cap(c,\infty))+2(E_{1}+E_{2}+E_{3}),$
(4.2)

where the error terms $E_{1}$, $E_{2}$ and $E_{3}$ are given by

$\displaystyle E_{1}$ $\displaystyle:=\int_{{\mathbb{H}}^{d}}\mathds{1}\\{x\in
B_{R}\\}{\mathbb{P}}[\mathcal{H}^{d}(B_{r(x,k,\eta)})>c+v_{k}(R),\,r(x,k,\eta)>r_{c^{\prime}}]\,\mathcal{H}^{d}(\mathrm{d}x),$
$\displaystyle E_{2}$
$\displaystyle:=\int_{{\mathbb{H}}^{d}}\int_{{\mathbb{H}}^{d}}\mathds{1}\\{x,z\in
B_{R},B(x,r_{c^{\prime}})\cap B(z,r_{c^{\prime}})\neq\emptyset\\}$
$\displaystyle\qquad\qquad\times{\mathbb{P}}[\mathcal{H}^{d}(B_{r(x,k,\eta)})>c+v_{k}(R)]{\mathbb{P}}[\mathcal{H}^{d}(B_{r(z,k,\eta)})>c+v_{k}(R)]\,\mathcal{H}^{d}(\mathrm{d}z)\mathcal{H}^{d}(\mathrm{d}x),$
$\displaystyle E_{3}$
$\displaystyle:=\int_{{\mathbb{H}}^{d}}\int_{{\mathbb{H}}^{d}}\mathds{1}\\{x,z\in
B_{R},B(x,r_{c^{\prime}})\cap B(z,r_{c^{\prime}})\neq\emptyset\\}$
$\displaystyle\qquad\qquad\times{\mathbb{P}}[\mathcal{H}^{d}(B_{r(x,k,\eta+\delta_{z})})>c+v_{k}(R),\,\mathcal{H}^{d}(B_{r(z,k,\eta+\delta_{x})})>c+v_{k}(R)]\,\mathcal{H}^{d}(\mathrm{d}z)\mathcal{H}^{d}(\mathrm{d}x).$

The remaining parts of the proof bound individually the four terms on the
right hand side of (4.2).

#### Bounding the total variation distance.

In a first step we investigate the intensity measure of $\xi_{R}^{(k)}$. Let
$r_{c}>0$ be such that $\mathcal{H}^{d}(B_{r_{c}})=v_{k}(R)+c$ and note that
$\mathcal{H}^{d}(B_{r(x,k,\mu)})>v_{k}(R)+c$ if and only if
$\mu(B(x,r_{c}))\leq{k-1}$. From the Mecke equation for Poisson processes [23,
Theorem 4.1] we obtain that for all $u>c$,

$\displaystyle\mathbb{E}\xi_{R}^{(k)}((u,\infty))$
$\displaystyle=\int_{{\mathbb{H}}^{d}}\mathds{1}\\{x\in
B_{R}\\}{\mathbb{P}}[\mathcal{H}^{d}(B_{r(x,k,\eta)})>u+v_{k}(R)]\,\mathcal{H}^{d}(\mathrm{d}x).$

Since the distribution of $\eta$ is invariant under hyperbolic isometries and
$\mathcal{H}^{d}(B_{r(x,k,\eta)})>u+v_{k}(R)$ if and only if there are at most
${k-1}$ points of $\eta$ in a ball with $\mathcal{H}^{d}$-measure $u+v_{k}(R)$
around $x$, the expression is equal to

$\displaystyle\mathcal{H}^{d}(B_{R})e^{-u-v_{k}(R)}\sum_{\ell=0}^{{k-1}}\frac{(u+v_{k}(R))^{\ell}}{\ell!}$
$\displaystyle\,=e^{-u}\mathcal{H}^{d}(B_{R})\frac{(k-1)!2^{d-1}(d-1)e^{-R(d-1)}}{\omega_{d}(R(d-1))^{k-1}}\sum_{\ell=0}^{{k-1}}\frac{\big{(}u+v_{k}(R)\big{)}^{\ell}}{\ell!},$

where we used the definition (2.1) of $v_{k}(R)$ and the Poisson property of
$\eta$. Hence, the Lebesgue density $\varrho$ of
$\mathbb{E}\xi_{R}^{(k)}\cap(c,\infty)$ is given by

$\displaystyle\varrho(u)=e^{-u}\mathcal{H}^{d}(B_{R})\frac{2^{d-1}(d-1)e^{-R(d-1)}(u+v_{k}(R))^{k-1}}{\omega_{d}(R(d-1))^{k-1}},\qquad
u>c.$

Now, (3.3) gives for $R>0$ large enough

$\displaystyle\textbf{d}_{\textbf{TV}}(\mathbb{E}\xi_{R}^{(k)}\cap(c,\infty),\mathbb{E}\zeta\cap(c,\infty))$
$\displaystyle\leq\int_{c}^{\infty}|\varrho(u)-e^{-u}|\,\mathrm{d}u$
$\displaystyle=\int_{c}^{\infty}\Big{|}e^{-u}\mathcal{H}^{d}(B_{R})\frac{2^{d-1}(d-1)e^{-R(d-1)}(u+v_{k}(R))^{k-1}}{\omega_{d}(R(d-1))^{k-1}}-e^{-u}\Big{|}\,\mathrm{d}u$
$\displaystyle\leq
e^{-c}\Big{[}(1+\beta_{1}e^{-2R})\Big{(}1+\frac{c+\beta_{2}\log
R}{R(d-1)}\Big{)}^{k-1}-1\Big{]}$
$\displaystyle\leq\begin{cases}\beta_{1}e^{-c}e^{-2R}&:k=1\\\
\beta_{3}e^{-c}(c+\beta_{2}\log R)R^{-1}&:k\geq 2\end{cases}$ (4.3)

for constants $\beta_{1},\beta_{2},\beta_{3}>0$ only depending on $k$ and $d$.

#### Bounding $\mathbf{E}_{1}$.

Note that $r(x,k,\eta)>r_{c^{\prime}}$ if and only if
$\mathcal{H}^{d}(B_{r(x,k,\eta)})>c^{\prime}+v_{k}(R)$. Hence, we find from
(3.3) similarly to the estimate leading to (4.3) that for all
$c^{\prime}>\max(c,0)$,

$\displaystyle
E_{1}=\mathbb{E}\xi_{R}^{(k)}\cap(c^{\prime},\infty)\leq\begin{cases}(1+\beta_{1}e^{-2R})e^{-c^{\prime}}&:k=1\\\
\beta_{4}e^{-c^{\prime}}\Big{(}\frac{c^{\prime}+v_{k}(R)}{R(d-1)}\Big{)}^{k-1}&:k\geq
2\end{cases}$ (4.4)

where the constant $\beta_{4}>0$ depends only on $k$.

#### Bounding $\mathbf{E}_{2}$.

Since $B(x,r_{c})\cap B(z,r_{c})\neq\emptyset$ if and only if the hyperbolic
distance of $x$ and $z$ is at most $2r_{c}$, we obtain from the invariance of
$\eta$ under hyperbolic isometries that $E_{2}$ is bounded by

$\displaystyle\mathbb{E}\xi_{R}^{(k)}((c,\infty))\mathbb{P}[\mathcal{H}^{d}(B_{r(p,k,\eta)})>c+v_{k}(R)]\int_{{\mathbb{H}}^{d}}\mathds{1}\\{B(p,r_{c^{\prime}})\cap
B(z,r_{c^{\prime}})\neq\emptyset\\}\,\mathcal{H}^{d}({\rm d}z)$
$\displaystyle\quad\leq(\mathbb{E}\xi_{R}^{(k)}((c,\infty)))^{2}\frac{\mathcal{H}^{d}(B_{2r_{c^{\prime}}})}{\mathcal{H}^{d}(B_{R})}.$

From (4.1) we have for $R>0$ large enough that

$r_{c^{\prime}}\leq\frac{1}{d-1}\log(c^{\prime}+v_{k}(R))+\beta_{5}$ (4.5)

for some constant $\beta_{5}>0$ only depending on $d$. Hence, using the
definition of $v_{k}(R)$, there is another constant $\beta_{6}>0$ only
depending on $d$ such that for all $R>0$,

$\displaystyle E_{2}\leq\beta_{6}(c^{\prime}+R)^{2}e^{-R(d-1)}.$ (4.6)

#### Bounding $\mathbf{E}_{3}$.

First we consider the case $k=1$. Note that $(\eta-\delta_{x})(B(x,r_{c}))=0$
and $(\eta-\delta_{z})(B(z,r_{c}))=0$ implies for $x,z\in\eta$ that
$d_{h}(z,x)\geq r_{c}$. Hence, we obtain for $E_{3}$ the upper bound

$\displaystyle E_{3}\leq$
$\displaystyle\int_{{\mathbb{H}}^{d}}\int_{{\mathbb{H}}^{d}}\mathds{1}\\{z\in
B_{R},x\in B_{R}\cap B(z,2r_{c^{\prime}})\setminus
B(z,r_{c})\\}{\mathbb{P}}[\eta(B(z,r_{c}))=0]$
$\displaystyle\qquad\qquad\times{\mathbb{P}}[\eta(B(z,r_{c})\setminus
B(x,r_{c}))=0]\,\mathcal{H}^{d}(\mathrm{d}z)\mathcal{H}^{d}(\mathrm{d}x)$
$\displaystyle\leq\mathbb{E}\xi_{R}^{(1)}((c,\infty))\int_{\mathbb{H}^{d}}\mathds{1}\\{x\in
B(p,2r_{c^{\prime}})\setminus
B(p,r_{c})\\}{\mathbb{P}}[\eta(B(z,r_{c})\setminus
B(x,r_{c}))=0]\,\mathcal{H}^{d}(\mathrm{d}x).$ (4.7)

Let $z\in{\mathbb{H}}^{d}\setminus B(p,r_{c})$ and observe that
$\mathcal{H}^{d}(B(z,r_{c})\setminus
B(p,r_{c}))\geq\mathcal{H}^{d}(B_{r_{c}})/2$. Thus,

$\displaystyle{\mathbb{P}}[\eta(B(z,r_{c})\setminus
B(x,r_{c}))=0]=\exp(-\mathcal{H}^{d}(B(z,r_{c})\setminus
B(x,r_{c})))\leq\exp(-\mathcal{H}^{d}(B_{r_{c}})/2).$

This gives for (4.7) the bound

$\displaystyle\mathbb{E}\xi_{R}^{(1)}((c,\infty))\mathcal{H}^{d}(B_{2r_{c^{\prime}}})e^{-R(d-1)/2}\leq\beta_{7}e^{2(d-1)r_{c^{\prime}}}e^{-R(d-1)/2}$
(4.8)

for some constant $\beta_{7}>0$ only depending on $d$.

Next we consider the case $k\geq 2$ and let $a\in(0,1]$, its precise value
will be specified later. In order to bound $E_{3}$ we distinguish the
situations that $z\in B(x,ar_{c^{\prime}})$ and $z\notin
B(x,ar_{c^{\prime}})$. In the first case we have that
$\mathcal{H}^{d}(B_{r(x,k,\eta+\delta_{z})})>v_{k}(R)+c$ and
$\mathcal{H}^{d}(B_{r(z,k,\eta+\delta_{x})})>v_{k}(R)+c$ if and only if
$\eta(B(x,r_{c}))\leq{k-2}$ and $\eta(B(z,r_{c}))\leq{k-2}$. This allows us to
bound $E_{3}$ by

$\displaystyle\int_{{\mathbb{H}}^{d}}\int_{{\mathbb{H}}^{d}}\mathds{1}\\{x\in
B_{R},z\in B_{R}\cap
B(x,{ar_{c^{\prime}}})\\}{\mathbb{P}}[\eta(B(x,r_{c}))\leq{k-2}]$
$\displaystyle\qquad\qquad\times{\mathbb{P}}[\eta(B(z,r_{c})\setminus
B(x,r_{c}))\leq{k-2}]\,\mathcal{H}^{d}({\rm d}z)\mathcal{H}^{d}({\rm d}x)$
(4.9)
$\displaystyle+\int_{{\mathbb{H}}^{d}}\int_{{\mathbb{H}}^{d}}\mathds{1}\\{x\in
B_{R},z\in B_{R}\cap(B(x,2r_{c^{\prime}})\setminus
B(x,{ar_{c^{\prime}}}))\\}{\mathbb{P}}[\eta(B(x,r_{c}))\leq{k-1}]$
$\displaystyle\qquad\qquad\times{\mathbb{P}}[\eta(B(z,r_{c})\setminus
B(x,r_{c}))\leq{k-1}]\,\mathcal{H}^{d}({\rm d}z)\mathcal{H}^{d}({\rm d}x).$
(4.10)

By invariance of $\eta$ under hyperbolic isometries, the first term (4.9) is
equal to

$\displaystyle\mathbb{E}\xi_{R}^{(k)}((c,\infty))\frac{{\mathbb{P}}[\eta(B_{r_{c}})\leq{k-2}]}{{\mathbb{P}}[\eta(B_{r_{c}})\leq{k-1}]}\int_{{\mathbb{H}}^{d}}\mathds{1}\\{z\in
B_{ar_{c^{\prime}}}\\}{\mathbb{P}}[\eta(B(z,r_{c})\setminus
B_{r_{c}})\leq{k-2}]\,\mathcal{H}^{d}(\mathrm{d}z)$
$\displaystyle\quad=\omega_{d}\mathbb{E}\xi_{R}^{(k)}((c,\infty))\frac{{\mathbb{P}}[\eta(B_{r_{c}})\leq{k-2}]}{{\mathbb{P}}[\eta(B_{r_{c}})\leq{k-1}]}$
$\displaystyle\quad\quad\times\sum_{\ell=0}^{{k-2}}\int_{0}^{{ar_{{c^{\prime}}}}}\sinh^{d-1}(s)\exp\left(-{\mathcal{H}}^{d}(B(z,r_{c})\setminus
B_{r_{c}})\right)\frac{({\mathcal{H}}^{d}(B(z,r_{c})\setminus
B_{r_{c}}))^{\ell}}{\ell!}\,\mathrm{d}s,$ (4.11)

where we applied the polar integration formula (3.2) with $z:=\exp_{p}(sv)$
for some arbitrary $v\in\mathbb{S}_{p}^{d-1}$ in (4.11). Note that we always
have that $2(r_{c}-2)>r_{c}\geq s$ if $r_{c}\geq 5$, which can be achieved by
choosing the parameter $R$ large enough. Thus, for such $R$ we find by Lemma 4
and Lemma 5(i) that

$\displaystyle\int_{0}^{{ar_{{c^{\prime}}}}}\sinh^{d-1}(s)\exp\left(-{\mathcal{H}}^{d}(B(z,r_{c})\setminus
B_{r_{c}})\right)\frac{({\mathcal{H}}^{d}(B(z,r_{c})\setminus
B_{r_{c}}))^{\ell}}{\ell!}\,\mathrm{d}s$ $\displaystyle\leq{1\over
2^{d-1}}\int_{0}^{{ar_{c^{\prime}}}}e^{(d-1)s}\,\exp\Big{(}-{\alpha_{1}se^{(d-1)(r_{c}-{s\over
2})}}\Big{)}{(\alpha_{2}se^{(d-1)r_{c}})^{\ell}\over\ell!}\,\mathrm{d}s$
$\displaystyle\leq{1\over
2^{d-1}}{(\alpha_{2}e^{(d-1)r_{c}})^{\ell}\over\ell!}\int_{0}^{\infty}\exp\Big{(}-s\Big{(}\alpha_{1}e^{{(d-1)(2r_{c}-ar_{c^{\prime}})\over
2}}-(d-1)\Big{)}\Big{)}s^{\ell}\,\mathrm{d}s$
$\displaystyle\leq\beta_{8}{\exp\Big{(}-(d-1)\frac{2r_{c}-a(\ell+1)r_{c^{\prime}}}{2}\Big{)}},$
(4.12)

which takes the maximum value
$\beta_{8}e^{-(d-1)\frac{2r_{c}-a(k-1)r_{c^{\prime}}}{2}}$ for $\ell=k-2$,
where $\beta_{8}>0$ is a constant only depending on $d$.

Next, we recall that for a general Poisson random variable $X$ with mean
$\lambda>0$ and any bounded function $h:\\{0,1,\ldots\\}\to{\mathbb{R}}$ one
has that

$\displaystyle{\mathbb{E}}[\lambda h(X+1)]={\mathbb{E}}[Xh(X)];$

in fact, this is the famous Chen-Stein characterization of the Poisson
distribution. In particular, applying this to the function $h(n)={\bf
1}\\{n\leq k\\}$ we see that

$\displaystyle\lambda{\mathbb{P}}(X\leq k-1)={\mathbb{E}}[X{\bf 1}\\{X\leq
k\\}]\leq k{\mathbb{P}}(X\leq k).$

Applying this to the Poisson random variable $\eta(B_{r_{c}})$ which has mean
${\mathcal{H}}^{d}(B_{r_{c}})=v_{k}(R)+c$ we find from (2.1) that

$\frac{{\mathbb{P}}[\eta(B_{r_{c}})\leq{k-2}]}{{\mathbb{P}}[\eta(B_{r_{c}})\leq{k-1}]}\leq{k-1\over
v_{k}(R)+c}.$

Hence, we obtain from (4.11) and (4.12) that (4.9) is bounded by

$\displaystyle\beta_{9}\mathbb{E}\xi_{R}^{(k)}((c,\infty)){k-1\over
v_{k}(R)+c}e^{-(d-1)\frac{2r_{c}-a(k-1)r_{c^{\prime}}}{2}}\leq\beta_{10}R^{-1}e^{-(d-1)\frac{2r_{c}-a(k-1)r_{c^{\prime}}}{2}},$

where $\beta_{9},\beta_{10}>0$ are constants only depending on $d$ and $k$.

Now, we consider (4.10). By invariance of $\eta$ under hyperbolic isometries,
(4.10) is the same as

$\displaystyle\mathbb{E}\xi_{R}^{(k)}((c,\infty))\int_{{\mathbb{H}}^{d}}\mathds{1}\\{z\in
B_{2r_{c^{\prime}}}\setminus
B_{ar_{c^{\prime}}}\\}{\mathbb{P}}[\eta(B(z,r_{c})\setminus
B_{r_{c}})\leq{k-1}]\,\mathcal{H}^{d}(\mathrm{d}z)$
$\displaystyle=\omega_{d}\mathbb{E}\xi_{R}^{(k)}((c,\infty))\sum_{\ell=0}^{{k-1}}\int_{ar_{c^{\prime}}}^{2r_{c^{\prime}}}\sinh^{d-1}(s)\exp\left(-{\mathcal{H}}^{d}(B(z,r_{c})\setminus
B_{r_{c}})\right)\frac{({\mathcal{H}}^{d}(B(z,r_{c})\setminus
B_{r_{c}}))^{\ell}}{\ell!}\,\mathrm{d}s,$

where we again applied the polar integration formula (3.2) with
$z:=\exp_{p}(sv)$ for some arbitrary $v\in\mathbb{S}_{p}^{d-1}$. Using that
${\mathcal{H}}^{d}(B(z,r_{c^{\prime}})\setminus
B_{r_{c}})\geq{\mathcal{H}}^{d}(B_{s/2})$ for $d_{h}(p,z)=s$, we find from
Lemma 4 that ${\mathcal{H}}^{d}(B(z,r_{c})\setminus
B_{r_{c}})\geq\gamma_{d}e^{(d-1)s/2}$ (the result can indeed be applied, since
$s/2\geq ar_{c}/2\geq 2$ for large enough $R$), see Figure 3. Since the
function $u\mapsto e^{-u}u^{\ell}$ is decreasing for $u\geq u_{0}$ large
enough, we obtain for $M>0$ to be specified below and $R$ large enough,

$\displaystyle\int_{ar_{c^{\prime}}}^{2r_{c^{\prime}}}\sinh^{d-1}(s)\exp\left(-{\mathcal{H}}^{d}(B(z,r_{c})\setminus
B_{r_{c}})\right)\frac{({\mathcal{H}}^{d}(B(z,r_{c})\setminus
B_{r_{c}}))^{\ell}}{\ell!}\,\mathrm{d}s$
$\displaystyle\quad\leq\frac{1}{2^{d-1}}\int_{ar_{c^{\prime}}}^{2r_{c^{\prime}}}e^{(d-1)s}\exp\Big{(}-\gamma_{d}se^{(d-1)s/2}\Big{)}\frac{(\gamma_{d}se^{(d-1)s/2})^{\ell}}{\ell!}\mathrm{d}s$
$\displaystyle\quad\leq\frac{e^{-r_{c^{\prime}}\frac{M(d-1)a}{2}}}{2^{d-2}(d-1)(ar_{c^{\prime}})^{M+2}}\int_{ar_{c^{\prime}}}^{2r_{c^{\prime}}}\Big{(}1+{s(d-1)\over
2}\Big{)}e^{(d-1)s/2}\exp\Big{(}-\gamma_{d}se^{(d-1)s/2}\Big{)}\frac{\gamma_{d}^{\ell}(se^{(d-1)s/2})^{\ell+1+M}}{\ell!}\mathrm{d}s$
$\displaystyle\quad\leq\beta_{11}r_{c^{\prime}}^{-M-2}e^{-r_{c^{\prime}}\frac{M(d-1)a}{2}},$

where we have used the substitution $t=se^{(d-1)s/2}$ for the last inequality
and $\beta_{11}>0$ is a constant depending on $d$ and on $M$ .

Figure 3: Illustration in the Beltrami-Klein model for the hyperbolic plane of
the argument used in the estimate of (4.10). Shown are the balls $B(p,r_{c})$
and $B(z,r_{c})$, where $d_{h}(p,z)=r_{c}$. The blue ball has radius
$r_{c}/2$.

Summarizing, we have shown that for large enough $R$,

$\displaystyle
E_{3}\leq\begin{cases}\beta_{7}e^{2(d-1)r_{c^{\prime}}}\exp(-\alpha_{1}r_{c}e^{(d-1)r_{c}/2})&:k=1\\\
\beta_{10}R^{-1}{\exp\left(-(d-2)\frac{2r_{c}-a(k-1)r_{c^{\prime}}}{2}\right)}+{\beta_{11}e^{-r_{c^{\prime}}\frac{M(d-1)a}{2}}}&:k\geq
2.\end{cases}$ (4.13)

Thus, choosing $a:={r_{c}\over{kr_{c^{\prime}}}}\in(0,1]$ and $M>0$ so large
that $Ma>2$, we find for $k\geq 2$ that

$\displaystyle E_{3}\leq\beta_{12}R^{-1},$ (4.14)

where $\beta_{12}$ is a constant depending on $c$, $d$ and on $k$.

#### Completing the proof.

After having estimated all terms in (4.2) we first conclude for $k=1$ from
(4.3), (4.4), (4.6) and (4.13) that for all sufficiently large $R$,

$\displaystyle\mathbf{d_{KR}}(\xi_{R}^{(1)}\cap(c,\infty),\zeta\cap(c,\infty))$
$\displaystyle\quad\leq\beta_{1}e^{-c}e^{-2R}+2\Big{\\{}(1+\beta_{1}e^{-2R})e^{-c^{\prime}}+\beta_{9}(c^{\prime}+R)^{2}e^{-R(d-1)}+\beta_{7}e^{2(d-1)r_{c^{\prime}}}e^{-R(d-1)/2}\Big{\\}}.$

Choosing $c^{\prime}=c+\log R$ and using (4.5) we see that

$\displaystyle\mathbf{d_{KR}}(\xi_{R}^{(1)}\cap(c,\infty),\zeta\cap(c,\infty))\leq\begin{cases}C_{1,d}\,Re^{-R(d-1)/2}&:d\leq
5\\\ C_{1,d}\,e^{-2R}&:d\geq 6.\end{cases}$

where $C_{1,d}>0$ is some constant only depending on $c$ and on $d$, which is
the result of Theorem 1(i).

On the other hand, for $k\geq 2$ and using again (4.3), (4.4), (4.6) and this
time (4.14) we arrive at the bound

$\displaystyle\mathbf{d_{KR}}(\xi_{R}^{(k)}\cap(c,\infty),\zeta\cap(c,\infty))$
$\displaystyle\quad\leq\beta_{3}e^{-c}(c+\beta_{1}\log
R)R^{-1}+2\Big{\\{}\beta_{4}e^{-c^{\prime}}\Big{(}\frac{c^{\prime}+v_{k}(R)}{R(d-1)}\Big{)}^{k-1}+\beta_{9}(c^{\prime}+R)^{2}e^{-R(d-1)}+\beta_{12}R^{-1}\Big{\\}}.$

Choosing again $c^{\prime}=c+\log R$ we conclude that

$\displaystyle\mathbf{d_{KR}}(\xi_{R}^{(k)}\cap(c,\infty),\zeta\cap(c,\infty))\leq
C_{k,d}\,{\log R\over R},$

where $C_{k,d}>0$ is some constant depending on $c$, $k$ and on $d$. This is
the assertion of Theorem 1(ii) and completes the proof. ∎

## Acknowledgements

The authors wish to thank two anonymous referees for their helpful suggestions
which led to an improvement of this article. This work has been supported by
the DFG priority program SPP 2265 Random Geometric Systems.

## References

  * [1] Barbour, A.D. and Brown, T.C. (1992). Stein’s method and point process approximation. Stochastic Process. Appl. 43 9–31.
  * [2] Besau, F., Rosen, D. and Thäle, C. (2021). Random inscribed polytopes in projective geometries. Math. Ann. 381, 1345–1372.
  * [3] Besau, F. and Thäle, C. (2020). Asymptotic normality for random polytopes in non-Euclidean geometries. Trans. Amer. Math. Soc. 373, 8911–8941.
  * [4] Bobrowski, O., Schulte, M. and Yogeshwaran, D. (2022). Poisson process approximation under stabilization and Palm coupling. Ann. Henri Lebesgue 5, 1489–1534.
  * [5] Bode, M., Fountoulakis, N. and Müller, T. (2015). On the largest component of a hyperbolic model of complex networks. Electron. J. Combinat. 22, 3–24.
  * [6] Bonnet, G. and Chenavier, N. (2020). The maximal degree in a Poisson-Delaunay graph. Bernoulli 26, 948–979.
  * [7] Calka, P. and Chenavier, N. (2014). Extreme values for characteristic radii of a Poisson-Voronoi tessellation. Extremes 17, 359–385.
  * [8] Cannon, J.W., Floyd, W.J., Kenyon, R. and Parry, W.R. (1997). Hyperbolic geometry. In Flavors of Geometry, MSRI Publications.
  * [9] Chavel, I. (1993). Riemannian Geometry - A Modern Introduction. Cambridge University Press.
  * [10] Chenavier, N. and Hemsley, R. (2016). Extremes for the inradius in the Poisson line tessellation. Adv. in Appl. Probab. 48, 544–573.
  * [11] Chenavier, N., Henze, N. and Otto, M. (2022). Limit laws for large $k$th-nearest neighbor balls. J. Appl. Probab. 59(3), 880–894.
  * [12] Chenavier, N. and Robert, C.Y. (2018). Cluster size distributions of extreme values for the Poisson-Voronoi tessellation. Ann. Appl. Probab. 28, 3291–3323.
  * [13] Decreusefond, L., Schulte, M. and Thäle, C. (2016). Functional Poisson approximation in Kantorovich–Rubinstein distance with applications to U-statistics and stochastic geometry. Ann. Prob. 44, 2147–2197.
  * [14] Fountoulakis, N. van der Hoorn, P., Müller, T. and Schepers, M. (2021). Clustering in a hyperbolic model of complex networks. Electron. J. Probab. 26, article 13.
  * [15] Fountoulakis, N. and Müller, T. (2018). Law of large numbers for the largest component in a hyperbolic model of complex networks. Ann. Appl. Probab. 28, 607–650.
  * [16] Godland, T., Kabluchko, Z. and Thäle, C. (2022). Beta-star polytopes and hyperbolic stochastic geometry. Adv. Math. 404, article 108382, 69 pages.
  * [17] Györfi, L., Henze, N. and Walk, H. (2019). The limit distribution of the maximum probability nearest-neighbour ball. J. Appl. Probab. 56, 574–589.
  * [18] Henze, N. (1982). The limit distribution for maxima of ’weighted’ rth-nearest-neighbour distances. J. Appl. Probab. 19, 344–354.
  * [19] Henze, N. (1983). Ein asymptotischer Satz über den maximalen Minimalabstand von unabhängigen Zufallsvektoren mit Anwendung auf einen Anpassungstest im $\mathbb{R}^{p}$ und auf der Kugel. Metrika 30, 245–259.
  * [20] Herold, F., Hug, D. and Thäle, C. (2021). Does a central limit theorem hold for the k-skeleton of Poisson hyperplanes in hyperbolic space? Probab. Theory Relat. Fields 179, 889–968.
  * [21] Isokawa, Y. (2000). Poisson-Voronoi tessellations in three-dimensional hyperbolic spaces. Adv. Appl. Probab. 32, 548–662.
  * [22] Jammalamadaka, S.R. and Janson, S. (2015). Asymptotic distribution of the maximum interpoint distance in a sample of random vectors with a spherically symmetric distribution. Ann. Appl. Probab. 25, 3571–3591.
  * [23] Last, G. and Penrose, M. (2017). Lectures on the Poisson Process. Cambridge University Press, Cambridge.
  * [24] Otto, M. (2020). Poisson approximation of Poisson-driven point processes and extreme values in stochastic geometry. arXiv:2005.10116.
  * [25] Owada, T. and Yogeshwaran, D. (2022). Sub-tree counts on hyperbolic random geometric graphs. Adv. in Appl. Probab. 54, 1032–1069.
  * [26] Pianoforte, F. and Schulte, M. (2021). Poisson approximation with applications to stochastic geometry. Electron. J. Probab. 26, article 149, 36 pages.
  * [27] Pianoforte, F. and Schulte, M. (2022). Criteria for Poisson process convergence with applications to inhomogeneous Poisson–Voronoi tessellations. Stochastic Processes Appl. 147, 388–422.
  * [28] Ratcliffe, J.C. (2019). Foundations of Hyperbolic Manifolds. 3rd edn., Springer, Berlin.
  * [29] Schrempp, M. (2019). Limit laws for the diameter of a set of random points from a distribution supported by a smoothly bounded set. Extremes 22, 167–191.
  * [30] Schulte, M. and Thäle, C. (2012). The scaling limit of Poisson-driven order statistics with applications in geometric probability. Stochastic Process. Appl. 122, 4096–4120.
  * [31] Schulte, M. and Thäle, C. (2016). Poisson point process convergence and extreme values in stochastic geometry. In Stochastic Analysis for Poisson Point Processes Bocconi Springer Ser. 7, Bocconi Univ. Press.


# Rotationally invariant formulation of spin-lattice coupling in multi-scale
modeling

Markus Weißenhofer<EMAIL_ADDRESS>Department of Physics,
University of Konstanz, DE-78457 Konstanz, Germany  Hannah Lange  Department
of Chemistry/Phys. Chemistry, LMU Munich, Butenandtstrasse 11, D-81377 Munich,
Germany  Akashdeep Kamra  Condensed Matter Physics Center (IFIMAC) and
Departamento de Física Teórica de la Materia Condensada, Universidad Autónoma
de Madrid, E-28049 Madrid, Spain  Sergiy Mankovsky  Department of
Chemistry/Phys. Chemistry, LMU Munich, Butenandtstrasse 11, D-81377 Munich,
Germany  Svitlana Polesya  Department of Chemistry/Phys. Chemistry, LMU
Munich, Butenandtstrasse 11, D-81377 Munich, Germany  Hubert Ebert  Department
of Chemistry/Phys. Chemistry, LMU Munich, Butenandtstrasse 11, D-81377 Munich,
Germany  Ulrich Nowak  Department of Physics, University of Konstanz, DE-78457
Konstanz, Germany

(September 3, 2024)

###### Abstract

In the spirit of multi-scale modeling, we  develop a theoretical framework for
spin-lattice coupling that connects, on the one hand, to ab initio
calculations of spin-lattice coupling parameters and, on the other hand, to
the magneto-elastic continuum theory.  The derived Hamiltonian describes a
closed system of spin and lattice degrees of freedom  and explicitly conserves
the total momentum, angular momentum and energy. Using a new numerical
implementation that corrects earlier Suzuki-Trotter decompositions we perform
simulations on the basis of the resulting equations of motion to investigate
the combined magnetic and mechanical motion of a ferromagnetic nanoparticle,
thereby validating our developed method. In addition to the ferromagnetic
resonance mode of the spin system we find another low-frequency mechanical
response and a rotation of the particle according to the Einstein-de-Haas
effect.  The framework developed herein will enable the use of multi-scale
modeling for investigating and understanding a broad range of magneto-
mechanical phenomena from slow to ultrafast time scales.

SD

    spin dynamics
SLD

    spin-lattice dynamics
FMR

    ferromagnetic resonance
DMI

    Dzyaloshinskii-Moriya interaction
ME

    magneto-elastic
angular momentum

    angular momentum
SLC

    spin-lattice coupling
SSC

    spin-spin coupling
EdH

    Einstein-de-Haas

The spin-orbit interaction is a relativistic  effect at the heart of modern
spintronics Manchon  _et al._ (2015). It couples the electron’s spin to its
orbital motion and  plays a central role in quantum materials  bearing high
potential for future nanoelectronic devices.  Its manifestations include
phenomena like magneto-crystalline anisotropy and Dzyaloshinskii-Moriya
interaction (DMI) Dzyaloshinsky (1958); Moriya (1960). While the focus in
spintronics has long been on  electrons or magnons as  carriers of angular
momentum, newer lines of research include even circularly polarized phonons to
fully understand and control the flow of angular momentum in a material
Hirohata  _et al._ (2020); Garanin and Chudnovsky (2015); Rückriegel  _et al._
(2020); Mentink  _et al._ (2019). Recently, it was demonstrated that even on
ultrashort time scales angular momentum can be transferred from the spin
system to the lattice Tauchert  _et al._ (2022). In the lattice, the spin
angular momentum is absorbed  by phonons carrying angular momentum till — on
larger times scales — the macroscopic Einstein-de-Haas (EdH) effect sets in
Dornes  _et al._ (2019). A coupling between spin and lattice degrees of
freedom that - beside the exchange of energy - includes the exchange of
angular momentum must be based on spin-orbit coupling, the effect  of which
has to be taken into account for a complete description of spin-lattice
dynamics (SLD).

Descriptions of spintronic phenomena are often based on spin models, which
treat the lattice degrees of freedom as a heat bath and define the spin
Hamiltonian with its magneto-crystalline anisotropy and DMI for a rigid
lattice. Consequently, the spin angular momentum is not conserved. First
attempts, to develop a framework for the calculation of coupled SLD — also
referred to as molecular and spin dynamics simulations — suffer from an
incomplete formulation of the spin lattice interaction Ma and Woo (2009);
Perera  _et al._ (2016); Strungaru  _et al._ (2021); Aßmann and Nowak (2019).
The works by Aßmann Aßmann and Nowak (2019) and Strungaru Strungaru  _et al._
(2021) assume a pseudo dipolar coupling that conserves the total angular
momentum \- a prerequisite for a well-defined spin-lattice coupling (SLC).
However, it suffers from the fact that it is not linked to first principles
calculations of SLC terms, which always rest on an expansion of the spin
Hamiltonian with respect to small distortions of the lattice. The works by
Hellsvik et al. Hellsvik  _et al._ (2019), Sadhukan et al. Sadhukhan  _et al._
(2022) and Mankovsky et al. Mankovsky  _et al._ (2022), use exactly these SLC
terms, that can be derived from first principles. However, these terms do not
conserve the total angular momentum since they are not rotationally invariant.
This inconsistency and the need for rotational invariance has already been
pointed out 50 years ago in the context of the magneto-elastic (ME) theory
Melcher (1970, 1972), a continuum theory that  approximates a microscopic spin
lattice model Hamiltonian on larger length scales.

In this  Letter, we  develop a rotationally invariant description of spin
lattice interaction for multi-scale modeling relaxing the assumption of a
rigid lattice with fixed orientation. The resulting Hamiltonian is
translationally and rotationally invariant, keeping total energy, momentum,
and angular momentum constant. All terms can be linked to the recently
developed ab initio methods that allow for a first principles calculation of
model parameters Hellsvik  _et al._ (2019); Sadhukhan  _et al._ (2022);
Mankovsky  _et al._ (2022), opening perspectives for multi-scale modeling of
SLD. We also demonstrate that our spin-lattice Hamiltonian represents the
discrete formulation of magneto-elastic theory and we link the microscopic
parameters with the magneto-elastic constants. We show that even terms that -
in a spin model - do not include any lattice distortions must transfer angular
momentum to the lattice. Furthermore, we derive the equations of motion for
spin and lattice degrees of freedom and solve them numerically with an
appropriate Suzuki-Trotter decomposition. Finlly, we present first simulations
of the precession of a magnetized body and spin dynamics including the
resulting response of the lattice.

A complete Hamiltonian accounting simultaneously for the spin and lattice
degrees of freedom contains contributions from the lattice degrees of freedom
(kinetic energy and pair potentials) as well as contributions which include
the spin degrees of freedom. The latter can be expressed as an expansion of
relativistic spin-spin interactions for small distortions Hellsvik  _et al._
(2019); Mankovsky  _et al._ (2022),

$\displaystyle\mathcal{H}_{\mathrm{SLC}}\approx\sum_{ij,\alpha\beta}\mathcal{J}^{\alpha\beta}_{ij}S_{i}^{\alpha}S_{j}^{\beta}+\sum_{ijk,\alpha\beta\mu}\mathcal{J}^{\alpha\beta,\mu}_{ij,k}S_{i}^{\alpha}S_{j}^{\beta}\left(u_{k}^{\,\mu}-u_{i}^{\,\mu}\right)$
$\displaystyle+\sum_{ijkl,\alpha\beta\mu\nu}\mathcal{J}^{\alpha\beta,\mu\nu}_{ij,kl}S_{i}^{\alpha}S_{j}^{\beta}\left(u_{k}^{\,\mu}-u_{i}^{\,\mu}\right)\left(u_{l}^{\,\nu}-u_{i}^{\,\nu}\right),$
(1)

where the summation  runs over the lattice (latin indices) and Cartesian
coordinates (greek indices)  111In this notation Eq. (1) also includes on-site
terms (where $i=j$).. $\bm{S}_{i}$ are  unit vectors representing the
direction of magnetic moments at sites $i$, and
$\bm{u}_{i}=\bm{r}_{i}-\bm{R}_{i}$ are displacement vectors of atoms $i$ at
position $\bm{r}_{i}$ (and equilibrium position $\bm{R}_{i}$ in a reference
configuration, see  Fig. 1). The spin-spin coupling (SSC)
$\mathcal{J}^{\alpha\beta}_{ij}$ and SLC tensors
$\mathcal{J}^{\alpha\beta,\mu}_{ij,k}=\partial\mathcal{J}^{\alpha\beta}_{ij}/\partial
u_{k}^{\mu}$ are defined with respect to a chosen coordinate system. As shown
by Mankovsky et al., these tensors can be calculated quantitatively from first
principles Mankovsky  _et al._ (2022).

The relative displacements $\left(u_{k}^{\,\mu}-u_{i}^{\,\mu}\right)$ w.r.t. a
reference atom $i$ take into account deformations of the lattice. 222Note that
in the work of Mankovsky et al. Mankovsky  _et al._ (2022) these tensors are
calculated taking only one displacement at site $k$ into account while all the
other atoms are in their equilibrium position. In this case the relative
displacement is $u_{k}^{\,\mu}-u_{i}^{\,\mu}=u_{k}^{\,\mu}$. However, for a
system that is displaced as a whole (equally for all sites $i$) the relative
displacement will vanish and there is no additional contribution to the
potential energy. They are the discrete lattice representation of the strain
and rotation tensor elements of elasticity theory. As such, Eq. (1) represents
the discrete formulation of  the ME theory Kittel (1949), from which we can
derive an extended expression for the ME energy density,

$\displaystyle\epsilon$
$\displaystyle=\sum_{\alpha\beta\mu\nu}S^{\alpha}S^{\beta}\Big{(}{B}^{\mathrm{s}}_{\alpha\beta,\mu\nu}\varepsilon_{\mu\nu}+{B}^{\mathrm{as}}_{\alpha\beta,\mu\nu}\omega_{\mu\nu}\Big{)}$
$\displaystyle+\sum_{\alpha\beta\gamma\mu\nu}\partial_{\beta}S^{\alpha}\partial_{\gamma}S^{\alpha}\Big{(}A^{\mathrm{s}}_{\beta\gamma,\mu\nu}\varepsilon_{\mu\nu}+A^{\mathrm{as}}_{\beta\gamma,\mu\nu}\omega_{\mu\nu}\Big{)}$
(2)
$\displaystyle+\sum_{\alpha\beta\gamma\delta\mu\nu}\varepsilon^{\alpha\beta\gamma}S^{\alpha}\partial_{\delta}S^{\beta}\Big{(}D^{\mathrm{s}}_{\gamma\delta,\mu\nu}\varepsilon_{\mu\nu}+D^{\mathrm{as}}_{\gamma\delta,\mu\nu}\omega_{\mu\nu}\Big{)},$

where $\bm{S}$ is the continuous magnetization, $\varepsilon_{\mu\nu}$ the
strain tensor, $\varepsilon^{\alpha\beta\gamma}$ the Levi-Civita symbol, and
$\omega_{\mu\nu}$ the rotation tensor. The  important role of the latter  in
the ME theory has been addressed before by Melcher Melcher (1970, 1972),  and
reaffirmed in recent experiments Xu  _et al._ (2020); Küß  _et al._ (2020,
2022). The terms in Eq. (2) model anisotropy, Heisenberg exchange, and DMI due
to lattice distortions and  the corresponding symmetric/antisymmetric ME
tensors ${B}^{\mathrm{s/as}}_{\alpha\beta,\mu\nu}$,
$A^{\mathrm{s/as}}_{\beta\gamma,\mu\nu}$, and
$D^{\mathrm{s/as}}_{\gamma\delta,\mu\nu}$ can be obtained from the microscopic
SLC tensors.  A detailed derivation can be found in the Supplemental Material
Sup  as well as the connection between the ME constants Kittel (1949) and the
microscopic SLC tensors.

Looking at Eq. (1) one finds immediately that this Hamiltonian does not
conserve the total (spin and lattice) angular momentum, since it is not
rotationally invariant. It is, hence, not capable of describing spins  plus
lattice as a closed system.  To understand this, we examine an isotropic
Heisenberg model with a uniaxial on-site anisotropy,

$\displaystyle\mathcal{H}_{ani}=-\sum_{ij}J_{ij}\bm{S}_{i}\cdot\bm{S}_{j}-d_{z}\sum_{i}\left(S_{i}^{z}\right)^{2},$
(3)

for a system with the $z$ axis being the easy axis of the magnetization.
Here, the Heisenberg exchange interaction term is rotationally invariant and
conserves the total spin angular momentum. The anisotropy term, however, is
not rotationally invariant and the total spin angular momentum is, hence, not
conserved. To keep the total angular momentum conserved, the spin angular
momentum would have to go to the lattice but since this term does not include
any lattice degrees of freedom it cannot.

Figure 1: Rotation and translation of a magnetized sample. The reference
configuration at $t=0$ is denoted by $\bm{R}_{i}$ (left), for $t>0$ by
$\bm{r}_{i}$ (right). During its motion, the easy axis for an atom $i$ at
position $\bm{r}_{i}$ or $\bm{R}_{i}$ can be defined via its upper and lower
neighbors at position ${r}_{i}^{z\pm}$.

The situation is shown in Fig. 1. Let us assume the atoms are at time $t=0$ in
equilibrium positions $\bm{R}_{i}$ in a reference  configuration with the easy
axis along $z$ (defined in the lab frame). When the sample starts moving the
lattice the atom positions at later times are
$\bm{r}_{i}(t)=\bm{R}_{i}+\bm{u}_{i}(t)$ and the easy axis may no longer be
aligned with the $z$-axis of the  lab frame. Consequently, the anisotropy term
in Eq. (3) has to be transformed. This can be done by projecting the spin
orientation using unit vectors that are  defined in terms of the respective
neighbor atoms.

In a cubic lattice an appropriate unit vector $\bm{e}^{z}$ can be defined via
the upper ($+$) and lower ($-$) nearest neighbors of atom $i$ at positions
$\bm{r}_{i}^{z\pm}=\bm{R_{i}}^{z\pm}+\bm{u_{i}}^{z\pm}$, yielding

$\bm{e}_{i}^{z(\pm)}=\frac{\bm{r_{i}}^{z\pm}-\bm{r}_{i}}{|\bm{r}_{i}^{z\pm}-\bm{r}_{i}|}\,.$
(4)

Now we can  write the anisotropy term of Eq. (3)  without reference to a
specific coordinate system and we obtain

$\displaystyle\mathcal{H}_{ani}$
$\displaystyle=-\sum_{ij}J_{ij}\bm{S}_{i}\cdot\bm{S}_{j}-\frac{d_{z}}{2}\sum_{i}\Big{[}\Big{(}\bm{S}_{i}\cdot\frac{\bm{r}_{i}^{z+}-\bm{r}_{i}}{|\bm{r}_{i}^{z+}-\bm{r}_{i}|}\Big{)}^{2}$
$\displaystyle+\Big{(}\bm{S}_{i}\cdot\frac{\bm{r}_{i}^{z-}-\bm{r}_{i}}{|\bm{r}_{i}^{z-}-\bm{r}_{i}|}\Big{)}^{2}\Big{]}.$
(5)

The resulting Hamiltonian contains only scalar products of the spins
$\bm{S}_{i}$ and differences of position vectors $\bm{r}_{i}$. It is hence
translationally and rotationally invariant and will keep the total momentum
and angular momentum constant. Most  importantly, the transformed Hamiltonian
contains the lattice degrees of freedom  explicitly, even though the original
one did not. Only this makes it possible to transfer angular momentum from the
spins to the lattice, keeping the total angular momentum constant.

The  microscopic origin of the  anisotropy can be crucial in determining the
exact form of the definition of the $z$ direction. It is however sufficient to
take only two neighbors into account to capture the majority of effects. Note
that this local definition of an easy axis does not only work for a global
rotation but also for deformations of the sample. Furthermore, both, upper and
lower neighbors, are used to define the local easy axis for spins $i$, a
definition that  holds also at surfaces with a reduced number of neighbors.
Nevertheless, it should be stressed that this definition is neither unique nor
trivial, since the choice of neighbors for the definition of the unit vectors
will affect the equations of motion and  the atoms the angular momentum is
transferred to.

The transformation above can be extended to other contributions of the spin-
lattice Hamiltonian (1). In order to do so, the local definition for the unit
vector in $z$ direction for atom $i$ from Eqs. (4) can be generalized to a set
of three orthogonal directions $\alpha$,

$\displaystyle\bm{e}_{i}^{\alpha(\pm)}=\frac{\bm{r_{i}}^{\alpha\pm}-\bm{r}_{i}^{\alpha}}{|\bm{r}_{i}^{\alpha\pm}-\bm{r}_{i}^{\alpha}|}.$
(6)

Similar to the case of a uniaxial anisotropy, these unit vectors can be used
to transform the first term of Eq. (1),

$\displaystyle\mathcal{H}_{\mathrm{SS}}$
$\displaystyle=\sum_{ij}\sum_{\alpha\beta}\mathcal{J}^{\alpha\beta}_{ij}S_{i}^{\alpha}S_{j}^{\beta}=\sum_{ij}\sum_{\alpha\beta}\mathcal{J}^{\alpha\beta}_{ij}(\bm{S}_{i}\cdot\bm{e}_{i}^{\,\alpha})(\bm{S}_{j}\cdot\bm{e}_{j}^{\,\beta})$
(7)

Again, this Hamiltonian consists of  scalar products of  spins and differences
of position vectors and is hence rotationally invariant. Analogously, the
spin-lattice Hamiltonian becomes

$\displaystyle\mathcal{H}_{\mathrm{SLC}}$
$\displaystyle=\mathcal{H}_{\mathrm{SS}}+\sum_{ijk,\alpha\beta\mu}\mathcal{J}^{\alpha\beta,\mu}_{ij,k}(\bm{S}_{i}\cdot\bm{e}_{i}^{\,\alpha})(\bm{S}_{j}\cdot\bm{e}_{j}^{\,\beta})\times$
$\displaystyle\times\big{[}\left(\bm{r}_{k}-\bm{r}_{i}\right)\cdot\bm{e}_{k}^{\,\mu}-R_{ki}\big{]}+\dots\,,$
(8)

where $R_{ki}$ is the equilibrium distance between atoms $k$ and $i$ in the
reference  configuration. This Hamiltonian consists of the rotationally
invariant spin-spin term ($\mathcal{H}_{\mathrm{SS}}$) and a second spin-spin-
lattice term.

Each term in the initial formulation (1) that breaks rotational symmetry now
depends on the spins and the lattice positions and, hence, can transfer
angular momentum between the two subsystems. Thus, the dominating terms for
angular momentum transfer may vary for different materials, depending on the
specific values of the SSC and SLC tensors.  For example, in Fe the transfer
is mainly via the spin-lattice DMI Mankovsky  _et al._ (2022), whereas in FePt
the dominating terms are two-site anisotropy terms (see Table 1).

Material |  Contribution to SSC | $|J_{ij}^{\mathrm{iso}}|$ | $|J_{ij}^{\mathrm{dia-a}}|$ | $|J_{ij}^{\mathrm{off-s}}|$ | $|\bm{D}_{ij}|$
---|---|---|---|---|---
Fe | $J_{ij}^{SSC}$ | 11.389 | 0.019 | 0.017 | 0.0
$\Delta J_{ij}^{SLC}(u_{j}^{x})$ | 1.587 | 0.002 | 0.003 | 0.062
FePt | $J_{ij}^{SSC}$ | 9.590 | 0.320 | 0.209 | 0.0
$\Delta J_{ij}^{SLC}(u_{j}^{x})$ | 1.960 | 0.023 | 0.024 | 0.089
Table 1: Maximal absolute SSC $J_{ij}^{SSC}$ and maximal modification of SSC
due to SLC $\Delta J_{ij}^{SLC}=J_{ij,j}\cdot u_{j}^{x}$ (in meV) in the
presence of a displacement $u_{j}^{x}=0.03\,a_{\mathrm{lat}}$ in Fe and FePt
for different SSC contributions. For both materials, we consider sites $i$ and
$j$ being Fe atoms with different distances $r_{ij}$ and list the values for
the pair $ij$ with the largest contribution to the respective parts of the SSC
tensor. In Fe, the largest contribution which can transfer angular momentum is
the spin-lattice DMI $|\Delta\bm{D}_{ij}|=|\bm{D}_{ij,j}^{x}\cdot u_{j}^{x}|$
for $r_{ij}=1\,a_{\mathrm{lat}}$, in FePt it is the spin-spin anti-symmetric
diagonal part $|J_{ij}^{\mathrm{dia-a}}|=\frac{1}{2}|J_{ij}^{xx}-J_{ij}^{zz}|$
for $r_{ij}=1.414\,a_{\mathrm{lat}}$. The symmetrized off-diagonal elements
are defined as
$J_{ij}^{\mathrm{off-s}}=\frac{1}{2}\left(J_{ij}^{xy}+J_{ij}^{yx}\right)$. For
details see the Supplemental  Material Sup .

As a first application of our formulation and to test its validity, we perform
combined SLD simulations using the following Hamiltonian for a simple cubic
lattice,

$\displaystyle\mathcal{H}=$
$\displaystyle\mathcal{H}_{ani}+\sum_{i}\frac{\bm{p}_{i}^{2}}{2m}+V_{0}\sum_{ij}\frac{(r_{ij}-R_{ij})^{2}}{R_{ij}}$
(9)

that extends our rotationally invariant formulation of the spin Hamiltonian of
Eq. (5) by terms describing the interaction and the kinetic energy of the
lattice, with $m$ being the mass of the atoms and $V_{0}$ describing the
strength of the lattice interactions in  the harmonic approximation. For the
sake of simplicity, we assume that these interactions are restricted to the
first three shells of neighbors and that they scale inversely with the
equilibrium distance.

Evaluating the dynamics of spin and lattice degrees of freedom
$\\{\bm{r}_{i},\bm{p}_{i},\bm{S}_{i}\\}$ requires the concurrent solution of
the coupled equations of motion,

$\displaystyle\dot{\bm{r}}_{i}=\frac{\partial\mathcal{H}}{\partial\bm{p}_{i}},\hskip
6.60004pt\dot{\bm{p}}_{i}=-\frac{\partial\mathcal{H}}{\partial\bm{r_{i}}},\hskip
6.60004pt\mathrm{and}\hskip
6.60004pt\dot{\bm{S}}_{i}=\frac{\gamma}{\mu_{\mathrm{s}}}\bm{S}_{i}\times\frac{\partial\mathcal{H}}{\partial\bm{S}_{i}},$
(10)

with $\gamma$ and $\mu_{\mathrm{s}}$ being the absolute values of the
gyromagnetic ratio and the magnetic spin moment, respectively. Conservation of
energy, momentum and angular momentum can be ensured by using a symplectic
algorithm. Here, we use a scheme based on the Liouville formalism Frenkel and
Smit (2001) and the Suzuki-Trotter decomposition Suzuki (1976) that was
initially proposed in Ref. Omelyan  _et al._ (2001) and has proven reliable
for the simulation of combined SLD Tsai  _et al._ (2005); Ma and Woo (2009);
Ma  _et al._ (2016); Aßmann and Nowak (2019); Strungaru  _et al._ (2021). Note
that the presence of a uniaxial on-site anisotropy term, which is quadratic in
the spins, requires a further decomposition of the integration scheme that has
not been discussed in literature so far. Details, tests of the conservation of
the total angular momentum and the energy of the system and a comparison of
the temperature dependence of the magnetization with spin dynamics (SD)
simulations based on the stochastic Landau-Lifshitz-Gilbert equation of motion
can be found in the Supplemental  Material Sup .

As a first application we study the coupled magnetization and lattice dynamics
of a free cubic nanoparticle. For this simulation we assume that initially the
cube is oriented such that the easy axis is aligned with the $z$-axis and all
spins point along $\bm{m}_{0}=(0.1,0,\sqrt{1-0.1^{2}})^{\mathrm{T}}$. This
gives rise to a coherent precession of the magnetization along with mechanical
motion of the cube. Fig. 2 displays the time evolution of the magnetization
$\bm{m}=\frac{1}{N}\sum_{i}^{N}\bm{S}_{i}$ and the Fourier transform of its
$y$-component for a nanoparticle consisting of $4^{3}$ atoms. The light curves
are obtained by pure SD simulations, for which we keep the position of the
atoms  fixed. In contrast to pure SD, the SLD simulations produce oscillations
at two characteristic frequencies $\omega_{\bm{n}_{3}}\approx 3.89\times
10^{-3}\gamma J/\mu_{\mathrm{s}}$ and $\omega_{\mathrm{FMR}}\approx
0.152\gamma J/\mu_{\mathrm{s}}$. The peak at $\omega_{\mathrm{FMR}}$ can be
attributed to the usual ferromagnetic resonance (FMR) frequency and is close
to the value predicted by linear spin wave theory
$\omega_{\mathrm{FMR}}=2\bar{d}\gamma/\mu_{\mathrm{s}}=0.15\gamma
J/\mu_{\mathrm{s}}$, where $\bar{d}$ is the averaged uniaxial magnetic
anisotropy 333Given the form of the anisotropy in Eq. (5), the effective
uniaxial anisotropy of spins at two faces of the cube is reduced by a factor
of two. Thus, for a cube consisting of $4^{3}$ atoms, we calculate
$\bar{d}=0.75d_{z}$.. The SD value ($\omega_{\mathrm{FMR}}\approx 0.148\gamma
J/\mu_{\mathrm{s}}$) is slightly smaller due to finite size effects: the spins
at the edges of the cube lag behind, since their anisotropy field is weaker
[cf. Eq. (5)], slowing down the overall precession frequency. Surprisingly,
the FMR frequency of SLD is shifted to higher values as compared to the SD
value.

The emergence of the peak at $\omega_{\bm{n}_{3}}$ as well as the shift of the
FMR frequency are a result of the mechanical motion of the nanoparticle, which
can be characterized  by the vectors $\bm{n}_{\alpha}$ ($\alpha\in{1,2,3}$)
normal to the faces of the cube.  As the easy axis is firmly attached to one
of these vectors (without loss of generality, $\bm{n}_{3}$), we can compare
the dynamics of the easy axis and the magnetization, see Fig. 3a. We find that
the magnetization precesses around the easy axis, which itself revolves around
$\bm{m}_{0}$ exactly at $\omega_{\bm{n}_{3}}$, giving rise to the second peak
in Fig. 2 and the shift of the FMR frequency.

Figure 2: Coherent magnetization dynamics of a free cubic nanoparticle
obtained from SLD in comparison with SD simulations (light curves). Top:
magnetization vector components $m_{\alpha}$ versus time. Bottom: Fourier
transform of $m_{y}$.

This emergence of the easy-axis-precession was predicted in Ref. Usov and
Liubimov (2015) based on a simple rigid-body-macrospin model Usov and Liubimov
(2012); Usadel and Usadel (2015); Keshtgar  _et al._ (2017). There, the
magnetic nanoparticle is described by the normalized magnetization $\bm{m}$
and the vectors $\bm{n}_{\alpha}$ introduced above. The dynamics of these
vectors $\bm{n}_{\alpha}$ are given by
$\dot{\bm{n}}_{\alpha}=\bm{\omega}\times\bm{n}_{\alpha}$, where $\bm{\omega}$
is the angular velocity of the nanoparticle in the laboratory frame, which is
related to its angular momentum via $\bm{L}=\Theta\bm{\omega}$. For a cubic
nanoparticle, the moment of inertia is given by $\Theta=\frac{1}{6}Nml^{2}$,
$l$ being the cube size and $N$ being the number of atoms. Conservation of
angular momentum requires that
$-N\frac{\mu_{\mathrm{s}}}{\gamma}\dot{\bm{m}}+\dot{\bm{L}}=0$. If the
nanoparticle is initially at rest, we get
$\bm{\omega}(t)=\frac{6\mu_{\mathrm{s}}}{ml^{2}\gamma}(\bm{m}(t)-\bm{m}_{0})$,
with $\bm{m}_{0}$ being the initial orientation of the magnetization. This
yields
$\dot{\bm{n}}_{\alpha}=\frac{6\mu_{\mathrm{s}}}{ml^{2}\gamma}(\bm{m}(t)\times\bm{n}_{\alpha}-\bm{m}_{0}\times\bm{n}_{\alpha})$.
The first contribution to the torque depends on the current value of the
magnetization and amounts to zero for $\bm{n}_{3}$, due to the rapid
oscillations of $\bm{m}(t)$ around $\bm{n}_{3}$. The second term describes a
simple rotation of $\bm{n}_{3}$ around the initial direction of the
magnetization $\bm{m}_{0}$ with frequency
$\omega_{\bm{n}_{3}}=\frac{6\mu_{\mathrm{s}}}{ml^{2}\gamma}$. For the
parameters used here, we calculate $\omega_{\bm{n}_{3}}\approx 3.75\times
10^{-3}\gamma J/\mu_{\mathrm{s}}$, which is in close agreement with the
simulation results.

In addition to the precession of $\bm{n}_{3}$, we find an EdH-type rotation of
$\bm{n}_{1,2}$ with $\omega^{\mathrm{EdH}}\approx 1.81\times 10^{-5}\gamma
J/\mu_{\mathrm{s}}$ around the $z$-axis of the laboratory frame (see Fig. 3b).
This rotation occurs when the average magnetization differs from its initial
value, since the average angular velocity of the cubic nanoparticle is given
by
$\langle\bm{\omega}(t)\rangle=\frac{6\mu_{\mathrm{s}}}{ml^{2}\gamma}\langle(\bm{m}(t)-\bm{m}_{0})\rangle$.
This leads to a nonzero value only for the $z$-component of the angular
velocity (cf. Fig. 2).

Fig. 3c displays all three characteristic frequencies versus cube size. The
FMR frequency is the highest and approaches the bulk value for large cubes.
The other two (mechanical) frequencies scale with $l^{-2}$ and have a constant
ratio over the range considered here.  Besides testing the validity of the
analytical expression for $\omega_{\bm{n}_{3}}$, this allows us to estimate
the mechanical frequencies for real materials. E.g., for an FePt nano-cube
with edge length of $100\text{\,}\mathrm{nm}$ we get
$\omega_{\bm{n}_{3}}\approx$250\text{\,}\mathrm{kHz}$$ and
$\omega_{\mathrm{EdH}}\approx$1\text{\,}\mathrm{kHz}$$ using
$\mu_{\mathrm{s}}\approx 3.23\mu_{\mathrm{B}}$,
$\gamma=$1.76\text{\times}{10}^{11}\text{\,}\mathrm{s}^{-1}\mathrm{T}^{-1}$$
and
$m=m_{\mathrm{Fe}}+m_{\mathrm{Pt}}\approx$4.17\text{\times}{10}^{-25}\text{\,}\mathrm{kg}$$
(the FMR frequency for a highly coercive material such as FePt is of the order
of $100\text{\,}\mathrm{GHz}$) Barker  _et al._ (2010).

Figure 3: Mechanical motion of a cubic nanoparticle excited by coherent
precession of the magnetization around $\bm{n}_{3}$ (easy axis). a) Spiraling
magnetization dynamics and easy axis precession in the time interval
$[0,10^{3}]\mu_{\mathrm{s}}/\gamma J$. b) Sketch of the two characteristic
mechanical modes. c) Characteristic frequencies versus cube size $l$. Dotted
lines correspond to theory curves as explained in the text.

In summary, we have developed a rotationally invariant formulation of coupled
spin-lattice dynamics for multi-scale modeling of magneto-mechanical motion.
It successfully integrates first principles evaluation of SLC parameters, ME
continuum theory, and spin lattice dynamics simulations. Employing our
developed framework and a new numerical implementation that corrects earlier
Suzuki-Trotter decompositions we simulate combined magneto-mechanical dynamics
of a ferromagnetic nanoparticle, thereby validating our formulation. Our
simulations demonstrate that in addition to the ferromagnetic resonance mode
of the spin system there are two low-frequency mechanical modes describing the
precession of the easy axis and a rotation of the particle according to the
EdH effect. By incorporating total angular momentum conservation, our work
provides the tools for simulation of a broad range of magneto-mechanical
phenomena. Therefore it is crucial to the understanding of recent and ongoing
intriguing experiments, e.g. on magnon-phonon coupling or ultrafast
magnetization dynamics.

_Acknowledgements.–_ Work in Konstanz is supported by the Deutsche
Forschungsgemeinschaft (DFG) via SFB 1432 and Project No. NO 290/5-1. A.K.
acknowledges financial support from the Spanish Ministry for Science and
Innovation – AEI Grant CEX2018-000805-M (through the “Maria de Maeztu”
Programme for Units of Excellence in R&D). Work in Munich is supported by the
DFG via SFB 1277.

## References

  * Manchon  _et al._ (2015) A. Manchon, H. C. Koo, J. Nitta, S. M. Frolov, and R. A. Duine, Nature Materials 14, 871 (2015).
  * Dzyaloshinsky (1958) I. Dzyaloshinsky, Journal of Physics and Chemistry of Solids 4, 241 (1958).
  * Moriya (1960) T. Moriya, Phys. Rev. 120, 91 (1960).
  * Hirohata  _et al._ (2020) A. Hirohata, K. Yamada, Y. Nakatani, I.-L. Prejbeanu, B. Diény, P. Pirro, and B. Hillebrands, Journal of Magnetism and Magnetic Materials 509, 166711 (2020).
  * Garanin and Chudnovsky (2015) D. A. Garanin and E. M. Chudnovsky, Phys. Rev. B 92, 024421 (2015).
  * Rückriegel  _et al._ (2020) A. Rückriegel, S. Streib, G. E. W. Bauer, and R. A. Duine, Phys. Rev. B 101, 104402 (2020).
  * Mentink  _et al._ (2019) J. H. Mentink, M. I. Katsnelson, and M. Lemeshko, Phys. Rev. B 99, 064428 (2019).
  * Tauchert  _et al._ (2022) S. R. Tauchert, M. Volkov, D. Ehberger, D. Kazenwadel, M. Evers, H. Lange, A. Donges, A. Book, W. Kreuzpaintner, U. Nowak, and P. Baum, Nature 602, 73 (2022).
  * Dornes  _et al._ (2019) C. Dornes, Y. Acremann, M. Savoini, M. Kubli, M. J. Neugebauer, E. Abreu, L. Huber, G. Lantz, C. A. F. Vaz, H. Lemke, E. M. Bothschafter, M. Porer, V. Esposito, L. Rettig, M. Buzzi, A. Alberca, Y. W. Windsor, P. Beaud, U. Staub, D. Zhu, S. Song, J. M. Glownia, and S. L. Johnson, Nature 565, 209 (2019).
  * Ma and Woo (2009) P.-W. Ma and C. H. Woo, Phys. Rev. E 79, 046703 (2009).
  * Perera  _et al._ (2016) D. Perera, M. Eisenbach, D. M. Nicholson, G. M. Stocks, and D. P. Landau, Phys. Rev. B 93, 060402(R) (2016).
  * Strungaru  _et al._ (2021) M. Strungaru, M. O. A. Ellis, S. Ruta, O. Chubykalo-Fesenko, R. F. L. Evans, and R. W. Chantrell, Phys. Rev. B 103, 024429 (2021).
  * Aßmann and Nowak (2019) M. Aßmann and U. Nowak, Journal of Magnetism and Magnetic Materials 469, 217 (2019).
  * Hellsvik  _et al._ (2019) J. Hellsvik, D. Thonig, K. Modin, D. Iuşan, A. Bergman, O. Eriksson, L. Bergqvist, and A. Delin, Phys. Rev. B 99, 104302 (2019).
  * Sadhukhan  _et al._ (2022) B. Sadhukhan, A. Bergman, Y. O. Kvashnin, J. Hellsvik, and A. Delin, Phys. Rev. B 105, 104418 (2022).
  * Mankovsky  _et al._ (2022) S. Mankovsky, S. Polesya, H. Lange, M. Weißenhofer, U. Nowak, and H. Ebert, Phys. Rev. Lett. 129, 067202 (2022).
  * Melcher (1970) R. L. Melcher, Phys. Rev. Lett. 25, 1201 (1970).
  * Melcher (1972) R. L. Melcher, Phys. Rev. Lett. 28, 165 (1972).
  * Note (1) In this notation Eq. (1) also includes on-site terms (where $i=j$).
  * Note (2) Note that in the work of Mankovsky et al. Mankovsky  _et al._ (2022) these tensors are calculated taking only one displacement at site $k$ into account while all the other atoms are in their equilibrium position. In this case the relative displacement is $u_{k}^{\,\mu}-u_{i}^{\,\mu}=u_{k}^{\,\mu}$. However, for a system that is displaced as a whole (equally for all sites $i$) the relative displacement will vanish and there is no additional contribution to the potential energy.
  * Kittel (1949) C. Kittel, Rev. Mod. Phys. 21, 541 (1949).
  * Xu  _et al._ (2020) M. Xu, K. Yamamoto, J. Puebla, K. Baumgaertl, B. Rana, K. Miura, H. Takahashi, D. Grundler, S. Maekawa, and Y. Otani, Science Advances 6, eabb1724 (2020).
  * Küß  _et al._ (2020) M. Küß, M. Heigl, L. Flacke, A. Hörner, M. Weiler, M. Albrecht, and A. Wixforth, Phys. Rev. Lett. 125, 217203 (2020).
  * Küß  _et al._ (2022) M. Küß, M. Albrecht, and M. Weiler, Front. Phys. 10, 981257 (2022).
  * (25) See Supplemental Material for (i) a detailed derivation of the continuum magnetoelastic theory relating its parameters to those obtained from first principles calculations, (ii) details on the numerical simulation, and (iii) evaluation of the magnetoelastic tensors using first principles calculation.
  * Frenkel and Smit (2001) D. Frenkel and B. Smit,  _Understanding molecular simulation: from algorithms to applications_ , Vol. 1 (Elsevier, 2001).
  * Suzuki (1976) M. Suzuki, Communications in Mathematical Physics 51, 183 (1976).
  * Omelyan  _et al._ (2001) I. P. Omelyan, I. M. Mryglod, and R. Folk, Phys. Rev. Lett. 86, 898 (2001).
  * Tsai  _et al._ (2005) S.-H. Tsai, H. K. Lee, and D. P. Landau, American Journal of Physics 73, 615 (2005), https://doi.org/10.1119/1.1900096 .
  * Ma  _et al._ (2016) P.-W. Ma, S. Dudarev, and C. Woo, Computer Physics Communications 207, 350 (2016).
  * Usov and Liubimov (2015) N. Usov and B. Y. Liubimov, Journal of Magnetism and Magnetic Materials 385, 339 (2015).
  * Note (3) Given the form of the anisotropy in Eq. (5), the effective uniaxial anisotropy of spins at two faces of the cube is reduced by a factor of two. Thus, for a cube consisting of $4^{3}$ atoms, we calculate $\bar{d}=0.75d_{z}$.
  * Usov and Liubimov (2012) N. A. Usov and B. Y. Liubimov, Journal of Applied Physics 112, 023901 (2012), https://doi.org/10.1063/1.4737126 .
  * Usadel and Usadel (2015) K. D. Usadel and C. Usadel, Journal of Applied Physics 118, 234303 (2015), https://doi.org/10.1063/1.4937919 .
  * Keshtgar  _et al._ (2017) H. Keshtgar, S. Streib, A. Kamra, Y. M. Blanter, and G. E. W. Bauer, Phys. Rev. B 95, 134447 (2017).
  * Barker  _et al._ (2010) J. Barker, R. F. L. Evans, R. W. Chantrell, D. Hinzke, and U. Nowak, Applied Physics Letters 97, 192504 (2010), https://doi.org/10.1063/1.3515928 .

  *[DMI]: Dzyaloshinskii-Moriya interaction
  *[EdH]: Einstein-de-Haas
  *[SLD]: spin-lattice dynamics
  *[SLC]: spin-lattice coupling
  *[ME]: magneto-elastic
  *[SSC]: spin-spin coupling
  *[SD]: spin dynamics
  *[FMR]: ferromagnetic resonance