solutionSolutionsolutionfile solutionfile[MultipleRISwithCorrelation]

#  Coverage Probability of Distributed IRS Systems Under Spatially Correlated
Channels

Anastasios Papazafeiropoulos, Cunhua Pan, Ahmet Elbir, Pandelis Kourtessis,
Symeon Chatzinotas, John M. Senior A. Papazafeiropoulos is with the
Communications and Intelligent Systems Research Group, University of
Hertfordshire, Hatfield AL10 9AB, U. K., and with SnT at the University of
Luxembourg, Luxembourg. C. Pan is with the School of Electronic Engineering
and Computer Science at Queen Mary University of London, London E1 4NS, U.K.
A. Elbir is with the EE department of Duzce University, Duzce, Turkey. P.
Kourtessis and John M. Senior are with the Communications and Intelligent
Systems Research Group, University of Hertfordshire, Hatfield AL10 9AB, U. K.
S. Chatzinotas is with the SnT at the University of Luxembourg, Luxembourg.
E-mails<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

This paper suggests the use of multiple distributed intelligent reflecting
surfaces (IRSs) towards a smarter control of the propagation environment.
Notably, we also take into account the inevitable correlated Rayleigh fading
in IRS-assisted systems. In particular, in a single-input and single-output
(SISO) system, we consider and compare two insightful scenarios, namely, a
finite number of large IRSs and a large number of finite size IRSs to show
which implementation method is more advantageous. In this direction, we derive
the coverage probability in closed-form for both cases contingent on
statistical channel state information (CSI) by using the deterministic
equivalent (DE) analysis. Next, we obtain the optimal coverage probability.
Among others, numerical results reveal that the addition of more surfaces
outperforms the design scheme of adding more elements per surface. Moreover,
in the case of uncorrelated Rayleigh fading, statistical CSI-based IRS systems
do not allow the optimization of the coverage probability.

###### Index Terms:

Intelligent reflecting surface (IRS), coverage probability, deterministic
equivalents, beyond 5G networks.

##  I Introduction

The advancements on metasurfaces have enabled the development of intelligent
reflecting surface (IRS), being a planar array that includes a large number of
nearly passive reflecting elements [1, 2]. IRS provides a smart radio
environment by realizing reflecting beamforming through its elements, which
can introduce phase adjustments on the impinging wave with different
objectives such as an increase of coverage and avoidance of obstacles. Also,
its construction principles allow affordable and green transmission due to its
low-cost hardware and low energy consumption, respectively.

Many works have approached the concept of IRS from the wireless communication
point of view due to its appealing advantages to achieve various tasks by
adjusting the phase shifts of the reflecting surface elements, e.g., see [2,
3, 4, 5, 6, 7, 8, 9] and references therein. Among them, in [4], we observe a
maximization of the sum-rate with a transmit power constraint, in [3], authors
achieved maximization of the energy efficiency with signal-to-interference-
plus-noise ratio (SINR) constraints, and in [7], channel estimation, being an
interesting research area in IRS-assisted systems due to their special
characteristics, by using a deep learning approach was proposed.

In particular, in the communication-theoretic direction, the study of the
coverage probability in IRS-assisted systems has attracted significant
attention [8, 10, 11]. However, all previous works assumed only one IRS while
the performance of systems aided simultaneously by multiple IRSs, offering
extended advantages such as a more robust avoidance of obstacles and improved
coverage has not been investigated except [12, 13, 14].

In parallel, the assumption of independent Rayleigh fading, which is commonly
assumed for tractable performance analysis is unrealistic for IRS-assisted
systems [15]. Although several works have accounted for the IRS correlation by
acknowledging its importance, they relied on conventional correlation models
[6], which are not directly applicable in IRSs as mentioned in [15], where a
practical correlation model was suggested.

Against the above background, we present the only work providing the coverage
probability in closed-form for single-input and single-output (SISO) systems
assisted simultaneously by multiple IRSs while accounting for the inevitable
correlated Rayleigh fading requiring suitable mathematical manipulations. In
particular, we consider two insightful design cornerstones: i) A finite number
of large IRSs (the number of elements per IRS grows large) and ii) a large
number of IRSs with each IRS having finite dimensions. Hence, contrary to [8],
we establish the theoretical framework incorporating correlated fading into
the analysis to identify the realistic potentials of IRSs before their final
implementation and we also study the performance when the IRSs number becomes
large. Compared to [9], which also assumed distributed IRSs and correlated
fading, we focus on the coverage probability instead of the achievable rate
and we rely on a more realistic correlation model while we account for the
scenario of a large number of IRSs, which has not been addressed before.
Moreover, we provide a methodology to optimize the reflect beamforming matrix
based on statistical channel state information (CSI) that enables optimization
at every several coherence intervals instead of frequent optimization at every
coherence interval as in works relying on instantaneous CSI.

Notation: Vectors and matrices are denoted by boldface lower and upper case
symbols, respectively. The notations $(\cdot)^{\scriptscriptstyle\mathsf{T}}$,
$(\cdot)^{\scriptscriptstyle\mathsf{H}}$, and $\tr\\!\left({\cdot}\right)$
represent the transpose, Hermitian transpose, and trace operators,
respectively. The expectation operator is denoted by
$\mathbb{E}\left[\cdot\right]$ while
$\text{diag}\left(\right)\left({\mathbf{a}}\right)$ and
$\text{diag}\left(\right)\left({\mathbf{A}}\right)$ express diagonal matrices
with diagonal elements being the elements of vector ${\mathbf{a}}$ and the
diagonal elements of ${\mathbf{A}}$, respectively. Also, the notations
$\arg\left(\cdot\right)$ and $\mod(\cdot,\cdot)$ denote the argument function
and the modulus operation while $\lfloor\cdot\rfloor$ truncates the argument.
Given two infinite sequences $a_{n}$ and $b_{n}$, the relation $a_{n}\asymp
b_{n}$ is equivalent to
$a_{n}-b_{n}\xrightarrow[n\rightarrow\infty]{\mbox{a.s.}}0$. Finally,
${\mathbf{b}}\sim{\cal C}{\cal N}{({\mathbf{0}},\mathbf{\Sigma})}$ represents
a circularly symmetric complex Gaussian vector with zero mean and covariance
matrix $\mathbf{\Sigma}$.

##  II System Model

We consider the smart connectivity between a single-antenna transmitter (TX)
and a single-antenna receiver (RX) enabled by means of a set of $M$
independent IRSs uniformly distributed in the intermediate space. We assume
that each IRS, controlled through a perfect backhaul link by the transmitter,
consists of a two-dimensional rectangular grid of
$N=N_{\mathrm{H}}N_{\mathrm{V}}$ passive unit elements with $N_{\mathrm{H}}$
elements per row and $N_{\mathrm{V}}$ elements per column that can modify the
phase shifts of impinging waves. Severe blockage effects make any direct
channel unavailable. To focus on the impact of the multitude of IRSs, we rely
on the assumption of perfect CSI, and thus, the results act as upper bounds of
practical implementations.

Let a block-fading model with independent realizations across different
coherence blocks for the description of all channels. In particular, we assume
the existence of a direct link and $M$ cascaded channels. The former is
described by $h_{\mathrm{d}}\sim\mathcal{CN}\left(0,\beta_{\mathrm{d}}\right)$
with $\beta_{\mathrm{d}}$ expressing the path-loss. Regarding the IRS-assisted
links,
${\mathbf{h}}_{m,1}=\left[h_{mn,1},\ldots,h_{mN,1}\right]^{{\scriptscriptstyle\mathsf{T}}}\in\mathbb{C}^{N\times
1}$ expresses the channel fading vector between the TX and the $m$th IRS while
${\mathbf{h}}_{m,2}=\left[h_{mn,2},\ldots,h_{mN,2}\right]^{{\scriptscriptstyle\mathsf{T}}}\in\mathbb{C}^{N\times
1}$ corresponds to the link between the $m$th IRS and the RX. Contrary to
existing works, relying on independent Rayleigh and Rician fading models, we
consider correlated Rayleigh fading111The extension to correlated Rician
fading, having a LoS component, is the topic of future work.. Hence, by
accounting for both small-scale fading and path-loss, we have

$\displaystyle{\mathbf{h}}_{m,1}$
$\displaystyle\sim\mathcal{CN}\left({\mathbf{0}},\beta_{m,1}{\mathbf{R}}_{m,1}\right),$
(1) $\displaystyle{\mathbf{h}}_{m,2}$
$\displaystyle\sim\mathcal{CN}\left({\mathbf{0}},\beta_{m,2}{\mathbf{R}}_{m,2}\right),$
(2)

where $\beta_{m,1}$, $\beta_{m,2}$ describe the path-losses while
${\mathbf{R}}_{m,1}\in\mathbb{C}^{N\times N}$,
${\mathbf{R}}_{m,2}\in\mathbb{C}^{N\times N}$ describe the spatial covariance
matrices of the respective links222The path-losses and the covariance matrices
are assumed known by applying practical methods, e.g., see [16].. Herein, we
account for the correlation model proposed in [15] as suitable for IRSs under
the conditions of rectangular IRSs and isotropic Rayleigh fading. Let the size
of each IRS element be $d_{\mathrm{H}}\times d_{\mathrm{V}}$, where
$d_{\mathrm{V}}$ and $d_{\mathrm{H}}$ express its vertical height and its
horizontal width, respectively. Then, the $\left(i,j\right)$th element of the
correlation matrix ${\mathbf{R}}_{m,k},$ $k\in\\{1,2\\}$ is given by

