A Derivation of linearized LLG equation

# Instability analysis of spin-torque oscillator with an in-plane magnetized free layer and a perpendicularly magnetized pinned layer

## Abstract

We study the theoretical conditions to excite a stable self-oscillation in a spin-torque oscillator with an in-plane magnetized free layer and a perpendicularly magnetized pinned layer in the presence of magnetic field pointing in an arbitrary direction. The linearized Landau-Lifshitz-Gilbert (LLG) equation is found to be inapplicable to evaluate the threshold between the stable and self-oscillation states because the critical current density estimated from the linearized equation is considerably larger than that found in the numerical simulation. We derive a theoretical formula of the threshold current density by focusing on the energy gain of the magnetization from the spin torque during a time shorter than a precession period. A good agreement between the derived formula and the numerical simulation is obtained. The condition to stabilize the out-of-plane self-oscillation above the threshold is also discussed.

###### pacs:
75.78.Jp, 75.76.+j, 85.75.-d

## I Introduction

A spin polarized current injected into a nanostructured ferromagnet creates spin torque through the spin-transfer effect [Slonczewski, 1996; Berger, 1996; Slonczewski, 2005]. The spin torque provides a rich variety of magnetization dynamics such as switching or self-oscillation [Katine et al., 2000; Kiselev et al., 2003; Rippard et al., 2004; Kubota et al., 2005; Krivorotov et al., 2005; Kubota et al., 2013; Tamaru et al., 2014]. In particular, a spin-torque oscillator consisting of an in-plane magnetized free layer and a perpendicularly magnetized pinned layer has been an attractive research subject in the field of magnetism. [Kent et al., 2004; Lee et al., 2005; Zhu and Zhu, 2006; Houssameddine et al., 2007; Firastrau et al., 2007; Ebels et al., 2008; Silva and Keller, 2010; Suto et al., 2012; Lacoste et al., 2013; Kudo et al., 2014; Bosu et al., 2016]. In this type of spin-torque oscillator, the spin torque forces the magnetization of the free layer into out of plane, and excites a large amplitude oscillation around the perpendicular axis. A high symmetry along the perpendicular direction in this system makes it easy to investigate the oscillation properties theoretically [Silva and Keller, 2010]. In order to observe the oscillation experimentally through magnetoresistance effect, however, the symmetry breaking should occur since the change of the relative angle between the magnetizations of the free and pinned layers in time is necessary. The linear analysis in the presence of an in-plane anisotropy [Ebels et al., 2008] or the perturbation approach to the system additionally having an in-plane magnetized reference layer [Lacoste et al., 2013] have been made to develop practical theory.

The application of an external magnetic field tilted from the perpendicular axis also breaks the symmetry and enables us to measure the oscillation experimentally. In other geometries, the experimental studies have shown that the oscillation properties such as the threshold current to excite the self-oscillation strongly depend on the field direction [Rippard et al., 2004; Tamaru et al., 2014]. On the other hand, the role of the magnetic field on the self-oscillation properties in this geometry has not been fully understood yet. For example, it is still unclear how much current is necessary to excite the out-of-plane self-oscillation in the presence of the magnetic field pointing in an arbitrary direction, while it is known that infinitesimal current can excite the self-oscillation for the highly symmetric case [Lee et al., 2005; Silva and Keller, 2010].

In this paper, we investigate theoretical conditions to excite the self-oscillation in a spin-torque oscillator with an in-plane magnetized free layer and a perpendicularly magnetized pinned layer in the presence of an external magnetic field. We solve the Landau-Lifshitz-Gilbert (LLG) equation both numerically and analytically. The main findings in this paper are as follows. First, we find that the linearized LLG equation is no longer useful to evaluate the instability threshold in the present system. The critical current density evaluated from the linearized LLG equation is two orders of magnitude larger than the instability threshold estimated from the numerical simulation. Second, we derive the theoretical formula determining the instability threshold. The main difference between the linear analysis and our result is that when a periodic precession around the stable state is assumed in the linear analysis, while we focus on the transition of the magnetization from the stable state to the out-of-plane self-oscillation state during a time shorter than the precession period. A good agreement between the numerical simulation and our formula is obtained. Third, we derive theoretical conditions to guarantee the present results, i.e., the condition that our formula of the threshold current density works better than the linear analysis to estimate the instability threshold, and the condition to stabilize the out-of-plane self-oscillation.

