Naoufel Ben Abdallah

Elise Fouassier

Clément Jourdana

David Sanchez

On a model of magnetization switching driven by a spin current: a multiscale approach *

Keywords: Spin transfer, spintronics, micromagnetics, Landau-Lifshitz equation, numerical simulations, multiscale analysis, asymptotic expansions. AMS subject classifications. 35B40, 35Q60, 41A60, 65Z05, 78A25

We study a model of magnetization switching driven by a spin current: the magnetization reversal can be induced without applying an external magnetic field. We first write our one dimensional model in an adimensionalized form, using a small parameter ε. We then explain the various time and space scales involved in the studied phenomena. Taking into account these scales, we first construct an appropriate numerical scheme, that allows us to recover numerically various results of physical experiments. We then perform a formal asymptotic study as ε tends to 0, using a multiscale approach and asymptotic expansions. We thus obtain approximate limit models that we compare with the original model via numerical simulations.

Introduction

In this article, we are concerned with a model of magnetization reversal induced by a spin polarized current. The idea is of switching the magnetization of a ferromagnetic material without applying an external magnetic field. The magnetization is rather reversed by an additional spin transfer torque. This concept was discovered by Slonczewski [START_REF] Slonczewski | Current-driven excitation of magnetic multilayers[END_REF] and Berger [START_REF] Berger | Emission of spin waves by a magnetic multilayer traversed by a current[END_REF] in 1996. Spin transfer appears to be a turning point in spintronics and is the subject of an extensive research in physics. The particular phenomenon we study here is of great interest to construct magnetic memories. The physical device proposed by [START_REF] Slonczewski | Current-driven excitation of magnetic multilayers[END_REF] and [START_REF] Berger | Emission of spin waves by a magnetic multilayer traversed by a current[END_REF] is a magnetic multilayer mainly composed of two ferromagnetic layers separated by a nonmagnetic spacer layer (see Figure 1). In a schematic way, it works as follows: the first ferromagnet F -, being thick, acts as a polarizer, the second ferromagnet F + , being thin, is to be switched.

Initially, the two magnetizations m -(t = 0) and m + (t = 0) respectively in F -and F + are polarized in different directions (θ is the angle between m -(t = 0) and m + (t = 0), see Figure 1). This is done by taking two ferromagnetic materials with different anisotropy directions. Roughly, we can say that a ferromagnetic material is a material having a preferred direction for the magnetization (the anisotropy direction): when there is no external field applied, then at equilibrium, the magnetization is parallel to this anisotropy direction.

In order to induce the switching of the magnetization in the thin layer, an electric current is injected in the device at the left end of the thick layer F -, perpendicularly to the layers (along the x axis). One can see [START_REF] Myers | Current-induced switching of domains in magnetic multilayer devices[END_REF] for experimental details for instance.

The first ferromagnetic layer F -is thick (around 100 nm), so that it acts as a polarizer. Its magnetization does not move under the spin current influence, whereas the spin density is polarized in the direction of the magnetization when passing through the F -layer. On the contrary, the second ferromagnetic layer F + is thin (around 1-5 nm). Its magnetization can move driven by the spin current. Indeed, the spin density, which is polarized along m -crossing F -, reaches F + with a transverse component with respect to m + . Then, a torque is established between this transverse component and the magnetization m + . When the spin transfer is intense enough (that is to say when the injected current is enough important), m + can be completely switched. Physical experiments request a very intense current to observe the switching.

Last, another non ferromagnetic layer is at the right end of the device. It allows to reduce the spin density (acting like a "sandbox").

The mathematical model we study here is a slightly simplified version of that introduced by Zhang, Levy and Fert [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF] (see also Shpiro, Levy and Zhang [START_REF] Shpiro | Self-consistent treatment of nonequilibrium spin torques in magnetic multilayers[END_REF] for details). In each layer, we write a system of coupled PDEs on the spin density and the local magnetization. The local magnetization satisfies a Landau-Lifshitz equation with an additional term in the effective field due to the spin current. For the spin density, the main particularity of this model lies in the description of spin transport using macroscopic spin diffusion in the entire structure, the treatment of interfaces scattering being done via boundary conditions. The equation satisfied by the spin density is thus a diffusion equation with an additional torque term where the magnetization appears. All the details are written in the next section.

In [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF], [START_REF] Shpiro | Self-consistent treatment of nonequilibrium spin torques in magnetic multilayers[END_REF], the authors considered a framework in which there are various scales for the different phenomena involved in the magnetization reversal process. Rather than making a priori approximations such as in [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF], in order to study these different scales, we introduce a small parameter, and we write adimensionalized equations. The definition of this parameter and the adimensionalization are detailed in Section 2.2. Let us just say here that, once the adimensionalization is done, the thick layer is of thickness 1 and the thin layer of thickness ε.

Our choice to introduce only one parameter could be discussed and is of course a simplification. But the approximate model we obtain with this hypothesis has good properties with respect to the full model in the limit ε → 0, which is a way to validate our choice a posteriori.

The several aspects of our results on the adimensionalized equations (see ( 14)-( 15)-( 16)-( 17)) are the following:

1. Description of the various scales. We first explain in schematic frameworks the phenomena appearing at the different scales in the model. The first point concerns space scales: it is the existence of a boundary layer of thickness ε at the right end of the thick layer. The second important point concerns time scales: it is the existence of very fast temporal oscillations at scale t/ε 4 for the spin density and the presence of time boundary layer at scale t/ε 2 . More precisely, the spin density in the thick layer has a component that is not parallel to the magnetization of this layer and that is oscillating (at scale t/ε 4 ), whereas the thick layer should act as a polarizer. But since this transverse component has zero mean value, it will be possible to avoid it in numerical calculations and in the asymptotic model.

2.

Numerical calculations for the full model. First, using a very small time step, we observe the very fast oscillations at scale t/ε 4 on the transverse part of the spin density in the thick layer as well as the decay at scale t/ε 2 . But, for the phenomenon we are interested in -the magnetization reversal-, these oscillations are not relevant, and of course calculations with a very small time step are not very efficient. So we construct a numerical scheme, adapted to the different scales of the model, that does not resolve the very fast oscillations. This will be done by using an implicit scheme for the temporal discretization of the equation on the spin density (with a time step much bigger than ε 4 ) that allows us to "filter" the oscillations.

Moreover, since there is a boundary layer near the interface in the thick layer, it is more efficient to use an adapted space mesh, refined near the interface only. Concerning the space discretization, we also take into account the fact that the behaviours of the parts of the spin density that are parallel or transverse to the magnetization are different by using adapted basis functions for the finite element method.

3. Formal asymptotic study. In order to handle very small values of the parameter ε, it is interesting to derive an asymptotic model that will be in particular much easier (and quicker) to simulate. So, we then perform a formal asymptotic study as ε tends to 0 in order to recover a simplified asymptotic model. We use a multiscale approach, introducing the variable x/ε to treat the boundary layer. Again, we do not want to compute the oscillating part or the time decay, so we omit the time scale t/ε 4 and we work with well-prepared initial data to avoid the study of the time-boundary layer at scale t/ε 2 . The first idea is to use a basic asymptotic expansion ansatz: we develop each quantity as a sum of powers of ε, and we compute the first profiles. Our limit model shows the "right" properties with respect to physical experiments, up to magnetization reversal. After it, whereas in physical experiments it is possible to recover the original configuration by changing the sign of the injected current, it is not possible with this first asymptotic model.

So, we then look for a better model. One quite usual idea is to allow the profiles in the previous ansatz to depend on the parameter ε (roughly, one can think that we replace a first profile of type m 0 by m ε 0 = m 0 + εm 1 ). The parameter ε still appears in the second asymptotic model thus obtained, but the numerical calculations for it are very quick, even for small values of ε.

4.

Numerical comparison between the approximate models and the full model. Once obtained the asymptotic models, we discuss their validity by making numerical comparisons with the full system.

The model we are interested in in this article has already been studied in previous works. There are several papers concerning theoretical results (existence of solutions), see for instance [START_REF] Hamdache | On a model of magnetization switching by spin-polarized current[END_REF], [START_REF] García-Cervera | Spin-polarized transport: existence of weak solutions[END_REF]. On numerical aspects, we would like to mention the works done by Garcia-Cervera, Wang and E, [START_REF] García-Cervera | Spin-polarized currents in ferromagnetic multilayers[END_REF], [START_REF] Wang | A Gauss-Seidel projection method for micromagnetics simulations[END_REF]. In [START_REF] García-Cervera | Spin-polarized currents in ferromagnetic multilayers[END_REF], they proposed a numerical scheme for the model of magnetization reversal we consider here. They gave three dimensional simulations but they did not study the case considered by [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF]. They treated Permalloy ferromagnetic multilayers for which ε is not very small (ε = 0.4). As far as we know, our work is the first one using a multiscale approach and studying asymptotic properties of this in the same framework as [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF].