$\displaystyle
r_{nm,k}=d_{\mathrm{H}}d_{\mathrm{V}}\mathrm{sinc}\left(2\|\mathbf{u}_{n,k}-\mathbf{u}_{m,k}\|/\lambda\right),$
(3)

where
$\mathbf{u}_{\epsilon,k}=[0,\mod(\epsilon-1,N_{\mathrm{H}})d_{\mathrm{H}},\lfloor(\epsilon-1)/N_{\mathrm{H}}\rfloor
d_{\mathrm{V}}]^{\scriptscriptstyle\mathsf{T}}$, $\epsilon\in\\{i,j\\}$ and
$\lambda$ is the wavelength of the plane wave.

Based on a slowly varying flat-fading channel model, the complex-valued
received signal at the RX through the network of $M$ IRSs is described by

$\displaystyle
y=\bigg{(}\sum^{M}_{m=1}{\mathbf{h}}_{m,1}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Phi}}_{m}{\mathbf{h}}_{m,2}+h_{\mathrm{d}}\\!\bigg{)}x+n,$
(4)

where $n\sim\mathcal{CN}\left(0,N_{0}\right)$ is the additive white Gaussian
noise (AWGN) sample and $x$ is the transmitted data symbol satisfying
$\mathbb{E}[|x|^{2}]=P$, i.e., $P$ denotes the average power of the symbol.
The diagonal matrix
${\bm{\Phi}}_{m}=\mathrm{diag}\left(\alpha_{m1}\exp\left(j\theta_{m1}\right),\ldots,\alpha_{mN}\exp\left(j\theta_{mN}\right)\right)\in\mathbb{C}^{N\times
N}$ expresses the response of the elements of the $m$th IRS. Note that
$\theta_{mn}\in\left[0,2\pi\right],n=1,\ldots,N$ and $\alpha_{mn}\in(0,1]$
express the phase shifts and the fixed amplitude reflection coefficients
introduced by the corresponding IRS element. The progress on loss-less meta-
surfaces allows to set $\alpha_{mn}=1$, which ensures maximum reflection [5].

##  III Performance Analysis

In this section, we present the derivation of the coverage probability when
multiple IRSs are subject to correlated Rayleigh fading by means of the
deterministic equivalent (DE) analysis. We focus on two interesting scenarios:
a) a finite set of large IRSs ($N\to\infty$); and b) a large number of IRSs
($M\to\infty$) with finite size.

###  III-A Main Results

The coverage probability $\bar{P}_{\mathrm{c}}$ is defined as the probability
that the effective received SNR at the RX is larger than a given threshold
$T$, i.e., $\bar{P}_{\mathrm{c}}=\mathrm{Pr}\left(\gamma>T\right)$, where

$\displaystyle\gamma={\gamma_{0}}{\bigg{|}\displaystyle\sum^{M}_{m=1}{\mathbf{h}}_{m,1}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Phi}}_{m}{\mathbf{h}}_{m,2}+h_{\mathrm{d}}\bigg{|}^{2}}{}$
(5)

is the received SNR in the general case with correlated fading that is
obtained by using (4) and assuming coherent communication. Also,
${\gamma_{0}}=P/N_{0}$ is the average transmit SNR. Under independent Rayleigh
fading with instantaneous CSI, it is known that the phase configuration
$\phi_{m,n}=\arg\left(h_{\mathrm{d}}\right)-\arg\left(h_{mn,1}^{*}\right)\arg\left(h_{mn,2}\right)$
provides the optimal SNR [2, 5]. However, in the practical case of correlated
Rayleigh fading, where only statistical CSI is available, we cannot directly
obtain the solution of the phase shifts. Moreover, since correlated fading
renders the exact derivation of the SNR intractable, we resort to the
application of the DE analysis to derive the approximated SNR 333Note that the
majority of works, deriving the coverage probability in IRS-assisted systems,
result in approximations since they are based on CLT.. In Section IV, we show
that the corresponding $P_{\mathrm{c}}$ provides a tight match with
$\mathrm{Pr}\left(\gamma>T\right)$.

###  III-B Finite $M$ and large $N$ analysis

In this part, we assume large IRSs, as usually considered in the existing
literature to obtain the coverage probability, e.g., see [10, 12].

###### Lemma 1

The SNR of a SISO transmission, enabled by $M$ large IRSs with correlated
Rayleigh fading is approximated by

$\displaystyle\\!{\gamma}$
$\displaystyle\asymp\gamma_{0}\\!\left(B_{M}+|h_{\mathrm{d}}|^{2}\right)\\!,\\!$
(6)

where
$B_{M}=\displaystyle\sum^{M}_{m=1}\\!\beta_{m}\tr\\!\left({\mathbf{R}}_{m,1}{\bm{\Phi}}_{m}{\mathbf{R}}_{m,2}{\bm{\Phi}}_{m}^{{\scriptscriptstyle\mathsf{H}}}\right)$
with $\beta_{m}\\!=\\!\beta_{m,1}$ $\beta_{m,2}$.

###### Proof:

See Appendix A. ∎

###### Proposition 1

The coverage probability of a SISO transmission, enabled by $M$ large IRSs
with correlated Rayleigh fading for arbitrary phase shifts, is tightly
approximated by

$\displaystyle\\!\\!P_{\mathrm{c}}=\left\\{\begin{array}[]{ll}\mathrm{exp}\\!\left(\\!-\frac{1}{\beta_{\mathrm{d}}}\\!\\!\left(\\!\frac{T}{\gamma_{0}}-B_{M}\\!\right)\\!\\!\right)&B_{M}<\frac{T}{\gamma_{0}}\\\
1&B_{M}\geq\frac{T}{\gamma_{0}}\end{array}\right..$ (9)

###### Proof:

See Appendix B. ∎

###### Remark 1

If the aggregate contribution from the IRS-assisted channels is larger than
$T/\gamma_{0}$, no outage is detected during the communication. Also, the
weaker the direct signal $(\beta_{\mathrm{d}}\to 0)$, the less severe its
impact is on $P_{\mathrm{c}}$ and the influence of the cascaded channels
becomes more pronounced. Moreover, when the path-losses of the cascaded links
increase, i.e., $\beta_{m}$ decreases, the coverage decreases.

###### Remark 2

From (9), we observe that when the number of surfaces $M$ increases, the
coverage probability is improved. In addition, by increasing the size of each
IRS in terms of $N$, $P_{\mathrm{c}}$ is enhanced. Hence, the use of more IRSs
or larger IRSs is proved to be beneficial for coverage.

###### Remark 3

Obviously, the coverage probability depends on the phase shifts, which could
be optimized. However, in the case of uncorrelated fading, i.e.,
${\mathbf{R}}_{m,1}={\mathbf{R}}_{m,2}={\bm{\mathrm{I}}}_{N}$,
$P_{\mathrm{c}}$ becomes independent of the reflect beamforming matrices
${\bm{\Phi}}_{m}$. In such case, the phase shifts of the IRS cannot be
optimized to improve the coverage.

###  III-C Large $M$ and finite $N$ analysis

The previous analysis does not allow to examine the coverage when $M\to\infty$
but $N$ is finite. To address this scenario, let
${\mathbf{g}}_{n,1}=\left[h_{1n,1},\ldots,h_{Mn,1}\right]^{{\scriptscriptstyle\mathsf{T}}}\in\mathbb{C}^{M\times
1}$ denote the channel fading vector between the TX and the $n$th elements of
all IRSs (first link). Also,
${\mathbf{g}}_{n,2}=\left[h_{1n,2},\ldots,h_{Mn,2}\right]^{{\scriptscriptstyle\mathsf{T}}}\in\mathbb{C}^{N\times
1}$ expresses the channel between the $n$th elements of all IRSs and the RX
(second link). Given that IRSs are reasonably far apart each other, we assume
no correlation among them444Note that not only this assumption is quite
reasonable but the modeling of a potential correlation would require the
conduct of measurement campaigns, which are currently unavailable.. In other
words, we have
$\mathbb{E}[{\mathbf{h}}_{m,1}{\mathbf{h}}_{l,1}^{{\scriptscriptstyle\mathsf{H}}}]=\mathbb{E}[{\mathbf{h}}_{m,2}{\mathbf{h}}_{l,2}^{{\scriptscriptstyle\mathsf{H}}}]={\mathbf{0}}_{N}~{}\forall~{}m\neq
l$ with $m=1,\ldots,M$ and $l=1,\ldots,M$.

Notably, a correlation appears between different channel vectors at each link.
Specifically, regarding the first link, let ${\mathbf{Q}}_{np,1}$ describe the
correlation between the $n$th and $p$th elements across all surfaces. It can
be written as

$\displaystyle{\mathbf{Q}}_{np,1}$
$\displaystyle=\mathbb{E}[{\mathbf{g}}_{n,1}{\mathbf{g}}_{p,1}^{{\scriptscriptstyle\mathsf{H}}}]$
$\displaystyle={\bm{\beta}}_{1}\text{diag}\left(\right)\left(r_{np,1}^{1},\ldots,r_{np,M}^{1}\right)$

where the matrix
${\bm{\beta}}_{1}=\text{diag}\left(\right)\left(\beta_{1,1},\ldots,\beta_{M,1}\right)\in\mathbb{C}^{M\times
M}$ is diagonal with elements expressing the path-losses between the TX and
the $M$ surfaces. Note that ${\bm{\beta}}_{1}$ does not depend the index $m$
but it includes the corresponding path-losses from all IRSs. The matrix
${\mathbf{Q}}_{np,1}$ is diagonal due to the independence among the IRSs.
Also, $r_{np,i}^{1}$ with $i=1,\ldots,M$ expresses the $(n,p)$th element of
the correlation matrix of the $i$th IRS of the first link, i.e.,
${\mathbf{R}}_{i,1}$. Similarly, for the second link, we have

$\displaystyle{\mathbf{Q}}_{np,2}$
$\displaystyle={\bm{\beta}}_{2}\text{diag}\left(\right)\left(r_{np,1}^{2},\ldots,r_{np,M}^{2}\right)$

where
${\bm{\beta}}_{2}=\text{diag}\left(\right)\left(\beta_{1,1},\ldots,\beta_{M,1}\right)\in\mathbb{C}^{M\times
M}$ is the diagonal matrix expressing the path-losses among the IRSs and the
RX, and
${\mathbf{Q}}_{np,2}=\text{diag}\left(\right)\left(r_{np,1}^{2},\ldots,r_{np,M}^{2}\right)\in\mathbb{C}^{M\times
M}$ describes the corresponding spatial correlation. Notably, in the case of
independent Rayleigh fading,
${\mathbf{Q}}_{np,1}={\mathbf{Q}}_{np,2}={\mathbf{O}}$ for $n\neq p$. As a
result, the corresponding channel vectors of the first and and second links
are formulated as

$\displaystyle{\mathbf{g}}_{n,1}$
$\displaystyle\sim\mathcal{CN}\left({\mathbf{0}},{\bm{\beta}}_{1}{\mathbf{Q}}_{nn,1}\right),$
(10) $\displaystyle{\mathbf{g}}_{n,2}$
$\displaystyle\sim\mathcal{CN}\left({\mathbf{0}},{\bm{\beta}}_{2}{\mathbf{Q}}_{nn,2}\right).$
(11)

The SNR in (5) can be rewritten in terms of a summation over the number of
elements of each IRS as

$\displaystyle\gamma={\gamma_{0}}{\bigg{|}\displaystyle\sum^{N}_{n=1}{\mathbf{g}}_{n,1}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Psi}}_{n}{\mathbf{g}}_{n,2}+h_{\mathrm{d}}\bigg{|}^{2}}{},$
(12)

