∎

e1Corresponding Author<EMAIL_ADDRESS>
11institutetext: Department of Physics, The University of South Dakota,
Vermillion, SD 57069, USA 22institutetext: School of Physics and Astronomy,
University of Minnesota, Minneapolis, MN 55455, USA 33institutetext:
Department of Physics and Astronomy, Texas A$\&$M University, College Station,
TX 77843, USA

# Observation of Time-Dependent Internal Charge Amplification in a Planar
Germanium Detector at Cryogenic Temperature

P. Acharyaaddr1  M. Frittsaddr2 D.-M. Meie1,addr1  V. Mandicaddr2  C.-J.
Wangaddr1  R. Mahapatraaddr3 and M. Plattaddr3

(Received: date / Accepted: date)

###### Abstract

For the first time, time-dependent internal charge amplification through
impact ionization has been observed in a planar germanium (Ge) detector
operated at cryogenic temperature. In a time period of 30 and 45 minutes after
applying a bias voltage, the charge energy corresponding to a baseline of the
59.54 keV $\gamma$ rays from a 241Am source is amplified for a short period of
time and then decreases back to the baseline. The amplification of charge
energy depends strongly on the applied positive bias voltage with drifting
holes across the detector. No such phenomenon is visible with drifting
electrons across the detector. We find that the observed charge amplification
is dictated by the impact ionization of charged states, which has a strong
correlation with impurity level and applied electric field. We analyze the
dominant physics mechanisms that are responsible for the creation and the
impact ionization of charged states. Our analysis suggests that the
appropriate level of impurity in a Ge detector can enhance charge yield
through the impact ionization of charged states to achieve extremely low-
energy detection threshold ($<$ 10 meV) for MeV-scale dark matter searches if
the charge amplification can be stabilized.

###### Keywords:

Charge Collection Efficiency (CCE); Time-Dependent Internal Charge
Amplification; Cluster Dipole States

††journal: Eur. Phys. J. C

##  1 Introduction

Dark matter (DM) is believed to be ubiquitous. It makes up 85% of the mass of
the universe fzw ; ghi ; jlf ; mwg . Many candidates including axions, low-
mass DM, and weakly interacting massive particles (WIMPs) have been postulated
axion ; lmdm ; wimp . WIMPs are expected to generate observable nuclear recoil
energy through elastic scattering off nuclei rjg . Despite great efforts made
in searching for axions and WIMPs ardm ; cd09 ; cd ; cd1 ; cd2 ; cog ; cre12 ;
cou ; bar ; dama ; dar ; dri ; ede ; kim ; lux ; pan ; pic ; cd14 ; xe11 ;
xe15 ; xe17 ; xma ; zep , DM remains undetected. Recently, low-mass DM in the
MeV-scale has become an exciting DM candidate ess2012 ; ess2016 ; ho ; ste .
To directly detect MeV-scale DM, a detector with sensitivity of measuring a
single electron-hole (e-h) pair is required, since the energy deposition
induced by MeV-scale DM through elastic scattering off electrons or nuclei is
in the range of sub-eV to 100 eV ess2012 ; mei . In 2018, Mei et al. proposed
to detect MeV-scale DM utilizing germanium internal charge amplification
(GeICA) geia  for the charge created by the ionization of impurities mei .
GeICA can potentially achieve a detection energy threshold of $\sim$0.1 eV
(100 meV), allowing a large portion of both electronic recoils and nuclear
recoils in the range of sub-eV to 100 eV induced by MeV-scale DM to be
accessible mei .

GeICA amplifies internal charge through impact ionization, which is a process
first observed in Ge diodes a few decades ago imp1 ; imp . In this process, a
charge carrier, electron or hole, with sufficient kinetic energy can knock a
bound electron out of its valence state and elevate it to a state in the
conduction band, creating an electron-hole pair. Carriers gain sufficient
kinetic energy through applying a strong electric field. Impact ionization of
Ge atoms requires higher electric field than that of impurities, since the
bandgap of Ge is about 0.7 eV (700 meV) while the ionization potential of
neutral impurity atoms in Ge is around 0.01 eV (10 meV). When a Ge detector is
cooled down to below 10 K, the residual impurities in Ge start to freeze out
from the conduction or valence band to localized states mei2022evidence . At
below 6 K, the localized states become thermally stable and form electric
dipole states mei2022evidence , which are excited neutral impurity states with
a binding energy of less than 10 meV. The dipole states can trap charge to
form cluster dipole states mei2022evidence  with even smaller binding energy
depending on the operational temperature.

The formation of excited dipole states and cluster dipole states in p-type Ge
is depicted in Figure 1 mei2022evidence . The phase space of an immobile
negative ion for trapping charge carriers is smaller than that of movable
bound holes, whose motion is restricted by the Onsager radius,
$R=\frac{1}{4\pi\varepsilon\varepsilon_{0}K_{B}T}$, where $\varepsilon$= 16.2
is the relative permittivity for Ge, $\varepsilon_{0}$ is the permittivity of
free space, $K_{B}$ is the Boltzmann constant, and $T$ is temperature.
Therefore, the probability of forming $A^{-^{*}}$ states is higher than
forming $A^{+^{*}}$ states in a p-type detector. This is why electrons are
trapped more severely than holes in a p-type detector.

Figure 1: Shown are the processes involved in the formation of the excited
dipole states and the cluster dipole states in a p-type Ge detector operated
at low temperatures, where $\vec{p}$ and $\vec{q}$ are the corresponding
dipole moments.

##  2 The Experimental Methods and the Observed Physical Phenomenon

The impact ionization of impurities in Ge specimens has been reported for a
range of temperatures ($\sim$ 4.2 K to 298 K) by many authors from 1950s -
1970s sclar ; pickin ; smith ; palm . The most recent impact ionization of
impurities was reported by Phipps et al. with SuperCDMS-style detectors at 40
milliKelvin (mK) phipps1 ; phipps2 , which is similar to this work with a
detector made from a USD crystal wang1 ; wang2 . The detector was fabricated
at Texas A $\&$ M University with four channels for charge readout,
geometrically similar to SuperCDMS style detectors supercdms , as shown in
Figure 2.

Figure 2: Shown is the detector studied in this work with four electrodes
read-out ($Q_{1},Q_{2},Q_{3},andQ_{4}$). Four ${}^{241}Am$ sources were
arranged in such a way that one source was over each channel. The opposite
side of the detector is grounded with uniform Al electrodes so that the
electric field points along the z-axis of the crystal. The size of the
detector is 10 cm in diameter and 3.3 cm in thickness with a mass of $\sim$1.4
kg.

The detector was wire-bonded, mounted in a dilution refrigerator and tested at
the K100 Detector Testing Facility at the University of Minnesota (UMN).

Figure 3: Shown is the orientation of the Am-241 source mover before it was
installed.

The detector was cooled down to $\sim$40 mK and then tested in two separate
refrigerator runs, Run 67 in 2018 and Run 74 in 2021. In Run 67 four 241Am
sources were mounted above each channel on the detector (see Figure 2). The
lead collimators with 0.2 mm holes allow 59.54 keV $\gamma$ rays through to
the detector. Alpha particles from the sources were blocked by the source
encapsulation. We observed 59.54 keV peaks in spectra from each channel and
performed different measurements over a course of two weeks. In Run 74 a
single 241Am source was mounted on a carriage that could be moved by a
superconducting stepper motor (see Figure 3). The source was of a different
design with a 0.5 mm collimator hole and a lower rate of $\gamma$ rays
incident on the detector (about 75% as much as the sources used in Run 67).

Figure 4: Shown is the time-dependent charge response from Q4 at positive
biases for Run 67. The baseline of 59.54 keV gamma ray from Q4 initially rose
linearly for a few minutes, then quickly transitioned to an $\sim$exponential
fall off over 10s of minutes. Figure 5: Shown is the time-dependent charge
response from Q4 at positive biases for Run 74. The baseline of 59.54 keV
gamma ray from Q4 initially rose linearly for a few minutes, then quickly
transitioned to an $\sim$exponential fall off over 10s of minutes.

##  3 Experimental Data Analysis and Results

We used the data measured at various bias voltages to characterize the charge
collection efficiency (CCE), defined as the ratio of the measured charge
energy to 59.54 keV gamma rays, as a function of bias for each channel, and
studied how the CCE varies with time under a given bias voltage. This work
presents the results of an analysis of the 59.54 keV calibration line
behavior, as shown in Figure 4 & Figure 5. The observation in particular had a
rich set of time dependent behaviors, which are the amplitude of the signals
from the 59.54 keV line to initially rise with time after bias, then fall to
some steady-state value. This phenomenon was only present with positive
biases. Since the 241Am sources are positioned at the biased side of the
crystal and the mean free path of 59.54 keV gamma rays in Ge is about 0.09 cm,
this means that this phenomenon is present for events in which electrons are
collected almost immediately and holes are drifted through the full detector
thickness. This indicates that the observed behavior has a polarity dependence
for all events.

In Run 67 the observed phenomenon was much stronger in the center channel (Q4)
and weakest in the outer ring channel (Q1), suggesting a radial dependence. In
Run 74, with only one source present with a weaker event rate, the radial
dependence is not evident, indicating an overall event rate dependence. In
this work, we only study the 59.54 keV line in the center channel (Q4) at
positive biases, since this channel outputs the clearest results.

The observed features of the time-dependent charge collection for the 59.54
keV line can be summarized as: (i) time-dependent impact ionization of
drifting holes kicks off at a positive bias of 4.5 volts (1.36 V/cm) and a
mixing chamber (MC) temperature in a range from (30-35) mK during Run 67 and
no such a phenomenon is observed with a negatively biased detector when
drifting electrons; similarly, we observed the time-dependent impact
ionization phenomenon for holes in Run 74 kicking off at a positive bias
around 7 volts (2.12 V/cm) and a MC temperature in a range from (35-40) mK
(ii) the hole impact ionization initially increases linearly with a rate
dependent on the bias and overall event rate; and (iii) after reaching a peak
the increased charge signal falls exponentially with a fall time dependent on
bias. The loss of both types of signal is easily explained as the loss of CCE
due to the breakdown of bulk field from trapped charges. The linear increase
of the hole impact ionization signal and lack of electron impact ionization
indicate that the population of possible hole impact ionization sites starts
small and is created by drifting electrons which are captured. This suggests a
possible physical mechanism involving the combination of the following three
processes: $e^{-}+A^{0^{*}}\rightarrow A^{-^{*}}$, $h^{+}+A^{-^{*}}\rightarrow
e^{-}+2h^{+}+A^{-}$, and $h^{+}+A^{-}\rightarrow A^{0^{*}}$. The first process
represents the creation of cluster dipole states ($A^{-^{*}}$) through
drifting electrons across the detector induced by background radiation as
shown in Figure 1; the second stands for hole impact ionization of cluster
dipole states; and the third is trapping of $h^{+}$, which determines the loss
of neutralization at the end of series.

There is also a possible process, $E_{ph}+A^{-^{*}}\rightarrow
e^{-}+A^{0{*}}$, which is the impact ionization of cluster dipole states
through absorbing Neganov-Luke phonons luke , where $E_{ph}$ is the energy of
phonons created by drifting holes under a given electric field. This process
is not a significant contribution in the experiment since the Neganov-Luke
phonons have energy smaller than 1 meV mei . However, the creation of cluster
dipole charge states ($e^{-}+A^{0^{*}}\rightarrow A^{-^{*}}$) and the impact
ionization of cluster dipole states ($h^{+}+A^{-^{*}}\rightarrow e^{-}$ \+
2$h^{+}+A^{-}$) inside the detector generate a dynamic electric field, which,
in turn, impacts charge transport and charge creation. As a result, in this
model the observed time-dependent impact ionization of 59.54 keV $\gamma$ rays
involves the growth of charge states as a function of time, the impact
ionization of time-dependent cluster dipole states, and the loss of CCE due to
the breakdown of bulk field from trapped charges, which is again a function of
time.

Thus, the analysis presented here attempts to quantify the observed behavior
in an empirical way. This will help lend insight into likely physical models
that dominate the observed behavior. We assume that the detected charge
energy, $E(t)$, is related to the input 59.54 keV $\gamma$ rays through the
following equation:

$E(t)=E_{\gamma}\\{p_{0}+p_{1}exp[\frac{p_{2}}{p_{3}}(1-exp(-p_{3}t))]\\}exp(-p_{4}t),$
(1)

where $E_{\gamma}$ = 59.54 keV. Other terms in Eq. 1 are explained below:

(1) $p_{0}+p_{1}$ represents the average CCE ($\epsilon_{0}$) at $t$ = 0. If
it is greater than 1, it means that charge energy is gained through impact
ionization. It is expected that the reaction, $h^{+}+A^{-^{*}}\rightarrow
2h^{+}+e^{-}+A^{-}$, dominates the hole impact ionization at $t$=0 in the
range of the applied field sund . The parameter $p_{0}+p_{1}$ is the direct
measurement of CCE at $t$ = 0 in equation 1, as depicted in Figure 4 & Figure
5.

(2) $\\{p_{0}+p_{1}exp[\frac{p_{2}}{p_{3}}(1-exp(-p_{3}t))]\\}exp(-p_{4}t)$ is
the average CCE ($\epsilon_{t})$ for the charge created by the impact
ionization at $t>0$. We can write $\epsilon_{t}$ = $\epsilon_{0}\times
M(t)\times d(t)$, where M(t) stands for the charge gained through time-
dependent impact ionization while $d(t)$ is a time-dependent charge damping
factor that describes the charge trapping due to a dynamic process described
in points (3) and (4) which creates more charge states. Note that
$\epsilon_{0}$, $d(t)$, and $M(t)$ are correlated with applied electric field.