This paper is organized as follows. In Sec. II, we show the numerical simulation results near the instability of the initial state. We also solve the linearized LLG equation analytically. In Sec. III, we derive a theoretical formula of the threshold current, and confirm its validity by comparing the results obtained from the formula with the numerical simulation. The conclusion is summarized in Sec. IV.

## Ii Numerical simulation and linear analysis

In this section, we investigate the threshold current density which is necessary to destabilize the magnetization in the stable state by numerically solving the LLG equation. We also compare the numerical result with the analytical values of the critical current density estimated from the linearized LLG equation. Throughout this paper, the term ”threshold current” indicates the current destabilizing the stable state calculated in the numerical simulation or from the formula which is also well consistent with the numerical simulation, while the term ”critical current” is a current estimated from the linearized LLG equation.

### ii.1 System description

The system we consider is schematically shown in Fig. 1(a). The axis is perpendicular to the film plane. The unit vectors pointing in the magnetization direction of the free and pinned layers are denoted as and , respectively. The magnetization of the pinned layer points to the positive direction, . The positive current is defined as the electrons flowing from the free layer to the pinned layer. We use the macrospin approximation to the free layer. The magnetization dynamics is described by the LLG equation

 dmdt=−γm×H−γHsm×(p×m)+αm×dmdt, (1)

where and are the gyromagnetic ratio and the Gilbert damping constant, respectively. We use the approximation because the damping constant for typical ferromagnets is on the order of [Oogane et al., 2006; Tsunegi et al., 2014]. The spin-torque strength is

 Hs=ℏηj2eMd, (2)

where is the spin polarization of the electric current density , while and are the saturation magnetization and the thickness of the free layer, respectively. We neglect the asymmetry of the spin torque described by the term [Slonczewski, 1996] here, for simplicity. The critical current density in the presence of this factor, as well as its role, is briefly summarized in Appendix A. The magnetic field consists of the demagnetization field along the direction, , and the applied expressed as

 H=Happl−4πMmzez. (3)

The applied field is tilted from the axis and assumed to lie in the plane for convention, i.e.,

 Happl=HapplsinθHex+HapplcosθHez, (4)

where and are the amplitude and the tilted angle from the axis of the applied field, respectively. The magnetic field relates to the energy density via , which in the present system is

 E=−MHappl(sinθHmx+cosθHmz)+2πM2m2z. (5)

Here, we assume that , and therefore, the stable state, i.e., the minimum of Eq. (5), locates close to the axis. Note that the magnetization dynamics described by Eq. (1) can be regarded as the motion of a point particle on a unit sphere.

The values of the parameters used in this section are brought from typical experiments [Hiramatsu et al., 2016], emu/c.c., rad/(Oe s), , nm, and . The magnitude of the applied field is Oe, while the field angle is . Figure 1(b) shows the constant energy curves of Eq. (5) with these parameters. Note that the stable (minimum energy) state, the saddle point, and the unstable (local maximum) states of the energy density all exist in the plane. The stable state locates in the positive region, while the saddle point exists in the negative region. Also, the unstable states slightly shift from the axis due to the applied field. We denote the energies corresponding to the stable state, the saddle point, and the unstable states as , , and , where the subscript distinguishes the unstable states in the positive () and negative () region. For , . The constant energy curves in Fig. 1(b) are classified to the ellipses around the and axes. The energy density corresponding to the curves around the axis is in the region of , while that for the curves around the axis is in the region of .

### ii.2 Linear analysis

The conventional method to estimate the minimum current density to destabilize the stable state is linearizing the LLG equation and investigating the oscillating solution of the magnetization with a complex frequency [Sun, 2000; Grollier et al., 2003; Morise and Nakamura, 2005]. In this section, we derive the theoretical formula of the critical current density and estimate its value.

We introduce the zenith and azimuth angles as to identify the magnetization direction. In particular, the angles corresponding to the stable state are denoted as . In the present case, , and is determined by the condition ,

 Happlsin(θH−θ0)+4πMsinθ0cosθ0=0. (6)