where
${\bm{\Psi}}_{n}=\mathrm{diag}\left(\exp\left(j\theta_{1n}\right),\ldots,\exp\left(j\theta_{Mn}\right)\right)\in\mathbb{C}^{M\times
M}$.

###### Lemma 2

The SNR of a SISO transmission, enabled by a large number of finite size IRSs
with correlated Rayleigh fading is approximated by

$\displaystyle{\gamma}$
$\displaystyle\asymp{\gamma_{0}}\left(B_{N}+|h_{\mathrm{d}}|^{2}\right)\\!,$
(13)

where
$B_{N}=\displaystyle\sum^{N}_{n=1}\\!\sum^{N}_{p=1}\\!\tr\\!\left({\mathbf{Q}}_{np,1}{\bm{\Psi}}_{n}{\mathbf{Q}}_{np,2}{\bm{\Psi}}_{p}^{{\scriptscriptstyle\mathsf{H}}}\right)$.
Note that ${\bm{\Psi}}=\text{diag}\left(\right)({\bm{\Psi}}_{1},$
$\ldots,{\bm{\Psi}}_{N})\in\mathbb{C}^{MN\times MN}$, i.e., ${\bm{\Psi}}$ is a
block diagonal matrix.

###### Proof:

See Appendix C. ∎

###### Proposition 2

The coverage probability of a SISO transmission, enabled by a large number of
finite size IRSs with correlated Rayleigh fading for arbitrary phase shifts,
is tightly approximated by

