# Effect of hydrostatic pressure on dynamic dielectric characteristics of
CsH2PO4 ferroelectric

A. S. Vdovych  $|$ ${}^{\textsf{\footnotesize{\mbox{?}}}}$, R. R.
Levitskii${}^{\textsf{\footnotesize{\mbox{?}}}}$, I. R.
Zachek${}^{\textsf{\footnotesize{\mbox{?}}}}$ Corresponding author:
<EMAIL_ADDRESS>
(Received March 15, 2022, in final form April 15, 2022)

###### Abstract

Based on the pseudospin model of the deformed CsH2PO4 crystal within the
Glauber method, the equation for the time-dependent mean value of the
pseudospin is obtained, which is solved in the case of small deviations from
the equilibrium state. Using the solution of the equation, we find expressions
for the longitudinal dynamic dielectric constant and relaxation time. Based on
the proposed parameters of the theory, the temperature and frequency
dependences of the dynamic dielectric constant and the temperature dependence
of the relaxation time are calculated and investigated. A detailed numerical
analysis of the obtained results was performed. The influence of hydrostatic
pressure on the dynamic characteristics of CsH2PO4 is studied.

Key words: ferroelectricity, ferroelectric phase transition, permittivity,
mechanical deformation

###### Abstract

Íà îñíîâ ïñåâäîñïíîâî¿ ìîäåë äåôîðìîâàíîãî êðèñòàëà CsH2PO4 â ðàìêàõ ìåòîäó
Ãëàóáåðà îòðèìàíî ðâíÿííÿ äëÿ çàëåæíîãî âä àñó ñåðåäíüîãî çíàåííÿ ïñåâäîñïíà,
ÿêå ðîçâ’ÿçàíî ó âèïàäêó ìàëèõ âäõèëåíü âä ñòàíó ðâíîâàãè. Âèêîðèñòîâóþè
ðîçâ’ÿçîê ðâíÿííÿ çíàéäåíî âèðàçè äëÿ ïîçäîâæíüî¿ äèíàìíî¿ äåëåêòðèíî¿
ïðîíèêíîñò òà àñó ðåëàêñàö¿. Íà îñíîâ çàïðîïîíîâàíèõ ïàðàìåòðâ òåîð¿
ðîçðàõîâàíî äîñëäæåíî òåìïåðàòóðíó àñòîòíó çàëåæíñòü äèíàìíî¿ äåëåêòðèíî¿
ïðîíèêíîñò òà òåìïåðàòóðíó çàëåæíñòü àñó ðåëàêñàö¿. Ïðîâåäåíî äåòàëüíèé
èñëîâèé àíàëç îòðèìàíèõ ðåçóëüòàòâ. Âèâåíî âïëèâ ãäðîñòàòèíîãî òèñêó íà äèíàìí
õàðàêòåðèñòèêè CsH2PO4.

Ключов слова: ñåãíåòîåëåêòðèêè, ñåãíåòîåëåêòðèíèé ôàçîâèé ïåðåõä, ïðîíèêíñòü,
ìåõàíí äåôîðìàö¿.

##  1 Introduction

A ferroelectric with hydrogen bonds CsH2PO4 (CDP) is an example of a crystal
where the effect of pressure is significant. In this crystal there are two
structurally equivalent types of hydrogen bonds of different lengths (figure
1b). Longer bonds have one equilibrium position for protons, and shorter bonds
have two equilibrium positions. They connect the groups PO4 in chains along
the $b$ -axis (figure 1a); therefore, the crystal is quasi-one-dimensional.

a

b

Figure 1: (Colour online) Primitive cell of CDP crystal in the ferroelectric
phase [25].

At room temperature in the absence of pressure, the crystal is in the
paraelectric phase and has monoclinic symmetry (space group P21/m) [1, 2]. In
this case, the protons on the short bonds are in two equilibrium positions
with equal probability. Below $T_{c}=153$ K, the crystal passes into the
ferroelectric phase (space group P21) [3, 4] with spontaneous polarization
along the crystallographic $b$ axis, and protons with a higher probability
occupy the upper position (figure 1a). Based on dielectric studies [5, 6] it
was established that at pressures $p=p_{c}=0.33$ GPa and
$T_{c}^{\text{cr}}=124.6$ K, double hysteresis loops appear, i.e., there is a
transition to the antiferroelectric phase. Using the neutron diffraction study
[7], it was found that in the antiferroelectric phase, the unit cell of the
CDP crystal doubles along the a-axis, since there are two sublattices in the
form of planes bc, which are polarized antiparallel along $b$-axis and
alternate along a-axis. The symmetry remains monoclinic (space group P21).
Protons on hydrogen bonds are arranged in adjacent sublattices in an
antiparallel manner.

The effect of hydrostatic pressure on the phase transition temperature and
dielectric properties of ferroelectrics Cs(H1-xD${}_{x})_{2}$PO4 was studied
in [8, 5, 6, 9, 10, 11, 12]. The molar heat capacity of CDP was measured in
[13], and was also calculated on the basis of lattice dynamics modelling in
[14, 15]. The important role of proton tunneling on bonds was shown based on
the first-principle calculations [16] and using theoretical calculations on
the basis of quasi-one-dimensional pseudospin model [17].

Using ab initio calculations, piezoelectric coefficients, elastic constants,
and molar heat capacity of CDP [18, 19] were obtained.

The theoretical description of the dielectric properties of CDP at different
values of hydrostatic pressure was carried out in [20, 21] on the basis of the
pseudospin model. However, in these works, the interaction parameters do not
depend on the deformations of the lattice. As a result, it is impossible to
obtain piezoelectric and elastic characteristics of the crystal, and the
critical pressure does not depend on temperature.

In [22], the temperature dependences of lattice strains $u_{1}$, $u_{2}$,
$u_{3}$, $u_{5}$ were measured. It also proposes a quasi-one-dimensional Ising
model for the CDP crystal, in which the interaction parameters are linear
functions of these strains. Based on this model, the temperature behavior of
$u_{j}(T)$ was explained. However, this model does not consider the crystal as
two sublattices and does not make it possible to describe the ferroelectric-
antiferroelectric transition at high pressures.

In [23, 24, 25, 26], there is proposed a two-sublattice pseudospin model of
the deformed crystal CDP, in which the interactions between the nearest
pseudospins in the chain are considered in the mean field approximation. In
this case, the interaction parameters are linear functions of the strains
$u_{j}$. As a result, the temperature dependences of spontaneous polarization,
dielectric constant, piezoelectric coefficients and elastic constants are
calculated, and the influence of hydrostatic and uniaxial pressures and
longitudinal electric field on these characteristics is studied.

In this paper, based on the model [25], the dynamic dielectric characteristics
of CDP in the presence of hydrostatic pressure are calculated.

##  2 The model of CDP crystal

To calculate the dynamic characteristics of CDP, we use the model [25], which
considers a system of protons on O-H … O bonds with two-minimum potential as a
system of pseudospins. The primitive cell has one chain, marked in figure 1 as
“A”. To describe the transition to antiferroelectric phase at high pressures,
paper [25] considers an extended primitive cell formed by two chains (“A” and
“B”). All chains “A” form a sublattice “A”, and all chains “B” form a
sublattice “B”. Each chain in the primitive cell contains two adjacent
tetrahedra PO4 (type “I” and “II”) together with two short hydrogen bonds
(respectively, “1” and “2”). Dipole moments $\vec{d}_{q1}^{\text{A}}$,
$\vec{d}_{q2}^{\text{A}}$, $\vec{d}_{q1}^{\text{B}}$,
$\vec{d}_{q2}^{\text{B}}$, are ascribed to the protons on the bonds.
Pseudospin variables ${\sigma_{q1}^{\text{A}}}/{2}$,
${\sigma_{q2}^{\text{A}}}/{2}$, ${\sigma_{q1}^{\text{B}}}/{2}$,
${\sigma_{q2}^{\text{B}}}/{2}$ describe reorientation of the respective dipole
moments of the base units:
$\vec{d}_{q1,2}^{\text{A,B}}=\vec{\mu}_{q1,2}^{\text{A,B}}{\sigma_{q1,2}^{\text{A,B}}}/{2}$.

Under the stresses maintaining the symmetry of crystal
$\sigma_{1}=\sigma_{xx}$, $\sigma_{2}=\sigma_{yy}$, $\sigma_{3}=\sigma_{zz}$,
$\sigma_{5}=\sigma_{xz}$ ($X$ $\perp$ (b,c), $Y$ $\parallel$ b, $Z$
$\parallel$ c), and in the presence of electric field $E_{2}=E_{y}$
Hamiltonian of the model of CDP is given by [25]:

$\displaystyle\hat{H}=NU_{\text{seed}}+\hat{H}_{\text{short}}+\hat{H}_{\text{long}}+\hat{H}_{E}+\hat{H}^{\prime}_{E},$
(2.1)

where $N$ is the total number of restricted primitive cells.

The first term in (2.1) is “seed” energy, which relates to the heavy ion
sublattice and does not explicitly depend on the configuration of the proton
subsystem. It includes elastic, piezolectric and dielectric parts, expressed
in terms of the electric field $E_{2}$ and strains maintaining the symmetry of
crystal $u_{1}=u_{xx}$, $u_{2}=u_{yy}$, $u_{3}=u_{zz}$, $u_{5}=2u_{xz}$:

$\displaystyle
U_{\text{seed}}=v\left\\{\frac{1}{2}\sum\limits_{j,j^{\prime}}c_{jj^{\prime}}^{E0}u_{j}u_{j}^{\prime}-\sum\limits_{j}e_{2j}^{0}E_{2}u_{j}-\frac{1}{2}\chi_{22}^{\varepsilon
0}E_{2}^{2}\right\\},~{}~{}~{}~{}j,j^{\prime}=1,2,3,5,$ (2.2)

Parameters $c_{jj^{\prime}}^{E0}$, $e_{2j}^{0}$, $\chi_{22}^{u0}$ are the so-
called “seed” elastic constants, “seed” coefficients of piezoelectric stresses
and “seed” dielectric susceptibility, respectively; $v$ is the volume of a
restricted primitive cell. In the paraelectric phase, all coefficients
$e_{ij}^{0}\equiv 0$.

The other terms in (2.1) describe the pseudospin part of Hamiltonian. In
particular, the second term in (2.1) is Hamiltonian of short-range
interactions:

$\hat{H}_{\text{short}}=-2w\sum\limits_{qq^{\prime}}\left(\frac{\sigma_{q1}^{\text{A}}}{2}\frac{\sigma_{q^{\prime}2}^{\text{A}}}{2}+\frac{\sigma_{q1}^{\text{B}}}{2}\frac{\sigma_{q^{\prime}2}^{\text{B}}}{2}\right)\bigl{(}\delta_{{\vec{R}}_{q}{\vec{R}}_{q^{\prime}}}+\delta_{{\vec{R}}_{q}+{\vec{R}}_{b},{\vec{R}}_{q^{\prime}}}\bigr{)}.$
(2.3)

In (2.3), $\sigma_{q1,2}^{\text{A,B}}$ are $z$-components of pseudospin
operator that describe the state of the bond “1” or “2” of the chain “A” or
“B”, in the $q$-th cell, ${\vec{R}}_{b}$ is the lattice vector along
$OY$-axis. The first Kronecker delta corresponds to the interaction between
protons in the chains near the tetrahedra PO4 of type “I”, where the second
Kroneker delta is near the tetrahedra PO4 of type “II”. Contributions into the
energy of interactions between pseudospins near tetrahedra of different type
are identical. Parameter $w$, which describes the short-range interactions
within the chains, is expanded linearly into a series with respect to strains
$u_{j}$:

$w=w^{0}+\sum\limits_{j}\delta_{j}u_{j},(j=1,2,3,5).$ (2.4)

The term $\hat{H}_{\text{long}}$ in (2.1) describes long-range dipole-dipole
interactions and indirect (through the lattice vibrations) interactions
between protons which are taken into account in the mean field approximation:

$\hat{H}_{\text{long}}=NH^{0}+\hat{H}_{2},$ (2.5)

where such notations are used:

$\displaystyle\hat{H}^{0}=\nu_{1}(\eta_{1}^{2}+\eta_{2}^{2})+2\nu_{2}\eta_{1}\eta_{2},$
(2.6)
$\displaystyle\hat{H}_{2}=\sum\limits_{q}\left\\{-\left(2\nu_{1}\eta_{1}+2\nu_{2}\eta_{2}\right)\\!\left(\\!\frac{\sigma_{q1}^{\text{A}}}{2}+\frac{\sigma_{q2}^{\text{A}}}{2}\\!\right)\\!-(2\nu_{2}\eta_{1}+2\nu_{1}\eta_{2})\\!\left(\\!\frac{\sigma_{q1}^{\text{B}}}{2}+\frac{\sigma_{q2}^{\text{B}}}{2}\\!\right)\right\\},$
$\displaystyle\nu_{1}=\nu_{1}^{0}+\sum\limits_{j}\psi_{j1}u_{j},~{}\nu_{2}=\nu_{2}^{0}+\sum\limits_{j}\psi_{j2}u_{j},~{}~{}~{}~{}\langle\sigma_{q1}^{\text{A}}\rangle=\langle\sigma_{q2}^{\text{A}}\rangle=\eta_{1},\quad\langle\sigma_{q1}^{\text{B}}\rangle=\langle\sigma_{q2}^{\text{B}}\rangle=\eta_{2}.$
(2.7)

The parameter $\nu_{1}$ describes the effective long-range interaction of the
pseudospin that with the pseudospins inside the sublattice, and parameter
$\nu_{2}$ \- with the pseudospins of the other sublattice.

The fourth term in (2.1) describes the interactions of pseudospins with the
external electric field:

$\displaystyle\hat{H}_{E}=-\sum\limits_{q}\mu_{y}E_{2}\left(\frac{\sigma_{q1}^{\text{A}}}{2}+\frac{\sigma_{q2}^{\text{A}}}{2}+\frac{\sigma_{q1}^{\text{B}}}{2}+\frac{\sigma_{q2}^{\text{B}}}{2}\right),$
(2.8)

where $\mu_{y}$ is y-component of effective dipole moments per one pseudospin.

The term $\hat{H}^{\prime}_{E}$ in Hamiltonian (2.1) takes into account the
dependence of effective dipole moments on the mean value of pseudospin
$s_{f}$:

$\displaystyle\hat{H}^{\prime}_{E}=-\sum\limits_{qf}s_{f}^{2}\mu^{\prime}E_{2}\frac{\sigma_{qf}}{2}=-\sum\limits_{qf}\left(\frac{1}{N}\sum\limits_{q^{\prime}}\sigma_{q^{\prime}f}\right)^{2}\mu^{\prime}E_{2}\frac{\sigma_{qf}}{2},$
(2.9)

where $\sigma_{qf}$ (f=1, 2, 3, 4) are a brief notation of pseudospins
$\sigma_{q1}^{\text{A}}$, $\sigma_{q2}^{\text{A}}$, $\sigma_{q1}^{\text{B}}$,
$\sigma_{q2}^{\text{B}}$, respectively. Here, we use corrections to dipole
moments $s_{f}^{2}\mu^{\prime}$ instead of $s_{f}\mu^{\prime}$ because of the
symmetry considerations, and the energy should not change when the field and
all pseudospins change their signs.

The term $\hat{H}^{\prime}_{E}$, as well as long-range interactions, are taken
into account in the mean field approximation:

$\displaystyle\hat{H}^{\prime}_{E}=-3\sum\limits_{q}\mu^{\prime}E_{2}\left(\frac{\eta_{1}^{2}\sigma_{q1}^{\text{A}}}{2}+\frac{\eta_{1}^{2}\sigma_{q2}^{\text{A}}}{2}+\frac{\eta_{2}^{2}\sigma_{q1}^{\text{B}}}{2}+\frac{\eta_{2}^{2}\sigma_{q2}^{\text{B}}}{2}\right)+2N(\eta_{1}^{3}+\eta_{2}^{3})\mu^{\prime}E_{2}.$
(2.10)

Our theory does not explicitly take into account the tunneling effects because
the tunneling of protons on the bonds and strains of the lattice being taken
into account simultaneously greatly complicates the calculations. In our
model, the tunneling parameter is partially taken into account in the form of
the other renormalized parameters, such as the interaction parameter $w^{0}$
and the effective dipole moment $\mu_{y}$. Explicit consideration of tunneling
could somewhat improve the agreement of the calculated characteristics with
the experimental data.

In the two-particle cluster approximation for short-range interactions, the
thermodynamic potential is given by:

$\displaystyle
G=NU_{\text{seed}}+NH^{0}+2N(\eta_{1}^{3}+\eta_{2}^{3})\mu^{\prime}E_{2}-Nv\sum\limits_{j}\sigma_{j}u_{j}$
$\displaystyle-k_{\text{B}}T\sum\limits_{q}\left\\{2\ln{\rm
Sp}\,\mathrm{e}^{-\beta\hat{H}^{(2)}_{q}}\left.-\ln{\rm
Sp}\,\mathrm{e}^{-\beta\hat{H}^{(1)\text{A}}_{q}}-\ln{\rm
Sp}\,\mathrm{e}^{-\beta\hat{H}^{(1)\text{B}}_{q}}\right\\}\right.,$ (2.11)

where $\beta=\frac{1}{k_{\text{B}}T}$, $k_{\text{B}}$ is Boltzmann constant,
$\hat{H}^{(2)}_{q}$, $\hat{H}^{(1)\text{A}}_{q}$, $\hat{H}^{(1)\text{B}}_{q}$
are two-particle and one-particle Hamiltonians:

$\displaystyle\hat{H}^{(2)}_{q}=-2w\left(\frac{\sigma_{q1}^{\text{A}}}{2}\frac{\sigma_{q2}^{\text{A}}}{2}+\frac{\sigma_{q1}^{\text{B}}}{2}\frac{\sigma_{q2}^{\text{B}}}{2}\right)-\frac{y_{1}}{\beta}\left(\frac{\sigma_{q1}^{\text{A}}}{2}+\frac{\sigma_{q2}^{\text{A}}}{2}\right)-\frac{y_{2}}{\beta}\left(\frac{\sigma_{q1}^{\text{B}}}{2}+\frac{\sigma_{q2}^{\text{B}}}{2}\right),$
(2.12)
$\displaystyle\hat{H}^{(1)\text{A}}_{q}=-\frac{\bar{y}_{1}}{\beta}\left(\frac{\sigma_{q1}^{\text{A}}}{2}+\frac{\sigma_{q2}^{\text{A}}}{2}\right),~{}~{}\hat{H}^{(1)\text{B}}_{q}=-\frac{\bar{y}_{2}}{\beta}\left(\frac{\sigma_{q1}^{\text{B}}}{2}+\frac{\sigma_{q2}^{\text{B}}}{2}\right),$
(2.13)

where such notations are used:

$\displaystyle
y_{1}=\beta\Delta_{1}+2\beta\nu_{1}\eta_{1}+2\beta\nu_{2}\eta_{2}+\beta(\mu_{y}E_{2}+3\eta_{1}^{2}\mu^{\prime}E_{2}),$
(2.14) $\displaystyle
y_{2}=\beta\Delta_{2}+2\beta\nu_{2}\eta_{1}+2\beta\nu_{1}\eta_{2}+\beta(\mu_{y}E_{2}+3\eta_{2}^{2}\mu^{\prime}E_{2}),$
$\displaystyle\bar{y}_{1}=\beta\Delta_{1}+y_{1},~{}~{}\bar{y}_{2}=\beta\Delta_{2}+y_{2}.$

Symbols $\Delta_{l}$ are the effective fields created by the neighboring bonds
from outside the cluster. From the condition of the minimum of the
thermodynamic potential $\partial G/\partial\Delta_{l}=0$, the system of
equations for these fields, as well as for the order parameters $\eta_{l}$ are
obtained:

$\displaystyle\eta_{1}=\frac{1}{D}\left[\sinh(y_{1}+y_{2})+\sinh(y_{1}-y_{2})+2a\sinh
y_{1}\right]=\tanh\frac{\bar{y}_{1}}{2},$
$\displaystyle\eta_{2}=\frac{1}{D}\left[\sinh(y_{1}+y_{2})-\sinh(y_{1}-y_{2})+2a\sinh
y_{2}\right]=\tanh\frac{\bar{y}_{2}}{2},$ (2.15)

Here, the following notations are used:

$\displaystyle D=\cosh(y_{1}+y_{2})+\cosh(y_{1}-y_{2})+2a\cosh y_{1}+2a\cosh
y_{2}+2a^{2},~{}~{}~{}~{}a=\mathrm{e}^{-\frac{w}{k_{\text{B}}T}}.$

Minimizing the thermodynamic potential with respect to the strains, an
additional system of equations for the strains $u_{j}$ was obtained in [25]:

$\displaystyle\sigma_{j}$ $\displaystyle=$ $\displaystyle
c_{j1}^{E0}u_{1}+c_{j2}^{E0}u_{2}+c_{j3}^{E0}u_{3}+c_{j5}^{E0}u_{5}-e_{2j}^{0}E_{2}-\frac{2\delta_{j}}{v}+\frac{4\delta_{j}}{vD}M$
(2.16) $\displaystyle-$
$\displaystyle\frac{1}{v}\psi_{j1}(\eta_{1}^{2}+\eta_{2}^{2})-\frac{2}{v}\psi_{j2}\eta_{1}\eta_{2}.$

where

$M=\bigl{[}a\,\cosh y_{1}+a\,\cosh y_{2}+2\,a^{2}\bigr{]}.$

In the case of ferroelectric ordering $\eta_{1}=\eta_{2}=\eta$,
$y_{l}=y_{2}=y$. Then, (2.14), (2.15) are given by:

$\displaystyle
y=\beta\Delta+2\beta(\nu_{1}+\nu_{2})\eta+\beta(\mu_{y}E_{2}+3\eta_{1}^{2}\mu^{\prime}E_{2}),~{}~{}~{}~{}\bar{y}=\beta\Delta+y,$
(2.17) $\displaystyle\eta=\frac{\sinh(2y)+2a\sinh y}{\cosh(2y)+1+4a\cosh
y+2a^{2}}=\tanh\frac{\bar{y}}{2}=\frac{\sinh y\cosh y+a\sinh
y}{\cosh^{2}y+2a\cosh y+a^{2}}=\frac{\sinh y}{\cosh y+a}.$ (2.18)

In the presence of hydrostatic pressure $\sigma_{1}=\sigma_{2}=\sigma_{3}=-p$,
$\sigma_{4}=\sigma_{5}=\sigma_{6}=0$. In [25], the expression for the
longitudinal component of polarization $P_{2}$ is also obtained, which in the
case of ferroelectric ordering is given by:

$\displaystyle
P_{2}=\sum\limits_{j}e_{2j}^{0}u_{j}+\chi_{22}^{u0}E_{2}+\frac{2\mu_{y}}{v}\eta+\frac{2\mu^{\prime}}{v}\eta^{3}.$
(2.19)

##  3 Dynamic dielectric properties of mechanically clamped CDP crystal.
Analytical results

Dynamic dielectric properties are studied in the absence of the electric field
and shear stresses. To study the dynamic properties of the CDP crystal, an
approach is used which is based on the ideas of Glauber’s stochastic model.
Based on this approach, we obtain a system of equations for time-dependent
unary correlation functions of pseudospins:

$-\alpha\frac{\mathrm{d}}{\mathrm{d}t}\langle\sigma_{qf}\rangle=\left\langle\sigma_{qf}\left[1-\sigma_{qf}\tanh\frac{1}{2}\beta{\varepsilon}_{qf}(t)\right]\right\rangle,$
(3.1)

where the parameter $\alpha$ sets the time scale of the dynamic processes in
the sytem; $\varepsilon_{qf^{\prime}}(t)$ is the local field acting on the
$f^{\prime}$-th pseudospin in the $q$-th cell; this is the factor at
$\sigma_{qf^{\prime}}/2$ in the initial Hamiltonian. In order to obtain a
closed system of equations for unary correlation functions, we use the two-
particle cluster approximation. In this approximation, the local fields
$\varepsilon_{qf}(t)$ are the coefficients at $\sigma_{qf}/2$ in two-particle
and one-particle Hamiltonians (2.13), (2.12). In the two-particle
approximation, they are given by:

$\displaystyle\varepsilon_{q1}=w\sigma_{q2}+\frac{y}{\beta},\quad\varepsilon_{q2}=w\sigma_{q1}+\frac{y}{\beta},$
(3.2)

and in the two-particle approximation they are given by:

$\displaystyle\varepsilon_{qf}=\frac{\bar{y}}{\beta}.$ (3.3)

Here, the effective fields $y$, $\bar{y}$ are given by expressions (2.14).

As a result of (3.1), we obtain a system of equations for the mean values of
pseudospins $\langle\sigma_{qf}\rangle=\eta_{f}$ in the two-particle
approximation

$\displaystyle\alpha\frac{\mathrm{d}}{\mathrm{d}t}\eta=-(1-P)\eta+L,$ (3.4)

and in the one-particle approximation

$\displaystyle\alpha\frac{\mathrm{d}}{\mathrm{d}t}\eta=-\eta+\tanh\frac{\bar{y}}{2},$
(3.5)

where the following notations are used:

$\displaystyle\\!\\!\\!\\!P=\frac{1}{2}\left[\tanh\left(\frac{\beta
w}{2}\\!+\\!\frac{y}{2}\right)\\!-\\!\tanh\left(-\frac{\beta
w}{2}\\!+\\!\frac{y}{2}\right)\right],~{}~{}L=\frac{1}{2}\left[\tanh\left(\frac{\beta
w}{2}\\!+\\!\frac{y}{2}\right)\\!+\\!\tanh\left(-\frac{\beta
w}{2}\\!+\\!\frac{y}{2}\right)\right].$

We limit ourselves to solving the equations (3.4) and (3.5) to the case of
small deviations from the equilibrium state. To do this, we present $\eta_{f}$
and the effective fields $y_{f}$, $\bar{y}_{f}$ as the sum of two terms –
equilibrium values and their deviations from the equilibrium state
(mechanically clamped crystal):

$\displaystyle\eta=\tilde{\eta}+\eta_{t},$ $\displaystyle
y=\tilde{y}+y_{t}=\beta[\Delta+2(\nu_{1}+\nu_{2})\tilde{\eta}+\Delta_{t}+2(\nu_{1}+\nu_{2})\tilde{\eta}_{t}+(\mu_{y}+3\tilde{\eta}_{1}^{2}\mu^{\prime})E_{2t}],~{}~{}E_{2t}=E_{2}\mathrm{e}^{\mathrm{i}\omega
t}.$ (3.6)

Here, $\Delta$ is the effective cluster field, and $\Delta_{t}$ is the
deviation from its equilibrium value, $\nu_{i}$ is the parameters of long-
range interactions.

We expand the coefficients $P_{f}$ and $L_{f}$ in series over $y_{ft}/{2}$
limited to linear terms:

$\displaystyle
P=P^{(0)}+\frac{y_{t}}{2}P^{(1)},~{}L=L^{(0)}+\frac{y_{t}}{2}L^{(1)},~{}\tanh\frac{\bar{y}}{2}=\tanh\frac{\tilde{\bar{y}}}{2}+\varphi\frac{\bar{y}_{t}}{2},$
(3.7)

where the following notations are used:

$\displaystyle
P^{(0)}=\frac{1-a^{2}}{Z},~{}~{}P^{(1)}=\frac{-4a(1-a^{2})\sinh\tilde{y}}{Z^{2}},~{}~{}~{}~{}~{}L^{(0)}=\frac{2a\sinh\tilde{y}}{Z},~{}L^{(1)}=\frac{4a[2a+(1+a^{2})\cosh\tilde{y}]}{Z^{2}},$
$\displaystyle
Z=1+a^{2}+2a\cosh\tilde{y};~{}~{}~{}~{}\varphi=(1-\tilde{\eta}^{2}),~{}~{}a=\mathrm{e}^{-\frac{w}{k_{\text{B}}T}},~{}~{}w=w^{0}+\sum\limits_{j}\delta_{j}u_{j},~{}~{}j=1,2,3,5.$

Substituting (3.6), (3.7) in (3.4), (3.5), we obtain a system of differential
equations for unary distribution functions:

$\displaystyle\alpha\frac{\mathrm{d}}{\mathrm{d}t}\eta_{t}=-\left(1-P^{(0)}-\frac{y_{t}}{2}P^{(1)}\right)(\tilde{\eta}+\eta_{t})+L^{(0)}+\frac{y_{t}}{2}L^{(1)},$
$\displaystyle\alpha\frac{\mathrm{d}}{\mathrm{d}t}\eta_{t}=-(\tilde{\eta}+\eta_{t})+\tanh\frac{\tilde{\bar{y}}}{2}+\varphi\frac{\bar{y}_{t}}{2}.$
(3.8)

Since the equilibrium state holds (2.18), it is easy to see that
$(1-P^{(0)})\tilde{\eta}=L^{(0)}$. Then, the system of equations (3.8) is
given by

$\displaystyle\alpha\frac{\mathrm{d}}{\mathrm{d}t}\eta_{t}=-\left(1-P^{(0)}\right)\eta_{t}+Y\frac{y_{t}}{2},~{}~{}~{}~{}~{}~{}Y=\left(P^{(1)}\tilde{\eta}+L^{(1)}\right),$
$\displaystyle\alpha\frac{\mathrm{d}}{\mathrm{d}t}\eta_{t}=-\eta_{t}+\varphi\frac{\bar{y}_{t}}{2}.$
(3.9)

Let us write it more in detail, taking into account (2.17):

$\displaystyle\alpha\frac{\mathrm{d}}{\mathrm{d}t}\eta_{t}=-\left(1-P^{(0)}\right)\eta_{t}+\frac{1}{2}\beta\left[\Delta_{t}+2(\nu_{1}+\nu_{2})\tilde{\eta}_{t}+(\mu_{y}+3\tilde{\eta}^{2}\mu^{\prime})E_{2t}\right]Y,$
$\displaystyle\alpha\frac{\mathrm{d}}{\mathrm{d}t}\eta_{t}=-\eta_{t}+\varphi\frac{1}{2}\beta\left[2\Delta_{t}+2(\nu_{1}+\nu_{2})\tilde{\eta}_{t}+(\mu_{y}+3\tilde{\eta}^{2}\mu^{\prime})E_{2t}\right].$
(3.10)

Excluding the parameter $\Delta_{ft}$, we obtain the equation for the time-
dependent order parameter:

$\displaystyle\frac{\mathrm{d}}{\mathrm{d}t}\eta_{t}=A\eta_{t}+B\beta(\mu_{y}+3\tilde{\eta}^{2}\mu^{\prime})E_{2t},$
(3.11)

where

$\displaystyle
A=\left[-\frac{(1-P^{(0)})2\varphi}{Y}+1+\varphi\beta(\nu_{1}+\nu_{2})\right]\Big{/}\left[\alpha\left(\frac{2\varphi}{Y}-1\right)\right],~{}~{}~{}~{}B=\frac{1}{2}\varphi\Big{/}\left[\alpha\left(\frac{2\varphi}{Y}-1\right)\right].$
(3.12)

Solving this equation, we obtain the time-dependent mean values of
pseudospins, and thus we find longitudinal dynamic dielectric susceptibility
of mechanically clamped CDP crystal:

$\displaystyle\chi_{22}(\omega)=\chi_{22}^{0}+\lim\limits_{E_{2t}\to
0}\frac{2}{v}(\mu_{y}+3\tilde{\eta}^{2}\mu^{\prime})\frac{\mathrm{d}\eta_{t}}{\mathrm{d}E_{2t}}=\chi_{22}^{0}+\frac{\chi}{1+\mathrm{i}\omega\tau},$
(3.13)

where $\chi$ is the static susceptibility of the pseudo-spin subsystem:

$\displaystyle\chi=\frac{2\beta(\mu_{y}+3\tilde{\eta}^{2}\mu^{\prime})^{2}}{v}\tau
B,$

$\tau$ is relaxation time:

$\tau=-1/A.$ (3.14)

Dynamic dielectric constant CDP is given by:

$\displaystyle\varepsilon_{22}(\omega)=1+4\piup\chi_{22}(\omega).$ (3.15)

##  4 Comparison of theoretical results with experimental data. Discussion of
the obtained results

The theory parameters are determined in [25] from the condition of agreement
of calculated characteristics with experimental data for temperature
dependences of spontaneous polarization $P_{2}(T)$ and dielectric permittivity
$\varepsilon_{22}(T)$ at different values of hydrostatic pressure [6],
spontaneous strains $u_{j}$ [22], molar heat capacity [13] and elastic
constants [27]; as well as the agreement with ab initio calculations of the
lattice contributions to molar heat capacity [18] and dielectric permittivity
at zero temperature [19].

It should be noted that the temperature dependences of the dielectric constant
$\varepsilon_{22}$ at different values of hydrostatic pressure were also
measured in [12]. However, they do not agree with experimental data [6]. It is
quite possible that another crystal sample was used there, which was grown
under different conditions. In addition, in [12] there are no data for the
temperature dependences of spontaneous polarization at different pressures, as
well as there are no data for dielectric characteristics at zero pressure.
Therefore, we used experimental data [6] to determine the model parameters. To
describe the dynamic dielectric permittivity, we used the data [28] because
they are measured at frequencies corresponding to the range of dispersion. In
the case of [12], the frequencies were much lower (up to 1MHz). At these
frequencies, the dielectric constant in the paraelectric phase behaves as
static, and in the ferroelectric phase there is a large contribution to the
permittivity associated with the reorientation of the domain walls. However,
our theory does not take into account the reorientation of the domain walls
and cannot describe the dielectric permittivity in the ferroelectric phase at
low frequencies.

Parameters of short-range interactions $w_{0}$ and long-range interactions
$\nu_{1}^{0}$ (“intra-sublattice”), $\nu_{2}^{0}$ (“inter-sublattice”) mainly
fix the phase transition temperature from paraelectric to ferroelectric phase
at the absence of external pressure and field, the order of phase transition
and the shape of curve $P_{2}(T)$. Their optimal values are:
$w_{0}/k_{\text{B}}=650$ K, $\nu_{1}^{0}/k_{\text{B}}=1.50$ K,
$\nu_{2}^{0}/k_{\text{B}}=0.23$ K.

To determine deformational potentials $\delta_{j}$ [see (2.4)] and $\psi_{j1}$
(2.7), $\psi_{j2}$ (2.7), it is necessary to use experimental data for the
shift of the phase transition temperature under hydrostatic and uniaxial
pressures as well as the data for temperature dependences of spontaneous
strains $u_{j}$, piezoelectric coefficients and elastic constants.
Unfortunately, only the data for the spontaneous strains and hydrostatic
pressure effect on the dielectric characteristics are available. As a result,
the experimental data for strains and dielectric characteristics can be
described using a great number of combinations of parameters $\psi_{j1}$,
$\psi_{j2}$. Therefore, for the sake of simplicity, we chose $\psi_{j2}$ to be
proportional to $\psi_{j1}$. Optimal values of deformational potentials are:
$\delta_{1}/k_{\text{B}}=1214$ K, $\delta_{2}/k_{\text{B}}=454$ K,
$\delta_{3}/k_{\text{B}}=1728$ K, $\delta_{5}/k_{\text{B}}=-131$ K;
$\psi_{11}/k_{\text{B}}=92.2$ K, $\psi_{21}/k_{\text{B}}=23.2$ K,
$\psi_{31}/k_{\text{B}}=139.7$ K, $\psi_{51}/k_{\text{B}}=5.5$ K; $\psi_{j2}$
= $\frac{1}{3}\psi_{j1}$.

The effective dipole moment in the paraelectric phase is found from the
condition of agreement of the calculated curve $\varepsilon_{22}(T)$ with
experimental data. We consider it to be dependent on the value of hydrostatic
pressure p, that is $\mu_{y}=\mu_{y}^{0}(1-k_{p}p)$, where
$\mu_{y}^{0}=2.63\cdot 10^{-18}$esu$\cdot$cm, $k_{p}=0.4\cdot
10^{-10}$cm2/dyn. The correction to the effective dipole moment
$\mu^{\prime}=-0.43\cdot 10^{-18}$esu$\cdot$cm is found from the condition of
the agreement of calculated saturation polarization with experimental data. In
[28], the experimental data for the dielectric permittivity $\varepsilon_{22}$
are almost twice as large as in [6]. This may be due to different growing
conditions of the samples and different features of measurement the dielectric
permittivity. Therefore, to describe the experimental data [28] for
$\varepsilon_{22}$, we assume that these samples are characterized by a more
effective dipole moment per pseudospin. Namely, $\mu_{y}^{0}=5.2\cdot
10^{-18}$esu$\cdot$cm, $\mu^{\prime}=-1.8\cdot 10^{-18}$esu$\cdot$cm.

The “seed” dielectric susceptibility $\chi_{22}^{u0}$, coefficients of
piezoelectric stress $e_{2j}^{0}$ and elastic constants $c_{ij}^{E0}$ are
found from the condition of the agreement of theory with experimental data in
the temperature regions far from the phase transition temperature $T_{c}$.
Their values are obtained as follows: $\chi_{22}^{u0}=0.443$; $e_{2j}^{0}=0.0$
esu/cm2; $c_{11}^{0E}=28.83\cdot 10^{10}$dyn/cm2, $c_{12}^{E0}=11.4\cdot
10^{10}$dyn/cm2, $c_{13}^{E0}=42.87\cdot 10^{10}$dyn/cm2,
$c_{22}^{E0}=26.67\cdot 10^{10}$dyn/cm2, $c_{23}^{E0}=14.5\cdot
10^{10}$dyn/cm2, $c_{33}^{E0}=65.45\cdot 10^{10}$dyn/cm2,
$c_{15}^{E0}=5.13\cdot 10^{10}$dyn/cm2, $c_{25}^{E0}=8.4\cdot 10^{10}$dyn/cm2,
$c_{35}^{E0}=7.50\cdot 10^{10}$dyn/cm2, $c_{55}^{E0}=5.20\cdot
10^{10}$dyn/cm2.

The volume of restricted primitive cell is $\upsilon$ = 0.467$\cdot 10^{-21}$
cm3 [7].

As can be seen from (3.13), the dynamic dielectric susceptibility is
determined by the behavior of the static dielectric susceptibility $\chi$ of
the pseudospin subsystem and the relaxation time $\tau$ in the system. Their
temperature dependences at different values of pressure are shown in figure 2.

a b

Figure 2: (Colour online) The temperature dependence of the inverse relaxation
time $\tau^{-1}$ and static dielectric susceptibility $\chi$ of the pseudospin
subsystem at different values of hydrostatic pressure. In the notation
$h_{p}$, the lower index means the value of pressure $p$ (kbar).

The calculated relaxation time $\tau$ tends to infinity at $T=T_{c}$. It is
associated with the relaxation frequency characteristic of this crystal (soft
relaxation mode) $\nu_{s}=(2\piup\tau)^{-1}$, which conditionally separates
the region of low-frequency and high-frequency dispersion. The inverse time
$\tau^{-1}$, as well as the relaxation frequency $\nu_{s}$ decrease while
approaching the phase transition temperature and tending to zero at a
temperature $T=T_{c}$.

At frequencies $\nu\ll\nu_{s}$, the real part of the dynamic dielectric
constant $\varepsilon_{22}^{\prime}$ behaves as static, and the imaginary part
$\varepsilon_{22}^{\prime\prime}$ is close to zero at all temperatures, except
for a narrow area near $T_{c}$. This can be seen in the frequency dependences
$\varepsilon_{22}(\nu)$ for different $\Delta T=T-T_{c}$ in the frequency
range $\nu<10^{7}$ Hz (figure 3, 4), as well as on the temperature dependences
$\varepsilon_{22}(T)$ at low frequencies and temperatures far from $T_{c}$
(figure 5).

Figure 3: (Colour online) The frequency dependences of the real
$\varepsilon^{\prime}_{22}$ and imaginary $\varepsilon^{\prime\prime}_{22}$
parts of the dielectric constant of CDP at different $\Delta T$ in
paraelectric phase and at different values of hydrostatic pressure $p$. In the
notation $h^{i}_{p}$, the superscripts number the temperature: $\Delta T=1.4$
K, $\circ$(1); $\Delta T=2.4$ K, $\square$(2); $\Delta T=6.3$ K,
$\triangledown$(3), and the lower indices indicate the value of pressure
(kbar). The symbols $\circ$, $\square$, $\diamond$ are experimental data [28].

Figure 4: (Colour online) The frequency dependences of the real
$\varepsilon^{\prime}_{22}$ and imaginary $\varepsilon^{\prime\prime}_{22}$
parts of the dielectric constant of CDP at different $\Delta T$ in
ferroelectric phase and at different values of hydrostatic pressure $p$. In
the notation $h^{i}_{p}$, the superscripts number the temperature: $\Delta
T=-1$ K (1); $\Delta T=-2$ K(2); $\Delta T=-5$ K(3), and the lower indices
indicate the value of pressure (kbar).

Figure 5: (Colour online) The temperature dependences of
$\varepsilon^{\prime}_{22}$ and $\varepsilon^{\prime\prime}_{22}$ of CDP at
different frequencies $\nu$ (MHz): 1.0 – 1,$\circ$,[28]; 12.0 –
2,$\square$,[28]; 72.4 – 3, $\diamond$, [28]; 251.2 – 4, $\vartriangle$, [28];
1000.0 – 5, $\triangledown$, [28]; 5000.0 – 6.