We introduce a new coordinate where the axis is parallel to the magnetization in the stable state . A small amplitude oscillation of the magnetization around a stationary point is described by the following linearized LLG equation (the detail of the derivation is shown in Appendix A)

 1γddt(mXmY)+M(mXmY)=Hs(sinθ00), (7)

where

 M=(αHX−Hscosθ0HY−HXαHY−Hscosθ0) (8)

with and . The solution of Eq. (7) has a form of . The critical current density is defined as the current density satisfying . For a small , this condition is approximated to because . Therefore, the critical current density becomes

 jc=2αeMdℏηcosθ0[Happlcos(θH−θ0)−4πMcos2θ0+cos2θ02]. (9)

Substituting the above parameters, we find that and A/cm.

### ii.3 Numerical simulation

Figures 2(a)-(d) show the magnetization dynamics on the unit sphere and time developments of the components of , obtained by numerically solving the LLG equation, Eq. (1). The current density is (a) 7.2, (b) 7.3, (c) -7.2, and (d) -7.3 A/cm. As shown, when the current magnitude is smaller than A/cm, the magnetization finally moves to another point and stops its dynamics. On the other hand, the magnetization shows the self-oscillation for A/cm. The component of the magnetization moves to the positive (negative) direction for the negative (positive) current because the negative (positive) current prefers to be parallel (antiparallel) to the magnetization of the pinned layer, .

Three important conclusions are obtained from Fig. 2. First, the threshold current density to destabilize the initial stable state, A/cm, is two orders of magnitude smaller than the critical current density, A/cm, estimated from the linearized LLG equation. Second, both positive and negative currents can destabilize the initial state, while the sign of is fixed (positive for ). Third, the magnetization precesses around the axis above the threshold. Note that the self-oscillation occurs on the trajectory close to the constant energy curve. Although the energy landscape has the constant energy curves around the axis, as shown in Fig. 1(b), an in-plane precession around the axis does not appear. In the next section, we explain the physical meanings of such behavior.

## Iii Theoretical formula of threshold current

The results discussed in the previous section indicates that the linear analysis is no longer applicable to evaluate the instability threshold, although the linear analysis has been widely used to analyze the spin torque induced magnetization dynamics [Sun, 2000; Grollier et al., 2003; Morise and Nakamura, 2005]. In this section, we clarify the reason for the breakdown of the linear analysis, and derive a theoretical formula of the threshold current density by focusing on the energy gain of the free layer generated from the work done by spin torque.

### iii.1 LLG equation averaged over constant energy curves

Here, let us discuss the averaging technique of the LLG equation on the constant energy curves. This method has been used in several works to analyze the self-oscillation and the thermally activated magnetization switching induced by spin torque, the microwave assisted magnetization reversal, and so on [Bertotti et al., 2004, 2005; Serpico et al., 2005; Bertotti et al., 2009; Apalkov and Visscher, 2005; Hillebrands and Thiaville, 2006; Bazaliy and Arammash, 2011; Dykman, 2012; Newhall and Eijnden, 2013; Taniguchi et al., 2013; Taniguchi, 2014, 2015; Pinna et al., 2014]. As will be discussed below, the critical current density introduced above corresponds to a special limit of this averaging technique. Therefore, by reviewing the derivation of the averaged LLG equation, the reason why the linearized LLG equation does not work to estimate the instability condition accurately will be clarified.

The self-oscillation is a steady precession on a constant energy curve of excited by the magnetic field torque (). To maintain the precession, the spin torque should balance with the damping torque. Note however that the spin torque and the damping torque have different angular dependences. Therefore, strictly speaking, the spin torque may overcome the damping torque at certain points on the precession trajectory, the damping torque may however overcome the spin torque at other points. The self-oscillation is maintained when the shift from the constant energy curve due to the imbalance between the spin torque and the damping torque is sufficiently small. In that case, the magnetization can return back to the original constant energy curve during the precession. When this condition is satisfied, we obtain the following averaged LLG equation,

 ∮dtdEdt=Ws(E)+Wα(E), (10)

where the integral range is a precession period on a constant energy curve of . The work done by spin torque and the dissipation due to the damping during the precession are

 Ws=∮dtγMHs[p⋅H−(m⋅p)(m⋅H)], (11)
 Missing \left or extra \right (12)

