Spin-spiral formalism based on the multiple scattering Green's function technique with applications to ultrathin magnetic films and multilayers

TL;DR

This paper introduces a Green's function-based computational method for analyzing spin-spiral magnetic states in layered systems, enabling detailed study of magnetic interactions and effects of spin-orbit coupling in ultrathin films and multilayers.

Contribution

The paper develops a novel computational scheme combining non-relativistic and relativistic Green's function techniques to accurately model spin-spiral states and Dzyaloshinskii-Moriya interactions in layered magnetic systems.

Findings

Successfully describes magnetic ground states of Mn/W(001) and Co/Pt(111).

Reveals changes in magnetic interactions when capping Co/Pt(111) with Ru.

Identifies long-wavelength spin-spiral ground states in multilayer systems.

Abstract

Based on the Korringa-Kohn-Rostoker Green's function technique we present a computational scheme for calculating the electronic structure of layered systems with homogeneous spin-spiral magnetic state. From the self-consistent non-relativistic calculations the total energy of the spin-spiral states is determined as a function of the wave vector, while a relativistic extension of the formalism in first order of the spin-orbit coupling gives an access to the effect of the Dzyaloshinskii-Moriya interactions. We demonstrate that the newly developed method properly describes the magnetic ground state of a Mn monolayer on W(001) and that of a Co monolayer on Pt(111). The obtained spin-spiral energies are mapped to a classical spin model, the parameters of which are compared to those calculated directly from the relativistic torque method. In case of the Co/Pt(111) system we find that the isotropic interaction between the Co atoms is reduced and the Dzyaloshinskii-Moriya interaction is increased when capped by a Ru layer. In addition, we perform spin-spiral calculations on Ir/Fe/Co/Pt and Ir/Co/Fe/Pt multilayer systems and find a spin-spiral ground state with very long wavelength due to the frustrated isotropic couplings between the Fe atoms, whereas the Dzyaloshinskii-Moriya interaction strongly depends on the sequence of the Fe and Co layers.

Tables

Table 1. Table 1: Isotropic exchange couplings for the first five nearest neighbors in a Mn monolayer on W(001) obtained by fitting of the spin-spiral energies to a Heisenberg model and from the relativistic torque method. Note that only J 2 + 2 ​ J 3 subscript 𝐽 2 2 subscript 𝐽 3 J_{2}+2J_{3} could be determined from the fitting procedure. For comparison, the corresponding parameters reported in Ref. Ferriani  et al. , 2008 are also shown as multiplied by two according to the spin model Eq. ( 49 ).

method	$J_1$	$J_2+2J_3$	$J_4$	$J_5$
spin spiral	59.62	-29.27	-1.54	-0.55
RTM	51.38	-26.23	-2.87	0. 81
spin spiral Nandy et al. (2016)	39.4	-11.0	-1.0	-0.30

Table 2. Table 2: Isotropic nearest-neighbor interaction J 1 subscript 𝐽 1 J_{1} and next nearest-neighbor interaction J 2 subscript 𝐽 2 J_{2} between the Co moments, as well as the effective nearest-neighbor interaction J 𝐽 J , the spin stiffness J eff subscript 𝐽 eff J_{\rm eff} and the micromagnetic exchange parameter 𝒥 𝒥 \mathcal{J} for Co/Pt(111) with hcp stacking of the Co monolayer obtained from non-relativistic self-consistent spin-spiral calculations and from the RTM.

method	$J_1$	$J_2$	$J$	$J_{\mathrm{eff}}$	$𝒥$
	(meV)	(meV)	(meV)	(meV$\cdot$ Å2)	(pJ/m)
spin spiral	58.81	2.08	65.30	376.87	47.80
RTM	44.43	2.41	50.15	289.45	36.71

Table 3. Table 3: Effective nearest-neighbor DM interaction D 𝐷 D , spiralization D eff subscript 𝐷 eff D_{\rm eff} and micromagnetic DM parameter 𝒟 𝒟 \mathcal{D} for a Co monolayer deposited on Pt(111) with hcp stacking obtained from fitting to E DM subscript 𝐸 DM E_{\rm DM} in Fig.  3 (b) and directly from the RTM.

method	$D$ (meV)	$D_{\mathrm{eff}}$ (meV$\cdot$ Å)	$𝒟$ (mJ/m2)
spin spiral	2.84	11.82	14.99
RTM	1.84	7.65	9.71

Table 4. Table 4: Interlayer distances for the two considered multilayer structures, Pt/Fe/Co/Ir and Pt/Co/Fe/Ir, after structural relaxations by using the VASP code. All distances are given in units of Å.

	dPt-Pt	dPt-Fe	dFe-Co	dCo-Pt	dIr-Ir
Pt/Fe/Co/Ir	2.56	2.14	2.02	2.10	2.21
	dPt-Pt	dPt-Co	dCo-Fe	dFe-Pt	dIr-Ir
Pt/Co/Fe/Ir	2.54	2.15	1.99	2.12	2.22

Table 5. Table 5: Effective nearest neighbor DM interaction D 𝐷 D and spiralization D eff subscript 𝐷 eff D_{\rm eff} of the CoFe bilayers in the Pt/Fe/Co/Ir and Pt/Co/Fe/Ir multilayers obtained from spin-spiral calculations and from the relativistic torque method with different stackings of the Co layer.

	stacking	method	$D$ (meV)	$D_{\mathrm{eff}}$ (meV$\cdot$ Å)
	Co-fcc	spin spiral	-1.78	7.25
Pt/Fe/Co/Ir	Co-fcc	RTM	-0.45	-1.83
	Co-hcp	RTM	1.11	4.52
	Co-fcc	spin spiral	4.57	-18.60
Pt/Co/Fe/Ir	Co-fcc	RTM	2.66	10.83
	Co-hcp	RTM	3.28	13.35

Equations

$${\mathbf{M}}_{\nu n}=M_{\nu }{\mathbf{S}}_{\nu n}.$$

$${\mathbf{S}}_{\nu n}=\mathcal{R}(\mathbf{n},{\varphi }_n){\mathbf{S}}_{\nu 0},$$

$${\mathbf{S}}_{\nu 0}=\mathbf{n}\cos {\theta }_{\nu }+$$

$${\mathbf{S}}_{\nu n}=\mathbf{n}\cos {\theta }_{\nu }+$$

$$+{\mathbf{e}}_{t2}\sin \left({\mathbf{qT}}_n+{\phi }_{\nu }\right)]\sin {\theta }_{\nu }.$$

$$V_{\nu n}\left(\mathbf{r}\right)=V_{\nu }\left(\mathbf{r}\right),$$

$${\mathbf{B}}_{\nu n}\left(\mathbf{r}\right)={\mathbf{S}}_{\nu n}{B}_{\nu }\left(\mathbf{r}\right)={\mathcal{R}}_n{\mathbf{S}}_{\nu 0}{B}_{\nu }\left(\mathbf{r}\right).$$

$$\displaystyle\left(-\Delta+V_{\nu}\left(\mathbf{r}\right)+\mu_{B}\,\left(\mathcal{R}_{n}^{-1}\mbox{\boldmath$\sigma$}\right)\cdot\mathbf{S}_{\nu 0}\,B_{\nu}\left(\mathbf{r}\right)\right)\psi_{\nu n}\left(\mathbf{r}\right)$$

$$=\epsilon {\psi }_{\nu n}\left(\mathbf{r}\right),$$

$${\mathcal{R}}_n^{-1}\mathbf{\sigma }=U_n^{+}\mathbf{\sigma }{U}_n,$$

$$U_n=\exp \left(\frac{i}{2}\mathbf{n}\mathbf{\sigma }{\varphi }_n\right)=\exp \left(\frac{i}{2}\left(\mathbf{n}\mathbf{\sigma }\right)\left({\mathbf{qT}}_n\right)\right).$$