(3) $p_{2}$ represents the rate of creating $A^{-^{*}}$ states. It depends on
the density of the dipole state, the overall event rate, the drift velocity,
and the charge trapping cross section. Since both the drift velocity and the
charge trapping cross section are field dependent, it is expected $p_{2}$ has
a correlation with applied electric field.

(4) $p_{3}$ represents the rate of decreasing $A^{-^{*}}$ states, which is
proportional to the applied field and the overall event rate.

(5) $exp(-p_{4}t)$ represents the loss of neutralization. If we let
$p_{4}=1/\tau$, the factor $exp(-t/\tau)$ describes the collection of charge
created by impact ionization falling exponentially as a function of time
depending on bias. The parameter $\tau$ measures the effective time constant
for the loss of charge signal due to the loss of CCE resulting from the
breakdown of bulk field from trapped charges.

Using this empirical model expressed in Eq. 1, we fit the data from Run 67 and
Run 74. The observed trends and the fits are shown for different applied bias
voltages in Figures 6 and 7. Since the two runs were taken at two different
periods at which the tower temperature of detector housing was slightly
different. Thus, we expect that the detector responds slightly different to
the impact ionization as it depends on the detector conditions including the
thermal effect. For accurate reading in a temperature, the mixing chamber (MC)
temperature was taken. The MC temperature in Run 67 was in a range between
(30-35) mK for the data series shown in Figure 6. Similarly in Run 74, the MC
temperature was in a range between (35-40) mK for the data series shown in
Figure 7, being higher compared to Run 67.

At a given bias the observed phenomenon is stronger in Run 67 than in Run 74.
This can be attributed to differences in the overall event rate in the
detector between the two runs. By introducing an additional external source it
was observed that the impact ionization effect seen in 59.54 keV events was
stronger when the overall event rate in the detector was increased. In Run 74
the overall rate was lower because the 241Am source used produced events at a
lower rate, only one source was used rather than four, and a lead shield was
erected around the cryostat which reduced the background rate from
radioactivity in the lab environment.

Figure 6: Shown is the detected charge energy created by impact ionization as
a function of applied bias voltage fitted by Eq. 1 using data from Run 67. The
error bars represent the vertical width of the band shown in Figure 4.  Figure
7: Shown is the detected charge energy created by impact ionization as a
function of applied bias voltage from Run 74 fitted by Eq. 1. The error bars
represent the vertical width of the band shown in Figure 5.

We summarize the fitting parameters in Table 1. The fitting to the impact
ionization curve using Eq. 1 shows that all the parameters have a good
correlation with the electric field. This indicates that the parameters have
their own effect under the electric field. The variation of those parameters
with the electric field is shown in Figure 8. Since each fitting parameter is
a combination of at least two physical processes, it is difficult to figure
out the impact of the applied field on any single process.

Table 1: A summary of the fit parameters using equation 1. Note that the
values of $p_{2}$ and $p_{3}$ are omitted for 5.4 V in Run 74 due to the rate
of creating and decreasing $A^{-*}$ states was very small, leading to abnormal
values for $p_{2}$ and $p_{3}$, which do not follow the trend observed at
higher bias voltages, where the dependency on the applied bias field is
evident. Run 67
---
Bias |  $p_{0}$ |  $p_{1}$ |  $p_{2}$ |  $p_{3}$ |  $p_{4}$
4.5 V |  1.05 |  0.015 |  0.67 |  0.22 |  0.0080
5.0 V |  1.07 |  0.018 |  0.97 |  0.29 |  0.013
5.5 V |  1.09 |  0.019 |  1.31 |  0.34 |  0.020
6.0 V |  1.13 |  0.020 |  1.41 |  0.35 |  0.026
7.0 V |  1.22 |  0.028 |  1.80 |  0.44 |  0.039
8.0 V |  1.28 |  0.032 |  2.79 |  0.67 |  0.046
9.0 V |  1.39 |  0.038 |  3.98 |  0.94 |  0.056
Run 74
Bias |  $p_{0}$ |  $p_{1}$ |  $p_{2}$ |  $p_{3}$ |  $p_{4}$
5.4 V |  0.65 |  0.024 |  - |  - |  0.0058
6.4 V |  0.78 |  0.026 |  0.67 |  0.26 |  0.0098
6.9 V |  0.92 |  0.027 |  0.80 |  0.27 |  0.012
7.2 V |  0.95 |  0.029 |  1.002 |  0.32 |  0.013
7.4 V |  0.96 |  0.031 |  1.13 |  0.35 |  0.015
9.4 V |  1.22 |  0.033 |  2.50 |  0.63 |  0.032
Figure 8: Shown are the fitting parameters, $p_{0}$ to $p_{4}$ from equation 1
as a function of the applied field, $E$, and fits a linear regression model,
which demonstrates the correlation of the parameters with the electric field
using data from each run. The fitting functions for those parameters in Run 67
(solid lines fitted to the circles) are: $p_{0}$=0.25$E$+0.69,
$p_{1}$=0.017$E$-0.0087, $p_{2}$=2.27$E$-2.57, $p_{3}$=0.49$E$-0.49, and
$p_{4}$=0.036$E$-0.039. Similarly, the fitting functions for those parameters
in Run 74 (dotted lines fitted to the squares ) are: $p_{0}=0.47E-0.11$,
$p_{1}=0.0078E+0.011$, $p_{2}=2.1E-3.5$, $p_{3}=0.43E-0.61$, and
$p_{4}=0.022E-0.033$.

##  4 Discussion and Perspective of the Physics Application

As can be seen from Figure 8, although the parameters, $p_{0}$, $p_{1}$,
$p_{2}$, $p_{3}$, and $p_{4}$ obtained from Run 67 and Run 74 have similar
tendencies as a function of applied bias field, the values of these parameters
differ between the two runs from a minimum of a few percent to a maximum of
$\sim$40%, depending on the applied bias field. The cause of this difference
is likely due to the overall event rates used in the two runs. Run 67 has a
higher event rate compared to Run 74.

For bias 4.5 V and above in Run 67, the value of $p_{0}+p_{1}$ is greater than
1.0, indicating impact ionization at $t=0$. In Run 74 the bias required to
produce impact ionization at $t=0$ was greater, 9.4 V. The difference between
the runs is apparently attributable to the higher event rate in Run 67
mentioned previously.

For stable operation at which the average CCE is not a function of time for a
given electric field, the parameters $p_{3}$ and $p_{4}$ from the empirical
model (Eq.1) should be zero. Using the fitted functions obtained from Run 67,
$p_{3}$ = 0 and $p_{4}$ = 0 at electric fields less than $\sim$1 V/cm. At such
a small electric field, the measured CCE, $p_{0}$ \+ $p_{1}$ = 0.9483 and is
constant in time. This means that at this electric field, the detector is a
normal detector, with no visible impact ionization. The charge collection
efficiency is consistent with typical SuperCDMS detectors operated at a
similar electric field phipps2  and suggests that detectors made from USD-
grown crystals are suitable as SuperCDMS-style detectors in terms of charge
collection.

At a higher electric field, for example at the applied bias of 9 V, the CCE is
about 1.4 at $t=0$ from Run 67. This means a higher electric field is needed
to generate impact ionization. If one assumes the charge collection time per
event is less than 2 microseconds in this large-size detector, we can predict
the CCE as a function of applied electric field using the fitted parameters
from Run 67 and the empirical model (Eq. 1). Figure 9 shows that the CCE can
be increased by a factor of $\sim$100 when the detector is operated at a field
of $\sim$ 400 V/cm with a charge collection time of 2 microseconds. This means
that the charge can be amplified by a factor of $\sim$100 at a bias voltage of
$\sim$1300 volts. Note that this extrapolation was obtained by extending the
fitted empirical trend as a function of applied electric field by more than
two orders of magnitude beyond the data. Therefore, the discussion above needs
to be verified using experimental data in the future.

Figure 9: Shown is the predicted CCE as a function of applied electric field.

##  5 Conclusion

We conclude that the observed time-dependent behavior of charge amplification
for the 59.54 keV line is mainly due to the impact ionization of charged
states by drifting holes across the detector. We attempted to set out possible
mechanisms for the impact ionization, which is mainly due to the creation of
cluster dipole states, and impact ionization of cluster dipole states. These
two processes dominate the production of additional charge that contributes to
the amplification of charge while drifting holes across the detector. Our
results suggest that a detector with an appropriate density of charged states,
created by controlled radiation from the appropriate density of neutral
states, can be developed for searching for MeV-scale dark matter with a
detection threshold as low as $<$10 meV, which corresponds to recoil energy
induced by sub-MeV dark matter particles. The challenge is to test such a
detector at a higher bias voltage such as 1300 volts without triggering
avalanche breakdown. The detector avalanche breakdown occurs at about 12 volts
in Run 67 and about 15 volts in Run 74. This could be increased by improving
the fabrication of the electrical contacts. But the overall event rate, from
the 241Am source plus background radiation, is also relevant to breakdown. Run
74 has a somewhat lower event rate compared to Run 67, and the avalanche
breakdown voltage was somewhat higher. We expect the avalanche breakdown
voltage can be much higher if events of only a few electron-hole pairs are
generated by a well-controlled optical source with the detector being well-
shielded from background radiation. More research and development is needed.

##  6 Acknowledgement

This work was supported in part by NSF OISE 1743790, NSF PHYS 1902577, DOE
grant DE-FG02-10ER46709, DE-SC0004768, and a research center supported by the
State of South Dakota.