This paper is organized as follows. We first explain in details how is derived the model introduced in [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF], and we write adimensionalized equations using the small parameter ε. In Section 3 and 4, we present our numerical results for the model we described in the previous section. Then, in section 5 we present our formal derivation of two asymptotic models, and we give some numerical comparisons between the original model and the limit models.

2 The spin transfer model

The model proposed by Zhang, Levy and Fert

In this section, we detail the model proposed by Zhang, Levy and Fert [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF] for the phenomenon of magnetization-switching we have described in the introduction. More precisely, we detail a slightly simplified model. We do not discuss here other kinds of models that have been proposed for the magnetization reversal, we refer to [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF], [START_REF] Shpiro | Self-consistent treatment of nonequilibrium spin torques in magnetic multilayers[END_REF] and the references therein. Let us just mention a more kinetic approach by Stiles and Zangwill [START_REF] Stiles | Anatomy of spin-transfer torque[END_REF]. Our aim is to perform both a numerical study and an asymptotic study of the model proposed by [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF] in order to recover the experimental results. Let us consider a magnetic multilayer with the current perpendicular to the plane of the layer (defined as the x-direction). Due to this geometry, we assume in the sequel that all the quantities only depend on the x-space variable (the direction perpendicular to the layers). Hence, the model we present is a one dimensional model (with respect to the space variable).

As suggested by Zhang, Levy and Fert [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF], we make the two following simplifications. First, we assume that the non-magnetic spacer layer (NF) is infinitely thin so that the quantities related to the spin current are conserved in this region. That is to say that we replace this layer by an interface, and put some continuity conditions on the spin density at this interface. Second, we assume that the thin ferromagnetic layer F + is backed by an ideal paramagnetic layer (some adequate boundary conditions at the right end of the device will be imposed), and we neglect the spin reflection at the interfaces.

Remark. The non-magnetic spacer layer (NF) normally induces a decay of the spin density s between the two magnetic layers. Without the infinitely thin assumption we have to introduce and study another thin layer or to replace it by equivalent transmission conditions obtained by solving a model equation for the non-magnetic spacer layer. As long as this layer remains very thin the attenuation factor is not too small and our study remains valid up to the modification of some coefficients to take into account the decay.

Notations. Let L be the thickness of the thick ferromagnetic layer F -and l the thickness of the thin ferromagnetic layer F + . We assume that F -corresponds to x ∈ (-L, 0), F + to x ∈ (0, l), and that the interface separating the layers is at x = 0.

The quantities we are interested in are the charge density ρ : R + × [-L, l] → R, (t, x) → ρ(t, x), the spin density s : R + × [-L, l] → R 3 , (t, x) → s(t, x) and the local magnetization m : R + × [-L, l] → R 3 , (t, x) → m(t, x).

For u and v two vectors in R 3 , we denote by u • v their scalar product, u × v their vector product and | u|, | v| their euclidian norms. We recall that u

× v • w = v × w • u = w × u • v.
In the sequel, ∂ t , ∂ x , ∂ 2 x denote the first derivative with respect to the t-variable, the first and second derivatives with respect to the x-variable respectively.

The equations in one ferromagnetic layer

Let us first focus on the equations written in one ferromagnetic layer, in order to describe the interaction between the spin accumulation and the local magnetization. The main aspect, compared to previous models, is that the model proposed by Zhang, Levy and Fert takes into account the effects of spin diffusion. In this model, on one hand, the spin density satisfies a diffusion equation with a term corresponding to a precession phenomenon around the magnetization. On the other hand, the magnetization is solution of a Landau-Lifshitz equation with an additional spin torque.

The equations on the charge density and the spin density

The system of equations on the charge density ρ and the spin density vector m is derived from a diffusion equation satisfied by the 2 × 2 matrix distribution function n(t, x) depending on the time variable t and the space variable x. The matrix n(t, x) is hermitian. If we denote by I 2 the 2 × 2 identity matrix and by σ, the vector of Pauli matrices

σ 1 = 0 1 1 0 , σ 2 = 1 0 0 -1 et σ 3 = 0 -i i 0 ,
we decompose the matrix distribution function n in the basis (I 2 , σ 1 , σ 2 , σ 3 ) of the space of 2×2 hermitian matrices as n = ρI 2 + s

• σ = ρI 2 + s 1 σ 1 + s 2 σ 2 + s 3 σ 3 ,
where ρ(t, x) ∈ R denotes the charge density and s(t, x) ∈ R 3 the spin density vector. Let , Ĉ, D be the matrices representing the current, the conductivity and the diffusion constant. We write

Ĉ = C 0 I 2 + C • σ, D = D 0 I 2 + D • σ,  = j e I 2 +  m • σ,
where j e and j m are respectively the electric current and the magnetic current. Then for diffusive transport, the current is given by

 = ĈE(x) -D∂ x n,
where E is the electric field. From this, we get

j e = 1 2 Re Tr() = C 0 E(x) -D 0 ∂ x ρ -D • ∂ x s,  m = 1 2 Re Tr(-σ) = CE(x) -D∂ x ρ -D 0 ∂ x s,
where the second line is to be read component by component as

(  m ) k = 1 2 Re Tr(-σ k ) , for k = 1, 2, 3.
The equation on n writes

∂ t n + ∂ x  + i 2 Ω • σ, n = Tr(n) 2 I 2 - n τ sf . (1)
where

i 2 Ω • σ, n = -Ω × s • σ.
This term describes the rotation effect of s around the effective field Ω (precession phenomenon). Now, we have to precise the interaction with the local magnetization m(t, x). It appears by several ways: through the diffusion constant and through the effective field. Following [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF], we introduce two spin polarization parameters β, β ′ and a parameter J quantifying the interaction between the spin accumulation and the local magnetization and we write

C = βC 0 m, D = β ′ D 0 m, Ω = J m,
where D 0 is the diffusive constant of the metal and is the reduced Plank constant. Using all these relations, we get the equations on the charge density ρ and the spin density vector:

         ∂ t ρ + ∂ x j e = 0 with j e = C 0 E -D 0 ∂ x ρ -β ′ D 0 m • ∂ x s, ∂ t s + ∂ x j m + J ( s × m) = - s τ sf with j m = βC 0 E m -β ′ D 0 m∂ x ρ -D 0 ∂ x s, (2)
where • = h 2π with h the Plank constant : h = 6, 62.10 -34 J.s, • J quantifies the interaction between s and m : J = 0.1 -0.4eV = 1.6 -6.4.10 -20 Joules,

• τ sf is the relaxation time of spin switching, evaluated in the order of 10 -12 s,

• D 0 is the diffusive constant of the metal in the order of 10 -3 m 2 .s -1 ,

• β ′ is a parameter of spin polarization.

Modeling simplification

As it is done in [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF], [START_REF] Shpiro | Self-consistent treatment of nonequilibrium spin torques in magnetic multilayers[END_REF], dropping the terms containing ∂ x ρ, we first rewrite

j m = βj e m -D 0 ∂ x s -ββ ′ ( m • ∂ x s) m , (3)
and we then neglect the first term in this expression of the magnetic current.

Since the system is quite simplified by putting ββ ′ = 0, without changing its structure, we chose to make one more simplification here: we only keep

j m = -D 0 ∂ x s.
All the calculations we present in the sequel will be easier to write (in particular the asymptotic study) in that framework.

Moreover, with these assumptions, the equation on the charge density ρ is not coupled anymore with the other quantities, so we omit this first equation in the sequel, and the equation on s now writes

∂ t s -D 0 ∂ 2 x s + J ( s × m) = - s τ sf . (4)

The equation on the local magnetization

The local magnetization m(t, x) satisfies the following Landau-Lifshitz equation (LL)

∂ t m = -γ m × ( H e + J s) + α m × ∂ t m, (5)
where parameters γ > 0 and α > 0 are respectively the gyromagnetic coefficient and the Gilbert damping parameter.

The coupling with the spin current appears via an additional term in the effective field: J s. The usual effective field appearing in the Landau-Lifshitz equation is H e . It includes contributions from the external, anisotropy and demagnetizing field, and an exchange term:

H e = H ext + ∇ m ψ( m) + H demagn + ν∂ 2 x m. (6)
In our study, we make the following assumptions:

• First, we assume that there is no external magnetic field applied: H ext = 0.

• Second, we take an anisotropy energy ψ( m) of the form:

ψ( m) = c 2 ( m • u) 2
, where c is a positive constant, and u is the anisotropy direction, a unit vector chosen orthogonal to the x-direction (i.e. in the plane of the layer).

• Third, since we study a one-dimensional model, the demagnetizing field has a simple form. Indeed, it fulfills the following equations in the whole space:

div H demagn + m = 0, curl H demagn = 0,
where m denotes the extension of m by 0 outside the ferromagnetic domain. In our one dimensional case, the demagnetizing field fulfills:

∂ x ( H demagn + m) • e x = 0, ∂ x e x × H demagn = 0.
Assuming moreover that H demagn is in L 2 (R), this implies that the demagnetizing field is a local field which writes