$$\left(-\Delta +V_{\nu }\left(\mathbf{r}\right)+{\mu }_B{\sigma \cdot \mathbf{S}}_{\nu 0}{B}_{\nu }\left(\mathbf{r}\right)\right)U_n{\psi }_{\nu n}\left(\mathbf{r}\right)$$

$$=\epsilon U_n{\psi }_{\nu n}\left(\mathbf{r}\right)$$

$${\mathcal{Z}}_{\nu n}\left(\epsilon ,\mathbf{r}\right)=\underline{J}\left(\epsilon ,\mathbf{r}\right){\underline{P}}_{\nu n}\left(\epsilon \right)+\sqrt{\epsilon }\mathcal{N}\left(\epsilon ,\mathbf{r}\right),$$

$${\underline{J}}_{\nu n}\left(\epsilon ,\mathbf{r}\right)=\underline{J}\left(\epsilon ,\mathbf{r}\right),$$

$${\underline{P}}_{\nu n}\left(\epsilon \right)={\underline{t}}_{\nu n}{\left(\epsilon \right)}^{-1}-i\sqrt{\epsilon }\underline{I},$$

$$U_n{\mathcal{Z}}_{\nu n}\left(\epsilon ,\mathbf{r}\right)$$

$${\underline{P}}_{\nu n}\left(\epsilon \right)={\underline{U}}_n^{+}{\underline{P}}_{\nu 0}\left(\epsilon \right){\underline{U}}_n,$$

$${\underline{t}}_{\nu n}\left(\epsilon \right)={\underline{U}}_n^{+}{\underline{t}}_{\nu 0}\left(\epsilon \right){\underline{U}}_n.$$

$$\underline{\underline{t}}\left(\epsilon \right)=\{{\delta }_{\nu \mu }{\delta }_{nm}{\underline{t}}_{\nu n}\left(\epsilon \right)\}$$

$${\underline{\underline{G}}}_0\left(\epsilon \right)=\{{\underline{G}}_{0,\nu n,\mu m}\left(\epsilon \right)\}$$

$$\underline{\underline{\tau }}\left(\epsilon \right)={\left[\underline{\underline{t}}{\left(\epsilon \right)}^{-1}-{\underline{\underline{G}}}_0\left(\epsilon \right)\right]}^{-1}.$$

$${\underline{\tau }}_{\nu n,\mu m}\left(\epsilon \right)$$

$$\underline{\underline{\tilde{\tau }}}\left(\epsilon \right)={\left[\underline{\underline{\tilde{t}}}{\left(\epsilon \right)}^{-1}-{\underline{\underline{\tilde{G}}}}_0\left(\epsilon \right)\right]}^{-1},$$

$$\underline{\underline{\tilde{t}}}\left(\epsilon \right)=\{{\delta }_{\nu \mu }{\delta }_{nm}{\underline{t}}_{\nu 0}\left(\epsilon \right)\},$$

$${\underline{\underline{\tilde{G}}}}_0\left(\epsilon \right)=\{{\underline{U}}_n{\underline{G}}_{0,\nu n,\mu m}\left(\epsilon \right){\underline{U}}_m^{+}\}.$$

$${\underline{U}}_n{\underline{G}}_{0,\nu n,\mu m}\left(\epsilon \right){\underline{U}}_m^{+}=$$

$$\frac{1}{{\Omega }_{BZ}}\underset{BZ}{\int }{e}^{i({\mathbf{T}}_m-{\mathbf{T}}_n)\mathbf{k}}{\tilde{\underline{G}}}_{0,\nu \mu }(\mathbf{k},\epsilon )d^{d}k,$$

$$\underline{\widetilde{G}}_{0,\nu\mu}(\mathbf{k},\varepsilon)=\frac{1}{2}\left(\begin{array}[c]{cc}\underline{G}_{0,\nu\mu}(\mathbf{k+}\frac{\mathbf{q}}{2},\varepsilon)(1-n_{3})+\underline{G}_{0,\nu\mu}(\mathbf{k-}\frac{\mathbf{q}}{2},\varepsilon)(1+n_{3})&(n_{1}+in_{2})\left[\underline{G}_{0,\nu\mu}(\mathbf{k+}\frac{\mathbf{q}}{2},\varepsilon)-\underline{G}_{0,\nu\mu}(\mathbf{k-}\frac{\mathbf{q}}{2},\varepsilon)\right]\\ \\ (n_{1}-in_{2})\left[\underline{G}_{0,\nu\mu}(\mathbf{k+}\frac{\mathbf{q}}{2},\varepsilon)-\underline{G}_{0,\nu\mu}(\mathbf{k-}\frac{\mathbf{q}}{2},\varepsilon)\right]&\underline{G}_{0,\nu\mu}(\mathbf{k+}\frac{\mathbf{q}}{2},\varepsilon)(1+n_{3})+\underline{G}_{0,\nu\mu}(\mathbf{k-}\frac{\mathbf{q}}{2},\varepsilon)(1-n_{3})\end{array}\right)\,,$$

$${\underline{\tilde{\tau }}}_{\nu n,\mu m}\left(\epsilon \right)=\frac{1}{{\Omega }_{BZ}}\underset{BZ}{\int }{e}^{i({\mathbf{T}}_m-{\mathbf{T}}_n)\mathbf{k}}{\tilde{\underline{\tau }}}_{\nu \mu }(\mathbf{k},\epsilon )d^{d}k,$$

$$\underline{\underline{\tilde{M}}}(\mathbf{k},\epsilon )=\left\{{\delta }_{\nu \mu }{\underline{t}}_{\nu 0}(\epsilon {)}^{-1}-{\underline{\tilde{G}}}_{0,\nu \mu }(\mathbf{k},\epsilon )\right\}.$$

$$G\left(\mathbf{r}+{\mathbf{R}}_{\nu n},{\mathbf{r}}^{\prime }+{\mathbf{R}}_{\mu m};\epsilon \right)=$$

Full text

I Introduction

Applications of complex magnetic structures in modern information technology are for several decades in the focus of broad research interest Parkin et al. (2008); Fert et al. (2013); Duine et al. (2018). Advanced numerical simulatons based on first principles methods Bihlmayer (2007) play an essential role for understanding the magnetic phenomena on a broad scale and also to design new devices. To this purpose classical spin models Nowak (2007) are also used extensively where the parameters of the spin models can be derived from first principles calculations for collinear or non-collinear magnetic states.

Spin-spiral states form a reach subset of non-collinear magnetic configurations that exist in nature, in particular, in thin magnetic films von Bergmann et al. (2014) and can be studied theoretically with analytical and computational tools. The first ab initio approach for calculating the electronic structure in the presence of a spin-spiral magnetic state in terms of the generalized Bloch theorem was introduced by Sandratskii Sandratskii (1986, 1991) and its implementation in the linearized muffin-tin orbital (LMTO) method was successfully used to calculate the energies of spin-spiral states in bulk Fe Sandratskii and Guletskii (1986); Mryasov et al. (1991). Using the augmented spherical wave (ASW) technique it was possible to study the static non-uniform spin-susceptibiity of various bulk systems Sandratskii (1998). The calculation of spin-spiral states was implemented in the full-potential linearized augmented plane-wave (FLAPW) method as combined with a constrained local moment treatment and used to determine the electronic structure of spin-spiral states in bcc Fe Kurz et al. (2004). In the spin-spiral calculations based on the LMTO, ASW or FLAPW band structure methods an appropriate basis function set is used to solve the Kohn-Sham equations. On the contrary, a non-relativistic multiple scattering Green’s function formalism was developed for spin-spiral configurations in Ref. Mankovsky et al., 2011 and applied to ordered and disordered solids with spiral magnetic order.