## References

  * (1) F. Zwicky, “Die Rotverschiebung von extragalaktischen Nebeln,” Helv. Phys. Acta 6 (1933) 110-127.
  * (2) G. Hinshaw et al., (WMAP Collaboration), “Nine-Year Wilkinson Microwave Anisotropy Probe (WMAP) Observations: Cosmological Parameter Results,” ArXiv: 1212.5226, APJS, 208 (2013) 19.
  * (3) J. L. Feng, “Supersymmetry and Cosmology,” Annals of Physics 315, Issue 1 (2005) 2-15, arXiv:hep-ph/0405215.
  * (4) M. W. Goodman and E. Witten, “Detectability of certain dark-matter candidates,” Phys. Rev. D 31 (June, 1985) 3059-3063.
  * (5) Peccei, Roberto D.; Quinn, Helen R. (20 June 1977). "CP Conservation in the Presence of Pseudoparticles". Physical Review Letters. 38 (25): 1440–1443. Bibcode:1977PhRvL..38.1440P. doi:10.1103/PhysRevLett.38.1440
  * (6) James Bateman, Ian McHardy, Alexander Merle, Tim R. Morris, and Hendrik Ulbricht, "On the Existence of Low-Mass Dark Matter and its Direct Detection," Sci Rep. 2015: 5: 8058.
  * (7) Ungman, Gerard; Kamionkowski, Marc; Griest, Kim (1996). "Supersymmetric dark matter". Physics Reports. 267 (5–6): 195–373. arXiv:hep-ph/9506380.
  * (8) R. J. Gaitskell, “Direct Detection of Dark Matter,” Ann. Rev. Nucl. Part. Sci. 54 (2004) 315- 359\.
  * (9) The ADMX Collaboration; Asztalos, S.J.; Carosi, G.; Hagmann, C.; Kinion, D.; van Bibber, K.; Hotz, M.; Rosenberg, L.; Rybka, G.; Hoskins, J.; Hwang, J.; Sikivie, P.; Tanner, D. B.; Bradley, R.; Clarke, J. (28 January 2010). "A SQUID-based microwave cavity search for dark-matter axions". Physical Review Letters. 104 (4): 041301.
  * (10) Z. Ahmed et al. (CDMS Collaboration), _Search for Weakly Interacting Massive Particles with the First Five-Tower Data from the Cryogenic Dark Matter Search at the Soudan Underground Laboratory_ , _Phys. Rev. Lett._ 102 (2009) 011301.
  * (11) S.K. Liu et al. (CDEX Collaboration), _Limits on light WIMPs with a germanium detector at 177 eVee threshold at the China Jinping Underground Laboratory_ , _Phys. Rev. D_ 90 (2014) 032003.
  * (12) H. Jiang et al. (CDEX Collaboration), _Limits on Light Weakly Interacting Massive Particles from the First 102.8 kg $\times$ day Data of the CDEX-10 Experiment_, _Phys. Rev. Lett._ 120 (2018) 241301.
  * (13) L. T. Yang et al. (CDEX Collaboration), _Search for Light Weakly-Interacting-Massive-Particle Dark Matter by Annual Modulation Analysis with a Point-Contact Germanium Detector at the China Jinping Underground Laboratory_ , _Phys. Rev. Lett._ 123 (2019), 221301.
  * (14) Aalseth et al. (CoGeNT Collaboration), _Experimental constraints on a dark matter origin for the DAMA annual modulation effect_ , _Phys. Rev. Lett._ 101 (2008) 251301.
  * (15) G. Angloher et al. (CRESST Collaboration), _Results from 730 kg days of the CRESST-II Dark Matter Search_ , _Eur. Phys. J. C_ 72 (2012) 1971.
  * (16) E. Behnke et al. (COUPP Collaboration), _First dark matter search results from a 4-kg CF 3I bubble chamber operated in a deep underground site_, _Phys. Rev. D_ 86 (2012) 052001.
  * (17) J. Barreto et al., _Direct Search for Low Mass Dark Matter Particles with CCDs_ , _Phys. Lett. B_ 711 (2012) 264.
  * (18) R. Bernabei et al. (DAMA/LIBRA Collaboration), _New results from DAMA/LIBRA_ , _Eur. Phys. J. C_ 67 (2010) 39.
  * (19) M. Bossa et al. (DarkSide Collaboration), _DarkSide-50, a background free experiment for dark matter searches_ , _JINST_ 9 (2014) C01034.
  * (20) J.B.R. Battat et al. (DRIFT Collaboration), _First background-free limit from a directional dark matter experiment: results from a fully fiducialised DRIFT detector_ , _Phys. Dark Univ._ 9 (2014) 1.
  * (21) E. Armengaud et al. (EDELWEISS Collaboration), _Background studies for the EDELWEISS dark matter experiment_ , _Astropart. Phys._ 47 (2013) 1.
  * (22) S.C. Kim et al. (KIMS Collaboration), _The recent results from KIMS experiment_ , _J. Phys. Conf. Ser._ 384 (2012) 012020.
  * (23) D. S. Akerib et al. (LUX Collaboration), _Results from a search for dark matter in the complete LUX exposure_ , _Phys. Rev. Lett._ 118 (2017) 021303.
  * (24) M. Xiao et al. (PandaX Collaboration), _First dark matter search results from the PandaX-I experiment_ , _Sci. China Phys. Mech. Astron._ 57 (2014) 2024.
  * (25) C. Amole et al. (PICO Collaboration), _Dark Matter Search Results from the PICO-2L C 3F8 Bubble Chamber_, _Phys. Rev. Lett._ 114 (2015) 231302.
  * (26) R. Agnese et al. (SuperCDMS Collaboration), _Search for Low-Mass WIMPs with SuperCDMS_ , _Phys. Rev. Lett._ 112 (2014) 241302.
  * (27) E. Aprile et al. (XENON Collaboration), _Design and Performance of the XENON10 Dark Matter Experiment_ , _Astropart. Phys._ 34 (2011) 679.
  * (28) E. Aprile et al. (XENON Collaboration), _Exclusion of leptophilic dark matter models using XENON100 electronic recoil data_ , _Science_ 349 (2015) 851.
  * (29) E. Aprile et al. (XENON Collaboration), _First Dark Matter Search Results from the XENON1T Experiment_ , _Phys. Rev. Lett._ 119 (2017) 181301.
  * (30) K. Abe et al. (XMASS Collaboration), _Distillation of Liquid Xenon to Remove Krypton_ , _Astropart. Phys._ 31 (2009) 290.
  * (31) D.Y. Akimov et al. (ZEPLIN-III Collaboration), _Limits on inelastic dark matter from ZEPLIN-III_ , _Phys. Lett. B_ 692 (2010) 180.
  * (32) R. Essig, J. Mardon and T. Volansky, _Direct detection of sub-GeV dark matter_ , _Phys. Rev. D_ 85 (2012) 076007.
  * (33) R. Essig et al., _Direct detection of sub-GeV dark matter with semiconductor targets_ , _J. High Energ. Phys._ 2016 (2016) 46.
  * (34) C.M. Ho and R.J. Scherrer, _Limits on MeV dark matter from the effective number of neutrinos_ , _Phys. Rev. D_ 87 (2013) 023505.
  * (35) G. Steigman, _Equivalent neutrinos, light WIMPs, and the chimera of dark radiation_ , _Phys. Rev. D_ 87 (2013) 103517.
  * (36)
  * (37) D.-M. Mei et al., _Direct detection of MeV-scale dark matter utilizing germanium internal amplification for the charge created by the ionization of impurities_ , _Eur. Phys. J. C_ 78 (2018) 187.
  * (38)
  * (39) A. S. Starostin and A. G. Beda, Phys. Atom. Nucl. 63 (2000) 1297-1300, arXiv:hepex/0002063v1.
  * (40) W. Pickin, "Impact ionization in n-germanium, 4 - 9.5 K," Phys. Rev. B V20 (1979) 2451.
  * (41) N. Sclar and E. Burstein, "Impact ionization of impurity in germanium," J. of Phys. and Chem. of Solids, V2 (1957) 1-23.
  * (42) Mei, D-M and Panth, R and Kooi, K and Mei, H and Bhattarai, S and Raut, M and Acharya, P and Wang, G-J, "Evidence of cluster dipole states in germanium detectors operating at temperatures below 10 K", AIP Advances 12, 065113 (2022).
  * (43) N. Sclar and E. Burstein, "Impact Ionization of Impuriteis in Germanium," J. Phys. Chem. Solids. Pergamon Press 1957. Vol. 2 pp. 1-23.
  * (44) W. Pickin, "Impact ionization in n-germanium, 4.5 - 9.5 K," Phys. Rev. B, V 20 (1979) 2451.
  * (45) D. L. Smith, D. S. Pan, and T. C. McGill, "Impact ionization of excitons in Ge and Si," Phys. Rev. B, V 12 (1975) 4360.
  * (46) J. F. Palmier, "Impact ionization in n-type germanium," Phys. Rev. B, V 6 (1972) 4557.
  * (47) A. Phipps, B. Sadoulet, K. M. Sundqvist, "Observation of Impact Ionization of Shallow States in Sub-Kevin, High-Purity Germanium," J. Low Temp. Phys. (2016) 184: 336 - 343.
  * (48) Arran Thomas James Phipps, "Ionization Collection in Detectors of the Cryogenic Dark Matter Searches," PhD Thesis, University of California, Berkeley, Spring of 2016.
  * (49) G.-J. Wang et al., _Development of large size high-purity germanium crystal growth_ , _Journal of Crystal Growth_ 352 (2012) 27-30.
  * (50) G.-J. Wang et al., _High purity germanium crystal growth at the University of South Dakota_ , _J. Phys.: Conf. Ser._ 606 (2015) 012012.
  * (51) SuperCDMS Collaboration, R. Agnese et al., Phys. Rev. Lett. 112 (2014) 241302, arXiv:1402.7137.
  * (52) P. Luke, "Voltage-Assisted Calorimetric Ionization Detector," Journal of Applied Phsyics 64, 6858-6860 (1988).
  * (53) K. Sundqvist, "Carrier Transport and Related Effects in Detectors of Cryogenic Dark Matter Searches," PhD thesis, University of California, Berkeley, 2012.


# Theory of inertial spin dynamics in anisotropic ferromagnets

Mikhail Cherkasskii<EMAIL_ADDRESS>Faculty of Physics and Center
of Nanointegration (CENIDE), University of Duisburg-Essen, Duisburg 47057,
Germany  Igor Barsukov<EMAIL_ADDRESS>Department of Physics and Astronomy,
University of California, Riverside, California 92521, USA  Ritwik Mondal
Present address: Department of Physics, Indian Institute of Technology (ISM)
Dhanbad, IN-826004, Dhanbad, India  Department of Spintronics and
Nanoelectronics, Institute of Physics ASCR, v.v.i., Cukrovarnická 10, Prague
6, 162 53, Czech Republic  Department of Physics and Astronomy, Uppsala
University, Box 516, SE-75120 Uppsala, Sweden  Michael Farle  Faculty of
Physics and Center of Nanointegration (CENIDE), University of Duisburg-Essen,
Duisburg 47057, Germany  Anna Semisalova<EMAIL_ADDRESS>Faculty
of Physics and Center of Nanointegration (CENIDE), University of Duisburg-
Essen, Duisburg 47057, Germany

###### Abstract

Recent experimental observation of inertial spin dynamics calls upon holistic
reevaluation of the theoretical framework of magnetic resonance in
ferromagnets. Here, we derive the secular equation of an inertial spin system
in analogy to the ubiquitous Smit-Beljers formalism. We find that the
frequency of precessional ferromagnetic resonances is decreased as compared to
the noninertial case. We also find that the frequency of nutational resonances
is generally increased due to the presence of magnetic anisotropy and applied
magnetic field. We obtain exact solutions of the secular equation and
approximations that employ the terminology of non-inertial theory and thus
allow for convenient estimates of the inertial effects.

##  I Introduction

Inertial effects of spin dynamics have mostly been neglected in the analysis
of experimental spin resonance and spin transport results, since they mainly
manifest themselves at hardly accessible terahertz frequencies in
ferromagnets. With recent advances in spectroscopic techniques [1, 2], this
fundamental phenomenon has been observed in permalloy, cobalt and cobalt-iron-
boron ferromagnetic films [3, 4]. Inertial spin dynamics offers novel avenues
for ultrahigh-frequency spintronic applications using well-established
ferromagnetic materials and is thus becoming a prominent field of research.

The mathematical concept of inertia of magnetization was introduced in the
context of magnetoelastic coupling in ferromagnets by Suhl [5]. It was
followed by an extension of the breathing Fermi surface model which
demonstrated the emergence of a damping contribution linked to inertia [6, 7].
It was noticed that the Landau-Lifshitz-Gilbert equation, which describes
magnetization motion in analogy with a spinning top, required an inertial
tensor of a rigid body. A revision [8] of this analogy within a macroscopic
Lagrangian approach suggested that inertia originates from generalization of
gyromagnetic ratio –- the magnetic moment is non-collinear to the angular
momentum. The Landau-Lifshitz-Gilbert equation was extended by including the
second-order time derivative of magnetization, which now resembles Newton’s
equation of motion for a massive point particle.

This inertial Landau-Lifshitz-Gilbert (ILLG) equation reads

$\begin{split}{\partial_{t}}{\bf{M}}&=-\left|\gamma\right|{\mu_{0}}{\bf{M}}\times{{\bf{H}}_{{\rm{eff}}}}+\dfrac{\alpha}{{{M_{0}}}}{\bf{M}}\times{\partial_{t}}{\bf{M}}\\\
&+\dfrac{\eta}{{{M_{0}}}}{\bf{M}}\times{\partial_{tt}}{\bf{M}},\end{split}$
(1)

where $\gamma=g{\mu_{\rm{B}}}/\hbar$ is the gyromagnetic ratio,
${\mu_{\rm{B}}}$ is the Bohr magneton, $g$ is the g-factor, $\mu_{0}$ is the
permeability of free space, ${\bf{M}}$ is the magnetization with the magnitude
${M_{0}}$, ${{\bf{H}}_{{\rm{eff}}}}$ is the effective magnetic field, $\alpha$
is the Gilbert damping, and $\eta$ is the inertial parameter. The inertial
parameter has been discussed to be correlated to the Gilbert damping within
the breathing Fermi surface model [6, 7] and torque-torque correlation model
[9], whereas inertia has been considered independent of damping within the
classical Lagrangian approach [8]. The ILLG equation was derived within the
framework of mesoscopic nonequilibrium thermodynamics [10] and the Dirac-Kohn-
Sham theory [11]. In an atomistic method, exchange interaction, damping, and
inertia were calculated from first principles [12]. The microscopic origin of
inertia has been asserted in the relativistic spin–orbit coupling [13, 14, 9].

Inertia leads mainly to nutation – a terahertz-frequency motion of
magnetization superimposed on the regular gigahertz-frequency precession [15]
(Fig. 1). Nutational resonances have been discussed in ferromagnets [16] and
antiferromagnets [17, 18, 19, 20]. Moreover, traveling nutational spin waves
[21, 22, 23] have been proposed. Besides nutational motion, inertia has been
found to result in a frequency shift of the uniform magnetization precession
[15, 24] and spin waves [22] at gigahertz frequencies. Previous studies of
inertial spin dynamics treated various ferromagnetic systems including
nanoparticles and nanostructures [25, 26, 27]; however, a general approach
based on ILLG for a ferromagnet with an arbitrary magnetic anisotropy has not
yet been proposed.

Here, we develop a holistic theoretical framework for a ferromagnet with an
arbitrary magnetic anisotropy energy landscape in analogy to the Smit-Beljers
(SB) approach. For the last six decades [28], the Smit-Beljers formalism has
been an indispensable tool for routinely predicting and analyzing macrospin
ferromagnetic resonance and, with extensions, spin-wave resonances. Now, our
theoretical framework allows for deriving the frequencies of the nutational
and precessional resonances of anisotropic ferromagnets in the presence of
inertia. We formulate the secular equation of inertial spin dynamics, provide
its exact and approximate solutions, and discuss the interplay of inertia and
magnetic anisotropy.

Figure 1: The magnetization ${\bf{M}}$ precesses around the effective magnetic
field and ${{\bf{H}}_{{\rm{eff}}}}$ due to the precessional torque (blue). The
inertial torque (red) causes the magnetization to undergo a concurrent
nutational motion. The total magnetization motion is subject to the damping
torque (green).

##  II Precessional and nutational resonance frequencies

The Smit-Beljers approach was initially developed via linearizing the Landau-
Lifshitz equation in spherical coordinates and using small-angle approximation
around the equilibrium direction of magnetization [29, 30]. Employing the
periodic solution ansatz, the ferromagnetic resonance (FMR) frequency was
derived as

$\omega^{2}_{\rm{SB}}=\frac{|\gamma|^{2}\left(1+\alpha^{2}\right)}{{M_{0}^{2}{{\sin}^{2}}{\theta_{0}}}}\left({{\partial_{\theta\theta}}F{\partial_{\phi\phi}}F-{{\left({{\partial_{\theta\phi}}F}\right)}^{2}}}\right),$
(2)

where ${\theta_{0}}$ is the equilibrium polar angle of magnetization. The
equation gives a closed-form relation between the FMR frequency and magnetic
anisotropy energy $F$ and allows for efficient and convenient numerical
analysis of experimental data [31].

In our approach, we similarly write the ILLG equation in spherical coordinates
which now contains first and second derivatives of the angles. As detailed in
Appendix A, the relative magnitude of the terms in this resulting equation can
be analyzed by their prefactors expressed in the orders of the Gilbert damping
parameter $\alpha$. Since the latter is typically $10^{-3}-10^{-2}$, we omit
higher-order terms and arrive at

