Micromagnetic simulations of spin-torque driven magnetisation dynamics with spatially resolved spin transport and magnetisation texture

TL;DR

This paper introduces a fast simulation method coupling spin transport and micromagnetics to accurately model spin-torque driven magnetisation dynamics in nano-structures, aligning well with experimental results.

Contribution

A novel coupling approach between spin transport and micromagnetics for realistic, spatially resolved simulations of spin-torque phenomena in nano-pillar devices.

Findings

Accurately reproduces spin-wave mode excitation and threshold currents.

Matches experimental observations of magnetoresistance and switching.

Provides insights into spin-caloritronics device behavior.

Abstract

We present a simple and fast method to simulate spin-torque driven magnetisation dynamics in nano-pillar spin-valve structures. The approach is based on the coupling between a spin transport code based on random matrix theory and a micromagnetics finite-elements software. In this way the spatial dependence of both spin transport and magnetisation dynamics is properly taken into account. Our results are compared with experiments. The excitation of the spin-wave modes, in- cluding the threshold current for steady state magnetisation precession and the nonlinear frequency shift of the modes are reproduced correctly. The giant magneto resistance effect and the magnetisa- tion switching also agree with experiment. The similarities with recently described spin-caloritronics devices are also discussed.

Figures

Tables

Table 1. Table 1: Parameters of thin (a) and thick (b) layer used in micromagnetics simulations.

$4\pi M_a$ (G)	${\alpha }_a$	$4\pi M_b$ (G)	${\alpha }_b$	$\gamma$ (rad.s${}_{}{}^{-}1$.G)
$8.2\times {10}^3$	$1.5\times {10}^{-2}$	$9.6\times {10}^3$	$9\times {10}^{-3}$	$1.87\times {10}^7$

Table 2. Table 2: Comparative table of the resonance fields of the ℓ = 0 ℓ 0 \ell=0 SW modes. Top are the peak locations measured experimentally in Ref.[ Naletov et al., 2011 ] Bottom are the eigen-frequencies extracted from the simulation with applied field H ext = 1 subscript 𝐻 ext 1 H_{\rm{ext}}=1  T along z 𝑧 z .

Exp. $f$ (GHz)	6.08	8.95	9.82	11.98
Sim. $f$ (GHz)	6.08	8.94	9.83	12.00
SW modes	$B_{00}$	$B_{01}$	$A_{00}$	$A_{01}$

Table 3. Table 4: Bulk VF parameters for the materials which constitute the multilayer. The VF bulk resistivity ρ ∗ superscript 𝜌 \rho^{*} is expressed in units 10 − 9 ​ Ω ​ m superscript 10 9 Ω 𝑚 10^{-9}\Omega m and the spin-flip l sf subscript 𝑙 sf l_{\mbox{sf}} length in n ​ m 𝑛 𝑚 nm .

Material	$\rho \ast$	$\beta$	$1/l_{\text{sf}}$
Cu	5	0	0.002
Au	20	0	0.033
Py	291	0.76	0.182

Equations

$$⟨M_z⟩=\frac{1}{V}{\int }_V{M}_z(\mathbf{r})d^3\mathbf{r}.$$

$$\frac{1}{\gamma }\dot{{\mathbf{M}}^j}={\mathbf{M}}^j\times {\mathbf{H}}_{\text{eff}}^j+\alpha {\mathbf{M}}^j\times {\dot{\mathbf{M}}}^j+{\mathbf{M}}^j\times {\mathbf{S}}^j.$$

$${\mathbf{d}}^j={\alpha }^j{\dot{\mathbf{M}}}^j$$

$${\mathbf{S}}^j=\frac{I_{dc}}{2\pi \lambda }\left[{\mathbf{m}}^j\times {\mathbf{m}}^{j\prime }\right].$$

$${\lambda }^j=\frac{2e{M}_s^{j}V}{\eta h}$$

$$d=\alpha (1-I_{\text{dc}}/I_{\text{th}}),$$

$$c_j=\frac{M_x^j+i{M}_y^j}{\sqrt{2M_s^j(M_s^j+M_z^j})},$$

$${\dot{c}}_a$$

$${\dot{c}}_b$$

$${\dot{p}}_a=-2({\Gamma }_{a-}-{\Gamma }_{a+})p_a+j_{ab}^p.$$

$$j_{ab}^p=2h_{ab}\sqrt{p_a{p}_b}\sin [{\varphi }_a(t)-{\varphi }_b(t)+\beta ].$$

$$c_{\ell ,m}(r,\varphi ,t)=J_{\ell }(k_{\ell ,m}r)e^{+i\ell \varphi }{e}^{-i{\omega }_{\ell m}t},$$

$${\dot{\mathbf{m}}}_i^j=-\gamma ({\mathbf{m}}_i\times {\mathbf{H}}_{\mathrm{eff}i}^j)+\frac{\alpha }{M_s}({\mathbf{m}}_i^j\times {\dot{\mathbf{m}}}_i^j)$$

$${\mathbf{m}}^j=\frac{1}{N}\sum _{i=1}^N{\mathbf{m}}_i^j.$$

$${\bm{\psi}}_{i\pm}=\left(\begin{array}[]{c}{\bm{\psi}}_{i\pm\uparrow}\\ {\bm{\psi}}_{i\pm\downarrow}\end{array}\right).$$

$$\left(\begin{array}[]{c}{\bm{\psi}}_{0-}\\ {\bm{\psi}}_{1+}\end{array}\right)=S\left(\begin{array}[]{c}{\bm{\psi}}_{0+}\\ {\bm{\psi}}_{1-}\end{array}\right),$$

$$S=\left(\begin{array}[]{cc}r^{\prime}&t\\ t^{\prime}&r\end{array}\right).$$

$$t=\left(\begin{array}[]{cc}t_{\uparrow\uparrow}&t_{\uparrow\downarrow}\\ t_{\downarrow\uparrow}&t_{\downarrow\downarrow}\end{array}\right)$$

$$G=\frac{e^2}{h}\text{mbox}Tr\left[t^{†}t\right],$$

$$\frac{\partial {\vec{J}}_0}{\partial \mu }=\frac{1}{4\pi }\text{mbox}Tr\left[t\vec{\sigma }{t}^{†}\right],$$

$${\hat{t}}_{\sigma \eta ,{\sigma }^{\prime }{\eta }^{\prime }}=\frac{1}{N_{ch}}\text{mbox}{Tr}_{N_{ch}}[t_{\sigma {\sigma }^{\prime }}{t}_{\eta {\eta }^{\prime }}^{†}].$$