For systems with broken space-inversion symmetry, spin orbit coupling (SOC) leads to the appearence of Dzyaloshinskii-Moriya interactions (DMI) Dzyaloshinsky (1958); Moriya (1960) that can stabilize non-collinear chiral magnetic structures, such as spin spirals and magnetic skyrmions. The DMI can be determined from a collinear magnetic structure in terms of the relativistic torque method (RTM) Udvardi et al. (2003); Ebert and Mankovsky (2009); Katsnelson et al. (2010); Mankovsky et al. (2019) or from the spin-cluster expansion (SCE) technique based on the disordered paramagnetic state Szunyogh et al. (2011); Deák et al. (2011). The energy due to the DMI can also be obtained using a first order perturbation treatment of the spin-orbit coupling on top of a non-relativistic spin-spiral calculation Heide et al. (2009) or, at least for commensurate spin spirals, employing supercell calculations Yang et al. (2015); Sandratskii (2017).

The total energy as a function of the wave vector of the spin spirals can be mapped to Heisenberg model giving thus an accurate access to the isotropic coupling parameters in the system. The main advantage of this procedure is that the longitudinal fluctuations of the magnetic moments, including the induced moments, are taken into account Ležaić et al. (2013), while in case of non self-consistent approaches based on the magnetic force theorem Liechtenstein et al. (1987); Jansen (1999), such as in the RTM os SCE, these longitudinal fluctuations are neglected. Morever, the self-consistent spin-spiral calculations include, in principle, all higher order magnetic exchange interactions. There are strong indications that these couplings can stabilize exotic complex magnetic states Heinze et al. (2011); Krönlein et al. (2018).

In this work we present the spin-spiral formalism within the multiple scattering Green’s function technique for both the non-relativistic and the relativistic cases and, as implemented with the screened Korringa-Kohn-Rostoker (SKKR) method Szunyogh et al. (1994); Zeller et al. (1995); Szunyogh et al. (1995a), its applications to ultrathin magnetic films and multilayers. In Sec. II a detailed description of the non-relativistic theory is presented, together with a first-order perturbation technique to include spin-orbit induced effects, while we also give details of the spin model we use for thin magnetic films. The applications are presented in Sec. III. We calculate the spin-spiral dispersion for a Mn monolayer on W(001) and find that the ground state is a right handed cycloidal spin-spiral state according to other theoretical and experimental results Ferriani et al. (2008). For a Co monolayer on Pt(111) we obtain that the ground state is ferromagnetic due to the large isotropic exchange coupling between the Co atoms and that the preferred rotational sense of the in-plane DM vector is left-handed, in agreement with previous theoretical works Dupé et al. (2014); Vida et al. (2016); Zimmermann et al. (2019) and with experiment Hrabec et al. (2014). We show that a Ru overlayer reduces the Co-Co exchange coupling and increases the in-plane DM interaction, thus the Ru/Co/Pt layer sequence can be an important component of the novel functional multilayer structures Karayev et al. (2019); Han et al. (2019). Room temperature skyrmions were observed in Ir/Fe/Co/Pt multilayers Soumyanarayanan et al. (2017); Yagil et al. (2018); Raju et al. (2019) and stable skyrmionic states in $4d$/Fe2/$5d$ multilayers were predicted theoretically Dupé et al. (2016). Motivated by these experimental and theoretical works we investigate Ir/Fe/Co/Pt and Ir/Co/Fe/Pt multilayer systems and highlight that the DMI energy strongly depends on the sequence of the Fe and Co layers. Finally in Sec. IV we summarize our results and draw possible conclusions.

II Theoretical background

II.1 Non-relativistic Green’s function technique for spin-spiral states

In this Section we present the non-relativistic multiple scattering formalism for homogeneous spin spirals along the lines of Ref. Mankovsky et al., 2011. Let us consider a lattice with translation vectors $\mathbf{T}_{n}$ and basis vectors $\mathbf{a}_{\nu}$ ($\nu=1,\dots,M$) defining the inequivalent sublattices of the lattice. For layered systems $\mathbf{T}_{n}$ are two-dimensional (2D) vectors, while $\mathbf{a}_{\nu}$ point to different layers. Note that for complex lattices the third ($z$) component of $\mathbf{a}_{\nu}$ should not necessarily be different for each of the layers. In a homogeneous spin spiral the spin-magnetic moments,

[TABLE]

differ only in their orientations $\mathbf{S}_{\nu n}$ as

[TABLE]

where ${\cal{R}}(\mathbf{n},\phi_{n})$ (later on denoted by ${\cal{R}}_{n}$) stands for a proper rotation around the axis defined by the unit vector $\mathbf{n}$ by an angle $\phi_{n}=\mathbf{qT}_{n}$, where $\mathbf{q}$ is the propagation vector of the spin spiral. If at the reference site of sublattice $\nu$, corresponding to $\mathbf{T}_{0}=\mathbf{0}$, the direction of the spin is

[TABLE]

where $\mathbf{e}_{t1}\perp\mathbf{n}$ and $\mathbf{e}_{t2}=\mathbf{n}\times\mathbf{e}_{t1}$, $\mathbf{S}_{\nu n}$ can be expressed as

[TABLE]

For layered magnetic systems the propagation vector $\mathbf{q}$ lies in the plane of the layers, usually chosen the $(x,y)$ plane of the global frame of reference, and the rotational axis $\mathbf{n}$ is an arbitrary unit vector. Various types of spin spirals are distinguished according to the relative direction of $\mathbf{q}$ and $\mathbf{n}$. Fig. 1 illustrates a Bloch-type spin spiral with $\mathbf{q}\parallel\mathbf{n}$, as well as an in-plane and an out-of-plane Néel-type spin spiral with $\mathbf{q}\perp\mathbf{n}$. For these spin configurations we choose $\theta_{\nu}=\pi/2$, consequently, $\mathbf{S}_{\nu n}\perp\mathbf{n}$, therefore, they represent flat spin spirals.

As the charge density and the magnitude of the magnetization densitiy are identical for each sites in a sublattice, within the local density approximation (LDA) of the density functional theory, the effective potential $V\left(\mathbf{r}\right)$ and exchange field $\mathbf{B}\left(\mathbf{r}\right)$ also preserves the translational-rotational symmetry,

[TABLE]

[TABLE]

The Kohn-Sham equation in cell $n$ of sublattice $\nu$ can then can be written as:

[TABLE]

where we used atomic (Rydberg) units ($\hbar=1$, $2m=1$) and \mbox{\boldmath\sigma}=(\sigma_{x},\sigma_{y},\sigma_{z}) denote the Pauli matrices. The transformation of the Pauli matrices can be expressed as

[TABLE]

where the unitary 2$\times$2 matrix $U_{n}$ is given by

[TABLE]

Thus the Kohn-Sham equation for cell $\nu n$ takes the form,

[TABLE]

which implies that if $\psi_{\nu n}\left(\mathbf{r}\right)$ is a solution of the Schrödinger equation in cell $\nu n$, then $U_{n}\psi_{\nu n}\left(\mathbf{r}\right)$ is a solution for cell $\nu 0$.

Within the non-relativistic multiple scattering theory (MST) the $\left(\ell\,m\,s\right)=\left(L\,s\right)$ angular momentum-spin representation is used with the free-space solutions $J^{Ls}(\varepsilon,\mathbf{r})=j_{\ell}(\sqrt{\varepsilon}r)Y_{\ell m}(\hat{r})\phi_{s}$ and $N^{Ls}(\varepsilon,\mathbf{r})=n_{\ell}(\sqrt{\varepsilon}r)Y_{\ell m}(\hat{r})\phi_{s}$, where $j_{\ell}(x)$ and $n_{\ell}(x)$ are the spherical Bessel- and Neumann functions, respectively, $Y_{\ell m}(\hat{r})$ are spherical harmonics, while $\phi_{s}$ denote spinor basis functions. Introducing the vector notation, $\underline{J}\left(\varepsilon,\mathbf{r}\right)=\left\{J^{Ls}\left(\varepsilon,\mathbf{r}\right)\right\}$ and $\mathcal{N}\left(\varepsilon,\mathbf{r}\right)=\left\{N^{Ls}\left(\varepsilon,\mathbf{r}\right)\right\}$, the regular and irregular scattering solutions of the Kohn-Sham equation Eq. (10) are normalized beyond the radius of the atomic (muffin-tin) sphere as