At frequencies $\nu\approx\nu_{s}$, there is a relaxation dispersion, which is
manifested in a rapid decrease of the real part of the dielectric constant
$\varepsilon_{22}^{\prime}$ with increasing frequency, as well as in large
values of the imaginary part $\varepsilon_{22}^{\prime\prime}$; in which case
the maximum value of $\varepsilon_{22}^{\prime\prime}$ corresponds to the
frequency $\nu_{s}$. This can be seen in the frequency dependences
$\varepsilon_{22}(\nu)$ for different $\Delta T=T-T_{c}$ in the frequency
range $10^{7}<\nu<10^{10}$ Hz (figure 3,4). At the temperature dependences
$\varepsilon_{22}(T)$, the relaxation dispersion is manifested in the decrease
of $\varepsilon_{22}^{\prime}$ near $T_{c}$ in the form of a depression, as
well as in the presence of the peak $\varepsilon_{22}^{\prime\prime}$ near
$T_{c}$ (figure 5). At frequencies $\nu\gg\nu_{s}$, the dielectric constant
behaves as a pure lattice contribution. This corresponds to the frequency
range $\nu>10^{10}$ Hz on the frequency dependences $\varepsilon_{22}(\nu)$ in
figure 3,4.

The increase in the relaxation time $\tau$ and the decrease in the relaxation
frequency $\nu_{s}$ as we approach the temperature $T=T_{c}$, is manifested in
the shift of the region of dispersion to lower frequencies on the frequency
dependence $\varepsilon_{22}(\nu)$ (figure 3,4), as well as in the expansion
of the “depression” of $\varepsilon_{22}^{\prime}$ and of the peak on
$\varepsilon_{22}^{\prime\prime}$ when approaching the temperature $T=T_{c}$
(figure 5). Since $\nu_{s}\rightarrow 0$ at $T=T_{c}$, this “depression” [as
well as the peak on $\varepsilon_{22}^{\prime\prime}(T)$] exists at all
frequencies; at low frequencies, it is very narrow. The value of permittivity
$\varepsilon_{22}^{\prime}$ at the minimum point (at $T=T_{c}$) is equal to
the lattice contribution $\varepsilon_{22}^{0}$.

Hydrostatic pressure $p_{h}$ lowers the phase transition temperature $T_{c}$
and reduces the value of $\tau^{-1}$ at a fixed $\Delta T=T-T_{c}$, (figure 2,
curves h1, h2, h3) compared to $\tau^{-1}$ without pressure (curve h0). As a
result, the dispersion region shifts to lower frequencies (figure 3,4, curves
h${}^{i}_{1}$, h${}^{i}_{2}$), compared to the case of no pressure (curves
h${}^{i}_{0}$). While on the temperature dependence of the permittivity, the
“depression” on $\varepsilon_{22}^{\prime}$ and the peak on
$\varepsilon_{22}^{\prime\prime}$ (figure 6, curves h${}^{i}_{1}$,
h${}^{i}_{2}$) broaden under pressure, compared with the case of no pressure
(curves h${}^{i}_{0}$).

Figure 6: (Colour online) The temperature dependences of
$\varepsilon^{\prime}_{22}$ and $\varepsilon^{\prime\prime}_{22}$ of CDP at
different frequencies $\nu$ and at different values of hydrostatic pressure
$p_{h}$. In the notation $h^{i}_{p}$, the superscripts number the frequency
$\nu$ (MHz): 72.4 – 3, $\diamond$, [28]; 251.2 – 4, $\vartriangle$, [28], and
the lower indices indicate the value of pressure (kbar).

In addition, at low frequencies in the paraelectric phase (figure 3) and at
temperatures far from $T_{c}$ the permittivity $\varepsilon_{22}^{\prime}$ in
figure 6 in the presence of pressure is greater than without pressure, because
it is close to static in this area. The static permittivity at $\Delta T=$
const increases with pressure (figure 2b). In the ferroelectric phase at low
temperatures, the pressure reduces the dielectric constant.

##  5 Conclusions

It is established that the dynamic dielectric constant of CDP at low
frequencies behaves as static; at frequencies commensurate with the inverse
relaxation time, there is a relaxation dispersion; and at high frequencies,
only the lattice contribution is manifested in permittivity. The region of the
longitudinal dispersion in CDP shifts to low frequencies when the temperature
approaches the phase transition point, which is associated with a significant
increase in relaxation time when approaching the temperature $T_{c}$. A
satisfactory agreement of theoretical results with experimental data is
obtained.

The effect of hydrostatic pressure on the dielectric properties is manifested
in a decrease of the phase transition temperature, in an increase of the
static dielectric constant in the paraelectric phase and in an increase of the
relaxation time. This leads to an increase of the dynamic permittivity at pre-
relaxation frequencies and to the shift of the dispersion region to lower
frequencies. In the ferroelectric phase at low temperatures, the pressure
reduces the dielectric constant. The dynamic permittivity is monodisperse.

As you can see, the dielectric constant in the dispersion region is quite
sensitive to the pressure. Thus, the CDP crystal can be used as a high-
frequency filter in which the absorption and transmission spectra can be
adjusted by pressure.

## References

  * [1] Matsunaga H., Itoh K., Nakamura E., J. Phys. Soc. Jpn., 1980, 48, No. 6, 2011–2014, doi:10.1143/JPSJ.48.2011.
  * [2] Itoh K., Hagiwara T., Nakamura E., J. Phys. Soc. Jpn., 1983, 52, No. 8, 2626–2629, doi:10.1143/JPSJ.52.2626.
  * [3] Iwata Y., Koyano N., Shibuya I., J. Phys. Soc. Jpn., 1980, 49, No. 1, 304–307, doi:10.1143/JPSJ.49.304.
  * [4] Iwata Y., Deguchi K., Mitani S., Shibuya I., Onodera Y., Nakamura E., J. Phys. Soc. Jpn., 1994, 63, No. 11, 4044–4050, doi:10.1143/JPSJ.63.4044.
  * [5] Yasuda N., Okamoto M., Shimizu H., Fujimoto S., Yoshino K., Inuishi Y., Phys. Rev. Lett., 1978, 41, No. 19, 1311–1314, doi:10.1103/PhysRevLett.41.1311.
  * [6] Yasuda N., Fujimoto S., Okamoto M., Shimizu H., Yoshino K., Inuishi Y., Phys. Rev. B, 1979, 20, No. 7, 2755–2764, doi:10.1103/PhysRevB.20.2755.
  * [7] Schuele P. J., Thoma R. A., Jpn. J. Appl. Phys., 1985, 24, Suppl. 24-2, 935–937, doi:10.7567/JJAPS.24S2.935.
  * [8] Gesi K., Ozawa K., Jpn. J. Appl. Phys., 1978, 17, No. 2, 435–436, doi:10.1143/JJAP.17.435.
  * [9] Brandt N. B., Zhukov S. G., Kulbachinskii V. A., Smirnov P. S., Strukov B. A., Fiz. Tverd. Tela, 1986,
28, 3159, (in Russian).

  * [10] Kobayashi Yu., Deguchi K., Azuma Sh., Suzuki E., Ming Li Ch., Endo Sh., Kikegawad T., Ferroelectrics, 2003, 285, No. 7, 83–89, doi:10.1080/00150190390205924.
  * [11] Magome E., Tomiaka S., Tao Y., Komukae M., J. Phys. Soc. Jpn., 2010, 79, No. 2,
025002 (2 pages), doi:10.1143/JPSJ.79.025002.

  * [12] Deguchi K., Azuma S., Kobayashi Y., Endo S., Tokunaga M., Phys. Rev. B, 2004, 69, No. 2,
024106 (5 pages), doi:10.1103/PhysRevB.69.024106.

  * [13] Imai K., J. Phys. Soc. Jpn., 1983, 52, No. 11, 3960–3965, doi:10.1143/JPSJ.52.3960.
  * [14] Shchur Ya., Phys. Rev. B, 2006, 74, 054301(8 pages), doi:10.1103/PhysRevB.74.054301.
  * [15] Shchur Ya., Phys. Status Solidi B, 2007, 244, No. 2, 569–577, doi:10.1002/pssb.200642176.
  * [16] Lasave J., Abufager P., Koval S., Phys. Rev. B, 2016, 93, No. 13, 134112 (11 pages),
doi:10.1103/PhysRevB.93.134112.

  * [17] Kojyo N., Onodera Y., J. Phys. Soc. Jpn., 1988, 57, No. 12, 439–4402, doi:10.1143/JPSJ.57.4391.
  * [18] Shchur Ya., Bryk T., Klevets I., Kityk A.V., Comput. Mater. Sci., 2016, 111, 301309,
doi:10.1016/j.commatsci.2015.09.014.

  * [19] Van Troeye B., van Setten M. J., Giantomassi M., Torrent M., Rignanese G.-M., Gonze X., Phys. Rev. B, 2017, 95, No. 2, 024112 (9 pages), doi:10.1103/PhysRevB.95.024112.
  * [20] Blinc R., Sa Baretto F. C., J. Chem. Phys., 1980, 72, No. 11, 6031–6034, doi:10.1063/1.439058.
  * [21] Stasyuk I. V., Levytsky R. R., Zachek I. R., Shchur Ya. Y., Kutny J. V., Miz E. V., Preprint of the Institute for Condensed Matter Physics, ICMP914R, Lviv, 1991, (in Russian).
  * [22] Deguchi K., Okaue E., Ushio S., Nakamura E., Abe K., J. Phys. Soc. Jpn., 1984, 53, No. 9, 3074–3080, doi:10.1143/JPSJ.53.3074.
  * [23] Levitskii R. R., Zachek I. R., Vdovych A. S., Phys. Chem. Solid State, 2012, 13, No. 1, 40–47.
  * [24] Levitskii R. R., Zachek I. R., Vdovych A. S., J. Phys. Stud., 2012, 16, No. 4, 4702,
doi:10.30970/jps.16.4702, (in Ukrainian).

  * [25] Vdovych A. S., Zachek I. R., Levitskii R. R., Condens. Matter Phys., 2020, 23, No. 3, 33702, doi:10.5488/CMP.23.33702.
  * [26] Vdovych A. S., Levitskii R. R., Zachek I. R., Moina A. P., J. Phys. Stud., 2021, 25, No. 3, 3702,
doi:10.30970/jps.25.3702.

  * [27] Prawer S., Smith T. F., Finlayson T. R., Aust. J. Phys., 1985, 38, No. 1, 63–83, doi:10.1071/PH850063.
  * [28] Deguchi K., Nakamura E., Okaue E., Aramaki N., J. Phys. Soc. Jpn., 1982, 51, No. 11, 3575–3582,
doi:10.1143/JPSJ.51.3575.

Ukrainian 0

Âïëèâ ãäðîñòàòèíîãî òèñêó íà äèíàìí äåëåêòðèí õàðàêòåðèñòèêè ñåãíåòîåëåêòðèêà
CsH2PO4 À. Ñ. Âäîâè ${}^{\textsf{\footnotesize{\mbox{?}}}}$, Ð. Ð. Ëåâèöüêèé
${}^{\textsf{\footnotesize{\mbox{?}}}}$, . Ð. Çàåê
${}^{\textsf{\footnotesize{\mbox{?}}}}$