$\displaystyle\\!\\!P_{\mathrm{c}}=\left\\{\begin{array}[]{ll}\mathrm{exp}\\!\left(\\!-\frac{1}{\beta_{\mathrm{d}}}\\!\\!\left(\\!\frac{T}{\gamma_{0}}-B_{N}\\!\right)\\!\\!\right)&B_{N}<\frac{T}{\gamma_{0}}\\\
1&B_{N}\geq\frac{T}{\gamma_{0}}\end{array}\right.\\!.$ (16)

###### Proof:

The proof follows similar lines with the proof of Proposition 1. ∎

###### Remark 4

We observe that (16) has a similar expression with (9). However, the main
characteristic of (16) is that it is written in terms of a double summation
expressing the correlation among the IRS elements instead of one summation in
(9). Notably, if the correlation matrix is identical across different IRSs,
${\mathbf{Q}}_{np,1}$ and ${\mathbf{Q}}_{np,2}$ are scaled identity matrices
but the coverage will always be dependent on the phase shifts due to the
contributions from the off-diagonal terms of the IRSs. Furthermore, we observe
a similar dependence from the path-loss of the direct signal and the number of
surfaces and their elements, i.e., their increase improves the coverage. Also,
under uncorrelated Rayleigh fading conditions, $P_{\mathrm{c}}$ does not
depend on the phases, and thus, cannot be optimized.

###  III-D Reflecting beamforming optimization

Both Propositions (1) and (2) are described by a similar expression in terms
of a trace that includes the reflecting beamforming matrices. Hence, their
optimization follows similar steps up to a point. Specifically, to achieve
maximum $P_{\mathrm{c}}$, we formulate the optimization problem, relying on
the common assumption of infinite resolution phase shifters, as

$\displaystyle\begin{split}&\\!\\!\\!(\mathcal{P}1)~{}\max_{{\bm{\Phi}}}~{}~{}P_{\mathrm{c}}\\\
&~{}~{}~{}~{}~{}~{}\mathrm{s.t}~{}|\phi_{mn}|\\!=\\!1,~{}~{}m\\!=\\!1,\dots,M~{}\mathrm{and}~{}n\\!=\\!1,\dots,N,\end{split}$
(17)

where $P_{\mathrm{c}}$ is given by (9) or (16) and
$\phi_{mn}=\exp\left(j\theta_{mn}\right)$.

The optimization problem $(\mathcal{P}1)$ is non-convex with respect to the
reflect beamforming matrix while having a unit-modulus constraint regarding
$\phi_{mn}$. Use of projected gradient ascent until converging to a stationary
point can provide a direct solution. In particular, since each surface has a
similar solution, we focus on the $m$th IRS. At the $i$th step, we assume the
vectors
${\mathbf{s}}_{m,i}=[\phi_{m1}^{i},\ldots,\phi_{mN}^{i}]^{{\scriptscriptstyle\mathsf{T}}}$,
which include the phases at this step. The next iteration increases
$P_{\mathrm{c}}$ until its convergence by projecting the solution onto the
closest feasible point based on
$\min_{|\phi_{mn}|=1,n=1,\ldots,N}\|{\mathbf{s}}_{m}-\tilde{{\mathbf{s}}}_{m}\|^{2}$
satisfying the unit-modulus constraint concerning $\phi_{mn}$ with

$\displaystyle\tilde{{\mathbf{s}}}_{m,i+1}$
$\displaystyle={\mathbf{s}}_{m,i}+\mu{\mathbf{q}}_{m.i},$ (18)
$\displaystyle{\mathbf{s}}_{m,i+1}$
$\displaystyle=\exp\left(j\arg\left(\tilde{{\mathbf{s}}}_{m,i+1}\right)\right).$
(19)

Note that $\mu$ expresses the step size computed at each iteration by means of
the backtracking line search [17] while ${\mathbf{q}}_{m,i}$ denotes the
adopted ascent direction at step $i$ with
${\mathbf{q}}_{m,i}=\partialderivative{P_{\mathrm{c}}}{\phi_{mn}}$, obtained
by Lemma 3 below. Algorithm III-D provides an outline of the proposed
algorithm for Proposition 1 and 2 by setting
$\tilde{{\bm{\Phi}}}={\bm{\Phi}}_{m}$ and
$\tilde{{\bm{\Phi}}}={\bm{\Psi}}_{n}$, respectively. 2\. Iteration $i$: for
$i=0,1,\dots,$ do 3\.
$[{\mathbf{q}}_{m,i}]_{n}=\partialderivative{P_{\mathrm{c}}}{{\mathbf{s}}_{m,i}^{*}}$,
where $\partialderivative{P_{\mathrm{c}}}{{\mathbf{s}}_{m,i}^{*}}$ is given by
Lemma 3; 4\. Find $\mu$ by backtrack line
search$(f\left(\tilde{{\bm{\Phi}}}_{0}\right),{\mathbf{q}}_{m,i},{\mathbf{s}}_{m,i})$
[17]; 5\.
$\tilde{{\mathbf{s}}}_{m,i+1}={\mathbf{s}}_{m,i}+\mu{\mathbf{q}}_{m,i}$; 6\.
${\mathbf{s}}_{m,i+1}=\exp\left(j\arg\left(\tilde{{\mathbf{s}}}_{m,i+1}\right)\right)$;
$\tilde{{\bm{\Phi}}}_{i+1}=\text{diag}\left(\right)\left({\mathbf{s}}_{m,i+1}\right)$;
7\. $P_{\mathrm{c}}^{i+1}=f\left(\tilde{{\bm{\Phi}}}_{i+1}\right)$; 8\. Until
$\|P_{\mathrm{c}}^{i+1}-P_{\mathrm{c}}^{i}\|^{2}<\epsilon$; Obtain
$\tilde{{\bm{\Phi}}}^{*}=\tilde{{\bm{\Phi}}}_{i+1}$; 9\. end for

###### Lemma 3

The derivative of the coverage probability with respect to
${\mathbf{s}}_{m,i}^{*}$ is given by

$\displaystyle\\!\partialderivative{{P_{\mathrm{c}}}}{{\mathbf{s}}_{m,i}^{*}}\\!\\!=\\!\frac{\beta_{m}}{\beta_{\mathrm{d}}}\\!\left\\{\\!\\!\\!\\!\begin{array}[]{ll}{P_{\mathrm{c}}}\,\text{diag}\left(\right)\left({\mathbf{R}}_{m,1}{\bm{\Phi}}_{m}{\mathbf{R}}_{m,2}\right),&\mathrm{Prop.}~{}\ref{GeneralPDF}\\\
{P_{\mathrm{c}}}\displaystyle\\!\sum^{N}_{p=1}{\mathbf{c}}_{p},&\mathrm{Prop.}~{}\ref{GeneralPDF1}\end{array}\right.\\!\\!$
(22)

when $B_{i}<\frac{T}{\gamma_{0}},~{}i=M,N$. Otherwise, it is zero. Note that
${\mathbf{c}}_{p}=\left[r_{1p,m}^{1}r_{1p,m}^{2}\phi_{m,1n},\ldots,r_{Np,m}^{1}r_{Np,m}^{2}\phi_{m,1N}\right]^{{\scriptscriptstyle\mathsf{T}}}$.

###### Proof:

See Appendix D. ∎

##  IV Numerical Results

Relying on a Cartesian coordinate system, we consider a cell, where the TX is
located at the origin and the RX at $\left(60,0\right)$. Also, we assume
$M=15$ IRSs being uniformly distributed between them while each IRS is
deployed with $N=225$ elements unless otherwise specified. The size of each
IRS element is given by
$d_{\mathrm{H}}\\!=\\!d_{\mathrm{V}}\\!=\\!\lambda/32$. The spatial
correlation matrix is given by (III-D). The large-scale fading coefficients
between the TX and the RX are given by
$\beta_{m,i}=G_{t}+G_{r}+10\nu_{e}\log_{10}(d_{i}/1\mathrm{m})-27.5$, where
$i\in\\{1,2\\}$ while $\beta_{\mathrm{d}}$ is given similarly. Also, the path-
loss exponents are $\nu_{1}=\nu_{2}=2$ and $\nu_{\mathrm{d}}=3.5$. Moreover,
we have $G_{t}=3.2$ dBi and $G_{r}=1.3$ dBi, the system bandwidth is $10$ MHz;
the carrier frequency is $3$ GHz and the noise variance is $-94$ dBm with the
noise figure being $10$ dB. Note that the transmitter power is $10$ dBm.
Monte-Carlo (MC) simulations verify the analytical results and corroborate
that the DE analysis provides tight approximations as has been already shown
in the literature, e.g., see [18].

Fig. 1 shows the coverage probability versus the target rate in the scenario
described by Proposition 1555The theoretical analysis of this proposition
relies on finite $N$ but we consider $N\geq 225$, which is common for
practical IRS implementations.. By increasing the number of elements in each
IRS, $P_{\mathrm{c}}$ increases. Also, the addition of more IRSs (an increase
of $M$) through $B_{M}$ contributes to the observation of less outage during
communication. Moreover, if no correlation is assumed, $P_{\mathrm{c}}$ is
lower because it becomes independent of the reflect beamforming matrix and
cannot be optimized. Notably, if the impact from the direct signal through
$\beta_{\mathrm{d}}$ becomes weaker, the coverage decreases, and the variation
regarding the number of elements as well as correlated Rayleigh fading have a
greater impact on $P_{\mathrm{c}}$ since the relevant gaps are larger.

Figure 1: Coverage probability of a SISO system with correlated Rayleigh
fading assisted by $M$ IRSs each having a large number of elements
($N\to\infty$) versus the target rate $T$ (analytical results and MC
simulations).

In Fig. 2, we depict the coverage probability versus the target rate by
accounting for a large number of surfaces, i.e., we shed light on the setting
referring to Proposition 2. Herein, we notice that the coverage is improved as
the number of IRSs increases. Similarly, if we increase the number of elements
per IRS, we observe a further improvement. Notably, a comparison between Figs.
1 and 2 reveal that an analogous increase concerning the number of IRSs
results in a larger improvement of the coverage compared to increasing the
number of elements per IRS. Specifically, in Fig. 2, it is shown that when $M$
increases from $15$ to $34$ b/s/Hz ($225\%$ increase), $P_{\mathrm{c}}$ starts
decreasing from full coverage when $T=1.47~{}\mathrm{dB}$ and
$2.2~{}\mathrm{dB}$, respectively. On the other hand, in Fig. 1, we observe
that for a similar increase concerning the number of elements per IRS, i.e.,
when $N$ increases from $100$ to $225$, the coverage is much lower.

Figure 2: Coverage probability of a SISO system with correlated Rayleigh
fading assisted by a large number of IRSs ($M\to\infty$) each having $N$
elements versus the target rate $T$ (analytical results and MC simulations).

##  V Conclusion

In this paper, we derived the coverage probability of a SISO system assisted
with multiple IRSs under the unavoidable conditions of correlated Rayleigh
fading. We considered two distinct scenarios: a finite multitude of large IRSs
and a large number of finite IRSs. Especially, we managed to derive and
optimize the coverage probability with respect to the phase shifts of the IRS
elements in both cases. The results enabled us to show that it is more
beneficial to increase the number of IRSs instead of increasing their
elements. Future works on coverage of distributed IRSs should take into
account the design of multi-user transmission with multiple antennas, and
possibly, under Rician fading conditions.

##  Appendix A Proof of Lemma 1

The proof starts by dividing (5) with $\frac{1}{N^{2}}$. Then, we have

$\displaystyle\\!\\!\frac{1}{N^{2}}{\gamma}=\gamma_{0}\frac{1}{N^{2}}\bigg{(}\displaystyle\sum^{M}_{m=1}\bigg{|}{\mathbf{h}}_{m,1}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Phi}}_{m}{\mathbf{h}}_{m,2}\bigg{|}^{2}+|h_{\mathrm{d}}|^{2}$
$\displaystyle\\!\\!+2\mathrm{Re}\left(h_{\mathrm{d}}^{*}\displaystyle\sum^{M}_{m=1}{\mathbf{h}}_{m,1}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Phi}}_{m}{\mathbf{h}}_{m,2}\right)$
$\displaystyle\\!\\!+\displaystyle\sum^{M}_{m=1}\sum^{M}_{\begin{subarray}{c}j=1\\\
n\neq
m\end{subarray}}{\mathbf{h}}_{m,1}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Phi}}_{m}{\mathbf{h}}_{m,2}{\mathbf{h}}_{n,2}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Phi}}_{n}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{h}}_{n,1}\bigg{)}$
(23)
$\displaystyle\\!\\!\asymp\gamma_{0}\frac{1}{N^{2}}\\!\\!\left(\displaystyle\sum^{M}_{m=1}\\!\beta_{m,1}\beta_{m,2}\\!\tr\\!\left({\mathbf{R}}_{m,1}{\bm{\Phi}}_{m}{\mathbf{R}}_{m,2}{\bm{\Phi}}_{m}^{{\scriptscriptstyle\mathsf{H}}}\right)\\!+\\!|h_{\mathrm{d}}|^{2}\\!\\!\right)\\!\\!,$
(24)

where, in (24), we have used [18, Lem. 4] for the first, third, and the fourth
terms. Especially, the third and fourth terms in (23) vanish as $N\to\infty$
due to the independence among the $M$ IRSs and between the two links,
respectively.

##  Appendix B Proof of Proposition 1

The coverage probability is written as

$\displaystyle
P_{\mathrm{c}}\\!=\\!\mathrm{Pr}\\!\left(\\!|h_{\mathrm{d}}|^{2}>\frac{T}{\gamma_{0}}\\!-\\!B_{M}\\!\\!\right)\\!,$
(25)

where in (25), we have used the SNR from (6). Given that
$|h_{\mathrm{d}}|^{2}$ is exponentially distributed with rate parameter
$\beta_{\mathrm{d}}$, i.e.,
$|h_{\mathrm{d}}|^{2}\sim\mathrm{Exp}(1/\beta_{\mathrm{d}})$, we obtain the
first branch in (9), if $B_{M}<\frac{T}{\gamma_{0}}$. Otherwise,
$P_{\mathrm{c}}=1$, and we conclude the proof.

##  Appendix C Proof of Lemma 2

By dividing $\gamma$ with $\frac{1}{M^{2}}$, we have