H demagn = -( m • e x ) e x .
This demagnetizing field is also the one obtained in thin layer models (see [START_REF] Carbou | Thin layers in micromagnetism[END_REF][START_REF] Sanchez | Behaviour of the landau-lifschitz equation in a periodic thin layer[END_REF]).

Hence, in the sequel, we work with an effective magnetic field ( 6) containing only the contributions from the anisotropy field, the demagnetizing field and an exchange term.

H e = c( m • u) u -( m • e x ) e x + ν∂ 2 x m. (7)

Boundary and interface conditions

To complete the description of the model, we have to talk about the boundary and interface conditions (at x = -L, x = 0 and x = l) both for the spin density and the magnetization.

For the Landau-Lifshitz equation, we choose, as usual, homogeneous Neumann conditions on the boundaries of each layer (both for F -and F + ). Consequently, we have

   ∂ x m(t, -L) = 0, ∂ x m(t, 0 -) = 0, ∀t ≥ 0. ∂ x m(t, 0 + ) = 0, ∂ x m(t, l) = 0, (8)
Then, for the spin density, the conditions we impose are the following:

• at x = -L, we put a (non homogeneous) Neumann condition. This value corresponds to the injected current:

∂ x s(t, -L) = j L (t) ∀t ≥ 0, (9)
where j L is a given function.

• at the interface x = 0, we preserve the continuity of the spin density s

s(t, 0 -) = s(t, 0 + ) ∀t ≥ 0, (10)
and the continuity of the current j m :

∂ x s(t, 0 -) = ∂ x s(t, 0 + ) ∀t ≥ 0, (11)
• at x = l, we want the system to have a free evolution as if a non ferromagnetic layer is at the end of the device ("sandbox"). A first alternative is to take a homogeneous Neumann condition

∂ x s(t, l) = 0 ∀t ≥ 0. ( 12
)
A second alternative is to establish a more relevant Fourier-Robin condition by solving the stationary problem on s with a infinite layer thickness and without magnetic field, that is:

-∂ 2 x s + s = 0, x ∈ (0, +∞). Keeping only L 2 solutions, we get the following boundary condition at x = l:

∂ x s(t, l) = -s(t, l).
But, since the numerical results (for the scaled system, this condition becomes ∂ x s(t, l) = -ε s(t, l), see next paragraph) with these two conditions are very similar, we chose to use the homogeneous Neumann condition.

The scaled model

In order to write a dimensionless system, we scale the space variable as x L in the thick material F - and as x l in the thin material F + . For time scales, the characteristic time scale for the magnetization corresponds to t 0 = γ -1 H -1 (with H = | H e + J s|). This term is of order of the nanosecond. We scale the time variable as t t0 . Using these notations, equation ( 4) rewrites, for instance in the thin layer, as

1 λ 2 0 ∂ t s - 1 l 2 ∂ 2 x s + 1 λ 2 J s × m + 1 λ 2 sf s = 0, where λ 2 sf = D 0 τ sf , λ 2 J = D 0 /J, λ 2 0 = D 0 t 0 , (
the equation obtained in the thick layer is the same with l replaced by L).

As it is indicated in [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF] or [START_REF] Shpiro | Self-consistent treatment of nonequilibrium spin torques in magnetic multilayers[END_REF], we assume that we have the following orders:

λ 2 sf = O(10 -15 ), λ 2 J = O(10 -18
), and λ 2 0 = O (10 -12 ). This order of magnitude of λ sf is typical for a transition-metal ferromagnet, e.g. Co (but is not valid anymore for Py for instance), see [START_REF] Shpiro | Self-consistent treatment of nonequilibrium spin torques in magnetic multilayers[END_REF]. We then define ε as the ratio

ε = λ J λ sf , (13)
so that ε 2 is of the order of 10 -3 and that we also have that λ sf λ 0 is of the order of ε.

Then we assumed that the thickness of the thin layer is of the order of 10 -9 m so that l 2 is of the same order as λ 2 J in the thin layer, whereas for the thick layer, we assume that L 2 = O(10 -15 ) is of the same order as λ 2

sf . Hence, we also have

ε = l L .
With that scaling, we denote s -: R + ×(-1, 0) → R 3 and m -: R + ×(-1, 0) → R 3 respectively the spin density vector and the local magnetization in the thick ferromagnetic material F -, and s + : R + ×(0, 1) → R 3 and m + : R + × (0, 1) → R 3 respectively the spin density vector and the local magnetization in the thin ferromagnetic material F + .

We still denote t and x the new time and space variables. We get the following scaled equations for the spin densities in the thick and thin layers:

ε 2 ∂ t s --∂ 2 x s -+ s -× m - ε 2 + s -= 0, x ∈ (-1, 0), ε 2 ∂ t s + - 1 ε 2 ∂ 2 x s + + s + × m + ε 2 + s + = 0, x ∈ (0, 1).
Moreover, rather than keeping track of the fact that the injected current must be very intense, of order 1/ε, through a Dirichlet data of that order at x = -L, we scale the spin density by a factor J/H. Hence, we will keep track of the fact that the injected current must be very intense through the factor 1 ε in front of s in the LL equation.

We reduce to the following dimensionless equations, written for t ≥ 0,

       ε 2 ∂ t s --∂ 2 x s -+ s -× m - ε 2 + s -= 0, x ∈ (-1, 0), ∂ t m -= -m -× c( m -• u -) u --( m -• e x ) e x + s - ε + ν∂ 2 x m -+ α m -× ∂ t m -, (14)
       ε 2 ∂ t s + - 1 ε 2 ∂ 2 x s + + s + × m + ε 2 + s + = 0, x ∈ (0, 1), ∂ t m + = -m + × c( m + • u + ) u + -( m + • e x ) e x + s + ε + ν ε 2 ∂ 2 x m + + α m + × ∂ t m + , (15)
       ∂ x s -(t, -1) = j -1 (t), s -(t, 0 -) = s + (t, 0 + ), ε∂ x s -(t, 0 -) = ∂ x s + (t, 0 + ), ∂ x s + (t, 1) = 0, (16)
∂ x m -(t, -1) = ∂ x m -(t, 0) = 0, ∂ x m + (t, 0) = ∂ x m + (t, 1) = 0. ( 17
)
Moreover, it is easy to see that the Landau-Lifshitz equation keeps the norm of m(t, x) constant through time evolution. We assume that

| m -(0, x)| = 1 for all x ∈ [-1, 0], and | m + (0, x)| = 1 for all x ∈ [0, 1], (18)
so that this property stays true for all t, x.

Without going further in the analysis, one can notice immediately that ε orders are different in the two ferromagnetic materials. One can guess that, in the thick ferromagnetic material F -, s -× m -= 0 when ε → 0. So, the spin density is polarized in the direction of the magnetization, and we recover, for the magnetization, a non-coupled Landau-Lifshitz equation. On the contrary, in the thin ferromagnetic layer, the two quantities stay coupled in the limit ε → 0, and a more precise study is needed in order to guess the asymptotic behaviour.

A first description of the various scales

Looking at the system (14)-( 15)-( 16)-(17), we can easily see that there are various scales (both time and space scales). In this paragraph, we make a first description of these different scales. As we will see in the next section, the understanding of these different scales is important in order to construct an appropriate numerical scheme for the previous system (for small values of the parameter ε). A more precise asymptotic study will be performed in Section 4 in order to construct an approximate model for the full system ( 14)-( 15)-( 16)-(17) as ε → 0.

Time scales

There are three time scales for the spin density evolution: one for the behaviour of s in the direction parallel to the magnetization m, another for the behaviour of the part of s that is orthogonal to m and a last one for a time-boundary layer linked to the initial data. The first time scale is the scale t given in the injected current j -1 (t). In order to exhibit the two others in a simple framework, let us take the diffusive equation on the spin density with a given magnetic field m, in the thick layer F -. Assume moreover that m -is constant (independent of t and x, this will be the case if we consider the layer at equilibrium from the beginning, m -(t, x) = u -for all (t, x)) and that s is homogeneous in space. The spin density s -, solution of the ordinary differential equation

ε 2 ∂ t s -+ s -× m - ε 2 + s -= 0,
highlights two time-scales: a damping at scale t/ε 2 arising from the equation ε 2 ∂ t s -+ s -= 0 and fast oscillations at scale t/ε 4 due to L t ε 4 where L(τ ) is the group associated with the equation

∂ τ v + v × m -= 0.
Since m -does not depend on τ , the group L(τ ) can be explicitly computed (calculating separately the parallel and transverse components), so that

v(τ ) = v |τ =0 • m -m -+ sin (τ ) m -× v |τ =0 + cos (τ ) m -× v |τ =0 × m -. (19)
The decay at time scale t/ε 2 only affects the initial data of the problem, thus leading to a time-boundary layer. On the other hand, one can see in the expression of L that very fast oscillations at scale t/ε 4 appear for the part of the spin density that is orthogonal to m -. But, we remind here that the first thick layer should act as a polarizer, that is to say that s -should be collinear to m -. This property will be satisfied by the local mean in time with respect to the fast time variable. Indeed, the oscillating part has zero mean: when w does not depend on τ ,