[TABLE]

and

[TABLE]

respectively, where $\underline{P}_{\nu n}\left(\varepsilon\right)$ is the inverse reactance matrix in the $(L\,s)$ representation, which is connected to the single-site $t$-matrix $\underline{t}_{\nu n}\left(\varepsilon\right)$ as

[TABLE]

with $\underline{I}$ denoting the unit matrix. According to Eq. (10) $U_{n}\,Z^{L\,s}_{\nu n}\left(\varepsilon,\mathbf{r}\right)$ is a solution of the Schrödinger equation at site $\nu 0$, therefore it must be a linear combination of the functions $Z^{L\,s}_{\nu 0}\left(\varepsilon,\mathbf{r}\right)$. Taking into account the boundary condition Eq. (11) this implies

[TABLE]

and similar for the functions $\mathcal{J}_{\nu n}\left(\varepsilon,\mathbf{r}\right)$, together with the transformation for the inverse reactance matrices,

[TABLE]

where we introduced the matrix $\underline{U}_{n}=\{\delta_{LL^{\prime}}U_{n}^{ss^{\prime}}\}$. Consequently, the same transformation applies to the single-site $t$-matrices,

[TABLE]

The following matrices in composite site-angular momentum-spin space,

[TABLE]

and

[TABLE]

are used to calculate the matrix of the scattering-path operator (SPO),

[TABLE]

In Eq. (18) $\underline{G}_{0,\nu n,\mu m}\left(\varepsilon\right)=\{\delta_{ss^{\prime}}{G}^{LL^{\prime}}_{0,\nu n,\mu m}\left(\varepsilon\right)\}$ denote the spin-independent free-space (bare) structure constants.

Utilizing the transformation of the $t$-matrices Eq. (16), the SPO matrix can be expressed as

[TABLE]

where the $\underline{\widetilde{\tau}}_{\nu n,\mu m}\left(\varepsilon\right)$ matrices can be calculated as

[TABLE]

with $\underline{\underline{\widetilde{t}}}\left(\varepsilon\right)$ comprising $t$-matrices that are identical within each sublattice

[TABLE]

and

[TABLE]

Exploiting the lattice Fourier transform of the free-space structure constants and Eq. (9), the above quantity can be expressed

[TABLE]

where $d$=2 or 3 for two-dimensional or three-dimensional translational invariance, respectively, and $\Omega_{BZ}$ is the volume of the Brillouin zone. The matrix $\widetilde{\underline{G}}_{0,\nu\mu}(\mathbf{k},\varepsilon)$ can be written in spin space as

[TABLE]

where in the above expression $\underline{G}_{0,\nu\mu}(\mathbf{k},\varepsilon)$ stands for the lattice Fourier transform of the free-space structure constants that can be obtained via Ewald summation Ewald (1921); Ham and Segall (1961); Kambe (1967a, b, 1968). For the case of a rotational axis parallel to the $z$ axis ($n_{1}=0$, $n_{2}=0$, $n_{3}=1$) $\underline{\widetilde{G}}_{0,\nu\mu}(\mathbf{k},\varepsilon)$ takes the diagonal form as reported in Ref. Mankovsky et al., 2011.

Utilizing also the translational invariance of the $\widetilde{t}$-matrices in Eq. (22), the matrices $\underline{\widetilde{\tau}}_{\nu n,\mu m}\left(\varepsilon\right)$ can be evaluted as

[TABLE]

where $\widetilde{\underline{\tau}}_{\nu\mu}(\mathbf{k},\varepsilon)$ are the blocks of the inverse of the following matrix in sublattice-angular momentum-spin space,

[TABLE]

The SPO-matrix can then be obtained from the transformation Eq. (20).

The Green’s function within the MST is given as

[TABLE]

with $\mathbf{r}_{<}$ and $\mathbf{r}_{>}$ denoting $\mathbf{r}$ or $\mathbf{r}^{\prime}$ with the smaller and larger magnitude, respectively, while the superscript $\times$ stands for the functions when replacing $Y_{L}(\hat{r})\phi_{s}$ by $Y_{L}(\hat{r})^{*}\phi_{s}^{+}$ in Eqs. (11) and (12). It is straightforward to show from Eqs. (14) and (20) that the site-diagonal Green’s function transforms as

[TABLE]

which immediately implies that the charge and magnetization densities are the same in each atomic cell of a sublattice and the spin-magnetic moments rotate from site to site according to Eq. (2). This means that the homogeneous spin-spiral state treated in terms of generalized Bloch theorem is consistent with the non-relativistic density functional theory and a self-consistent calculation for a given wave-vector ${\mathbf{q}}$ can be performed on the cost of a calculation of a periodic collinear magnetic state. It should be noted that if ${\mathbf{S}}_{\nu 0}\perp{\mathbf{n}}$ does not apply, i.e. in case of conical spin spirals, the spin configuration does not correspond to a stationary state, therefore, the self-consistent electronic structure should be determined by exerting appropriate transversal constraint to the local moments.

II.2 Relativistic correction to the spin-spiral energy

Within a relativistic theory the above formalism can not be used for calculating the electronic structure of a spin-spiral state. Formally, one can see this by using the fact that, even neglecting orbital polarization effects, the operator $W_{n}=\exp\left(\frac{i}{\hbar}\left(\mathbf{n}\bm{J}\right)\left(\mathbf{qT}_{n}\right)\right)$, $\bm{J}=\bm{L}+\bm{S}$ being the total angular momentum operator, must be used to describe the rotation of the spin-moments, which precludes to express the structure constants, $\widetilde{\underline{G}}_{0,\nu n,\mu m}(\varepsilon)$, via lattice Fourier transformation as in Eq. (24). In more feasible terms, the magnetic anisotropy induced by the spin-orbit coupling (SOC) differentiates between different directions of the magnetic moments in the spin spiral, thus alters the electronic states from site to site. It is well-known also from spin-model simulations that in the presence of magnetic anisotropy an inhomogeneous spin spiral is formed Hasselberg et al. (2015). Consequently, when including spin-orbit coupling, supercell calculations as based on constrained local moments are needed to treat homogeneous spin-spiral states self-consistently Yang et al. (2015).

It is appropriate to treat the SOC as perturbation in spin-spiral calculations Heide et al. (2009). Here we present a numerical scheme based on the magnetic force theorem to include the effect of SOC in first order. Instead of calculating the total energy of the spin spiral, we will consider the zero-temperature grand potential of the electronic system

[TABLE]

where $E_{F}$ denotes the Fermi energy, $n(\varepsilon)$ is the density of states and $N(\varepsilon)$ is the integrated density of states. Employing Lloyd’s formula Lloyd (1967), the grand potential can be expressed within MST as

[TABLE]

Introducing the difference between the inverse $t$-matrices in the presence of SOC, $\underline{\underline{t}}^{\prime}\left(\varepsilon\right)$, and in the absence of SOC, $\underline{\underline{t}}\left(\varepsilon\right)$,

[TABLE]

the change of the grand potential to first order in $\Delta\underline{\underline{m}}\left(\varepsilon\right)$ can be expressed as, see Ref. Udvardi et al. (2003),