$$\hat{t}=\frac{1}{N_{ch}}\mbox{Tr}_{N_{ch}}\left(\begin{array}[]{cccc}t_{\uparrow\uparrow}t_{\uparrow\uparrow}^{\dagger}&t_{\uparrow\uparrow}t_{\uparrow\downarrow}^{\dagger}&t_{\uparrow\downarrow}t_{\uparrow\uparrow}^{\dagger}&t_{\uparrow\downarrow}t_{\uparrow\downarrow}^{\dagger}\\ t_{\uparrow\uparrow}t_{\downarrow\uparrow}^{\dagger}&t_{\uparrow\uparrow}t_{\downarrow\downarrow}^{\dagger}&t_{\uparrow\downarrow}t_{\downarrow\uparrow}^{\dagger}&t_{\uparrow\downarrow}t_{\downarrow\downarrow}^{\dagger}\\ t_{\downarrow\uparrow}t_{\uparrow\uparrow}^{\dagger}&t_{\downarrow\uparrow}t_{\uparrow\downarrow}^{\dagger}&t_{\downarrow\downarrow}t_{\uparrow\uparrow}^{\dagger}&t_{\downarrow\downarrow}t_{\uparrow\downarrow}^{\dagger}\\ t_{\downarrow\uparrow}t_{\downarrow\uparrow}^{\dagger}&t_{\downarrow\uparrow}t_{\downarrow\downarrow}^{\dagger}&t_{\downarrow\downarrow}t_{\downarrow\uparrow}^{\dagger}&t_{\downarrow\downarrow}t_{\downarrow\downarrow}^{\dagger}\end{array}\right),$$

$$\hat{t}=\left(\begin{array}[]{cccc}T_{\uparrow\uparrow}&0&0&T_{\uparrow\downarrow}\\ 0&T_{mx}&0&0\\ 0&0&T_{mx}^{*}&0\\ T_{\downarrow\uparrow}&0&0&T_{\downarrow\downarrow}\end{array}\right),$$

$$\hat{S}=\left(\begin{array}[]{cc}\hat{r}^{\prime}&\hat{t}\\ \hat{t}^{\prime}&\hat{r}\end{array}\right).$$

$$R_{\theta ,\vec{n}}=\exp (-i\vec{\sigma }\cdot \vec{n}\theta /2)$$

$$\hat{\tilde{S}}={\hat{R}}_{\theta ,\vec{n}}\hat{S}{\hat{R}}_{\theta ,\vec{n}}^{†},$$

$$G=\frac{1}{{\mathcal{R}}_{sh}}\left(T_{\uparrow \uparrow }+T_{\uparrow \downarrow }+T_{\downarrow \uparrow }+T_{\downarrow \downarrow }\right)$$

$${\bm{P}}_{i\pm}=\left(\begin{array}[]{l}{P}_{i\pm,\uparrow}\\ {P}_{i\pm,mx}\\ {P}^{*}_{i\pm,mx}\\ {P}_{i\pm,\downarrow}\\ \end{array}\right).$$

$$\displaystyle\left(\begin{array}[]{c}{\bm{P}}_{0-}\\ {\bm{P}}_{1+}\end{array}\right)$$

$$\displaystyle\left(\begin{array}[]{c}{\bm{P}}_{1-}\\ {\bm{P}}_{2+}\end{array}\right)$$

$$\left(\begin{array}[]{c}{\bm{P}}_{0-}\\ {\bm{P}}_{2+}\end{array}\right)=\hat{S}_{a+b}\left(\begin{array}[]{c}{\bm{P}}_{0+}\\ {\bm{P}}_{2-}\end{array}\right).$$

Full text

I Introduction

The orientation of the magnetization in a magnetic film can be influenced using a spin-polarized current. Consequently, a direct current can transfer spin angular momentum between magnetic layers, separated by either a normal metal or a thysichin insulating layer. This effect is called spin transfer torque (STT) and was first discussed in the 1970s in the context of moving magnetic domain walls Berger (1978) and fully understood in the 1990s Berger (1996); Slonczewski (1996). STT has been of profound importance for the development of spintronic devices such as read-heads based on the giant magnetoresitive (GMR) effect Baibich et al. (1988); Binasch et al. (1989), the spin-transfer torque magnetic random-access memory (STT-MRAM) Diao et al. (2007) and spin-torque nano-oscillators (STNO)Katine et al. (2000).

Until very recently C. et al. (2014) the approaches to theoretically describe the magnetization dynamics induced by a spin torque usually greatly simplified or neglected the description of either the spatial inhomogeneity of the spin torque, or the three-dimensional magnetization texture K. J. Lee (2004); Xiao et al. (2005); Berkov and Miltat (2008).

In the present work, we go beyond such approaches by coupling a finite element micromagnetic method Fischbacher et al. (2007) to a numerical solver for spin transport, based on continuous random matrix theory (CRMT) Rychkov et al. (2009); Borlenghi et al. (2011). In this way the effect of spin torque is described, with the transport and magnetic degrees of freedom treated on an equal footing. The spatial inhomogeneity of both spin transport and magnetization dynamics is thus explicitly included. In our implementation, CRMT is parametrised by the same set of experimentally accessible parameters as in Valet-Fert theory Valet and Fert (1993), so that our numerical simulations contain no free parameters.

We demonstrate the capabilities of our computational method by addressing theoretically the effect of spin torque on the magnetization dynamics in the perpendicularly magnetized circular spin-valve nanopillar experimentally investigated in Ref. [Naletov et al., 2011]. This configuration is obtained by saturating the device with a large applied field perpendicular to the layers. The device setup is kept very simple and as such serves as a prototype for spin-valve structures, which find application also in the emerging field of spin-caloritronics Borlenghi et al. (2014a, b, 2015a, 2015b). In particular, our configuration corresponds to a circular precession of the magnetisation, allowing for a precise identification of the spin wave (SW) modes.This is crucial if one wants to couple the system to an external rf signal, since only signal with the same symmetries of the SW modes can excite the magnetisation. Breaking the axial symmetry results in a more complicated configuration where modes with different symmetries mix up Naletov et al. (2011).

Moreover, the system can be described using the language of coupled oscillators Slavin and Tiberkevich (2009) and in particular of the discrete nonlinear Schrödinger equation Borlenghi et al. (2014a) (DNLS). Indeed, this setup corresponds to the simplest realisation of the DNLS, containing only two elements. Here spin transfer torque physically corresponds to a magnon chemical potential Borlenghi et al. (2014b) that controls the propagation of the energy and magnetisation currents between the two layers. The DNLS appears in many branches of Physics, including Bose-Einstein condensates, photonics waveguides and photosynthetic reactions. Understanding the dynamics in our setup can therefore shed light on a very general oscillator model.