respectively. Since the energy density averaged over the precession is conserved in the self-oscillation state, the self-oscillation is described by the equation . Therefore, the current density necessary to excite a self-oscillation on a certain constant energy curve of is

 j(E)=2αeMdℏη∮dt[H2−(m⋅H)2]∮dt[p⋅H−(m⋅p)(m⋅H)]. (13)

The explicit form of for an arbitrary is obtained, in principle, by substituting the solution of the precession trajectory on a constant energy curve, which is described by . However, the solution is hardly obtained because the equation is a nonlinear equation. Therefore, we evaluate the integrals in Eq. (13) numerically, except for special cases mentioned below (see also Appendix B). The technique to evaluate the integrals in Eq. (13) is shown, for example, in Ref. [Taniguchi, 2014]. The damping constant is assumed to be scalar in the above formulation. On the other hand, a tensor damping was proposed in Ref. [Safonov, 2002]. The presence of the tensor damping was also suggested in the spin-torque problem [Zhang and Zhang, 2009]. The effect of the tensor damping can be taken into account by replacing in Eq. (12) with the tensor damping; see Appendix C of Ref. [Taniguchi, 2014].

### iii.2 Derivation of threshold current

Note that the critical current density , Eq. (9), obtained from the linearized LLG equation relates to Eq. (13) via

 jc=limE→Eminj(E). (14)

Therefore, the fact that the critical current density is quite larger than the threshold current density found in the numerical simulation indicates the breakdown of applying averaged LLG equation.

An important assumption in the averaged LLG equation is that the magnitudes of the spin torque and the damping torque are sufficiently small. Thus, a shift of the magnetization from a constant energy curve due to the imbalance between these torques is also small. However, this assumption is not satisfied in the present case. Figure 3 shows the trajectory of the magnetization dynamics obtained from the numerical simulation, where the current density A/cm is the threshold value found in Fig. 2(d). We also show the constant energy curve including the saddle point . We should remind the readers that there are two kinds of constant energy curves, as shown in Fig. 1(b), i.e., the curves around the axis corresponding to and the curves around the axis corresponding to . The constant energy curve of separates these in-plane and out-of-plane regions. As shown in Fig. 3, while the magnetization moves from the initial state to a point close to the saddle point, the magnetization crosses the constant energy curves of , and transfers from the in-plane region to the out-of-plane region. A periodic oscillation around the stable state ( axis) is not excited. This result is the evidence that the assumption used in Eq. (13), as well as Eq. (9), is broken. Therefore, the critical current density does not work to estimate the instability of the magnetization around the stable state accurately.

The inapplicability of the linearized LLG equation also relates to the value of the damping constant . Note that both the spin torque and the damping torque move the magnetization from a constant energy curve, by either supplying or dissipating the energy from the free layer. Therefore, the averaging technique of the LLG equation, as well as the linearization of the LLG equation, works well for low damping case. The fact that the linearized LLG equation could not be applied in the above numerical simulation indicates that the value of the damping constant in the present system is high and that the precession around the stable state is not stabilized. The range of the damping constant where the linearized LLG equation will be applicable is discussed in Sec. III.3 below.

Figure 3 suggests that the magnetization can climb up the energy barrier by absorbing energy due to the work done by the spin torque during a time shorter than a precession period around the stable state. Therefore, the threshold current density can be defined as a current density satisfying the following equation,

where corresponds to the initial stable state. Strictly speaking, the exact solution of the LLG equation is necessary to evaluate the threshold current density from Eq. (15). However, the LLG equation is a nonlinear equation, and it is difficult to obtain the exact solution. Instead, we approximate Eq. (15) as

where are the points on the constant energy curve of and are located in the plane; see Fig. 3. We note that Eq. (15) is well approximated by Eq. (16) when locate close to , which means that . Note that the left hand side of Eq. (16) can be evaluated in a similar manner to calculating Eq. (13) because the integral range is on the constant energy curve. However, the integral range is limited to in Eq. (16), while the range is over a periodic precession in Eq. (13). The values of the integrals for these different regions are, in general, different. Since the value of the integral in Eq. (16) is determined by the energy landscape, and the time-dependent solution of Eq. (1) is unnecessary, the integral in Eq. (16) is more easily evaluated than that in Eq. (15) [Bertotti et al., 2009].