$\displaystyle\\!\\!\\!\frac{1}{M^{2}}\gamma\\!=\\!{\gamma_{0}}\frac{1}{M^{2}}\\!\bigg{(}\displaystyle\\!\sum^{N}_{n=1}\\!\bigg{|}{\mathbf{g}}_{n,1}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Psi}}_{n}{\mathbf{g}}_{n,2}\bigg{|}^{2}\\!\\!\\!+\\!2\mathrm{Re}\\!\left(\\!\\!h_{\mathrm{d}}^{*}\\!\displaystyle\sum^{N}_{n=1}\\!{\mathbf{g}}_{n,1}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Psi}}_{n}{\mathbf{g}}_{n,2}\\!\\!\right)$
$\displaystyle+|h_{\mathrm{d}}|^{2}+\displaystyle\sum^{N}_{n=1}\sum^{N}_{\begin{subarray}{c}p=1\\\
p\neq
n\end{subarray}}{\mathbf{g}}_{n,1}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Psi}}_{n}{\mathbf{g}}_{n,2}{\mathbf{g}}_{p,2}^{{\scriptscriptstyle\mathsf{H}}}{\bm{\Psi}}_{p}^{{\scriptscriptstyle\mathsf{H}}}{\mathbf{g}}_{p,1}\bigg{)}$
(26)
$\displaystyle={\gamma_{0}}\frac{1}{M^{2}}\bigg{(}\displaystyle\sum^{N}_{n=1}\sum^{N}_{p=1}\tr\left({\mathbf{Q}}_{np,1}{\bm{\Psi}}_{n}{\mathbf{Q}}_{np,2}{\bm{\Psi}}_{p}^{{\scriptscriptstyle\mathsf{H}}}\right)+|h_{\mathrm{d}}|^{2}\bigg{)},$
(27)

where the second term in (26) vanishes as $M\to\infty$ due to the independence
between the two links. Application of [18, Lem. 4] at the first and fourth
terms in (26) gives (27) after a direct combination of the two resultant
traces.

##  Appendix D Proof of Lemma 3

First, by applying the chain rule, we obtain

$\displaystyle\partialderivative{{P_{\mathrm{c}}}}{{\mathbf{s}}_{m,i}^{*}}$
$\displaystyle=\partialderivative{{P_{\mathrm{c}}}}{B_{i}}\partialderivative{B_{i}}{{\mathbf{s}}_{m,i}^{*}},$
(28)

where $i=M$ and $i=N$ correspond to Propositions 1 and Propositions 2,
respectively. The first-order derivative in (28) becomes

$\displaystyle\partialderivative{{P_{\mathrm{c}}}}{B_{i}}$
$\displaystyle=\frac{1}{\beta_{\mathrm{d}}}{P_{\mathrm{c}}}.$ (29)

By taking into account the expression of $B_{M}$, we have

where we have used the property
$\tr\left({\mathbf{A}}\text{diag}\left(\right)({\mathbf{s}}_{m})\right)=\left(\text{diag}\left(\right)(A)\right)^{{\scriptscriptstyle\mathsf{T}}}{\mathbf{s}}_{m}$.
In the case of $\partialderivative{B_{N}}{{\mathbf{s}}_{m,i}^{*}}$, we obtain

(30)

By substituting (LABEL:deriv4) or (30) together with (29) into (28), we obtain
the desired results.

## References

  * [1] Q. Wu and R. Zhang, “Towards smart and reconfigurable environment: Intelligent reflecting surface aided wireless network,” _IEEE Commun. Mag._ , vol. 58, no. 1, pp. 106–112.
  * [2] E. Basar _et al._ , “Wireless communications through reconfigurable intelligent surfaces,” _IEEE Access_ , vol. 7, pp. 116 753–116 773, 2019\.
  * [3] C. Huang _et al._ , “Reconfigurable intelligent surfaces for energy efficiency in wireless communication,” _IEEE Trans. Wireless Commun._ , vol. 18, no. 8, pp. 4157–4170, 2019.
  * [4] C. Pan _et al._ , “Multicell MIMO communications relying on intelligent reflecting surfaces,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 8, pp. 5218–5233, 2020.
  * [5] E. Björnson, Ö. Özdogan, and E. G. Larsson, “Intelligent reflecting surface versus decode-and-forward: How large surfaces are needed to beat relaying?” vol. 9, no. 2, pp. 244–248.
  * [6] Q. U. A. Nadeem _et al._ , “Asymptotic max-min SINR analysis of reconfigurable intelligent surface assisted MISO systems,” _IEEE Trans. Wireless Commun._ , vol. 19, no. 12, pp. 7748–7764, 2020.
  * [7] A. M. Elbir _et al._ , “Deep channel learning for large intelligent surfaces aided mm-Wave massive MIMO systems,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 9, pp. 1447–1451, 2020.
  * [8] C. Guo _et al._ , “Outage probability analysis and minimization in intelligent reflecting surface-assisted MISO systems,” _IEEE Commun. Lett._ , vol. 24, no. 7, pp. 1563–1567, 2020.
  * [9] Y. Gao _et al._ , “Distributed IRS with statistical passive beamforming for MISO communications,” _IEEE Wireless Commun. Lett._ , pp. 1–1, 2020.
  * [10] L. Yang _et al._ , “Coverage, probability of SNR gain, and DOR analysis of RIS-aided communication systems,” _IEEE Wireless Commun. Lett._ , vol. 9, no. 8, pp. 1268–1272, 2020.
  * [11] T. Van Chien _et al._ , “Coverage probability and ergodic capacity of intelligent reflecting surface-enhanced communication systems,” _IEEE Commun. Lett._ , 2020.
  * [12] Z. Zhang _et al._ , “Analysis and optimization of outage probability in multi-intelligent reflecting surface-assisted systems,” _arXiv preprint arXiv:1909.02193_ , 2019.
  * [13] S. Sun _et al._ , “Towards reconfigurable intelligent surfaces powered green wireless networks,” in _IEEE Wireless Communications and Networking Conference (WCNC)_ , 2020, pp. 1–6.
  * [14] Z. Zhang and L. Dai, “Capacity improvement in wideband reconfigurable intelligent surface-aided cell-free network,” in _IEEE 21st International Workshop on Signal Processing Advances in Wireless Communications (SPAWC)_ , 2020, pp. 1–5.
  * [15] E. Björnson and L. Sanguinetti, “Rayleigh fading modeling and channel hardening for reconfigurable intelligent surfaces,” _IEEE Wireless Commun. Lett._ , 2020.
  * [16] D. Neumann, M. Joham, and W. Utschick, “Covariance matrix estimation in massive MIMO,” _IEEE Signal Process. Lett._ , vol. 25, no. 6, pp. 863–867, 2018.
  * [17] S. Boyd, S. P. Boyd, and L. Vandenberghe, _Convex optimization_. Cambridge university press, 2004.
  * [18] A. K. Papazafeiropoulos and T. Ratnarajah, “Deterministic equivalent performance analysis of time-varying massive MIMO systems,” _IEEE Trans. Wireless Commun._ , vol. 14, no. 10, pp. 5795–5809, 2015.

solutionfile


# Diversity of Ultrafast Spin Dynamics Near the Tricritical Point in a
Ferrimagnetic Gd/FeCo Multilayer

T.G.H. Blank  Radboud University, Institute for Molecules and Materials, 6525
AJ Nijmegen, the Netherlands  Department of Applied Physics, Eindhoven
University of Technology, P.O. Box 513, Eindhoven 5600 MB, the Netherlands
B.D. Muis  Kavli Institute of Nanoscience, Delft University of Technology,
2600 GA Delft, the Netherlands  T. Lichtenberg  Department of Applied Physics,
Eindhoven University of Technology, P.O. Box 513, Eindhoven 5600 MB, the
Netherlands  B. Koopmans  Department of Applied Physics, Eindhoven University
of Technology, P.O. Box 513, Eindhoven 5600 MB, the Netherlands  A.V. Kimel
Radboud University, Institute for Molecules and Materials, 6525 AJ Nijmegen,
the Netherlands

###### Abstract

It is found that subtle changes in the external magnetic field and temperature
result in dramatic changes in the ultrafast response of spins to a femtosecond
laser excitation in a ferrimagnetic Gd/FeCo multilayer. A total of six
distinct types of spin dynamics were observed and explained by considering the
spin-flop transition to the noncollinear phase and the concept of a
tricritical point in the $H$-$T$ phase diagram. A particularly interesting
type of dynamics is the exchange-driven reversal. These exchange-driven
dynamics provide new insights into the tricritical point, which is shown to
separate two thermodynamically distinct noncollinear phases with the
transition-metal magnetization pointing on adjacent sides of the anisotropy
plane.

††preprint: APS/123-QED

Ultrafast magnetism is a rapidly developing field of physics that explores
spin dynamics launched in magnetically ordered materials by ultrashort stimuli
with a duration of just a few ps or less. Such stimuli bring the magnetic
media into a strongly non-equilibrium state where the conventional
thermodynamic description of magnetic phenomena fails, often resulting in
counter-intuitive spin dynamics [1, 2, 3]. Better fundamental understanding
and exploring the practical limits on the timescales to control the magnetic
ordering is expected to have a great impact on the development of magnetic
information processing and data storage technologies [4, 5].

Among all classes of magnetic materials studied in ultrafast magnetism
(ferro-, ferri-, and antiferromagnets), ferrimagnets are most promising for
the ultrafast and efficient control of magnetic properties. In these
materials, the exchange interaction aligns the spins antiparallel as in
antiferromagnets, but the spins are not equivalent such that the net
magnetization is nonzero. Examples of ultrafast phenomena discovered in
ferrimagnets include magnetization reversal via a strongly non-equilibrium
state initiated by fs laser pulses [6, 7], current-induced magnetization
reversal [8], as well as the record fast and least dissipative magnetic
recording [9]. A central ingredient that enriches ultrafast magnetism in
ferrimagnets is the fact that the balance between the spins and therefore
their magnetic structure changes with temperature. In particular, a
compensation temperature may exist where the spins cancel out and the net
magnetization is zero. Ferrimagnets in the vicinity of their compensation
temperature combine the properties of ferro- and antiferromagnets and
therefore offer a highly intriguing playground in ultrafast magnetism [10, 11,
12].