The remainder of this paper is organized as follows: in Sec. II we describe the geometry of the nanopillar and also briefly review the classification of SW modes in a perpendicularly magnetised nanopillar. Those modes are then identified by means of micromagnetic simulations at zero current (thus without spin torque). This study was previously performed in Ref. [Naletov et al., 2011], to which we refer for a thorough discussion. Sec. III reviews key concepts of scattering approach and CRMT to describe transport in magnetic multilayers, and consitutes an extension of the material presented in Refs. [Waintal et al., 2000] and [Rychkov et al., 2009]. In Sec. IV we describe how to couple the CRMT transport code to Nmag Fischbacher et al. (2007), in order to simulate the effect of spatially inhomogeneous spin transfer torque and magnetisation dynamics on the same footing. Sec. V contains our micromagnetic simulations of current driven spin dynamics. Here we identify the SW modes excited by spin transfer torque, the critical current for auto oscillations, and the frequency shift beyond the critical threshold Slavin and Tiberkevich (2009). Using CRMT we provide a precise characterisation of the magnetoresistance. Our simulations are then compared with the experimental results found in Ref. Naletov et al., 2011. Finally, in the conclusion we summarise the main results of this paper and point out further possible developments.

II Physical system and model

II.1 Spin valve structure

The nanopillar studied here is displayed in Fig.1(a). It consists of a trilayer structure made of two Permalloy ($\rm{Ni_{80}Fe_{20}}$ alloy) disks Pya and Pyb separated by a 10 nm Cu spacer. The disks have diameter $d=200$ nm and thicknesses $t_{a}=4$ nm (upper disk) and $t_{b}=15$ nm (lower disk). The upper disk is connected to a 25 nm Au contact and the lower disk to a 60 nm Cu contact. In the experimentNaletov et al. (2011), the sample was mounted inside a magnetic resonance force microscope (MRFM) and the whole apparatus was placed inside a vacuum chamber operated at room temperature. The external magnetic field $\bm{H}_{\rm{ext}}$, was oriented along the pillar axis $\bm{z}$, which corresponds to the precession axis of the magnetisation.

The MRFM consists of an ultra-soft cantilever with a 100 nm diameter magnetic sphere glued to its tip. The sphere is positioned precisely above the center of the nano-pillar, so as to retain the axial symmetry. The mechanical-FMR spectroscopy consists in recording, by optical means, the vibration amplitude of the cantilever as a function of the bias out-of-plane magnetic field in the presence of a RF field with fixed frequency excitation Naletov et al. (2011). The dynamics can be excited also by injecting a dc current along $\bm{z}$, which excites the dynamics in the thin layer due to spin transfer torque. In the steady state, the combined torques exerted by the RF excitations and spin transfer compensate the damping, so that the local magnetization vector precesses regularly at the Larmor frequency along a circular orbit, see Fig.1a).

The instantaneous magnetization $\bm{M}$ can be decomposed into a large static component $\bm{M}_{\parallel}$ and a small oscillating component $\bm{M}_{\perp}$, where $\bm{M}_{\parallel}$ is parallel to the local precession axis (i.e., the direction of the external magnetic field), and $\bm{M}_{\perp}$ is perpendicular to that axis.

An essential feature of our system is that different SW modes have different spatial distribution of the phase of $\bm{M}_{\perp}$ inside the magnetic disks. Those modes can be excited only by RF fields with the same rotational symmetries, giving selection rules for the excitation of the SW modes Naletov et al. (2011).

The dipolar force acting on the cantilever is proportional to the spatially averaged value of the longitudinal (static) component of the magnetization inside the whole nanopillar,

[TABLE]

The latter is not subject to any selection rule, so that the the mechanical-FMR setup detects all possible SW modes. Experimentally, it has been observed that the presence of the cantilever introduces a shift of $+0.57$ GHz in the SW spectrum. This has be taken into account in our simulations.

The dynamics of the magnetization can be excited also by means of an rf current flowing along the axis of the pillar, which generates an rf orthoradial Oersted field, and by mean of an homogeneous in-plane RF magnetic field.

We remark that, although in experiments the dynamics driven by a combination of RF fields and STT, in our simulations the dynamics is excited by STT only. The effect of the RF field is modelled by setting different initial condition for the magnetisation.

In our circuit, a positive current corresponds to a flow of electrons from the bottom Pyb thick layer to the top Pya thin layer, and stabilizes the parallel configuration due to the spin transfer effect. Vice-versa, a negative current stabilizes the thick layer and destabilizes the thin one. At low current, thick and thin layer are thus the fix and free layer correspondingly. However, at high enough current, both layers precess due to the repulsive dipolar interaction, as will be discussed later.

II.2 Magnetization dynamics

In this section, we briefly review the magnetization dynamics in our system. For a more comprehensive discussion, see Ref. [Naletov et al., 2011].

The local dynamics of the magnetisation $\bm{M}^{j}$, which depends continuously on the position $\bm{r}^{j}$ in the layer $j=(a,b)$, is described by the Landau-Lifshitz-Gilbert equation Gilbert (2004); Landau and Lifshitz (1965):

[TABLE]

Here, $\gamma<0$ is the gyromagnetic ratio in the magnetic layer. The first term on the right-hand side of Eq.(2) describes the adiabatic torque, that accounts for the precession of the magnetization vector around the local equilibrium direction. This precession axis is defined by the effective magnetic field experienced locally by the magnetization, $\bm{H}^{j}_{\text{eff}}=-\frac{\partial F}{\partial{\bm{M}}^{j}}$, which contains all the static contributions to the free energy $F$ of the layers Gurevich and Melkov (1996). In particular, the effective field contains contributions from applied field, exchange interaction, and dipolar interaction between the layers. We refer to Refs.Gurevich and Melkov, 1996; Naletov et al., 2011; Borlenghi et al., 2015b; Borlenghi et al., 2011 for the explicit expressions.

The second terms on the right-hand side of Eq.(2) is the damping torque

[TABLE]

proportional to the Gilbert damping parameter $\alpha^{j}$. We introduce here the notation $\bm{M}^{j}\equiv M_{s}^{j}\bm{m}^{j}$, with $M_{s}^{j}$ the norm of the magnetization (a constant of the motion) and $\bm{m}^{j}$ the unit vector along the magnetization direction. The dissipative term $\bm{d}^{j}$ is responsible for the finite linewidth (full width at half height) of the resonance peaks, $\Delta H=2d$. For a normally magnetized nanopillar, with circular precession of the magnetization, the simple relation $\alpha=|\gamma|\Delta H/(2\omega)$ holds Slavin and Tiberkevich (2009).

If a charge current $I_{dc}$ is flowing between two layers $j$ and $j^{\prime}$, the Slonczewski-Berger Slonczewski (1996); Berger (1996) spin transfer torque reads

[TABLE]

The latter depends on the relative angle between the magnetization $\bm{m}^{j}$ in the layer $j$ and the spin polarization of the current, which coincides with the direction of the magnetization $\bm{m}^{j^{\prime}}$ of the polarizer (here the thick layer). The term

[TABLE]

has the dimension of a distance. Here, $h$ is the Planck constant, $e$ the absolute value of the electron charge. $\eta$ is the spin polarisation of the current and $V$ the volume of the thin layer.

