# Electric-field driven stability control of skyrmions in an ultrathin
transition-metal film

Souvik Paul  s.paul<EMAIL_ADDRESS>Peter
Gr$\ddot{u}$nberg Institut and Institute for Advanced Simulation,
Forschungszentrum J$\ddot{u}$lich and JARA, 52425 J$\ddot{u}$lich, Germany.
Institute of Theoretical Physics and Astrophysics, Christian-Albrechts-
Universit$\ddot{a}$t zu Kiel, Leibnizstrasse 15, 24098 Kiel, Germany  Stefan
Heinze<EMAIL_ADDRESS>Institute of Theoretical Physics and
Astrophysics, Christian-Albrechts-Universit$\ddot{a}$t zu Kiel, Leibnizstrasse
15, 24098 Kiel, Germany

(August 27, 2024)

###### Abstract

To realize future spintronic applications with magnetic skyrmions –
topologically nontrivial swirling spin structures – it is essential to achieve
efficient writing and deleting capabilities of these quasi-particles.
Electric-field assisted nucleation and annihilation is a promising route,
however, the understanding of the underlying microscopic mechanisms is still
limited. Here, we show how the stability of individual magnetic skyrmions in
an ultrathin transition-metal film can be controlled via external electric
fields. We demonstrate based on density functional theory that it is important
to consider the changes of all interactions with electric field, i.e., the
pair-wise exchange, the Dzyaloshinskii-Moriya interaction, the
magnetocrystalline anisotropy energy, and the higher-order exchange
interactions. The energy barriers for electric-field assisted skyrmion writing
and deleting obtained via atomistic spin simulations vary by up to a factor of
three more than the variations of the interactions calculated from first-
principles. This surprising effect originates from the electric-field
dependent size of metastable skyrmions at a fixed magnetic field. The large
changes of lifetimes allow the possibility of electric-field assisted
thermally activated writing and deleting of skyrmions.

Introduction

Magnetic skyrmions Bogdanov1994 ; tokura13  show great promises as information
carriers in future magnetic memory, logic devices, and neuromorphic computing
due to their nanoscale size and ultralow current-driven manipulation, achieved
by spin transfer torque (STT) tokura13 ; tomasello14 ; zhou14 ; junichi13 ;
sampaio13 ; fert13  and spin orbit torque alex2009 ; miron2010 ; miron2011 ;
liu2012 ; woo2017 ; montoya2018 ; callum2020 . However, the main drawback of
these techniques to realize a low-energy-dissipation device is Joule heating,
which destabilizes the skyrmionic bits. The electric-field-induced
manipulation offers an efficient route for creating, deleting and controlling
skyrmions avoiding the heating problem. The energy dissipation can be reduced
by a factor of 100 by using the electric field as compared to STT matsukura15
. Another important benefit of using an electric field is that it can be
applied locally and it does not displace the skyrmionic bits, which is
desirable for encoding information at a particular position.

In spite of recognizing the potential of electric-field-induced manipulation,
only a few experimental studies have been reported on the electric-field-
induced switching of skyrmions and skyrmion bubbles in transition-metal
multilayers hsu17 ; schott17 ; ma18 ; titiksha18  and recently in multiferroic
heterostructures Wang2020 ; Ba2021 . Theoretical studies have addressed either
the variation of magnetic anisotropy upadhyaya15 ; fook16 ; nakatani16  or the
Dzyaloshinskii-Moriya interaction (DMI) yang2018 ; titiksha18 ; desplat2021
directly by the electric field or indirectly due to the electric-field-induced
strain Wang2020 ; Ba2021 . However, it is the interplay of the exchange
interaction, the DMI, and the anisotropy which is responsible for the
properties of magnetic skyrmions Bogdanov1994 . Therefore, in a study on the
influence of the electric field, one needs to take the variation of all these
interactions into account. Moreover, the importance of higher-order exchange
interactions (HOI) beyond the conventional Heisenberg pair-wise exchange
mechanism in ultrathin films for the stability of skyrmions has recently been
demonstrated hoi2020 .

Here, we investigate the stability of isolated magnetic skyrmions in an
ultrathin film from first-principles electronic structures theory. We take an
atomic Fe/Rh bilayer on the Re(0001) surface, Fe/Rh/Re(0001), as a model
system, which is capable of stabilizing isolated skyrmions at external
magnetic fields in absence of an electric field hoi2020 . From density
functional theory (DFT) calculations, we find that the Heisenberg pair-wise
exchange interactions vary by about 15$\%$, the DMI by only 8$\%$, while the
magnetocrystalline anisotropy energy (MAE) varies by about 30$\%$ for an
electric field difference of 1 V/Å. Among the HOI, only the four-site four
spin interaction shows a significant variation of about 20$\%$ for an electric
field difference of 1 V/Å. We analyze the changes of the magnetic interactions
based on the spin-dependent screening of the electric field at the surface.

We study the formation and collapse of individual skyrmions under electric
fields by atomistic spin simulations using the DFT parameters. The energy
barriers preventing the collapse of individual skyrmions varies by about
60$\%$ and the barriers for skyrmion creation by about 30$\%$ for an electric
difference of 1 V/Å. The enhanced energy barriers with respect to the
variation of the magnetic interactions can be explained by a shift of the
critical magnetic field by about 0.6 T of the phase boundary between the
skyrmion and field-polarized state. Thereby, the skyrmion radius varies
significantly with electric field at a given magnetic field value which leads
to the large change of the energy barriers and thereby of the lifetimes which
makes electric-field driven switching feasible.

Figure 1: Effect of electric fields on the spin spiral energy dispersion and
multi-$Q$ states. a Illustration of the Fe monolayer (ML) on a Rh ML on the
Re(0001) surface, denoted as Fe/Rh/Re(0001), exposed to a perpendicular
uniform electric field. In the DFT calculation, the electric field is created
by a charged sheet located at 5.3 Å above the Fe layer. The direction of
electric field lines for $\mathcal{E}$= $-$0.5 V/Å is shown. b Energy
dispersion $E(\mathbf{q})$ of flat spin spirals along two high symmetry
directions ($\overline{\Gamma\mathrm{KM}}$ and $\overline{\Gamma\mathrm{M}}$)
without SOC at $\mathcal{E}$= $+$0.5 V/Å (red), $\mathcal{E}$= 0.0 V/Å
(green), and $\mathcal{E}$= $-$0.5 V/Å (blue). The filled circles represent
DFT data and the solid lines are fits to the Heisenberg model. The filled
diamonds, highlighted by the yellow border, represent the multi-$Q$ (3$Q$ and
$uudd$) states without SOC, which are shown at the q points of the
corresponding single-$Q$ states (see Fig. 6 for spin structures). Inset shows
$E(\mathbf{q})$ of flat cycloidal spin spirals including DMI and MAE along
$\overline{\Gamma\mathrm{M}}$ for zero and two finite electric fields. The
filled circles represent DFT data and the solid lines are fits to the
Heisenberg plus DMI spin model (see text for details). Grey (cyan) shaded
region indicates right (left) rotating spin spirals. Note that the MAE shifts
$E(\mathbf{q})$ by $K$/2 with respect to the FM state.

Results

First-principles calculations. Fig. 1a shows the setup used to perform the DFT
calculations including an electric field based on the FLEUR code fleur  (see
methods for computational details). The electric field is introduced by
placing a charged sheet in the vacuum region of the Fe/Rh/Re(0001) film
weinert2009 ; oba15 . Charge neutrality of the whole system is maintained by
adding or removing the same amount of opposite charge to or from the film. In
this way, a uniform electric field perpendicular to the film plane is
generated. For the electric field strength, we chose values of $\mathcal{E}$=
$\pm$0.5 V/Å as in the experimental work of Ref. hsu17 , which demonstrated
switching of skyrmions in Fe films at these electric field strengths.

We first discuss the energy dispersion $E(\mathbf{q})$ of homogeneous flat
spin spirals without spin-orbit coupling (SOC) obtained via DFT for
Fe/Rh/Re(0001) along two high-symmetric directions of the two-dimensional
Brillouin zone (2DBZ) (Fig. 1b). The magnetic moment $\mathbf{M}_{i}$ at
lattice site $\mathbf{R}_{i}$, of a spin spiral is given by
$\mathbf{M}_{i}=M(\cos{\mathbf{q}\mathbf{R}_{i}},\sin{\mathbf{q}\mathbf{R}_{i}},0)$,
where $M$ is the size of the moment (about $2.9\mu_{\rm B}$ per Fe atom) and
$\mathbf{q}$ is the spin spiral vector. At $\mathcal{E}$= $0$, the
ferromagnetic (FM) state ($\overline{\Gamma}$ point) is energetically lowest.
The Néel state with an angle of 120∘ between adjacent spins
($\overline{\mathrm{K}}$ point) and the row-wise antiferromagnetic (AFM) state
($\overline{\mathrm{M}}$ point) are significantly higher in energy. The
dispersion calculated for $\mathcal{E}$= $\pm$0.5 V/Å shows the same trend as
of the zero field. The electric field induced modification of $E(\mathbf{q})$
is not visible on this scale at small $\mathbf{q}$ (see Supplementary Fig. 1
for a close-up), however, it is significant at large $\mathbf{q}$ with an
energy rise (drop) at $\mathcal{E}$$>$ $0$ ($\mathcal{E}$$<$ $0$). This kind
of modification has been observed earlier in a freestanding Fe monolayer (ML)
and for a Co ML on Pt(111), which can be explained by the spin spiral
($\mathbf{q}$) dependent screening of the electric field oba15 .

Electric field | $J_{1}$ | $J_{2}$ | $J_{3}$ | $J_{4}$ | $J_{5}$ | $J_{6}$ | $J_{7}$ | $B_{1}$ | $Y_{1}$ | $K_{1}$ | $D_{\rm eff}$ | $K$
---|---|---|---|---|---|---|---|---|---|---|---|---
$+$0.5 V/Å | 9.37 |  $-$1.03 | 0.07 |  $-$0.24 | 0.27 |  $-$0.01 |  $-$0.12 |  $-$0.34 | 1.05 |  $-$1.23 | 0.87 | $-0.16$
0.0 V/Å | 8.85 |  $-$0.77 |  $-$0.05 |  $-$0.22 | 0.27 | 0.05 |  $-$0.16 |  $-$0.39 | 1.00 |  $-$1.36 | 0.89 | $-0.20$
$-$0.5 V/Å | 8.08 |  $-$0.55 | 0.03 |  $-$0.20 | 0.23 | 0.01 |  $-$0.12 |  $-$0.33 | 1.00 |  $-$1.53 | 0.94 | $-0.22$
Table 1: Variation of interaction constants with electric fields. Constants
for the $i$-th nearest-neighbor exchange ($J_{i}$), biquadratic exchange
($B_{1}$), three-site four spin exchange ($Y_{1}$), four-site four spin
exchange ($K_{1}$), effective DMI ($D_{\rm eff}$), and the MAE ($K$) obtained
via DFT for Fe/Rh/Re(0001) at zero and two finite electric fields. The
positive sign of $D_{\rm eff}$ indicates a clockwise rotational sense and
negative value of $K$ indicates an in-plane easy magnetization axis (easy
plane). Data of $\mathcal{E}$= 0.0 V/Å are taken from Ref. hoi2020 . Note that
the first three exchange interactions are modified according to Eqs. (8a-c).
For more details, see Supplementary Table 1. All energies are given in meV.

SOC adds two contributions: the MAE and the DMI (inset of Fig. 1b). The easy
magnetization axis of the Fe/Rh bilayer is in the film plane (easy plane) and
the value of the MAE is $K$= $-0.2$ meV/Fe atom that leads to an energy offset
of $K/2$ for spin spirals with respect to the FM state ($\overline{\Gamma}$
point). The MAE changes only slightly upon applying an electric field which is
barely visible in the inset. The DMI arises due to breaking of the inversion
symmetry at the surface and here it favors cycloidal spin spirals with a
clockwise rotational sense (inset of Fig. 1b). For $\mathcal{E}$= $+$0.5 V/Å,
the spin spiral energy minimum is $-$0.05 meV/Fe below the FM state and the
spiral exhibits a pitch of 18.4 nm. For $\mathcal{E}$= $-$0.5 V/Å, the minimum
is $-$0.20 meV/Fe lower and the pitch significantly decreases to 10.2 nm. Note
that the field-induced change of the energy minimum by 0.15 meV/Fe atom is
significant since the Zeeman energy due to an applied magnetic field of 1 T
amounts to about 0.2 meV/Fe atom.