Lesser known is that, next to temperature, the magnetic structure of
ferrimagnets can also be tuned with an external magnetic field, resulting in a
rich $H$-$T$ phase diagram [13]. When a critical field is reached, the spins
undergo a spin-flop transition to a noncollinear spin configuration [14].
Previous studies of ferrimagnetic rare-earth (RE) transition-metal (TM) alloys
have indeed shown that sufficiently high external magnetic fields unlock a new
dimension in the ultrafast magnetism of ferrimagnets related to this
noncollinear phase [15, 16]. However, due to relatively strong exchange
interaction between the sublattices, the studies required exceptionally high
magnetic fields up to $30$ T. Here, we overcome this obstacle by fabricating a
high-quality synthetic ferrimagnet – a Gd/FeCo multilayer. In such a
heterostructure, the amount of nearest neighbors from different species
(RE/TM) decreases as compared to the alloy, meaning that the RE-TM exchange
interaction, as well as the critical field, are expected to drop significantly
[17]. In this Letter, we show in Gd/FeCo that by changing the temperature or
magnetic field in relatively narrow ranges, from $320$ K to $340$ K and from
$0.5$ T to $1$ T, we can reach all possible magnetic phases hosted by
ferrimagnets and trigger largely diverse ultrafast spin dynamics upon
femtosecond laser pulse excitation. We propose to explain this diversity with
the notion of a tricritical point.

The studied multilayer structure Glass/Ta($4$)/Al($10$)/
Ta($4$)/[Fe87.5Co12.5($0.5$)/Gd($0.5$)]×20/Ta($4$), with layer thicknesses (in
nm) indicated between parentheses, was grown by magnetron sputtering at an Ar
pressure of $10^{-2}$ mbar. Further details on the growth procedure can be
found in Ref. [18]. The atomic magnetic moments (to which we simply refer as
“spins”) in the ferromagnetic FeCo layers with total magnetization
$\mathbf{M}_{\mathrm{TM}}$ couple antiferromagnetically to the spins in the
ferromagnetic Gd layers with magnetization $\mathbf{M}_{\mathrm{RE}}$, as is
depicted schematically in Fig. 1(a). In the absence of an external magnetic
field, the spins align along the out-of-plane anisotropy axis. Three
magnetostatic phases can be distinguished in this material (see Fig. 1(b)).
For low temperatures, the absolute Gd magnetization $M_{\mathrm{RE}}$ exceeds
that of the FeCo layers $M_{\mathrm{TM}}$ (Phase I). As $M_{\mathrm{RE}}(T)$
decreases faster with increasing temperature than $M_{\mathrm{TM}}(T)$, a
compensation temperature $T_{\mathrm{M}}$ exists at which
$M_{\mathrm{TM}}=M_{\mathrm{RE}}$, and above which
$M_{\mathrm{TM}}>M_{\mathrm{RE}}$ (phase II). Moreover, when an external
magnetic field $\mathbf{B}_{\mathrm{ext}}$ is applied along the anisotropy
axis, a critical spin-flop field $B_{\mathrm{SF}}$ can be reached at which the
gain in Zeeman energy due to an increase of net magnetization
$\mathbf{M}=\mathbf{M}_{\mathrm{RE}}+\mathbf{M}_{\mathrm{TM}}$ when the spins
would arrange noncollinearly exceeds its associated loss in exchange and
anisotropy energy. In this case, the noncollinear phase (phase III) is
attained, described by the polar angles $\theta_{\mathrm{TM}}$ and
$\theta_{\mathrm{RE}}$ of the sublattice magnetizations. This spin-flop phase-
transition (PT) from phases I/II to phase III can be either of the first, or
second-order [19]. In terms of Landau’s theory of PTs, in a first-order PT,
the changes of the order parameter ($\theta_{\mathrm{TM}}$) are discontinuous
and the transition between the phases often exhibits hysteresis. In a second-
order PT, the order parameter changes continuously and without hysteresis [19,
16, 7].

Figure 1: (a) Illustration of the spins in the multilayer. The Gd moments
(purple, hollow) align antiparallel to those of the FeCo layers (red, filled).
(b) The arrows depict the sublattice magnetizations of Gd and FeCo in the
three static phases: two collinear (I and II) and one noncollinear phase
(III).

To pinpoint $T_{\mathrm{M}}$ and $B_{\mathrm{SF}}(T)$ experimentally, we
measured the magneto-optical Faraday effect of Gd/FeCo using linearly
polarized light from a HeNe laser ($\lambda=632.8$ nm) that transmitted
perpendicularly through the multilayer. As commonly accepted, we assume that
the Faraday rotation in GdFeCo at the wavelengths $600$-$800$ nm is dominated
by the out-of-plane projection of the TM magnetization
$M_{\mathrm{TM},z}=-|\mathbf{M}_{\mathrm{TM}}|\cos\theta_{\mathrm{TM}}$, and
we will take this assumption throughout the article [20, 21]. The magneto-
optical dominance of $\mathbf{M}_{\mathrm{TM}}$ is corroborated by the fact
that the magneto-optical signal does not decrease/disappear at the
magnetization compensation point and moreover remains unchanged in the
temperature range $293-350$ K (see Supplemental Material [22]). A selection of
hysteresis loops measured using external magnetic fields up to $\sim 1$ T
applied along the sample normal is shown in Fig. 2(a). The step around
$B_{\mathrm{ext}}=0$ corresponds to a $180^{\circ}$ magnetization reversal.
The diverging coercive field and the flip of the loops between $329$ K and
$330$ K are clear indications of a compensation temperature
$T_{\mathrm{M}}\approx 330$ K. Moreover, the data measured at 310 K reveal a
transition from a fully saturated state, which is observed at small fields, to
a state with a reduced magneto-optical signal, which is observed at larger
fields i.e. above $\sim 0.6$ T. This transition sets in abruptly, but is
continuous. Thus, the corresponding change in the order parameter
$\theta_{\mathrm{TM}}$ can be classified thermodynamically as a second-order
PT. The critical field depends on temperature and reaches a minimum at the
compensation point. When the temperature is slightly above the compensation
point, the continuous change has become a discontinuous jump, while also
hysteresis can be observed. These features are characteristic of a first-order
PT, and it persists for about $5$ K above $T_{\mathrm{M}}$ after which the PT
changes back to the second-order. The point with temperature
$T_{\mathrm{3P}}\approx 335$ K at which the first and second-order PTs meet is
by convention called a tricritical point.

Figure 2: (a) Faraday rotation measurements of light at the wavelength
$\lambda=632.8$ nm. The bending and discontinuities of the loops correspond to
the spin-flop transition. (b) Critical field $B_{\mathrm{SF}}$ and order of
the spin-flop transition deduced from the isothermal susceptibility
$\chi_{T}$.  Figure 3: (a) Illustration of the experimental pump-probe setup
(abbreviations: beamsplitter (BS), and long-pass filter (LP)). (b) The six
different types of magnetization dynamics that were distinguished in the pump-
probe measurements. The rotation data divided by $\theta_{0}$ can be
interpreted as $\theta_{\mathrm{F}}(t)/\theta_{0}=\Delta
M_{\mathrm{TM,z}}(t)/|M_{\mathrm{TM}}(0)|$. The arrows provide a pictorial
description of $\mathbf{M}_{\mathrm{TM}}$, based on our interpretation of the
data. (c) Overview of where each type of dynamics (see (b)) occurs in a
$H$-$T$ diagram. Note that $T_{\mathrm{M}}$ shifted $\sim 5$ K compared to the
results in Fig. 2 due to heat accumulation in the sample. The purple area
denotes the region where the amount of UD decreased below $20\%$ of the
maximal demagnetization measured at low fields ($B_{\mathrm{ext}}<0.5$ T),
where the system started in the collinear state. It can be seen that the
purple region terminates near the tricritical point.

The exact locations of the second-order PTs were determined by taking the
center of a sigmoid function that was fitted to the isothermal susceptibility
$\chi_{T}=\left(\frac{\partial M_{\mathrm{TM}}}{\partial B}\right)_{T}$.
Analogously, a first-order PT is marked by a divergence of this
susceptibility. The resulting critical field $B_{\mathrm{SF}}$ as a function
of $T$ is summarized in Fig. 2(b). The temperature dependence of the critical
field is in perfect agreement with the spin-flop transition which has been
observed and predicted for GdFeCo [15, 19] and TbFeCo [16] alloys, except that
the typically required field decreased nearly an order of magnitude. The
results are supported by measurements using a SQUID-VSM, which probes net
magnetization $\mathbf{M}$ instead of $M_{\mathrm{TM},z}$ (see Supplemental
Material 111See Supplemental Material at [URL] for the full static magnetic
characterization of the Gd/FeCo sample using magneto-optics and the SQUID-VSM,
details of the optical pump-probe experimental method, all of the measured
dynamical transients, estimations of the temperature increase due to laser-
induced heating, measurements of UD with $100$ fs time-resolution and an
example of fluence dependence of the data, which includes Refs. [26, 27].)).