The current density satisfying Eq. (16) is given by

Equation (17) is the theoretical formula of the threshold current density and is the main result in this paper. This equation provides the estimation of the threshold current density with high accuracy. For example, the values of with the parameters used in Fig. 2 are A/cm and A/cm, which show good agreement with the numerical results in Fig. 2. These values are estimated for . Below, we show that the agreement between Eq. (17) and the numerical simulation is obtained also for different values of ; see Fig. 5. Note that because the magnetic field pointing in the positive direction breaks the symmetry between the magnetization dynamics moving to the positive and negative directions, although the difference is small. We emphasize that Eq. (17) consists of two parts. One is proportional to because this term arises from the energy dissipation due to the damping. The other is, on the other hand, independent of but proportional to the energy barrier .

Equation (17) can be simplified into a different form for (see also Appendix C). In this case, , , and with and . Then, we find that

 ∫mdmd±dt[p⋅H−(m⋅p)(m⋅H)]=∓πγ(1−h)2, (18)
 ∫mdmd±dt[H2−(m⋅H)2]=16πM3γ√h(1−h)(3−5h+2h2). (19)

Therefore, Eq. (17) becomes

 jth±(θH=90∘)=∓2eM2dℏη×[16α3√h1−h(3−2h)+8h(1−h)2]. (20)

Note that in the limit of , indicating that infinitesimal current can destabilize the stable state in the absence of the applied field.

We note that both the positive and negative currents can destabilize the stable state in our picture, contrary to having a fixed sign (positive for ). The physical meaning of this difference is as follows. Since the damping torque always dissipates energy from the free layer, positive energy should be supplied from the work done by spin torque to destabilize the stable state. In the derivation of , a steady precession around the stable state is assumed. On the precession trajectory, the spin torque has a component antiparallel to the damping torque when and has a component parallel to the damping torque when , for a positive current. The spin torque supplies energy to the free layer in the former case, but dissipates energy from the free layer in the latter case. Note that the trajectory slightly shifts to the positive direction due to the magnetic field having the positive component, i.e., the trajectory is not symmetric with respect to the plane. Then, the work done by spin torque during the precession becomes finite and positive. The spin torque overcomes the damping torque when the current density becomes larger than . When the current direction is reversed, the work done by spin torque becomes negative, and thus, the spin torque cannot overcome the damping. As a result, the sign of is positive. However, as emphasized above, a periodic precession around the easy axis assumed in the derivation of is not excited in the present case. Instead, we focused on the magnetization dynamics from to . The work done by spin torque during becomes positive when the current has the positive sign. Similarly, the work during is positive when the current sign is negative. Therefore, both the positive and negative currents can destabilize the stable state by compensating the damping torque. Note also that the magnetization crosses the constant energy curve of during a time shorter than a precession period around the axis. Therefore, an in-plane self-oscillation on a constant energy curve of around the axis cannot be excited in the present case.

### iii.3 Applicability of the present theory

There are two characteristic current scales, and , related to the magnetization dynamics, as discussed above. These two currents are defined from different mechanisms of the instability of the stable state. The instability condition of a precession around the stable state gives . On the other hand, was derived by the condition that the energy gain by the spin torque during a time shorter than the precession period becomes larger than the energy barrier between the stable state and the saddle point. The initial state is destabilized when the current magnitude becomes larger than . For the present parameters, is smaller than , and therefore, determines the instability threshold. The condition that works well to estimate the instability of the stable state can be expressed as

 jth±jc<1. (21)

This is another important equation in this paper, guaranteeing the validity of our approach. Whether Eq. (21) is satisfied or not depends on the material parameters, as well as the applied field magnitude and angle. If Eq. (21) is unsatisfied, determines the instability threshold, the magnetization moves to the out-of-plane region after the magnetization precesses around the in-plane axis.