$\left\\{\begin{array}[]{l}{\partial_{tt}}\theta=\dfrac{{\left|\gamma\right|{\mu_{0}}{H_{\theta}}}}{\eta}-\dfrac{{\alpha\left|\gamma\right|{\mu_{0}}{H_{\varphi}}}}{\eta}+\dfrac{{{\partial_{t}}\varphi\sin\theta}}{\eta}\\\
+{\left({{\partial_{t}}\varphi}\right)^{2}}\sin\theta\cos\theta,\\\
{\partial_{tt}}\varphi\sin\theta=\dfrac{{\left|\gamma\right|{\mu_{0}}{H_{\varphi}}}}{\eta}+\dfrac{{\alpha\left|\gamma\right|{\mu_{0}}{H_{\theta}}}}{\eta}-\dfrac{{{\partial_{t}}\theta}}{\eta}\\\
-2{\partial_{t}}\varphi{\partial_{t}}\theta\cos\theta.\end{array}\right.$ (3)

Using the small-angle approximation, we develop the equations around the
equilibrium direction of magnetization which introduces second-order
derivatives of the energy $F$. The system of equations can be further
linearized employing a Jacobian matrix of the angles (Appendix A). Using the
periodic solution ansatz, we arrive at a fourth-order characteristic
polynomial constituting the secular equation of the inertial spin system:

$\begin{split}&\left[\frac{\omega^{2}}{\left|\gamma\right|^{2}}-\frac{{\left({1+{\alpha^{2}}}\right)}}{{M_{0}^{2}{{\sin}^{2}}{\theta_{0}}}}\left({{\partial_{\theta\theta}}F{\partial_{\phi\phi}}F-{{\left({{\partial_{\theta\phi}}F}\right)}^{2}}}\right)\right]\\\
&-{\eta^{2}\omega^{2}}\left[\frac{{\omega^{2}}}{\left|\gamma\right|^{2}}-\frac{1}{\eta\left|\gamma\right|M_{0}}\left({{\partial_{\theta\theta}}F+\frac{{{\partial_{\varphi\varphi}}F}}{{{{\sin}^{2}}{\theta_{0}}}}}\right)\right]\\\
&-i\omega\frac{{\alpha}}{\left|\gamma\right|M_{0}}\left({{\partial_{\theta\theta}}F+\frac{{{\partial_{\varphi\varphi}}F}}{{{{\sin}^{2}}{\theta_{0}}}}}\right)=0.\end{split}$
(4)

The first group of terms corresponds to Eq. 2. The second group of terms
introduces inertia of magnetization. The third group of terms corresponds to
the frequency-domain linewidth of the ferromagnetic resonance

$\Delta\omega_{\rm{SB}}=\frac{\left|\gamma\right|\alpha}{M_{0}}\left({{\partial_{\theta\theta}}F+\frac{{{\partial_{\varphi\varphi}}F}}{{{{\sin}^{2}}{\theta_{0}}}}}\right)$
(5)

as it does in the noninertial case [28, 32]. The presented approach has the
advantage to converge to the Smit-Beljers secular equation when the inertial
parameter vanishes and can be written as

$\left(\omega^{2}-\omega_{\mathrm{SB}}^{2}\right)-\eta^{2}\omega^{2}\left(\omega^{2}-\frac{1}{\eta\alpha}\Delta\omega_{\rm{SB}}\right)-i\omega\Delta\omega_{\rm{SB}}=0$
(6)

Equation 6 has two physical solutions: precessional resonance
$\omega_{\mathrm{p}}$ at lower frequency and nutational resonance
$\omega_{\mathrm{n}}$ at higher frequency. In Appendix B, we calculate the
explicit (exact but complex) solutions, shown in Fig. 2, as a benchmark for
the consecutive approximations.

First, similarly to the original Smit-Beljers formalism, we can omit the
imaginary term in the secular Equation 6 – an approximation that we mark with
’a’ – and derive an analytical form of the resonance frequencies:

$\omega_{\rm{p}}^{\left({\rm{a}}\right)}={\left({p-\sqrt{{{p}^{2}}-q}}\right)^{1/2}},$
(7)
$\omega_{\rm{n}}^{\left({\rm{a}}\right)}={\left({p+\sqrt{{{p}^{2}}-q}}\right)^{1/2}},$
(8)

where

$\begin{array}[]{l}p=\dfrac{1}{{{2\eta^{2}}}}+\dfrac{\Delta\omega_{\rm{SB}}}{2\alpha\eta},\\\
q=\dfrac{\omega_{\rm{SB}}^{2}}{\eta^{2}}.\end{array}$ (9)

The analytical form can be used conveniently to calculate the resonance
frequencies via $\omega_{\rm{SB}}$ and $\Delta\omega_{\rm{SB}}$, thus adding
just a few extra steps compared to the Smit-Beljers formalism. However, the
analytical form is still too complex and the effect of magnetic anisotropy on
precessional and nutational behavior is not immediately clear.

Figure 2: Frequency-field relation of ferromagnetic resonance and nutational
resonances. Explicit solutions of Eq. 6 for the precessional resonance
(orange) show a red-shift compared to the non-inertial Smit-Beljers case
(dashed blue). Explicit solutions of Eq. 6 for the nutational resonance (red)
show a blue-shift compared to the zeroth-order approximation $1/\eta$. (a) The
calculation parameters for a thin film with cubic magnetocrystalline
anisotropy are: $g=2.09$, ${\mu_{0}}{M_{0}}=2.1\>\rm{T}$, $\alpha=0.002$,
$\eta=75\>\rm{fs}\cdot\rm{rad}^{-1}$, $K_{\rm{cub1}}=4.9\times
10^{4}\>\rm{J}\cdot\rm{m}^{-3}$, Ref. [33, 3, 4]. (b)-(d) The following
parameters for a thin film with uniaxial magnetocrystalline anisotropy have
been used: $g=2.17$, ${\mu_{0}}{M_{0}}=1.8\>\rm{T}$, $\alpha=0.10$,
$\eta=75\>\rm{fs}\cdot\rm{rad}^{-1}$,
$K_{\rm{u1}}=4.1\times{10^{5}}\>\rm{J}\cdot\rm{m}^{-3}$, Ref. [34, 3, 4].

We thus implement another approximation – marked with ’b’ – by expanding the
analytical form into a Taylor series (Appendix B) and neglecting higher-order
terms in $\Delta\omega_{\rm{SB}}\eta\ll 1$:

$\begin{split}{\omega_{\rm{p}}^{(\rm{b})}}=\omega_{\rm{SB}}\sqrt{1-\eta\frac{\Delta\omega_{\rm{SB}}}{\alpha}},\end{split}$
(10)
${\omega}_{\rm{n}}^{(\rm{b})}=\frac{1}{\eta}+\frac{\Delta\omega_{\rm{SB}}}{2\alpha}.$
(11)

Here, we immediately see a systematic red-shift of the precessional frequency
as compared to the non-inertial case of $\omega_{\rm{SB}}$. As shown in Fig.
1, the inertial torque vector has a notable component that is antiparallel to
the precessional torque, thus effectively reducing the latter.

The nutational frequency ${\omega}_{\rm{n}}^{(\rm{b})}$, on the other hand,
shows a substantial frequency increase (a blueshift) as compared to an earlier
estimation $\omega_{\mathrm{n}}\sim 1/\eta$ for an isotropic ferromagnet [3,
4]. Another approximation, for instance employed in Ref. [15], accounts for a
frequency shift due to an applied magnetic field

${\bar{\omega}_{\rm{n}}}=\frac{{\sqrt{1+\eta\left|\gamma\right|{\mu_{0}}{H_{0}}}}}{\eta}=\frac{1}{\eta}+\frac{{\left|\gamma\right|{\mu_{0}}{H_{0}}}}{2}+...$
(12)

but neglects the effects of magnetic anisotropy. While the nutational
frequency obtained in our model converges to the estimate
$\bar{\omega}_{\rm{n}}$ in the case of vanishing anisotropy, it demonstrates
that magnetic anisotropy (both, shape and magnetocrystalline) shifts the
nutation resonance frequency as compared with ${\bar{\omega}_{\rm{n}}}$, and
must be accounted for according to the characteristic polynomial (Eq. 6) and
its solutions.

##  III Effect of magnetic anisotropy on inertial spin dynamics

We calculate the effect of magnetic anisotropy on precessional and nutational
resonances for four concrete examples of magnetic samples that have been and
may likely be used in an experiment probing inertial spin dynamics. We
consider single-crystal ferromagnetic thin films with cubic magnetocrystalline
anisotropy (iron) and uniaxial magnetocrystalline anisotropy (hexagonal-close-
packed cobalt) with magnetic parameters obtained from experimental data of
Refs. [3, 4]. The explicit (exact) solutions of the characteristic polynomial
of Eq. (6) are plotted in Fig. 2 for various configurations of applied
magnetic field with respect to the film surface and crystal symmetry axes. We
use free energy density and equilibrium angles defined in Appendices C and D.
As shown in Fig. 2, the effect of inertia is consistent in all calculated
scenarios. The precessional frequency experiences a redshift due to inertia as
compared to resonance frequency $\omega_{\rm{SB}}$ for noninertial Smit-
Beljers case. Both aligned precessional modes (above the saturation field) and
nonaligned precessional modes (below the saturation field in hard-axes
configurations) [31] experience a redshift which increases with increasing
precessional frequency. For the aligned modes, the redshift thus becomes
stronger with increasing magnetic field. For nonaligned modes, on the other
hand, decreasing magnetic field can result in increasing redshift.

The nutational frequency experiences a blue shift due to magnetic anisotropy
and magnetic field. The blueshift typically increases with increasing magnetic
field. However, below the saturation field in hard-axes configurations, the
blue-shift of the nutational frequency can become stronger with decreasing
magnetic field (see the red line in Fig. 2(a-c)). The exact behavior is
dominated by the $\Delta\omega_{\mathrm{SB}}$ term in Eqs. (10) and (11). Both
shifts can be substantial (reaching up to 12% in the magnetic field of 10T)
for all calculated scenarios.

While we use the explicit solutions of the secular equation in Fig. 2 to
visualize the effects of magnetic anisotropy in inertial spin systems, the
observed frequency shifts are in qualitative agreement with the
approximations. To assess the quantitative validity of our approximations ’a’
and ’b’, we compare them with the benchmark of the explicit solutions. We find
that the analytical form ’a’ (Eqs. (7-8)) accrues less than 0.5% error for
aligned modes when compared to the exact solution of Eq. 4; however, it should
be stressed that the characteristic polynomial Eq. 4 itself has been derived
for inertial parameters $\eta\ll 1/\Delta\omega_{\rm{SB}}$. The next-step
approximation by the Taylor series ’b’ (Eqs. (10-11)) introduces an error less
than 6% for aligned modes with $\eta<100\>\rm{fs}\cdot\rm{rad}^{-1}$, while at
higher values of the inertial parameter, the Taylor series causes a
substantial error of the frequencies. While the approximation ’b’ in Eqs.
(10-11) should thus be treated with caution, for comparison, we calculate the
explicit forms of precessional and nutational frequencies for aligned modes in
the configurations displayed in Fig. 2:

(a) The out-of-plane cubic case:

$\begin{split}&{\omega_{\rm{p}}^{(\rm{b})}}^{2}=\left|\gamma\right|^{2}\left(1+\alpha^{2}\right)\left(-\mu_{0}M_{0}+\mu_{0}H_{0}+\dfrac{2K_{\rm{cub1}}}{M}\right)^{2}\\\
&\times\left[1-\eta\left|\gamma\right|\left(-2\mu_{0}M_{0}+2\mu_{0}H_{0}+\dfrac{4K_{\rm{cub1}}}{M}\right)\right],\end{split}$
(13)

${\omega}_{\rm{n}}^{(\rm{b})}=\frac{1}{\eta}+\left|\gamma\right|\left({-{\mu_{0}}{M_{0}}+{\mu_{0}}{H_{0}}+\frac{{2{K_{{\rm{cub1}}}}}}{{{M_{0}}}}}\right)$
(14)

(b) The out-of-plane uniaxial case:

$\begin{split}{\omega_{\rm{p}}^{(\rm{b})}}^{2}=\left(1+\alpha^{2}\right)|\gamma|^{2}\left(-\mu_{0}M_{0}+\mu_{0}H_{0}+\dfrac{2K_{\rm{u1}}}{M_{0}}\right)^{2}\times\\\
\left(1-\eta|\gamma|\left(-2\mu_{0}M_{0}+2\mu_{0}H_{0}+\dfrac{4K_{\text{u1}}}{M_{0}}\right)\right),\end{split}$
(15)

$\omega_{\rm{n}}^{\left({\rm{b}}\right)}=\frac{1}{\eta}+\left|\gamma\right|\left({-{\mu_{0}}{M_{0}}+{\mu_{0}}{H_{0}}+\frac{{2{K_{{\rm{u1}}}}}}{{{M_{0}}}}}\right).$
(16)

(c) The in-plane perpendicular uniaxial case:

$\begin{split}\begin{array}[]{l}{\omega_{\rm{p}}^{(\rm{b})}}^{2}=\left(1+\alpha^{2}\right)|\gamma|^{2}\left(\mu_{0}M_{0}+\mu_{0}H_{0}\right)\left(\mu_{0}H_{0}-\dfrac{2K_{\rm{u1}}}{M_{0}}\right)\times\\\
\left(1-\eta|\gamma|\left(\mu_{0}M_{0}+2\mu_{0}H_{0}-\dfrac{2K_{\rm{u1}}}{M_{0}}\right)\right),\end{array}\end{split}$
(17)

$\omega_{\rm{n}}^{\left({\rm{b}}\right)}=\frac{1}{\eta}+\left|\gamma\right|\left({\frac{{{\mu_{0}}{M_{0}}}}{2}+{\mu_{0}}{H_{0}}-\frac{{{K_{{\rm{u1}}}}}}{{{M_{0}}}}}\right).$
(18)

(d) The in-plane parallel uniaxial case:

$\begin{split}\begin{array}[]{l}{\omega_{\rm{p}}^{(\rm{b})}}^{2}=\left(1+\alpha^{2}\right)|\gamma|^{2}\left(\mu_{0}H_{0}+\dfrac{2K_{\rm{u1}}}{M_{0}}\right)\times\\\
\left(\mu_{0}M_{0}+\mu_{0}H_{0}+\dfrac{2K_{\rm{u1}}}{M_{0}}\right)\times\\\
\left(1-\eta|\gamma|\left(\mu_{0}M_{0}+2\mu_{0}H_{0}+\dfrac{4K_{\rm{u1}}}{M_{0}}\right)\right),\end{array}\end{split}$
(19)

$\omega_{\rm{n}}^{\left({\rm{b}}\right)}=\frac{1}{\eta}+\left|\gamma\right|\left({\frac{{{\mu_{0}}{M_{0}}}}{2}+{\mu_{0}}{H_{0}}+\frac{{2{K_{{\rm{u1}}}}}}{{{M_{0}}}}}\right).$
(20)

It is a common practice to analyze experimentally determined dependences of
ferromagnetic resonance frequency on applied magnetic field for evaluating
magnetic parameters such as magnetic anisotropy and g-factor [35, 36, 31, 37,
38, 39, 40]. Our work demonstrates that such evaluation needs to be adjusted
by taking into account the inertial redshift. In particular, measurements at
higher fields/frequencies have been considered to result in more accurate
determination of magnetic parameters [41, 42, 43, 44, 45]. Our model, however,
shows that especially at high magnetic fields, inertial redshift is strong and
needs to be taken into account.

It should be noted that in the framework of the extended breathing Fermi
surface model [6, 7], the inertial term with negative sign was derived. Such
negative inertial term would formally result in a blueshift of the
precessional frequencies. However, since the origin of inertia is still under
discussion, we consider here only the effects of the positive inertial term
suggested in Ref. [8].

##  IV Summary

In summary, we derived the secular equation for an inertial spin system with
an arbitrary magnetic anisotropy energy in analogy with the Smit-Beljers
approach. We find that ferromagnetic resonance experiences a substantial
redshift due to the inertia, while nutational resonance experiences a
blueshift due to magnetic anisotropy and field. For an accurate evaluation of
magnetic parameters from magnetic resonance measurements, inertia needs to be
taken into account. Our model (Eq. 6) allows for convenient calculation of
precessional and nutational resonances of an inertial spin system using
parameters ($\omega_{\mathrm{SB}}$ and $\Delta\omega_{\mathrm{SB}}$) obtained
from noninertial models.

## ACKNOWLEDGMENTS

M.Ch., M.F. and A.S. acknowledge funding by Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) – Project No. 392402498 (SE 2853/1-1) and
Project No. 405553726 CRC/TRR 270. R.M. acknowledges the funding from the
Swedish Research Council via VR 2019-06313. I.B. acknowledges funding by the
National Science Foundation through Grant No. ECCS-1810541.