Since $\bm{d}^{j}$ and $\bm{S}^{j}$ are collinear, spin transfer torque can compensate the damping torque, as shown in Fig.1. When the dc current through the nano-pillar reaches the threshold current $I_{\rm{th}}=-2\pi\lambda\alpha H_{\rm{eff}}$, the thin layer starts auto-oscillating. Combining Eqs.(3) and (4), it is possible to define an effective damping for the thin layer,

[TABLE]

which depends linearly on the spin polarized dc current Sankey et al. (2006); Chen et al. (2008). The critical thresohold corresponds to the value of the current $I_{dc}$ at which the effective damping vanishes and the system starts auto-oscillating.

II.3 Coupled oscillator model and classification of the SW modes

Since our layers are thinner than 15 nm, one can assume that the magnetization dynamics is uniform along the thickness. In this approximation, the linearized LLG equation simplifies to two equations describing the circular precession of the transverse magnetization projections $M_{x}^{j}$ and $M_{y}^{j}$ around the $z$ axis, which depends only on the two spatial variables $(x,y)$ in the layer $j$ Gurevich and Melkov (1996); Naletov et al. (2011). The two real equations of each layer can be rewritten as one complex equation for the dimensionless spin-wave amplitude

[TABLE]

that depends on the polar coordinates $(r_{j},\phi_{j})$ of disk $j$.

The dynamics of the two disks, written in terms of the $c_{j}$s, is described by the equations Slavin and Tiberkevich (2009); Naletov et al. (2011); Borlenghi et al. (2014a)

[TABLE]

which are the equations of motion of two coupled nonlinear oscillators with resonance frequency $\omega_{j}(p_{j})$ and damping rates $\Gamma_{j}(p_{j})$. Both depend on the SW power $p_{j}=|c_{j}|^{2}$, which describes the amplitude of the oscillations in each disk. From hereon, to keep the notation simple we do not write explicitly the dependence on $p_{j}$. The frequencies $\omega_{j}=\gamma|H^{j}_{\rm{eff}}|$ are proportional to the local magnetic field, while the damping rates $\Gamma_{j-}$ are proportional to $\alpha_{j}\omega_{j}$. Both can be therefore controlled by means of the applied field along $\bm{z}$.

The terms $\Gamma_{j+}$, proportional to $I_{dc}$ are due to spin transfer torque, which can compensate the damping and lead to auto oscillations of the layers. In the present case, those terms do not have the same sign. At positive current, $\Gamma_{a+}$ is positive, while $\Gamma_{b+}$ is negative, favouring the auto oscillations in the thin layer $a$, and stabilising the thick layer $b$.

The dipolar coupling strength $h_{jj^{\prime}}$ is an effective term obtained by averaging the dipolar field over the volumes of the samples, see Refs. Naletov et al., 2011 for the explicit expression.

Eqs. (8) and (9) describe the dynamics of a nonlinear Schrödinger dimer, the simplest realisation of the discrete nonlinear Schrödinger equation (DNLS) Iubini et al. (2013); Borlenghi et al. (2014a, b). Upon multiplying Eqs.(8) and (9) respectively by $c_{a}^{*}$ and $c_{b}^{*}$ and summing them with their complex conjugate equations, one has the following continuity equation for the SW power

[TABLE]

and a similar equation for $p_{b}$. The magnetisation current $j^{p}_{ab}=2\mbox{Im}[h_{ab}c_{a}c^{*}_{b}]$ describe the transfer of $M_{z}$ between the two layers, and is essentially the SW current written for a discrete systems with only two spins Borlenghi et al. (2014a). Upon writing $c_{a}=\sqrt{p_{a}(t)}e^{i\phi_{a}(t)}$, the current reads

[TABLE]

The quantity $\beta$ comes from the condition of dissipative coupling between the oscillators Iubini et al. (2013); Borlenghi et al. (2015a). When the two oscillators are synchronised, $\phi_{a}\approx\phi_{b}$ and the magnetisation current approaches the constant value $j^{p}_{ab}\propto\sin\beta$. On the other hand, if the oscillators are not synchronised, the magnetisation current oscillates around zero and vanishes in average. Within this DNLS formulation, the spin transfer torque that appears in Eq.(10) plays the role of a magnon chemical potential, that by controlling the lifetime of the excitations, controls also the SW current between them.

The diagonalization of the LLG equation in a confined geometry leads to a discrete series of normal modes having each a different eigen-value, $\omega/(2\pi)$, the so-called Larmor precession frequency. The normal modes of the system are numbered according to the number of half waves in the vibration. In the case of a 2D axially symmetric structure, the normal modes are identified by two integers: $\ell$ and $m$, respectively the mode number in the azimuthal and radial directions. The analytical expression of the normal modes of a perpendicularly magnetized disk is found in Refs. [Damon and Eshbach, 1961; Naletov et al., 2011]

[TABLE]

where $J_{\ell}$ are the Bessel functions of the first kind and $k_{\ell,m}$ is the modulus of the in-plane SW wave-vector, which depends on the boundary conditions.

The above labelling can be extended to the case of two different magnetic disks coupled by dipolar interaction. In the perpendicular geometry, the strength of the dynamical dipolar coupling is attractive (lower in energy) when both layers vibrate in antiphase, because the dynamical dipolar charges in each layer are alternate Belmeguenai et al. (2007); Gubbiotti et al. (2004). Thus the binding state $B$ corresponds to a collective motion where the two layers vibrate anti-symmetrically and the anti-binding state $A$ to a collective motion where the two layers vibrate symmetrically. The $B$ modes correspond to a precession amplitude that is larger in the thick layer $b$, while the $A$ modes correspond to a precession amplitude that is larger in the thin layer $a$ Naletov et al. (2011), see Fig.1 for a cartoon.

The dynamical dipolar coupling does not modify the nature of the modes hence, in order to describe the dynamics of the bi-layer system, we shall just add a new index $B$ or $A$ indicating if the precession occurs in antiphase (mostly in the thick layer) or in phase (mostly in the thin layer), respectively. There are thus three indices to label the observed eigen-modes: the usual azimuthal and radial indices for a single disk ($\ell,m$), plus an additional index referring to the symmetrical or anti-symmetrical ($A$ or $B$) coupling between both layers.

The identification of the SW modes and their symmetry is essential to couple the oscillator to an external source, since SW modes can couple only to a source with the same symmetry. Here the $\ell$ index determines the rotational symmetry of the SW mode. The $\ell=0$ modes correspond to SW that do not rotate in the $x$-$y$ plane and can be excited only by a spatially uniform in-plane RF field, while the $\ell=1$ modes correspond to SW that rotate around the disk in the same direction as the Larmor precession, and can be excited only by an RF Oersted field with orthoradial symmetry. Thus, exciting the system with these different means gives two different spectra.

III Micromagnetic simulations at zero current