Note that the first term on the right hand side of Eq. (17) is proportional to the damping constant , while the second term is independent of . On the other hand, Eq. (9) is proportional to . Therefore, Eq. (21) is not satisfied when becomes sufficiently small. When Eq. (21) is unsatisfied, determines the instability of the stable state. Then, we can discuss the minimum value of guaranteeing the applicability of Eq. (21) [Taniguchi et al., 2013]. The value which falls off from the condition in Eq. (21) for the parameters used in Fig. 2 is . This value of is at least one to two orders of magnitude smaller than the experimentally reported values for conventional ferromagnets used in spin-torque oscillator, such as CoFeB [Oogane et al., 2006; Tsunegi et al., 2014]. Therefore, we consider that determines the instability of the stable state for typical experiments.

### iii.4 In the presence of the angular dependence of the spin torque

When the applied field points to the in-plane direction, , the stable state corresponds to , and the critical current density in Eq. (9) diverges. This is because the work done by spin torque during the precession around the stable state becomes zero. Therefore, Eq. (21) is always satisfied for .

The divergence of appears at a different field angle when the angular dependence of the spin torque is taken into account, although this term is neglected in the above calculation, for simplicity. In this case, Eq. (2) is replaced by

 Hs=ℏηj2e(1+λm⋅p)Md. (22)

Then, the critical current density becomes

 jc=2αeMdℏηP(θ0)[Happlcos(θH−θ0)−4πMcos2θ0+cos2θ02], (23)

where is given by

 P(θ0)=cosθ01+λcosθ0+λsin2θ02(1+λcosθ0)2, (24)

see Appendix A. The divergence of the critical current density, Eq. (23), occurs at the angle satisfying . In particular, when , Eqs. (23) becomes

 jc(θH=90∘)=4αeMdℏηλ(Happl+2πM). (25)

On the other hand, Eq. (17) is generalized for finite as

Equation (26) for is

 jth±(θH=90∘)=∓2eMdℏη4πMND±, (27)

where and are

 N=4λ2[2α(3−2h)(1−h)√h(1−h)+3h]×√1−4λ2h(1−h) (28)
 D±=3{√1−4λ2h(1−h)[π∓4λ√h(1−h)]−2[1−2λ2(1−h)]cos−1[±2λ√h(1−h)]}, (29)

see Appendix B. In the presence of a finite , even for . Eq. (27) reproduces Eq. (20) in the limit of . The currents, and , in Eq. (21) should be replaced by Eqs. (23) and (26) in the presence of the angular dependence of the spin torque.

### iii.5 Validity of Eq. (17) and condition to excite out-of-plane self-oscillation

In this section, we confirm the validity of Eq. (17) for a wide range of by comparing with the numerical simulation of the LLG equation.

Before the comparison, we briefly discuss the definition of the threshold current density estimated from the numerical simulation. We emphasize that just determines the instability of the stable state, and does not guarantee the existence of the out-of-plane self-oscillation. The out-of-plane self-oscillation is excited when a condition,

 j(E)jth±>1, (30)

is satisfied [Taniguchi, 2015], where the range of is . Note that the reason why the out-of-plane self-oscillations appear in Figs. 2(b) and 2(d) is that there exists a certain satisfying Eq. (30). On the other hand, when Eq. (30) is not satisfied for any value of , the magnetization moves to the point close to for a positive (negative) current above because the spin-torque magnitude becomes sufficiently strong, and the magnetization eventually becomes parallel or antiparallel to the magnetization of the pinned layer, . Figure 4 shows an example of such dynamics, where and the current density is close to the threshold value, A/cm, for this . As shown, the magnetization finally becomes almost parallel to the axis. Such magnetization dynamics was observed experimentally [Houssameddine et al., 2007]. The threshold current density evaluated from the numerical simulation should be defined as the current density above which the magnetization shows a stable out-of-plane self-oscillation or the magnetization moves to the points close to . The detail of the method numerically defining the threshold current density is summarized in Appendix D.