Since higher-order interactions (HOI) beyond pair-wise Heisenberg type
exchange can be important for the magnetic ground state in transition-metal
films pkurz ; rhfeir111 ; 3spin ; spethmann2020 ; hoffmann ; li2020  as well
as for the stability of individual skyrmions hoi2020  and skyrmion lattices
tf1 , we have also calculated their variation upon applying electric fields.
To calculate the biquadratic interaction as well as the three-site and the
four-site four spin interactions hoffmann , we consider three multi-$Q$
states: the two collinear up-up-down-down ($uudd$) states hardrat  and a
three-dimensional noncollinear 3Q-state pkurz  (see methods for details). As
seen in Fig. 1b, the multi-$Q$ states, calculated at zero and two finite
electric fields, are lower in energy compared to the corresponding spin spiral
states (for energy differences see Supplementary Table 2). The energies of the
multi-Q states shift with field in a similar way as the energy dispersion,
i.e., their energy increases (decreases) for positive (negative) electric
field.

The total energy from DFT calculations are used to parametrize an atomistic
spin model which is given by:

$\displaystyle\mathcal{H}=-\sum_{ij}J_{ij}(\textbf{m}_{i}\cdot\textbf{m}_{j})-\sum_{<ij>}\textbf{D}_{ij}\cdot(\textbf{m}_{i}\times\textbf{m}_{j})$
$\displaystyle-\sum_{i}K(m^{z}_{i})^{2}-\sum_{i}\mu_{s}\textbf{B}\cdot\textbf{m}_{i}-B_{1}\sum_{<ij>}(\textbf{m}_{i}\cdot\textbf{m}_{j})^{2}$
$\displaystyle-2\leavevmode\nobreak\
Y_{1}\sum_{<ijk>}(\textbf{m}_{i}\cdot\textbf{m}_{j})(\textbf{m}_{j}\cdot\textbf{m}_{k})$
$\displaystyle-
K_{1}\sum_{<ijkl>}(\textbf{m}_{i}\cdot\textbf{m}_{j})(\textbf{m}_{k}\cdot\textbf{m}_{l})+(\textbf{m}_{i}\cdot\textbf{m}_{l})(\textbf{m}_{j}\cdot\textbf{m}_{k})$
$\displaystyle-(\textbf{m}_{i}\cdot\textbf{m}_{k})(\textbf{m}_{j}\cdot\textbf{m}_{l})$
(1)

where the magnetic moment of Fe at site $i$ is denoted by $\textbf{M}_{i}$ and
$\textbf{m}_{i}$= $\textbf{M}_{i}/M_{i}$. $J_{ij}$, $\textbf{D}_{ij}$,
$\mu_{s}$ and $K$ denote the pair-wise exchange constants, the DMI vectors,
the magnetic moment and the MAE constant, respectively. $B_{1}$ is the
biquadratic constant, $Y_{1}$ and $K_{1}$ are the three-site and four-site
four spin constant, respectively. The higher-order interactions are taken into
account in nearest-neighbor approximation, indicated in the summation by
$<..>$, since they arise from fourth-order perturbation theory hoffmann . The
DMI has been calculated in the effective nearest-neighbor approximation, i.e.,
determined from the slope of the energy contribution due to SOC close to
$\mathbf{q}=0$ Malottki2017a .

The DFT values of the interaction parameters are given in table 1 for zero and
two finite electric fields. At $\mathcal{E}$= $+0.5$ V/Å, the ferromagnetic
nearest-neighbor exchange constant, $J_{1}$, is enhanced by about 6$\%$, while
at $\mathcal{E}$= $-0.5$ V/Å, it decreases by about 9$\%$ with respect to zero
electric field. The absolute change of $J_{1}$ amounts to about 1.3 meV for a
field change by 1 V/Å, which is similar to the value of 1.2 meV, reported from
a DFT study for a Co ML on Pt(111) oba15 . $J_{2}$ becomes more (less)
antiferromagnetic at $\mathcal{E}$$>$ $0$ ($\mathcal{E}$$<$ $0$) with a
variation almost linear with electric field. The pair-wise exchange constants
beyond second neighbors remain fairly constant with electric fields. Among the
HOI, only the four-site four spin interaction, $K_{1}$, varies significantly
with electric field. It shows a decrease of about 10$\%$ at $\mathcal{E}$=
$+0.5$ V/Å and an increase of about 13$\%$ at $\mathcal{E}$= $-0.5$ V/Å.

The effective nearest-neighbor DMI constant, $D_{\rm eff}$, varies by only
0.07 meV upon changing the electric field by 1 V/Å, which amounts to a
relative change of $\Delta D_{\rm eff}/D_{\rm eff}\approx 0.08$. The field-
induced variation of the DMI is still twice larger than the previously
reported value for a MgO/Co/Pt trilayer yang2018 . The MAE changes by 30$\%$
from $\mathcal{E}$= $+0.5$ V/Å to $\mathcal{E}$= $-0.5$ V/Å and it contributes
0.03 meV to the change of the spin spiral energy minimum (Fig. 1b).

To obtain insight into the electric-field induced changes of the magnetic
interactions, we analyze the spin-dependent screening of the electric field in
the film. The modification of the charge density, $\Delta\rho$, at positive
and negative electric fields with respect to zero field is displayed for the
ferromagnetic state in Fig. 2a,b. $\Delta\rho$ is spin-dependent and
therefore, it creates surface electric dipoles which screen the electric field
close to the surface of the film. The direction of the electric dipole at the
surface switches upon changing the sign of the electric field. We observe that
$\Delta\rho$ is sizable at the Fe layer, which quickly decays to small values
in the film and it displays Friedel oscillations, as expected due to screening
at a metal surface mitsui2020 .

Figure 2: Spin-dependent screening of electric fields at the surface. a,b
Electric-field induced spin-dependent charge density difference
$\Delta\rho_{\sigma}(z,\mathcal{E})$=
$\rho_{\sigma}(z,\mathcal{E})-\rho_{\sigma}(z,0)$ along the $z$-direction and
averaged over the film ($xy$) plane at $\mathcal{E}$= $\pm$0.5 V/Å. Solid blue
and red lines denote $\Delta\rho_{\sigma}(z,\mathcal{E})$ for spin-up
($\sigma=\uparrow$) and spin-down ($\sigma=\downarrow$) electrons,
respectively. The position of the charged sheet, Fe, Rh and first two Re
layers of Fe/Rh/Re(0001) film are indicated. c-f Spin-polarized local density
of states (LDOS) of Fe $d_{z^{2}}$, $d_{xz}$, $p_{z}$, $s$, Rh $p_{z}$ and Re
$p_{z}$ orbitals, respectively. Each panel displays LDOS for electric fields
of $+$0.5 V/Å (red), 0 V/Å (green) and $-$0.5 V/Å (blue). The Fermi energy is
indicated as EF.

The largest charge differences, observed in front of the surface (Fig. 2a,b),
originate from the extended $p_{z}$ and $d_{z^{2}}$ states, which can be
understood from the orbitally decomposed spin-resolved LDOS of the Fe atom
(Fig. 2c-f). The decrease (increase) of the spin down density at positive
(negative) electric field in the Fe layer and the vacuum region is due to
minority $p_{z}$ and $d_{z^{2}}$ states, just below the Fermi energy (Fig.
2c,e). A similar field effect arises for the spin up channel in the vacuum
region, which can be attributed to the majority $p_{z}$ states. On the other
hand, the increase of the spin up charge density at the Fe atom for
$\mathcal{E}>0$ (Fig. 2a) is due to an increase of the majority $d$ states as
can be concluded from the change of the integrated charge density within the
muffin-tin spheres (see Supplementary Table 3).

In contrast, the electric-field induced effect is very small for the Rh and Re
atoms, for which only $p_{z}$ orbitals are slightly affected in the LDOS (Fig.
2g-h). Note that the minority $p_{z}$ peak of Re and Rh just below the Fermi
energy is strongly hybridized with the Fe $p_{z}$ state, as can be concluded
from the same energy position and peak shape as well as from the band
structures (see Supplementary Figures 2 to 5). Thereby, the $p_{z}$ state of
Re and Rh is affected as the Fe $p_{z}$ state is filled or emptied due to an
electric field.

The spin-dependent imbalance of surface charge caused by the electric field
leads to a small linear change of the Fe magnetic moment with electric field,
while the induced moments of Rh and Re remain almost unchanged (Supplementary
Figure 6). Note that the Fe moment for spin spirals at zero and finite
electric fields remains fairly constant upon varying q.

The increase of the nearest-neighbor exchange interaction (cf. Table 1) in the
Fe ML can be explained based on the change of the charge density. Upon
increasing the electric field from $\mathcal{E}$= $-0.5$ V/Å to $\mathcal{E}$=
$+0.5$ V/Å, there is a decrease of the total $sp$ charge by 2.5$\times
10^{-3}$ electrons and an increase of the total $d$ charge by 6.8$\times
10^{-4}$ electrons (see Supplementary Table 3). The Slater-Pauling curve
relates the Curie temperature to the number of $d$ electrons in an alloy and
here, we find that such an increase of $d$ electrons enhances the Curie
temperature of Fe, which can be related to an enhanced nearest-neighbor
exchange interaction Takahashi2007 ; oba15 .

The electric-field induced change of the DMI can be understood based on the
electric dipoles at the interfaces of the film. The DMI constant can be
related to the electric dipole moment at the interfaces of magnetic
multilayers jia2020 . Since the screening charge penetrates into the ultrathin
film (Fig. 2a,b), it contributes to the electric dipole moment. Upon changing
from a positive to a negative electric field, switching of the dipole moment
causes the observed increase and decrease of the DMI strength, respectively.
The Re layer is quite well screened from the electric field, therefore, its
contribution to the change of the DMI is quite small despite being a $5d$
transition metal with a large SOC constant (Supplementary Figure 7). A larger
effect is expected from Rh due to the non-negligible screening charge density
at the Fe/Rh interface (Fig. 2a,b). However, its contribution to the total DMI
is small since Rh is a $4d$ transition metal. As a result, the total change of
DMI constant is relatively small compared to that of the exchange constant.
Nevertheless, it has a profound effect on skyrmion switching as discussed
below.

The hybridization between the $p_{z}$ and $d_{z^{2}}$ states of Fe is
significantly affected by the electric field. It has been shown previously
that the electric-field dependent $p$-$d$ hybridization can explain the change
of the magnetocrystalline anisotropy energy Nakamura2009 .

Summarizing the DFT results, we find that the electric field influences the
ground state through the variation of exchange interactions, DMI, and MAE,
which leads to the change of the spin spiral period and the depth of the
energy minimum.

Figure 3: Effect of electric fields on the phase diagram. Zero temperature
phase diagram of Fe/Rh/Re(0001) at a $\mathcal{E}$= $+0.5$ V/Å , b
$\mathcal{E}$= $0$ V/Å , and c $\mathcal{E}$= $-0.5$ V/Å. Energies of the
relaxed spin spiral (SS, black circles), skyrmion lattice (SkX, red circles)
and field-polarized ferromagnetic (FM, green circles) states are shown with
reference to the homogeneous spin spiral (dashed line). The SS, SkX and FM
phases are denoted with blue, red and green background colors, respectively.
Inset of c shows the variation of the critical field $B_{\mathrm{c}}$ (onset
of FM phase) with the electric field. Note that the critical field is
calculated where the energy of the skyrmion lattice (SkX) phase becomes equal
to the FM phase. Figure 4: Effect of electric fields on radius and energy
barriers of isolated skyrmions. a Radius and b energy barriers of isolated
skyrmions in Fe/Rh/Re(0001) with applied magnetic fields at $\mathcal{E}$=
$+0.5$ V/Å (red), $\mathcal{E}$= 0 V/Å (green) and $\mathcal{E}$= $-0.5$ V/Å
(blue). In b, filled circles show collapse barriers, $\Delta E_{\rm col}$,
while filled diamonds mark creation barriers, $\Delta E_{\rm crea}$. Inset of
a shows the skyrmion profiles at $B$= $3$ T for zero and two finite electric
fields and b shows the variation of the critical field $B_{\mathrm{c}}$ (onset
of FM phase) with the electric field. Note that the critical field is
calculated where the energy of isolated skyrmions is equal to the FM phase.