In this section we describe the SW spectra by means of micromagnetic simulations without spin-transfer torque. Those simulations were performed with the NMAG micromagnetic software Fischbacher et al. (2007), where the sample is described by finite element tetrahedral mesh. The latter has a maximum intersite distance of 6 nm, of the order of the Py exchange length. The micromagnetic parameters are the same used in Ref.[Naletov et al., 2011] and are reported in Tab.I for convenience. The dynamics at each of the $i=1,...,N$ nodes of the mesh of disk $j=(a,b)$ is described by the following LLG equation for the unit magnetisation vector, $\bm{m}_{i}^{j}=\bm{M}_{i}^{j}/M_{s}$:

[TABLE]

The integration of the LLG equation at each mesh site is performed by the Sundials ODE solver Hindmarsh et al. (2005), which is based on variable steps multistep methods. The field $\bm{H}^{j}_{i}$ at each mesh node has contributions from applied field, first neighbour exchange interaction and long range dipolar interaction, responsible for the coupling between the layers. The quantity of interest is the space-averaged magnetisation $\left\langle{\bm{m}}^{j}(t)\right\rangle=\frac{1}{V^{j}}\int_{V^{j}}{\bm{m}}^{j}(\bm{r}^{j},t)d^{3}{\bm{r}}$, which for our finite-elements mesh reduces to

[TABLE]

From this quantity, the SW amplitudes $c_{j}(t)$ are calculated. The power spectrum, shown in Fig. 2, is given by the Fourier transform of the time series of the collective SW amplitude averaged over the sample thicknesses: $c=(c_{a}t_{a}+c_{b}t_{b})/(t_{a}+t_{b})$.

The modes with $\ell=0$ (displayed in blue tones) are excited starting from an initial condition where the magnetization uniformly tilted $8^{\circ}$ in the x direction with respect to the precession axis $\bm{z}$. Instead, the the modes with $\ell=+1$ (displayed in red tones) are excited by applying to the magnetization aligned with the $\bm{z}$ axis the orthoradial vector field perturbation $\theta(r,z)=\epsilon\hat{\bm{z}}\times\hat{\bm{\rho}}$. Here $\epsilon=0.01$ and $\hat{\bm{\rho}}$ is the unit vector in the radial direction. Starting from these conditions, we have computed the time evolution of the system for 120 ns, with an integration time step of 5 ps. From the Fourier transform, the maximum frequency is 100 GHz and the frequency resolution 15 MHz. The frequencies of the SW modes are displayed in Tab.3 and 3 and compared with the experimental values.

IV Continuous Random Matrix Theory (CRMT) for spin transport

This section contains a thorough review of the CRMT semi-classical theory of spin-dependent transport in magnetic multilayers. We follow closely the material presented in Refs. [Waintal et al., 2000; Rychkov et al., 2009; Borlenghi et al., 2011] and we extend it by providing an explicit formula for spin torque, that will be used in micromagnetics simulations of next section.

IV.1 Scattering matrix approach

Within the scattering matrix formulation developed by Landauer and Buttiker Büttiker et al. (1985), a sample is defined by the scattering matrix $S$ which expresses the outgoing propagating modes in term of the incoming ones. The incoming modes are filled according to the Fermi-Dirac distribution of the leads to which they are connected. From the elements of the scattering matrix various physical quantities can be calculated, such as the conductance $G$, the spin currents $\bm{J}$ and charge current $I$.

The system contains $N_{ch}\gg 1$ propagating modes (or channels) per spin. In particular, one has $N_{ch}\approx A/\lambda_{F}^{2}$, where $A$ is the transverse area of the electrode and $\lambda_{F}$ is the Fermi wavelength. The amplitude of the wave-function on the different modes is given by the vector ${\bm{\psi}}_{i\pm}$ with $N_{ch}$ elements

[TABLE]

The latter contains the amplitudes for the right (+) and left (-) moving electron direction with spin $\sigma=\uparrow,\downarrow$ along the $z$-axis in region $i=0,2$ of the multilayer, see Fig.3 for a cartoon. The $S$ matrix is a $4N_{ch}\times 4N_{ch}$ unitary matrix that relates the outgoing modes to the ingoing ones:

[TABLE]

see Fig.3 (a) for a schematic of the system.

and consists of $2N_{ch}\times 2N_{ch}$ transmission $t,t^{\prime}$ and reflection $r,r^{\prime}$ sub-blocks.

[TABLE]

Here the $(r,t)$ and $(r^{\prime},t^{\prime})$ describe reflection and transmission respectively from left to right and from right to left. The transmission and reflection matrices have an internal spin structure:

[TABLE]

where $t_{\sigma\sigma^{\prime}}$ are $N_{ch}\times N_{ch}$ matrices containing amplitudes for transmission between $\sigma^{\prime}$ and $\sigma$ spin states, capturing both spin preserving and spin-flip phenomena.

The conductance of the system is given by the Landauer formula Büttiker et al. (1985)

[TABLE]

while the spin current in region 0 reads

[TABLE]

$\mu$ being the difference of chemical potential between the two electrodes Waintal et al. (2000).

IV.2 Random Matrix Theory (RMT)

The scattering matrix approach is fully quantum and contains interference effects such as weak localization or universal conductance fluctuations Waintal et al. (2000); Akkermans and Montambaux (2007). The system studied in this paper is a nanopillar of 250 nm of diameter connected to top and bottom electrodes. Those contain $\approx 10^{4}-10^{5}$ propagative channels. Here the scattering is not perfectly ballistic (mismatch at the interfaces, surface roughness or impurity scattering) so that channels get mixed up. Random Matrix Theory (RMT) Beenakker (1997); Mehta (2004); Waintal et al. (2000) assumes that this mixing is ergodic: an electron entering the system in a given mode will leave it in an arbitrary mode, acquiring a random phase in the process.

In this case, the transmission and reflection probabilities of an electron are well caracterized by their average over the propagative channels. This average is obtained by taking the trace over the $N_{ch}$ of the original reflection and transmission matrices. For instance, the hat matrix $\hat{t}$ Rychkov et al. (2009); Borlenghi et al. (2011) is defined as

[TABLE]

Explicitely, this reads

[TABLE]

with the same structure for the reflection hat matrix.

The elements of this matrix correspond to the probability for an electron with a given spin to be transmitted ($\hat{t}$) or reflected ($\hat{r}$) by the system. In particular, the terms $(1/N_{ch}){\rm Tr_{N_{ch}}}(t_{\uparrow\uparrow}t_{\uparrow\uparrow}^{\dagger})\equiv T_{\uparrow\uparrow}$ and $(1/N_{ch}){\rm Tr_{N_{ch}}}(t_{\downarrow\downarrow}t_{\downarrow\downarrow}^{\dagger})\equiv T_{\downarrow\downarrow}$ correspond to the probability to transmit an electron with up and down spin correspondingly. The terms at the corners, $(1/N_{ch}){\rm Tr_{N_{ch}}}(t_{\uparrow\downarrow}t_{\downarrow\uparrow}^{\dagger})\equiv T_{\uparrow\downarrow}$ and $T_{\downarrow\uparrow}$ correspond to transmission probabilities with spin flip. The so-called ”mixing transmission”, $(1/N_{ch}){\rm Tr_{N_{ch}}}(t_{\uparrow\uparrow}t_{\downarrow\downarrow}^{\dagger})=T_{mx}$ is a complex number whose amplitude measures how much of a spin transverse to the magnetic layer can be transmitted through the system, while its phase amounts for the corresponding precession. $T_{mx}$ decays exponentially with the size of the ferromagnet, and accounts for spin transfer effect Waintal et al. (2000); Borlenghi et al. (2011).