[TABLE]

where $\underline{\tau}_{\nu n,\nu n}\left(\varepsilon\right)$ are the SPO matrices related to the spin spiral within the non-relativistic theory.

In the presence of SOC, we write the Hamiltonian in cell $\nu n$ as

[TABLE]

where $H_{\nu n}$ is the non-relativistic Kohn-Sham Hamiltonian and, for spherical potentials, the Hamiltonian of the spin-orbit coupling is given by

[TABLE]

with

[TABLE]

where $\mathbf{L}$ and $\mathbf{S}$ denote the operators of the electron’s angular momentum and spin, respectively. In order to evaluate the change of the $t$-matrix due to SOC we use the regular solutions,

[TABLE]

normalized beyond the radius of the atomic sphere $S_{\nu}$ as

[TABLE]

where $\mathcal{H}\left(\mathbf{r},\varepsilon\right)$ are the Hankel-type solutions of the free-space Schrödinger equation. The first-order (Born) approximation to the Lippmann-Schwinger equation for $\mathcal{R}^{\prime}_{\nu n}\left(\mathbf{r},\varepsilon\right)$ then reads

[TABLE]

where $G_{\nu n}\left(\mathbf{r},\mathbf{r}^{\prime},\varepsilon\right)$ is the single-site Green’s function, which for $r>S_{\nu}$ and $r^{\prime}<S_{\nu}$ can be expressed as

[TABLE]

Inserting Eq. (41) into Eq. (40) we obtain a form like Eq. (39) for $\mathcal{R}^{\prime}_{\nu n}\left(\mathbf{r},\varepsilon\right)$ from which the change of the $t$-matrix can be read off,

[TABLE]

Using the relationship Eq. (38), $\Delta\underline{m}_{\nu n}\left(\varepsilon\right)$ can finally be expressed as

[TABLE]

We have to emphasize that the regular solutions $\mathcal{Z}_{\nu n}(\mathbf{r},\varepsilon)$ and the SPO matrices $\underline{\tau}_{\nu n,\nu n}(\varepsilon)$ entering Eqs. (43) and (34), respectively, refer to the non-relativistic case, therefore, the transformations (14) and (20) apply. From these it follows that

[TABLE]

where on the right-hand side we explicitely marked that at the reference site $\nu 0$ the orientation of the magnetization is changed from $\mathbf{S}_{\nu 0}$ to $\mathbf{S}_{\nu n}$. Here we have to emphasize that due to the first-order perturbation treatment of the SOC there is no magnetic anisotropy included in the energy $\Delta\Omega$. Consequently, such a change of the spin vector at the reference site shouldn’t affect the contribution associated with site $\nu 0$. We tested numerically that, within a relative accuracy of about $10^{-5}$, the different sites in a sublattice add the same amount to the change of the energy (grand potential) due to the SOC, thus, it is sufficient to evaluate only the term for $\nu 0$ in Eq. (34) to obtain the energy correction per site.

We also note that it is possible to give an explicit expression for the SOC-induced energy per site, $\Delta\overline{\Omega}$. This is based on the reformulation of Eq. (44)

[TABLE]

where $\Delta\underline{\widehat{m}}_{\nu n}\left(\varepsilon\right)$ is defined by taking

[TABLE]

instead of $H^{\rm SOC}_{\nu}\left(\mathbf{r}\right)$ in Eq. (43) with the scattering solutions $\mathcal{Z}_{\nu 0}\left(\mathbf{r},\varepsilon\right)$. Averaging $\mathcal{R}_{n}^{-1}\mathbf{S}$ over $n$ for a wavelength of the spin spiral yields $\mathbf{n}\left(\mathbf{n}\mathbf{S}\right)$, where $\mathbf{n}$ refers to the axis of the rotation. Thus we arrive at

[TABLE]

with

[TABLE]

Despite of the above closed expression for $\Delta{\overline{\Omega}}$, we used the $\nu 0$ term in Eq. (34) to calculate the DM energy per site of the spin spiral, since the evaluation of the matrix in Eq. (43) is fairly simple when working in the relativistic $(j,\ell,m_{j})$ representation.

Because of the missing magnetic anisotropy $\Delta\Omega$ vanishes at $\mathbf{q}=0$ and it changes sign when reversing the wavevector, $\mathbf{q}\rightarrow-\mathbf{q}$, or, equivalently, the sense of rotation, $\mathbf{n}\rightarrow-\mathbf{n}$. $\Delta\Omega$ can therefore be identified with the DM energy of the spin spiral, $E_{\rm DM}$, which we add to the self-consistent non-relativistic total energy associated with the energy of isotropic spin-spin interactions. In addition, the expression Eq. (48) should be correlated with the fact that only the components of the DM vectors being parallel to ${\mathbf{n}}$ contribute to the DM energy of a spin spiral.

II.3 Spin-model parameters

For comparison with the results of the spin-spiral calculations as outlined above, we will use the extended classical Heisenberg model,

[TABLE]

where $\underline{J}_{ij}$ is the exchange coupling tensor Udvardi et al. (2003), and $\underline{K}_{i}$ is the on-site anisotropy matrix. The exchange coupling tensor can be decomposed into an isotropic, an antisymmetric and a traceless symmetric parts,

[TABLE]

where $\underline{I}$ stands for the $3\times 3$ unit matrix. The isotropic part $J_{ij}=\frac{1}{3}\textrm{Tr}\underline{J}_{ij}$ represents the Heisenberg couplings between the magnetic moments. According to the sign convention of Eq. (49), $J_{ij}>0$ and $J_{ij}<0$ indicate ferromagnetic (FM) and antiferromagnetic (AFM) couplings, respectively. The antisymmetric part of the exchange tensor can be identified with the DM vector $\mathbf{D}_{ij}$ as follows

[TABLE]

The traceless symmetric part of the exchange tensor is related to the two-site magnetic anisotropy, while the second term on the right-hand side of Eq. (49) to the one-site magnetic anisotropy. In all cases considered in this work the symmetry of the system implied uniaxial on-site anisotropy $-\sum_{i}K_{i}S_{i,z}^{2}$.

The spin-model parameters in Eq. (49) were determined by using the relativistic torque method as outlined in Ref. Udvardi et al., 2003. Note that, in order to obtain all the matrixelements of the $3\times 3$ exchange coupling matrices, ferromagnetic reference states oriented along different crystallographic directions should be used Udvardi et al. (2003). In order to produce coupling matrices that respect the symmetry of the lattice, for these orientations we considered the out-of plane ($z$) direction and, in case of $C_{4v}$ and $C_{3v}$ point-group symmetry, two and three independent in-plane directions, respectively.

In order to facilitate a comparison between the spin-spiral calculations and the spin model, Eq. (49), in case of fcc(111) and hcp(0001) surfaces we also determined the effective interaction parameters from the atomic interaction parameters derived for a large number of neighbors using RTM. According to Ref. Simon et al. (2018), the effective spin-model parameters for $C_{3v}$ point group symmetry are defined as

[TABLE]

known also as the spin stiffness Schweflinghaus et al. (2016) and spiralization Freimuth et al. (2014), respectively. The relationship between the micromagnetic and the effective parameters is given by

[TABLE]

where $V_{a}=\frac{\sqrt{3}}{2}a_{\rm 2D}^{2}t$ is the atomic volume, $a_{\rm 2D}$ and $t$ being the in-plane lattice constant and the film thickness, respectively. Beyond the effective and micromagnetic parameters we can also define effective nearest neighbor interactions, $J$ and $D$, which are related to the effective parameters as:

[TABLE]