Atomistic spin simulations. Next we show by atomistic spin dynamics
simulations based on Eq. (1) with the DFT parameters that the field-dependent
spin spiral minimum allows to tune the onset of the FM phase, where magnetic
skyrmions are metastable. From the zero-temperature phase diagram (Fig. 3a-c),
we find that, independent of the electric field, a spin spiral ground state
occurs at zero and small magnetic fields, consistent with the spin spiral
minima in Fig. 1b. At $B$$>$ $0.7$ T, the skyrmion lattice phase becomes
energetically favorable. With further increase of the magnetic field, the
field-polarized or ferromagnetic (FM) phase becomes the lowest energy state.
Note that the onset of the FM phase, i.e., the critical field
$B_{\mathrm{c}}$, changes with electric field (inset of Fig. 3c). The change
of $B_{\mathrm{c}}$ (0.56 T) with the electric field (1 V/Å) is close to the
estimated value from the spin spiral energy minimum in Fig. 1b ($\approx$ 0.75
T).

Now we consider individual magnetic skyrmions in the field-polarized phase. We
find that the skyrmion radius increases significantly for $\mathcal{E}$$<$
$0$, while it decreases for $\mathcal{E}$$>$ $0$ (Fig. 4a). At first sight, it
seems surprising that the order of skyrmion radius with the electric field is
reversed with respect to the change of the spin spiral period at the energy
minimum found from DFT (Fig. 1b). However, this is due to the fact that the
critical magnetic field, at which the transition from the skyrmion to the
field-polarized or ferromagnetic (FM) phase occurs, also changes with the
electric field. Since, in an experiment, switching of a skyrmion is performed
at a fixed magnetic field tf2 ; hsu17 , the variation of the skyrmion radius
(see inset of Fig. 4a) strongly affects the energy barriers for skyrmion
creation or annihilation as shown below footnote1 .

Figure 5: Minimum energy path, collapse and creation barriers at zero and
finite electric fields. a Minimum energy paths (MEP) from the initial
(skyrmion) to the final (FM) state in Fe/Rh/Re(0001) at $B$= $3$ T for zero
and two finite electric fields ($\mathcal{E}$= $\pm 0.5$ V/Å). The energy
difference between the saddle point and the skyrmion (FM) state defines the
collapse (creation) barrier. Total and individual energy contributions at $B$=
$3$ T for b the collapse barrier and c for the creation barrier. For energy
decomposition of the full MEP, see Supplementary Fig. 8. Combined exchange (C.
exchange) indicates the sum of pair-wise exchange, biquadratic and three-site
four spin interactions. The four-site four spin interaction is denoted as
4-Spin for brevity. The electric fields $\mathcal{E}$= $+0.5$ V/Å 0 V/Å and
$-0.5$ V/Å are shown as red, green and blue, respectively.

We used the geodesic nudged elastic band (GNEB) method Bessarab2015  to
calculate the minimum energy path (MEP) between an isolated skyrmion and the
FM background (Fig. 5a). We observe that the skyrmions annihilate by the well-
known radial collapse mechanism muckel2021 . From the MEP, we find the energy
barrier $\Delta E_{\rm col}$ protecting the skyrmion from collapsing into the
FM state, i.e., the difference between the saddle point (SP), maximum energy
point on the path, and the initial (skyrmion) state (Fig. 5a). The creation
barrier, $\Delta E_{\rm crea}$, is obtained from the difference between the SP
and the final (FM) state.

It is apparent that $\Delta E_{\rm col}$ decreases as a function of applied
magnetic field similar to the decrease of the skyrmion radius (Fig. 4b). This
is due to the fact that the energy terms which contribute to the barrier, in
particular the DMI, scale with the number of spins in the skyrmion. On the
other hand, $\Delta E_{\rm crea}$ displays only a small nearly linear rise
with applied magnetic field (Fig. 4b). This can be understood from the fact
that it depends on the energy difference of the FM state and the saddle point.
We can relate $\Delta E_{\rm crea}$ to the energy difference between the
skyrmion and the FM state, $\Delta E_{\rm sk-FM}$, and the collapse barrier by
$\Delta E_{\rm crea}=\Delta E_{\rm sk-FM}+\Delta E_{\rm col}$ (Fig. 5a). As
discussed above, $\Delta E_{\rm col}$ decreases with the magnetic field.
However, the FM state becomes more favorable with increasing magnetic field
and $\Delta E_{\rm sk-FM}$ increases (Fig. 5a). The gradual rise of $\Delta
E_{\rm crea}$ with magnetic field is due to these two opposing contributions.

The magnetic field at which both barriers are equal (green curve in Fig. 5a)
marks the transition from the isolated skyrmion to the FM phase. As seen in
the inset of Fig. 4b, this critical magnetic field shifts with the electric
field. At $\mathcal{E}$= $0$ V/Å, $B_{\rm c}$ is about 2.9 T. If one changes
the electric field to $+0.5$ V/Å, at a fixed magnetic field, the creation
barrier is enhanced, while the collapse barrier decreases since $B_{c}$ shifts
to a lower value (Fig. 3a). Therefore, the FM state becomes more favorable
(red curve in Fig. 5a). For $\mathcal{E}$$<$ $0$, the opposite effect occurs
since one moves into the skyrmion phase at a fixed magnetic field of 2.9 T.
The skyrmion state is therefore lower (blue curve in Fig. 5a) than the FM
state and the collapse barrier rises. This explains why collapse and creation
barriers show opposite trends upon variation of the electric field. Note that
the critical fields in Fig. 4b is slightly higher than the fields in Fig. 3c.
The reason is that, in the first case, the critical fields are obtained from
the energy equality of the isolated skyrmion to the FM phase, and in the
second case, from the energy equality of the skyrmion lattice phase to the FM
phase.

Now we quantify the energy barrier at $B$= $3$ T in terms of the magnetic
interactions (Fig. 5b,c). At this value, we find $\Delta E_{\rm col}$$\approx$
$125$ meV for zero electric field. $\Delta E_{\rm col}$ increases by $\approx
35$ meV at $\mathcal{E}$= $-0.5$ V/Å (Fig. 5b). Interestingly, the total
energy barrier arises due to the increase of DMI ($\approx 110$ meV) and the
MAE ($\approx 90$ meV), while both the Zeeman term ($\approx-140$ meV) and the
four-site four spin interaction ($\approx-20$ meV) lead to opposite
contributions and thus reduce the electric field effect footnote3 . The
combined exchange interaction, i.e., the sum of the contributions from pair-
wise exchange, biquadratic and three-site four spin interaction, does not
contribute to the change of energy barriers. For $\mathcal{E}$= $+0.5$ V/Å, we
find a decrease of $\Delta E_{\rm col}$ by $\approx 35$ meV and the magnetic
interactions contribute in an analogous way.

The variation of $\Delta E_{\rm col}$ by about $\pm 30\%$ for $\mathcal{E}$=
$\pm$0.5 V/Å (Fig. 5b) cannot be directly understood from the electric-field
induced changes of the magnetic interactions. Since the magnetic moment does
not change much with the electric field, the large variation of the Zeeman
energy must be directly linked to the change of the skyrmion radius with the
electric field (Fig. 4a). For $\mathcal{E}$= $\pm$0.5 V/Å, the skyrmion radius
varies by 1.3 nm at $B$= $3$ T. We can estimate the relative change of the
number of moments in the skyrmion from $\Delta N_{\rm sk}/N_{\rm
sk}=(R^{2}_{{\rm sk}}(\mathcal{E})-R^{2}_{\rm sk}(0))/R^{2}_{\rm sk}(0)\approx
0.3$, which closely matches the variation of the Zeeman contribution to the
energy barrier $\Delta E_{\rm Zeeman}/E_{\rm Zeeman}\approx 0.3$. The
electric-field dependent skyrmion radius also explains the large relative
change of the DMI term, $\Delta E_{\rm DMI}/E_{\rm DMI}\approx 0.25$, and the
anisotropy term, $\Delta E_{\rm MAE}/E_{\rm MAE}\approx 0.45$, which are about
four times larger than the change of $D_{\rm eff}$ and $K$, respectively
(Table 1). For the four-spin interaction, $K_{1}$, this effect is much reduced
since the barrier contribution of this term depends on the saddle point
structure hoi2020  which is less affected by the electric field.

The creation barrier (Fig. 5c) varies by about 15% for $\mathcal{E}$= $\pm$0.5
V/Å, but displays an opposite trend with respect to $\Delta E_{\rm col}$,
i.e., the barrier decreases for $\mathcal{E}$$<$ $0$ and rises for
$\mathcal{E}$$>$$0$. As stated above, this can be explained based on the MEP
(Fig. 5a). The energy decomposition (Fig. 5c) shows that the DMI plays a minor
role, while MAE and four-spin exchange determine the electric field dependence
footnote2 . The Zeeman term shows an opposite trend and again reduces the
electric field effect. The scaling of anisotropy and Zeeman term with electric
field is also much larger than expected from the electric-field induced
changes of the magnetic interactions and is due to the change of skyrmion
radius as explained above for $\Delta E_{\rm col}$.

Discussion

Based on our results, we can discuss the electric-field assisted switching of
skyrmions in a scanning tunneling microscopy (STM) experiment hsu17 . The
electric current injected from the STM tip results in magnon excitations which
eventually trigger the skyrmion collapse. In Ref. muckel2021 , an effective
temperature of the magnon bath due to single hot electron events was estimated
and it is around $50$ K. The skyrmion collapse and creation rates calculated
based on the Arrhenius law compared well with experimental values muckel2021 .
Therefore, we use it to estimate the effect of the electric field on the
switching.

The Arrhenius law is given by $\tau$= $\tau_{0}\exp{(\Delta E_{\rm col}/k_{\rm
B}T)}$. Since the lifetime $\tau$ is dominated by the energy barrier $\Delta
E_{\rm col}$ at low temperatures, we neglect the variation of the prefactor
$\tau_{0}$. Then, the ratio of lifetimes at $B$= $3$ T is given by
$\tau(\mathcal{E}$= $+0.5\leavevmode\nobreak\ {\rm
V/{\AA}})/\tau(\mathcal{E}$= $0)$= $\exp{(-35\leavevmode\nobreak\ {\rm
meV}/k_{\rm B}T)}$. This leads to a change of the skyrmion lifetime by a
factor of about $2$$\times$$10^{-4}$ at a temperature of $T$= $50$ K.
Therefore, a significant effect in deleting individual skyrmions can be
obtained by a local electric field in the tunnel junction. A similar estimate
for the creation leads to a factor of $8$$\times$$10^{-3}$ which shows that
writing of skyrmions can also be greatly manipulated by an electric field.

We have demonstrated that the stability of isolated skyrmions in an ultrathin
film can be changed significantly by external electric field. The electric-
field induced changes of the magnetic interactions lead to a shift of the
critical magnetic field for the onset of the field-polarized phase, which
exhibits metastable skyrmions. As a consequence, changes in the creation and
collapse barriers are much larger than expected from the variations of the
interactions which lead to the possibility of writing and deleting skyrmions.
Our study shows that it is indispensable to consider all magnetic interactions
to evaluate the electric-field effect for skyrmion stability.

Methods