The other off diagonal elements in the hat matrix Eq.(21) can be ignored Waintal et al. (2000); Borlenghi et al. (2011). In the basis parallel to the magnetisation, Eq.(22) therefore becomes:

[TABLE]

with the same structure for the reflection matrix $\hat{r}$. The hat-matrix $\hat{S}$ has a form similar to Eq.(17),

[TABLE]

In order to describe transport in non-collinear multilayers, where the orientation of the magnetization changes inside the system, one needs to rotate the original $S$ matrix as $\tilde{S}=R_{\theta,\vec{n}}SR_{\theta,\vec{n}}^{\dagger}$ in the chosen working basis. Here the matrix

[TABLE]

is the rotation matrix of angle $\theta$ around the unit vector $\vec{n}$ that brings the magnetization onto the z-axis of the working basis. In term of hat matrices, this translates directly into

[TABLE]

with $\hat{R}_{\sigma\eta,\sigma^{\prime}\eta^{\prime}}=R_{\sigma\sigma^{\prime}}R^{*}_{\eta\eta^{\prime}}$ a unitary matrix.

From Eq.(23), the conductance is given by

[TABLE]

where the Sharvin resistance $\mathcal{R}_{sh}=\frac{h}{N_{ch}e^{2}}$, is a material property that can be experimentally measured, related to the number $N_{ch}$ of transverse propagative channels for the electrons crossing the system. Eq.(27) is analogous to the Landauer formula Eq.(19) and consists of the sum of all the possible transmission processes (spin preserving and spin flipping) for an electron.

In analogy with the modes $\psi_{i\pm\sigma}$ for the scattering matrix, the $4$-vectors ${\bm{P}}_{i\pm}$ are introduced

[TABLE]

The components o${\bm{P}}_{i\pm\uparrow,\downarrow}$ have interpretation in term of probabilities for an electron to propagate in the region $i$ of the system. The ”mixing” components, ${P}_{mx}$ are complex numbers which correspond to probability to find the electron with spin along the $x$ (real part) or $y$ (imaginary part) axis. Inside magnetic layers where the $z$ axis xill correspond to the direction of the magnetisation, they will correspond to the probability for the spin to have a part transverse to the magnetisation.

For the following discussion, it is convenient to consider a system made of two conductors connected in series, described by the two hat matrices $\hat{S}_{a}$ and $\hat{S}_{b}$, see Fig. 3 (b). The space is thus divided into three regions: region 0 and 2 respectively at the leftmost and rightmost part of the system, and region 1 in between the two hat matrices. In analogy with Eq.(17) which expresses the amplitudes of the outgoing modes in term of the incoming ones, in subsystems $a$ and $b$ one has respectively

[TABLE]

while for the total system $a+b$

[TABLE]

Here we have used the addition law of hat matrices, a fundamental property that will be useful in the derivation of spin currents and spin torque inside the system. According to this rule, given the hat matrices for separate systems $a$ and $b$, the hat matrices of the composed system $a+b$ read Rychkov et al. (2009); Borlenghi et al. (2011)

[TABLE]

with $\hat{1}$ the $4\times 4$ identity matrix. Similar expression holds for $r_{a+b}^{\prime}$ and $t_{a+b}^{\prime}$.

The main result of Refs.[Rychkov et al., 2009; Borlenghi et al., 2011] is that the spin current in region $i$ of the system can be expressed in terms of probability vectors as follows:

[TABLE]

where $\vec{\bm{\sigma}}=(\vec{\sigma}_{\uparrow\uparrow},\vec{\sigma}_{\uparrow\downarrow},\vec{\sigma}_{\downarrow\uparrow},\vec{\sigma}_{\downarrow\downarrow})$ is the vector of components of Pauli matrices. To express the spin current in region 1 in between the two conductors as a function of the hat matrices, we use Eqs.(33)-(42) to eliminate $P_{2+}$ and $P_{2-}$ and obtain

[TABLE]

The theory is completed imposing boundary conditions on the incoming electrons on both sides of the system. For normal electrodes one has,

[TABLE]

where $\mu_{0}$ and $\mu_{2}$ are the respective chemical potentials of the two electrodes. The generalization to magnetic electrodes is done imposing different chemical potentials for majority and minority electrons in the leads. By taking as boundary conditions $\mu_{0}=eU$, with $U$ the potential difference between the two sides and $\mu_{2}=0$, The spin current finally reads

[TABLE]

where the spin current per channel and per unit of potential difference reads

[TABLE]

Eq.(54) allows one to compute the spin current in the region between two bulk materials of arbitrary thicknesses. In the multilayer considered here, the currents $\bm{J}_{a/b\pm\delta}$ are calculated at positions $\delta=\pm 1$ nm before and after the two magnetic layers, as shown in Fig.4. The torque is the spin current absorbed by each layer, i.e. the quantity

[TABLE]

with $f_{j}=\bm{J}_{j-\delta}-\bm{J}_{j+\delta}$, $j=a,b$.

IV.3 From scattering matrices to CRMT

In order to calculate the current in different regions of the multilayer, one needs to calculate the matrices $\hat{S}(L)$ as a function of the position $L$ inside the system. The main result of Refs.[Rychkov et al., 2009; Borlenghi et al., 2011] is that the matrix $\hat{S}(L+\delta L)$ for an infinitesimal increment of the position $\delta L$ is entirely characterised by two matrices $\Lambda^{t}$ and $\Lambda^{r}$, defined as

[TABLE]

Once $\hat{S}(\delta L)$ is known, one can make use of the addition law Eqs.(41) and (42) to obtain a differential equation to compute $\hat{S}(L+\delta L)$. By taking the limit $\delta L\rightarrow 0$ one gets the two CRMT differential equations that describe hat matrices as a function of the length of the system:

[TABLE]

A bulk magnetic material is characterised by four independent parameters $\Gamma_{\uparrow}$, $\Gamma_{\downarrow}$ $\Gamma_{\rm sf}$ and $\Gamma_{\rm mx}$, which describe spin preserving and spin flip phenomena. With this parametrisation, one has for the transmission

[TABLE]

and for the reflection

[TABLE]