1 2π 2π 0 L(τ ) w dτ = ( w • m -) m -. (20)
Since these very fast oscillations at scale t/ε 4 are not significant for the phenomenon we want to study here, in the sequel, in particular for numerical simulations, we will use techniques that allow us to avoid them, mainly filtering.

Space scales

In order to understand the various scales with respect to the x-variable in the thick layer, let us consider again a simple setting. We consider the stationary equation on the spin density, with a given constant magnetization m -:

-∂ 2 x σ -+ σ -× m - ε 2 + σ -= 0.
Now, we can write the equations on the components of σ -that are parallel or transverse to m -. For the transverse part, we write σ ⊥ = σ 1 w 1 + σ 2 w 2 with ( w 1 , w 2 ) an orthonormal basis of ( m -) ⊥ , and let

Z -= σ 1 + iσ 2 . Then, we get      -∂ 2 x ( σ -• m -) + ( σ -• m -) = 0, -∂ 2 x Z -+ - i ε 2 + 1 Z -= 0. ( 21
)
When solving the last equation, the roots

± - i ε 2 + 1 appear (where √ z is chosen such that Re( √ z) > 0). Since, - i ε 2 + 1 ≃ √ -i ε when ε → 0, we get Z -(x) ≃ Z -(0)e - √ -i x ε .
Hence, we can see that the orthogonal part of σ -depends on the variable x/ε. This term corresponds to a boundary layer near the interface x = 0 in the thick material, which we will describe precisely in Section 3. Thus, it will be important to take into account the two space scales in the thick layer: x and x/ε.

Numerical calculations

As we explained in the previous section, when the parameter ε is small, we have to handle two different space scales: x and x/ε, and three different time scales: t, t/ε 2 and t/ε 4 . For the space discretization, the two ideas are the following. First, to treat the existence of a boundary layer of size ε, we will use a non-uniform space mesh, refined near the interface in the thick material. Then, as we have described in the previous section, the behaviour of the spin density is quite different in the directions parallel or transverse to the magnetization. Hence, we will use a finite element method for the equation on s, using basis functions adapted to this phenomenon.

Concerning the time discretization, the time scale t/ε 4 is not relevant for the phenomena we are interested in (namely the magnetization switching). More precisely, since this scale only appears as oscillations with zero mean value on the transverse part of the spin density that are not observed in physical experiments, we do not want to simulate them. As a first step, in order to observe numerically the oscillations at the fast time scale, we will present calculations with a short final time, and a very small time step. Then, we will use an implicit time scheme with a bigger time step that will act as a filter with respect to these oscillations, and will allow us to recover only the non-oscillating part. With that implicit scheme, it is easy to make calculations for longer final times, and to observe the magnetization switching.

Construction of the numerical scheme

In this part, we describe in details the numerical scheme we have constructed for the system ( 14)-( 15)-( 16)-(17).

Notations. Let N be the spatial mesh element number and (x k ) 1≤k≤N +1 the coordinates of mesh points, with x 1 = -1 and x N +1 = 1. We denote ∆t the time step, and for all n ∈ N, t n = n∆t. We denote s n , and m n the discrete spin density and magnetization vectors obtained from s(t n , x k ) 1≤k≤N +1 and m(t n , x k ) 1≤k≤N +1 .

Let ( e 1 , e 2 , e 3 ) be the canonical basis in R 3 .

Let us now explain how we compute ( s n+1 , m n+1 ) from ( s n , m n ). We first compute s n+1 by solving the spin density equation with m n . Then, we compute m n+1 by solving the LL equation with s n+1 .

Time discretization of the diffusive equation on the spin density

As we will show in paragraph 3.2, using a very small time step, we will be able to observe the oscillations at scale t/ε 4 . But, we have already explained that actually we do not want to observe them, that is why we will filter them by using an implicit Euler scheme for the spin-density equation time-discretization.

It is easy to see how the filtering occurs with an implicit time scheme with a "large" time step, at least in a schematic framework. Let us consider the differential equation

u ′ (t) + λ u(t) × e 3 = 0, ( 22
)
where u : R → R 3 and λ is a real parameter. Then, the components of u satisfy

   u ′ 1 (t) + λu 2 (t) = 0, u ′ 2 (t) -λu 1 (t) = 0, u ′ 3 (t) = 0. ( 23
)
Using an implicit Euler scheme in time, denoting by u n the value calculated by the scheme at time t n = n∆t, we get for all n ≥ 0,

   u n+1 1 u n+1 2 = M u n 1 u n 2 , u n+1 3 = u n 3 , (24)
where

M = 1 1 + (λ∆t) 2 1 -λ∆t λ∆t 1 .
Now, when λ is very large and when ∆t is chosen such that λ∆t is still very large, we get that M behaves like

M ≃ 1 (λ∆t) 2 -1 λ∆t 1 λ∆t 1 (λ∆t) 2 .
Hence, each term u n 1 , u n 2 computed by the implicit scheme in that case is less than 1 λ∆t for n ≥ 1 : the oscillating part is filtered when λ∆t is large.

In the problem we consider here, the parallel with this schematic framework can be done with λ = 1/ε 4 (isolate the part ε 2 ∂ t s+ s× m ε 2 ). Hence, when choosing a time step such that ∆t/ε 4 is large, the oscillations at scale t/ε 4 are filtered when using an implicit scheme.

Space discretization of the diffusive equation on the spin density

Let n ∈ N be fixed, and assume that ( s n , m n ) is known. To determine s n+1 , we use an implicit Euler scheme for the time discretization (as described above) and a finite element method on the non-uniform space mesh (x k ) 1≤k≤N +1 for the space discretization. We first describe the basis functions, and then give the appropriate variational formulation.

Basis functions

For each time t n , each 1

≤ k ≤ N + 1, we construct a triplet of basis functions θ k (t n ) = ( θ 1 k , θ 2 k , θ 3 k )(t n ) supported in [x k-1 ; x k+1 ] ( [x 1 , x 2 ] for k = 1 and [x N , x N +1 ] for k = N + 1
). These basis functions are defined as being solutions of the following stationary problems (with suitable modifications for the cases k = 1 and k = N + 1): for 1

