# Refining magnetic interactions from the magnetic field dependence of spin-
wave excitations in magnetoelectric LiFePO4

L. Peedu  National Institute of Chemical Physics and Biophysics, Akadeemia tee
23, 12618 Tallinn, Estonia  V. Kocsis  RIKEN Center for Emergent Matter
Science (CEMS), Wako, Saitama 351-0198, Japan  D. Szaller  Institute of Solid
State Physics, TU Wien, 1040 Vienna, Austria  B. Forrai  Department of
Physics, Institute of Physics, Budapest University of Technology and
Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary  S. Bordács  Department
of Physics, Institute of Physics, Budapest University of Technology and
Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary  I. Kézsmárki
Department of Physics, Institute of Physics, Budapest University of Technology
and Economics, Műegyetem rkp. 3., H-1111 Budapest, Hungary  Experimental
Physics 5, Center for Electronic Correlations and Magnetism,Institute of
Physics, University of Augsburg, 86159 Augsburg, Germany  J. Viirok  National
Institute of Chemical Physics and Biophysics, Akadeemia tee 23, 12618 Tallinn,
Estonia  U. Nagel  National Institute of Chemical Physics and Biophysics,
Akadeemia tee 23, 12618 Tallinn, Estonia  B. Bernáth  D. L. Kamenskyi  High
Field Magnet Laboratory (HFML-EMFL), Radboud University, Toernooiveld 7, 6525
ED Nijmegen, The Netherlands  A. Miyata  O. Portugall  Laboratoire National
des Champs Magnétiques Intenses (LNCMI-EMFL), CNRS-UGA-UT3-INSA,143 Avenue de
Rangueil, 31400 Toulouse, France  Y. Tokunaga  RIKEN Center for Emergent
Matter Science (CEMS), Wako, Saitama 351-0198, Japan  Department of Advanced
Materials Science, University of Tokyo, Kashiwa 277-8561, Japan  Y. Tokura
RIKEN Center for Emergent Matter Science (CEMS), Wako, Saitama 351-0198, Japan
Department of Applied Physics, University of Tokyo, Hongo, Tokyo 113-8656,
Japan  Y. Taguchi  RIKEN Center for Emergent Matter Science (CEMS), Wako,
Saitama 351-0198, Japan  T. Rõõm  National Institute of Chemical Physics and
Biophysics, Akadeemia tee 23, 12618 Tallinn, Estonia

###### Abstract

We investigated the spin excitations of magnetoelectric LiFePO4 by THz
absorption spectroscopy in magnetic fields up to 33 T. By studying their
selection rules, we found not only magnetic-dipole, but also electric-dipole
active (electromagnons) and magnetoelectric resonances. The magnetic field
dependence of four strong low-energy modes is reproduced well by our four-
sublattice spin model for fields applied along the three orthorhombic axes.
From the fit, we refined the exchange couplings, single-ion anisotropies, and
the Dzyaloshinskii-Moriya interaction parameters. Additional spin excitations
not described by the mean-field model are observed at higher frequencies. Some
of them shows a strong shift with magnetic field, up to 4
$\,\mathrm{cm}^{-1}$/T, when the field is applied along the easy axis. Based
on this field dependence, we attribute these high frequency resonances to
excitation of higher spin multipoles and of two magnons, which become THz-
active due to the low symmetry of the magnetically ordered state.

##  1 Introduction

Recent optical studies of multiferroic materials have revealed non-reciprocal
directional dicroism, which is the light absorption difference for unpolarized
counter-propagating beams [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
16, 17, 18, 19, 20, 21, 22, 23, 24, 25]. This unusual phenomenon is the
finite-frequency manifestation of the magnetoelectric (ME) effect, which
emerges at simultaneously electric- and magnetic-dipole allowed excitations,
that we term as ME resonance111Usually, magnons couple to the magnetic
component of electromagnetic radiation, i.e. they are magnetic-dipole active.
If the magnons are electric-dipole active, the term electromagnon is often
used [26]. Magnetoelectric resonance is a spin-wave excited by both components
of electromagnetic radiation [3, 27]. For the rest of the paper we classify
the spin-waves, based on their coupling to the electromagnetic radiation,
using magnetic-dipole active, electric-dipole active and magnetoelectric spin-
wave. We use “magnon” for the spin-wave excitation described by the mean-field
model without making a difference in its coupling to the electromagnetic
radiation.. Since the relative orientation of the electric and magnetic fields
is different for counter-propagating beams ME coupling generates an absorption
difference between the two beams and may even lead to one-way transparency
[7]. This non-reciprocal absorption may gain applications in photonics and
spintronics [3, 13]. For example, materials with ME resonances can be used as
optical diodes where the direction of transparency for the THz radiation can
be switched by magnetic fields [4, 3, 5, 7, 9, 8, 17], electric fields [28,
20, 24], or both [13]. From the fundamental science point of view, the
spectroscopy of the ME excitations promotes the understanding of the static ME
response which is linked to the non-reciprocal directional dichroism spectrum
via the Kramers-Kronig relations [29, 16]. Moreover, a THz absorption study,
combined with magnetization, inelastic neutron scattering measurements [30,
31, 32], and theoretical modeling [33, 34, 35, 36] can resolve realistic spin
Hamiltonians of ME compounds.

The relativistic spin-orbit coupling plays an essential role for ME spin
excitations. It establishes a coupling between spins and electric dipoles and
also introduces single-ion anisotropy for $S>1/2$. The single-ion anisotropy
expands the frequency scale of spin excitations as it separates the $\pm
m_{s}$ doublets from each other in zero field, where $m_{s}$ is the spin
quantum number. In addition to conventional spin waves, spin-quadrupolar
excitations corresponding to $\Delta m_{s}=\pm 2$ may appear in systems with
strong single-ion anisotropy and spin $S>1/2$ [37, 38, 39, 40, 41], broadening
the frequency range for possible applications of ME materials. In general, if
there are $N$ spins in the magnetic unit cell we expect $2NS$ spin
excitations, which can be described by the multi-boson spin-wave theory [37,
38, 13, 42] or by single-ion spin Hamiltonian with added molecular field to
take into account spin-spin interactions [43, 32, 44].

The Li$M$PO4 ($M=$ Mn, Co, Fe, Ni) orthophosphate compounds become ME as their
magnetic order breaks the inversion symmetry [45]. This, together with their
large single-ion anisotropy [31, 32, 46], makes them appealing candidates to
explore unconventional spin excitations. Among these compounds, LiFePO4 has
the highest Néel temperature, $T_{\mathrm{N}}$ $=50$ K below which an
antiferromagnetic (AFM) order develops, as depicted in Fig. 1. The spins of
the four magnetic ions of the unit cell are nearly parallel to the $y$ axis
[47]. Detailed neutron diffraction experiments showed that the spins are
slightly rotated in the $xy$ plane and canted toward the $z$ axis [31].
LiFePO4 has one of the largest spins in the orthophosphate family but the
number of spin-wave modes detected in the magnetically ordered phase has been
substantially less than $2NS=16$, allowed for a $S=2$ spin system. In zero-
field inelastic neutron scattering (INS) studies two spin-wave branches [30,
31, 32] and a dispersionless mode were observed below 10 meV [32]. Whereas, a
recent high-frequency electron spin resonance study detected two modes in the
vicinity of the spin-flop field, 32 T [46]. Therefore, further experimental
data is needed to understand better the spin dynamics and spin Hamiltonian of
LiFePO4.

In this work, we studied the magnetic field dependence of the spin excitations
using THz absorption spectroscopy in the low temperature AFM phase of LiFePO4.
The spectral range of our experiments extending up to 175 $\,\mathrm{cm}^{-1}$
(22 meV) covers two and five times larger energy window compared to former INS
[30, 31, 32] and electron spin resonance studies [48, 46], respectively. The
broader spectral range allowed us to observe 17 spin excitations and to
determine their selection rules. The absorption spectra were measured with
magnetic field applied along all three principal crystallographic axes up to
33 T in the Faraday configuration (light propagates along the field,
$\mathbf{k}\\!\parallel\\!\mathbf{H}$) and up to 17 T in Voigt geometry (light
propagates perpendicular to the field, $\mathbf{k}\\!\perp\\!\mathbf{H}$).
Beside THz spectroscopy, we measured high-field magnetization up to 120 T
along the easy-axis from which we determined the spin-flop and the saturation
fields. Finally, we successfully employed a mean-field model to describe the
field dependence of the magnetization and the resonance frequencies of the
four strongest low-frequency spin-wave modes in the AFM state. By fitting the
magnetic field dependence of four magnons we have refined the values of the
exchange couplings, the single-ion anisotropies, and the Dzyaloshinskii-Moriya
interaction.

##  2 Experimental

The LiFePO4 single crystals were grown by the floating zone method [49]. The
quality of the crystals was verified by powder diffraction and Laue XRD, which
confirmed the orthorhombic structure with the same lattice constants as
reported in Ref. [50].

The low field magnetization measurements were done using a 14 T PPMS with VSM
option (Quantum Design). High-field magnetization measurements were carried
out up to 120 T using ultra-high semidestructive pulses at the Laboratoire
National des Champs Magnétiques Intenses in Toulouse [51, 52]. The maximum
field of a semidestructive pulse was reached in $\sim$2.5 $\mu$s.

For THz spectroscopy studies the single crystal was cut into three 1 mm thick
slabs each with a large face normal to one of the three principal
crystallographic directions. The slabs were wedged by two degrees to suppress
the fringes in the spectra produced by the internal reflections in the
crystal.

The THz measurements up to 17 T were performed with a polarizing Martin-
Puplett interferometer and a 0.3 K silicon bolometer in Tallinn. High field
spectra from 17 T up to 33 T were measured using a Bruker IFS 113v
spectrometer and a 1.6 K silicon bolometer in High Field Magnet Laboratory in
Nijmegen. The experiments above 17 T were done in Faraday configuration, while
below 17 T both Faraday and Voigt configuration experiments were performed.
All spectra were measured with an apodized spectral resolution of 0.3 or 0.5
$\,\mathrm{cm}^{-1}$. The polarizer angle with respect to the crystal axes was
determined by evaluating the intensity change of the strongest modes in the
THz absorption spectra as the function of rotation angle of the polarizer.
This information was also used to mount the polarizer in the high field
experiments in Nijmegen where the in situ polarizer rotation was not possible.
Absorption was determined by using a reference spectrum of an open hole,
sample spectrum in the paramagnetic phase or by statistically calculating the
baseline from the magnetic field dependence of sample spectra. In the first
method, the absorption was calculated as

$\alpha=-d^{-1}\ln(I/I_{r}),$ (1)

where $I_{r}$ is the intensity through the reference hole with the area equal
to the sample hole area and $d$ is the sample thickness. In the second method,
the absorption difference was calculated,

$\displaystyle\Delta\alpha(H,T)$ $\displaystyle=$
$\displaystyle\alpha(H,T)-\alpha(0\,\text{T},55\,\text{K})$ (2)
$\displaystyle=$
$\displaystyle-d^{-1}\ln\left[I(H,\textit{T})/I(0\,\text{T},55\,\text{K})\right],$

where $I(0\,\text{T},55\,\text{K})$ is the intensity through the sample
measured at 0 T and 55 K in the paramagnetic phase. In the third method, the
statistically calculated baseline, $\alpha(0\,\text{T})$, was found as a
minimum of differential absorption,