##  Appendix A SMIT-BELJERS APPROACH WITH INERTIA

First, we transform the ILLG equation into a spherical coordinate system and
find equilibrium angles of magnetization. Second, we linearize the system of
equations describing magnetization dynamics at the equilibrium point. Finally,
we derive the eigenfrequencies corresponding to the resonances.

In a spherical coordinate system one writes
${\bf{M}}={M_{0}}{{\bf{e}}_{\rm{r}}},$
${{\bf{H}}_{{\rm{eff}}}}={H_{\rm{r}}}{{\bf{e}}_{\rm{r}}}+{H_{\rm{\theta}}}{{\bf{e}}_{\rm{\theta}}}+{H_{\rm{\varphi}}}{{\bf{e}}_{\rm{\varphi}}},$
where the magnitude of the magnetization vector persists over time. Using

$\begin{array}[]{l}{\partial_{t}}{\bf{M}}={M_{0}}\left({{\partial_{t}}\theta{{\bf{e}}_{\rm{\theta}}}+\sin\theta{\partial_{t}}\varphi{{\bf{e}}_{\rm{\varphi}}}}\right),\\\
{\partial_{tt}}{\bf{M}}={M_{0}}\left\\{{\left[{-{{\left({{\partial_{t}}\theta}\right)}^{2}}-{{\left({{\partial_{t}}\varphi}\right)}^{2}}{{\sin}^{2}}\theta}\right]{{\bf{e}}_{\rm{r}}}}\right.\\\\[4.0pt]
+\left[{{\partial_{tt}}\theta-{{\left({{\partial_{t}}\varphi}\right)}^{2}}\sin\theta\cos\theta}\right]{{\bf{e}}_{\rm{\theta}}}\\\\[5.0pt]
\left.{+\left[{{\partial_{tt}}\varphi\sin\theta+2{\partial_{t}}\varphi{\partial_{t}}\theta\cos\theta}\right]{{\bf{e}}_{\rm{\varphi}}}}\right\\},\end{array}$
(21)

one transforms the ILLG equation into