These four parameters correspond in turn to 5 different lengths. The two most important one are the mean free paths for majority ($\sigma=\uparrow$) and minority ($\sigma=\downarrow$) electrons defined as $\ell_{\sigma}=1/\Gamma_{\sigma}$. Next comes the spin diffusion length $l_{\rm{sf}}=[4\Gamma_{\rm{sf}}(\Gamma_{\uparrow}+\Gamma_{\downarrow})]-1/2$. Last come the complex number $\Gamma_{\rm{mx}}=1/l_{\perp}+i/l_{L}$ where $l_{\perp}$ is the penetration length of transverse spin current inside the magnet while $l_{L}$ is the Larmor precession length. Upon integrating Eq. (59) one obtains for the mixing transmission.

[TABLE]

which shows the exponential decay of the transverse spin current, absorbed by the layer, giving the phenomenon of spin transfer torque. The two lengths $l_{\perp}$ and $l_{L}$ are of the order of the nm Borlenghi et al. (2011), since spin transfer torque is a phenomenon that occurs close to the interface. For bulk materials, another important quantity is the current polarisation, $P_{\sigma}=\frac{T_{\uparrow\sigma}-T_{\downarrow\sigma}}{T_{\uparrow\sigma}+T_{\downarrow\sigma}}$. The latter decays exponentially as a function of the distance in the material, with spin-flip length scale of the order of 15 nm for Py, and reaches a constant value. The behaviour of spin and charge currents in a multilayer for different magnetic configurations have been extensively studied in Ref.[Borlenghi et al., 2011], to which we refer for a thorough discussion.

The CRMT parameters ($\Gamma_{\uparrow}$, $\Gamma_{\downarrow}$ and $\Gamma_{sf}$) are in 1-1 correspondence with the Valet-Fert (VF) parameters Valet and Fert (1993) [$\rho_{\uparrow}$, $\rho_{\downarrow}$ (resistivities for majority and minority electrons) and $l_{sf}$ (spin-flip diffusion length)]. Using the standard notations for the average resistivity $\rho^{*}$ and polarization $\beta$ [$\rho_{\uparrow(\downarrow)}=2\rho^{*}(1\mp\beta)]$ one has Rychkov et al. (2009):

[TABLE]

This parametrization does not fix the mixing coefficients $\Gamma_{mx}$, which only play a role in non collinear configurations Waintal et al. (2000).

In the VF theory, the interfaces are characterised by the average resistance $r^{b*}$ and polarisation $\gamma$ with $r_{\uparrow,\downarrow}=2r^{b*}(1\pm\gamma)]$. In the CRMT formalism, interfaces are modelled as a virtual material with their own transmission and reflection hat matrices. Those are related to the VF parameters as follows

[TABLE]

while for the reflection coefficient one has

[TABLE]

and $R_{\uparrow\downarrow}=R_{\downarrow\uparrow}=0$.

At this point the parametrisation of CRMT is complete. Eqs.(58) and (59) have been solved numerically to simulate a one dimensional magnetic multilayer. Each material and interface are simulated separately, and the whole system is recovered by applying the addition law Eqs.(41) and (42).

In order to calculate magnetoresistance and spin torque as a function of the magnetic configuration, we simulate transport the nano-pillar depicted in Fig.4 by keeping the magnetisation $M_{b}$ fixed and rotating the magnetisation $M_{a}$ of an angle $\theta$ in spin space. Such rotation is obtained by applying to the hat matrices of the thin layer the transformation Eq.(26).

Fig.5a) and b) shows respectively the angular dependence of spin torque and resistance inside the nanopillar for a single propagative channel and with $eU=1$. The resistance difference between parallel and antiparallel configuration is 27 m$\Omega$. We remark that this calculation, performed without adjustable parameters, reproduces experimental data within the $10\%$.

We note also that in our system, the current is spin-polarized by the thick layer and impinges the thin layer, exerting a torque that destabilizes its magnetization. Because of multiple reflections of spin polarized electrons between the two Py layers, STT tends to stabilize the thick layer, increasing its effective damping. Thus, at low current, one can consider the thick layer as ”fixed” and the thin layer as ”free”. However, at high current both layers undergo a coupled precession (see next section).

V Simulations of current driven dynamics

This section contains the main results of the paper. Here solve simultaneously the transport equations coupled to the LLG equation with spin transfer torque, by coupling CRMT to Nmag. The approach described here is valid for one dimensional systems, where the magnetisation varies only along the direction of propagation of electrons. In the present case the magnetisation is uniform along $z$ inside the material, but its dynamics is different in each layer. A generalization of CRMT for a fully three dimensional spin transport has been recently developed Petitjean et al. (2012). However, for the case considered here the one dimensional CRMT approach captures the physics well and is extremely fast. The other advantage of CRMT is that it is parametrised by the same set of experimentally accessible parameters as the VF theory (reported in Tabs.5 and 5), so that no free parameter is needed to characterise realistic systems and materials.

To include the effect of spin transfer torque into Nmag, the LLG equation Eq.(13) at node $i$ of disk $j=a,b$ needs to be modified as follows

[TABLE]

Here $g$ is the Landé factor, $\mu_{B}$ is the Bohr magneton and $C$ is the volume associated to each site of the mesh. $\bm{\tau}^{j}_{i}$ is the torque given by Eq.(56) For our relatively homogeneous mesh, we have taken this volume as the total volume $V_{j}$ of disk $j=a,b$ divided by the total number of sites $N_{j}$. For simplicity, we have defined the vector $\hat{\bm{w}}_{i}^{j}=\bm{m}_{i}^{j}\times(\bm{m}_{i}^{j}\times\bm{m}_{i}^{j^{\prime}})$.

To include STT in our micromagnetic simulations, we adopt the following self-consistent loop:

At time $t$ and site $i$ in disk $j$, Nmag computes the vector $\bm{m}_{i}^{j}$ 2. 2.

The magnetic configuration is passed to the CRMT solver which computes the torque $\bm{\tau}^{j}_{i}$ 3. 3.

Then the quantity $\bm{m}_{i}^{j}(t)+{\bm{\tau}_{i}^{j}}(t)\Delta t$ is set as new initial condition for the Nmag solver, which performs the time integration of the LLG equation between times $t$ and $t+\Delta t$.

At this point spin torque is re-calculated as a function of the new magnetic configuration and the loop starts again. The integration time step $\Delta t$ is of the order of the ps.

Since Nmag uses a finite element mesh while CRMT a finite difference discretisation of space, to implement the coupling, we divide the mesh into columns of sites, each column $k$ representing a CRMT pillar with cross section $S_{k}$. The whole system is then considered as an assembly of CRMT pillars connected in parallel, as shown in Fig.6.

Each column contains two sites, lying on surfaces of thin ($a$) and thick ($b$) layer facing each other, with magnetizations $\bm{m}_{ak}(t)$, and $\bm{m}_{bk}(t)$ correspondingly. From the angle $\theta_{k}$ between these two sites, one obtains the torques $\bm{\tau}^{j}_{k}(\theta)$ that act on the magnetisation vectors $\bm{m}^{j}_{k}$, and the resistance $R_{k}(\theta_{k})$.