Density functional theory calculations. All DFT calculations with applied
electric fields were performed using the fleur code fleur . We relaxed the top
three layers of Fe/Rh/Re(0001), i.e., the Fe, the Rh and the first Re layer,
along the $z$-direction by minimizing the atomic force on each atom down to
values of less than 0.04 eV/Å in the presence of $\mathcal{E}$= $\pm$0.5 V/Å.
We chose the generalized gradient approximation (GGA) exchange-correlation
functional as parametrized by Perdew, Burke and Ernzerhof ggapbe , 66 $k$
points in two-dimensional Brillouin zone (2DBZ) and $k_{max}$= $4.0$ a.u.-1
for relaxation. We found that there were almost no changes (below 0.01 a.u.)
of the top three inter-layer distances compared to the zero field values given
in Ref. paul2020 . Therefore, we take the zero electric field interlayer
distances of Ref. paul2020  for calculations including electric fields.
Similar to our result, no electric-field induced relaxation was observed in a
Co ML on the Pt(111) surface oba15 .

To check the existence of a noncollinear ground state and to extract the
Heisenberg pair-wise exchange parameters, we calculated the energy dispersion
of homogeneous flat spin spirals of the form
$\mathbf{M}_{i}=M(\cos{\mathbf{q}\mathbf{R}_{i}},\sin{\mathbf{q}\mathbf{R}_{i}},0)$,
where $\mathbf{M}_{i}$ is the magnetic moment at lattice site $\mathbf{R}_{i}$
and $\mathbf{q}$ is the wave vector in 2DBZ kurz2004 . We used the generalized
Bloch theorem gbt  to compute spin spiral energies within the chemical unit
cell. To study the surface, we chose an asymmetric film of two atomic
overlayers on nine Re substrate layers. Since we used an asymmetric film to
compute the magnetic interactions and our interest is on the electric-field
induced changes in magnetism, we only apply electric fields perpendicular to
the Fe surface. To be consistent with the zero electric field spin spiral
calculations paul2020 , we used the local density approximation (LDA)
exchange-correlation functional form given by Vosko, Wilk and Nusair vwn , a
dense mesh of 44$\times$44 $k$-points in the full 2DBZ and $k_{max}=4.0$
a.u.-1.

The DMI was computed within the first-order perturbation theory heide  on the
self-consistent spin spiral state. The MAE was calculated self-consistently as
described in Ref. Li1990  starting from a self-consistent scalar-relativistic
density. To obtain an accurate value of the MAE, we varied the substrate
layers from 13 to 17 of the asymmetric film.

Higher-order exchange interactions from DFT. The fourth-order perturbative
expansion in the hopping parameter over the Coulomb interaction parameter of
the Hubbard model hubbard1  leads to a two-site four spin (biquadratic,
$B_{1}$), a three-site four spin ($Y_{1}$) hubbard2 ; hubbard3  and a four-
site four spin ($K_{1}$) hoffmann  term (see Eq. (1)). These terms arise due
to hopping of the electron among two (two-site four spin or biquadratic),
three (three-site four spin) and four (four-site four spin) lattice sites. We
evaluate these three higher-order interaction constants directly from DFT
within nearest-neighbor approximation.

Figure 6: Formation and spin structures of the multi-$Q$ states. a
Superposition of two 90∘ spin spirals (1$Q$ states) with clockwise (CW) and
counterclockwise (CCW) rotational sense results in an $uudd$ state ($2Q$
state). b 2D Brillouin zone of hexagonal lattice and reciprocal vectors b1 and
b2. Two high-symmetry directions $\overline{\Gamma\mathrm{K}}$ and
$\overline{\Gamma\mathrm{M}}$ are shown. The q vectors corresponding to the
$uudd$ state along $\overline{\Gamma\mathrm{K}}$ (green circle) and along
$\overline{\Gamma\mathrm{M}}$ (blue circle) as well as the 3$Q$ state at the
$\overline{\mathrm{M}}$ point (red circle) are shown. c Spin structure of the
$uudd$ state along $\overline{\Gamma\mathrm{K}}$, which is formed by a
superposition of two 90∘ spin spirals at q=
$\pm(3/4)\overline{\Gamma\mathrm{K}}$, d spin structure of the $uudd$ state
along $\overline{\Gamma\mathrm{M}}$, which is formed by a superposition of two
90∘ spin spirals at q= $\pm(1/2)\overline{\Gamma\mathrm{M}}$ and e spin
structure of the 3$Q$ state, which is formed by a superpostion of three spin
spirals at the $\overline{\mathrm{M}}$ point. Note that the spin structure of
the $uudd$ states is collinear, whereas it is non-collinear for the 3$Q$
state.

To compute the constants, we consider three multi-$Q$ states (Fig. 6c-e). The
multi-$Q$ state is constructed from the superposition of spin spirals
(single-$Q$) related to symmetry equivalent q vectors of the 2DBZ (Fig. 6a,b).
The single-$Q$ and multi-$Q$ states are energetically degenerate within the
Heisenberg pair-wise spin model. The degeneracy is lifted when the mentioned
three higher-order interactions are taken into account which provides a way to
compute their strength.

We consider two collinear $uudd$ or 2$Q$ states hardrat  and a three-
dimensional noncollinear 3$Q$ state pkurz  to determine the three higher-order
exchange constants (for spin structures see Fig. 6c-e). The three higher-order
exchange interaction constants without SOC are determined from the energy
differences of the single-$Q$ and multi-$Q$ states using the following
equations:

$\displaystyle B_{1}=\frac{3}{32}\Delta
E_{\overline{\mathrm{M}}}^{3Q}-\frac{1}{8}\Delta
E_{\overline{\mathrm{M}}/2}^{uudd}$ (2) $\displaystyle
Y_{1}=\frac{1}{8}(\Delta E_{3\overline{\mathrm{K}}/4}^{uudd}-\Delta
E_{\overline{\mathrm{M}}/2}^{uudd})$ (3) $\displaystyle
K_{1}=\frac{3}{64}\Delta E_{\overline{\mathrm{M}}}^{3Q}+\frac{1}{16}\Delta
E_{3{\overline{\mathrm{K}}}/4}^{uudd}$ (4)

where the total energy differences between multi-$Q$ and corresponding
single-$Q$ (spin spiral) states given by

$\displaystyle\Delta
E_{\overline{\mathrm{M}}}^{3Q}=E_{\overline{\mathrm{M}}}^{3Q}-E_{\overline{\mathrm{M}}}^{SS}$
(5) $\displaystyle\Delta
E_{\overline{\mathrm{M}}/2}^{uudd}=E_{\overline{\mathrm{M}}/2}^{uudd}-E_{\overline{\mathrm{M}}/2}^{SS}$
(6) $\displaystyle\Delta
E_{3\overline{\mathrm{K}}/4}^{uudd}=E_{3{\overline{\mathrm{K}}}/4}^{uudd}-E_{3{\overline{\mathrm{K}}}/4}^{SS}$
(7)

see Fig. 2 and Supplementary Table 2 for energy differences.

Upon taking the higher-order exchange constants into account, we need to
modify the first three Heisenberg pair-wise exchange interaction parameters
obtained from fitting of the homogeneous spin spirals neglecting higher-order
interactions as follows (for a derivation see Ref. hoi2020 ):

$\displaystyle J^{\prime}_{1}$ $\displaystyle=J_{1}-Y_{1}$ (8a) $\displaystyle
J^{\prime}_{2}$ $\displaystyle=J_{2}-Y_{1}$ (8b) $\displaystyle
J^{\prime}_{3}$ $\displaystyle=J_{3}-B_{1}/2$ (8c)

Note that the four-site four spin interaction does not adjust any exchange
interaction parameters, since its contribution to the homogeneous spin spirals
is a constant term of $-$12$K_{1}$.

Due to the large 2D unit cell, the energy of the multi-$Q$ states were
evaluated from asymmetric films consisting of 8 layers in total. However, we
have checked for zero electric field that the energy difference between the 8
and 11 layer film calculations is less than 1 meV.

Atomistic spin dynamics simulations. To relax and calculate the energy of the
spin spirals, skyrmion lattice, FM phases as well as the isolated skyrmions,
we performed spin dynamics simulations based on the Landau-Lifshitz equation:

$\displaystyle\hslash\frac{d\textbf{m}_{i}}{dt}=\frac{\partial
H}{\partial{\textbf{m}}_{i}}\times\textbf{m}_{i}-\alpha\left(\frac{\partial
H}{\partial{\textbf{m}}_{i}}\times\textbf{m}_{i}\right)\times\textbf{m}_{i}$
(9)

where $\hslash$ is the reduced Planck constant, $\alpha$ is the damping
parameter and the Hamiltonian $H$ is defined in Eq. (1). We used a time step
of 0.1 fs, $\alpha$ is varied from 0.05 to 0.1 and the simulations were
carried out over 4 to 6 million steps for relaxation. We solve Eq. (9) by
semi-implicit method as proposed by Mentink et al. mentink2000 .

Geodesic nudged elastic band method. We first create isolated skyrmions in the
field-polarized background, i.e., at a magnetic field above $B_{\mathrm{c}}$
from the theoretical profile Bogdanov1994  and then relax the spin structure
using spin dynamics with the full set of DFT parameters (Table I). We computed
the collapse and creation barriers of isolated skyrmions using the GNEB method
gneb . The method finds the minimum energy path (MEP) connecting the initial
state (IS), i.e., skyrmions, and final state (FS), i.e., FM state, on a
multidimensional energy surface. Within GNEB, an initial path connecting the
IS and FS is created by a chain of images of the system. The objective of the
method is to bring the initial path to MEP via relaxing the intermediate
images. The relaxation is achieved by a force projection scheme. For this, the
effective field is calculated at each images and its local tangent to the path
is replaced by a spring force which maintains a uniform distribution of
images. The maximum energy of MEP corresponds to the saddle point (SP) which
defines the energy barrier separating two stable states. We compute the energy
of the SP accurately using a climbing image (CI) method on top of GNEB.

Data availability

The authors declare that the data supporting the findings of this study are
available within the article and its Supplementary Information files.

Code availability

The atomistic spin dynamics code used in this work is available from the
authors upon a reasonable request.