We study the validity of Eq. (17) by comparing with the threshold current estimated by numerically solving Eq. (1) for several values of and . The threshold current density estimated from the numerical simulation of the LLG equation is shown by dots in Fig. 5(a), where the magnetic field angle varies in the range of while the magnitude is fixed to 650 Oe. We also shows the value of evaluated from Eq. (17) by solid lines. We find a good agreement between the numerical and theoretical results, supporting the validity of Eq. (17). The comparison between the numerically evaluated instability threshold and the analytical formula, Eq. (20), for several values of the field magnitude is shown in Fig. 5(b), where the field angle is fixed to . The theoretical formula agrees with the numerical result when , while the numerical result becomes different with the theoretical formula for relatively large magnetic field. This is because the derivation of the theoretical formula, Eq. (17), assumes that , as mentioned below Eq. (16). The current magnitude above which the difference between the theoretical and numerical results appears is on the order of A/cm, which is one to two orders of magnitude larger than the current magnitude used in typical experiments [Houssameddine et al., 2007; Suto et al., 2012; Bosu et al., 2016; Hiramatsu et al., 2016]. Therefore, we consider that the present formula works well to analyze experiments for wide range of the applied field angle and magnitude.

We notice that is an increasing function of for the out-of-plane self-oscillation when , as in the zero field case [Lee et al., 2005; Silva and Keller, 2010]. Then, there is a certain satisfying Eq. (30) if

 ju±jth±>1, (31)

is satisfied, where are Eq. (13) at the unstable states, ,

 ju±≡limE→Emax±j(E). (32)

The dependences of and on the field angle and the magnitude are also shown in Figs. 5(a) and 5(b), respectively. It is shown that is almost independent of and , while increases with increasing these parameters. For example, we find that for . This result indicates that the out-of-plane self-oscillation can be excited for for the present parameters. This finding is consistent with the numerical results shown in Figs. 2 and 4, supporting the validity of our argument. We notice that the linearized LLG equation is useful to estimate by replacing with the zenith angle corresponding to the maximum point; see Eq. (48). In particular, when , the unstable states locate at and . Then, we find that (see also Appendix B)

 ju±(θH=90∘)=∓2αeMdℏη√1−h24πM(1−h22). (33)

This equation indicates that for , i.e., is almost independent of , which is consistent with the result shown in Fig. 5(b).

## Iv Conclusion

In conclusion, we studied the theoretical conditions to excite the self-oscillation in a spin-torque oscillator consisting of an in-plane magnetized free layer and a perpendicularly magnetized pinned layer in the presence of an external magnetic field pointing in an arbitrary direction. The numerical simulation in Fig. 2 showed that the initial stable state is destabilized by current density much smaller than the critical current density estimated from the linearized LLG equation, Eq. (9). The fact implies that the linearized LLG equation is no longer applicable to evaluate the instability threshold in the present system. Then, we derived the theoretical formula of the threshold current density, Eq. (17), by focusing on the transition of the magnetization from the stable state to the out-of-plane precession during a time shorter than a precession period around the stable state. The derived formula consists of two parts, where one is proportional to the damping constant , while the other is independent of but proportional to the energy barrier for the transition. A good agreement between the numerical simulation and our formula, Eq. (17), is obtained in Fig. 5, indicating the validity of the formula. The condition that our formula of the threshold current density works better than the linear analysis to instability threshold is Eq. (21). We also derived the theoretical condition, Eq. (30), to stabilize the out-of-plane self-oscillation.

## Acknowledgement

The authors express gratitude to Shinji Yuasa, Kay Yakushiji, Akio Fukushima, Shingo Tamaru, Sumito Tsunegi, Ryo Hiramatsu, Yoichi Shiota, Takehiko Yorozu, Hirofumi Suto, and Kiwamu Kudo for valuable discussion they had with us. This work is supported by Japan Society and Technology Agency (JST) strategic innovation promotion program ”Development of new technologies for 3D magnetic recording architecture”.

## Appendix A Derivation of linearized LLG equation

In this Appendix, we show the detail of the derivation of Eq. (7). For generality, we consider a ferromagnet having uniaxial anisotropies along the , , and axes with an external magnetic field applied in an arbitrary direction. The magnetic field is given by

 H=⎛⎜ ⎜⎝HapplsinθHcosφH−4πM~NxmxHapplsinθHsinφH−4πM~NymyHapplcosθH−4πM~Nzmz⎞⎟ ⎟⎠. (34)