In the present case we use the same sign convention as in Ref. Simon et al. (2018), namely, $D>0$ corresponds to the left-handed (counterclockwise) rotational direction, while $D<0$ to the right-handed (clockwise) rotational direction. The effective spin-model parameters can obviously determined from the long-wavelength (small-wavenumber) limit of the spin-wave spectrum: $J_{\rm eff}$ is related to the curvature of the non-relativistic energy dispersion, while $D_{\rm eff}$ to the slope of the SO induced contribution at $q=0$. From the calculated self-consistent total energy of the spin spirals we also determined the isotropic exchange parameters for several neighbors using least-squares fit, where the isotropic couplings from the RTM served as initial parameters for the fitting procedure.

In our applications we considered out-of-plane cycloidal spin spirals implying that all the spin vectors are in the ($\mathbf{q},\mathbf{z}$) plane. It is well-known that in this case the uniaxial anisotropy gives a contribution of $-K_{i}/2$ to the energy per site of a homogeneous spin spiral, while it adds $-K_{i}$ to the energy of the ferromagnetic state along the $z$ direction with respect to an in-plane direction of the magnetization. We then approximated the missing magnetic anisotropy part to the energy of the spin spirals by adding half of the magnetoscrystalline anisotropy energy (MAE) calculated for the ferromagnetic system. We derived the MAE in the spirit of the magnetic force theorem as a difference of the grand potential based on the band energy between the $x$ and $z$ direction of the magnetization, MAE = $\Omega_{x}-\Omega_{y}$ Jansen (1999); Daalderop et al. (1990); Szunyogh et al. (1995b).

In all calculations we used the atomic sphere approximation (ASA) for the effective potential with an angular momentum cut-off $\ell_{max}=2$ and the local spin-density approximation as parametrized by Vosko et al. Vosko et al. (1980). The energy integrals were performed along a semicircle contour in the upper complex energy semiplane. In case of the self-consistent calculations we used $3300-3600$ $k$-points in the full Brillouin-zone, while for the calculation of the DM energy more than $12000$ $k$-points were necessary to achieve a reliable accuracy.

III Results

III.1 Mn monolayer on W(001)

First we investigated the magnetic ground state of a Mn monolayer on W(001) in terms of spin-spiral calculations as implemented within the SKKR code for layered system. The model system consisted of four W layers, one Mn monolayer and three layers of empty spheres between a semi infinite W substrate and a semi infinite vacuum region. For the W and Mn layers epitaxial growth was assumed on a bcc(001) surface with the in-plane lattice constant of W(001), $a_{2\rm D}=3.165$ Å. The interlayer distance between the Mn layer and the topmost W layer was optimized by VASP calculations Kresse and Furthmüller (1996a, b); Hafner (2008). Relative to the interlayer distance in bulk W we found an inward relaxation of 12.6 % for the Mn monolayer.

The total energy of homogeneous flat spin spirals, Etot propagating along the (110) direction is shown in Fig. 2 (a) as a function of the wavenumber $q$. The magnitude of the Mn magnetic moment remained nearly constant with a value of 3.15 $\mu_{\rm B}$, while the induced moments of the W atoms changed in magnitude as a function of the wavenumber of the spin spiral. Note that at the $\overline{\rm M}$ point of the Brillouin zone, i.e. for a row-wise AFM configuration of the Mn moments, the induced moments of W vanishes. Apparently, as the state with lowest energy we obtained a spin spiral with a wavenumber of $q=0.3$ Å*-1*, thus, with a wavelength of $\lambda=2.1$ nm. This is in relatively good agreement with the corresponding wavelength of $\lambda=3.1$ nm reported in Ref. Ferriani et al., 2008 for the same system in terms of FLAPW spin-spiral calculations. The energy gain of this spin spiral is about 3 meV/Mn atom with respect to the ferromagnetic state.

We fitted the calculated total energies of the spin spirals to an isotropic Heisenberg model containing the first five nearest neighbor (NN) interactions. We found that the fitted curve matched the calculated points very accurately along the whole $\overline{\Gamma\rm M}$ line. The fitted interactions are shown in Table 1, together with those calculated by using the RTM. It should be mentioned that along this propagation direction of the spin spiral, the second and third NN couplings ($J_{2}$ and $J_{3}$) can not be obtained independently from a fitting up to the fifth NN interactions, since only $J_{2}+2J_{3}$ can uniquely be determined. We therefore present this value in Table 1. Obviously, the dominant coupling is the strong ferromagnetic nearest-neighbor interaction, while the farther interactions are antiferromagnetic and considerably smaller in size than the NN interaction. The interactions obtained from the two methods compare well to each other, except for the fifth NN coupling which is ferromagnetic for the RTM and antiferromagnetic in case of the spin-spiral fit. The isotropic interactions fitted to spin-spiral energies calculated by the FLAPW method Nandy et al. (2016) are also listed in Table 1. Note that due to the definition of the spin Hamiltonian in Eq. (49), the parameters given in Ref. Nandy et al., 2016 were multiplied by a factor of two. While there is an overall good agreement between the parameters obtained from the two spin-spiral calculations, $J_{1}$ and $J_{2}+2J_{3}$ from the FLAPW method are clearly smaller in magnitude than our calculated values. This can partly be attributed to the considerably smaller inward relaxation of the Mn monolayer used in the FLAPW calculations (4.7 %) Ferriani et al. (2005) as compared to our calculations (12.6 %).

The DM energy, EDM calculated from Eq. (34) adds an antisymmetric term to the spin-spiral dispersion as shown in Fig. 2(b). As a result we find the energy minimum for a right-handed cycloidal spin-spiral with a period of $\lambda=1.92$ nm which agrees well both with the experimental value of $\lambda=2.1$ nm and with the value of $\lambda=2.3$ nm obtained from FLAPW spin-spiral calculations including SOC Ferriani et al. (2008). Reassuringly, the W atoms exhibiting large SOC and considerable spin-polarization have an overwhelming contribution to the DM energy, while the contribution of the Mn atoms is negligible, see Fig. 2(b). For small $q$ values, the DM energy of the spin spirals is proportional with $q$ and from the slope of the curve an effective nearest neighbor DM interaction of $11.6$ meV can be fitted, which is in good agreement with the corresponding value of $9.2$ meV reported in Ref. Ferriani et al. (2008) using our convention for the exchange interactions, see above. The spin-spiral energies based on the spin model Eq. (49) with the tensorial interactions from the RTM but excluding the anisotropy terms are also presented in Fig. 2(a). In accordance with the spin-spiral calculations this also prefers a right-handed spin-spiral, nevertheless with a somewhat smaller period of $1.6$ nm.

We can account for the magnetic anisotropy as explained in Section II. C. Within the framework of magnetic force theorem, we found an out-of plane anisotropy with value of $K=4.15$ meV. This implies that the energy of the spin spirals should be shifted upwards with respect to the energy of the FM state by $K/2=2.08$ meV, close to the value of $1.8$ meV reported in Ref. Ferriani et al., 2008. Since this energy shift is much smaller than the energy gain of the spin spiral due to isotropic and DM interactions ($>$15 meV, see Fig. 2(a)), this spin-spiral state remains lower in energy than the FM state.

III.2 Co monolayer on Pt(111)

Next we performed spin-spiral calculations for a Co monolayer on Pt(111). The self-consistently treated layer structure considered in the SKKR method consisted of five Pt atomic layers, one Co monolayer and three layers of empty spheres between a semi-infinite Pt substrate and a semi-infinite vacuum region. For modelling the geometry of the system we used the in-plane lattice constant of Pt(111), $a_{2\rm D}=2.774$ Å, fcc growth was assumed for the Pt layers and hcp stacking was used for the Co monolayer. The distance between the atomic layers were optimized in terms of VASP calculations. Relative to the interlayer distance in bulk Pt, for the Co monolayer we found an inward relaxation of $11$% .