Theoretical work showed that the tricritical point should appear either above
or below $T_{\mathrm{M}}$ depending on whether the magnetic anisotropy of the
RE sublattice is larger or smaller than that of the TM sublattice [19]. Yet,
experimental studies of the noncollinear state in GdFeCo and TbFeCo alloys
have not revealed the role of the tricritical point, due to the experimentally
challenging requirement of high magnetic fields. Having manufactured a Gd/FeCo
multilayer whose tricritical point appears well below $1$ T, we could study
the laser-induced spin dynamics in detail using a table-top optical pump-probe
setup with external field applied out of the sample plane as depicted
schematically in Fig. 3(a) (see Supplementary Material [22] for details of the
method). The pump-induced changes of the Faraday rotation
$\theta_{\mathrm{F}}(t)$ is again assumed to be proportional to
$M_{\mathrm{TM},z}$. We divided $\theta_{\mathrm{F}}(t)$ by the static
magneto-optical rotation $\theta_{0}$ associated with
$\mathbf{M}_{\mathrm{TM}}$ pointing completely out of the plane, i.e.
$2\theta_{0}\approx 0.9^{\circ}$ equals the Faraday rotation contrast (at
$\lambda=800$ nm) in the case of a full magnetization reversal (see Fig.
2(a)). This ratio equals the normalized changes of the out-of-plane projection
of the FeCo magnetization $\Delta
M_{\mathrm{TM,z}}(t)/|\mathbf{M}_{\mathrm{TM}}(0)|=\theta_{\mathrm{F}}(t)/\theta_{0}$.
For example, $\theta_{\mathrm{F}}(t)/\theta_{0}=2$ implies a full
$180^{\circ}$ reversal of $\mathbf{M}_{\mathrm{TM}}$.

Figure 4: Several quantities derived from the pump-probe measurements: (a) The
initial polar angle of the FeCo magnetization $\theta_{\mathrm{TM}}(t<0)$. The
dotted line indicates where $\theta_{\mathrm{TM}}=\pi/2$. (b) The final value
of $\theta_{\mathrm{TM}}$ after $1.2$ ns. The outcome is prone to several
measurement and calibration uncertainties. (c) The typical timescale
associated with the sub-ns dynamics of $\mathbf{M}_{\mathrm{TM}}(t)$. No
reversal was observed below the dotted line on the left and in the blue region
on the right. (d) Frequencies of the oscillations measured at $T=325$ K as a
function of the external magnetic field. The slope can be related to the
effective $g$-factor, excluding the soft mode [23] near the spin-flop PT at
$B_{\mathrm{ext}}=0.6$ T.

Using this method, we observed a remarkable diversity in magnetization
dynamics for different external fields and temperatures. In particular, we
were able to distinguish six different types of dynamics depicted in Fig.
3(b). The curves of type $1$a and $1$b correspond to a sub-ps laser-induced
partial loss of magnetic order, known as ultrafast demagnetization (UD) [1],
of about $\sim 30~{}\%$ followed by a recovery of the FeCo-sublattice,
measured below and above the magnetization compensation point for the Gd-
dominated and FeCo-dominated phases I and II, respectively. Type $2$ can be
assigned to UD of FeCo followed by a slow (sub-ns) reversal of the
magnetization. The dynamics of types $3$ and $5$ have a reduced
demagnetization compared to the types $1$ and $2$, followed by an ultrafast
reorientation of the spins accompanied by spin oscillations. Type $3$ and $5$
could only be observed for high fields ($>0.5$ T), therefore we attribute the
decrease in demagnetization to the reduced $z$-projection
($\propto\cos\theta_{\mathrm{TM}}$) of the magnetization in a noncollinear
configuration. Also, the UD has an opposite sign in type $3$ and $5$, which
indicates that the $\mathbf{M}_{\mathrm{TM}}$ pointed originally below and
above the sample plane, i.e. $\theta_{\mathrm{TM}}\lessgtr\frac{\pi}{2}$,
respectively (where “above” points in the direction of the external magnetic
field). Naturally, an intermediate regime of dynamics (type $4$) exists which
is similar to $3$ and $5$, but without observable UD. We assigned type $4$
when the UD decreased below the arbitrarily chosen border of $<20~{}\%$ of the
original degree of UD (in type $1$ and $2$), in which case the degree of UD is
comparable to the experimental noise.

The field and temperature where each type of spin dynamics was observed (the
dynamical $H$-$T$ phase diagram) are shown in Fig. 3(c). All the raw curves
associated with each point are plotted in the Supplemental Material [22]. We
compared this dynamical $H$-$T$ phase diagram to the static diagram of Fig.
2(b), and discuss the conclusions that can be drawn from that. To aid this
discussion, we calculated several quantities from the dynamical measurements
shown in Fig. 4. The calculation of the initial polar angle of the FeCo
magnetization $\theta_{\mathrm{TM}}(t<0)$ in Fig. 4(a) was based on the
reduced out-of-plane projection of UD in the noncollinear phase as compared to
the collinear phase (see Supplemental Material [22] for details). Using these
initial angles, we determined the final state angle at $1.2$ ns as is depicted
in Fig. 4(b). In addition, we fitted the growth of the signal right after UD
with a sigmoidal function and extracted a timescale in which $\sim 96\%$ of
the growth takes place (see Supplementary information). The results are
depicted in Fig. 4(c).

By comparing the dynamical phase-diagram of Fig. 3(c) with the results of Fig.
4, we first conclude that the flip in the sign of the UD between types $1$a
and $1$b marks the compensation temperature $T_{\mathrm{M}}$, which has
shifted with about $5$ K compared to the static phase diagram likely due to
heat accumulation in the sample. Furthermore, for $T<T_{\mathrm{M}}$ and for
fields $0.15$ T $\leq B_{\mathrm{ext}}<B_{\mathrm{SF}}$ (type $2$), the spins
undergo a sub-ns magnetization reversal after the UD. The end-state after
$1.2$ ns is independent of the applied external field (see Fig. 4(b)), which
implies that FeCo undergoes a full reversal to phase II. An exception might be
the results just above the threshold field of $0.15$ T, seen as the yellow
line in Fig. 4(b). Here, the final state could correspond to a non-collinear
configuration, but it might also correspond simply to a partially reversed
spot. In Fig. 4(c), one can see that above the threshold field of $0.15$ T
(see dashed line), the timescale of reversal gradually reduces for larger
external fields from $\sim 1$ ns to about $100$ ps upon approaching phase III.
We estimated that the sample temperature is elevated by laser-heating to about
$140$ K after $1.2$ ns (see Supplemental Material [22]). Therefore, these
observations let us conclude that reversal of type $2$ occurs via transverse
motion after the system is laser-heated above the compensation point.

Next, we address the dynamics of types $3$-$5$, which are only triggered when
the system starts in the noncollinear state (phase III). Figure 4(b) shows
that the final angle in this entire region reaches
$\theta_{\mathrm{TM}}\sim\pi$, suggesting that after laser-induced heating the
system ends up in phase II. Most interestingly, the typical timescale of this
spin reorientation is $\sim 20$ ps, which is about an order of magnitude
faster than the spin reversal starting from the collinear phase I. This
ultrafast timescale is roughly independent of the applied magnetic field (see
Fig. 4(c)), indicating it is exchange-driven. One should note, however, that
despite the reorientation being fast, it is followed by spin precession around
the new equilibrium and the magnetization will eventually arrive at this new
equilibrium only after the oscillations are damped. The frequency of
precession $\nu_{\mathrm{SF}}$ is faster than expected for the ferromagnetic
resonance mode in a collinear state [24, 12], hence this mode should
correspond to the spin-flop resonance [25, 15]. The mode frequency follows a
slope defined by the effective $g$-factor, which is found to be $1.21\pm 0.11$
near $T_{M}$ (see Fig. 4(d)). Similar values of $g_{\mathrm{eff}}$ near
$T_{M}$ of the spin-flop resonance mode in GdFeCo were reported in high fields
[15].

Finally, the other key observation in our work regards the novel insight about
the tricritical point, as can be deduced from the dynamics in the purple
region (type $4$ in Fig. 3(c)), where the demagnetization is nearly invisible
or completely absent and therefore corresponds to the situation where the
$\mathbf{M}_{\mathrm{TM}}$ started in the plane
$\theta_{\mathrm{TM}}\approx\pi/2$. This region is seen to terminate at the
spin-flop transition, exactly where the tricritical point is expected. This
allows us to establish a novel empirical rule, which states that the
tricritical point separates the regions with
$\theta_{\mathrm{TM}}\lessgtr\frac{\pi}{2}$, i.e. when
$\mathbf{M}_{\mathrm{TM}}$ initially points towards different sides of the
anisotropy plane. This rule can be explained intuitively: when initially
$\theta_{\mathrm{TM}}<\frac{\pi}{2}$, $\mathbf{M}_{\mathrm{TM}}$ has to jump
over a potential barrier at $\theta_{\mathrm{TM}}=\frac{\pi}{2}$ due to
magnetic anisotropy to reach the TM-dominated collinear phase II (where
$\theta_{\mathrm{TM}}=\pi$). This jump is discontinuous and can thus be
classified as a first-order PT. Naturally, when $T<T_{\mathrm{M}}$ or when
$T>T_{\mathrm{3P}}$, $\mathbf{M}_{\mathrm{TM}}$ remains to point to the same
side of the plane during the spin-flop transition. Therefore it has no
potential barrier to pass, and the associated spin-flop PT can be continuous
and thus of the second order.

In conclusion, our results highlight the unique diversity and tunability of
magnetization dynamics in RE-TM ferrimagnets, and the established phase
diagram provides new insights into the role and origin of the so far illusive
tricritical point present in these materials. The extensive dataset provides
many parameters associated with the statics and dynamics of the noncollinear
phase, which could not be resolved before. We note that despite intense
studies of ultrafast magnetism in ferrimagnets during the last 15 years, the
role of the tricritical point in the diversity of the spin dynamics has not
been treated either theoretically or computationally. Finally, our results
demonstrate the possibility of switching between the noncollinear and
collinear spin configurations on ultrafast $\mathcal{O}(10~{}\text{ps})$
timescales. This puts the non-collinear phase forward as a potentially
interesting host of bits in magnetic recording.