≤ j ≤ 3,                  -∂ 2 x θ j k + θ j k × m n k ε 2 + θ j k = 0, if x ∈] -1, 0[ -∂ 2 x θ j k + θ j k × m n k + ε 2 θ j k = 0, if x ∈]0, 1[ θ j k (x k ) = e j , θ j k (x k±1 ) = 0,
that is to say that θ j k (t n ) is the solution of two Dirichlet problems, the first one on [x k-1 , x k ], the second one on [x k , x k+1 ]. The equation satisfied by

θ j k (t n ) on (x k-1 , x k ) (respectively (x k , x k+1 )) is the first equation above if [x k-1 , x k ] ⊂ [-1, 0] (respectively [x k , x k+1 ] ⊂ [-1, 0]) and the second one above if [x k-1 , x k ] ⊂ [0, 1] (respectively [x k , x k+1 ] ⊂ [0, 1]
). These functions are thus piecewise exponential type functions (on [x k-1 , x k ] and [x k , x k+1 ]). They can be written as

θ j k (t n )(x) = T P ( m n k )R k (x)P ( m n k ) e j ,
where R k is a matrix that does not depend on t n , so that it is computed only once (and not at each step of the numerical calculation), and P ( m n k ) is the matrix of one orthonormal basis adapted to m n k with respect to the canonical basis. To obtain this formula, we separate the parallel and the transverse components compared to m n k (same procedure as in section 2.4). To illustrate the procedure and to give a more precise formula, let us solve the problem on

[x k-1 , x k ], if [x k-1 , x k ] ⊂ [0, 1].
The adaptation is straightforward for the other cases.

Let us denote by η

ε + = √ ε 2 -i = 1 √ 2 ε 2 + √ 1 + ε 4 -i -ε 2 + √ 1 + ε 4 the square root with positive real part and r k (x)e iφ k (x) = sinh η ε + (x -x k-1 ) sinh η ε + (x k -x k-1 )
.

Then the corresponding matrix R k is

R k (x) =     sinh ε(x-x k+1 ) sinh ε(x k -x k+1 ) 0 0 0 r k (x) cos φ k (x) -r k (x) sin φ k (x) 0 r k (x) sin φ k (x) r k (x) cos φ k (x)     .
Once the construction of these appropriate functions is done, the implementation of the finite element method is as usual. We first give the variational formulation we used, and then the related matrix system that is actually solved numerically.

Variational formulation

In order to implement the finite element method, we have to write a variational formulation of the equation on the spin density. Let us define for all (t, x),

s(t, x) = s -(t, x) for x ∈ (-1, 0), s + (t, x) for x ∈ (0, 1). ( 25
)
In order to obtain the variational problem satisfied by ω, we combine the equation on s -with the equation on s + multiplied by ε, together with the boundary and interface conditions on s -and s + . We thus get that s is defined as the solution in

H 1 (R + × [-1, 1], R 3 ) of the problem : for all θ ∈ H 1 (R + × [-1, 1], R 3 ), for all t > 0, ε 2 0 -1 ∂ t s • θ dx + ε 3 1 0 ∂ t s • θ dx + 0 -1 ∂ x s • ∂ x θ dx + 1 ε 1 0 ∂ x s • ∂ x θ dx + 1 ε 2 0 -1 s × m -• θ dx + 1 ε 1 0 s × m + • θ dx + 0 -1 s • θ dx + ε 1 0 s • θ dx = -j -1 (t) • θ(t, -1),
which we rewrite in the following concise form

∂ t s(t), (ε 2 1 (-1,0) + ε 3 1 (0,1) ) θ(t) L 2 x (-1,1) + a ε s(t), θ(t) = L ε θ(t) . ( 26
)
Matrix system Now, we write the discretized variational problem related to the previous variational formulation (26) and the basis ( θ n i ) 1≤i≤N +1 . This discretized problem is obtained in the following way. For the temporal discretization, we use an implicit Euler scheme, so that we write for n ∈ N,

s(t n+1 ) -s(t n ) ∆t , (ε 2 1 (-1,0) + ε 3 1 (0,1) ) θ(t n ) L 2 x (-1,1) + a ε s(t n+1 ), θ(t n ) = L ε θ(t n ) . (27)
Let n ∈ N be fixed. We decompose s(t n ) on the basis functions derived above:

s(t n , x) = N +1 k=1 s n k • θ k (t n )(x), (28)
where s n k ∈ R 3 for all k = 1 . . . N + 1. Next, we construct the vector sn as being the column vector made up from the vectors s n k for 1 ≤ k ≤ N +1. The contributions of the two first terms of the variational formulation (containing the time derivatives) are inserted in a matrix B ε . The matrix A ε contains the terms corresponding to a ε . Finally, the L ε term (which allows to treat the non homogeneous Neumann condition) is in a vector L ε . Finally, we write an implicit Euler scheme for the time discretization. Thus, the problem consists to solve the following matrix system

A ε + ε 2 ∆t B ε sn+1 = L ε + ε 2 ∆t B ε sn . (29)
We would like to indicate that the coefficients of the matrices A ε and B ε are of the type θ j k (t n )θ j ′ k ′ (t n )dx (or with derivatives). These integrals can be calculated analytically using the explicit formulas we obtained for the basis functions. We use these exact values in our numerical scheme. This method of finite element with exponential basis functions allows to model quite accurately the problem with an acceptable space step. For example, taking ε = 1.10 -2 , we need a space step of 5.10 -4 with the finite differences method to take in account correctly the limit layer of thickness ε. With the finite element method, we just need to use a refined mesh in the asymptotic region. The number of point is divided by 10 and the time computation divided by 5 to obtain the same accuracy.

Landau-Lifshitz equation discretization

For the Landau-Lifshitz equation, we use the scheme proposed by Wang, Garcia-Cervera, E [START_REF] Wang | A Gauss-Seidel projection method for micromagnetics simulations[END_REF]. We do not rewrite it here in full details. Let us just say that it is an implicit Gauss-Seidel projection scheme. It is a fractional step method adapted to the LL equation ∂ t m = -m × H + α∂ t m × m, that can be rewritten (at least for smooth solutions)

(1 + α 2 )∂ t m = -m × H -αm × m × H .
The first step consists of an implicit Gauss-Seidel method. The second step corresponds to performing the heat flow, the third and last step is just the projection on S 2 (in order to keep the norm of m constant). Let us just mention that a priori the effective field in our coupled LL equation contains singular terms:

m -× s - ε in the thick layer, m + × s + ε and 1 ε 2 ∂ 2
x m + in the thin layer. These terms do not have to be treated in a specific way. Indeed, even if they seem singular, the two terms in the thin layer are actually of order O(1) when ε → 0 (as shown in the asymptotic study in Section 4). In the thick layer, the term m -× s - ε contains only a very fast oscillating part that is filtered by the implicit scheme used for the equation on the spin density. So, a posteriori nothing has to be done with respect to these terms.

Numerical observation of the very fast time oscillations

First, we present here calculations with an explicit time scheme with a very small time step. Since the calculations are quite time demanding, we chose a not too small value of ε and a short final time: ε = 10 -1 , T = 2.10 -2 . We took a time step δt = 10 -8 , a space step of order h = 10 -2 (since there is a CFL type condition δt ≤ ε 4 h 2 ). It gives Figure 2 in which the torque s × m is presented in logarithmic scale (right) and for a given x inside F -(left). We clearly observe the fast oscillations at scale t/ε 4 and the decay at scale t/ε 2 .

Figure 2: The oscillations on the spin density in the thick layer at scale t/ε 4 .

Numerical simulations of the magnetization switching

In this part, we present our numerical simulations obtained with the scheme we have described in Section 3.1 (using the implicit time discretization). In particular, we show the magnetization switching in the thin layer, and we study the impact of the intensity of the injected current on this switching phenomenon. As we will see, our numerical simulations are in accordance with physical experiments. As initial conditions, we assume that the thick ferromagnetic material F -is near the equilibrium and that the thin ferromagnetic material F + is at equilibrium before injecting the spin current, that is to say that we assume that for all x (respectively in (-1, 0) and (0, 1)),

u -=   0 -sin θ cos θ   , m -(t = 0, x) =   0 -sin θ ′ cos θ ′   and m + (t = 0, x) = u + =   0 0 1   , ( 30
)
where θ = 30 • is the angle between the two anisotropy directions u -and u + and θ ′ = 40 • i.e. the angle between the initial condition for m -and the anisotropy direction u -is 10 • . As we begin the injection at t = 0 + , we assume that s(t = 0) = 0.

The injected current is chosen in the form

j -1 = j -1   0 -sin ϕ cos ϕ   with ϕ = 10 • .
The following numerical simulations are done with a parameter ε = 10 -2 . The time step is ∆t = 10 -2 , and the space step is 10 -2 away from the interface, and 5 • 10 -4 in the boundary layers of the thick layer.

Observation of the magnetization switching

To start with numerical results, we present the evolution of the magnetization in a case where we observe a switching (we inject an intense enough current here, we will discuss that point in the next paragraph). In Figures 3 and4, one can see the evolution of the magnetization components with respect to time, both in the thick and in the thin layers (the left part between -1 and 0 describes F -and the right part between 0 and 1 describes F + ). On one hand, as expected, we observe that the magnetization in the thick ferromagnet F -remains at equilibrium. On the other hand, we see that the magnetization moves in F + under the spin current influence, and that around time 6, it is almost completely switched. It is obvious that the end point reached by m + is not the equilibrium point -u + , that would correspond to m + z = -1. We identified this point as being -u -: that is to say that the magnetization in the thin material is almost reversed, but the end point point is not in the anisotropy direction of the F + ferromagnet but in the anisotropy direction of the thick ferromagnet F -. This point is confirmed by Figure 5. The reader can notice on Figures 3 and4 that the magnetization m + is constant with respect to the space variable through the thin layer. We will recover this result in our asymptotic study.

Study of the torque s × m

It is also interesting to study the torque s × m, since it is responsible of the magnetization switching. We remind the switching is due to interactions between the transverse component of the spin density and the magnetization in the thin material. The torque is presented in Figure 6. One can notice that it is indeed important in F + when the switching starts, but | s × m| is never greater than O(ε).

Current intensity influence

In physical experiments, it is known that there is a threshold for the intensity of the injected current above which it is possible to observe the magnetization reversal. In Figure 7, we plotted the results we obtain with our model about the impact of the injected current. We recover the existence of such

Current sign influence

The second experiment is to take a current whose intensity allows a switching and to change its sign during the time. We obtain Figure 8. It represents the evolution of the component m + z during the time. In this case, we start with a negative sign (such that j -1 • u -< 0). This sign is necessary if we want to observe a magnetization switching in F + . At time t = 10, the magnetization has converged to -u - and we decide to invert the current sign. Then, the switching is reversed and the magnetization goes up to u -. At time 20, we change the sign again coming back to the initial current. The magnetization is switched again and returns to -u -. At the times 30 and 40, we make the same manipulations. The magnetization follows exactly the same way around the unit sphere as previously. Finally, at time 50, we stop the current taking j -1 = 0. The magnetization converges slowly to the inferior pole of the sphere -u + . We recall here that ±u + are known to be stable equilibrium states for the non coupled Landau-Lifshitz equation

∂ t m + = -m + × c( m + • u + ) u + -( m • e x ) e x + ν ε 2 ∂ 2 x m + + α m + × ∂ t m + .
This succession of switchings allows to see the role of the current sign. It shows that the magnetization switching is reversible. We can choose to guide the magnetization up to the superior part or down to the inferior part of the sphere.

Derivation of an asymptotic model and numerical comparison

In this section, we are interested in the derivation of an approximate model for the full system ( 14)-( 15)-( 16)-(17) when ε goes to 0. Our derivation here is purely formal, and the validation is done numerically, by presenting numerical comparisons between our asymptoic model and the full system. The mathematical study of existence and convergence will be addressed in a following paper by E. Fouassier and D. Sanchez.

A first model

The first idea is to obtained an approximate model that does not depend on ε anymore. In order to find such a model, we use a multiscale approach. We introduce the various scales appearing in the phenomena that we described in the previous section. We recall that there were three time scales, t, t/ε 2 and t/ε 4 , and two space scales in the thick material, x and x/ε. We have already mentioned that the oscillations at the fast time scale t/ε 4 are not significant here, so we will omit that scale in the asymptotic study. Indeed, these oscillations have zero mean, so if we do not take into account this scale in our study, we recover the behaviour of the local mean in time of s with respect of this fast time scale, which is what we are interested in here (one can have a look at the formula (19) obtained when the three scales are distinguished). Moreover we assume that the initial data is well prepared and that we are looking for the steady-state model, thus avoiding the time scale t/ε 2 .

Because of the boundary layer near the interface in the thick material, we decompose the solutions in this material as

s ε-(t, x) = s ε -(t, x) + s ε -t, -x ε , m ε-(t, x) = m ε -(t, x) + m ε -t, -x ε .
Then, the first idea is to make the following ansatz for U = s ε -, m ε -, s ε+ or m ε+ , and U = s ε -or m ε -: we expand U and U as a sum of powers of ε with profiles that are independent of ε, in the form:

U (t, x) = U 0 (t, x) + εU 1 (t, x) + ε 2 U 2 (t, x) + . . . U (t, z) = U 0 (t, z) + ε U 1 (t, z) + ε 2 U 2 (t, z) + . . . (31)
In these ansatz, we assume that the profiles U j (t, z) together with all their derivatives go to 0 when z → +∞: they describe a boundary layer at the right end of the thick ferromagnetic layer (near the interface x = 0). To simplify the notations we let s - j = s - j + s - j and m - j = m - j + m - j for all j. Our aim is then to obtain the equations fulfilled by the first orders profiles. To do this, we plug the asymptotic expansion into the equations ( 14)-( 15)-( 16)-( 17), and identify powers of ε. Each equation can then be rewritten in the following form j≥-2

ε j P j (s ± k , s ± k , m ± k , m ± k ) = 0,
where P j 's are operators acting on some profiles

s ± k , s ± k , m ± k , m ± k (with several indices k).
From this, as usual, we now identify P j (s ± k , s± k , m± k , m± k ) = 0 for all j ≥ -2. In order to separate the boundary layer profiles U j from the other part U j , we take the limit z → +∞ in the previous equations. Thus, we first obtain equations on U j profiles, and we then make the difference between the first equations and the last ones to obtain equations on the boundary layer part.

Finally, we plug the ansatz into the boundary and interface conditions and into the condition | m ± | 2 = 1, and we again identify powers of ε. For instance, we get |m 0 | 2 + 2εm 0 • m 1 + . . . = 1, from which we deduce |m 0 | 2 = 1, m 0 • m 1 = 0 and so on.

To calculate the first profiles, and obtain an asymptotic model, we need the equations satisfied at orders ε -2 , ε -1 , and ε 0 .

Equations in the thick ferromagnet

At order ε -2 in the two equations in the thick material, we first obtain the following equations. First, s - 0 × m - 0 = 0, so there exists a profile a 0 such that

s - 0 (t, x) = a 0 (t, x)m - 0 (t, x). (32)
Second, we get

m - 0 × ∂ 2 z m - 0 = 0. (33) Morerover, since |m - 0 | 2 = 1, we have m - 0 • ∂ z m - 0 = 0.
Using this property and (33), we write

∂ 2 z m - 0 = ∂ 2 z s - 0 • m - 0 m - 0 .
We then multiply by ∂ z m - 0 and we obtain ∂ 2 z m - 0 • ∂ z m - 0 = 0. Now, we integrate twice and we use that m - 0 , ∂ z m - 0 → 0 as z → ∞, to conclude that m - 0 = 0. Third, we obtain the equation fulfilled by the boundary layer profile s -

0 -∂ 2 z s - 0 + s - 0 × m - 0 = 0. (34)
Since, m - 0 does not depend on z, we can explicitely solve this equation (as it has been done before). The property that we will use in the sequel on s-0 , is that this equation allows us to write a Dirichlet to Neumann relation on s-0 ∂ z s-0 (t, 0 -) = Ds - 0 (t, 0 -), where D is a positive matrix.

At order ε -1 , we then obtain

s - 0 × m - 1 + s - 1 × m - 0 = 0, (36)
and

m - 0 × ∂ 2 z m - 1 = 0. ( 37
)
As above, we combine (37) with m - 0 • m - 1 = 0 and it gives m - 1 = 0. We also have

-∂ 2 z s - 1 + s - 1 × m - 0 + s - 0 × m - 1 = 0. ( 38
)
At order ε 0 , we get

-∂ 2 x s - 0 + s - 0 × m - 2 + s - 1 × m - 1 + s - 2 × m - 0 + s - 0 = 0, (39)
which gives, taking the scalar product with m - 0 , and using (36)

-∂ 2 x a 0 + (1 + |∂ x m - 0 | 2 )a 0 = 0. ( 40
)
Here again, we used that the norm of the magnetization m - 0 is preserved:

|m - 0 | 2 = 1 so that m - 0 •∂ x m - 0 = 0 and m - 0 • ∂ 2 x m - 0 = -|∂ x m - 0 | 2 .
We then write the equation on m - 0 . Using (36), we get the non coupled LL equation satisfied by m -

0 ∂ t m - 0 = -m - 0 × c(m - 0 • u -)u --(m - 0 • e x )e x + ν∂ 2 x m - 0 + αm - 0 × ∂ t m - 0 . (41)

Equations in the thin ferromagnet

At order ε -2 , we first obtain the following equations

-∂ 2 x s + 0 + s + 0 × m + 0 = 0, (42)
and

m + 0 × ∂ 2 x m + 0 = 0. ( 43
)
Combined with |m + 0 | 2 = 1 and the Neumann conditions (52), equation (43) gives that m + 0 (t, x) is independent of the space variable x.

At order ε -1 , we then obtain

-∂ 2 x s + 1 + s + 1 × m + 0 + s + 0 × m + 1 = 0, (44)
-m + 0 × (ν∂ 2 x m + 1 + s + 0 ) = 0. ( 45
)
At order ε 0 , we get

∂ t m + 0 = -m + 0 × c(m + 0 • u + )u + -(m + 0 • e x )e x + ν∂ 2 x m + 2 -(m + 0 × s + 1 + m + 1 × s + 0 ) + αm + 0 × ∂ t m + 0 .
Since m + 0 is independent of x, an integration over x ∈ [0, 1] of the previous equation gives

∂ t m + 0 = -m + 0 × c(m + 0 • u + )u + -(m + 0 • e x )e x + αm + 0 × ∂ t m + 0 -m + 0 × ν∂ x m + 2 (t, 1) -ν∂ x m + 2 (t, 0) + ∂ x s + 1 (t, 1) -∂ x s + 1(t, 0) .
Let γ(t, x) = m - 0 (t, 0 -) • m + 0 (t, x). The vector field m + 0 then fulfills:

∂ t m + 0 = -m + 0 × C(m + 0 • u + )u + -(m + 0 • e x )e x + ∂ x a 0 (t, 0 -) 1 -γ 2 -(pγ + q)m - 0 + m - 0 × m + 0 +αm + 0 × ∂ t m + 0 .
Since |m + 0 (t)| 2 = 1 for all t, as it is usual when dealing with LL equations, this last equation can be rewritten as

(1 + α 2 )∂ t m + 0 = -m + 0 × C(m + 0 • u + u + ) -(m + 0 • e x )e x -∂ x a 0 (t, 0) pγ(t) + q -α 1 -γ 2 (t) m - 0 (t, 0) -m + 0 × m + 0 × αC(m + 0 • u + u + ) -α(m + 0 • e x )e x -∂ x a 0 (t, 0) αpγ(t) + αq + 1 1 -γ 2 (t) m - 0 (t, 0) .( 55
)
where a 0 satisfies (54).

Let us explain how to see on this limit equation that the magnetization m + 0 will be switched. First, let us mention that when there is no spin-current (a 0 (t, -1) = 0 which implies a 0 (t, x) = 0 for all (t, x)), the magnetization m + 0 satisfies a non coupled Landau-Lifshitz equation, and it classically converges to the nearest stable equilibrium state (±u + ). On the contrary, when a 0 (t, -1) = f (t) = 0, we can study the behaviour and the switching of m + 0 thanks to the quantity γ(t) = m - 0 (t, 0) • m + 0 (t), that is a solution of

(1 + α 2 )∂ t γ(t) = -m + 0 × C(m + 0 • u + )u + -(m + 0 • e x )e x + αm + 0 × C(m + 0 • u + )u + -(m + 0 • e x )e x • m - 0 (t, 0) +m + 0 • ∂ t m - 0 (t, 0) + m + 0 × m + 0 × ∂ x a 0 (t, 0) αpγ(t) + αq + 1 1 -γ 2 (t) m - 0 (t, 0) • m - 0 (t, 0) = F (m - 0 , m + 0 ) -∂ x a 0 (t, 0)(αpγ + αq + 1),
where F (m - 0 , m + 0 ) is bounded. In this equation, ∂ x a 0 (t, 0) writes -f (t)u(t) where u(t) is bounded, greater than a strictly positive constant and only depends on m - 0 . Thus, for a 0 (t, 0) = f (t) large enough, γ(t) converges in finite time to ±1, which means that m + 0 converges to ±m - 0 . Here again, we see that our limit model is valid only up to convergence, when m + 0 and m - 0 are not collinear. Indeed, this limit model has a singularity when γ(t) converges to ±1.

Numerical comparison between the full system and the first asymptotic expansion

In this part, we compare the asymptotic expansion obtained in Section 4.1 (when m + and m -are never collinear) with the general system. All the following calculations are done with ε = 1.10 -2 . The errors are calculated in the L 2 norm with respect to x at each instant.

Discretization of the limit system

We will not give many details on the discretization of the limit system, since we use a simple finite difference method. Let us just say that for the LL equation (55) satisfied by m + 0 (that does not contain any laplacian term anymore), we use the scheme proposed by Joly and Vacus [6]. This scheme has the nice property to keep the norm of the magnetization constant during time evolution, as the LL equation does. It is written as follows, m is approximated by m n+1 + m n 2 :

m n+1 -m n ∆t = - m n + m n+1 2 × H e ( m n ) + α m n + m n+1 2 × m n+1 -m n ∆t , (56)
where H e is the total effective field.

Numerical comparison

As we expected, our limit model is a good approximation up to order O(ε) only up to the magnetization switching in the thin material. The three following figures show the L 2

x error between the solutions of the full system and the solutions of the limit model. When the magnetization is switched in the thin material, the errors between the magnetization of the full system m ε,+ and m + 0 and between s ε,-and s - 0 grow up to the order 0.1 whereas the small parameter ε equals 10 -2 . on the spin density increases even more and never decreases again (Figure 10). This error comes from the behaviour of s -: in the full system, the norm of s -increases brutally when the magnetization m + and m -become collinear (Figure 11) whereas in the asymptotic model the norm of the spin density s - 0 does not vary.

Modified asymptotic expansion

The previous asymptotic model is not satisfying from two points of view. First, it is valid only up to convergence. Second, due to the behaviour of the solution of the limit ordinary differential equation satisfied by the magnetization m + 0 , the limit model is not reversible, in the sense that if we change the sign of the injected current the magnetization will not come back to the upper position. In order to obtain a limit model with a bigger domain of validity, we use a slightly different approach, The idea is of keeping a dependence on the parameter ε for the profiles in our expansion to better track the change of behaviour when m + 0 and m - 0 are almost collinear. Indeed when 1 -γ 2 (t) is of order ε a part of the term m + 1 is of order 1/ε and not taken into account in the expression of m + 0 thus leading to the loss of accuracy of the previous asymptotic expansion.

Rather than using the previous ansatz (formal)

U ε (t, x) = ε k U k t, x, x ε , we use U ε (t, x) = ε 2k U ε k t, x, x ε .
The new profiles now depend on ε and the approximate model (i.e. the equations on the new first profiles U ε 0 ) we will get will still depend on ε. Our approach could be described by two (equivalent) ways. We look for an approximate model that is a good approximation of the full model up to order O(ε), when ε → 0. The idea is of putting together the profiles with indices 0 and 1 obtained previously. One possible approach to do this is to write the equations satisfied by s 0 + εs 1 , and so on . . ., and to truncate them at order O(ε). Or, equivalently, we can repeat the previous procedure with an Ansatz containing only even powers of ε, and keep the orders ε 2k and ε 2k-1 together to obtain the equations fulfilled by the profiles (now depending on ε). In the sequel, we present our results using the second approach.

The second model

We look for a formal limit model approximating the equations up to first order. As before, since we expect the formation of a boundary layer in the thick ferromagnet near x = 0 -, we make the following ansatz:

s ε-(t, x) = s ε-(t, x) + s ε-t, -x ε , m ε-(t, x) = m ε-(t, x) + m ε-t, -x ε .
With these notations, the boundary conditions write:

• ∂ x s ε-(t, -1) = j -1 (t) (given, independent of ε), • ∂ x s ε+ (t, 1) = 0, • ∂ x m ε-(t, -1) = 0, ∂ x m ε-(t, 0) + 1 ε ∂ z m ε-(t, 0) = 0, • ∂ x m ε+ (t, 0) = ∂ x m ε+ (t, 1) = 0.
and the transmission conditions now write:

• s ε-(t, 0) + s ε-(t, 0) = s ε+ (t, 0), • - 1 ε ∂ z s ε-(t, 0) + ∂ x s ε-(t, 0) = 1 ε ∂ x s ε+ (t, 0).

Equations in the thick ferromagnet

By taking the limit as z goes to +∞ in the equation in the thick ferromagnet we obtain:

ε 2 ∂ t s ε--∂ 2 x s ε-+ s ε-× m ε- ε 2 + s ε-= 0, (57)
∂ t m ε-= -m ε-× C(m ε-• u -)u --(m ε-• e x )e x + s ε- ε + ν∂ 2 x m ε-+ αm ε-× ∂ t m ε-. (58)
To obtain an approximate solution exact up to order 1 in ε we perform an asymptotic expansion of the equation in ε 2 by keeping the orders ε 2k and ε 2k-1 together for k ≥ -1. Let us denote by s ε- 0 , m ε- 0 , . . ., the first profiles we are looking for. In this way we do not have to perform further approximations on the Landau-Lifshitz equation: m ε- 0 satisfies (58). Then, thanks to this policy the equation on s ε-gives:

s ε- 0 × m ε- 0 = 0 and -∂ 2 x s ε- 0 + s ε- 0 = 0.
In the same way, we now obtain equations on m ε- 0 and s ε- 0 : keeping the terms of order -2 and -3 in ε we have:

m ε- 0 × ∂ 2 z m ε- 0 = (m ε- 0 + m ε- 0 ) × ∂ 2 z m ε- 0 = 0, -∂ 2 z s ε- 0 + s ε- 0 × (m ε- 0 + m ε- 0 ) = 0. Since m ε-is a solution to the Landau-Lifschitz equation, |m ε-| = 1 which implies that ∂ 2 z m ε- 0 = ∂ 2 z m ε- 0 = 0 and then m ε- 0 = 0.
This implies that m ε- 0 fulfills the classical Landau-Lifschitz equation in [-1, 0] with Neumann boundary conditions in {-1, 0}:

   ∂ t m ε- 0 = -m ε- 0 × C(m ε- 0 • u -)u --(m ε- 0 • e x )e x + ν∂ 2 x m ε- 0 + αm ε- 0 × ∂ t m ε- 0 , m ε- 0 (0, x) given, ∂ x m ε- 0 (t, -1) = ∂ x m ε- 0 (t, 0) = 0.
We also have the existence of a function a -: R +

t × [-1, 0] → R such that s ε- 0 (t, x) = a -(t, x)m ε- 0 (t, x).
By taking the scalar product with m ε- 0 of the equation on s ε- 0 we obtain the following equation on a

-: -∂ 2 x a -+ (1 + |∂ x m ε- 0 | 2 )a -= 0, x ∈ [-1, 0], ∂ x a -(t, -1) = j -1 (t) • m ε- 0 (t, -1). The equation on s ε- 0 becomes: -∂ 2 z s ε- 0 + s ε- 0 × m ε- 0 = 0. Solving this equation, we get s ε- 0 (t, z) = β 1 (z)m ε- 0 (t, x) × s ε- 0 (t, 0) × m ε- 0 (t, x) -β 2 (z) s ε- 0 (t, 0) × m ε- 0 (t, x), where (β 1 + iβ 2 )(z) = e - √ -iz with √ -i = e -iπ/4 .

Equations in the thin ferromagnet

We study in the same way the equations in the thin ferromagnet by performing the same asymptotic expansion. We get at orders ε -3 and ε -2 :

m ε+ 0 × ∂ 2 x m ε+ 0 = 0, -∂ 2
x s ε+ 0 + s ε+ 0 × m ε+ 0 = 0. As before, this implies that m ε+ 0 is independent from x and we have an explicit formula for s ε+ 0 :

s ε+ 0 (x) = (s ε+ 0 (0) • m ε+ 0 )m ε+ 0 + β ′ 1 (x)s ε+ 0 (0) × m ε+ 0 + β ′ 2 (x)m ε+ 0 × s ε+ 0 (0) × m ε+ 0 , where (β ′ 1 + iβ ′ 2 )(x) = i cosh( √ -i(x -1)) cosh( √ -i) .
At orders ε -1 and ε 0 we obtain: .

We now only need to determine a condition on a -(t, 0) and ∂ x a -(t, 0) to solve the problem. We compute the scalar product of the transmission condition with m ε- 0 (t, 0) and we obtain a -(t, 0) = s ε+ 0 (t, 0) • m ε- 0 (t, 0) which gives the following Dirichlet-to-Neumann condition:

∂ x a -(t, 0) = - a -(t, 0) ε a 2 + b 2 + a √ 2 1 -γ(t) 2 (a 2 + b 2 ) √ 2 + a γ(t) 2 -2bγ(t) + a + √ 2 .
As one can easily see, this model is a generalization of the model we obtained thanks to the first asymptotic expansion. If we perform the extra assumption that M ε+ 0 and M ε- 0 are never collinear we come down to the previous asymptotic model.

Numerical comparison with the modified model

Comparison

In this part, we compare the modified asymptotic expansion obtained in Section 4.3 with the general system. We use the same parameters as in the previous simulations (ε = 1.10 -2 ). With this new asymptotic model we obtain a better qualitative behaviour of the magnetization m ε+ 0 in the thin ferromagnet F + (no more spinning around the expected limit) and the L 2 -error always remains (even in the worst case, when m ε- 0 and m ε+ 0 line up) of order at most ε 2 (Figure 12). Moreover, Figure 13 shows that this better behaviour of the model also appears on the spin density s ε- 0 with a better L 2 error and a better qualitative behaviour of s ε- 0 .

Computation time

To compute the general problem, we have to solve a large linear system at each time iteration and we have to take a refined space step in regions with an asymptotic behaviour (specially when ε is small). Consequently, the computation time is important. With the asymptotic expansion, we only solve the first order. It saves a considerable amount of time. For example, on the same machine with ε = 1.10 -2 , the simulation of first term of the asymptotic expansion is 250 times faster than the simulation of the general problem.

Conclusion

In this paper, we studied the model of magnetization reversal proposed by [START_REF] Zhang | Mechanisms of spin-polarized current-driven magnetization switching[END_REF], in a one dimensional framework. We emphasized the different time and space scales that exist in the model and we constructed a numerical scheme that is appropriate to treat these scales. This scheme allowed us to recover the results observed in the experiments lead by the physicists on the magnetic switching thanks to a spin current.

The second part of our work consisted in the derivation of an asymptotic model as ε → 0. We ended up with a limit model that is very easy to handle numerically and is good approximation of the original model (up to O(ε 2 )) . It could be a good alternative to the original model when dealing with very small values of ε.

As we mentioned in the introduction, spin transfer appears to be a turning point in spintronics and is the subject of an extensive research for applications. The spin transfer torque seems to be involved in several different phenomena: the magnetization switching, that we studied here, magnetization excitations (or magnetic oscillators), or domain wall motion. The very fast oscillations we pointed out in this paper could appear to be a key point in the magnetic oscillators applications, and the study of that phenomenon, both theoretically and numerically, will be our next target.

Aknowledgements. Elise Fouassier, Clément Jourdana and David Sanchez would like to thank Naoufel Ben Abdallah who proposed them to work on this subject. They would have liked him to read this completed work. They also want to thank Gilles Carbou for fruitful discussions on their work. This work was partially supported by the ANR project QUATRAIN (QUAntum TRAnsport In Nanostructures).

Figure 1 :

 1 Figure 1: A multilayered ferromagnetic device proposed by Slonczewski.

Figure 3 :

 3 Figure 3: Evolution of the m z component during the time (3D view and projection).

Figure 4 :

 4 Figure 4: Evolution of the m x and m y components.

Figure 5 :

 5 Figure 5: Evolution of the magnetization m + during the time around the unit sphere (the starting point is the north pole, and the arrow points -u -) (left), evolution of m + . u -during the time (right).

Figure 6 :

 6 Figure 6: Evolution of the torque s × m 2 during the time: 3D view (left), and with a logarithmic scale in the thin material (right).

Figure 7 :

 7 Figure 7: Evolution of the component m + z for different values of the injected current intensity j -1 .

Figure 8 :

 8 Figure 8: A succesion of switchings during the time.

Figure 9 :

 9 Figure 9: Comparison of m + and m + 0 .

Figure 10 :

 10 Figure 10: Comparison of m -and m - 0 (left), comparison of s -and s - 0 (right).The problems observed on the magnetization m + are linked to the variations of the spin density s - 0 (the asymptotic expansion gave s + 0 = 0) and we find when m + 0 goes to m - 0 that the already big L 2 -error

∂ t m ε+ 0 = 0 ε + αm ε+ 0 × ∂ t m ε+ 0 .∂ t m ε+ 0 = 4 . 3 . 4

 0000434 -m ε+ 0 × C(m ε+ 0 • u + )u + -(m ε+ 0 • e x )e x + s ε+Since m ε+ 0 does not depend on x ∈ [0, 1], we perform an integration over x ∈ (0, 1):-m ε+ 0 × C(m ε+ 0 • u + )u + -(m ε+ 0 • e x )e x dx +αm ε+ 0 × ∂ t m ε+ 0 .Determination of s ε+ 0 (x = 0) We now use the transmission conditions at x = 0. Here is the only difference with the first asymptotic model we obtained. They can be written now as

2 2 2 + 1 - 2 + a √ 2 2 ( 1 -

 221221 s ε- 0 (t, 0) + a -(t, 0)m ε- 0 (t, 0) = s ε+ 0 (t, 0), -1 ε ∂ z s ε- 0 (t, 0) + ∂ x a -(t, 0)m ε- 0 (t, 0) = 1 ε ∂ x s ε+ 0 (t, 0).Thanks to the previous results we obtained on s ε- 0 and s ε+ 0 , we gets ε+ 0 (t, 0) = √ a -(t, 0) -ε∂ x a -(t, 0) a 2 + b 2 + a √ bγ(t, 0) √ γ(t, 0) 2 ) × -b + γ(t, 0)(a 2 √ 2 + b 2 √ 2 + a) m ε+ 0 (t, 0) + -bγ(t, 0) + a + √ 2 m ε- 0 (t, 0) + (aγ(t, 0) + b) m ε- 0 (t, 0) × m ε+ 0 (t, 0) ,with a + ib = √ -i tanh √ -i and γ(t, x) = m ε+ 0 (t, x) • m ε-(t, 0).

4. 3 . 5 2 2 2 + 1 -bγ √ 2 + a √ 2 2 ( 1 -

 35221221 Equation on m ε+ 0 and Dirichlet to Neumann condition on a - We finally obtain an ordinary differential equation that rules the behaviour of m ε+ 0 :∂ t m ε+ 0 = -m ε+ 0 × C(m ε+ 0 • u + )u + -(m ε+ 0 • e x )e x + D (a 2 + b 2 )γ -b √ 2 m ε- 0 (t, 0) -D(a 2 + b 2 + a √ 2)m ε- 0 (t, 0) × m ε+ + αm ε+ 0 × ∂ t m ε+ 0 , a -(t, 0) -ε∂ x a -(t, 0) a 2 + b 2 + a √ γ 2 )

Figure 12 :

 12 Figure 12: Comparison between m + and m ε+ 0 .

Figure 13 :

 13 Figure 13: Comparison between s -and s ε-.

  The model proposed by Zhang, Levy and Fert . . . . . . . . . . . . . . . . . . . . . . . . . 2.1.1 The equations in one ferromagnetic layer . . . . . . . . . . . . . . . . . . . . . . . 2.1.2 Boundary and interface conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The scaled model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 A first description of the various scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

	1 Introduction
	2 The spin transfer model
	2.1

2.3.1 Time scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Space scales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Boundary and interface conditions

We will use (or we have already used) the following conditions with indices i = 0 and i = 1:

Solving the equations on the spin density

We first use the equation on s + 0 and the Dirichlet to Neumann equation ( 35) to write the following energy estimate

Since D is a positive matrix, we deduce the two following properties

• s-0 (t, 0 -) = 0, so that s-0 (t, z) = 0 for all t, z (see equation ( 34)),

• ∂ x s + 0 = 0, so that s + 0 is independent of x. Using the transmission condition at x = 0 (48), and we write s + 0 (t) = s - 0 (t, 0 -) = a 0 (t, 0 -) m-0 (t, 0 -). We plug this into the equation on s + 0 , and we get

From this result, we can deduce all the spin density profiles at order ε 0 , while the condition m-0 (t, 0 -)× m + 0 (t) = 0 is satisfied:

• there is no boundary layer at order ε 0 : s-0 (t, z) = 0,

• s-0 is polarized along m-0 : s-0 = a 0 m-0 where a 0 is entirely determined by the equation

• the spin density in F + satisfies s + 0 = 0.

Equation on m + 0

To obtain the equation fulfilled by m + 0 (t) we also need to compute the profile s + 1 (at least its values at the boundaries). As in the previous paragraph, this can be done using the equations fulfilled by the three profiles s - 1 , s - 1 and s + 1 and the transmission and boundary conditions that link them. After these calculations, we get the equation satisfied by the magnetization m + 0 . In order to write it, we need the following notations. We let A + iB = e -iπ/4 tanh e -iπ/4 and