## References

  * (1) Bogdanov, A. N. & Hubert, A.  The properties of isolated magnetic vortices.  _phys. stat. sol. (b)_ 186, 527–543 (1994).
  * (2) Nagaosa, N. & Tokura, Y.  Topological properties and dynamics of magnetic skyrmions FeRe.  _Nat. Nanotechnol._ 8, 899–911 (2013).
  * (3) Tomasello, R. _et al._ A strategy for the design of skyrmion racetrack memories.  _Sci. Rep._ 4, 6784 (2014).
  * (4) Zhou, Y. & Ezawa, M.  A reversible conversion between a skyrmion and a domain-wall pair in a junction geometry.  _Nat. Commun._ 5, 4652 (2014).
  * (5) Iwasaki, J., Mochizuki, M. & Nagaosa, N.  Universal current-velocity relation of skyrmion motion in chiral magnets.  _Nat. Commun._ 4, 1463 (2013).
  * (6) Sampaio, J., Cros, V., Rohart, S., Thiaville, A. & Fert, A.  Nucleation, stability and current-induced motion of isolated magnetic skyrmions in nanostructures.  _Nat. Nanotechnol._ 8, 839–844 (2013).
  * (7) Fert, A., Cros, V. & Sampaio, J.  Skyrmions on the track.  _Nat. Nanotechnol._ 8, 152–156 (2013).
  * (8) Chernyshov, A. _et al._ Evidence for reversible control of magnetization in a ferromagnetic material by means of spin-orbit magnetic field.  _Nat. Phys._ 5, 656 – 659 (2009).
  * (9) Miron, I. M. _et al._ Current-driven spin torque induced by the rashba effect in a ferromagnetic metal layer.  _Nat. Mater._ 9, 230 – 234 (2010).
  * (10) Miron, I. M. _et al._ Perpendicular switching of a single ferromagnetic layer induced by in-plane current injection.  _Nature_ 476, 189 – 193 (2011).
  * (11) Liu, L. _et al._ Spin-torque switching with the giant spin hall effect of tantalum.  _Science_ 336, 555 – 558 (2012).
  * (12) Woo, S. _et al._ Spin-orbit torque-driven skyrmion dynamics revealed by time-resolved X-ray microscopy.  _Nat. Commun._ 8, 15573 (2017).
  * (13) Montoya, S. A. _et al._ Spin-orbit torque induced dipole skyrmion motion at room temperature.  _Phys. Rev. B_ 98, 104432 (2018).
  * (14) MacKinnon, C. R., Lepadatu, S., Mercer, T. & Bissell, P. R.  Role of an additional interfacial spin-transfer torque for current-driven skyrmion dynamics in chiral magnetic layers.  _Phys. Rev. B_ 102, 214408 (2020).
  * (15) Matsukura, F., Tokura, Y. & Ohno, H.  Control of magnetism by electric fields.  _Nat. Nanotechnol._ 8, 209–220 (2015).
  * (16) Hsu, P.-J. _et al._ Electric-field-driven switching of individual magnetic skyrmions.  _Nat. Nanotechnol._ 12, 123–126 (2017).
  * (17) Schott, M. _et al._ The skyrmion switch: Turning magnetic skyrmion bubbles on and off with an electric field.  _Nano Lett._ 17, 3006–3012 (2017).
  * (18) Ma, C. _et al._ Electric field-induced creation and directional motion of domain walls and skyrmion bubbles.  _Nano Lett._ 19, 353–361 (2019).
  * (19) Srivastava, T. _et al._ Large-voltage tuning of dzyaloshinskii-moriya interactions: A route toward dynamic control of skyrmion chirality.  _Nano Lett._ 18, 4871–4877 (2018).
  * (20) Wang, Y. _et al._ Electric-field-driven non-volatile multi-state switching of individual skyrmions in a multiferroic heterostructure.  _Nat. Commun._ 11, 3577 (2020).
  * (21) Ba, Y. _et al._ Electric-field control of skyrmions in multiferroic heterostructure via magnetoelectric coupling.  _Nat. Commun._ 12, 322 (2021).
  * (22) Upadhyaya, P., Yu, G., Amiri, P. K. & Wang, K. L.  Electric-field guiding of magnetic skyrmions.  _Phys. Rev. B_ 92, 134411 (2015).
  * (23) Fook, H. T., Gan, W. L. & Lew, W. S.  Gateable skyrmion transport via field-induced potential barrier modulation.  _Sci. Rep._ 6, 21099 (2016).
  * (24) Nakatani, Y., Hayashi, M., Kanai, S., Fukami, S. & Ohno, H.  Electric field control of skyrmions in magnetic nanodisks.  _Appl. Phys. Lett._ 108, 1524403 (2016).
  * (25) Yang, H., Boulle, O., Cros, V., Fert, A. & Chshiev, M.  Controlling Dzyaloshinskii-Moriya interaction via chirality dependent atomic-layer stacking, insulator capping and electric field.  _Sci Rep_ 8, 12356 (2018).
  * (26) Desplat, L. _et al._ Mechanism for ultrafast electric-field driven skyrmion nucleation.  _Phys. Rev. B_ 104, L060409 (2021).
  * (27) Paul, S., Haldar, S., von Malottki, S. & Heinze, S.  Role of higher-order exchange interactions for skyrmion stability.  _Nat. Commun._ 11, 4756 (2020).
  * (28) _https://www.flapw.de_ .
  * (29) Weinert, M., Schneider, G., Podloucky, R. & Redinger, J.  FLAPW: applications and implementations.  _J. Phys.: Condens. Matter_ 21, 084201 (2009).
  * (30) Oba, M. _et al._ Electric-field-induced modification of the magnon energy, exchange interaction, and Curie temperature of transition-metal thin films.  _Phys. Rev. Lett._ 114, 107202 (2015).
  * (31) Kurz, P., Bihlmayer, G., Hirai, K. & Blügel, S.  Three-dimensional spin structure on a two-dimensional lattice: Mn/Cu(111).  _Phys. Rev. Lett._ 86, 1106–1109 (2001).
  * (32) Romming, N. _et al._ Competition of dzyaloshinskii-moriya and higher-order exchange interactions in Rh/Fe atomic bilayers on Ir(111).  _Phys. Rev. Lett._ 120, 207201 (2018).
  * (33) Krönlein, A. _et al._ Magnetic ground state stabilized by three-site interactions: $\mathrm{Fe}/\mathrm{Rh}(111)$.  _Phys. Rev. Lett._ 120, 207202 (2018).
  * (34) Spethmann, J. _et al._ Discovery of magnetic single- and triple-$\mathbf{q}$ states in $\mathrm{Mn}/\mathrm{Re}(0001)$.  _Phys. Rev. Lett._ 124, 227203 (2020).
  * (35) Hoffmann, M. & Blügel, S.  Systematic derivation of realistic spin models for beyond-heisenberg solids.  _Phys. Rev. B_ 101, 024418 (2020).
  * (36) Li, W., Paul, S., von Bergmann, K., Heinze, S. & Wiesendanger, R.  Stacking-dependent spin interactions in $\mathrm{Pd}/\mathrm{Fe}$ bilayers on Re(0001).  _Phys. Rev. Lett._ 125, 227205 (2020).
  * (37) S. Heinze _et al._ Spontaneous atomic-scale magnetic skyrmion lattice in two dimensions.  _Nat. Phys._ 7, 713 (2011).
  * (38) Hardrat, B. _et al._ Complex magnetism of iron monolayers on hexagonal transition metal surfaces from first principles.  _Phys. Rev. B_ 79, 094411 (2009).
  * (39) von Malottki, S., Dupé, B., F. Bessarab, P., Delin, A. & Heinze, S.  Enhanced skyrmion stability due to exchange frustration.  _Sci. Rep._ 7, 12299 (2017).
  * (40) Mitsui, T. _et al._ Magnetic friedel oscillation at the Fe(001) surface: Direct observation by atomic-layer-resolved synchrotron radiation ${}^{57}\mathrm{Fe}$ Mössbauer spectroscopy.  _Phys. Rev. Lett._ 125, 236806 (2020).
  * (41) Takahashi, C., Ogura, M. & Akai, H.  First-principles calculation of the Curie temperature Slater-Pauling curve.  _J. Phys.: Condens. Matter_ 19, 365233 (2007).
  * (42) Jia, H., Zimmermann, B., Michalicek, G., Bihlmayer, G. & Blügel, S.  Electric dipole moment as descriptor for interfacial Dzyaloshinskii-Moriya interaction.  _Phys. Rev. Materials_ 4, 024405 (2020).
  * (43) Nakamura, K. _et al._ Giant modification of the magnetocrystalline anisotropy in transition-metal monolayers by an external electric field.  _Phys. Rev. Lett._ 102, 187201 (2009).
  * (44) N. Romming _et al._ Writing and Deleting Single Magnetic Skyrmions.  _Science_ 341, 636 (2013).
  * (45) Note that the spin spiral periods obtained at zero magnetic field in the atomistic spin dynamics simulations are consistent with those from the DFT calculations.
  * (46) Bessarab, P. F., Uzdin, V. M. & Jónsson, H.  Method for finding mechanism and activation energy of magnetic transitions, applied to skyrmion and antivortex annihilation.  _Comput. Phys. Commun._ 196, 335 – 347 (2015).
  * (47) Muckel, F. _et al._ Experimental identification of two distinct skyrmion collapse mechanisms.  _Nat. Phys._ (2021).
  * (48) Note that for $K_{1}<0$, the four-site four spin interaction destabilizes skyrmions hoi2020 .
  * (49) Note that the four-spin term exhibits the same field dependence for creation and collapse barriers since it has a similar contribution to the fm state as to the skyrmion state.
  * (50) Perdew, J. P., Burke, K. & Ernzerhof, M.  Generalized gradient approximation made simple.  _Phys. Rev. Lett._ 77, 3865–3868 (1996).
  * (51) Paul, S. & Heinze, S.  Tailoring magnetic interactions in atomic bilayers of Rh and Fe on Re(0001).  _Phys. Rev. B_ 101, 104408 (2020).
  * (52) Kurz, P., Förster, F., Nordström, L., Bihlmayer, G. & Blügel, S.  Ab initio treatment of noncollinear magnets with the full-potential linearized augmented plane wave method.  _Phys. Rev. B_ 69, 024415 (2004).
  * (53) Sandratskii, L. M.  Energy Band Structure Calculations for Crystals with Spiral Magnetic Structure.  _Phys. Status Solidi B_ 136, 167 (1986).
  * (54) Vosko, S. H., Wilk, L. & Nusair, M.  Accurate spin-dependent electron liquid correlation energies for local spin density calculations: a critical analysis.  _Can. J. Phys._ 58, 1200–1211 (1980).
  * (55) Heide, M., , Bihlmayer, G. & Blügel, S.  Describing Dzyaloshinskii-Moriya spirals from first principles.  _Phys. B: Condens. Matter_ 404, 2678 (2009).
  * (56) Li, C., Freeman, A. J., Jansen, H. J. F. & Fu, C. L.  Magnetic anisotropy in low-dimensional ferromagnetic systems: Fe monolayers on Ag(001), Au(001), and Pd(001) substrates.  _Phys. Rev. B_ 42, 5433–5442 (1990).
  * (57) Hubbard, J.  Electron correlations in narrow energy bands.  _Proc. R Soc. London. Ser. A Mathematical and Physical Sciences_ 276, 238–257 (1963).
  * (58) MacDonald, A. H., Girvin, S. H. & Yoshioka, D.  $\frac{t}{U}$ expansion for the Hubbard model.  _Phys. Rev. B_ 37, 9753–9756 (1988).
  * (59) Takahashi, M.  Half-filled Hubbard model at low temperature.  _J. Phys. C Solid State Phys._ 10, 1289–7301 (1977).
  * (60) H. Mentink, J., V. Tretyakov, M., Fasolino, A., I. Katsnelson, M. & Rasing, T.  Stable and fast semi-implicit integration of the stochastic Landau–Lifshitz equation.  _J. Phys.: Condens. Matter_ 22, 176001 (2010).
  * (61) F. Bessarab, P., M. Uzdin, V. & Jónsson, H.  Method for finding mechanism and activation energy of magnetic transitions, applied to skyrmion and antivortex annihilation.  _Comput. Phys. Commun._ 196, 335–347 (2015).

Acknowledgements

We gratefully acknowledge computing time at the supercomputer of the North-
German Supercomputing Alliance (HLRN) and financial support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) via project no.
414321830 (HE3292/11-1). S.P. graciously acknowledge 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”. S.P.
and S.H. thank G. Bihlmayer and M. Hoffmann for useful discussion. S.H. thanks
M. Goerzen for insightful discussions.

Author contributions

S.H. devised the project. S.P. performed the calculations. S.P. and S.H.
analyzed the results and wrote the paper.

Competing interests

The authors declare no competing financial interests.

Additional information

Supplementary information is available for this paper at


††thanks: contributed to this work equally††thanks: contributed to this work
equally

# Cubic magneto-optic Kerr effect in Ni(111) thin films with and without
twinning

Maik Gaerner  Center for Spinelectronic Materials and Devices, Department of
Physics, Bielefeld University, 33615 Bielefeld, Germany  Robin Silber
IT4Innovations, VŠB-Technical University of Ostrava, Ostrava 70800, Czech
Republic  Nanotechnology Centre, VŠB-Technical University of Ostrava, Ostrava
70800, Czech Republic  Tobias Peters  Center for Spinelectronic Materials and
Devices, Department of Physics, Bielefeld University, 33615 Bielefeld, Germany
Jaroslav Hamrle  Faculty of Mathematics and Physics, Charles University,
Prague 12116, Czech Republic  Timo Kuschel<EMAIL_ADDRESS>Center for Spinelectronic Materials and Devices, Department of Physics,
Bielefeld University, 33615 Bielefeld, Germany

###### Abstract