${}^{\textsf{\footnotesize 1}}$ íñòèòóò ôçèêè êîíäåíñîâàíèõ ñèñòåì Íàöîíàëüíî¿
àêàäåì¿ íàóê Óêðà¿íè,
âóë. Ñâíööüêîãî, 1, 79011 Ëüââ, Óêðà¿íà ${}^{\textsf{\footnotesize 2}}$
Íàöîíàëüíèé óíâåðñèòåò ‘‘Ëüââñüêà ïîëòåõíêà’’, Óêðà¿íà, 79013, Ëüââ, âóë. Ñ.
Áàíäåðè, 12


# Going beyond skyrmions optically: Ultrafast generation of antiferromagnetic
meron-antimeron pairs

Sumit Ghosh<EMAIL_ADDRESS>Peter Grünberg Institut (PGI-1),
Forschungszentrum Jülich GmbH, 52428 Jülich, Germany  Institute of Physics,
Johannes Gutenberg-University Mainz, 55128 Mainz, Germany  Stefan Blügel
Peter Grünberg Institut (PGI-1), Forschungszentrum Jülich GmbH, 52428 Jülich,
Germany  Yuriy Mokrousov  Peter Grünberg Institut (PGI-1), Forschungszentrum
Jülich GmbH, 52428 Jülich, Germany  Institute of Physics, Johannes Gutenberg-
University Mainz, 55128 Mainz, Germany

###### Abstract

We propose a new mechanism to produce a meron-antimeron pair in a two
dimensional antiferromagnet with an ultrafast laser pulse via thermal
Schwinger mechanism. Unlike ultrafast skyrmion nucleation, this process
conserves total topological charge. We systematically show different stages of
the dynamics and define proper topological invariants to characterise the
configurations. The emergent structure can retain its topological structure
for up to 100 ps. By introducing a topological structure factor we show that
pair formation is robust against any random choice of initial magnetic
configuration and can survive against disorder. Our findings demonstrate that
the rich world of spin textures, which goes beyond conventional skyrmions, can
be reached optically.

In the quest for fast and efficient mechanism for manipulating magnetic order,
ultrafast optical manipulation of magnetism emerged one of the most promising
paradigms. Interaction between an ultrafast laser and magnetic moments takes
place at a much faster timescale and consumes far less energy Kirilyuk  _et
al._ (2010); Radu  _et al._ (2011) compared to electrical switching Yang  _et
al._ (2017). The underlying physical mechanism is, however, still in a mist.
From the initial experimental demonstration of ultrafast demagnetisation
Beaurepaire  _et al._ (1996) this effect was considered to be of thermal
origin conventionally modelled with the phenomenological three temperature
model. Zhang and Hübner Zhang and Hübner (2000) proposed a microscopic
mechanism mediated by spin-orbit coupling, whereas Koopmans and co-workers
proposed an alternative mechanism mediated by electron-phonon coupling
Koopmans  _et al._ (2005). These theoretical frameworks aim to understand the
transition from ordered magnetic states to strongly disordered and often
nonmagnetic states, which gives the impression that the governing mechanism
behind the ultrafast demagnetisation is of thermal origin Koopmans  _et al._
(2010). Based on time-dependent density functional theory approach, Dewhurst
et al. recently proposed an alternative mechanism Dewhurst  _et al._ (2018)
for transferring angular momentum between different sites,which has been
subsequently observed experimentally Willems  _et al._ (2020). This mechanism
relies on a coherent redistribution of spin angular momentum Dewhurst  _et
al._ (2021) and therefore can be exploited to drive a ferromagnet-to-
antiferromagnet transition Golias  _et al._ (2021).

In this context, an ability to generate any desired magnetic configuration,
i.e. a spin texture of desired properties, with an ultrafast optical laser,
remains a challenge. Recent experimental demonstration of skyrmion nucleation
with an ultrafast laser by Büttner et al. Büttner  _et al._ (2021) has proved
the feasibility of optical generation of spin textures. However, the proposed
explanation of the effect relies solely on classical magnetisation dynamics
which does not take into account electronic interactions. Besides, the working
model contains the anti-symmetric Dzyaloshinskii-Moriya interaction which is
known to drive the system into a skyrmionic state without any additional
excitation Fert  _et al._ (2013). Optical skyrmion generation has been also
demonstrated with a tight-binding electronic Hamiltonian with spin-orbit
coupling Viñas Boström  _et al._ (2022), which also makes its solely-optical
origin somewhat ambiguous. Such topological structures are known to arise
during the transition between two different magnetic configurations Heo  _et
al._ (2016); Ono and Ishihara (2019), however, their transient nature makes
them useless for any practical purpose. An alternative mechanism of laser
induced magnetisation dynamics was recently proposed by Ghosh and co-workers
Ghosh  _et al._ (2022), showing the emergence of a new chiral spin-mixing
interaction which can lead to chiral formation even in absence of any
intrinsic spin-orbit coupling. The emergent chirality is quite sensitive to
the laser parameters which emphasizes the importance of coherent interaction
of the laser with the electronic degrees of freedom for chiral state
generation. This mechanism leads to a quasi-stable chiral formation which can
survive for several picoseconds, and it is distinctly different from that
driven by thermal excitations.

Figure 1: Generation of a meron-antimeron pair by an ultrafast laser pulse.
(a) Initial antiferromagnetic configuration. (b) Emergence of a quasi-stable
meron-antimeron pair after the action of the laser pulse.

It is well-known that skyrmions are nontrivial magnetic configurations that
are equivalent to magnetic monopoles which can be characterised by a skyrmion
number Nagaosa and Tokura (2013). Due to their heavy mass Saha (1936),
creating isolated magnetic monopoles still remains quite challenging. An
alternative paradigm would lie in producing a monopole-antimonopole pair
Kobayashi (2021), which is also making its appearance in recent experimental
observations of magnetic textures. Indeed, topological spin textures often
emerge in pairs consisting of a texture and its anti-texture, keeping the
total topological charge to zero. Such pair formation has been observed both
at metallic interfaces Gao  _et al._ (2019) and predicted in two-dimensional
magnets Lu  _et al._ (2020); Augustin  _et al._ (2021). These pairs can be
further transmuted to a pair of different but opposite textures (for example
meron to skyrmion Yu  _et al._ (2018) conversion by an external magnetic
field) without violating the conservation of topological charge.

All these studies so far have been done on ferromagnetic materials which
naturally raises an alluring question $-$ is it possible to combine the
aforementioned ideas and stimulate such pair production in an antiferromagnet
via an optical excitation? Antiferromagnets are favorable materials due to
their natural abundance and faster response as compared to ferromagnets.
Besides, they are immune to any external magnetic field and also demonstrate
rich and topologically nontrivial electronic properties. However, regarding
the nontrivial real space textures in antiferromagnets, reported studies are
restricted to a narrow domain of skyrmions only Zhang  _et al._ (2016);
*Buhl2017; *Gobel2017. Here, we demonstrate that it is indeed possible to go
beyond skyrmions and nucleate meron-antimeron pairs in an antiferromagnet with
an ultrafast laser pulse. We present a systematic study of the different
stages of nucleation and demonstrate the gradual formation of meron-antimeron
pair starting from a collinear antiferromagnetic state. We analyze different
physical observables at different stages of the dynamics and define a proper
topological invariant to characterise their dynamical evolution. Finally we
reveal the robustness of the overall process versus disorder.

## Result and Discussion:

We define a two-dimensional antiferromagnet on a square lattice with a time
dependent double-exchange Hamiltonian Koshibae  _et al._ (2009); Ono and
Ishihara (2017); Ghosh  _et al._ (2022)

$\displaystyle H$ $\displaystyle=$
$\displaystyle-J\sum_{i;\mu,\nu}c^{\dagger}_{i,\mu}(\hat{\bm{m}}_{i}\cdot\bm{\sigma})_{\mu\nu}c_{j,\nu}-\sum_{i,j;\mu}\tau_{ij}(t)c^{\dagger}_{i,\mu}c_{j,\nu},\quad$
(1)

where $c^{\dagger}_{i,\mu},c_{i,\mu}$ are the creation (annihilation)
operators of an electron at site $i$ with spin $\mu$. $\hat{\bm{m}}_{i}$ is a
unit vector at site $i$ denoting the direction of the magnetisation and
$\bm{\sigma}$ is the vector of Pauli spin matrices. $J$ denotes the coupling
between the local magnetic moment and the itinerant moment, which, following
Hund’s rule is kept positive Ono and Ishihara (2017) and set to 1 eV.
$\tau_{ij}(t)$ is the time-dependent hopping parameter between site $i$ and
site $j$. The action of the laser is defined as a time varying electric field
with a Gaussian envelope modelled as a time varying vector potential
$\bm{A}(t)=-\frac{\mathcal{E}_{0}}{\omega}e^{-(t-t_{0})^{2}/2s^{2}}(\cos(\omega(t-t_{0}))\hat{\bm{x}}+\sin(\omega(t-t_{0}))\hat{\bm{y}})$,
where $\mathcal{E}_{0}$ is the peak amplitude of laser occurring at time
$t=t_{0}$ and $\omega$ is the angular frequency. $s$ denotes the broadening of
the pulse which we kept at 5 fs. In this work we keep $\mathcal{E}_{0}=0.1$
eV$\cdot a_{0}^{-1}$ ($a_{0}$ being the lattice constant), and
$\omega=2.02\,|J|/\hbar=3.07\times 10^{15}$ Hz. This is incorporated in the
Hamiltonian (Eq.1) via Pierel’s substitution $\tau_{ij}(t)=\tau
e^{i(e/\hbar)\bm{A}(t)\cdot\bm{r}_{ij}}$, where $r_{ij}$ is the vector
connecting site $i$ to site $j$. The static value of the hopping parameter is
kept at $\tau=0.4$ eV. For our study we consider a super-cell of dimension
$20\times 20$ with periodic boundary condition.

Figure 2: Time evolution of different physical observables. (a) Change of
occupation for different states with their colour denoting corresponding
eigenvalues. Here we show every 5th state. Gray shaded region shows the time
delay between the incident of laser and onset of magnetisation dynamics. (b)
Change of magnetisation vector (solid) and effective field (transparent)
components of the magnetic moment at $x$=1,$y$=10. Yellow line denotes the $x$
component of the laser electric field. (c) Change of average total energy
(gray dashed) and rate of change of energy (red) where the green portion shows
when the system absorbs energy. (d) Evolution of different components of
average torque.