Considering homogeneous flat spin spirals rotating in the $xz$ plane (the axis $x$ denoting an in-plane nearest neighbor direction), we calculated the total energy for ${\bf q}$-vectors along the $\overline{\Gamma\rm K}$ direction in the 2D Brillouin zone. The magnetic moment of the Co atoms proved to be fairly independent of $q$ with a value of $2.13\,\mu_{\rm B}$ and the Pt layers showed to have induced moments changing according to the spin-spiral wave vector. In Fig. 3(a) the calculated non-relativistic and relativistic dispersion of the spin spirals are shown, the latter one obtained by adding the DM term Eq. (34) to the non-relativistic spin spiral energy. The spin-spiral energy calculated from the spin model Eq. (49) with only isotropic and DM interactions determined by RTM is also presented in Fig. 3(a).

Within the investigated range of $q$, the non-relativistic spin-spiral dispersion turned out to be fairly parabolic, thus a fit to a Heisenberg spin model allowed us to determine the first two nearest-neighbor interactions between the Co moments, $J_{1}$ and $J_{2}$. We also calculated the effective and micromagnetic isotropic parameters defined in Eqs. (52), (54) and (55), $J_{\rm eff}$, $\mathcal{J}$ and $J$, respectively. These parameters are listed in Table 2. The interactions derived from the RTM were previously reported in Ref. Vida et al. (2016) for hcp stacking of the Co monolayer. For comparison, the corresponding values for $J_{1}$ and $J_{2}$, as well as for the effective and micromagnetic parameters are also presented in Table 2. We find that the self-consistent spin-spiral calculations give a NN isotropic coupling and effective parameters by about 30 % larger than the respective parameters from the torque method. This difference can be attributed to the ferromagnetic coupling between the Co moments and the induced moments of Pt that are included in the self-consistent spin-spiral calculations, but not taken into account in the Co-Co interactions obtained from the RTM.

By using Eq. (34) we calculated the DM contribution to the spin-spiral energy and presented it in Fig. 3(b). Due to the large SOC of Pt, the DM energy mainly originates from the topmost Pt layer, while the Co layer has a much smaller contribution. From the slope of the EDM curve at $q=0$ we determined the spiralization $D_{\rm eff}$, the microscopic DM parameter $\mathcal{D}$ from Eq. (54) and the effective nearest neighbor DM coupling $D$ from Eq. (55). We also derived these parameters based on the previously reported in-plane DM interactions calculated in terms of the RTM Vida et al. (2016) and summarized them in Table 3. Most likely again due to the strong interaction between the Co and Pt moments, the effective DM parameters turned out to be by about 50 % larger in case of the spin-spiral calculations as compared to the torque method and in both calculations the rotational direction of the DM vectors is left-handed (counterclockwise) in agreement with other theoretical results Yang et al. (2015); Vida et al. (2016); Freimuth et al. (2014); Dupé et al. (2014). Due to the large ferromagnetic isotropic coupling, the energy gain from the formation of spin-spiral states is negligible and the easy-axis MAE of $0.57$ meV (see in Ref. Vida et al. (2016) ) clearly stabilizes an out-of-plane ferromagnetic ground state.

III.3 Co/Pt(111) capped by a Ru overlayer

In Ref. Simon et al., 2018 we showed that the DM interaction is decreased when Co/Pt(111) is capped by a $5d$ monolayer. In this Section we examine the effect of the $4d$ Ru overlayer on top of Co/Pt(111) which was recently claimed to induce enhanced interfacial DMI Karayev et al. (2019). The thin film system was modelled by four monolayers of Pt, one Co monolayer, one Ru monolayer and three monolayers of empty sheres between a semi-infinite Pt substrate and a semi-infinite vacuum region. We used the in-plane lattice constant of the Pt(111) surface, $a_{2\rm D}=2.774$ Å. From VASP calculations we have found an inward relaxation for the Co monolayer and for the Ru overlayer of $-8$% and $-15$% relative to the interlayer distance in bulk Pt, respectively. We assumed fcc growth for the Pt layers fcc, while hcp stacking for the Co and Ru monolayers.

Figure 4(a) shows the spin-spiral dispersion as a function of the spin-spiral wave vector along the $\overline{\Gamma\rm K}$ direction for the Ru/Co/Pt(111) system. The magnitude of the Co magnetic moment of $1.90\mu_{B}$ was fairly the same for any $q$, while the Ru and Pt layers exhibited induced moments decreasing in size with increasing $q$. Similar to the uncapped Co/Pt(111) monolayer, the non-relativistic spin-spiral dispersion indicates the preference of a ferromagnetic order. From a parabolic fit of the non-relativistic dispersion we obtained an effective isotropic nearest-neighbor interaction of $J_{\rm SS}=44.43$ meV, implying that the Ru overlayer decreased the FM coupling between the Co atoms as compared to the uncapped Co/Pt(111) system (compare with the value of 65.30 meV in Table  2). From the RTM calculations we determined a value of $J_{\rm RTM}=35.72$ meV, which also reflects the reduced isotropic Co-Co interaction.

The calculated Dzyaloshinskii-Moriya energy and its layer-wise composition is presented in Fig. 4(b). Apparently, the Ru overlayer significantly contributes to EDM for larger $q$ values and it has the same sign as the contribution of the Pt layers, which is slightly enhanced as compared to the case of uncapped Co/Pt(111). In spite that the Co monolayer has now a remarkably reduced contribution, the Ru overlayer overall increase the DM energy. This is reflected also in the effective NN DM interaction of $D_{\rm SS}=3.87$ meV which is by about 1 meV larger than for the uncapped system, see in Table 3. From the RTM method we also found a larger effective DM parameter, $D_{\rm RTM}=3.14$ meV as compared to $D_{\rm RTM}=1.84$ meV for Co/Pt(111). The rotational direction of the in-plane DM vectors is left-handed as in the Co/Pt(111) system. The Ru overlayer also drastically modified the magnetic anisotropy, because the preferred magnetization direction became in-plane, while for the uncapped system we obtained an out-of plane magnetization. The obtained in-plane anisotropy energy, $E_{x}-E_{z}=-0.56$ meV, is larger in magnitude than the very small energy gain from the formation of a spin-spiral state, indicating that the magnetic ground state is in-plane ferromagnetic. Note that in Ref. Karayev et al., 2019 a perpendicular magnetic anisotropy has been observed for the Co/Pt/Ru multilayers, but this does not contradict our present result, which refers to an overlayer system with a free surface.

III.4 Pt/(FeCo)/Ir multilayer systems

Finally, we investigated periodic Pt/(FeCo)/Ir superlattices to model the [Ir(10 Å)/Fe(0-6 Å)/Co(4-6 Å)/Pt(10 Å)]20 multilayers in which room-temperature magnetic skyrmions have been found recently (Soumyanarayanan et al., 2017). In the calculations we considered an Fe and a Co monolayer in both possible sequences, sandwiched between an Ir and a Pt bilayer, and repeated this unit periodically along the direction normal to the planes. As what follows we will label these multilayers by Pt/Co/Fe/Ir and Pt/Fe/Co/Ir. For the hexagonal layers we used the in-plane lattice constant of Ir(111), $a_{2\rm D}=2.714$ Å. The Pt, Ir and Fe monolayers were stacked in fcc geometry. For the Co monolayer both hcp and fcc stackings were considered, however, the self-consistent spin-spiral calculations were performed only for the fcc stacking of the Co layer. The interlayer distances were optimized from VASP calculations, where we found that the interlayer distances were independent on the stacking of the Co layer. The calculated interlayer distances are summarized in Table 4. These interlayer distances were used in the self-consistent spin-spiral SKKR and RTM calculations and the Wigner-Seitz radii of the atomic spheres in the Co, Fe, Pt and Ir layers were modified according to the relaxations.