In most studies utilizing the magneto-optic Kerr effect (MOKE), the detected
change of polarized light upon reflection from a magnetized sample is supposed
to be proportional to the magnetization $\bm{M}$. However, MOKE signatures
quadratic in $\bm{M}$ have also been identified and utilized, e.g., to sense
the structural order in Heusler compounds, to detect spin-orbit torque or to
image antiferromagnetic domains. In our study, we observe a strong anisotropic
MOKE contribution of third order in $\bm{M}$ in Ni(111) thin films, attributed
to a cubic magneto-optic tensor $\propto$ $\bm{M}^{3}$. We further show that
the angular dependence of cubic MOKE (CMOKE) is affected by the amount of
structural domain twinning in the sample. Our detailed study on CMOKE will
open up new opportunities for CMOKE applications with sensitivity to twinning
properties of thin films, e.g. CMOKE spectroscopy and microscopy or time-
resolved CMOKE.

I. Introduction - The magneto-optic Kerr effect (MOKE) [1], being magnetic
circular dichroism and birefringence upon reflection from a magnetized sample,
serves as a powerful tool in the area of thin-film magnetic characterization
[2, 3]. So far, applications such as vectorial magnetometry [4, 5], MOKE
spectroscopy [6, 7] and microscopy [8, 9] as well as time-resolved MOKE [10,
11] mainly rely on linear MOKE (LinMOKE), being linearly proportional to the
magnetization $\bm{M}$.
In addition, quadratic-in-magnetization MOKE (QMOKE) [12, 13], being magnetic
linear dichroism and birefringence upon reflection proportional to
$\bm{M}^{2}$, can be employed to study antiferromagnets, for example in time-
resolved experiments [14, 15, 16]. QMOKE showed to be sensitive to the
structural ordering of Heusler compounds [17, 18] and to spin-orbit torque
effects [19]. Vectorial magnetometry including QMOKE [20, 21, 22, 23, 24],
QMOKE spectroscopy [18, 25, 26, 27, 28], as well as QMOKE microscopy [29, 30]
have been realized.
However, any higher-order MOKE, such as cubic-in-magnetization MOKE (CMOKE),
being proportional to $\bm{M}^{3}$, has only been mentioned very rarely so far
[31, 32, 33] and has not yet been systematically investigated in a combined
theoretic and experimental study. Furthermore, higher-order transport effects
have been detected in magnetoresistive and anomalous Hall effects [34, 35, 36]
and can be expected analogously for magnetothermal and anomalous Nernst
effects. Those effects possess equal symmetries as in the case of MOKE.
However, MOKE offers advantages for spectroscopy, microscopy and time-resolved
experiments that electrical measurements cannot provide such as spectral
range, spatial and time resolution.
In this letter, we report on the observation of third-order-in-magnetization
MOKE in Ni(111) thin films within a combined systematic theoretic and
experimental study. We identify this CMOKE as a threefold in-plane angular
dependence of the magnetically saturated longitudinal MOKE (LMOKE) response
which cannot be attributed to LinMOKE or QMOKE, solely. For the quantitative
description of our experimental data, we expand the Taylor series of the
permittivity tensor up to third order in $\bm{M}$. We analyze the impact of
structural domain twinning and find that the threefold angular dependence of
CMOKE is linearly decreasing with the degree of twinning. Thus, the observed
MOKE contributions will be important for future applications related to
structural domain twinning analysis, such as the investigation of huge
magnetostriction in Ni2MnGa [37].

II. Theoretic description - MOKE can be related to the magnetically perturbed
permittivity of the ferromagnetic layer. Hence, the Kerr angles $\Phi_{s/p}$
for $s$\- and $p$-polarized incident light write [38]

$\displaystyle\Phi_{s}$
$\displaystyle{}=\theta_{s}+i\epsilon_{s}=A_{s}\left(\varepsilon_{yx}-\frac{\varepsilon_{yz}\varepsilon_{zx}}{\varepsilon_{d}}\right)+B_{s}\varepsilon_{zx}\text{
,}$ (1) $\displaystyle\Phi_{p}$
$\displaystyle{}=\theta_{p}+i\epsilon_{p}=-A_{p}\left(\varepsilon_{xy}-\frac{\varepsilon_{xz}\varepsilon_{zy}}{\varepsilon_{d}}\right)+B_{p}\varepsilon_{xz}\text{
,}$ (2)

in which $\theta_{s/p}$ is the Kerr rotation and $\epsilon_{s/p}$ is the Kerr
ellipticity. $A_{s/p}$ and $B_{s/p}$ are the optical weighting factors which
are even and odd functions of the angle of incidence, respectively.
$\varepsilon_{d}$ represents the non-magnetic diagonal permittivity tensor
elements of a cubic crystal structure.
In order to describe QMOKE, the permittivity tensor $\bm{\varepsilon}$ is
typically developed up to second order in $\bm{M}$ [39, 13, 25]. However, to
successfully describe all experimental observations presented in our study, we
have to develop $\bm{\varepsilon}$ up to third order in $\bm{M}$. Using the
Einstein summation, the elements of $\bm{\varepsilon}$ are expressed as

$\displaystyle\varepsilon_{ij}$
$\displaystyle=\varepsilon_{ij}^{(0)}+K_{ijk}M_{k}+G_{ijkl}M_{k}M_{l}+H_{ijklm}M_{k}M_{l}M_{m}\text{
.}$ (3)

Here, $\varepsilon_{ij}^{(0)}$ describes the elements of the non-perturbed
permittivity tensor. $M_{k/l/m}$ are components of the normalized
magnetization $\bm{M}$, and $K_{ijk}$ and $G_{ijkl}$ are the components of the
linear and quadratic magneto-optic (MO) tensors $\bm{K}$ and $\bm{G}$ [40],
respectively. We further add the components $H_{ijklm}$ of the cubic MO five-
rank tensor $\bm{H}$ [33].
Using the Onsager relation
$\varepsilon_{ij}(\bm{M})=\varepsilon_{ij}(-\bm{M})$ as well as symmetry
arguments of cubic crystal structures, these tensors can be simplified in such
a way that $\bm{K}$ is solely determined by one independent parameter $K$,
while $\bm{G}$ is described by two independent parameters, namely
$G_{s}=(G_{11}-G_{12})$ and $2G_{44}$ [39, 28]. The difference of these two
parameters $\Delta G=G_{s}-2G_{44}$ denotes the anisotropic strength of
$\bm{G}$ [13, 21]. $\bm{H}$ possesses two independent parameters [32] which
can be named $H_{123}$ and $H_{125}$ following the derivation in Ref. [32].
Moreover, we introduce an anisotropy parameter $\Delta H=H_{123}-3H_{125}$
which describes the anisotropic strength of $\bm{H}$.
In order to yield the perturbed permittivity tensor for a (111)-oriented cubic
crystal structure, the permittivity tensor is rotated around the
$\hat{x}$-axis of our coordinate system (see Supp. Mat. [41]) by
$\vartheta_{[100]}=45^{\circ}$ and subsequently around the $\hat{y}$-axis by
$\vartheta_{[010]}=\arcsin(1/\sqrt{3})\approx 35.26^{\circ}$. To describe the
in-plane sample rotation in our experiment, we further apply a rotation to the
permittivity tensor by the general angle $\alpha$ around the $\hat{z}$-axis of
our coordinate system. The respective expressions for the permittivity tensor
elements are written down in the Supplemental Materials [41].

III. Eight-directional method - In order to analyze individual MOKE
contributions with different dependencies on $\bm{M}$ experimentally, a
separation algorithm, known as the eight-directional-method [42], can be used.
So far, this method or similar ones have been utilized to characterize (001)-
and (011)-oriented thin films of cubic crystal structure [28, 43] while the
QMOKE in (111)-oriented thin films so far has been identified only by
magnetization loop symmetrization [44] and theoretically [45]. During the
eight-directional method, the MOKE signal is measured for eight different in-
plane magnetization directions $\mu=0^{\circ}+k\cdot
45^{\circ},k=\\{0,1,...,7\\}$ (see Supp. Mat. [41]). From these measurements,
four different MOKE contributions are extracted. Analytical equations for each
of the four MOKE contributions can be derived by inserting the permittivity
tensor elements into the Eqs. (1) and (2). Note that we only consider
contributions up to third order in $\bm{M}$ and stick to in-plane
magnetization.
The first MOKE contribution, separated by the eight-directional method, is the
LinMOKE contribution ($\propto M_{L}$) together with the longitudinal CMOKE
contribution ($\propto M_{L}^{3}$) for the magnetically saturated longitudinal
magnetization component $M_{L}$ (in-plane and parallel to the plane of
incidence)

$\displaystyle\Phi_{M_{L},M_{L}^{3}}$
$\displaystyle=\frac{1}{2}\left(\Phi_{s/p}^{\mu=90^{\circ}}-\Phi_{s/p}^{\mu=270^{\circ}}\right)$
$\displaystyle=\pm B_{s/p}\left(K+\frac{H_{123}+3H_{125}}{2}\right)$
$\displaystyle\quad-A_{s/p}\frac{\sqrt{2}}{6}\left(\Delta H+\frac{K\Delta
G}{\varepsilon_{d}}\right)\sin{(3\alpha)}\text{ .}$ (4)

Here, we can identify an isotropic term $\propto$ $M_{L}$ and $\propto$
$M_{L}^{3}$ as well as an anisotropic threefold term $\propto$ $M_{L}^{3}$.
The amplitude of the anisotropic contribution is stemming from two MO
parameters, being $\Delta H$ and the product $K\Delta G$. The latter stems
from the intermixing of $\varepsilon_{yz}\varepsilon_{zx}$ and
$\varepsilon_{xz}\varepsilon_{zy}$ in Eqs. (1) and (2).
The second MOKE contribution is the transverse CMOKE contribution ($\propto
M_{T}^{3}$) of the magnetically saturated transverse magnetization component
$M_{T}$ (in-plane and perpendicular to the plane of incidence)

$\displaystyle\Phi_{M_{T}^{3}}$
$\displaystyle=\frac{1}{2}\left(\Phi_{s/p}^{\mu=0^{\circ}}-\Phi_{s/p}^{\mu=180^{\circ}}\right)$
$\displaystyle=A_{s/p}\frac{\sqrt{2}}{6}\left(\Delta H+\frac{K\Delta
G}{\varepsilon_{d}}\right)\cos{(3\alpha)}$ (5)

which only consists of an anisotropic threefold term $\propto$ $M_{T}^{3}$. In
contrast to $\Phi_{M_{L},M_{L}^{3}}$, $\Phi_{M_{T}^{3}}$ depends on
$\cos{(3\alpha)}$ instead of $\sin{(3\alpha)}$, but has the same amplitude.
The third MOKE contribution is the QMOKE contribution ($\propto M_{L}M_{T}$)

$\displaystyle\Phi_{M_{L}M_{T}}$
$\displaystyle=\frac{1}{2}\left(\Phi_{s/p}^{\mu=45^{\circ}}+\Phi_{s/p}^{\mu=225^{\circ}}-\Phi_{s/p}^{\mu=135^{\circ}}-\Phi_{s/p}^{\mu=315^{\circ}}\right)$
$\displaystyle=\pm A_{s/p}\left(2G_{44}+\frac{1}{3}\Delta
G-\frac{K^{2}}{\varepsilon_{d}}\right)$ $\displaystyle\quad-
B_{s/p}\frac{\sqrt{2}}{3}\Delta G\sin{(3\alpha)}$ (6)

which consists of an anisotropic threefold term and an isotropic term that
depends on the MO parameters of $\bm{G}$ as well as $K^{2}$.
The fourth and final MOKE contribution, separated by the eight-directional
method, is the QMOKE contribution ($\propto M_{T}^{2}-M_{L}^{2}$)

$\displaystyle\Phi_{M_{T}^{2}-M_{L}^{2}}$
$\displaystyle=\frac{1}{2}\left(\Phi_{s/p}^{\mu=0^{\circ}}+\Phi_{s/p}^{\mu=180^{\circ}}-\Phi_{s/p}^{\mu=90^{\circ}}-\Phi_{s/p}^{\mu=270^{\circ}}\right)$
$\displaystyle=-B_{s/p}\frac{\sqrt{2}}{3}\Delta G\cos{(3\alpha)}$ (7)