The generalized demagnetization coefficient () is defined as , where is the shape anisotropy (demagnetization) field with , while is the crystalline or interface anisotropy field. The energy density is

 EM=−Happl[sinθHsinθcos(φH−φ)+cosθHcosθ]+2πM~Nxsin2θcos2φ+2πM~Nysin2θsin2φ+2πM~Nzcos2θ. (35)

The system in the main text corresponds to the case of , , and .

Since we are interested in a small oscillation of the magnetization around the stable state, the zenith and azimuth angles corresponding to the stable state should be identified. The stable state is determined by the conditions that , which are explicitly given by

 Happl[sinθHcosθcos(φH−φ)−cosθHsinθ]−4πM~Nxsinθcosθcos2φ−4πM~Nysinθcosθsin2φ+4πM~Nzsinθcosθ=0, (36)
 HapplsinθHsinθsin(φH−φ)+4πM~Nxsin2θsinφcosφ−4πM~Nysin2θsinφcosφ=0. (37)

Let us denote the zenith and azimuth angles satisfying Eqs. (36) and (37) as . As mentioned in the main text, we introduce the coordinate where the axis is parallel to the stable state . The rotation to the coordinate to the coordinate is described by the rotation matrix

 R=⎛⎜⎝cosθ00−sinθ0010sinθ00cosθ0⎞⎟⎠⎛⎜⎝cosφ0sinφ00−sinφ0cosφ00001⎞⎟⎠. (38)

The relations between the components of in the and coordinates are , , . Also, the magnetic field in the coordinate is

 H=⎛⎜⎝HXXmX+HXYmYHYXmX+HYYmYHZXmX+HZYmY+HZZ⎞⎟⎠, (39)

where

 HXX=−4πM~Nxcos2θ0cos2φ0−4πM~Nycos2θ0sin2φ0−4πM~Nzsin2θ0, (40)
 HXY=HYX=−4πM(~Ny−~Nx)cosθ0sinφ0cosφ0, (41)
 HYY=−4πM~Nxsin2φ0−4πM~Nycos2φ0, (42)
 HZX=−4πM~Nxsinθ0cosθ0cos2φ0−4πM~Nysinθ0cosθ0sin2φ0+4πM~Nzsinθ0cosθ0, (43)
 HZY=−4πM(~Ny−~Nx)sinθ0sinφ0cosφ0, (44)
 HZZ=Happl[sinθHsinθ0cos(φH−φ0)+cosθHcosθ0]−4πM~Nxsin2θ0cos2φ0−4πM~Nysin2θ0sin2φ0−4πM~Nzcos2θ0. (45)

Similarly, the magnetization of the pinned layer in the coordinate transforms in the coordinate to

 p≡⎛⎜⎝pXpYpZ⎞⎟⎠=⎛⎜⎝sinθpcosθ0cos(φp−φ0)−cosθpsinθ0sinθpsin(φp−φ0)sinθpsinθ0cos(φp−φ0)+cosθpcosθ0⎞⎟⎠. (46)

Now we consider a small oscillation of the magnetization around the stable state. Using the approximations and , the LLG equation is linearized as

 Unknown environment '% (47)

where and . The terms proportional to are neglected because these terms are on the order of . The condition that the trace of the coefficient matrix is zero gives

 jc=2αeMdℏηpZ(HX+HY2). (48)

Substituting , , , and , Eq. (48) reproduces Eq. (9). On the other hand, in the case of the in-plane magnetized system considered in Ref. [Grollier et al., 2003], i.e., , , , , , , and , we find that , , and , where is the in-plane anisotropy. Then, the critical current density becomes , which is consistent with the result in Ref. [Grollier et al., 2003].

The angular dependence of the spin torque, characterized by the factor , can be taken into account as follows. As mentioned in the main text, in this case is given by Eq. (22). In this case, Eq. (2) is replaced by Eq. (22). The factor is linearized as

 11+λm⋅p=11+λmZpZ11+λ(mXpX+mYpY)1+λmZpZ≃11+λpZ[1−λ(mXpX+mYpY)1+λpZ]. (49)

We introduce the following notations,

 H(0)s=ℏηj2e(1+λpZ)Md, (50)
 Λ=λ1+λpZ. (51)

Then, Eq. (47) becomes

 1γddt(mXmY)+M(mXmY)=−H(0)s(pXpY),