The calculated spin-spiral dispersions for the Pt/Fe/Co/Ir and Pt/Co/Fe/Ir multilayer structures are shown in Fig. 5 for a propagation direction along the $\overline{\Gamma\rm K}$ line in the 2D Brillouin zone. The fairly stable moments of the Fe and Co were m_{\rm Fe}=2.84\,$$\mu_{\rm B} and m_{\rm Co}=1.59\,$$\mu_{\rm B} for Pt/Fe/Co/Ir, and m_{\rm Co}=1.83\,$$\mu_{\rm B} and m_{\rm Fe}=2.43\,$$\mu_{\rm B} for Pt/Co/Fe/Ir. The non-relativistic total energy dispersion curve shows a faster increase for Pt/Fe/Co/Pt than for Pt/Co/Fe/Pt, while in both cases a shallow minimum can be found at small $q$, which indicates the appearance of frustrated isotropic interactions.

We calculated the isotropic exchange interactions for the Fe-Fe, Fe-Co and Co-Co pairs using the relativistic torque method and plotted them in Fig. 6 as a function of the interatomic distances in case of an fcc-stacking of the Co layer. Note that interactions are presented only within a bilayer of CoFe or FeCo, since the interactions between the bilayers are negligible. In both multilayer structures the ferromagnetic NN Fe-Fe, Fe-Co and Co-Co interactions are dominating. The NN Fe-Co interaction in Pt/Fe/Co/Ir is clearly enhanced as compared to Pt/Co/Fe/Ir, being the main reason for the steeper spin-spiral energy dispersion for Pt/Fe/Co/Ir seen in Fig. 5. While the Fe-Co interactions remain ferromagnetic for larger distances, the second and third NN Fe-Fe interactions, as well as the third NN Co-Co interactions are antiferromagnetic, which gives rise to frustration and to the stabilization of long wavelength spin spirals as seen in Fig. 5.

The DM energy as calculated from Eq. (34) and its resolution into layer-wise contributions is shown in Fig. 7 for an fcc stacking of the Co monolayer. In case of Pt/Fe/Co/Ir the Pt layers have a dominating contribution to EDM, while the contribution of the Ir layers is about ten times smaller and, surprisingly, it is similar in magnitude as the contribution of the Fe layer. Moreover, the contribution of Pt is different in sign as those of Ir and Fe, and EDM prefers a right-handed spin-spiral state. In case of Pt/Co/Fe/Ir structure the contribution of the Pt layers is practically unchanged, but the contribution of the Ir layers enhances in size by a factor of about five and becomes dominant in EDM. The contribution of Fe considerably decreases and also reverses sign. Note that in both cases the DM energy related to Co is negligible. As a result, EDM for Pt/Co/Fe/Ir reverses sign as compared with Pt/Fe/Co/Ir, thus it favors left-handed spin spirals.

From the slope of EDM at $q=0$ we obtained the spiralization $D_{\rm eff}$ and the effective nearest-neighbor DM interaction and presented them in Table 5 with the corresponding parameters calculated in terms of the RTM for both the fcc and hcp stacking of the Co monolayer. In case of fcc stacking of the Co monolayer, by using both computational methods the sign of the effective DM parameters is different for the two multilayer systems indicating that the sense of rotation of the DM vectors depends on the sequence of the magnetic layers. On the contrary, for hcp stacking of the Co monolayer, the sign of the DMI does not change and similar to the case of together fcc stacking of Co, the effective DM coupling is three times larger for Pt/Co/Fe/Ir than for Pt/Fe/Co/Ir.

By using the magnetic force theorem we also determined the magnetic anisotropy energy of the multilayer systems for both kinds of Co stacking. We have found that the Pt/Fe/Co/Ir multilayer shows perpendicular magnetic anisotropy in case of both the fcc and the hcp stacking of Co, with $\Delta E=E_{x}-E_{z}=0.76$ meV and $0.36$ meV, respectively. In case of the Pt/Co/Fe/Ir multilayer with Co-fcc stacking we obtained an in-plane anisotropy with $\Delta E=-0.64$ meV, while for Co-hcp stacking an easy-axis anisotropy with $\Delta E=0.59$ meV. It is known, however, that the magnetic anisotropy of multilayers is very sensitive to the growth conditions den Broeder et al. (1991); Daalderop et al. (1990); Wu et al. (2012) and, similarly as in Ref. Soumyanarayanan et al., 2017, our calculated values of the MAE are much larger than the experimental values. Due to the frustrated isotropic couplings between the Fe and Co atoms and due to the large DM interaction, skyrmionic states can then likely be stabilized in Pt/Co/Fe/Ir multilayers.

IV Summary and conclusions

In summary, by using the multiple scattering Green’s function technique we presented a theoretical approach to calculate the electronic structure of layered systems with spiral magnetic structure. For layered systems the propagation direction of the spin-spiral is restricted to the plane of the layers, while an arbitrary rotational direction can be chosen. The non-relativistic allows for self-consistent calculations, from which the total energy of the system can be obtained as a function of the spin-wave vector. We employed a first-order perturbation technique to include the effect of spin-orbit coupling in the calculations. A particular advantage of this approach is that the energy related to the Dzyaloshinskii-Moriya interactions can be resolved into layer-wise contributions.

We performed ab initio calculations for ultrathin films and multilayers and demonstrated that the newly developed method gives an accurate access to the magnetism of these systems. We found that the magnetic ground state of a Mn monolayer on W(001) is a right-handed spin-spiral in good agreement with the experiment and calculations in Ref. Ferriani et al., 2008. From the spin-spiral dispersion we derived spin-model parameters, which compared well with the couplings calculated by the relativistic torque method. We investigated the spin-spiral states of a Co monolayer on Pt(111) with hcp stacking of the Co monolayer and concluded that the ground state of the system is ferromagnetic, similar to previous calculations Dupé et al. (2014); Zimmermann et al. (2019). The nearest-neighbor isotropic exchange interaction determined from the non-relativistic spin-spiral dispersion was found significantly larger than the corresponding interaction from the torque method which can be attributed to the effect of the induced moments of the Pt atoms, included inherently in the spin-spiral calculations. For the Co monolayer on Pt(111) capped by a Ru overlayer we showed that the isotropic coupling between the Co atoms is reduced, while the interfacial DMI was increased. These results correlate well with recent experiments on Pt/Co/Ru superlattices Karayev et al. (2019). We also investigated Pt/Fe/Co/Ir and Pt/Co/Fe/Ir multilayers and found that the non-relativistic dispersion implies the appearence of spin-spiral states with large wavelength due to the frustrated couplings between the Fe and Co atoms. Remarkably, the rotational sense of the DMI was opposite for the two multilayers, which could be attributed to the largely enhanced contribution of the Ir layer in the Pt/Co/Fe/Ir multilayer. Moreover, we found that the effective DMI is three times larger for Pt/Co/Fe/Ir than in case of Pt/Fe/Co/Ir independently on the stacking of the Co layer.

Our results obtained from the presented spin-spiral approach provide thus a theoretical support to the fine-tuning of the magnetic and non-magnetic layers in multilayer structures with the purpose of manipulating the interfacial DMI and designing new building elements for spintronics applications Soumyanarayanan et al. (2017). An obvious possibility to proceed on this way is to consider disordered alloys in these structures that is easily feasible within the Green’s function technique in terms of the coherent potential approximation Mankovsky et al. (2011).

Acknowledgements.

The authors are grateful for illuminating discussions with Manuel Dos Santos Dias and Laszlo Udvardi, and acknowledge NIIF and PRACE for awarding access to the resources based in Hungary at Debrecen and Szeged, and in Germany at Jülich. This work was supported by the National Research, Development and Innovation Office of Hungary under projects No. PD120917 and K115575, as well as by the BME Nanotechnology FIKP grant of EMMI (BME FIKP-NAT).