which only consists of an anisotropic threefold term. While
$\Phi_{M_{T}^{2}-M_{L}^{2}}$ depends on $\cos{(3\alpha)}$, the anisotropic
part of $\Phi_{M_{L}M_{T}}$ depends on $\sin{(3\alpha)}$, but with the same
amplitude determined by $\Delta G$.

Figure 1: XRD and MOKE data of the Ni(111) samples with almost no (9.8%)
twinning (left) and almost maximum (86%) twinning (right). (a,e) Off-specular
XRD scan (Euler’s cradle texture map) showing the Ni{200} peaks at
$2\theta=51.832^{\circ}$. (b,f) XRD $\theta-2\theta$ scans of the Ni{200}
peaks at selected $\alpha$ angles. (c,g) Eight-directional method performed
with a magnetic field strength of 230 mT using a wavelength of 635 nm. The
dependence of four MOKE contributions on the sample orientation $\alpha$ is
demonstrated. (d,h) Magnetization curves for exemplary sample orientations.
(i) Schematic displaying the eight directions in which the magnetic field is
applied to measure the different MOKE contributions.

IV. Experimental background - The fcc Ni(111) thin films were prepared by
magnetron sputtering onto MgO(111) substrates in an Ar atmosphere of
$2.1\times 10^{-3}$ mbar. Parts of the growth procedure have been adapted from
Ref. [46]. The substrates have been ultrasonically cleaned in acetone,
isopropanol and ethanol and then dried in $N_{2}$. Afterwards, they were
annealed in vacuum at 800∘C. The cleaning and annealing times were varied
throughout the samples series to achieve different degrees of twinning. The
substrates of the exemplarily displayed samples 1 and 2 were annealed for 1h
and 15 minutes, respectively (see Supp. Mat. [41] for more details). About 20
nm thick Ni layers were then deposited onto the substrates with a growth rate
of 0.066 nm/s at a substrate temperature of 350∘C. Afterwards, the Ni layers
were capped with approximately 3 nm of Si which was grown at room temperature.
The Ni layer thicknesses were confirmed by X-ray reflectivity measurements
(see Supp. Mat. [41]), using the Cu $K_{\alpha}$ source of a Phillips X’pert
Pro MPD PW3040-60 diffractometer. Specular X-ray diffraction (XRD)
$\theta$-$2\theta$ scans (see Supp. Mat. [41]) as well as off-specular XRD
texture mapping using an Euler cradle were employed to examine the crystalline
growth of the samples. The texture measurements were performed for 360∘ of
sample rotation $\alpha$ with a tilt of the samples $\Psi=\langle 50^{\circ}$,
$60^{\circ}\rangle$, effectively capturing the Ni{200} diffraction peaks.
Since the Ni layer is grown in a (111) orientation, there are two growth types
possible (structural domain twinning) which differ from each other by an in-
plane rotation of 60∘. Each (111) phase produces a diffraction pattern with
threefold symmetry, but with a $\Delta\alpha=60^{\circ}$ difference between
the patterns.
MOKE measurements were carried out using p-polarized light of 635 nm
wavelength and a 45∘ angle of incidence. Both the samples as well as the
magnetic field were rotated in the sample plane. More detailed descriptions of
the MOKE setup can be found in Refs. [47, 48].

V. Results - The off-specular XRD texture mapping of sample 1 is presented in
Fig. 1(a). The scan was performed at 2$\theta$=51.832∘ and thus shows the
Ni{200} peaks. Additional $\theta$-2$\theta$ scans at all Ni{200} peak
positions were made and are displayed exemplarily in Fig. 1(b). The appearance
of only three intense Ni{200} peaks indicates that almost no twinning occurs
in sample 1. The integrated peak intensities of the dominant phase make up
95.1% of the combined peak intensities which lead to a degree of twinning of
$100\%-(95.1\%-4.9\%)=9.8\%$.
The in-plane angular dependencies of the four MOKE contributions are measured
by the eight-directional method according to Eqs. (4)-(7) and presented in
Fig. 1(c). The experimentally observed dependencies are in good agreement with
the analytical equations. In all four MOKE contributions, pronounced threefold
angular dependencies are visible. While $\Phi_{M_{L},M_{L}^{3}}$ and
$\Phi_{M_{L}M_{T}}$ follow a $\sin{(3\alpha)}$ curve (see Eqs. (4) and (6)),
$\Phi_{M_{T}^{3}}$ and $\Phi_{M_{T}^{2}-M_{L}^{2}}$ describe a
$\cos{(3\alpha)}$ behaviour (see Eqs. (5) and (7)). In addition, the expected
offsets of $\Phi_{M_{L},M_{L}^{3}}$ and $\Phi_{M_{L}M_{T}}$ can be identified
(see Eqs. (4) and (6)), while $\Phi_{M_{T}^{2}-M_{L}^{2}}$ has a vanishing
offset (see Eq. (7)). $\Phi_{M_{T}^{3}}$ shows a finite offset which is not
predicted by theory (see Eq. (5)). Magnetic field misalignments as a possible
origin of this offset have been excluded (see Supp. Mat. [41]).
A small superimposed onefold angular dependence can be seen mainly in the
$\Phi_{M_{L},M_{L}^{3}}$ and $\Phi_{M_{T}^{3}}$ contributions. We attribute it
to a vicinal MOKE (VISMOKE) contribution [38] caused by a slight miscut of the
substrate (see Supp. Mat. [41] for further details). Note that the sample
possesses a uniaxial magnetic anisotropy with the magnetic easy axis oriented
parallel to the step edges of the vicinal surface (see Supp. Mat. [41]).
Therefore, it can be concluded that the three-fold anisotropies in
$\Phi_{M_{L},M_{L}^{3}}$ and $\Phi_{M_{T}^{3}}$ do not stem from the uniaxial
magnetic anisotropy of the sample.
The off-specular XRD texture mapping and $\theta$-2$\theta$ scans of the
Ni{200} peaks of sample 2 are displayed in Figs. 1(e,f). Here, six peaks of
similar intensity can be identified. The dominant phase makes up 57.0% of the
combined peak intensities related to 86% of twinning. In all four MOKE
contributions displayed in Fig. 1(g), the threefold angular dependencies have
undergone a clear reduction, compared to sample 1. This reduction of the Q-
and CMOKE amplitudes is due to the fact that the threefold angular
dependencies of both (111) phases cancel each other out due to their
$\Delta\alpha=60^{\circ}$ difference in orientation. Exemplary magnetization
curves at $\alpha=30^{\circ},90^{\circ}$ are displayed for both samples in
Figs. 1(d) and 1(h), respectively. A stronger anisotropy in the saturation
values can be identified in sample 1 compared to sample 2 due to the lower
degree of twinning.

Figure 2: Dependence of (a) amplitudes and (b) offsets of the individual MOKE
contributions on the degree of structural domain twinning. The dependence of
the amplitudes on the degree of twinning is fitted using a linear regression.

In order to analyze offsets, amplitudes and phase dependencies quantitatively,
fits of the form
$C+A_{1}\sin{(3\alpha+\alpha_{1})}+A_{2}\sin{(\alpha+\alpha_{2})}$ have been
performed (see dashed lines in Figs. 1(c) and 1(g)). Here, the onefold sine
contribution is related to VISMOKE and the phase shift $\alpha_{2}$ correlates
with the miscut direction of the respective Ni film (see Supp. Mat. [41]). The
phase shifts $\alpha_{1}$ of all contributions are in good agreement to each
other with additional 90∘ offset for $\alpha_{1}$ of $\Phi_{M_{T}^{3}}$ and
$\Phi_{M_{T}^{2}-M_{L}^{2}}$ due to the $\cos{(3\alpha)}$ dependence instead
of $\sin{(3\alpha)}$ (see Eqs. (4)-(7)).
As predicted by theory, both CMOKE contributions $\Phi_{M_{L},M_{L}^{3}}$ and
$\Phi_{M_{T}^{3}}$ have similar amplitudes of the threefold angular
dependencies. For sample 1, these are $(1.17\pm 0.01)$ mdeg and $(1.16\pm
0.02)$ mdeg. In sample 2, both of these amplitudes are reduced to $(0.33\pm
0.01)$ mdeg. The QMOKE contributions $\Phi_{M_{L}M_{T}}$ and
$\Phi_{M_{T}^{2}-M_{L}^{2}}$ also show similar amplitudes for each sample. In
sample 1, these are $(0.71\pm 0.03)$ mdeg and $(0.77\pm 0.02)$ mdeg. In sample
2, these amplitudes are reduced to $(0.10\pm 0.01)$ mdeg and $(0.11\pm 0.02)$
mdeg.
An overview of the amplitudes and offsets in all samples is displayed in Fig.
2. It is clearly visible that the amplitudes of the angular dependencies
stemming from CMOKE and QMOKE are directly proportional to the amount of
structural domain twinning in the sample. The anisotropic MOKE contributions
vanish in a fully twinned sample and are maximal in a single crystalline
sample. The offsets of $\Phi_{M_{L},M_{L}^{3}}$ and $\Phi_{M_{L}M_{T}}$ are
independent from the degree of twinning. Slight variations between the offsets
can be attributed to minor variations in layer thicknesses.
Measurements of the Kerr ellipticity and with $\lambda=406$ nm have also been
conducted for two samples using the eight-directional method (see Supp. Mat.
[41]). The same reduction of CMOKE and QMOKE amplitudes with the degree of
twinning has been observed, while the offsets stay constant. In the Supp. Mat.
[41] the dependencies of the MOKE contributions on the magnetic field strength
are also shown. While the offsets of all contributions are constant once the
sample is magnetically saturated, the CMOKE amplitudes slightly decrease with
increasing magnetic field.

VI. Numerical simulations - The CMOKE threefold angular dependence can be
described qualitatively just by $K$ $\Delta G$ (see Eqs. (4) and (5)). In
order to check whether the $\bm{H}$ tensor is needed to quantitatively
describe the amplitude of the CMOKE oscillations, numerical simulations using
Yeh’s transfer matrix formalism [49] with and without the tensor $\bm{H}$ have
been conducted. In the numerical model, the MO parameters are set as free
parameters and the model is fitted to the experimental data.

Figure 3: $\Phi_{M_{L},M_{L}^{3}}$ contribution of sample 1 (VISMOKE
subtracted). (a) Simulations with and without the permittivity tensor of third
order in $\bm{M}$ are shown. (b) Close-up of the simulation without the
$\bm{H}$ tensor.

In Fig. 3(a), simulations with and without $\bm{H}$ tensor are shown for
sample 1. Although the fit was performed for all contributions of the eight-
directional method simultaneously, we only show the results for
$\Phi_{M_{L},M_{L}^{3}}$ here (see Supp. Mat. [41] for the numerical
simulation of all MOKE contributions). The superimposed onefold angular
dependency due to VISMOKE has been subtracted prior to the numerical fitting.
The simulation of MOKE up to second order in $\bm{M}$ (i.e. without $\bm{H}$)
shows negligible angular dependencies of $\Phi_{M_{L},M_{L}^{3}}$ (see Fig.
3(b)). In order to describe $\Phi_{M_{L},M_{L}^{3}}$ solely by $K$ and $\Delta
G$ in the numerical model, the calculated amplitudes for $\Phi_{M_{L}M_{T}}$
and $\Phi_{M_{T}^{2}-M_{L}^{2}}$ would be largely above the experimental
values. However, with the $\bm{H}$ tensor included in the simulation, CMOKE
and the QMOKE contributions can be fitted well, simultaneously. Therefore, the
contribution which is stemming from $K\,\Delta G$ can be assumed to be small
compared to the contribution stemming from $\Delta H$. This shows that the
$\bm{H}$ tensor is indispensable for a quantitative description of CMOKE in
Ni(111) films.