$\displaystyle\Delta\alpha_{H}(H_{i})$ $\displaystyle=$
$\displaystyle\alpha(H_{i})-\alpha(0\,\text{T})$ (3) $\displaystyle=$
$\displaystyle-d^{-1}\ln\left[I(H_{i})/I(0\,\text{T}\right],$

at each frequency over several magnetic field values $H_{i}$. By adding
$\alpha(0\,\text{T})$ to the differential absorption we get the dependence of
absorption spectra on magnetic field. This method was used to obtain the
spectra measured above 17 T.

Figure 1:  The ground state spin configuration of LiFePO4 in zero magnetic
field. There are four Fe2+ spins, $S=2$, in the magnetic unit cell drawn as a
box. $S_{1}$ is in the corner, $S_{2}$ is on the face and $S_{3}$ and $S_{4}$
are inside the unit cell. The spins are rotated towards the $x$ axis and
canted towards the $z$ axis away from the $y$ axis [31]. The numbering of
spins and the labeling of exchange interactions corresponds to the spin
Hamiltonian described by Eq. (4).  Figure 2:  (a) Magnetic field dependence of
the magnetization $M$ (solid, red) and (b) the directly measured
$\mathrm{d}M/\mathrm{d}H$ (solid, red) at $T=5\,\mathrm{K}$ for increasing
pulsed field in H$\,\parallel\,$y. The inset of panel (a) shows VSM $M$-$H$
measurements in quasi-static fields at $T=2.4\,\mathrm{K}$, where the magnetic
field directions are H$\,\parallel\,$x (black), H$\,\parallel\,$y (red) and
H$\,\parallel\,$z (green). The dashed lines in panel (a) show the results of
the mean-field calculations with the parameters from Table 1. For comparison,
we show the pulsed field magnetization data from Ref. [46] with open red
circles. The AFM, spin-flopped and spin polarized state regions are shown for
H$\,\parallel\,$y.  Figure 3:  (a) THz absorption spectra of LiFePO4 at $55$ K
in the paramagnetic phase, and (b) the difference between the zero-field
absorption spectra recorded at $3.5$ K (magnetically ordered phase) and 55 K,
demonstrating spectral features associated with the onset of magnetic order.
Line colors correspond to the propagation direction of the THz radiation:
$k_{x}$ (blue), $k_{y}$ (green) and $k_{z}$ (red). Two orthogonal
polarizations $\\{E^{\omega}_{i},H^{\omega}_{j}\\}$ for the given propagation
direction, $k_{k}\sim E^{\omega}_{i}\times H^{\omega}_{j}$, are indicated by
the solid and dashed lines, according to the inset of panel (b). $F_{n}$ with
$n=1,\ldots,17$ labels the modes in the magnetically ordered phase and
$F_{\mathrm{os}}$ is an on-site magnetic excitation in the paramagnetic phase.
$H^{\omega}_{j}$ or $E^{\omega}_{i}$ indicate the magnetic- or electric-dipole
activity of the mode, respectively. The blue, green and red rectangles mark
the peaks with absorption above the upper detection limit, $F_{4}$, $F_{6}$,
$F_{11}$, and phonon.

##  3 Mean-field model

The mean-field theory of localized magnetic moments is a widely applied tool
to interpret the static and dynamic magnetic properties of systems with
periodic magnetic structures [53], e.g. ferro- [54], ferri- [34] and
antiferromagnetic- [33] insulators. Particularly, the microscopic spin
Hamiltonian of LiFePO4 has been discussed by several papers [30, 55, 31, 32,
56, 46].

Here we aim to describe the static magnetism and the infrared-active optical
magnetic resonances of LiFePO4. Thus, we use a simplified Hamiltonian where
the exchange coupling terms $J_{y}$ and $J_{z}$ have been omitted as they
connect spins at crystallographically equivalent sites, see Fig. 1. While the
one-magnon THz spectrum is insensitive to the same energy shift of all states
at the $\Gamma$ point of the Brillouin zone produced by $J_{y}$ and $J_{z}$,
these couplings are relevant when describing the dispersion of the magnon
modes [31]. The Hamiltonian of our study is

$\displaystyle\mathcal{H}$ $\displaystyle=$ $\displaystyle
4\left[J_{xz}\left(\mathbf{S}_{1}\cdot\mathbf{S}_{2}+\mathbf{S}_{3}\cdot\mathbf{S}_{4}\right)\right.$
(4)
$\displaystyle+J_{xy}\left(\mathbf{S}_{1}\cdot\mathbf{S}_{3}+\mathbf{S}_{2}\cdot\mathbf{S}_{4}\right)$
$\displaystyle+J_{yz}\left(\mathbf{S}_{1}\cdot\mathbf{S}_{4}+\mathbf{S}_{2}\cdot\mathbf{S}_{3}\right)$
$\displaystyle+\left.D_{14}\left(S_{1}^{y}S_{4}^{z}-S_{1}^{z}S_{4}^{y}+S_{3}^{y}S_{2}^{z}-S_{3}^{z}S_{2}^{y}\right)\right]$
$\displaystyle+$
$\displaystyle\sum_{i=1}^{4}\left[\Lambda_{x}\left(\vphantom{S_{i}^{y}}S_{i}^{x}\right)^{2}+\Lambda_{z}\left(S_{i}^{z}\right)^{2}+\Lambda_{xy}\left(S_{i}^{x}S_{i}^{y}\right)\right.$
$\displaystyle-\left.\mu_{\mathrm{B}}\mu_{0}\left(g_{x}\mathit{H}_{x}\mathit{S}^{x}_{i}+g_{y}\mathit{H}_{y}\mathit{S}^{y}_{i}+g_{z}\mathit{H}_{z}\mathit{S}^{z}_{i}\right)\right],$

where the terms, exchange interactions, Dzyaloshinskii-Moriya term, single-ion
anisotropy terms, and the Zeeman energy, have been considered in the earlier
works of LiFePO4 summarized in Table 1. The model is based on four Fe2+ spins,
here represented by classical vectors of $S\,=\,2$ length, that occupy
crystallographically non-equivalent positions of the unit cell. As shown in
Fig.1, the spins are coupled by three different exchange couplings with
parameters $J_{xz},J_{xy}$ and $J_{yz}$. There are two single-site hard-axis
anisotropies, $\Lambda_{x}$ and $\Lambda_{z}$, that effectively produce the
easy axis along $y$. The spins are slightly rotated away from the $y$ axis
towards the $x$ axis as observed by neutron scattering [31]. Extending
previous studies to reproduce this small deviation of the magnetic structure
from the collinear antiferromagnetic order we introduced an additional single-
site anisotropy term $\Lambda_{xy}S^{x}S^{y}$. The Dzyaloshinskii-Moriya
interaction is $\mathbf{D}=\\{D_{14},0,0\\}$. Since the spins are
predominantly along the $y$ direction, $D_{14}$ cants spins towards the $z$
direction. The last term in Eq. (4) is the interaction of the electron spin
with the applied magnetic field taking into account the $g$-factor anisotropy.

We model the spin dynamics using the Landau-Lifshitz-Gilbert equation, Ref.
[57], as used in Ref. [58], by assuming that the spins are oscillating about
their equilibrium orientations without changing their lengths. The equilibrium
orientation of the spins is found by minimizing the energy described in the
Hamiltonian of Eq. (4), with respect to the spin orientations. The magnetic-
dipole absorption of light by spin-waves is calculated assuming that the
magnetic field $\mathbf{H}^{\omega}$ of radiation couples to the total
magnetic moment of the spins [58]. Dielectric permittivity in the absorption
coefficient formula, Eq. (10) in Ref. [58], was assumed to be real and
frequency-independent with components $\epsilon_{x}=8.1$, $\epsilon_{y}=7.3$
and $\epsilon_{z}=7.6$ [59].

Table 1: The parameters of the mean-field model used to describe the static magnetic properties and spin-waves in LiFePO4: exchange couplings $J_{i}$ and $J_{ij}$, single-ion anisotropies $\Lambda_{i}$ and $\Lambda_{ij}$, Dzyaloshinskii-Moriya coupling $D_{14}$, and anisotropic $g$-factor $g_{i}$. All parameters are in units of meV except the dimensionless $g_{i}$.  $J_{xz}$ | $J_{xy}$ | $J_{yz}$ | $\Lambda_{x}$ | $\Lambda_{z}$ | $\Lambda_{xy}$ | $D_{14}$ | $g$ | Ref.
---|---|---|---|---|---|---|---|---
$-0.006$ | 0.086 | 0.51 | 0.52 | 1.52 | $-0.009$ | 0.027 | $g_{x}=2.04$ | 222This work
|  |  |  |  |  |  | $g_{y}=2.3$ |
|  |  |  |  |  |  | $g_{z}=2.1$ |
0.05 | 0.14 | 0.77 | 0.62 | 1.56 | - | 0.038333The Dzyaloshinskii-Moriya parameter $D_{14}=J_{\mathrm{DM}}/4$ where $J_{\mathrm{DM}}$ is from Ref. [46]. | $g_{x}=2.24$ | [46]
|  |  |  |  |  |  | $g_{y}=2.31$ |
|  |  |  |  |  |  | $g_{z}=1.99$ |
0.01 | 0.09 | 0.46 | 0.86 | 2.23 | - | - | - | [32]
0.05 | 0.14 | 0.77 | 0.62 | 1.56 | - | - | - | [31]

##  4 Results

###  4.1 Magnetization

We characterized LiFePO4 samples by measuring the magnetization at 2.4 K along
the principal axes up to 14 T. Along H$\,\parallel\,$y, the measurement is
extended up to 120 T at 5 K using pulsed fields, see Fig. 2. The $y$-axis
magnetization determined from the pulsed-field measurements was normalized to
the value of static field measurements in the range from 4 to 14 T, neglecting
a small hysteresis of magnetization between 0 and 4 T. In the AFM state the
spins are predominantly aligned along the easy axis, the $y$ axis in LiFePO4.
The magnetization grows approximately linearly in increasing field for
H$\,\parallel\,$x and H$\,\parallel\,$z. When H$\,\parallel\,$y is applied,
the spins maintain easy-axis alignment up to the spin-flop field marked by a
jump in the magnetization at 32$\pm$3 T. As the field further increases the
magnetization grows linearly and reaches saturation at 56$\pm$3 T. In the
field-polarized state, the saturation magnetization is estimated to 4.4
$\mu_{B}$ per iron. The spin-flop field deduced from our measurements is in
agreement with former high-field magnetization measurements [46].

###  4.2 THz absorption spectra in zero field

The zero-field THz absorption spectra of LiFePO4 are presented in Fig. 3 and
the mode parameters are collected in Table 2, while Fig. 3(b) features
absorption spectra in the AFM phase, relative to the paramagnetic phase.

The spectra in the paramagnetic phase show a broad but weak magnetic-dipole
active peak $F_{\mathrm{os}}$ at around $\sim$55 $\,\mathrm{cm}^{-1}$, Fig. 3
(a). The magnetic on-site excitation $F_{\mathrm{os}}$ is
$\mathit{H}^{\omega}_{x}$ active as it is seen in two polarization
configurations, {$E^{\omega}_{y}$,$\mathit{H}^{\omega}_{x}$} and
{$E^{\omega}_{z}$,$\mathit{H}^{\omega}_{x}$}. The frequency and the selection
rules of $F_{\mathrm{os}}$ are reproduced by exact diagonalization of a four-
spin cluster, see Fig. S6 in the Supplementary Material. Other features in the
paramagnetic phase spectra are $E^{\omega}_{x}$-active phonon at
140$\,\mathrm{cm}^{-1}$, with a strong absorption exceeding the detection
limit, and absorption rising towards higher frequencies due to the phonons
with resonance frequencies above 175$\,\mathrm{cm}^{-1}$.

To better resolve spectral features emerging in the magnetically ordered phase
we plot the difference spectra,
$\alpha(3.5\,\mathrm{K})-\alpha(55\,\mathrm{K})$, Fig. 3 (b). We observe a
diminished absorption at the tails of phonons at low $T$ appearing as negative
features in the difference spectra between 140 and 175$\,\mathrm{cm}^{-1}$.
The change of the 140$\,\mathrm{cm}^{-1}$ phonon mode is obscured by the
strong absorption and therefore the $E^{\omega}_{x}$-spectra, green and red
solid lines, are discontinued where the 140$\,\mathrm{cm}^{-1}$ phonon peaks.
The broad peak $F_{\mathrm{os}}$ from the high-$T$ paramagnetic phase appears
as a negative feature in the difference spectra in $\mathit{H}^{\omega}_{x}$
polarization.

All sharp modes, labeled $F_{1}$,$\ldots$, $F_{17}$, are absent above
$T_{\mathrm{N}}$ and we assign them to spin excitations. The seven
excitations, $F_{3}$, $F_{4}$, $F_{5}$, $F_{6}$, $F_{8}$, $F_{9}$ and
$F_{11}$, are identified as magnetic-dipole active modes. Six modes, $F_{10}$,
$F_{12}$, $F_{13}$, $F_{15}$, $F_{16}$ and $F_{17}$, are identified as
electric-dipole active resonances. The mode $F_{13}$ has a shoulder, thus, it
was fitted with two Gaussian lines with maxima at 124.4$\,\mathrm{cm}^{-1}$
and 127.6$\,\mathrm{cm}^{-1}$. Two modes, $F_{7}$ at 71.4$\,\mathrm{cm}^{-1}$
and $F_{14}$ at 137.1$\,\mathrm{cm}^{-1}$, are both electric- and magnetic-
dipole allowed, therefore, we identified them as ME resonances. $F_{7}$ is the
strongest in {$E^{\omega}_{y}$,$\mathit{H}^{\omega}_{x}$} polarization, red
dashed line in Fig. 3 (b), and its intensity is halved if only one of the
components, $E^{\omega}_{y}$ or $\mathit{H}^{\omega}_{x}$, is present. Thus,
$F_{7}$ is an example of a ME resonance which couples equally to the magnetic
and electric components of radiation. We detected $F_{14}$ in the same three
polarization configuration, thus, we also assigned it to a ME resonance with
the same selection rule as mode $F_{7}$,
{$E^{\omega}_{y}$,$\mathit{H}^{\omega}_{x}$}.

The three strongest magnetic-dipole active modes $F_{4}$, $F_{6}$, and
$F_{11}$ show only weak absorption in polarizations orthogonal to their main
magnetic dipole component. The weak absorption in other polarizations could be
explained by the imperfections of the polarizer. However, we can not
completely rule out that some of these modes are ME resonances with a weak
electric-dipole component which can be tested by further measurements of the
non-reciprocal directional dichroism on magneto-electrically poled samples
[13, 16]. We can not identify the selection rules for modes $F_{1}$ and
$F_{2}$ as they are too weak.

###  4.3 Magnetic field dependence of spin-waves

The magnetic field dependence of mode frequencies and intensities between 0
and 17 T is shown in Fig. 4 for Faraday, panels (a)-(c), and Voigt
configuration, (d) and (e). The modes mostly stay at constant frequency when
the magnetic field is applied along the hard axes, H$\,\parallel\,$x, Fig.
4(a) and H$\,\parallel\,$z, Fig. 4(c, e). However, most of the resonances
shift with the magnetic field for H$\,\parallel\,$y. We assigned a slope,
$b_{1}=\Delta E/\Delta B$, calculated between 15 and 17 T in units
$\,\mathrm{cm}^{-1}$T-1, to each of the modes and collected them in Table 2.
If the mode was not visible in this range, a lower magnetic field range was
used. One mode, $F_{17}$, has zero slope and $F_{9}$, $F_{13}$, and $F_{16}$
have a moderate value, $|b_{1}|<0.3$. Modes $F_{14}$ and $F_{15}$ have the
largest $|b_{1}|$ for H$\,\parallel\,$y but also a substantial $|b_{1}|$ for
H$\,\parallel\,$z.

Assuming $g\approx 2$ we estimated from the slopes the change of the spin
projection quantum number, $\Delta m_{s}$, upon the excitations. The results
are listed in Table 2. The spin-waves below 80$\,\mathrm{cm}^{-1}$ (zero-field
frequency) have $|\Delta m_{s}|=1$ while above 100$\,\mathrm{cm}^{-1}$
$|\Delta m_{s}|$ is 2, 3 or 4. $|\Delta m_{s}|$ was not assigned to $F_{1}$
and $F_{2}$ where $b_{1}\approx 1.5$$\,\mathrm{cm}^{-1}$T-1 below 8 T, which
is between $\Delta m_{s}=1$ and 2. We note that $b_{1}$ of $F_{1}$ changes
with field. It is 0.9$\,\mathrm{cm}^{-1}$T-1 above 8 T. This change of slope
could be due to the anti-crossing with $F_{4}$ but we do not have evidence for
that because the mode was too weak to be detected in the high-field magnet
set-up above 17 T.

The absorption spectra in high magnetic field H$\,\parallel\,$y up to 31.6 T
are presented in Fig. 5. The spin-wave excitations, $F_{6}$ and $F_{7}$, start
softening before reaching the spin-flop transition at 32 T, in accordance with
the mean-field calculation. Also, $F_{13}$ at about 125$\,\mathrm{cm}^{-1}$
shows softening close to 30 T. Spectra in other two field directions,
H$\,\parallel\,$x and H$\,\parallel\,$z, above 17 T are shown in Supplementary
Material, Figs. S1 and S3.

Table 2:  The excitation configurations and field dependence of LiFePO4 modes in the AFM phase. The selection rules were found by measuring polarization dependence of spin excitations in three principal directions without magnetic field. The absorption line energy and area in zero field were obtained from the fit to Gaussian lineshape, except $F_{13}$ where the sum of two Gaussians was used. The slopes of the modes were estimated from the linear field dependence between 15 and 17 T; if mode was not visible in this field range, the lower field range was used. From the slopes the $|\Delta m_{s}|$ values are proposed assuming $g\approx 2$. Modes $F_{4}$ to $F_{7}$ were observed by INS spectroscopy [31, 32] and are fitted to the mean-field model in this work.  Mode | Energy at 0 T | Area at 0 T | Selection | Magnetic field | Slope $b_{1}$ | $|\Delta m_{s}|$
---|---|---|---|---|---|---
| (cm-1) | (cm-2) | rules at 0 T | direction | (cm-1T-1) |
$F_{1}$ | $18.3$ | 4 |  | $z$ | $+1.4$ |
$F_{2}$ | $24.7$ | 2 |  | $z$ | $+1.5$ |
$F_{3}$ | $30.8$ | 2 | $\mathit{H}^{\omega}_{z}$ | $y$ |  $-0.9$, $+0.9$ | 1
$F_{4}$ | $46.2\,\,(5.7\,\mathrm{meV})$ |  $>$100 | $\mathit{H}^{\omega}_{z}$ | $y$ | $-1.1$ | $1$
$F_{5}$ | $58.0\,\,(7.2\,\mathrm{meV})$ | 6 | $\mathit{H}^{\omega}_{z}$ | $y$ | $-1.1$ | $1$
$F_{6}$ | $67.9\,\,(8.4\,\mathrm{meV})$ |  $>$200 | $\mathit{H}^{\omega}_{x}$ | $y$ | $+0.9$ | $1$
$F_{7}$ | $71.4\,\,(8.9\,\mathrm{meV})$ | 37 |  $\mathit{H}^{\omega}_{x}$, $E^{\omega}_{y}$ | $y$ | $+1.0$ | 1
$F_{8}$ | $76.2$ | 9 | $\mathit{H}^{\omega}_{x}$ | $y$ |  $-0.8$, $+1.0$ | 1
$F_{9}$ | $90.8$ | 2 | $\mathit{H}^{\omega}_{z}$ | $x$ | $+0.1$ |
$F_{10}$ | $102.2$ | 57 | $E^{\omega}_{z}$ | $y$ | $-3.3$ | 3
$F_{11}$ | $109.0$ | 74 | $\mathit{H}^{\omega}_{y}$ | $y$ | $+1.8$ | $2$
$F_{12}$ | $120.8$ | 50 | $E^{\omega}_{y}$ | $y$ | $-1.9$ | $2$
$F_{13}$ | 124.4, 127.6 | 185 | $E^{\omega}_{x}$ | $y$ | $-0.3$ |
$F_{14}$ | $137.1$ | 17 |  $\mathit{H}^{\omega}_{x}$, $E^{\omega}_{y}$ | $y$ |  $-3.0$, $+2.8$ | $3$
|  |  |  | $z$ | $-0.6$ |
$F_{15}$ | $146.3$ | 30 | $E^{\omega}_{z}$ | $y$ |  $-3.7$, $+3.8$ | $4$
|  |  |  | $z$ | $+0.7$ |
$F_{16}$ | $163.7$ | 2 | $E^{\omega}_{y}$ | $x$ | $-0.3$ |
$F_{17}$ | $164.8$ | 4 | $E^{\omega}_{x}$ | $y$ | $0.0$ |
Figure 4:  Magnetic field dependence of the spin-wave resonance frequencies
and absorption line areas at $T=3.5\,\mathrm{K}$ in LiFePO4. Panels (a), (b),
and (c) correspond to measurements in the Faraday (${\bf
k}\parallel\mathbf{H}$), while panels (d) and (e) correspond to experiments in
the Voigt (${\bf k}\perp\mathbf{H}$) configuration. The direction of the
magnetic field is (a) – H$\,\parallel\,$x, (b), (d) – H$\,\parallel\,$y, and
(c), (e) – H$\,\parallel\,$z. The symbols correspond to six combinations of
linear light polarization $\\{E^{\omega}_{i},H^{\omega}_{j}\\}$ as indicated
at bottom left of the figure. The symbol height is proportional to the square
root of experimental absorption line area with the same scaling as wavenumber
axis. To simplify the figure the larger symbols are not shown for every
measured field. The error bars (vertical green lines) from fitting the line
positions in most cases are too small to be seen in the figure. The black
lines are the results of the mean-field model calculations, modes $F_{4}$,
$F_{5}$, $F_{6}$, and $F_{7}$. Comparison of experimental and calculated
intensities is in the Supplementary Material.  Figure 5:  THz absorption
spectra of spin-wave excitations in magnetic field
$\mathbf{H}\parallel\mathbf{k}\parallel\mathbf{y}$ at $T=3.5$ K in two
orthogonal polarizations, (a) $\\{\mathit{H}^{\omega}_{z},E^{\omega}_{x}\\}$
and (b) $\\{\mathit{H}^{\omega}_{x},E^{\omega}_{z}\\}$. The broad stripes are
the results of the mean-field calculation, $F_{4}$, $F_{5}$, $F_{6}$ and
$F_{7}$, in order of increasing frequency; the width of the line is
proportional to the square root of the line area with the same scale as
wavenumber axis and calculated in the magnetic dipole approximation. The zero
field 55 K spectrum was used as a reference below 17 T and the low temperature
zero-field spectrum above 17 T. Filled circles are the spin-wave excitation
energies reproduced from Ref. 46.

###  4.4 Mean-field model results

The mean-field model parameters of Table 1 were obtained by fitting the
magnetic field dependence of frequencies of the spin-wave modes $F_{4}$,
$F_{5}$, $F_{6}$, and $F_{7}$. The magnetic field dependence of these modes
below 17 T is reproduced remarkably well by the model for all three magnetic
field directions, Fig. 5 and Fig. 4. The isotropic $g$-factor was not
sufficient to quantitatively describe the magnetic field dependence of mode
frequencies. The anisotropic $g$-factor values improved not only the magnetic
field dependence of spin-wave frequencies, but also reproduce the value of the
spin-flop field and the saturation field, Fig. 2. In addition, the calculated
magnetization as a function of $H_{x}$ and $H_{z}$ follows the measured $M(H)$
below 15 T, inset to Fig. 2.

The equilibrium spin configuration deviates in two ways from the perfect
collinear arrangement of spins along the $y$ axis. With the parameters from
Table 1 we get the canting away from the $y$ axis towards the $z$ axis by 0.86
degrees, driven by the Dzyaloshinskii-Moriya interaction $D_{14}$. The
rotation of spins, driven by $\Lambda_{xy}$, away from the $y$ axis towards
the $x$ axis is 0.95 degrees. Using the spin length $S=2$, the out-of-easy-
axis magnetic moments per spin are $|m_{z}|=0.063\mu_{\mathrm{B}}$ and
$|m_{x}|=0.067\mu_{\mathrm{B}}$ reproducing the experimentally determined
deviations, $m_{z}$ = 0.063(5)$\mu_{B}$ and $m_{x}$ = 0.067(5)$\mu_{B}$ [31].

The saturation value of the magnetization for H$\,\parallel\,$y calculated
from the mean-field model is 4.5% higher than the experimentally observed,
Fig. 2 (a). Reason for the failure to reproduce the saturation magnetization
and the spin-flop field with the same set of magnetic-field independent
parameters could be magnetostriction [56]. Magnetostriction, as was proposed
in Ref. [46], could also be the reason why the mean-field model does not
reproduce the frequency of $F_{4}$ close to the spin-flop field, 32 T in Fig.
5.

##  5 Discussion

###  5.1 Spin-wave excitations from the mean-field model

We found that the the mean-field model quantitatively describes the magnetic
field dependence of the frequencies of spin-waves $F_{4}$, $F_{5}$, $F_{6}$,
and $F_{7}$, Fig. 4. The modes $F_{4}$, $F_{5}$, $F_{6}$, and $F_{7}$ have a
linear field dependence with the slope close to $\pm 1$$\,\mathrm{cm}^{-1}$T-1
when the field is along the easy axis $y$. This slope corresponds to a spin-
wave excitation with $\Delta m_{s}=\pm 1$, assuming $g\approx 2$. Other
studies also found a g-factor close to 2 [46]. Other candidates for the
$\Delta m_{s}=\pm 1$ spin-wave excitations are $F_{3}$ and $F_{8}$. However,
both of these modes have two branches degenerate in zero field. The
magnetization measurements, inset of Fig. 2(a), indicate bi-axial magnetic
anisotropy in LiFePO4 which lifts the degeneracy of magnetic resonances in
zero field. Therefore, $F_{3}$ and $F_{8}$ cannot be consistently included
into the mean-field description.

The spin-waves of the mean-field model have oscillating spin components,
$\delta\mathbf{S}_{i}=\mathbf{S}_{i}-\bar{\mathbf{S}}_{i}$, perpendicular to
the equilibrium direction of the $i$-th spin, $\bar{\mathbf{S}}_{i}$. The
spin-wave couples to the magnetic field of radiation if the oscillating spin
component of the whole magnetic unit cell is finite,
$\mathbf{H}^{\omega}\cdot\left(\sum_{i=1}^{4}\delta\mathbf{S}_{i}\right)$. The
equilibrium direction of the spins is aligned to the easy axis $y$ within few
degrees in LiFePO4. The selection rules, Table 2, show that $F_{4}$ and
$F_{5}$ are excited by the $\mathit{H}^{\omega}_{z}$ component of radiation
and modes $F_{6}$ and $F_{7}$ by the $\mathit{H}^{\omega}_{x}$ component,
which both are perpendicular to $\bar{\mathbf{S}}_{i}$. The magnetic field
dependence of intensities of the strongest modes $F_{4}$ and $F_{6}$ is well
described by the mean-field model. Firstly, $F_{4}$ is
$\mathit{H}^{\omega}_{z}$\- and $F_{6}$ is $\mathit{H}^{\omega}_{x}$ -active
in zero field, Table 2. Secondly, as $H_{y}$ increases, $F_{4}$ becomes
$\mathit{H}^{\omega}_{x}$-active and $F_{6}$ becomes
$\mathit{H}^{\omega}_{z}$-active, Fig. S4 in the Supplementary Material. Thus,
for modes $F_{4}$ and $F_{6}$ the agreement between theory and experiment is
good.

The experimental and theoretical selection rules of magnetic dipole transition
for $F_{5}$ agree, it is $\mathit{H}^{\omega}_{z}$-active. For the spin-wave
$F_{7}$ the theory predicts $\mathit{H}^{\omega}_{y}$-activity, although it is
$\mathit{H}^{\omega}_{x}$-active in the experiment, see Fig. S4 in the
Supplementary Material. Overall, theory underestimates $F_{5}$ and $F_{7}$
magnetic dipole transition intensity by two orders of magnitude. It is not
surprising as the modes $F_{5}$ and $F_{7}$ are relatively weak as compared to
$F_{4}$ and $F_{6}$ and therefore they are sensitive to the composition of the
spin-wave state. If the mean-field model, as an approximation, does not give
the true spin-wave state, the weak intensities could be seriously affected. As
observed experimentally, Fig. 3, $F_{7}$ is electric-dipole active in
addition. The coupling of spins to the electric field was not included in the
mean-field model.

Similar to the quantum-mechanical formulas of the magnon dispersion relation
of former studies [30, 60, 31], the classical expressions for the zero-field
resonance frequencies of the magnon modes can be derived. For the two
strongest spin-waves $F_{4}$ and $F_{6}$

$\nu_{4/6}\approx
2S\sqrt{\Lambda_{x/z}\left(4(J_{yz}+J_{xy})+\Lambda_{z/x}\right)},$ (5)

while the zero-field frequencies of the weaker $F_{5}$ and $F_{7}$ are

$\nu_{5/7}\approx
2S\sqrt{\left(\frac{\Lambda_{x/z}}{1-\sqrt{J_{xy}/\Lambda_{x/z}}}-J_{xz}\right)\left(4J_{yz}+\Lambda_{z/x}\right)},$
(6)

where we neglected the weak Dzyaloshinskii-Moriya interaction and the single-
ion anisotropy $\Lambda_{xy}$ terms. While these two terms are necessary to
give finite magnetic dipole activity to the weak $F_{5}$ and $F_{7}$
resonances by breaking the equivalence of $S_{1}$ and $S_{2}$ ($S_{3}$ and
$S_{4}$, respectively) sublattices, they do not change the resonance
frequencies significantly.

As follows from Eq. (5) and Eq. (6), if $J_{xy}=J_{xz}=0$, $F_{4}$ and $F_{5}$
are degenerate in zero field, $\nu_{4}=\nu_{5}$, and also $\nu_{6}=\nu_{7}$.
In this case the nearest-neighbor $(100)$ planes of the $\\{S_{1},S_{4}\\}$
and $\\{S_{2},S_{3}\\}$ sub-lattices, separated by $a/2$, are decoupled from
each other, thus, their in-phase and out-of-phase excitations with respect to
each other are degenerate. Consequently, $F_{4}$ and $F_{6}$ can be considered
as the in-phase while $F_{5}$ and $F_{7}$ as the out-of-phase resonances of
the nearest-neighbor $(100)$ planes. Without Dzyaloshinskii-Moriya interaction
and $\Lambda_{xz}$ anisotropy the total oscillating magnetic dipole moment of
the unit cell produced by $F_{5}$ and $F_{7}$ is zero. This explains the weak
intensity of $F_{5}$ and $F_{7}$ compared to $F_{4}$ and $F_{6}$ in the THz
absorption spectrum. Furthermore, the correspondence between the INS magnon
dispersion interpreted in the two-spin unit cell scheme [30, 31, 32] and our
$\Gamma$-point optical experiments can also be formulated based on the mean-
field results. Namely, $F_{4}$ and $F_{6}$ correspond to the spin-waves
observed in the zone center, $\mathbf{Q}=(0,2,0)$ [31] or $\mathbf{Q}=(0,0,2)$
[32] while $F_{5}$ and $F_{7}$ are zone-boundary excitations of the two-spin
unit cell, seen at $\mathbf{Q}=(0,0,1)$ [32], $\mathbf{Q}=(1,1,0)$ [30, 31]
and $\mathbf{Q}=(0,1,1)$ [31] in the INS experiments [30, 31, 32].

###  5.2 Spin excitations beyond the mean-field model

Out of 17 lines appearing below $T_{\mathrm{N}}$ in the THz absorption
spectrum only four can be described by the classical four-spin mean-field
model. The rest can be (_i_) spin-stretching excitations captured only by
multi-boson spin-wave theory or alternatively by crystal-field schemes
including exchange fields, (_ii_) two-magnon excitations (two spin waves with
nearly opposite k vectors), or can even be (_iii_) excitations from impurity
spins. Assuming that the spins are aligned along the $y$-axis the magnetic
symmetry reduces to $Pnma^{\prime}$ [30]. Since all spatial symmetries of the
paramagnetic state remain in the AFM phase, at least in combination with time-
reversal operation, we do not expect new optical phonon modes to emerge below
$T_{\mathrm{N}}$.

We assign absorption lines $F_{1}$, $F_{2}$ and $F_{3}$ to impurities because
these very weak modes are located below the lowest magnon mode $F_{4}$. In
addition, the frequencies of $F_{1}$ and $F_{2}$ increase linearly in magnetic
field H$\,\parallel\,$z, not coinciding with easy-axis direction $y$. Previous
works have found that Fe2+ at Li+ site has zero field splitting
7.3$\,\mathrm{cm}^{-1}$ (220 GHz) [48]. The lowest impurity absorption line in
our spectrum is $F_{1}$ at 18$\,\mathrm{cm}^{-1}$ in zero field. This suggests
that we are observing different impurities than reported in Ref. [48].

The mean-field model does not describe spin excitations $F_{8}$–$F_{17}$.
Several of them are electric-dipole active and have a steep magnetic field
dependence of frequency, suggesting $|\Delta m_{s}|>1$ change of a spin
projection quantum number. This is unusual for a spin-wave excitation but can
be explained by a large single-ion anisotropy ($\Lambda$) which is comparable
or stronger than the exchange coupling ($J$) [37], see Table 1. In that case a
suitable approach is a multi-boson spin-wave theory, which describes more than
four spin-wave excitations in a four-sublattice magnet. Out from the ortho-
phosphate compounds, the multi-boson spin-wave theory has been only applied to
LiCoPO4, a $S=3/2$ spin system [13]. Developing a multi-boson spin-wave theory
for LiFePO4 is a tedious calculation, therefore, it is out of the scope of
this paper.

Some of the observed features can be explained qualitatively in the limit of
zero exchange and Dzyaloshinskii-Moriya coupling. Assuming rotational symmetry
about the $y$ axis in Eq. (4), $\Lambda_{z}=\Lambda_{x}$ , the spins are
parallel to the quantization axis $y$, and the energy levels $E_{m_{s}}$ of
spin $S=2$ are $E_{0}$, $E_{\pm 1}$ and $E_{\pm 2}$. When the
H$\,\parallel\,$y field is applied, the energy difference $E_{+2}-E_{-2}$
increases approximately at a rate 4 $\,\mathrm{cm}^{-1}$T-1, as observed for
the spin-wave excitation $F_{15}$. The electric dipole activity comes from the
on-site spin-induced polarization which in the lowest order of spin operators
is $P\propto\hat{S}_{\alpha}\hat{S}_{\beta}$ ($\alpha,\beta=x,y,z$) [13].
Although $P\propto\hat{S}_{x}^{2}$ and $\hat{S}_{z}^{2}$ (quantization axis is
$y$) couple states different by $\Delta m_{s}=\pm 2$ it does not explain the
$|\Delta m_{s}|\geq 3$ transitions, $F_{10}$, $F_{14}$ and $F_{15}$. However,
in LiFePO4 the single ion anisotropies are not equal,
$\Lambda_{z}\neq\Lambda_{x}$ and mix $E_{0}$ into $E_{\pm 2}$ states, see
Table I in Ref. [32]. Therefore, the selection rule for the electric-dipole
transition, $\Delta m_{s}=2$, and mixing of states gives finite electric-
dipole moment to the $\Delta m_{s}=4$ transition. In a similar manner,
$P\propto\hat{S}_{x}\hat{S}_{y}$ and $\hat{S}_{y}\hat{S}_{z}$ could give rise
to $\Delta m_{s}=\pm 1$ transitions and if the mixing of states is taken into
account, then to the electric-dipole allowed $\Delta m_{s}=\pm 3$ transitions.

Two spin-waves, $\omega_{1}(\mathbf{q}_{1})$ and $\omega_{2}(\mathbf{q}_{2})$,
can be excited by THz radiation of frequency $\omega=\omega_{1}+\omega_{2}$ if
$\mathbf{q}_{1}=-\mathbf{q}_{2}$, which is termed as two-magnon excitation.
The exact frequency dependence of this absorption depends on the coupling
mechanism between the radiation and the spin-wave and on the density of spin-
wave states [61, 62, 63, 64, 65, 66]. This leads to broad absorption bands
with peaks at the highest density of spin-wave states [61, 62, 63, 67, 64, 68,
58], mostly with spin-waves from the edge of the Brillouin zone. Since the
product of the two spin operators has the same time-reversal parity as the
electric dipole moment, the simultaneous excitation of two spin-waves by the
electric field is allowed and this mechanism usually dominates over the
magnetic-dipole active absorption [69]. A relatively broad electric-dipole
active absorption line is $F_{13}$. If
$\omega_{1}(\mathbf{q}_{1})=\omega_{2}(\mathbf{q}_{2})$, the spin-wave
frequency should be $\omega_{1}\approx 60$$\,\mathrm{cm}^{-1}$$=7.4$ meV. At
about the same energy two dispersion curves cross in the $[0,1.5,0]$ Brillouin
zone point of the two-spin unit cell [31, 32]. The $[0,0.5,0]$ point,
equivalent to $[0,1.5,0]$, is the Brillouin zone boundary of the four-spin
unit cell and therefore we expect anti-crossing of two dispersion curves which
leads to increase in the density of magnon states at this point. Thus,
considering the linewidths, energy scale, and the electric-dipole activity,
$F_{13}$ could be a two-magnon excitation. Another candidate for a two-magnon
excitation is the electric dipole active $F_{12}$. Although it is relatively
broad in zero field, it has a complicated field dependence in
H$\,\parallel\,$y, see Fig. S5 in the Supplementary Material, what can not be
explained within a simple model of two-magnon excitation.

##  6 Summary

We studied the magnetic ground state and the spin excitations of the
magnetoelectric antiferromagnet LiFePO4 by magnetization measurements in
magnetic fields up to 120 T and by THz absorption spectroscopy up to 33 T.
Magnetization measurements revealed a spin-flop transition at 32 T before
reaching the saturation at 56 T. We found 17 absorption lines below
175$\,\mathrm{cm}^{-1}$ (5.25 THz) appearing in the magnetically ordered
phase. Based on the magnetic field dependence of the resonance frequencies and
the intensities, we assigned four of them to magnon modes ($F_{4}$-$F_{7}$),
eight to multiboson spin-wave excitations ($F_{8}$-$F_{11}$,
$F_{14}$-$F_{17}$), two to two-magnon excitations ($F_{12}$, $F_{13}$) and the
rest to the absorption by impurity spins ($F_{1}$-$F_{3}$). We applied a mean-
field model, which describes well the four magnon modes ($F_{4}$-$F_{7}$). We
attribute the other spin-wave modes to excitations with $|\Delta m_{s}|>1$
arising due to the large, $S=2$, spin of octahedrally coordinated Fe2+ ions.
Such excitations may become electric-dipole active due to symmetry allowed
coupling between spin-quadrupolar fluctuations and electric polarization. Two
modes, $F_{7}$ and $F_{14}$, are magneto-electric resonances with significant
coupling to both, electric and magnetic field component of radiation.
Additional experiments on magneto-electrically poled samples are needed to
clarify if these two resonances show non-reciprocal directional dichroism [13,
16].

##  7 Acknowledgments

The authors acknowledge the valuable discussions with Karlo Penc and thank
Kirill Amelin and Jakub Vit for fruitful discussions and for the help with the
THz spectroscopy measurements. This project was supported by the Estonian
Research Council grant PRG736, institutional research funding IUT23-3 of the
Estonian Ministry of Education and Research, the European Regional Development
Fund project TK134, by the bilateral program of the Estonian and Hungarian
Academies of Sciences under Contract No. NMK2018-47, by the Hungarian National
Research, Development and Innovation Office—NKFIH Grants No. FK 135003. The
high magnetic field magnetization experiments were supported by LNCMI-CNRS and
HFML-RU/NWO-I, members of the European Magnetic Field Laboratory (EMFL). D. S.
acknowledges the FWF Austrian Science Fund I 2816-N27 and TAI 334-N. The
cooperation of Austrian and Hungarian partners was supported by Austrian
Agency for International Cooperation in Education Research Grant No. WTZ HU
08/2020 and by the Hungarian NKFIH Grant No. 2019-2.1.11-TÉT-2019-00029. V. K.
was supported by the RIKEN Incentive Research Project and B.B. acknowledges
the support by the European Research Council (Grant Agreement
No.835279-Catch-22). The data handling, calculations and figures were done in
Python programming language using libraries NumPy [70], Matplotlib [71], Scipy
[72] and Pandas [73].

L.P., V.K., and D.S. contributed equally to this work.

## References

  * Tokura [2007] Y. Tokura, Multiferroics—toward strong coupling between magnetization and polarization in a solid, J Magn Magn Mater 310, 1145 (2007).
  * Saito  _et al._ [2008] M. Saito, K. Taniguchi, and T. Arima, Gigantic optical magnetoelectric effect in CuB2O4, J. Phys. Soc. Jpn. 77, 013705 (2008).
  * Kézsmárki  _et al._ [2011] I. Kézsmárki, N. Kida, H. Murakawa, S. Bordács, Y. Onose, and Y. Tokura, Enhanced directional dichroism of terahertz light in resonance with magnetic excitations of the multiferroic Ba2CoGe2O7 oxide compound, Phys. Rev. Lett. 106, 057403 (2011).
  * Miyahara and Furukawa [2011] S. Miyahara and N. Furukawa, Theory of magnetoelectric resonance in two-dimensional $S=3/2$ antiferromagnet Ba2CoGe2O7 via spin-dependent metal-ligand hybridization mechanism, J. Phys. Soc. Jpn. 80, 073708 (2011).
  * Bordács  _et al._ [2012] S. Bordács, I. Kézsmárki, D. Szaller, L. Demkó, N. Kida, H. Murakawa, Y. Onose, R. Shimano, T. Rõõm, U. Nagel, S. Miyahara, N. Furukawa, and Y. Tokura, Chirality of matter shows up via spin excitations, Nat. Phys. 8, 734 (2012).
  * Takahashi  _et al._ [2013] Y. Takahashi, Y. Yamasaki, and Y. Tokura, Terahertz magnetoelectric resonance enhanced by mutual coupling of electromagnons, Phys. Rev. Lett. 111, 037204 (2013).
  * Kézsmárki  _et al._ [2014] I. Kézsmárki, D. Szaller, S. Bordács, V. Kocsis, Y. Tokunaga, Y. Taguchi, H. Murakawa, Y. Tokura, H. Engelkamp, T. Rõõm, and U. Nagel, One-way transparency of four-coloured spin-wave excitations in multiferroic materials, Nat Commun 5, 3203 (2014).
  * Kézsmárki  _et al._ [2015] I. Kézsmárki, U. Nagel, S. Bordács, R. S. Fishman, J. H. Lee, H. T. Yi, S.-W. Cheong, and T. Rõõm, Optical diode effect at spin-wave excitations of the room-temperature multiferroic BiFeO3, Phys. Rev. Lett. 115, 127203 (2015).
  * Bordács  _et al._ [2015] S. Bordács, V. Kocsis, Y. Tokunaga, U. Nagel, T. Rõõm, Y. Takahashi, Y. Taguchi, and Y. Tokura, Unidirectional terahertz light absorption in the pyroelectric ferrimagnet CaBaCo4O7, Phys. Rev. B 92, 214441 (2015).
  * Kuzmenko  _et al._ [2015] A. M. Kuzmenko, V. Dziom, A. Shuvaev, A. Pimenov, M. Schiebl, A. A. Mukhin, V. Y. Ivanov, I. A. Gudim, L. N. Bezmaternykh, and A. Pimenov, Large directional optical anisotropy in multiferroic ferroborate, Phys. Rev. B 92, 184409 (2015).
  * Toyoda  _et al._ [2015] S. Toyoda, N. Abe, S. Kimura, Y. H. Matsuda, T. Nomura, A. Ikeda, S. Takeyama, and T. Arima, One-way transparency of light in multiferroic CuB2O4, Phys. Rev. Lett. 115, 267207 (2015).
  * Iguchi  _et al._ [2017] Y. Iguchi, Y. Nii, and Y. Onose, Magnetoelectrical control of nonreciprocal microwave response in a multiferroic helimagnet, Nat Commun 8, 10.1038/ncomms15252 (2017).
  * Kocsis  _et al._ [2018] V. Kocsis, K. Penc, T. Rõõm, U. Nagel, J. Vít, J. Romhányi, Y. Tokunaga, Y. Taguchi, Y. Tokura, I. Kézsmárki, and S. Bordács, Identification of antiferromagnetic domains via the optical magnetoelectric effect, Phys. Rev. Lett. 121, 057601 (2018).
  * Yu  _et al._ [2018] S. Yu, B. Gao, J. W. Kim, S.-W. Cheong, M. K. L. Man, J. Madéo, K. M. Dani, and D. Talbayev, High-temperature terahertz optical diode effect without magnetic order in polar FeZnMo3O8, Phys. Rev. Lett. 120, 037601 (2018).
  * Tokura and Nagaosa [2018] Y. Tokura and N. Nagaosa, Nonreciprocal responses from non-centrosymmetric quantum materials, Nat Commun 9, 10.1038/s41467-018-05759-4 (2018).
  * Kocsis  _et al._ [2019] V. Kocsis, S. Bordács, Y. Tokunaga, J. Viirok, L. Peedu, T. Rõõm, U. Nagel, Y. Taguchi, Y. Tokura, and I. Kézsmárki, Magnetoelectric spectroscopy of spin excitations in LiCoPO4, Phys. Rev. B 100, 155124 (2019).
  * Viirok  _et al._ [2019] J. Viirok, U. Nagel, T. Rõõm, D. G. Farkas, P. Balla, D. Szaller, V. Kocsis, Y. Tokunaga, Y. Taguchi, Y. Tokura, B. Bernáth, D. L. Kamenskyi, I. Kézsmárki, S. Bordács, and K. Penc, Directional dichroism in the paramagnetic state of multiferroics: A case study of infrared light absorption in Sr2CoSi2O7 at high temperatures, Phys. Rev. B 99, 014410 (2019).
  * Yokosuk  _et al._ [2020] M. O. Yokosuk, H.-S. Kim, K. D. Hughey, J. Kim, A. V. Stier, K. R. O’Neal, J. Yang, S. A. Crooker, K. Haule, S.-W. Cheong, D. Vanderbilt, and J. L. Musfeldt, Nonreciprocal directional dichroism of a chiral magnet in the visible range, npj Quantum Materials 5, 10.1038/s41535-020-0224-6 (2020).
  * Kimura  _et al._ [2020a] K. Kimura, T. Katsuyoshi, Y. Sawada, S. Kimura, and T. Kimura, Imaging switchable magnetoelectric quadrupole domains via nonreciprocal linear dichroism, Communications Materials 1, 39 (2020a).
  * Kimura  _et al._ [2020b] S. Kimura, M. Matsumoto, and H. Tanaka, Electrical switching of the nonreciprocal directional microwave response in a triplon Bose-Einstein condensate, Phys. Rev. Lett. 124, 217401 (2020b).
  * Ogino  _et al._ [2020] M. Ogino, Y. Kaneko, Y. Tokura, and Y. Takahashi, Gyrotropic birefringence via electromagnon resonance in a multiferroic of spin origin, Phys. Rev. Research 2, 023345 (2020).
  * Kimura  _et al._ [2021] S. Kimura, N. Terada, M. Hagiwara, M. Matsumoto, and H. Tanaka, Electric dipole active magnetic resonance and nonreciprocal directional dichroism in magnetoelectric multiferroic materials in terahertz and millimeter wave regions, Appl. Magn. Reson. 10.1007/s00723-020-01307-w (2021).
  * Toyoda  _et al._ [2021] S. Toyoda, M. Fiebig, T. hisa Arima, Y. Tokura, and N. Ogawa, Nonreciprocal second harmonic generation in a magnetoelectric material, Sci. Adv. 7, 10.1126/sciadv.abe2793 (2021).
  * Vít  _et al._ [2021] J. Vít, J. Viirok, L. Peedu, T. Rõõm, U. Nagel, V. Kocsis, Y. Tokunaga, Y. Taguchi, Y. Tokura, I. Kézsmárki, P. Balla, K. Penc, J. Romhányi, and S. Bordács, In situ electric-field control of THz nonreciprocal directional dichroism in the multiferroic Ba2CoGe2O7, Phys. Rev. Lett. 127, 157201 (2021).
  * Reschke  _et al._ [2022] S. Reschke, D. G. Farkas, A. Strinić, S. Ghara, K. Guratinder, O. Zaharko, L. Prodan, V. Tsurkan, D. Szaller, S. Bordács, J. Deisenhofer, and I. Kézsmárki, Confirming the trilinear form of the optical magnetoelectric effect in the polar honeycomb antiferromagnet Co2Mo3O8, npj Quantum Materials 7, 10.1038/s41535-021-00417-3 (2022).
  * Pimenov  _et al._ [2006] A. Pimenov, A. A. Mukhin, V. Y. Ivanov, V. D. Travkin, A. M. Balbashov, and A. Loidl, Possible evidence for electromagnons in multiferroic manganites, Nat. Phys. 2, 97 (2006).
  * Takahashi  _et al._ [2012] Y. Takahashi, R. Shimano, Y. Kaneko, H. Murakawa, and Y. Tokura, Magnetoelectric resonance with electromagnons in a perovskite helimagnet, Nat. Phys. 8, 121 (2012).
  * Kuzmenko  _et al._ [2018] A. M. Kuzmenko, D. Szaller, T. Kain, V. Dziom, L. Weymann, A. Shuvaev, A. Pimenov, A. A. Mukhin, V. Y. Ivanov, I. A. Gudim, L. N. Bezmaternykh, and A. Pimenov, Switching of magnons by electric and magnetic fields in multiferroic borates, Phys. Rev. Lett. 120, 027203 (2018).
  * Szaller  _et al._ [2014] D. Szaller, S. Bordács, V. Kocsis, T. Rõõm, U. Nagel, and I. Kézsmárki, Effect of spin excitations with simultaneous magnetic- and electric-dipole character on the static magnetoelectric properties of multiferroic materials, Phys. Rev. B 89, 184419 (2014).
  * Li  _et al._ [2006] J. Li, V. O. Garlea, J. L. Zarestky, and D. Vaknin, Spin-waves in antiferromagnetic single-crystal LiFePO4, Phys. Rev. B 73, 024410 (2006).
  * Toft-Petersen  _et al._ [2015] R. Toft-Petersen, M. Reehuis, T. B. S. Jensen, N. H. Andersen, J. Li, M. D. Le, M. Laver, C. Niedermayer, B. Klemke, K. Lefmann, and D. Vaknin, Anomalous magnetic structure and spin dynamics in magnetoelectric LiFePO4, Phys. Rev. B 92, 024404 (2015).
  * Yiu  _et al._ [2017] Y. Yiu, M. D. Le, R. Toft-Petersen, G. Ehlers, R. J. McQueeney, and D. Vaknin, Hybrid excitations due to crystal field, spin-orbit coupling, and spin waves in LiFePO4, Phys. Rev. B 95, 104409 (2017).
  * Szaller  _et al._ [2017] D. Szaller, V. Kocsis, S. Bordács, T. Fehér, T. Rõõm, U. Nagel, H. Engelkamp, K. Ohgushi, and I. Kézsmárki, Magnetic resonances of multiferroic TbFe3(BO3)4, Phys. Rev. B 95, 024427 (2017).
  * Szaller  _et al._ [2020] D. Szaller, K. Szász, S. Bordács, J. Viirok, T. Rõõm, U. Nagel, A. Shuvaev, L. Weymann, A. Pimenov, A. A. Tsirlin, A. Jesche, L. Prodan, V. Tsurkan, and I. Kézsmárki, Magnetic anisotropy and exchange paths for octahedrally and tetrahedrally coordinated Mn2+ ions in the honeycomb multiferroic ${\mathrm{Mn}}_{2}{\mathrm{Mo}}_{3}{\mathrm{O}}_{8}$, Phys. Rev. B 102, 144410 (2020).
  * Rõõm  _et al._ [2020] T. Rõõm, J. Viirok, L. Peedu, U. Nagel, D. G. Farkas, D. Szaller, V. Kocsis, S. Bordács, I. Kézsmárki, D. L. Kamenskyi, H. Engelkamp, M. Ozerov, D. Smirnov, J. Krzystek, K. Thirunavukkuarasu, Y. Ozaki, Y. Tomioka, T. Ito, T. Datta, and R. S. Fishman, Magnetoelastic distortion of multiferroic BiFeO3 in the canted antiferromagnetic state, Phys. Rev. B 102, 214410 (2020).
  * Farkas  _et al._ [2021] D. G. Farkas, D. Szaller, I. Kézsmárki, U. Nagel, T. Rõõm, L. Peedu, J. Viirok, J. S. White, R. Cubitt, T. Ito, R. S. Fishman, and S. Bordács, Selection rules and dynamic magnetoelectric effect of the spin waves in multiferroic BiFeO3, Phys. Rev. B 104, 174429 (2021).
  * Penc  _et al._ [2012] K. Penc, J. Romhányi, T. Rõõm, U. Nagel, A. Antal, T. Fehér, A. Jánossy, H. Engelkamp, H. Murakawa, Y. Tokura, D. Szaller, S. Bordács, and I. Kézsmárki, Spin-stretching modes in anisotropic magnets: Spin-wave excitations in the multiferroic Ba2CoGe2O7, Phys. Rev. Lett. 108, 257203 (2012).
  * Romhányi and Penc [2012] J. Romhányi and K. Penc, Multiboson spin-wave theory for Ba2CoGe2O7: A spin-3/2 easy-plane Néel antiferromagnet with strong single-ion anisotropy, Phys. Rev. B 86, 174428 (2012).
  * Akaki  _et al._ [2017] M. Akaki, D. Yoshizawa, A. Okutani, T. Kida, J. Romhányi, K. Penc, and M. Hagiwara, Direct observation of spin-quadrupolar excitations in Sr2CoSi2O7 by high-field electron spin resonance, Phys. Rev. B 96, 214406 (2017).
  * Legros  _et al._ [2021] A. Legros, S.-S. Zhang, X. Bai, H. Zhang, Z. Dun, W. A. Phelan, C. D. Batista, M. Mourigal, and N. Armitage, Observation of 4- and 6-magnon bound states in the spin-anisotropic frustrated antiferromagnet FeI2, Phys. Rev. Lett. 127, 267201 (2021).
  * Bai  _et al._ [2021] X. Bai, S.-S. Zhang, Z. Dun, H. Zhang, Q. Huang, H. Zhou, M. B. Stone, A. I. Kolesnikov, F. Ye, C. D. Batista, and M. Mourigal, Hybridized quadrupolar excitations in the spin-anisotropic frustrated magnet FeI2, Nat. Phys. 17, 467 (2021).
  * Fishman  _et al._ [2018] R. S. Fishman, J. A. Fernandez-Baca, and T. Rõõm,  _Spin-Wave Theory and its Applications to Neutron Scattering and THz Spectroscopy_ (IOP Concise Physics, Morgan and Claypool Publishers, 1210 Fifth Avenue, Suite 250, San Rafael, CA, 94901, USA, 2018).
  * Chaix  _et al._ [2014] L. Chaix, S. de Brion, S. Petit, R. Ballou, L.-P. Regnault, J. Ollivier, J.-B. Brubach, P. Roy, J. Debray, P. Lejay, A. Cano, E. Ressouche, and V. Simonet, Magneto- to electroactive transmutation of spin waves in ErMnO3, Phys Rev Lett 112, 137201 (2014).
  * Strinić  _et al._ [2020] A. Strinić, S. Reschke, K. V. Vasin, M. Schmidt, A. Loidl, V. Tsurkan, M. V. Eremin, and J. Deisenhofer, Magnetoelectric properties and low-energy excitations of multiferroic FeCr2S4, Phys Rev B 102, 134409 (2020).
  * Santoro and Newnham [1967] R. P. Santoro and R. E. Newnham, Antiferromagnetism in LiFePO4, Acta Crystallogr. 22, 344 (1967).
  * Werner  _et al._ [2021] J. Werner, C. Neef, C. Koo, A. Ponomaryov, S. Zvyagin, and R. Klingeler, Exceptional field dependence of antiferromagnetic magnons in LiFePO4, Phys. Rev. B 103, 174406 (2021).
  * Creer and Troup [1970] J. Creer and G. Troup, The magnetic susceptibility of LiFePO4 and LiCoPO4, Phys. Lett. A 32, 439 (1970).
  * Werner  _et al._ [2020] J. Werner, C. Neef, C. Koo, S. Zvyagin, A. Ponomaryov, and R. Klingeler, Antisite disorder in the battery material LiFePO4, Phys. Rev. Materials 4, 115403 (2020).
  * Baker  _et al._ [2011] P. J. Baker, I. Franke, F. L. Pratt, T. Lancaster, D. Prabhakaran, W. Hayes, and S. J. Blundell, Probing magnetic order in Li$M$PO4 ($M$$=$ Ni, Co, Fe) and lithium diffusion in LixFePO4, Phys. Rev. B 84, 174403 (2011).
  * García-Moreno  _et al._ [2001] O. García-Moreno, M. Alvarez-Vega, F. Garcıá-Alvarado, J. García-Jaca, J. M. Gallardo-Amores, M. L. Sanjuán, and U. Amador, Influence of the structure on the electrochemical performance of lithium transition metal phosphates as cathodic materials in rechargeable lithium batteries: A new high-pressure form of Li$M$PO4 ($M$ $=$ Fe and Ni), Chem. Mater. 13, 1570 (2001).
  * Portugall  _et al._ [1999] O. Portugall, N. Puhlmann, H. U. Müller, M. Barczewski, I. Stolpe, and M. von Ortenberg, Megagauss magnetic field generation in single-turn coils: new frontiers for scientific experiments, J. Phys. D: Appl. Phys. 32, 2354 (1999).
  * Takeyama  _et al._ [2012] S. Takeyama, R. Sakakura, Y. H. Matsuda, A. Miyata, and M. Tokunaga, Precise magnetization measurements by parallel self-compensated induction coils in a vertical single-turn coil up to 103T, J. Phys. Soc. Japan 81, 014702 (2012).
  * Turov [1963] E. A. Turov,  _Physical properties of magnetically ordered crystals_ (Moscow: Izdat. Acad. Sci. SSSR, 1963).
  * Kittel [1948] C. Kittel, On the theory of ferromagnetic resonance absorption, Phys. Rev. 73, 155 (1948).
  * Liang  _et al._ [2008] G. Liang, K. Park, J. Li, R. E. Benson, D. Vaknin, J. T. Markert, and M. C. Croft, Anisotropy in magnetic properties and electronic structure of single-crystal LiFePO4, Phys. Rev. B 77, 064414 (2008).
  * Werner  _et al._ [2019] J. Werner, S. Sauerland, C. Koo, C. Neef, A. Pollithy, Y. Skourski, and R. Klingeler, High magnetic field phase diagram and failure of the magnetic Grüneisen scaling in LiFePO4, Phys. Rev. B 99, 214432 (2019).
  * Gilbert [2004] T. L. Gilbert, A phenomenological theory of damping in ferromagnetic materials, IEEE Trans. Magn. 40, 3443 (2004).
  * Peedu  _et al._ [2019] L. Peedu, V. Kocsis, D. Szaller, J. Viirok, U. Nagel, T. Rõõm, D. G. Farkas, S. Bordács, D. L. Kamenskyi, U. Zeitler, Y. Tokunaga, Y. Taguchi, Y. Tokura, and I. Kézsmárki, Spin excitations of magnetoelectric LiNiPO4 in multiple magnetic phases, Phys. Rev. B 100, 024406 (2019).
  * Stanislavchuk  _et al._ [2017] T. Stanislavchuk, D. S. Middlemiss, J. S. Syzdek, Y. Janssen, R. Basistyy, A. A. Sirenko, P. G. Khalifah, C. Grey, and R. Kostecki, Infrared-active optical phonons in LiFePO4 single crystals, J. Appl. Phys. 122, 045107 (2017).
  * Tian  _et al._ [2008] W. Tian, J. Li, J. W. Lynn, J. L. Zarestky, and D. Vaknin, Spin dynamics in the magnetoelectric effect compound LiCoPO4, Phys. Rev. B 78, 184429 (2008).
  * Halley and Silvera [1965] J. W. Halley and I. Silvera, Odd-exciton magnon interaction and explanation of anomalous far-infrared absorption in antiferromagnetic FeF2, Phys. Rev. Lett. 15, 654 (1965).
  * Allen  _et al._ [1966] S. J. Allen, R. Loudon, and P. L. Richards, Two-magnon absorption in antiferromagnetic MnF2, Phys. Rev. Lett. 16, 463 (1966).
  * Loudon [1968] R. Loudon, Theory of infra-red and optical spectra of antiferromagnets, Adv. Phys. 17, 243 (1968), https://doi.org/10.1080/00018736800101296 .
  * Tanaka and Nagasaka [1990] Y. Tanaka and K. Nagasaka, Far infrared activity due to two-magnon excitation in antiferromagnetic MEM(TCNQ)2, Solid State Commun 73, 735 (1990).
  * Tanabe  _et al._ [2005] Y. Tanabe, Y. Fujimaki, K. Kojima, S. Uchida, S. Onari, T. Matsuo, S. Azuma, and E. Hanamura, Direct optical excitation of two and three magnons in $\alpha$-Fe2O3, Low Temp Phys+ 31, 780 (2005).
  * Filho  _et al._ [2015] C. C. Filho, P. Gomes, A. García-Flores, G. Barberis, and E. Granado, Two-magnon raman scattering in LiMnPO4, J. Magn. Magn. Mater. 377, 430 (2015).
  * Fert  _et al._ [1978] A. Fert, D. Bertrand, J. Leotin, J. Ousset, J. Magariño, and J. Tuchendler, Excitation of two spin deviations by far infrared absorption in FeI2, Solid State Commun. 26, 693 (1978).
  * Hildebrand  _et al._ [1999] M. G. Hildebrand, A. Slepkov, M. Reedyk, G. Amow, J. E. Greedan, and D. A. Crandles, Far-infrared optical properties of antiferromagnetic SmTiO3, Phys Rev B 59, 6938 (1999).
  * Richards [1967] P. L. Richards, Far infrared absorption by two magnon excitations in antiferromagnets, J. Appl. Phys. 38, 1500 (1967), http://dx.doi.org/10.1063/1.1709683 .
  * Harris  _et al._ [2020] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. Fernández del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant, Array programming with NumPy, Nature 585, 357–362 (2020).
  * Hunter [2007] J. D. Hunter, Matplotlib: A 2D graphics environment, Computing in Science Engineering 9, 90 (2007).
  * Virtanen  _et al._ [2020] P. Virtanen, R. Gommers, T. E. Oliphant, M. Haberland, T. Reddy, D. Cournapeau, E. Burovski, P. Peterson, W. Weckesser, J. Bright, S. J. van der Walt, M. Brett, J. Wilson, K. J. Millman, N. Mayorov, A. R. J. Nelson, E. Jones, R. Kern, E. Larson, C. J. Carey, İ. Polat, Y. Feng, E. W. Moore, J. VanderPlas, D. Laxalde, J. Perktold, R. Cimrman, I. Henriksen, E. A. Quintero, C. R. Harris, A. M. Archibald, A. H. Ribeiro, F. Pedregosa, P. van Mulbregt, and SciPy 1.0 Contributors, SciPy 1.0: Fundamental algorithms for scientific computing in Python, Nat. Methods 17, 261 (2020).
  * Wes McKinney [2010] Wes McKinney, Data structures for statistical computing in Python, in  _Proceedings of the 9th Python in Science Conference_, edited by Stéfan van der Walt and Jarrod Millman (2010) pp. 56 – 61.


# Forensic analysis of the Turkey 2023 presidential election reveals extreme
vote swings in remote areas

Peter Klimek1,2,3∗, Ahmet Aykaç4, Stefan Thurner1,2,3,5 1Section for Science
of Complex Systems, CeDAS, Medical University of Vienna, Spitalgasse 23,
A-1090, Austria
2Complexity Science Hub Vienna, Josefstädter Strasse 39, A-1080 Vienna,
Austria
3Supply Chain Intelligence Institute Austria (ASCII), Josefstädter Strasse 39,
A-1080 Vienna, Austria
4Theseus International Management Institute, Sophia Antipolis, France
5Santa Fe Institute, 1399 Hyde Park Road, Santa Fe, NM 85701, USA
<EMAIL_ADDRESS>
###### Abstract

Concerns about the integrity of Turkey’s elections have increased with the
recent transition from a parliamentary democracy to an executive presidency
under Recep Tayyip Erdoğan. Election forensics tools are used to identify
statistical traces of certain types of electoral fraud, providing important
information about the integrity and validity of democratic elections. Such
analyses of the 2017 and 2018 Turkish elections revealed that malpractices
such as ballot stuffing or voter manipulation may indeed have played a
significant role in determining the election results. Here, we apply election
forensic statistical tests for ballot stuffing and voter manipulation to the
results of the 2023 presidential election in Turkey. We find that both rounds
of the 2023 presidential election exhibit similar statistical irregularities
to those observed in the 2018 presidential election, however the magnitude of
these distortions has decreased. We estimate that 2.4% (SD 1.9%) and 1.9% (SD
1.7%) of electoral units may have been affected by ballot stuffing practices
in favour of Erdoğan in the first and second rounds, respectively, compared to
8.5% (SD 3.9%) in 2018. Areas with smaller polling stations and fewer ballot
boxes had significantly inflated votes and turnout, again, in favor of
Erdoğan. Furthermore, electoral districts with two or fewer ballot boxes were
more likely to show large swings in vote shares in favour of Erdoğan from the
first to the second round. Based on a statistical model, we estimate that
these swings translate into 342,000 excess votes (SD 4,900) or 0.64% for
Erdoğan. Our results suggest that Turkish elections continue to be riddled
with statistical irregularities, that may be indicative of electoral fraud.

The first round of the 2023 presidential election was held on 14 May 2023 and
pitted incumbent President Recep Tayyip Erdoğan of the Justice and Development
Party (AKP) against opposition candidate Kemal Kiliçdaroğlu, who led an
alliance of six opposition parties. The election was seen as the first real
threat to Erdoğan’s presidency in a long time, as Turkey was reeling from a
prolonged economic crisis with inflation rates of up to 85% and the aftermath
of devastating earthquakes in February 2023 that killed more than 50,000
people, coupled with public outrage at the government’s slow response to these
crises WaPo23 . Despite polls to the contrary, Erdoğan won by a healthy margin
over Kiliçdaroğlu with 49 per cent to 45 per cent of the vote, sending Turkey
into a run-off between the two candidates on 28 May. Erdoğan won this second
round of the election with a margin of 4.36% (52.18% versus 47.82%) over
Kiliçdaroğlu.

While elections in Turkey are generally considered to be free and fair, the
playing field is clearly not level. For example, according to Turkish media
monitoring, Erdoğan received almost 33 hours of airtime on the main state
television channel, compared to 32 minutes for Kiliçdaroğlu APN . These and
other imbalances became possible after Erdoğan transformed Turkey from a
parliamentary democracy to an executive presidency following a 2017 referendum
in which a reform package narrowly won. That election, however, was marred by
allegations of fraud, as unverified videos and reports emerged on social media
showing various forms of electoral malpractice, such as ballot stuffing
(casting multiple votes for one candidate) and voter coercion (preventing
potentially opposing voters from casting their ballots) Observer17 ; NYT .
Also practices such as handing out cash to supporters have been reported
Guardian .

The emerging field of electoral forensics seeks to diagnose the extent to
which a particular type of malpractice may have affected the outcome of an
election, in order to identify electoral malpractice in a timely and fully
quantitative manner Levin12 . A disproportionate abundance of round numbers
was often the focus of early work in election forensics. The basic principle
of these tests is that humans have a particular tendency to favour round
numbers, or numbers with certain digits, when producing results. These
tendencies are at odds with the statistics of the expected number and digit
distributions of fair elections Cantu11 ; Beber12 , including violations of
Benford’s Law Mebane08 ; Pericchi11 ; Gueron22 . However, it has been shown
that such digit-based tests need to be combined with contextual information
such as country-specific risk factors, socio-economic inequalities or ethnic
fractionalisation Montgomery15 .

As a result, there has been a growing interest in forensic electoral testing
that attempts to provide “mechanistic” or generative models of the impact of
specific electoral malpractices on the expected distributions of votes and
turnout across polling stations, as well as on the correlations between votes
and turnout Myagkov09 ; Levin09 ; Klimek12 ; Jimenez17 ; Rozenas17 ; Jimenez11
; Lacasa19 . The basic rationale of such approaches is to consider elections
as large-scale natural experiments in which a population is divided into a
large number of electoral units in which each registered voter makes the
decision to (i) cast a valid ballot or not, and (ii) vote for a particular
candidate. The large number of electoral units in most countries means that
certain statistical regularities can be expected to hold. Election forensics
then tests whether deviations from these regularities are consistent with
specific types of fraud. Similar principles can be used to apply machine
learning models for election forensics Zhang19 . These statistical tools are
often complemented by analyses of secondary data, such as exit polls or survey
and sampling data Prado11 ; Hausmann11 .

In a forensic analysis of the 2017 election, we confirmed that the election
records indeed show specific statistical irregularities that point toward
ballot stuffing and voter coercion Klimek18 . More specifically, ballot
stuffing, when carried out on a large scale, typically leads to a correlation
between voter turnout and the vote for the candidate Klimek12 . Voter
manipulation or coercion tends to inflate votes and turnout in smaller or more
remote regions, where opponents are easier to identify and irregularities are
less likely to be observed or reported Jimenez17 . For the constitutional
referendum, we found that 11% of areas were potentially affected by ballot
stuffing and that removing the affected influences from the data would have
turned the overall vote of the referendum from “Yes” to “No”. There were also
small but significant traces of ballot stuffing. Similar statistical
irregularities were also observed in the 2018 presidential and parliamentary
elections.

Here we ask whether similar forensic patterns are also present in both rounds
of the 2023 presidential election YSK . We test whether statistical
fingerprints of ballot stuffing or voter coercion can also be found in the
2023 results. We also seek to clarify the contribution of small and/or remote
electoral units to such irregularities. In particular, we develop a
statistical test to compare the results at the polling station level for the
first and second rounds of the presidential election. With this test, we ask
whether boxes in such areas are more likely to show improbably large swings of
votes in favour of either candidate from one round to the next.

##  I Data and Methods

We analyse the official and final election results of both rounds of the 2023
presidential election, as provided by the election commission YSK . At the
finest level of aggregation available, the data for the first round consists
of 191,872 electoral units (here also referred to as “boxes”) for which we
consider the size of the electorate (number of eligible voters), the number of
valid votes and the number of votes for Recep T. Erdoğan. The same variables
were extracted from the data for the second round, which contained 192,214
electoral units. Of these, 12 and 18 units respectively were removed because
their number of valid votes was zero.

When applying the ballot-stuffing and voter-rigging tests, all electoral units
with an electorate of less than one hundred voters are excluded in order to
rule out that the results are driven by small number artefacts. This reduces
the number of electoral units used to 180,301 and 180,629 respectively. The
number of valid votes divided by the size of the electorate, $n_{i}^{(j)}$,
for a ballot box $i$ is called the “turnout”, $t_{i}^{(j)}$, for the first
($j=1$) or second ($j=2$) round of the election. The number of votes,
$V_{i}^{c,j}$, for candidate $c$ divided by the valid votes is called the
“vote share”, $v_{i}^{(c,j)}$ .

The election data includes multiple administrative levels for each ballot box,
namely 81 provinces, 948 districts and 28,268 counties. The number of ballot
boxes is very unevenly distributed across these districts. In particular, we
identify districts that are typically smaller and more remote as those that
have only two or fewer ballot boxes.

We apply ballot stuffing Klimek12  and voter coercion Jimenez17  to the data
as previously reported Klimek18 . In order to investigate whether there is
systematic bias in vote shifts, we propose the following test procedure. In
the first round, Sinan Oğan also ran and received 5.17% of the vote, after
which he supported Erdoğan in the run-off. Therefore, we compare the vote
shares of Erdoğan ($c=E$) or Oğan ($c=O$) in round 1 with Erdoğan’s shares in
round 2. To do this, we compute the “vote shift” at the ballot box level as
$\delta
v_{i}=\frac{V_{i}^{(E,2)}}{t_{i}^{(2)}}-\frac{V_{i}^{(E,1)}+V_{i}^{(O,1)}}{t_{i}^{(1)}}$.
We assume that there is a general population-level trend in how preferences
for a candidate may have shifted between the first and second rounds, which
can be obtained as the mode of the distribution of $\delta v_{i}$, denoted as
$\bar{\delta v}$. Algorithmically, we estimate $\bar{\delta v}$ using MatLab’s
kernel density estimation procedure. Our null hypothesis is that the vote
shifts at the ballot box level are symmetrically distributed around this mode,
i.e. the expectation $E$ for the deviations from $\bar{\delta v}$ is zero,
$H_{0}:E(\delta v_{i}\bar{\delta v})=0$.

Let $B^{+/-}$ be the set of ballot boxes for which this deviation is
greater/smaller than zero, $B^{+/-}=\\{i|\delta v_{i}-\bar{\delta v}>/<0\\}$.
A symmetrized vote shift distribution can then be constructed by replacing,
for example, $(\delta v_{i}-\bar{\delta v})\forall i\in B^{+}$ by values of
$-(\delta v_{j}-\bar{\delta v})$, where $j$ was randomly sampled from $B^{-}$.
For the case were vote shifts typically favor one candidate, say $E(|\delta
v_{i}\bar{\delta v}|)_{B^{+}}>E(|\delta v_{i}\bar{\delta v}|)_{B^{-}}$
(meaning that the expectation value is taken over all $i$ in $B^{+/-}$), one
replaces vales in $B^{+}$ with values sample from within $B^{-}$. This gives a
model estimate for a symmetrized vote shift distribution from which corrected
vote totals for individual candidates can be estimated. By comparing the
actual vote tallies with expectations from the model with symmetrized vot
shift distributions, one obtains the number of excess votes due to large vote
swings.

##  II Results

The cumulative percentage of votes for Erdoğan is shown as a function of
turnout in Figure 1. For each turnout level (x-axis), the share of votes from
boxes with that turnout level or lower is shown on the y-axis. In 2018, the
share of votes exceeds the 50% threshold only if we include voting boxes with
a turnout of more than 90%. In 2023, we observe a similarly shaped curve with
an overall higher turnout (the curve is shifted to the right) but without
crossing the 50% threshold in the first round of the election. For the second
round we find again a qualitatively similar curve as in 2018 that crosses the
50% threshold at a turnout of around 92%.

Figure 1: Votes for Erdoğan as a function of voter turnout for 2018 and 2023.
For a given turnout level, the cumulative vote share of ballot boxes with this
or lower turnout is shown. In 2018, a majority of more than 50% is achieved by
including boxes with a turnout of more than 90% (blue). In the first round in
2023 (red line) we observe similar characteristics with higher turnout and
lower vote shares (below 50%), whereas in the second round (yellow) Erdoğan
reached a majority of votes.

###  II.1 Ballot stuffing test.

In Figure 2, we test for the presence of electoral malpractices that lead to
vote-turnout correlations (such as ballot stuffing) using so-called electoral
fingerprints, i.e., a 2-d histogram of the vote-turnout distribution. The
fingerprint for the first round and second round of the 2023 Turkish
presidential election is shown in Figure 2(A,D), respectively. The colour
intensity (blue) indicates the number of boxes with the corresponding
percentage of votes (x-axis) and turnout (y-axis), together with a box plot of
the distribution of turnout for a given percentage of votes. If there were no
non-linear correlations between votes and turnout, the bulk of the
distribution in Figure 2(A,D) should be circular or elliptical symmetric.
Malpractices such as ballot stuffing would inflate turnout and simultaneously
increase vote shares, breaking the elliptical symmetry in the fingerprints if
the number of affected boxes is large enough.

Considering the region of high voting and turnout, in both rounds we see a
smearing of the bulk towards inflated votes and turnout, towards values of
100% votes for Erdoğan and 100% turnout. To assess whether such deviations
between symmetric and biased fingerprints are statistically significant, we
run a parametric test that was proposed previously Klimek18 . This model is
designed to test if the observed deviations from the normal distribution in
vote and turnout shares can be better explained by a model where ballot
stuffing occurs in a given fraction of ballot boxes (fraud parameter $f$). To
fit and evaluate this statistical model, we follow a previously described
strategy Klimek18  and further restrict the analysis to boxes with a vote and
turnout share of more than 25%.

For the first round we find a fraud parameter of $f=0.024$ with a standard
deviation (SD) of $0.019$ and for the second round $f=0.19$ (SD $1.7\%$). Note
that for 2018 we find higher values with $f=0.085$ (SD $0.039$). Therefore,
the test suggests that the number of ballot boxes affected by such statistical
irregularities has decreased from 2018 to 2023, to a point where the ballot
stuffing test does not detect statistically significant effects.

Figure 2: Forensic electoral fingerprints for the two rounds of the 2023
presidential election. The fingerprints for (A) first round and (B) second
round show joint vote-turnout distributions with color intensity encoding the
number of ballot boxes with a given vote (y-axis) and turnout (x-axis). For
both elections, we find a visible correlation in the region of high vote and
turnout (e.g. more than 80%), which can be associated with ballot stuffing. A
box plot (red horizontal boxes) shows the 25th, 50th and 75th percentiles of
voter turnout as a function of votes (whiskers indicate the 95% confidence
interval). (C) To adjust for regional characteristics, the fingerprints can be
adjusted by rescaling the vote and turnout shares by their typical levels in
the unit’s region, resulting in the standardised fingerprint shown. (D) Traces
of voter coercion can be identified by comparing the standardised fingerprints
of small (red lines) and large (blue) units, as voter coercion results in
their being shifted towards inflated votes and turnout, as observed here.
(E,F) The standardised fingerprints of the second round are similar to those
of the first round.

###  II.2 Voter rigging test.

The fingerprints shown in Figure 2(A) and (D) may also show deviations from
elliptical symmetry due to geographical effects. To account for such effects,
it has been suggested to compare the unit with other units in close
geographical proximity Jimenez17 . Here, we compare the vote and turnout
figures of a polling station with the averages observed in other polling
stations in the same constituency. We refer to these rescaled vote and turnout
shares as standardised votes and turnout, respectively. We call their joint
distribution (2D histogram) the “standardised fingerprint” . Standardised
fingerprints are shown in Figure 2(C) for round 1. For the voter coercion
test, we ask whether small and large units have different standardised
fingerprints. The underlying hypothesis of this test is that coercion is more
likely to occur in smaller units because they are more susceptible to coercion
tactics. Reasons for this include that (i) it is easier to identify opponents
in smaller units, (ii) fewer eyewitnesses can be expected, and (iii) election
observers are less likely to be present. In line with these assumptions, voter
manipulation suggests that the standardised fingerprints of small units are
biased towards increased voting and turnout compared to larger units.

We use different definitions of “small units”. Figure 2(B) shows the
standardised fingerprints for small (red) and large (blue) units, where small
units are those in the lowest $p=10$th percentile of all units. It can clearly
be seen that the fingerprints for small units are shifted towards the upper
right corner, see Figure 2(C), which is consistent with voter manipulation.
For the second round, we found almost identical standardised fingerprints, see
Figs. 2(E) and 2(F).

Figure 3: Results of the statistical test for voter manipulation. (A) The
displacement $\delta(p)$ between small and large units for the first (solid
dark magenta line) and second (solid light magenta line) round of the 2023
elections is very narrowly outside the accepted range for a restricted set of
size thresholds, similar as it was in 2018 (solid black line). These
displacements are much smaller displacements in the Russian or Venezuelan
elections (dashed lines); the reference elections are shown as dotted lines.
(B) Units in the 2023 and 2018 Turkish elections are ranked according to their
electorate size. We show the cumulative vote share, $cum_{i}(v)$, calculated
over all units with a size greater than the given rank. As in 2018, we observe
a characteristic “hockey stick” in 2023, meaning that units with high rank
(low electorate size) show a clear tendency to favour Erdoğan.

The magnitude of the displacement between the average standardised votes and
the turnout of small and large units depends on the size threshold p and is
denoted by $\delta(p)$, see Klimek18 ; Jimenez17  for methodological details.
To assess whether this shift is statistically significant, we apply the
Jimenez et al. voter rigging test Jimenez17 . The idea behind this test is to
estimate the expected shifts between small and large units based on a
reference set of trustworthy elections, yielding a range of “acceptable shift
sizes”. We obtain this acceptable region from 21 different reference elections
in ten countries, see Klimek18 ; Jimenez17 . For a given election, one can now
check whether the actual observed displacement between small and large units
for a size threshold $p$ lies within this region (“accepted region”) or not.

The displacement, $\delta(p)$ , is shown in Figure 3(A) for the reference set
of elections (solid lines), the elections in Russia and Venezuela (dashed
line) next to the two rounds in 2023 (solid magenta) and 2018 (solid black)
Turkish elections. For small size thresholds,$p$, all Turkish datasets show
shifts that are slightly outside the acceptable range. This indicates
statistically significant signs of voter manipulation, however, only for a
limited region of thresholds. The shifts in the first round in 2023 are
slightly smaller than in 2018, whereas in the second round they were slighly
larger.

To assess the potential impact of these voter rigging effects in the data, we
rank the units by their electorate size in descending order and calculate the
vote share over all units with smaller ranks (higher electorate size), see
Figure 3(B). In this plot, voter manipulation takes the form of a “hockey
stick”, i.e., a sharp increase for the smallest units. This signal is also
found in Russia and Venezuela, but is absolutely absent in the reference
elections, see the insets in Figure 3(B).

To further show the bias in small units in the first and second round in 2023,
respectively, we compare the fingerprint observed in all units with an
electorate size larger than 100, Figure 4(A,B), with the fingerprints of two
different definitions of “small”. First, we consider only boxes with an
electorate of one hundred or less, resulting in 11,689 units, see Figure
4(C,D). Alternatively, one can consider only boxes from areas with one or two
ballot boxes in total, resulting in 38,662 boxes in the first round, see
Figure 4(E,F). The small units show completely different fingerprints when
compared to the large ones. They show a bimodal distribution with a larger
mode in the high vote-high turnout region and a smaller mode in the low vote-
high turnout region.

Figure 4: The fingerprints in 2023 for small units (C-F) show very different
patterns than those with an electorate size of more than 100 (A,B).
Considering only areas with an electorate size of less than 100 (C,D), or
alternatively ballot boxes from areas with one or two ballot boxes (E,F),
results in a bimodal distribution with a large mode in the region with very
high turnout and votes for Erdoğan.

Comparing Figures 4(E) and (F), there is a tendency for the high vote/high
turnout mode to have shifted further to the upper right of the plot in the
second round compared to the first. To systematically compare such effects, we
examine the shift of vote shares between Erdoğan or Oğan in the first round
and Erdoğan in the second round (see Methods). A scatterplot of these vote
shares is shown in Figure 5(A), which shows that in the vast majority of cases
the vote shares for Erdoğan (or candidates who later endorsed him) were
similar in the first and second rounds; the points cluster around the $x=y$
line. The distribution of vote shifts, $\delta v_{i}$, in Figure 5(B) shows a
clear tendency for large shifts in vote shares from Kiliçdaroğlu in round 1
toward Erdoğan in round 2 to be more common than vice versa (the distribution
is skewed to the right, favouring Erdoğan). Figure 5(C) shows the percentage
of ballot boxes that are in regions with two or fewer ballot boxes in total
for each of the bins in the histogram (bins in the lowest or highest
percentile have been combined). Interestingly, larger districts (in the sense
of having more than two ballot boxes) have very small vote shifts. Strong
positive vote shifts in favour of Erdoğan occur mainly in districts with few
ballot boxes.

Figure 5: Statistical analysis of vote shifts from Kiliçdaroğlu to Erdoğan
from the first to the second round. (A) The vote shares for Erdoğan and Oğan
in the first round are strongly correlated with Erdoğan’s vote share in the
second round, although there are some outliers. (B) A histogram of these vote
shifts shows that these outliers typically favour Erdoğan in the sense that
large shifts from anti to pro Erdoğan are much more common than the other way
around. (C) Looking at the proportion of ballot boxes coming from districts
with two or fewer ballot boxes as a function of vote shifts, we see that these
extreme vote shifts favouring Erdoğan occur predominantly in districts with
fewer ballot boxes.

Similar observations apply, to a lesser extent, to vote shifts in favour of
Kiliçdaroğlu. To compare the magnitude of these vote shifts, we consider a
symmetrized distribution of vote shifts in which vote shifts in favour of
Erdoğan follow the same distribution as those in favour of Kiliçdaroğlu (see
Methods). Overall, we find that Erdoğan received 342,000 (SD 4,900) excess
votes, which corresponds to 0.64% of all valid votes.

##  III Discussion

An electoral forensic analysis of the first and second round of the 2023
presidential election in Turkey identifies statistical irregularities similar
to those observed in the 2018 election and the 2017 constitutional referendum.
However, the estimated magnitude of these irregularities has decreased in 2023
compared to the 2018 presidential election. For 2023, we observe trends in
turnout inflation, as would be expected in the presence of electoral
malpractices such as ballot stuffing. However, the percentage of electoral
units potentially affected by these distortions has fallen to 2.4% (SD 1.9%)
in the first round and to 1.9% (SD 1.7%), making the results statistically
insignificant. In both rounds in 2023, we also observe a tendency for areas
with small electorates to show different voting and turnout patterns compared
to other regions. Such biases are consistent with the presence of voter
coercion or intimidation techniques, to which smaller and more remote
electoral units are more susceptible. The effect size of these deviations is
statistically significant only by a small margin and for a limited range in
the parameter space (meaning these effects vanish when increasing the size
threshold of what constitutes “small” unit above the 20th percentile).

We found, however, suprising irregularities when considering ballot boxes in
areas with two or fewer ballot boxes in total, which typically hints at remote
areas. There are roughly 40,000 such boxes (out of 190,000 in total) that show
completely different trends in their vote–turnout distribution when compared
to the other areas. In such remote areas the vote–turnout distribution becomes
bimodal with high vote/high turnout bulks for both candidates (Erodgan and
Kiliçdaroğlu, respectively), though the Erdoğan mode contains much more ballot
boxes compared to the Kiliçdaroğlu mode. What is more, comparing the first and
the second round of the election, a much higher fraction of such boxes has
flipped from the Kiliçdaroğlu to the Erdoğan mode compared to the other way
around. This high number of districts (around 90%) that flipped from one
candidate to the other at very high turnout levels in the span of no more than
two weeks is certainly surprising. Correcting the election results for such
suprising vote shifts would reduce Erdoğan’s vote tally by about 342,000 votes
or 0.64% of valid votes.

When interpreting the results of forensic voting tests, several limitations
must be borne in mind. First and foremost, none of these tests can provide
incontrovertible evidence of electoral fraud; they provide correlation, not
causation. Although most of the tests make adjustments for regional
characteristics and rule out the possibility that these correlations are
spurious artefacts from small districts, it is not possible to control for all
potential confounders. A positive result in such a test indicates that the
data are compatible with certain types of fraud (ballot stuffing, etc.) and
typically gives a quantitative estimate of how many regions might have been
affected by such malpractice. Other non-fraudulent influences can never be
completely ruled out. In order to understand whether the observed
irregularities may be due to such non-fraudulent phenomena, such as
heterogeneous voter mobilisation through strategic voting, forensic tests must
always be evaluated in conjunction with external information Mebane16 ;
Mebane18 . Conversely, even if the tests give a negative result, one cannot
rule out the presence of other types of fraud, which would require a different
test.

Taken together, these results suggest that the presence of certain types of
electoral malpractice in both rounds of the 2023 presidential election cannot
be ruled out. However, these malpractices appear to have been less frequent in
2023 than in 2018, and less decisive in swinging the vote one way or the
other. There is a notion that Turkish elections have become “free and unfair”
because Turkey’s political playing field is known to be massively tilted in
Erdoğan’s favour WaPo23a . The authors are satisfied that Erdoğan or his
supporters have autocratised the media, the judiciary, civil society and
academia Coskun20 ; Esen20 . Our analysis suggests that it may be more
appropriate to describe Turkish elections as “mostly free and unfair”, as
consistent trends of small but discernible electoral irregularities can be
consistently found by electoral forensic tools. While these statistical
irregularities were not large enough to determine the outcome in 2023 on their
own, they certainly tilt Turkey’s political playing field even further towards
an illiberal democracy.

## References

  * (1) “After bitter campaign, a divided Turkey votes in critical runoff election”, Washington Post, May 28, 2023, https://www.washingtonpost.com/world/2023/05/28/turkey-election-runoff-Erdogan-Kilicdaroglu/.
  * (2) “Turkey’s opposition denounces fairness of vote under Erdoğan”, AP News, May 8, 2023, https://apnews.com/article/turkey-elections-opposition-43f8d2c0aba17958c8e7123cf39eeabf.
  * (3) “Turkey: Videos show electoral fraud and ballot stuffing”, The Observers, France 24, April 18 2017, http://observers.france24.com/en/20170418-turkey-videos-show-electoral-fraud-ballot-stuffing.
  * (4) “Videos Fuel Charges of Fraud in Erdoğan’s Win in Turkey Referendum”, The New York Times, April 18, 2017, https://www.nytimes.com/2017/04/18/world/europe/turkey-referendum-is-haunted-by-allegations-of-voter-fraud.html.
  * (5) “Erdoğan hands banknotes to supporters outside Turkish polling station”, The Guardian News, May 28, 2023, https://www.youtube.com/watch?v=LVVBYl-QcsY.
  * (6) Levin I, Cohn GA, Ordeshook PC, and Alvarez RM (2009) Detecting voter fraud in an electronic voting context: An analysis of the unlimited reelection vote in Venezuela, in EVT/WOTE’09 Proceedings of the 2009 conference on Electronic voting technology/workshop on trustworthy elections, USENIX Association, Berkeley.
  * (7) Cantú F and Saiegh SM (2011) Fraudulent Democracy? An Analysis of Argentina’s Infamous Decade Using Supervised Machine Learning. Political Analysis 19: 409–433.
  * (8) Beber B and Sacco A (2012) What the Numbers Say: A Digit-Based Test for Election Fraud. Political Analysis 20:211–234.
  * (9) Mebane W (2008) Election forensics: The second-digit Benford’s law test and recent American presidential elections, in Election Fraud: Detecting and Deterring Electoral Manipulation, eds Alvarez RM, Hall TE and Hyde SD (Brooking Press, Washington DC), pp 162–181.
  * (10) Pericchi L and Torres D (2011) Quick anomaly detection by the Newcomb-Benford Law, with applications to electoral processes data from the USA, Puerto Rico, and Venezuela. Statist Sci 26:513–527.
  * (11) Gueron E and Pellegrini J (2022) Application of Benford–Newcomb law with base change to electoral fraud detection. Physica A 607:128208.
  * (12) Montgomery JM, Olivella S, Potter JD, Crisp BF (2015) An informed forensics approach to detecting vote irregularities. Political Analysis 23(4):488–505.
  * (13) Lacasa L and Fernandez-Gracia J (2019) Election forensics: Quantitative methods for electoral fraud detection. Forensic Science International 294:e19-e22.
  * (14) Myakgov M, Ordeshook PC, and Shaikin D (2009) The Forensics of Election Fraud, Cambridge University Press.
  * (15) Levin I, Cohn GA, Ordeshook PC, and Alvarez RM (2009) Detecting voter fraud in an electronic voting context: An analysis of the unlimited reelection vote in Venezuela, in EVT/WOTE’09 Proceedings of the 2009 conference on Electronic voting technology/workshop on trustworthy elections, USENIX Association, Berkeley.
  * (16) Klimek P, Yegorov Y, Hanel R, and Thurner S (2012) Statistical detection of systematic election irregularities. Proc Natl Acad Sci USA 109:16469–16473.
  * (17) Jiménez R, Hidalgo M, Klimek P (2017) Testing for voter rigging in small polling stations. Science Advances 3:e1602363.
  * (18) Rozenas A (2017) Detecting election fraud from irregularities in vote-share distributions. Political Analysis 25(1):41–56.
  * (19) Jiménez R (2011) Forensic analysis of the Venezuelan recall referendum. Statist Sci 26:564–583.
  * (20) Zhan M, Alvarez M, Levin I (2019) Using machine learning and synthetic data for possible election anomaly detection. Plos ONE 14(10):e0223950.
  * (21) Prado R and Sansó B (2011) The 2004 Venezuelan presidential recall referendum: Discrepancies between two exit polls and official results. Statist Sci 26:502–512.
  * (22) Hausmann R and Rigobón R (2011) In search of the black swan: Analysis of the statistical evidence of fraud in Venezuela. Statist Sci 26:543–563.
  * (23) Klimek P, Jiménez R, Hidalgo M, Hinteregger A, Thurner S (2018) Forensic analysis of Turkish elections in 2017-2018. Plos ONE 13(10):e0204975.
  * (24) Supreme Board of Elections, https://www.ysk.gov.tr, accessed June 6, 2023.
  * (25) Mebane W (2016) Election Forensics: Frauds Tests and Observation-level Frauds Probabilities. In 2016 Annual Meeting of the Modwest Political Science Association, Chicago, April 7–10.
  * (26) Ferrari D, McAlister K, Mebane W (2018) Developments in Positive Empirical Models of Election Frauds: Varying Dimensions. In 2018 Annual Meeting of the Modwest Political Science Association, Chicago, April 5–18.
  * (27) “Turkey points to a global trend: Free and unfair elections”, Washington Post, May 19, 2023, https://www.washingtonpost.com/opinions/2023/05/19/Erdogan-turkey-autocrats-manipulation-elections/.
  * (28) Coşkun Gü and Kölemen A (2020) Illiberal Democracy or Electoral Autocracy: The Case of Turkey. In The Emergence of Illiberalism (Routledge, 2020), pp. 166-187.
  * (29) Esen B and Yardimci-Geyikci S (2020) The Turkish presidential elections of 24 June 2018. Mediterranean Politics 25(5):682–689.