$\left\\{\begin{array}[]{l}{\partial_{tt}}\theta=\dfrac{{\left|\gamma\right|{\mu_{0}}{H_{\rm{\theta}}}}}{\eta}-\dfrac{{\alpha{\color[rgb]{0.71875,0.05859375,0.2265625}{\partial_{t}}\theta}}}{\eta}+\dfrac{{{\partial_{t}}\varphi\sin\theta}}{\eta}\\\
+{\left({{\partial_{t}}\varphi}\right)^{2}}\sin\theta\cos\theta,\\\\[8.0pt]
{\partial_{tt}}\varphi\sin\theta=\dfrac{{\left|\gamma\right|{\mu_{0}}{H_{\rm{\varphi}}}}}{\eta}-\dfrac{{\alpha{\color[rgb]{0.71875,0.05859375,0.2265625}{\partial_{t}}\varphi\sin\theta}}}{\eta}-\dfrac{{{\partial_{t}}\theta}}{\eta}\\\
-2{\partial_{t}}\varphi{\partial_{t}}\theta\cos\theta.\end{array}\right.$ (22)

Here, we introduce the first approximation, i.e. we replace the ”red” terms in
the system (22) with the ”blue” ones in the system (25) as follows. Based on
the ILLG equation, we write

$\left\\{\begin{array}[]{l}\dfrac{{\alpha{\partial_{t}}\varphi\sin\theta}}{\eta}=-\dfrac{{\alpha\left|\gamma\right|{\mu_{0}}{H_{\theta}}}}{\eta}+\dfrac{{{\alpha^{2}}{\partial_{t}}\theta}}{\eta}\\\
{\rm{}}-\alpha{\left({{\partial_{t}}\varphi}\right)^{2}}\sin\theta\cos\theta+\alpha{\partial_{tt}}\theta,\\\\[8.0pt]
\dfrac{{\alpha{\partial_{t}}\theta}}{\eta}=\dfrac{{\alpha\left|\gamma\right|{\mu_{0}}{H_{\varphi}}}}{\eta}-\dfrac{{{\alpha^{2}}{\partial_{t}}\varphi\sin\theta}}{\eta}\\\
{\rm{}}-2\alpha{\partial_{t}}\varphi{\partial_{t}}\theta\cos\theta-\alpha{\partial_{tt}}\varphi\sin\theta\par\end{array}\right.$
(23)

and substitute these equations in (22) instead of the corresponding ”color”
terms. We obtain

$\left\\{\begin{array}[]{l}{\partial_{tt}}\theta=\dfrac{{\left|\gamma\right|{\mu_{0}}{H_{\theta}}}}{\eta}-\dfrac{{\alpha\left|\gamma\right|{\mu_{0}}{H_{\varphi}}}}{\eta}+\dfrac{{{\partial_{t}}\varphi\sin\theta}}{\eta}\left({1+{\alpha^{2}}}\right)\\\
+{\left({{\partial_{t}}\varphi}\right)^{2}}\sin\theta\cos\theta{\rm{}}+2\alpha{\partial_{t}}\varphi{\partial_{t}}\theta\cos\theta+\alpha{\partial_{tt}}\varphi\sin\theta,\\\\[8.0pt]
{\partial_{tt}}\varphi\sin\theta=\dfrac{{\left|\gamma\right|{\mu_{0}}{H_{\varphi}}}}{\eta}+\dfrac{{\alpha\left|\gamma\right|{\mu_{0}}{H_{\theta}}}}{\eta}-\dfrac{{{\partial_{t}}\theta}}{\eta}\left({1+{\alpha^{2}}}\right)\\\
-2{\partial_{t}}\varphi{\partial_{t}}\theta\cos\theta{\rm{}}+\alpha{\left({{\partial_{t}}\varphi}\right)^{2}}\sin\theta\cos\theta-\alpha{\partial_{tt}}\theta.\end{array}\right.$
(24)

The last two terms in both equations are negligible, since they are multiplied
by $\alpha<1$, while $\Delta\omega_{\rm{SB}}\eta\ll 1$. The ${\alpha^{2}}$
terms are much less than 1, hence the ILLG equation is converted to

$\left\\{\begin{array}[]{l}{\partial_{tt}}\theta=\dfrac{{\left|\gamma\right|{\mu_{0}}{H_{\rm{\theta}}}}}{\eta}-\dfrac{{\alpha{\color[rgb]{0.36328125,0.5,0.70703125}\left|\gamma\right|{\mu_{0}}{H_{\rm{\varphi}}}}}}{\eta}+\dfrac{{{\partial_{t}}\varphi\sin\theta}}{\eta}\\\
+{\left({{\partial_{t}}\varphi}\right)^{2}}\sin\theta\cos\theta,\\\\[8.0pt]
{\partial_{tt}}\varphi\sin\theta=\dfrac{{\left|\gamma\right|{\mu_{0}}{H_{\rm{\varphi}}}}}{\eta}+\dfrac{{\alpha{\color[rgb]{0.36328125,0.5,0.70703125}\left|\gamma\right|{\mu_{0}}{H_{\rm{\theta}}}}}}{\eta}-\dfrac{{{\partial_{t}}\theta}}{\eta}\\\
-2{\partial_{t}}\varphi{\partial_{t}}\theta\cos\theta.\end{array}\right.$ (25)

This transformation is commonly adopted and was performed in Ref. [32]. The
advantage of the first approximation is that the final result, which is to be
shown below, converges to the SB equation for $\eta=0$. Note that the
effective magnetic field
${{\bf{H}}_{{\rm{eff}}}}=-\mu_{0}^{-1}{\partial_{\bf{M}}}F$ in the spherical
coordinate system is given by

${H_{\rm{\theta}}}=-\frac{1}{{{\mu_{0}}{M_{0}}}}{\partial_{\theta}}F,\,{H_{\rm{\varphi}}}=-\frac{1}{{{\mu_{0}}{M_{0}}\sin\theta}}{\partial_{\varphi}}F.$
(26)

In order to find the eigenfrequencies from the nonlinear system of equations
(25), it is necessary to linearize it and to determine the equilibrium
orientation of magnetization. The equilibrium given by the angles
${\theta_{0}}$ and ${\varphi_{0}}$ is found from the extremum conditions

${\partial_{\theta}}F=0,\,{\partial_{\varphi}}F=0$ (27)

limited by the conditions for the minimum, namely, the determinant of a
Hessian matrix has to be positive

${\partial_{\theta\theta}}F{\partial_{\varphi\varphi}}F-{\partial_{\theta\varphi}}F{\partial_{\varphi\theta}}F>0$
(28)

and one of the second derivative has to be positive as well

${\partial_{\theta\theta}}F>0.$ (29)

In the excited state, magnetization is deflected from the equilibrium
orientation by the effective magnetic field changes over time. Here, we
introduce the second approximation, which corresponds to the standard SB
approach, that is the deflection from equilibrium is small

$\Delta\theta\left(t\right)=\theta\left(t\right)-{\theta_{0}},\,\Delta\varphi\left(t\right)=\varphi\left(t\right)-{\varphi_{0}}$
(30)

and it is sufficient to limit the expansion of free energy to the linear terms

$\begin{array}[]{l}{\partial_{\theta}}F={\partial_{\theta\theta}}F\Delta\theta+{\partial_{\theta\varphi}}F\Delta\varphi,\\\\[4.0pt]
{\partial_{\varphi}}F={\partial_{\theta\varphi}}F\Delta\theta+{\partial_{\varphi\varphi}}F\Delta\varphi,\end{array}$
(31)

where the second derivatives are evaluated at the equilibrium. Using the
small-angle approximation, one obtains from expressions (25)-(31)

$\begin{array}[]{l}{\partial_{tt}}\Delta\theta=\left({-\dfrac{{\left|\gamma\right|}}{{\eta{M_{0}}}}{\partial_{\theta\theta}}F+\dfrac{{\alpha\left|\gamma\right|}}{{\eta{M_{0}}\sin{\theta_{0}}}}{\partial_{\theta\varphi}}F}\right)\Delta\theta\\\\[8.0pt]
+\left({-\dfrac{{\left|\gamma\right|}}{{\eta{M_{0}}}}{\partial_{\theta\varphi}}F+\dfrac{{\alpha\left|\gamma\right|}}{{\eta{M_{0}}\sin{\theta_{0}}}}{\partial_{\varphi\varphi}}F}\right)\Delta\varphi\\\\[8.0pt]
+\dfrac{{{\partial_{t}}\Delta\varphi\sin{\theta_{0}}}}{\eta}+{\left({{\partial_{t}}\Delta\varphi}\right)^{2}}\sin{\theta_{0}}\cos{\theta_{0}},\\\\[8.0pt]
{\partial_{tt}}\Delta\varphi\sin{\theta_{0}}=\left({-\dfrac{{\left|\gamma\right|}}{{\eta{M_{0}}\sin{\theta_{0}}}}{\partial_{\theta\varphi}}F-\dfrac{{\alpha\left|\gamma\right|}}{{\eta{M_{0}}}}{\partial_{\theta\theta}}F}\right)\Delta\theta\\\\[8.0pt]
+\left({-\dfrac{{\left|\gamma\right|}}{{\eta{M_{0}}\sin{\theta_{0}}}}{\partial_{\varphi\varphi}}F-\dfrac{{\alpha\left|\gamma\right|}}{{\eta{M_{0}}}}{\partial_{\theta\varphi}}F}\right)\Delta\varphi\\\\[8.0pt]
-\dfrac{{{\partial_{t}}\Delta\theta}}{\eta}-2{\partial_{t}}\Delta\varphi{\partial_{t}}\Delta\theta\cos{\theta_{0}}.\end{array}$
(32)

In order to linearize this system of equations, the following notations are
introduced

$\begin{array}[]{l}{a_{41}}=-\dfrac{{\left|\gamma\right|{\partial_{\theta\varphi}}F}}{{\eta{M_{0}}{{\sin}^{2}}{\theta_{0}}}}-\dfrac{{\alpha\left|\gamma\right|{\partial_{\theta\theta}}F}}{{\eta{M_{0}}\sin{\theta_{0}}}},\\\\[8.0pt]
{a_{42}}=-\dfrac{1}{{\eta\sin{\theta_{0}}}},\\\\[8.0pt]
{a_{43}}=-\dfrac{{\left|\gamma\right|{\partial_{\varphi\varphi}}F}}{{\eta{M_{0}}{{\sin}^{2}}{\theta_{0}}}}-\dfrac{{\alpha\left|\gamma\right|{\partial_{\theta\varphi}}F}}{{\eta{M_{0}}\sin{\theta_{0}}}},\\\\[8.0pt]
{\nu_{4}}=-2\cot{\theta_{0}},\\\\[8.0pt]
{a_{21}}=\dfrac{\alpha{\left|\gamma\right|{\partial_{\theta\varphi}}F}}{{\eta{{M}_{0}}\sin{\theta_{0}}}}-\dfrac{{\left|\gamma\right|{\partial_{\theta\theta}}F}}{{\eta{{M}_{0}}}},\\\\[8.0pt]
{a_{23}}=\dfrac{{\alpha\left|\gamma\right|{\partial_{\varphi\varphi}}F}}{{\eta{M_{0}}\sin{\theta_{0}}}}-\dfrac{{\left|\gamma\right|{\partial_{\theta\varphi}}F}}{{\eta{M_{0}}}},\\\\[8.0pt]
{a_{24}}=\dfrac{{\sin{\theta_{0}}}}{\eta},\\\\[8.0pt]
{\nu_{2}}=\sin{\theta_{0}}\cos{\theta_{0}},\\\\[8.0pt]
{x_{1}}=\Delta\theta,\\\\[8.0pt]
{x_{2}}={\partial_{t}}\Delta\theta,\\\\[8.0pt]
{x_{3}}=\Delta\varphi,\\\\[8.0pt]
{x_{4}}={\partial_{t}}\Delta\varphi.\end{array}$ (33)

Employing (33), we rewrite the system (32) as

$\begin{array}[]{l}{\partial_{t}}{x_{1}}={x_{2}},\\\
{\partial_{t}}{x_{2}}={a_{21}}{x_{1}}+{a_{23}}{x_{3}}+{a_{24}}{x_{4}}+{\nu_{2}}x_{4}^{2},\\\
{\partial_{t}}{x_{3}}={x_{4}},\\\
{\partial_{t}}{x_{4}}={a_{41}}{x_{1}}+{a_{42}}{x_{2}}+{a_{43}}{x_{3}}+{\nu_{4}}{x_{2}}{x_{4}}.\end{array}$
(34)

At the fixed point
${{\bf{x}}^{*}}=\left({x_{1}^{*},x_{2}^{*},x_{3}^{*},x_{4}^{*}}\right)$ the
dynamics of the nonlinear system (34) are qualitatively similar to the
dynamics of a linear system (35) associated with the Jacobian matrix
$J\left({{{\bf{x}}^{*}}}\right)$ [46], i.e.,

$\begin{split}&\left({\begin{array}[]{*{20}{l}}{\partial_{t}}{x_{1}}\\\
{\partial_{t}}{x_{2}}\\\ {\partial_{t}}{x_{3}}\\\
{{\partial_{t}}{x_{4}}}\end{array}}\right)=\\\
&\left({\begin{array}[]{*{20}{c}}{{\partial_{{x_{1}}}}{f_{1}}\left({{{\bf{x}}^{*}}}\right)}&\ldots&{{\partial_{{x_{4}}}}{f_{1}}\left({{{\bf{x}}^{*}}}\right)}\\\
\vdots&\ddots&\vdots\\\
{{\partial_{{x_{1}}}}{f_{4}}\left({{{\bf{x}}^{*}}}\right)}&\ldots&{{\partial_{{x_{4}}}}{f_{4}}\left({{{\bf{x}}^{*}}}\right)}\end{array}}\right)\left({\begin{array}[]{*{20}{l}}{x_{1}}\\\
{x_{2}}\\\ {x_{3}}\\\ {x_{4}}\end{array}}\right),\end{split}$ (35)

where the right-hand sides of Eqs. (34) are denoted as ${f_{i}}$. The fixed
point is determined by equating the derivatives of the nonlinear system (34)
to zero, which gives the following equations

$\begin{array}[]{l}{x_{2}}=0,\\\ {a_{21}}{x_{1}}+{a_{23}}{x_{3}}=0,\\\
{x_{4}}=0,\\\ {a_{41}}{x_{1}}+{a_{43}}{x_{3}}=0,\end{array}$ (36)

with the solution $x_{1}^{*}=x_{2}^{*}=x_{3}^{*}=x_{4}^{*}=0$. The Jacobian
matrix of Eqs. (34) is

$J=\left[{\begin{array}[]{*{20}{c}}0&1&0&0\\\
{{a_{21}}}&0&{{a_{23}}}&{{a_{24}}+2{\nu_{2}}{x_{4}}}\\\ 0&0&0&1\\\
{{a_{41}}}&{{a_{42}}+{\nu_{4}}{x_{4}}}&{{a_{43}}}&{{\nu_{4}}{x_{2}}}\end{array}}\right]$
(37)

and at the point $x_{1}^{*}=x_{2}^{*}=x_{3}^{*}=x_{4}^{*}=0$ it provides the
linear system of equations

$\begin{array}[]{l}{\partial_{t}}{x_{1}}={x_{2}},\\\
{\partial_{t}}{x_{2}}={a_{21}}{x_{1}}+{a_{23}}{x_{3}}+{a_{24}}{x_{4}},\\\
{\partial_{t}}{x_{3}}={x_{4}},\\\
{\partial_{t}}{x_{4}}={a_{41}}{x_{1}}+{a_{42}}{x_{2}}+{a_{43}}{x_{3}}.\end{array}$
(38)

This third approximation goes beyond the SB approach and it is the
linearization of the system (34). The eigenvalues of these equations give
resonance frequencies, which are calculated from the characteristic polynomial

$\begin{split}&{\omega^{4}}+\left({{a_{21}}+{a_{24}}{a_{42}}+{a_{43}}}\right){\omega^{2}}\\\
&-i\left({{a_{24}}{a_{41}}+{a_{23}}{a_{42}}}\right)\omega-{a_{23}}{a_{41}}+{a_{21}}{a_{43}}=0.\end{split}$
(39)

Restoring the original variable notations, one finds the equation describing
eigenfrequencies of a ferromagnet with inertia

$\begin{split}&\left[\frac{\omega^{2}}{\left|\gamma\right|^{2}}-\frac{{\left({1+{\alpha^{2}}}\right)}}{{M_{0}^{2}{{\sin}^{2}}{\theta_{0}}}}\left({{\partial_{\theta\theta}}F{\partial_{\phi\phi}}F-{{\left({{\partial_{\theta\phi}}F}\right)}^{2}}}\right)\right]\\\
&-{\eta^{2}\omega^{2}}\left[\frac{{\omega^{2}}}{\left|\gamma\right|^{2}}-\frac{1}{\eta\left|\gamma\right|M_{0}}\left({{\partial_{\theta\theta}}F+\frac{{{\partial_{\varphi\varphi}}F}}{{{{\sin}^{2}}{\theta_{0}}}}}\right)\right]\\\
&-i\omega\frac{{\alpha}}{\left|\gamma\right|M_{0}}\left({{\partial_{\theta\theta}}F+\frac{{{\partial_{\varphi\varphi}}F}}{{{{\sin}^{2}}{\theta_{0}}}}}\right)=0.\end{split}$
(40)

Note that this equation can be converted to SB formula (2) if the inertial
parameter vanishes.

##  Appendix B EXACT AND APPROXIMATE EXPRESSIONS OF RESONANCE FREQUENCIES

The quartic equation (40) results in two pairs of roots. The first pair is
precessional frequency modified by inertia, one root of the pair is positive,
the second is negative. The same applies to the other pair corresponding to
the nutational frequency. Here we consider only positive roots. Let us use the
Ferrari’s solution for this quartic equation to write exact expressions of
resonance frequencies, and introduce the notations:

$\begin{split}&{A_{\rm{r}}}=M_{0}^{2}{\eta^{2}},\\\\[4.0pt]
&{C_{\rm{r}}}=-M_{0}^{2}-{M_{0}}\eta\left|\gamma\right|\left({{\partial_{\theta\theta}}F+\dfrac{{{\partial_{\varphi\varphi}}F}}{{{{\sin}^{2}}{\theta_{0}}}}}\right),\\\\[4.0pt]
&{D_{\rm{r}}}=i{M_{0}}\alpha\left|\gamma\right|\left({{\partial_{\theta\theta}}F+\dfrac{{{\partial_{\varphi\varphi}}F}}{{{{\sin}^{2}}{\theta_{0}}}}}\right)\\\\[4.0pt]
&{E_{\rm{r}}}=\dfrac{{{{\left|\gamma\right|}^{2}}\left({1+{\alpha^{2}}}\right)}}{{{{\sin}^{2}}{\theta_{0}}}}\left({{\partial_{\theta\theta}}F{\partial_{\phi\phi}}F-{{\left({{\partial_{\theta\phi}}F}\right)}^{2}}}\right),\\\\[4.0pt]
&{a_{\rm{r}}}=\dfrac{{{C_{\rm{r}}}}}{{{A_{\rm{r}}}}},\;{b_{\rm{r}}}=\dfrac{{{D_{\rm{r}}}}}{{{A_{\rm{r}}}}},\;{c_{\rm{r}}}=\dfrac{{{E_{\rm{r}}}}}{{{A_{\rm{r}}}}}.\end{split}$
(41)

In Ferrari’s method, one determines a root of the nested depressed cubic
equation. In our case, the root is written

${y_{\rm{r}}}=-\frac{{5{a_{\rm{r}}}}}{6}+{U_{\rm{r}}}+{V_{\rm{r}}},$ (42)

where

$\begin{array}[]{l}{U_{\rm{r}}}=\sqrt[3]{{-\sqrt{\dfrac{{P_{\rm{r}}^{3}}}{{27}}+\dfrac{{Q_{\rm{r}}^{2}}}{4}}-\dfrac{{{Q_{\rm{r}}}}}{2}}},\\\
{V_{\rm{r}}}=-\dfrac{{{P_{\rm{r}}}}}{{3{U_{\rm{r}}}}},\\\
{P_{\rm{r}}}=-\dfrac{{a_{\rm{r}}^{2}}}{{12}}-{c_{\rm{r}}},\\\
{Q_{\rm{r}}}=\dfrac{1}{3}{a_{\rm{r}}}{c_{\rm{r}}}-\dfrac{{a_{\rm{r}}^{3}}}{{108}}-\dfrac{{b_{\rm{r}}^{2}}}{8}.\end{array}$
(43)

Thus, the exact precessional angular frequency modified by inertia is given by

$\begin{split}&\omega_{\rm{p}}=\frac{{\sqrt{{a_{\rm{r}}}+2{y_{\rm{r}}}}}}{2}\\\
&-\frac{1}{2}\sqrt{-3{a_{\rm{r}}}-2{y_{\rm{r}}}-\frac{{2{b_{\rm{r}}}}}{{\sqrt{{a_{\rm{r}}}+2{y_{\rm{r}}}}}}},\end{split}$
(44)

The exact nutational angular frequency can be written as

$\begin{split}&\omega_{\rm{n}}=\frac{{\sqrt{{a_{\rm{r}}}+2{y_{\rm{r}}}}}}{2}\\\
&+\frac{1}{2}\sqrt{-3{a_{\rm{r}}}-2{y_{\rm{r}}}-\frac{{2{b_{\rm{r}}}}}{{\sqrt{{a_{\rm{r}}}+2{y_{\rm{r}}}}}}}.\end{split}$
(45)

Next, we write a few approximations allowing one to elucidate the physics
behind Eqs. (44)-(44). The approximation ”a” of resonance frequencies is
derived by taking into account the real part of the quartic equation (40),
which transforms this equation into a bi-quadratic one. Thus, the
approximation ”a” of precessional frequency reads

$\omega_{\rm{p}}^{\left({\rm{a}}\right)}={\left({p-\sqrt{{{p}^{2}}-q}}\right)^{1/2}}.$
(46)

The approximate nutational frequency is given by

$\omega_{\rm{n}}^{\left({\rm{a}}\right)}={\left({p+\sqrt{{{p}^{2}}-q}}\right)^{1/2}},$
(47)

where

$\begin{array}[]{l}p=\dfrac{1}{{{2\eta^{2}}}}+\dfrac{\Delta\omega_{\rm{SB}}}{2\alpha\eta},\\\
q=\dfrac{\omega_{\rm{SB}}^{2}}{\eta^{2}},\end{array}$ (48)

This approximation introduces an additional error, which does not exceed 0.5%
for the aligned modes for the parameters employed in the main part of the
paper. We thus find that the solution ’a’ can be considered sufficiently
accurate in the context of this work.

Expressions (46) and (47) can be further simplified by employing Taylor series
expansion and assumption that $\Delta\omega_{\rm{SB}}\eta\ll 1$. The
approximation ’b’ of precessional frequency is

$\begin{split}{\omega_{\rm{p}}^{(\rm{b})}}=\omega_{\rm{SB}}\sqrt{1-\eta\frac{\Delta\omega_{\rm{SB}}}{\alpha}}.\end{split}$
(49)

From Eq. (49), one can see that this expression of precessional frequency
modified by inertia converges to the conventional expression of FMR at
$\eta=0$. The approximation ”b” of the nutational frequency reads

${\omega}_{\rm{n}}^{(\rm{b})}=\frac{1}{\eta}+\frac{\Delta\omega_{\rm{SB}}}{2\alpha}.$
(50)

The series expansion leads to a further error of about 6% for the parameters
used for numerical calculations presented in the main part of the paper.

##  Appendix C FREE ENERGY DENSITY

We consider two geometrical configurations are of interest for the study of
resonances in magnetic materials. In the first configuration, a magnetic field
rotates tangentially through the plane of the film surface the _x0y_ plane
(Fig. 3(a)). In the second configuration, the magnetic field goes from the
tangential direction to the normal direction in the _x0y_ plane (Fig. 3(b)).
In the geometries selected here, the films are located differently relative to
the axes, thus the aforementioned planes do not coincide. Such choice of the
axes allows one to avoid the division by zero in the out-of-plane applied
field configuration (Fig. 3(b)). Otherwise, if one directs the magnetic field
along the normal to the film in the configuration of axes shown in Fig. 3(a),
one obtains singularity $\left({{\theta_{0}}=0}\right)$ in Eqs. (4) and (2).
Note that an alternative approach was derived in the past by Baselgia et al.
[47].

Figure 3: The orientation of coordinate frames in the case of magnetic field
applied in-plane (a) and out-of-plane (b).  Figure 4: The sequence of z-y-z
rotations to the $\tilde{\alpha},$ $\tilde{\beta}$ and $\tilde{\gamma}$ angles
respectively.

In order to write expressions of the energy contributions, let us define
coordinate frames. In the general case, the axes of magnetocrystalline
anisotropy (cubic, uniaxial, etc) may not coincide with the demagnetization
axes, therefore, one needs to make a transition from one axis to another to
calculate the energy. We indicate angles of magnetization vector as $\theta$
and $\varphi$ respectively to Cartesian coordinate system _xyz_ , defining the
demagnetizing energy. The polar and azimuthal angles of the vector of the
magnetic field are denoted by ${\theta_{\rm{H}}}$ and ${\varphi_{\rm{H}}}$
with respect to the same axes. The axes specifying the energy of
magnetocrystalline anisotropy are indicated by
${x_{\rm{a}}}{y_{\rm{a}}}{z_{\rm{a}}}$. For instance, we focus on uniaxial
anisotropy in the rotated coordinate system such that the c-axis
$\left({{z_{a}}}\right)$ is aligned with the y-axis. Thus, the free energy
density is given by

$F={F_{\rm{Z}}}+{F_{{\rm{dm}}}}+{F_{\rm{a}}},$ (51)

where ${F_{\rm{Z}}}$ is the Zeeman energy density, ${F_{{\rm{dm}}}}$ is the
demagnetizing energy density and ${F_{\rm{a}}}$ is related to
magnetocrystalline anisotropy. Using representation of vectors in a spherical
coordinate system, one writes the Zeeman energy as

$\begin{split}F_{\rm{Z}}&=-\mu_{0}\bf{MH}\\\
&=-\mu_{0}M_{0}H\left(\right.\sin\theta\sin\theta_{\rm{H}}\cos\left(\varphi-\varphi_{\rm{H}}\right)\\\
&+\cos\theta\cos\theta_{\rm{H}}\left.\right)\end{split}$ (52)

and the demagnetizing energy as

$\begin{split}F_{\rm{dm}}=\frac{1}{2}\mu_{0}M_{0}^{2}\left(\right.N_{\rm{x}}\sin^{2}\theta\cos^{2}\varphi&+N_{\rm{y}}\sin^{2}\theta\sin^{2}\varphi\\\
&+N_{\rm{z}}\cos^{2}\theta\left.\right),\end{split}$ (53)

where ${N_{\rm{x}}},$ ${N_{\rm{y}}},$ ${N_{\rm{z}}}$ are demagnetizing
factors. The magnetocrystalline anisotropy energy of a ferromagnet with cubic
symmetry is given by

$\begin{split}{F_{{\rm{cub}}}}&=K_{{\rm{cub1}}}\left({\alpha_{1}^{2}\alpha_{2}^{2}+\alpha_{2}^{2}\alpha_{3}^{2}+\alpha_{1}^{2}\alpha_{3}^{2}}\right)\\\
&+K_{{\rm{cub2}}}\alpha_{1}^{2}\alpha_{2}^{2}\alpha_{3}^{2}+...\\\\[4.0pt]
&=K_{{\rm{cub1}}}{\cos^{2}}{\varphi_{\rm{a}}}{\sin^{2}}{\theta_{\rm{a}}}\left[{{{\cos}^{2}}{\varphi_{\rm{a}}}}{+\left({1+{{\sin}^{2}}{\theta_{\rm{a}}}}\right){{\sin}^{2}}{\varphi_{\rm{a}}}}\right]\\\
&+K_{{\rm{cub2}}}{\sin^{4}}{\theta_{\rm{a}}}{\sin^{2}}{\varphi_{\rm{a}}}{\cos^{4}}{\varphi_{\rm{a}}}+...,\end{split}$
(54)

where ${\theta_{\rm{a}}}$ and ${\varphi_{\rm{a}}}$ are polar and azimuthal
angles of the magnetization vector in the
${x_{\rm{a}}}{y_{\rm{a}}}{z_{\rm{a}}}$ frame, and ${\alpha_{1}},$
${\alpha_{2}}$ and ${\alpha_{3}}$ are directional cosines with respect to the
same frame. Finally, the uniaxial anisotropy energy can be written as

$\begin{split}{F_{{\rm{uni}}}}&=K_{{\rm{u1}}}{\sin^{2}}{\theta_{\rm{a}}}+K_{{\rm{u2}}}{\sin^{4}}{\theta_{\rm{a}}}+K_{{\rm{u3}}}{\sin^{6}}{\theta_{\rm{a}}}\\\
&+K_{{\rm{u4}}}{\sin^{6}}{\theta_{\rm{a}}}\cos
6{\varphi_{\rm{a}}}+...,\end{split}$ (55)

where constants of anisotropy are denoted with ${K_{i}}$.

The magnetization vector can be specified in two equivalent ways

$\begin{array}[]{l}{\bf{M}}=M\left[\begin{array}[]{l}{\alpha_{1}}\\\
{\alpha_{2}}\\\
{\alpha_{3}}\end{array}\right]=M\left[\begin{array}[]{l}\sin{\theta_{\rm{a}}}\cos{\varphi_{\rm{a}}}\\\
\sin{\theta_{\rm{a}}}\sin{\varphi_{\rm{a}}}\\\
\cos{\theta_{\rm{a}}}\end{array}\right],\end{array}$ (56)

therefore, one can write

$\begin{array}[]{l}{\theta_{\rm{a}}}=\arccos{\alpha_{3}},\\\
{\varphi_{\rm{a}}}=\arctan\dfrac{{{\alpha_{2}}}}{{{\alpha_{1}}}}.\end{array}$
(57)

On the other hand, one can match the vector components in the
${x_{\rm{a}}}{y_{\rm{a}}}{z_{\rm{a}}}$ frame with the $xyz$ frame using the
Euler rotation matrix in the form

$\left[\begin{array}[]{l}{\alpha_{1}}\\\ {\alpha_{2}}\\\
{\alpha_{3}}\end{array}\right]={E_{\rm{u}}}{\left({\tilde{\alpha},\tilde{\beta},\tilde{\gamma}}\right)^{T}}\left[\begin{array}[]{l}\sin\theta\cos\varphi\\\
\sin\theta\sin\varphi\\\ \cos\theta\end{array}\right],$ (58)

$\begin{split}&{E_{\rm{u}}}\left({\tilde{\alpha},\tilde{\beta},\tilde{\gamma}}\right)\\\
&=\left[{\begin{array}[]{*{20}{c}}{{c_{\tilde{\alpha}}}{c_{\tilde{\beta}}}{c_{\tilde{\gamma}}}-{s_{\tilde{\alpha}}}{s_{\tilde{\gamma}}}}&{-{c_{\tilde{\gamma}}}{s_{\tilde{\alpha}}}-{c_{\tilde{\alpha}}}{c_{\tilde{\beta}}}{s_{\tilde{\gamma}}}}&{{c_{\tilde{\alpha}}}{s_{\tilde{\beta}}}}\\\
{{c_{\tilde{\beta}}}{c_{\tilde{\gamma}}}{s_{\tilde{\alpha}}}+{c_{\tilde{\alpha}}}{s_{\tilde{\gamma}}}}&{{c_{\tilde{\alpha}}}{c_{\tilde{\gamma}}}-{c_{\tilde{\beta}}}{s_{\tilde{\alpha}}}{s_{\tilde{\gamma}}}}&{{s_{\tilde{\alpha}}}{s_{\tilde{\beta}}}}\\\
{-{c_{\tilde{\gamma}}}{s_{\tilde{\beta}}}}&{{s_{\tilde{\beta}}}{s_{\tilde{\gamma}}}}&{{c_{\tilde{\beta}}}}\end{array}}\right],\end{split}$
(59)

thereby one rotates the $xyz$ axes to the
${x_{\rm{a}}}{y_{\rm{a}}}{z_{\rm{a}}}$ axes. Note that the coordinate system
rotates, not the vector, therefore the Euler matrix is transposed. Here we
introduce the short notations for trigonometric functions
${c_{\tilde{\alpha}}}=\cos\tilde{\alpha},$
${s_{\tilde{\alpha}}}=\sin\tilde{\alpha}$ and so on. The given Euler rotation
matrix describes a sequence of rotations to the angles $\tilde{\alpha},$
$\tilde{\beta}$ and $\tilde{\gamma}$ around the z, y and z local axes (Fig.
4). Thus, one can express directional cosines ${\alpha_{1}},$ ${\alpha_{2}}$
and ${\alpha_{3}}$ through the predetermined rotation angles $\tilde{\alpha},$
$\tilde{\beta},$ $\tilde{\gamma}$ and angles of magnetization $\theta$ and
$\varphi$ in the $xyz$ frame; then one can use formulas (57) and the
corresponding expression of energy density to calculate anisotropy energy in
the the $xyz$ coordinate system.

For example, we focus on ferromagnets with uniaxial symmetry and cases shown
in Fig. 3, then the consistency between the magnetization angles in $xyz$ and
${x_{\rm{a}}}{y_{\rm{a}}}{z_{\rm{a}}}$ frames is given by

$\begin{split}{E_{\rm{u}}}{\left({\frac{\pi}{2},\frac{\pi}{2},0}\right)^{T}}=\left[{\begin{array}[]{*{20}{c}}0&0&{-1}\\\
{-1}&0&0\\\ 0&1&0\end{array}}\right],\end{split}$ (60)
$\left[\begin{array}[]{l}{\alpha_{1}}\\\ {\alpha_{2}}\\\
{\alpha_{3}}\end{array}\right]=\left[\begin{array}[]{l}-\cos\theta\\\
-\sin\theta\cos\varphi\\\ \sin\theta\sin\varphi\end{array}\right],$ (61)
$\begin{split}{\theta_{\rm{a}}}=\arccos\left({\sin\theta\sin\varphi}\right),\\\
{\varphi_{\rm{a}}}=\arctan\left({\tan\theta\cos\varphi}\right).\end{split}$
(62)

The energy density is defined as

$\begin{split}F&=-\mu_{0}M_{0}H\left(\right.\sin\theta\sin\theta_{H}\cos\left(\varphi-\varphi_{H}\right)\\\
&+\cos\theta\cos\theta_{H}\left.\right)\\\
&+K_{\rm{u1}}\left(1-\sin^{2}\theta\sin^{2}\varphi\right)+F_{\rm{dm}},\end{split}$
(63)

where we neglect the high-order anisotropy terms. For the in-plane
configuration (Fig. 3(a)), one writes

${F_{{\rm{dm}}}}=\frac{{{\mu_{0}}M_{0}^{2}}}{2}{\cos^{2}}\theta,$ (64)

whereas for the out-of-plane case shown in Fig. 3(b), the demagnetizing energy
is given by

${F_{{\rm{dm}}}}=\frac{{{\mu_{0}}M_{0}^{2}}}{2}{\sin^{2}}\theta{\sin^{2}}\varphi.$
(65)

The presented expressions of energy density allow one to avoid the division by
zero in the out-of-plane magnetization configuration
$\left({{\theta_{0}}=0}\right)$. One can calculate the second derivatives of
the energy density and substitute the results in Eq. (4) to find the FMR
frequency modified by inertia or the nutation frequency.

##  Appendix D EQUILIBRIUM ANGLES OF MAGNETIZATION

Figure 5: Equilibrium angles of magnetization. (a-d) polar angles. (e-h)
azimuthal angles.

Based on the presented approach, we find equilibrium angles of magnetization
for all the investigated configurations and the results are plotted in Fig. 5.
Note that the angles for the in-plane cases are calculated for the geometry
shown in Fig. 3(a), while the angles for out-of-plane cases are given in other
axes (Fig. 3(b)).

## References

  * Fulop  _et al._ [2020] J. A. Fulop, S. Tzortzakis, and T. Kampfrath, Laser-driven strong-field terahertz sources, Advanced Optical Materials 8, 1900681 (2020).
  * Seifert  _et al._ [2016] T. Seifert, S. Jaiswal, U. Martens, J. Hannegan, L. Braun, P. Maldonado, F. Freimuth, A. Kronenberg, J. Henrizi, I. Radu, E. Beaurepaire, Y. Mokrousov, P. M. Oppeneer, M. Jourdan, G. Jakob, D. Turchinovich, L. M. Hayden, M. Wolf, M. Münzenberg, M. Kläui, and T. Kampfrath, Efficient metallic spintronic emitters of ultrabroadband terahertz radiation, Nature Photonics 10, 483–488 (2016).
  * Neeraj  _et al._ [2021] K. Neeraj, N. Awari, S. Kovalev, D. Polley, N. Zhou Hagstrom, S. S. P. K. Arekapudi, A. Semisalova, K. Lenz, B. Green, J. C. Deinert, I. Ilyakov, M. Chen, M. Bawatna, V. Scalera, M. Aquino, C. Serpico, O. Hellwig, J.-E. Wegrowe, M. Gensch, and S. Bonetti, Inertial spin dynamics in ferromagnets, Nature Physics 17, 245–250 (2021).
  * Unikandanunni  _et al._ [2021] V. Unikandanunni, R. Medapalli, M. Asa, E. Albisetti, D. Petti, R. Bertacco, E. Fullerton, and S. Bonetti, Inertial spin dynamics in epitaxial cobalt films, arXiv 2109.03076 (2021).
  * Suhl [1998] H. Suhl, Theory of the magnetic damping constant, IEEE Transactions on Magnetics 34, 1834–1838 (1998).
  * Fähnle  _et al._ [2011] M. Fähnle, D. Steiauf, and C. Illg, Generalized Gilbert equation including inertial damping: Derivation from an extended breathing Fermi surface model, Physical Review B 84, 172403 (2011).
  * Fähnle  _et al._ [2013] M. Fähnle, D. Steiauf, and C. Illg, Erratum: Generalized Gilbert equation including inertial damping: Derivation from an extended breathing Fermi surface model [phys. rev. b 84, 172403 (2011)], Physical Review B 88, 219905 (2013).
  * Wegrowe and Ciornei [2012] J.-E. Wegrowe and M.-C. Ciornei, Magnetization dynamics, gyromagnetic relation, and inertial effects, American Journal of Physics 80, 607–611 (2012).
  * Thonig  _et al._ [2017] D. Thonig, O. Eriksson, and M. Pereiro, Magnetic moment of inertia within the torque-torque correlation model, Scientific Reports 7, 931 (2017).
  * Ciornei  _et al._ [2011] M.-C. Ciornei, J. M. Rubí, and J.-E. Wegrowe, Magnetization dynamics in the inertial regime: Nutation predicted at short time scales, Physical Review B 83, 020410 (2011).
  * Mondal  _et al._ [2018] R. Mondal, M. Berritta, and P. M. Oppeneer, Generalisation of Gilbert damping and magnetic inertia parameter as a series of higher-order relativistic terms, Journal of Physics: Condensed Matter 30, 265801 (2018).
  * Bhattacharjee  _et al._ [2012] S. Bhattacharjee, L. Nordstrom, and J. Fransson, Atomistic spin dynamic method with both damping and moment of inertia effects included from first principles, Physical Review Letters 108, 057204 (2012).
  * Mondal  _et al._ [2015] R. Mondal, M. Berritta, K. Carva, and P. M. Oppeneer, Ab initio investigation of light-induced relativistic spin-flip effects in magneto-optics, Physical Review B - Condensed Matter and Materials Physics 91, 174415 (2015).
  * Mondal  _et al._ [2017] R. Mondal, M. Berritta, A. K. Nandy, and P. M. Oppeneer, Relativistic theory of magnetic inertia in ultrafast spin dynamics, Physical Review B 96, 024425 (2017).
  * Olive  _et al._ [2015] E. Olive, Y. Lansac, M. Meyer, M. Hayoun, and J.-E. Wegrowe, Deviation from the Landau-Lifshitz-Gilbert equation in the inertial regime of the magnetization, Journal of Applied Physics 117, 213904 (2015).
  * Cherkasskii  _et al._ [2020] M. Cherkasskii, M. Farle, and A. Semisalova, Nutation resonance in ferromagnets, Physical Review B 102, 184432 (2020).
  * Mondal  _et al._ [2021] R. Mondal, S. Grossenbach, L. Rozsa, and U. Nowak, Nutation in antiferromagnetic resonance, Physical Review B 103, 104404 (2021).
  * Mondal [2021] R. Mondal, Theroy of magnetic inertial dynamics in two-sublattice ferromagnets, Journal of Physics: Condensed Matter 33, 275804 (2021).
  * Mondal and Oppeneer [2021] R. Mondal and P. M. Oppeneer, Influence of intersublattice coupling on the terahertz nutation spin dynamics in antiferromagnets, Physical Review B 104, 104405 (2021).
  * Mondal and Kamra [2021] R. Mondal and A. Kamra, Spin pumping at terahertz nutation resonances, Phys. Rev. B 104, 214426 (2021).
  * Makhfudz  _et al._ [2020] I. Makhfudz, E. Olive, and S. Nicolis, Nutation wave as a platform for ultrafast spin dynamics in ferromagnets, Applied Physics Letters 117, 132403 (2020).
  * Cherkasskii  _et al._ [2021] M. Cherkasskii, M. Farle, and A. Semisalova, Dispersion relation of nutation surface spin waves in ferromagnets, Physical Review B 103, 174435 (2021).
  * Lomonosov  _et al._ [2021] A. M. Lomonosov, V. V. Temnov, and J.-E. Wegrowe, Anatomy of inertial magnons in ferromagnetic nanostructures, Physical Review B 104, 054425 (2021).
  * Li  _et al._ [2015] Y. Li, A.-L. Barra, S. Auffret, U. Ebels, and W. E. Bailey, Inertial terms to magnetization dynamics in ferromagnetic thin films, Phys. Rev. B 92, 140413 (2015).
  * Titov  _et al._ [2021a] S. V. Titov, W. T. Coffey, Y. P. Kalmykov, M. Zarifakis, and A. S. Titov, Inertial magnetization dynamics of ferromagnetic nanoparticles including thermal agitation, Physical Review B 103, 144433 (2021a).
  * Titov  _et al._ [2021b] S. V. Titov, W. T. Coffey, Y. P. Kalmykov, and M. Zarifakis, Deterministic inertial dynamics of the magnetization of nanoscale ferromagnets, Physical Review B 103, 214444 (2021b).
  * Rahman and Bandyopadhyay [2021] R. Rahman and S. Bandyopadhyay, An observable effect of spin inertia in slow magneto-dynamics: Increase of the switching error rates in nanoscale ferromagnets, Journal of Physics Condensed Matter 33, 355801 (2021).
  * Vonsovskii [1966] S. Vonsovskii, Chapter i - magnetic resonance in ferromagnetics, in  _Ferromagnetic Resonance_, edited by S. Vonsovskii (Pergamon, 1966) pp. 1–11.
  * Smit [1955] J. Smit, Ferromagnetic resonance absorption in $\rm{BaFe}_{12}\rm{O}_{19}$, a highly anisotropic crystal, Philips Res. Rep. 10, 113 (1955).
  * Suhl [1955] H. Suhl, Ferromagnetic resonance in nickel ferrite between one and two kilomegacycles, Physical Review 97, 555 (1955).
  * Farle [1998] M. Farle, Ferromagnetic resonance of ultrathin metallic layers, Reports on progress in physics 61, 755 (1998).
  * Skrotskii and Kurbatov [1959] G. V. Skrotskii and L. V. Kurbatov, Theory of the anisotropy of the width of ferromagnetic resonance absorption line, JETP 8, 148–151 (1959).
  * Goryunov  _et al._ [1995] Y. V. Goryunov, N. Garif’yanov, G. Khaliullin, I. Garifullin, L. Tagirov, F. Schreiber, T. Mühge, and H. Zabel, Magnetic anisotropies of sputtered Fe films on MgO substrates, Physical Review B 52, 13450 (1995).
  * Coey [2010] J. M. Coey,  _Magnetism and magnetic materials_ (Cambridge university press, 2010).
  * Finizio  _et al._ [2018] S. Finizio, S. Wintz, D. Bracher, E. Kirk, A. S. Semisalova, J. Förster, K. Zeissler, T. Weßels, M. Weigand, K. Lenz, A. Kleibert, and J. Raabe, Thick permalloy films for the imaging of spin texture dynamics in perpendicularly magnetized systems, Phys. Rev. B 98, 104415 (2018).
  * Ortiz  _et al._ [2021] V. H. Ortiz, B. Arkook, J. Li, M. Aldosary, M. Biggerstaff, W. Yuan, C. Warren, Y. Kodera, J. E. Garay, I. Barsukov, and J. Shi, First- and second-order magnetic anisotropy and damping of europium iron garnet under high strain, Phys. Rev. Materials 5, 124414 (2021).
  * Gonçalves  _et al._ [2018] A. Gonçalves, F. Garcia, H. Lee, A. Smith, P. Soledade, C. Passos, M. Costa, N. Souza-Neto, I. Krivorotov, L. Sampaio, and I. Barsukov, Oscillatory interlayer coupling in spin hall systems, Scientific reports 8, 1 (2018).
  * Etesamirad  _et al._ [2021] A. Etesamirad, R. Rodriguez, J. Bocanegra, R. Verba, J. Katine, I. N. Krivorotov, V. Tyberkevych, B. Ivanov, and I. Barsukov, Controlling magnon interaction by a nanoscale switch, ACS Applied Materials & Interfaces 13, 20288 (2021).
  * Barsukov  _et al._ [2019] I. Barsukov, H. K. Lee, A. A. Jara, Y.-J. Chen, A. M. Gonçalves, C. Sha, J. A. Katine, R. E. Arias, B. A. Ivanov, and I. N. Krivorotov, Giant nonlinear damping in nanoscale ferromagnets, Science Advances 5, eaav6943 (2019).
  * Navabi  _et al._ [2019] A. Navabi, Y. Liu, P. Upadhyaya, K. Murata, F. Ebrahimi, G. Yu, B. Ma, Y. Rao, M. Yazdani, M. Montazeri, L. Pan, I. N. Krivorotov, I. Barsukov, Q. Yang, P. Khalili Amiri, Y. Tserkovnyak, and K. L. Wang, Control of spin-wave damping in YIG using spin currents from topological insulators, Phys. Rev. Applied 11, 034046 (2019).
  * Shaw  _et al._ [2013] J. M. Shaw, H. T. Nembach, T. J. Silva, and C. T. Boone, Precise determination of the spectroscopic g-factor by use of broadband ferromagnetic resonance spectroscopy, Journal of Applied Physics 114, 243906 (2013).
  * Barsukov  _et al._ [2014] I. Barsukov, Y. Fu, A. M. Gonçalves, M. Spasova, M. Farle, L. C. Sampaio, R. E. Arias, and I. N. Krivorotov, Field-dependent perpendicular magnetic anisotropy in CoFeB thin films, Applied Physics Letters 105, 152403 (2014).
  * Barsukov  _et al._ [2015] I. Barsukov, Y. Fu, C. Safranski, Y.-J. Chen, B. Youngblood, A. M. Gonçalves, M. Spasova, M. Farle, J. A. Katine, C. C. Kuo, and I. N. Krivorotov, Magnetic phase transitions in Ta/CoFeB/MgO multilayers, Applied Physics Letters 106, 192407 (2015).
  * Ordonez-Romero  _et al._ [2013] C. L. Ordonez-Romero, M. A. Cherkasskii, N. Qureshi, B. A. Kalinikos, and C. E. Patton, Direct Brillouin light scattering observation of dark spin-wave envelope solitons in magnetic films, Physical Review B 87, 174430 (2013).
  * Wang  _et al._ [2014] Z. Wang, M. Cherkasskii, B. A. Kalinikos, L. D. Carr, and M. Wu, Formation of bright solitons from wave packets with repulsive nonlinearity, New Journal of Physics 16, 053048 (2014).
  * Verhulst [1996] F. Verhulst,  _Nonlinear Differential Equations and Dynamical Systems_ (Springer Berlin Heidelberg, 1996).
  * Baselgia  _et al._ [1988] L. Baselgia, M. Warden, F. Waldner, S. L. Hutton, J. E. Drumheller, Y. Q. He, P. E. Wigen, and M. Maryško, Derivation of the resonance frequency from the free energy of ferromagnets, Physical Review B 38, 2237–2242 (1988).