This procedure takes into account the three dimensional texture of the magnetization, while electronic transport is along ${\bm{z}}$ only, without considering the lateral diffusion of spins.

Since the system is an assembly of columns connected in parallel, the current $I_{k}$ flowing in each column is given by the total current $I_{\rm{dc}}$ divided by the number of columns $N_{k}$, which in our system corresponds to the number of sites lying at the surface of each disk. Each column has a cross section $S_{k}\approx S/N_{k}$. The current $I_{k}$ is related to the potential difference between the conductors via Ohm’s law: $eU=R_{k}(\theta)I_{dc}/N_{k}$, where the resistance $R_{k}$ depends on the transmission probability $T_{k}(\theta)$ for an electron to cross column $k$:

[TABLE]

Inserting Eq.(78) into Eq.(LABEL:eq:torque_nmag1), and recalling that $\mathcal{R}_{\rm sh}=h/(e^{2}N_{ch})$, the torque finally reads

[TABLE]

We remark that Eq.(79) is quite general and can be applied to systems with arbitrary geometry, provided that electronic transport is one dimensional only.

VI Results and discussion

We turn now to the discussion of the coupled Nmag-CRMT simulations. A qualitative description of the different dynamical regimes is given by the time evolution of $M_{z}$, displayed in Fig.7. Panel (a) shows the region beyond critical current, where both the magnetisations are fixed and aligned with the $\bm{z}$ axis. Panels (b) shows the region slightly beyond critical current, where the thin layer precesses and the thick one remains fixed. Finally, panel (c) displays the high current regime, where the thin layer is reversed and both layers undergo coupled precession. This behaviour is due to the repulsive character of dipolar interaction: as the magnetisation of the thin layer reverses, it repels the magnetisation of the thick one, causing it to precess.

Fig.8 shows the resistance as a function of the dc current, for values of the applied field between $1$T ands $1.3$T. One can recognise clearly the three dynamical regions discussed above. In particular, the resistance starts increasing around $-6$ mA, indicating the region of critical currents. The resistance increases with the current until it reaches a peak between $-13$ and $-16$ mA, corresponding to the reversal of the magnetisation ${\bm{m}}^{a}$ of the thin layer.

Note that at increasing field the resistivity graph moves rightwards towards higher currents. This is due to the fact that the precession frequency increases proportionally to the applied field, so that the damping $\Gamma_{a+}=\approx 2\alpha_{a}\omega_{a}$ also increases and a higher current is needed to compensate it.

The power spectrum of the system, with the modes $\ell=0$ and $+1$, is shown respectively in Fig.9 (a) and (b) in logarithmic scale. The modes $A_{00}$ and $A_{10}$, corresponding to the dynamics of the thin layer, increase with the current and dominate the spectrum around $I_{dc}=-5$ mA. Beyond the critical threshold, their frequency increases with the field. On the other hand, the modes $B_{00}$ and $B_{10}$, which dominate the spectrum at zero currents, decrease until they almost disappear near the critical threshold.

Fig.10 shows the linewidth and frequencies of the modes $A_{00}$ and $A_{10}$ as a function of the dc current. In panel (a), one can see that the linewidths of both modes decrease linearly of an order of magnitude at increasing current, and they vanish at the critical threshold, where the damping is compensated and the system begins to auto oscillate. Note that the two modes have slightly different critical currents, of about $5$ ($\ell=0$) and $5.5$ ($\ell=1$)mA. The behaviour of the spectrum and the critical currents agree with the experimental result of Refs.[Borlenghi, 2011; Naletov et al., 2011]. The frequencies of the excited modes remain constant until the critical threshold, and then starts increasing linearly with a rate around 2 GHz per mA. This behaviour also agrees with experimental data.

We conclude this section by noting that the physics discussed here is very similar to that of the spin-caloritronics diode Borlenghi et al. (2014a). In Fig.11, one can see that increasing the electrical dc current leads to an increase of the magnon current between the two disks. This current describes the transfer of magnetic moment $M_{z}$ between the two disks, and corresponds to the usual SW current written for a system of only two spins Borlenghi et al. (2014a, b, 2015a, 2015b).The increase of the SW current with the dc current is due to the fact that spin transfer torque excites only one mode, which eventually dominates the spectrum, inducing a phase locking between the two disks. In this kind of discrete systems, the magnon current is indeed a measure of the phase synchronisation of the system.

VII Conclusions

By coupling CRMT to micromagnetic simulations with the Nmag code we allow description, on an equal footing and without free parameters, of both transport and magnetic degrees of freedom. The results of this work have important consequences regarding the characterization and the optimization of the performance of STNOs.

Using the developed method, we have identified the nature of the modes that auto-oscillate when the current exceeds the critical threshold. In particular, we have predicted different critical currents for the modes $A_{00}$ and $A_{10}$, and a nonlinear frequency shift of the order of 2 GHz per mA. The precise determination of the SW mode symmetry of the auto-oscillating mode is important for the phase synchronization of a STNO to an external source. In fact, it will be successful only if the latter can couple efficiently to the spin transfer driven auto-oscillation modes, i.e., if it has the appropriate symmetry.

The flexibility of CRMT and of the procedure to couple it with Nmag allows one to use our numerical tool to simulate different geometries and materials. Further development and investigations are possible. A multiscale approach which combines systematically CRMT with a fully quantum approach has been implemented Borlenghi et al. (2011). This should allow to compute current driven dynamics in a large variety of systems, including multi-terminal devices and tunnel junctions.

We remark that in our simulations we have not taken into account the lateral diffusion of spins, since the system is described using one dimensional CRMT columns, where electrons propagate only along the $z$ direction.

This 1D model of transport is effective to describe selection rules into a perpendicularly magnetized nanopillar with magnetic field applied along ${\bm{e}}_{z}$, but more complicated configurations (such as magnetic vortexes and multi-terminal spin valves) may require a fully three dimensional description.

Finally, further investigation is necessary to understand the behaviour of the system at high current, where nonlinear effects (such as the dependence of the Gilbert damping on current Slavin and Tiberkevich (2009)) may play an important role. The present work can be considered as an intermediate step towards a fully 3D description of transport and magnetization dynamics in realistic systems.

VIII Acknowledgements

We wish to thank R. Lassalle-Balier and J. Dubois for fruitful discussions, X. Waintal, O. Klein and G. de Loubens for useful comments and assistance in the analysis of the experiments. Financial support from Vetenskapsrådet (VR), Carl Tryggers Stiftelse (CTS), Stiftelsen Olle Engqvist Byggmästare and Swedish Energy Agency (STEM) is gratefully acknowledged.

The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Center (NSC), Linköping University, the PDC Centre for High Performance Computing (PDC-HPC), KTH, and the High Performance Computing Center North (HPC2N), Umeå University.