The ground state is constructed by filling half of the lowest eigenvalue
states which in our case are all the states with negative eigenvalues. For our
choice of parameters, this corresponds to an antiferromagnetic alignment
Koshibae  _et al._ (2009); Ghosh  _et al._ (2022). We simulate the thermal
equilibrium at finite temperature by adding a small random fluctuation
(0.1$\pi$) to the polar angle while the azimuthal angle is chosen randomly
between 0 and 2$\pi$. This randomness is crucial to initiate the magnetisation
dynamics by generating initial spin-mixing interaction Ghosh  _et al._ (2022).
Note that our initial configuration does not posses any spin-orbit coupling or
any specific non-collinear order and therefore constitutes a trivial
antiferromagnetic insulator. At any instance of time $t$ any quantum state
$|\psi_{l}\rangle$ can be expressed as the linear combination of instantaneous
eigenstates $|n(t)\rangle$ of the Hamiltonian which is evolved within
Schrödinger picture Koshibae  _et al._ (2009); Ghosh  _et al._ (2022). This
allows us to evaluate the instantaneous effective field for any site $i$ as
$\frac{1}{\mu_{B}}\sum\langle\psi_{l}|-\nabla_{\bm{m}_{i}}H|\psi_{l}\rangle$,
$\mu_{B}$ being the Bohr magneton. This effective field is used in a set of
Landau-Lifsitz-Gilbert equations Ono and Ishihara (2017); Petrović  _et al._
(2018); Ghosh  _et al._ (2022)

$\displaystyle\frac{d\bm{m}_{i}}{dt}=-\gamma(\bm{m}_{i}\times\bm{B}_{i})-\lambda\bm{m}_{i}\times(\bm{m}_{i}\times\bm{B}_{i}),$
(2)

where $\gamma=\frac{g_{e}\mu_{B}}{\hbar\mu_{s}}\frac{1}{1+\alpha^{2}}$,
$\lambda=\frac{g_{e}\mu_{B}}{\hbar\mu_{s}}\frac{\alpha}{1+\alpha^{2}}$, with
$g_{e}$ being the gyromagnetic ratio, $\mu_{s}$ is the onsite magnetic moment
which we keep at 1 $\mu_{B}$ and $\alpha$ is the dimensionless damping
coefficient which we keep fixed at 0.2.

After the action of the ultrafast laser, one can observe three different
phases of dynamics originated from complex interactions between electronic and
magnetic degrees of freedom (Fig. 2). The peak amplitude of the laser occurs
at $t=25$ fs (Fig. 2b), marked by the left edge of the gray area in Fig. 2.
This triggers a redistribution of occupation of quantum states (Fig. 2a) which
in turns start building the torque. An initial randomness is imperative to
start this process which also determines how fast the system can respond to
the pulse Ghosh  _et al._ (2022). In our present setup after a latency period
of approximately 100 fs (right edge of gray region in Fig. 2) the magnetic
moments start re-organising themselves causing the fast relaxation phase Ghosh
_et al._ (2022) that lasts for approximately 0.5 ps (Fig. 2b). The system then
enters into the slow relaxation phase which eventually leads to the final
texture shown in Fig. 1. Note that the energy absorption takes place only
during the short duration defined by the Gaussian envelope of the pulse (green
segment in Fig. 2c). The steady state can be characterised by evaluating the
average torque Ghosh and Manchon (2019) defined as
$\langle\tau_{s}\rangle=1/N\sum_{ij}\nu_{ij}(\bm{m}_{ij}\times\bm{B}_{ij})$,
where the index $i,j$ correspond to the position, $\nu_{ij}=\pm 1$ is the
sublattice index, and $N$ is the total number of sites. Note that the optical
excitation predominantly produces the out of plane component of the torque
(Fig. 2d) leading to the final in-plane texture. Assuming that the average
energy gain per site is $\sim$ 0.1 eV with the dissipation rate being
$~{}10^{-6}$ eV fs-1, it would take the system about 100 ps to dissipate the
excess energy and regain its initial configuration.

Figure 3: Skyrmionic charge and winding number for antiferromagnetic meron-
antimeron pairs. (a) Distribution of skyrmionic charge is shown with red and
blue colour denoting positive and negative values, with winding centres
$\mathcal{W}_{1,2,3,4}$ marked with yellow spheres. Solid blue and red lines
show the contours along which winding number is calculated with red and blue
colours denoting a winding number $\pm$1\. (b) Time evolution of total
absolute skyrmion charge density and winding number.

## Dynamical evolution of topological charge:

From Fig. 1 one can see that the final configuration contains a pair of
antiferromagnetic meron and anti-meron. This is reminiscent of Schwinger
mechanism which produces a particle and anti-particle. Similar mechanism can
also generate monopole-antimonopole pairs Affleck and Manton (1982);
*Affleck1982a, that can survive even in presence of a thermal bath Gould  _et
al._ (2019). Such pair nucleation for topological defects was also predicted
in the context of bubble collision Digal and Srivastava (1996); *Copeland1996;
*Digal1997 following Kibble mechanism Kibble (1976). To characterise the pairs
we choose their skyrmionic charge. On a ferromagnetic background, each meron
can be characterised by a skyrmion number $\pm\frac{1}{2}$, resulting in a
total zero skyrmion number for the whole configuration. For an
antiferromagnet, one can define a skyrmion number for each sublattice Ono and
Ishihara (2019), however it does not take the inter-sublattice interaction
properly into account. To avoid this, we define a staggered skyrmionic charge
(Fig. 3)

$\displaystyle
S^{s}_{ij}=\nu_{ij}[\bm{m}_{ij}\cdot(\partial_{x}\bm{m}_{ij}\times\partial_{y}\bm{m}_{ij})],$
(3)

where $\nu_{ij}$=$\pm$1 for the opposite sublattices and
$\partial_{x}\bm{m}_{ij}=(\bm{m}_{i+1,j}-\bm{m}_{ij})$. Since the total
skyrmionic charge of the complete cell is zero, we use an absolute skyrmion
density $S^{s}_{Tot}=\sum|S^{s}_{ij}|$ to characterise the configuration. For
a well-isolated meron-antimeron pair this number should be 1. However due to
the finite size there is always some overlap between the different regions
causing a deviation of the total absolute skyrmion charge from 2 for the whole
super-cell (Fig. 3b). For a proper topological characterisation in this case,
we evaluate the winding number Rybakov (2021); *Kuchkin2020a for all four
meron/antimeron centres defined as

$\displaystyle\mathcal{W}=\frac{\mathcal{P}_{s}}{2\pi}\oint_{\mathcal{C}}d\bm{l}\cdot\bm{\nabla}\phi$
(4)

where $\mathcal{P}_{s}$ is the sub-lattice polarisation, given by
sgn$(m_{z})\cdot\nu$ with $m_{z}$ being the out of plane component of the
magnetisation at the centre and $\nu=\pm 1$ for the sub lattices. $\phi$ is
the azimuthal angle and the contour $\mathcal{C}$ encircles the winding centre
in anticlockwise direction passing through only one type of sublattice. We
choose four centres (Fig. 3a), each specifying a meron/antimeron and calculate
the corresponding winding number over the entire evolution period (Fig. 3).
The winding numbers, although initially showing some fluctuating non-integer
values originating from the lack of well defined winding centres during the
fast and slow relaxation process, eventually attain a steady magnitude of $\pm
1$ for $t\gtrsim 0.6$ ps. Compared to that, the total absolute staggered
skyrmion charge saturates slowly towards a magnitude of 2, indicating the
presence of four half instantons. It is worth mentioning that the winding
number remains unchanged under small changes in the winding center as long as
path encircles it. Therefore, the result remains consistent with the
skyrmionic charge distribution even when the winding centre differs from the
centre of topological charge.

Figure 4: Topological structure factor and its stability. (a) Topological
structure factor of the configuration shown in Fig. 3a. Inset shows the
corresponding real space distribution of the topological charge. (b) Average
topological structure factor averaged over 70 different configurations. (c)
Average topological structure factor for different strength of scalar impurity
averaged over 32 different configurations and scaled with its average value at
$V_{0}=0$. For each configuration we consider the mean value at
$q=0.1(\pm\hat{x}\pm\hat{y})\pi$.

## Topological structure factor:

The evolution of absolute staggered skyrmion charge and winding number ensures
the emergence of topological objects with opposite charges coming in pairs,
however it does not yet specify the distribution of the topological charge. To
determine that we define a topological structure factor

$\displaystyle\mathcal{S}_{\bm{q}}=\frac{1}{\mathcal{N}}\sum_{ij}S^{s}_{\bm{r}_{i}}S^{s}_{\bm{r}_{j}}e^{i\bm{q}\cdot(\bm{r}_{i}-\bm{r}_{i})}$
(5)

where $\mathcal{N}=\sum_{i}(S^{s}_{\bm{r}_{i}})^{2}$ is the normalisation
constant. One can readily see from Fig. 4a that the topological structure
factor is reminiscent of a checker-board pattern (Fig. 4a inset) which is
consistent with the skyrmionic charge distribution (Fig. 3). However, it is
also possible that the laser excited system moves into a spin-spiral
configuration or a combination of random clusters with no global long range
order. To verify that this is not the case, we consider 70 different
configurations with the same normalisation constant and evaluate the average
distribution of the topological structure factor (Fig. 4). The obtained maxima
in the distribution occur for $\bm{q}=0.1(\pm\hat{x}\pm\hat{y})\pi$ which
establishes the emergence of a meron-antimeron pair in our $20\times 20$ super
cell.

We finally study the robustness of these emergent configuration against
disorder. For that, we introduce random “impurities” modelled by on-site
energy deviations ranging from $-V_{0}$ to $V_{0}$, and measure the
topological structure factor at points $\bm{q}=0.1(\hat{x}\pm\hat{y})\pi$
(Fig.4). One can readily see that the structure manages to retain its basic
features for fairly strong values of the impurity potential, and then
gradually collapses into a featureless distribution. Note that a small random
impurity can enhance the formation as denoted by the initial rise in average
topological structure factor which is consistent with the behaviour against
small random fluctuation of magnetic moment. This establishes the robustness
of the pair formation process which would be helpful for its physical
realisation.

## Conclusions:

In this work we present a new paradigm for generating non-trivial magnetic
textures in antiferromagnets. We show that when excited by an ultrafast laser
a two-dimensional antiferromagnet can host meron-antimeron pairs via thermal
Schwinger mechanism. The topological texture can survive for up to 100 ps
which makes it possible to observe experimentally. The process is independent
of initial helicity of the laser and responsive within a moderate range of
frequency and peak amplitude of the laser. We define a staggered skyrmionic
charge, winding number and topological structure factor to demonstrate the
formation of topological charge and show that on average the system is more
likely to form a pair of meron and antimeron in a checker board pattern. This
formation is fairly robust against scalar impurity which makes it easier to
observe experimentally. Since the total topological charge of a texture /
anti-texture pairs is always zero, pairs of different genera can be connected
smoothly allowing more flexibility in exploring topological objects of higher
genus. Our results thus open a new route to generate out-of-equilibrium
topological features which remain inaccessible in their ground state
configuration. This takes us beyond the generation of conventional nontrivial
structures such as skyrmions and domain walls, and paves the way to generating
more exotic higher-order topological objects such as chiral bobbers Rybakov
_et al._ (2015) or hopfions Rybakov  _et al._ (2019).

We thank Nikolai Kiselev, Stefan Eisebitt, Frank Freimuth and Olena Gomonay
for discussions. We acknowledge financial support from Leibniz Collaborative
Excellence project OptiSPIN $-$ Optical Control of Nanoscale Spin Textures. We
acknowledge funding under SPP 2137 “Skyrmionics” of the DFG and financial
support from the European Research Council (ERC) under the European Union’s
Horizon 2020 research and innovation program (Grant No. 856538, project “3D
MAGiC”). The work was also supported also by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) $-$ TRR 173/2 $-$
268565370 (project A11), TRR 288 $-$ 422213477 (project B06). We also
gratefully acknowledge the Jülich Supercomputing Centre and RWTH Aachen
University for providing computational resources.

## References

  * Kirilyuk  _et al._ (2010) A. Kirilyuk, A. V. Kimel, and T. Rasing, Rev. Mod. Phys. 82, 2731 (2010).
  * Radu  _et al._ (2011) I. Radu, K. Vahaplar, C. Stamm, T. Kachel, N. Pontius, H. A. Dürr, T. A. Ostler, J. Barker, R. F. L. Evans, R. W. Chantrell, A. Tsukamoto, A. Itoh, A. Kirilyuk, T. Rasing, and A. V. Kimel, Nature 472, 205 (2011).
  * Yang  _et al._ (2017) Y. Yang, R. B. Wilson, J. Gorchon, C.-H. Lambert, S. Salahuddin, and J. Bokor, Sci. Adv. 3, 1 (2017), arXiv:1609.06392 .
  * Beaurepaire  _et al._ (1996) E. Beaurepaire, J.-C. Merle, A. Daunois, and J.-Y. Bigot, Phys. Rev. Lett. 76, 4250 (1996).
  * Zhang and Hübner (2000) G. P. Zhang and W. Hübner, Phys. Rev. Lett. 85, 3025 (2000).
  * Koopmans  _et al._ (2005) B. Koopmans, J. J. M. Ruigrok, F. D. Longa, and W. J. M. de Jonge, Phys. Rev. Lett. 95, 267207 (2005).
  * Koopmans  _et al._ (2010) B. Koopmans, G. Malinowski, F. Dalla Longa, D. Steiauf, M. Fähnle, T. Roth, M. Cinchetti, and M. Aeschlimann, Nat. Mater. 9, 259 (2010).
  * Dewhurst  _et al._ (2018) J. K. Dewhurst, P. Elliott, S. Shallcross, E. K. U. Gross, and S. Sharma, Nano Lett. 18, 1842 (2018).
  * Willems  _et al._ (2020) F. Willems, C. von Korff Schmising, C. Strüber, D. Schick, D. W. Engel, J. K. Dewhurst, P. Elliott, S. Sharma, and S. Eisebitt, Nat. Commun. 11, 871 (2020).
  * Dewhurst  _et al._ (2021) J. K. Dewhurst, S. Shallcross, P. Elliott, S. Eisebitt, C. v. K. Schmising, and S. Sharma, Phys. Rev. B 104, 054438 (2021).
  * Golias  _et al._ (2021) E. Golias, I. Kumberg, I. Gelen, S. Thakur, J. Gördes, R. Hosseinifar, Q. Guillet, J. K. Dewhurst, S. Sharma, C. Schüßler-Langeheine, N. Pontius, and W. Kuch, Phys. Rev. Lett. 126, 107202 (2021), arXiv:2006.16061 .
  * Büttner  _et al._ (2021) F. Büttner, B. Pfau, M. Böttcher, M. Schneider, G. Mercurio, C. M. Günther, P. Hessing, C. Klose, A. Wittmann, K. Gerlinger, L. M. Kern, C. Strüber, C. von Korff Schmising, J. Fuchs, D. Engel, A. Churikova, S. Huang, D. Suzuki, I. Lemesh, M. Huang, L. Caretta, D. Weder, J. H. Gaida, M. Möller, T. R. Harvey, S. Zayko, K. Bagschik, R. Carley, L. Mercadier, J. Schlappa, A. Yaroslavtsev, L. Le Guyarder, N. Gerasimova, A. Scherz, C. Deiter, R. Gort, D. Hickin, J. Zhu, M. Turcato, D. Lomidze, F. Erdinger, A. Castoldi, S. Maffessanti, M. Porro, A. Samartsev, J. Sinova, C. Ropers, J. H. Mentink, B. Dupé, G. S. Beach, and S. Eisebitt, Nat. Mater. 20, 30 (2021).
  * Fert  _et al._ (2013) A. Fert, V. Cros, and J. Sampaio, Nat. Nanotechnol. 8, 152 (2013).
  * Viñas Boström  _et al._ (2022) E. Viñas Boström, A. Rubio, and C. Verdozzi, npj Comput. Mater. 8, 62 (2022), arXiv:2010.16125 .
  * Heo  _et al._ (2016) C. Heo, N. S. Kiselev, A. K. Nandy, S. Blugel, and T. Rasing, Sci. Rep. 6, 1 (2016).
  * Ono and Ishihara (2019) A. Ono and S. Ishihara, J. Phys. Soc. Japan 88, 023703 (2019).
  * Ghosh  _et al._ (2022) S. Ghosh, F. Freimuth, O. Gomonay, S. Blügel, and Y. Mokrousov, Commun. Phys. 5, 69 (2022), arXiv:2011.01670 .
  * Nagaosa and Tokura (2013) N. Nagaosa and Y. Tokura, Nat. Nanotechnol. 8, 899 (2013).
  * Saha (1936) M. N. Saha, Indian J. Phys. 10, 141 (1936).
  * Kobayashi (2021) T. Kobayashi, Phys. Rev. D 104, 043501 (2021), arXiv:2105.12776 .
  * Gao  _et al._ (2019) N. Gao, S. G. Je, M. Y. Im, J. W. Choi, M. Yang, Q. Li, T. Y. Wang, S. Lee, H. S. Han, K. S. Lee, W. Chao, C. Hwang, J. Li, and Z. Q. Qiu, Nat. Commun. 10, 5603 (2019).
  * Lu  _et al._ (2020) X. Lu, R. Fei, L. Zhu, and L. Yang, Nat. Commun. 11, 4724 (2020), arXiv:2002.05208 .
  * Augustin  _et al._ (2021) M. Augustin, S. Jenkins, R. F. L. Evans, K. S. Novoselov, and E. J. G. Santos, Nat. Commun. 12, 185 (2021), arXiv:2012.03296 .
  * Yu  _et al._ (2018) X. Z. Yu, W. Koshibae, Y. Tokunaga, K. Shibata, Y. Taguchi, N. Nagaosa, and Y. Tokura, Nature 564, 95 (2018).
  * Zhang  _et al._ (2016) X. Zhang, Y. Zhou, and M. Ezawa, Sci. Rep. 6, 24795 (2016), arXiv:1504.01198 .
  * Buhl  _et al._ (2017) P. M. Buhl, F. Freimuth, S. Blügel, and Y. Mokrousov, Phys. status solidi - Rapid Res. Lett. 11, 1700007 (2017), arXiv:1701.03030 .
  * Göbel  _et al._ (2017) B. Göbel, A. Mook, J. Henk, and I. Mertig, Phys. Rev. B 96, 060406 (2017), arXiv:1707.05267 .
  * Koshibae  _et al._ (2009) W. Koshibae, N. Furukawa, and N. Nagaosa, Phys. Rev. Lett. 103, 266402 (2009).
  * Ono and Ishihara (2017) A. Ono and S. Ishihara, Phys. Rev. Lett. 119, 207202 (2017).
  * Petrović  _et al._ (2018) M. D. Petrović, B. S. Popescu, U. Bajpai, P. Plecháč, and B. K. Nikolić, Phys. Rev. Appl. 10, 054038 (2018).
  * Ghosh and Manchon (2019) S. Ghosh and A. Manchon, Phys. Rev. B 100, 014412 (2019), arXiv:1901.08314 .
  * Affleck and Manton (1982) I. K. Affleck and N. S. Manton, Nucl. Phys. B 194, 38 (1982).
  * Affleck  _et al._ (1982) I. K. Affleck, O. Alvarez, and N. S. Manton, Nucl. Phys. B 197, 509 (1982).
  * Gould  _et al._ (2019) O. Gould, S. Mangles, A. Rajantie, S. Rose, and C. Xie, Phys. Rev. A 99, 052120 (2019), arXiv:1812.04089 .
  * Digal and Srivastava (1996) S. Digal and A. M. Srivastava, Phys. Rev. Lett. 76, 583 (1996).
  * Copeland and Saffin (1996) E. J. Copeland and P. M. Saffin, Phys. Rev. D 54, 6088 (1996).
  * Digal  _et al._ (1997) S. Digal, S. Sengupta, and A. M. Srivastava, Phys. Rev. D 56, 2035 (1997).
  * Kibble (1976) T. W. B. Kibble, J. Phys. A. Math. Gen. 9, 1387 (1976).
  * Rybakov (2021) P. Rybakov,  _Topological excitations in field theory models of superconductivity and magnetism_ , Ph.D. thesis (2021).
  * Kuchkin  _et al._ (2020) V. M. Kuchkin, B. Barton-Singer, F. N. Rybakov, S. Blügel, B. J. Schroers, and N. S. Kiselev, Phys. Rev. B 102, 144422 (2020), arXiv:2007.06260 .
  * Rybakov  _et al._ (2015) F. N. Rybakov, A. B. Borisov, S. Blügel, and N. S. Kiselev, Phys. Rev. Lett. 115, 117201 (2015).
  * Rybakov  _et al._ (2019) F. N. Rybakov, N. S. Kiselev, A. B. Borisov, L. Döring, C. Melcher, and S. Blügel, Unpublished , 1 (2019), arXiv:1904.00250 .