VII. Conclusion - We observed MOKE contributions of third order in $\bm{M}$,
called CMOKE, in Ni(111) thin films. These contributions manifest themselves
in strong threefold angular dependencies of e.g. LMOKE which are even larger
than the ones of the QMOKE contributions. We were able to describe CMOKE
quantitatively by describing the permittivity tensor up to third order in
$\bm{M}$ and including the fifth-rank tensor $\bm{H}$. Furthermore, we showed
that the angular dependencies of the different MOKE contributions are
suppressed in case of structural domain twinning. Our findings can be used in
new applications (spectroscopy, microscopy, pump-probe) based on CMOKE, e.g.
for the analysis of structural domain twinning.

## Acknowledgments

This article has been produced with the financial support of the European
Union under the REFRESH - Research Excellence For Region Sustainability and
High-tech Industries project number CZ.10.03.01/00/22$\\_$003/0000048 via the
Operational Programme Just Transition. Additionally, the infrastructure used
was made available through project No. CZ02.01.01/00/22$\\_$008/0004631-
”Materials and Technologies for Sustainable Development,” funded by the
European Union and the state budget of the Czech Republic within the framework
of the Jan Amos Komensky Operational Program.
J.H. acknowledges support by the European Union project Matfun, Project No.
CZ.02.1.01/0.0/0.0/15$\\_$003/0000487. M.G. acknowledges support from the BMBF
project DiProMag (Plattform MaterialDigital, Project No. 13XP5 120B).

## References

  * J. Kerr [1877] J. Kerr, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science 3, 321 (1877).
  * E. R. Moog and S. D. Bader [1985] E. R. Moog and S. D. Bader, Superlattices and Microstructures 1, 543 (1985).
  * Qiu and Bader [2000] Z. Q. Qiu and S. D. Bader, Review of Scientific Instruments 71, 1243 (2000).
  * Florczak and Dahlberg [1990] J. M. Florczak and E. D. Dahlberg, Journal of Applied Physics 67, 7520 (1990).
  * Vavassori [2000] P. Vavassori, Applied Physics Letters 77, 1605 (2000).
  * Harp  _et al._ [1993] G. R. Harp, D. Weller, T. A. Rabedeau, R. F. C. Farrow, and M. F. Toney, Physical Review Letters 71, 2493 (1993).
  * Veis  _et al._ [2014] M. Veis, L. Beran, R. Antos, D. Legut, J. Hamrle, J. Pistora, C. Sterwerf, M. Meinert, J.-M. Schmalhorst, T. Kuschel, and G. Reiss, Journal of Applied Physics 115, 17A927 (2014).
  * Schäfer [2007] R. Schäfer, ed.,  _Investigation of domains and dynamics of domain walls by the magneto-optical Kerr-effect_ , Handbook of Magnetism and Advanced Magnetic Materials (John Wiley & Sons, Ltd., 2007).
  * J. McCord [2015] J. McCord, Journal of Physics D: Applied Physics 48, 333001 (2015).
  * Beaurepaire  _et al._ [1996] E. Beaurepaire, J.-C. Merle, A. Daunois, and J.-Y. Bigot, Physical Review Letters 76, 4250 (1996).
  * Kirilyuk  _et al._ [2010] A. Kirilyuk, A. V. Kimel, and T. Rasing, Reviews of Modern Physics 82, 2731 (2010).
  * Metzger  _et al._ [1965] G. Metzger, P. Pluvinage, and R. Torguet, Annale de Physique 13, 5 (1965).
  * J. Hamrle  _et al._ [2007] J. Hamrle, S. Blomeier, O. Gaier, B. Hillebrands, H. Schneider, G. Jakob, K. Postava, and C. Felser, Journal of Physics D: Applied Physics 40, 1563 (2007).
  * Saidl  _et al._ [2017] V. Saidl, P. Němec, P. Wadley, V. Hills, R. P. Campion, V. Novák, K. W. Edmonds, F. Maccherozzi, S. S. Dhesi, B. L. Gallagher, F. Trojánek, J. Kuneš, J. Železný, P. Malý, and T. Jungwirth, Nature Photonics 11, 91 (2017).
  * Zheng  _et al._ [2018] Z. Zheng, J. Y. Shi, Q. Li, T. Gu, H. Xia, L. Q. Shen, F. Jin, H. C. Yuan, Y. Z. Wu, L. Y. Chen, and H. B. Zhao, Phys. Rev. B 98, 134409 (2018).
  * Zhao  _et al._ [2021] H. C. Zhao, H. Xia, S. Hu, Y. Y. Lv, Z. R. Zhao, J. He, E. Liang, G. Ni, L. Y. Chen, X. P. Qiu, S. M. Zhou, and H. B. Zhao, Nature Communications 12, 5266 (2021).
  * Wolf  _et al._ [2011] G. Wolf, J. Hamrle, S. Trudel, T. Kubota, Y. Ando, and B. Hillebrands, Journal of Applied Physics 110, 043904 (2011).
  * Silber  _et al._ [2020] R. Silber, D. Král, O. Stejskal, T. Kubota, Y. Ando, J. Pištora, M. Veis, J. Hamrle, and T. Kuschel, Applied Physics Letters 116, 262401 (2020).
  * Montazeri  _et al._ [2015] M. Montazeri, P. Upadhyaya, M. C. Onbasli, G. Yu, K. L. Wong, M. Lang, Y. Fan, X. Li, P. Khalili Amiri, R. N. Schwartz, C. A. Ross, and K. L. Wang, Nature Communications 6, 8958 (2015).
  * Mewes  _et al._ [2004] T. Mewes, H. Nembach, M. Rickart, and B. Hillebrands, Journal of Applied Physics 95, 5324 (2004).
  * T. Kuschel  _et al._ [2011] T. Kuschel, H. Bardenhagen, H. Wilkens, R. Schubert, J. Hamrle, J. Pištora, and J. Wollschläger, Journal of Physics D: Applied Physics 44, 265003 (2011).
  * Tesařová  _et al._ [2012] N. Tesařová, P. Němec, E. Rozkotová, J. Šubrt, H. Reichlová, D. Butkovičová, F. Trojánek, P. Malý, V. Novák, and T. Jungwirth, Applied Physics Letters 100, 102403 (2012).
  * Tesařová  _et al._ [2013] N. Tesařová, P. Němec, E. Rozkotová, J. Zemen, T. Janda, D. Butkovičová, F. Trojánek, K. Olejník, V. Novák, P. Malý, and T. Jungwirth, Nature Photonics 7, 492 (2013).
  * Jang  _et al._ [2020] H. Jang, L. Marnitz, T. Huebner, J. Kimling, T. Kuschel, and D. G. Cahill, Phys. Rev. Applied 13, 024007 (2020).
  * Sepúlveda  _et al._ [2003] B. Sepúlveda, Y. Huttel, C. Martinez Boubeta, A. Cebollada, and G. Armelles, Physical Review B 68, 064401 (2003).
  * Hamrlová  _et al._ [2016] J. Hamrlová, D. Legut, M. Veis, J. Pištora, and J. Hamrle, Journal of Magnetism and Magnetic Materials 420, 143 (2016).
  * R. Silber  _et al._ [2018] R. Silber, M. Tomíčková, J. Rodewald, J. Wollschläger, J. Pištora, M. Veis, T. Kuschel, and J. Hamrle, Photonics and Nanostructures - Fundamentals and Applications 31, 60 (2018).
  * Silber  _et al._ [2019] R. Silber, O. Stejskal, L. Beran, P. Cejpek, R. Antos, T. Matalla-Wagner, J. Thien, O. Kuschel, J. Wollschläger, M. Veis, T. Kuschel, and J. Hamrle, Physical Review B 100, 064403 (2019).
  * Janda  _et al._ [2018] T. Janda, L. Nádvorník, J. Kuchařík, D. Butkovičová, E. Schmoranzerová, F. Trojánek, and P. Němec, Review of Scientific Instruments 89, 073703 (2018).
  * Xu  _et al._ [2019] J. Xu, C. Zhou, M. Jia, D. Shi, C. Liu, H. Chen, G. Chen, G. Zhang, Y. Liang, J. Li, W. Zhang, and Y. Wu, Phys. Rev. B 100, 134413 (2019).
  * Gridnev  _et al._ [1997] V. N. Gridnev, B. B. Krichevtsov, V. V. Pavlov, and R. V. Pisarev, Journal of Experimental and Theoretical Physics Letters 65, 68 (1997).
  * Petukhov  _et al._ [1998] A. V. Petukhov, T. Rasing, T. Katayama, N. Nakajima, and Y. Suzuki, Journal of Applied Physics 83, 6742 (1998).
  * K. Postava  _et al._ [2004] K. Postava, J. Pištora, and T. Yamaguchi, Sensors and Actuators A: Physical 110, 242 (2004).
  * Limmer  _et al._ [2006] W. Limmer, M. Glunk, J. Daeubler, T. Hummel, W. Schoch, R. Sauer, C. Bihler, H. Huebl, M. S. Brandt, and S. T. B. Goennenwein, Phys. Rev. B 74, 205205 (2006).
  * Limmer  _et al._ [2008] W. Limmer, J. Daeubler, L. Dreher, M. Glunk, W. Schoch, S. Schwaiger, and R. Sauer, Phys. Rev. B 77, 205210 (2008).
  * Meyer  _et al._ [2015] S. Meyer, R. Schlitz, S. Geprägs, M. Opel, H. Huebl, R. Gross, and S. T. B. Goennenwein, Applied Physics Letters 106, 132402 (2015).
  * Musiienko  _et al._ [2021] D. Musiienko, F. Nilsén, A. Armstrong, M. Rameš, P. Veřtát, R. H. Colman, J. Čapek, P. Müllner, O. Heczko, and L. Straka, Journal of Materials Research and Technology 14, 1934 (2021).
  * Hamrle  _et al._ [2003] J. Hamrle, J. Ferré, J. P. Jamet, V. Repain, G. Baudot, and S. Rousset, Physical Review B 67, 155411 (2003).
  * Višňovský [2006] Š. Višňovský,  _Optics in Magnetic Multilayers and Nanostructures_ (CRC Press, Boca Raton, 2006).
  * Višňovský [1986] Š. Višňovský, Czechoslovak Journal of Physics 36, 1424 (1986).
  * [41] See Supplemental Materials at [URL will be inserted by publisher] including Refs. [50, 51, 52].
  * Postava  _et al._ [2002] K. Postava, D. Hrabovský, J. Pištora, A. R. Fert, Š. Višňovský, and T. Yamaguchi, Journal of Applied Physics 91, 7293 (2002).
  * Liang  _et al._ [2016] J. H. Liang, Y. L. Chen, L. Sun, C. Zhou, Y. Yang, and Y. Z. Wu, Applied Physics Letters 108, 082404 (2016).
  * P. K. Muduli  _et al._ [2009] P. K. Muduli, W. C. Rice, L. He, B. A. Collins, Y. S. Chu, and F. Tsui, Journal of Physics: Condensed Matter 21, 296005 (2009).
  * Hamrlová  _et al._ [2013] J. Hamrlová, J. Hamrle, K. Postava, and J. Pištora, physica status solidi (b) 250, 2194 (2013).
  * P. Sandström  _et al._ [1999] P. Sandström, E. B. Svedberg, J. Birch, and J.-E. Sundgren, Surface Science 437, L767 (1999).
  * Kehlberger  _et al._ [2015] A. Kehlberger, K. Richter, M. C. Onbasli, G. Jakob, D. H. Kim, T. Goto, C. A. Ross, G. Götz, G. Reiss, T. Kuschel, and M. Kläui, Physical Review Applied 4, 014008 (2015).
  * Müglich  _et al._ [2016] N. D. Müglich, A. Gaul, M. Meyl, A. Ehresmann, G. Götz, G. Reiss, and T. Kuschel, Physical Review B 94, 184407 (2016).
  * P. Yeh [1980] P. Yeh, Surface Science 96, 41 (1980).
  * Parratt [1954] L. G. Parratt, Physical Review 95, 359 (1954).
  * Björck and Andersson [2007] M. Björck and G. Andersson, Journal of Applied Crystallography 40, 1174 (2007).
  * Mende  _et al._ [1983] G. Mende, J. Finster, D. Flamm, and D. Schulze, Surface Science 128, 169 (1983).