###### Acknowledgements.

The authors thank K. Saaedi, C. Berkhout, and S. Semin for their technical
support. The work was supported by the European Research Council ERC Grant
Agreement No.101054664 273 (SPARTACUS), and the research program “Materials
for the Quantum Age” (QuMat). This program (registration number 024.005.006)
is part of the Gravitation program financed by the Dutch Ministry of
Education, Culture and Science (OCW).

## References

  * Beaurepaire  _et al._ [1996] E. Beaurepaire, J.-C. Merle, A. Daunois, and J.-Y. Bigot, Ultrafast Spin Dynamics in Ferromagnetic Nickel, Phys. Rev. Lett. 76, 4250 (1996).
  * Eschenlohr and Bovensiepen [2017] A. Eschenlohr and U. Bovensiepen, Special issue on ultrafast magnetism, Journal of Physics: Condensed Matter 30, 030301 (2017).
  * Kimel  _et al._ [2020] A. V. Kimel, A. M. Kalashnikova, A. Pogrebna, and A. K. Zvezdin, Fundamentals and perspectives of ultrafast photoferroic recording, Physics Reports 852, 1 (2020), fundamentals and perspectives of ultrafast photoferroic recording.
  * Stanciu  _et al._ [2007a] C. D. Stanciu, F. Hansteen, A. V. Kimel, A. Kirilyuk, A. Tsukamoto, A. Itoh, and T. Rasing, All-Optical Magnetic Recording with Circularly Polarized Light, Phys. Rev. Lett. 99, 047601 (2007a).
  * Kimel and Li [2019] A. V. Kimel and M. Li, Writing magnetic memory with ultrashort light pulses, Nat Rev Mater 4, 189 (2019).
  * Radu  _et al._ [2011] I. Radu, K. Vahaplar, C. Stamm, T. Kachel, N. Pontius, H. A. Dürr, T. A. Ostler, J. Barker, R. F. L. Evans, R. W. Chantrell, A. Tsukamoto, A. Itoh, A. Kirilyuk, T. Rasing, and A. V. Kimel, Transient ferromagnetic-like state mediating ultrafast reversal of antiferromagnetically coupled spins, Nature 472, 205 (2011).
  * Davies  _et al._ [2020] C. S. Davies, T. Janssen, J. H. Mentink, A. Tsukamoto, A. V. Kimel, A. F. G. van der Meer, A. Stupakiewicz, and A. Kirilyuk, Pathways for single-shot all-optical switching of magnetization in ferrimagnets, Phys. Rev. Appl. 13, 024064 (2020).
  * Yang  _et al._ [2017] Y. Yang, R. B. Wilson, J. Gorchon, C.-H. Lambert, S. Salahuddin, and J. Bokor, Ultrafast magnetization reversal by picosecond electrical pulses, Science Advances 3, e1603117 (2017).
  * Stupakiewicz  _et al._ [2017] A. Stupakiewicz, K. Szerenos, D. Afanasiev, A. Kirilyuk, and A. V. Kimel, Ultrafast nonthermal photo-magnetic recording in a transparent medium, Nature 542, 71 (2017).
  * Mangin  _et al._ [2014] S. Mangin, M. Gottwald, C.-H. Lambert, D. Steil, V. Uhlíř, L. Pang, M. Hehn, S. Alebrand, M. Cinchetti, G. Malinowski, Y. Fainman, M. Aeschlimann, and E. E. Fullerton, Engineered materials for all-optical helicity-dependent magnetic switching, Nat. Mater. 13, 286 (2014).
  * Kim  _et al._ [2022] S. K. Kim, G. S. D. Beach, K.-J. Lee, T. Ono, T. Rasing, and H. Yang, Ferrimagnetic spintronics, Nat. Mater. 21, 24 (2022).
  * Hintermayr  _et al._ [2023] J. Hintermayr, Y. L. W. van Hees, and B. Koopmans, Exploring terahertz-scale exchange resonances in synthetic ferrimagnets with ultrashort optically induced spin currents, Phys. Rev. B 108, 024401 (2023).
  * Zvezdin  _et al._ [1998] A. K. Zvezdin, I. A. Lubashevsky, R. Z. Levitin, G. M. Musaev, V. V. Platonov, and O. M. Tatsenko, Spin-flop and metamagnetic transitions in itinerant ferrimagnets, in  _Itinerant Electron Magnetism: Fluctuation Effects_, edited by D. Wagner, W. Brauneck, and A. Solontsov (Springer Netherlands, Dordrecht, 1998) pp. 285–302.
  * Clark and Callen [1968] A. E. Clark and E. Callen, Néel Ferrimagnets in Large Magnetic Fields, J. Appl. Phys. 39, 5972 (1968).
  * Becker  _et al._ [2017] J. Becker, A. Tsukamoto, A. Kirilyuk, J. C. Maan, T. Rasing, P. C. M. Christianen, and A. V. Kimel, Ultrafast Magnetism of a Ferrimagnet across the Spin-Flop Transition in High Magnetic Fields, Phys. Rev. Lett. 118, 117203 (2017).
  * Pogrebna  _et al._ [2019] A. Pogrebna, K. Prabhakara, M. Davydova, J. Becker, A. Tsukamoto, T. Rasing, A. Kirilyuk, A. K. Zvezdin, P. C. M. Christianen, and A. V. Kimel, High-field anomalies of equilibrium and ultrafast magnetism in rare-earth–transition-metal ferrimagnets, Phys. Rev. B 100, 174427 (2019).
  * Antropov  _et al._ [2021] N. O. Antropov, E. A. Kravtsov, M. V. Makarova, V. V. Proglyado, T. Keller, I. A. Subbotin, E. M. Pashaev, G. V. Prutskov, A. L. Vasiliev, Y. M. Chesnokov, N. G. Bebenin, M. A. Milyaev, V. V. Ustinov, B. Keimer, and Y. N. Khaydukov, Tunable spin-flop transition in artificial ferrimagnets, Phys. Rev. B 104, 054414 (2021).
  * Blank  _et al._ [2022] T. G. H. Blank, S. Hermanussen, T. Lichtenberg, T. Rasing, A. Kirilyuk, B. Koopmans, and A. V. Kimel, Laser-induced transient anisotropy and large amplitude magnetization dynamics in a Gd/FeCo multilayer, Adv. Mater. Interfaces 9, 2201283 (2022).
  * Davydova  _et al._ [2019] M. D. Davydova, K. A. Zvezdin, J. Becker, A. V. Kimel, and A. K. Zvezdin, $H\text{$-$}T$ phase diagram of rare-earth–transition-metal alloys in the vicinity of the compensation point, Phys. Rev. B 100, 064409 (2019).
  * Stanciu  _et al._ [2007b] C. D. Stanciu, A. Tsukamoto, A. V. Kimel, F. Hansteen, A. Kirilyuk, A. Itoh, and T. Rasing, Subpicosecond Magnetization Reversal across Ferrimagnetic Compensation Points, Phys. Rev. Lett. 99, 217204 (2007b).
  * Mekonnen  _et al._ [2013] A. Mekonnen, A. R. Khorsand, M. Cormier, A. V. Kimel, A. Kirilyuk, A. Hrabec, L. Ranno, A. Tsukamoto, A. Itoh, and T. Rasing, Role of the inter-sublattice exchange coupling in short-laser-pulse-induced demagnetization dynamics of GdCo and GdCoFe alloys, Phys. Rev. B 87, 180406 (2013).
  * Note [1] See Supplemental Material at [URL] for the full static magnetic characterization of the Gd/FeCo sample using magneto-optics and the SQUID-VSM, details of the optical pump-probe experimental method, all of the measured dynamical transients, estimations of the temperature increase due to laser-induced heating, measurements of UD with $100$ fs time-resolution and an example of fluence dependence of the data, which includes Refs. [26, 27].
  * Venkataraman [1979] G. Venkataraman, Soft modes and structural phase transitions, Bulletin of Materials Science 1, 129 (1979).
  * Stanciu  _et al._ [2006] C. D. Stanciu, A. V. Kimel, F. Hansteen, A. Tsukamoto, A. Itoh, A. Kirilyuk, and T. Rasing, Ultrafast spin dynamics across compensation points in ferrimagnetic $\mathrm{GdFeCo}$: The role of angular momentum compensation, Phys. Rev. B 73, 220402 (2006).
  * Gurevich and Melkov [1996] A. Gurevich and G. Melkov,  _Magnetization Oscillations and Waves_ (Taylor & Francis, 1996).
  * Kools  _et al._ [2023] T. J. Kools, Y. L. W. van Hees, K. Poissonnier, P. Li, B. Barcones Campo, M. A. Verheijen, B. Koopmans, and R. Lavrijsen, Aging and passivation of magnetic properties in Co/Gd bilayers, Applied Physics Letters 123, 042406 (2023).
  * Ostler  _et al._ [2011] T. A. Ostler, R. F. L. Evans, R. W. Chantrell, U. Atxitia, O. Chubykalo-Fesenko, I. Radu, R. Abrudan, F. Radu, A. Tsukamoto, A. Itoh, A. Kirilyuk, T. Rasing, and A. Kimel, Crystallographically amorphous ferrimagnetic alloys: Comparing a localized atomistic spin model with experiments, Phys. Rev. B 84, 024407 (2011).