Bifurcation to complex dynamics in largely modulated voltage-controlled parametric oscillator

An experimental demonstration of a parametric oscillation of a magnetization in a ferromagnet was performed recently by applying a microwave voltage, indicating the potential to be applied in a switching method in non-volatile memories. In the previous works, the modulation of a perpendicular magnetic anisotropy field produced by the microwave voltage was small compared with an external magnetic field pointing in an in-plane direction. A recent trend is, however, opposite, where an efficiency of the voltage controlled magnetic anisotropy (VCMA) effect is increased significantly by material research and thus, the modulated magnetic anisotropy field can be larger than the external magnetic field. Here, we solved the Landau–Lifshitz–Gilbert equation numerically and investigated the magnetization dynamics driven under a wide range of the microwave VCMA effect. We evaluated bifurcation diagrams, which summarize local maxima of the magnetization dynamics. For low modulation amplitudes, the local maximum is a single point because the dynamics is the periodic parametric oscillation. The bifurcation diagrams show distributions of the local maxima when the microwave magnetic anisotropy field becomes larger than the external magnetic field. The appearance of this broadened distribution indicates complex dynamics such as chaotic and transient-chaotic behaviors, which were confirmed from an analysis of temporal dynamics.

Voltage controlled magnetic anisotropy (VCMA) effect 1 modulates a perpendicular magnetic anisotropy at a ferromagnetic metal/nonmagnetic insulating layer interface by modulating of electron states near the interface [2][3][4] and/or inducing magnetic moments 5 .It enables us to manipulate the direction of the magnetization in a ferromagnet electrically without the Joule heating, and thus, is expected to be a new writing method in magnetoresistive random access memory (MRAM), whose writing scheme currently relies on spin-transfer torque 6,7 .The material researches for highly efficient VCMA effect has been reported, where the perpendicular magnetic anisotropy is largely modulated by a small voltage [8][9][10][11][12][13][14][15][16][17] .The efficiency recently achieved reaches to about 300 fJ/(Vm) 18 , which corresponds to a magnetic anisotropy field on the order of kilo-Oersted for typical VCMA-based MRAM.At the same time, analyses on the magnetization dynamics driven by the VCMA effect have been investigated both experimentally and numerically [19][20][21][22][23][24][25][26][27][28] .It has been revealed that the switching is unstable when the pulse width of the voltage is short because the dynamics becomes very sensitive to the pulse shape in such condition 24 .For a long-pulse regime, however, the switching also becomes unstable due to noise 23 .
To overcome the issue, a parametric oscillation of the magnetization by applying microwave voltage was proposed in Ref. 29 .For the magnetization switching in the VCMA-based MRAM using a perpendicularly magnetized free layer, an external magnetic field H appl pointing in an in-plane direction is applied and induces the magnetization precession around H appl , which eventually relaxes to the direction of H appl .When the microwave voltage witn an oscillation frequency f = 2 f L , where f L = γH appl /(2π) (γ is the gyromagnetic ratio) is the Larmor frequency, is applied, however, a sustainable oscillation of the magnetization is excited.This oscillation is classified into the parametric oscillation, and for simplicity, we call the oscillator as voltage-controlled parametric oscillator.Since this oscillation is stable, it can be used to manipulate the magnetization direction even in a long-pulse regime, which results in a reliable magnetization switching.Note that the previous work 29 focuses on a parameter region of H Ka /H appl ≪ 1, where H Ka is the amplitude of the modulated magnetic anisotropy field generated by the micowave VCMA effect.The recent progress of the VCMA efficiency, however, makes the opposite limit, H Ka /H appl ≫ 1, available because the value of H Ka is growing rapidly, as mentioned above.The dynamical behavior of the magnetization in this limit has not been investigated yet.
In this work, we study the magnetization dynamics driven by the microwave VCMA effect by solving the Landau-Lifshitz-Gilbert (LLG) equation numerically.We change the value of H Ka in a wide range including the limit of H Ka /H appl ≫ 1.To clarify the change of the dynamical behavior, bifurcation diagrams are evaluated, which summarize local maxima of the oscillating magnetization.When the magnitude of H Ka is small, the parametric oscillation is induced, and the bifurcation diagram becomes single points because the dynamics is periodic.When H Ka becomes larger than H appl , however, the bifurcation

System description
In Fig. 1(a), we show a schematic illustration of a ferromagnetic/nonmagnetic/ferromagnetic trilayer.The top and bottom ferromagnets correspond to free and reference layer, respectively.We apply a macrospin assumption in the free layer, whose validity in dynamical state driven by the VCMA effect has been confirmed experimentally 29 .Thus, the unit vector m pointing in the magnetization direction in the free layer conserves its magnitude, i.e., |m| = 1 and d|m|/dt = 0.An external magnetic field H appl is applied to an in-plane direction.We assume that the shape of the free layer is a cylinder, and therefore, the free layer does not have an in-plane magnetic anisotropy.For convenience, we use a Cartesian coordinate, where the x axis is parallel to H appl while the z axis is normal to the film plane.The magnetic field in the free layer H , in the absence of an applied voltage, is given by where e k (k = x, y, z) is the unit vector in the k direction.The perpendicular magnetic anisotropy field includes the contribution from the interfacial magnetic anisotropy field H K 30-32 and the shape magnetic anisotropy field 4πM(N z − N x ) with the demagnetization coefficients N k (N x = N y due to the cylindrical symmetry).The net perpendicular magnetic anisotropy field, H K − 4πM(N z − N x ), determines the retention time of MRAM.In the conventional scheme of the writing in the VCMA-based MRAM, the direct voltage modulates this net perpendicular magnetic anisotropy close to zero to excite the magnetization precession around the external magnetic field 33 .If the perpendicular component largely remains finite, the magnetization just moves its direction to the direction of the external field with the angle sin −1 (H appl /H ′ K ) and the precession cannot be excited, where H ′ K is the reduced perpendicular magnetic anisotropy field by the VCMA effect.In fact, the numerical simulation in Ref. 33 assumes that the net perpendicular magnetic anisotropy field is completely canceled by the VCMA effect.Note that this assumption is important not only for making the simulation simple but also for experiments.If the direct (or intrinsic) component of the perpendicular magnetic anisotropy field remains finite during the precession, an instantaneous frequency becomes nonuniform.Such a nonuniform frequency will increase the switching error because the pulse width of the voltage for writing the bit is determined as a half of the Larmor precession period for VCMA-based MRAM.Therefore, it is preferable to make the direct component of the perpendicular magnetic anisotropy field zero during the switching for the conventional switching scheme.In the parametric oscillation state, both direct and microwave voltages are applied to the trilayer, and the direct voltage modulates the total perpendicular magnetic anisotropy so that it becomes close to zero 29 , while the microwave voltage provides an oscillating magnetic anisotropy field.Accordingly, the magnetic field used in the following calculation becomes, where H Ka and f are the magnitude and frequency of the magnetic anisotropy field due to the microwave voltage.The magnetization dynamics driven by this magnetic field is described by the LLG equation, where α is the damping constant.Throughout this paper, we use the values of γ = 1.764 × 10 7 rad/(Oe s) and α = 0.005.The preparation of the initial state of m by investigating thermal equilibrium is explained in Methods 34,35 .Recall that the LLG equation conserves the norm of m as |m| = 1.Therefore, although m is a three-dimensional vector, its dynamical degree of freedom is two due to this constraint.In fact, if we use a spherical coordinate, for example, the dynamics of m is described by two variables (zenith and azimuth angles).It should be noted that dynamical systems described by differential equations cannot show chaotic behavior when the dynamical degree of freedom is less than or equal to two, according to the Poincaré-Bendixson theorem 36 .The presence of the microwave voltage, however, makes the present system non-autonomous and provides a possibility to excite chaotic behavior, as shown below.

Results
Here, we study the change of the magnetization dynamics for various magnitude of H Ka .

Parametric oscillation
Let us first start by confirming the parametric oscillation studied previously 29 .It was shown in Ref. 29 that a sustainable oscillation of the magnetization is excited when the frequency of the microwave voltage, f , is twice the Larmor frequency, f L = γH appl /(2π).Therefore, in the following, we fix the value of f to be f = 2 f L . Figure 1(b) shows the dynamical trajectory of the magnetization in a steady state, where H appl = 300 Oe while H Ka = 100 Oe, i.e., H Ka /H appl ≪ 1, as in the case of the previous work 29 .Since |m| = 1 is satisfied in the LLG equation, it is useful to draw the dynamical trajectory on a unit sphere, as shown in this figure.Time evolution of m z in a steady state is also shown in Fig. 1(c).These results indicate an appearance of the sustainable oscillation of the magnetization mentioned above.In Fig. 1(d), the Fourier spectrum of m z is shown, where the inset shows it around the main peak in a linear scale.Its main peak appears at 0.84 GHz, which is the same with f L with H appl = 300 Oe.These results are consistent with the previous works 29 .Since the magnetization switches its direction between m z ≃ +1 and m z ≃ −1 periodically with the period 1/(2 f L ), this parametric oscillation can be used as a switching scheme in VCMA-based MRAM 29 .Recall that this oscillation is sustained by the microwave modulation of the magnetic anisotropy; if this time-dependent modulation is absent, the magnetization monotonically relaxes to the direction of the external magnetic field.Spin-wave propagation through a parametric excitation is another example of the magnetization dynamics caused by microwave voltage, which has been studied previously [37][38][39][40] .

Appearance of complex dynamics and bifurcation diagram
When the value of H Ka further increases, the magnetization dynamics becomes complex.In Fig. 2(a), we show the dynamical trajectory of the magnetization for H Ka = 620 Oe, while H appl = 300 Oe is the same with that used in Fig. 1 indicate that the application of the microwave voltage is no longer applicable to the switching method for the VCMA-based MRAM when the modulation of the magnetic anisotropy field becomes larger than the external magnetic field due to the breakdown of the simple parametric oscillation.
A way to distinguish these dynamics, such as the difference between Figs. 1(b) and 2(a), qualitatively is to draw a bifurcation diagram 36 .In Fig. 3(a), we summarize local maxima of m z for various H Ka , where H appl = 300 Oe.Recall that the parametric  oscillation is excited when H Ka is small.In this case, the dynamics is periodic and m z is similar to a simple trigonometric function, as shown in Fig. 1(c).The local maxima of m z for this case, thus, saturate to a single point.When the dynamics become complex, the bifurcation diagram shows broadened structure.For example, there are three points at H Ka = 620 Oe in Fig. 3(a), the validity of which is confirmed from Fig. 2(b).When the value of H Ka further increases, the bifurcation diagram shows largely broadened structure and finally shows simplified structure again.As discussed below, these correspond to chaotic and transient-chaotic behaviors [41][42][43] .Note that the appearance of complex but still periodic structure might weakly depend on the initial state (see Methods) while the region of the broadened structure might depend on the measurement time, as will be mentioned below.It is difficult to analytically estimate the value of H Ka at which the bifurcation from a simple parametric oscillation to a complex oscillation occurs; see Methods.However, the numerical simulations for various parameters imply that the complex oscillation appears when H Ka becomes larger than H appl , as discussed below.We also evaluated Lyapunov exponent 36 by Shimada-Nagashima method 44 , as shown in Fig. 3(b).The Lyapunov exponent is an inverse of a time scale characterizing an expansion of a difference between two solutions of the LLG equation with an infinitesimally different initial conditions; see Methods explaining the evaluation method of the Lyapunov exponent.A negative Lyapunov exponent means that the magnetization saturates to a fixed point.The Lyapunov exponent is zero when the magnetization dynamics are periodic, while it becomes positive when the dynamics are chaotic.We notice that the Lyapunov exponent in the present system is zero or positive, depending on the value of H Ka .Thus, the Lyapunov exponent can be used as an indicator to distinguish between periodic oscillation and chaotic dynamics.For example, we can conclude that the dynamics in Fig. 2(a) is periodic not only from the temporal dynamics shown in Fig. 2(b) but also from the fact that the Lyapunov exponent for H Ka = 620 Oe is zero.

Chaotic dynamics
In Fig. 4(a), we show the dynamical trajectory of the magnetization for H Ka = 1200 Oe.The dynamical trajectory covers almost the whole region of the unit sphere.The time evolution of m z becomes non-periodic, as shown in Fig. 4(b), and the Fourier spectrum shows a broad structure having several peaks.These results imply chaotic dynamics of the magnetization 41,42 .The appearance of chaos for this parameter can also be concluded from the fact that the Lyapunov exponent shown in Fig. 3(b) is positive.

4/14
As mentioned above, chaotic dynamics in the present system are excited because of the presence of the microwave voltage.Similar examples in spintronics devices have been found in spin-torque oscillators with time-dependent inputs [45][46][47] .The differences of the phenomena observed between the voltage-controlled parametric oscillator studied here and spin-torque oscillator are as follows.In the spin-torque oscillators, a sustainable oscillation of the magnetization is driven by direct currents, and an injection of time-dependent inputs is not a necessary condition for the oscillation.Spin-torque oscillator cannot show chaotic behavior when only the direct current is injected due to the constraint by the Poincaré-Bendixson theorem.A way to excite chaos in spin-torque oscillator is to inject time-dependent inputs such as alternating current and/or magnetic field [45][46][47] .When the magnitudes of these time-dependent inputs are relatively small, synchronization may occur.When their magnitudes are further increased, chaos might be induced.On the other hand, the sustainable oscillation of the magnetization in the present voltage-controlled parametric oscillator is driven by microwave voltage.In other words, this time-dependent input is a necessary condition for the oscillation.Chaotic dynamics appear when the magnitude of the microwave voltage becomes relatively large, as shown in Fig. 4.
Let us briefly mention applicability of the chaotic dynamics for practical devices.Chaos in spintronics devices has been studied both experimentally and theoretically using various methods [48][49][50][51][52][53][54][55][56][57] .An excitation of chaos in spintronics devices may evoke interest from a viewpoint of brain-inspired computing 58,59 .For example, it was found that a computational capability of physical reservoir computing is enhanced when spin-torque oscillators are near the edge of chaos 58 , where chaos was excited by adding another ferromagnet to the oscillator.An excitation of chaos in spin-torque oscillator by time-dependent inputs, however, seems to require large power consumption.6][47] is on the order of 10 7 − 10 8 A/cm 2 .The value is larger than the switching current density of the state of the art spin-transfer torque driven MRAM 60 , and thus, is hardly desirable.The large current also causes large power consumption, which is also unsuitable for practical applications.The voltage-controlled parametric oscillator, on the other hand, ideally reduces the power consumption significantly.In fact, this point has been a motivation for developing VCMA-based MRAM.However, these VCMA-based devices often require an external magnetic field for both the switching and parametric oscillation, which is not preferable in practical applications.A way to solve the issue might be to use an effective field, instead of applying external magnetic field, such as an interlayer exchange coupling field, as investigated in the study of spin-orbit torque driven MRAM 61 .Physical reservoir computing by the VCMA-based MRAM without an external magnetic field was investigated recently 62 , however, switching nor parametric oscillation was used there.A development of the voltage-controlled parametric oscillator without requiring an external magnetic field will be of interest as a future work for applying it to practical applications.

Transient chaos
In Fig. 5(a), we show the dynamical trajectory of the magnetization for H Ka = 1600 Oe.The time evolution of m z and its Fourier transformation are also shown in Figs.5(b) and 5(c), respectively.We note that the dynamics shown here corresponds to m z (t) in a long-time limit, i.e., a steady state.The results look similar to those shown in Figs. 1 and 2. Also, the local maxima of m z in the bifurcation diagram in Fig. 3(a) concentrates on a single point.However, there is a critical difference between the magnetization dynamics shown in Fig. 5 and those shown in Figs. 1 and 2.
To clarify their differences, it is necessary to focus on a process to reach the steady state.In Figs.6(a) and 6(b), we show the time evolution of m z (t) from t = 0 to the steady states for H Ka = 620 Oe and 1600 Oe, respectively.When H Ka is relatively   43 , where dynamical systems initially show chaotic behavior but it suddenly disappears.Time necessary to move from chaos to a steady state is sensitive to various parameters and initial conditions of systems 43 .For example, in Figs.7(a), 7(b), and 7(c), we show time evolution of m z for H Ka = 1900 Oe, 1980 Oe, and 2000 Oe, respectively, where the initial conditions are common.The results indicate that the time necessary to realize the steady state depends on the parameter H Ka and it does not show, for example, a monotonic change with respect to the change of the parameter.The transient chaos in spintronics devices was predicted in a spin-torque oscillator with a delayed-feedback circuit 54 .The phenomenon has not been verified experimentally yet, although chaos was confirmed recently 57 .
We note that the classification of chaos and transient chaos depends on a measurement time 43 .For example, in the present work, we solve the LLG equation from t = 0 to t = 5 µs and classify the magnetization dynamics.If we change this maximum time (5 µs) to, for example 2 µs, the dynamics for H Ka = 2000 Oe, shown in Fig. 7(c), will be classified as chaos.Another example can be seen in the bifurcation diagram and the Lyapunov exponent shown in Fig. 3, where a broadened structure in the bifurcation diagram and a positive Lyapunov for H Ka = 1920 Oe, indicating chaos for this parameter.If we measure the dynamics for this parameter longer, however, the dynamics might change to a steady state; in such a case, the dynamics will be classified as transient chaos.As a general knowledge on transient chaos 43 , we should remember that the classification of chaos and transient chaos has such an arbitrary property.

Bifurcation diagrams for various applied fields
As mentioned above, a bifurcation diagram is useful to classify the magnetization dynamics, although, for example, the difference between a simple steady oscillation and a transient chaos should be verified from temporal dynamics, as mentioned above.We also notice that the whole structure of the bifurcation diagram does not change qualitatively even when the initial condition is slightly changed.Recall that the present system includes only three parameters, H Ka , H appl , and α.Therefore, in Figs.8(a)-8(c), we show the bifurcation diagrams for (a) H appl = 200 Oe, (b) H appl = 300 Oe, and (c) H appl = 400 Oe.These figures indicate that the local maxima of m z for relatively small H Ka concentrate on single points, which indicate the excitation of the parametric oscillation.As H Ka increases, broadened structures appear, i.e., the local maxima of m z have various values, implying the appearance of complex oscillation and chaos.The figures also imply that this boundary between the parametric oscillation and the complex dynamics, i.e., the boundary between the single and multiple points in the bifurcation diagram, locates near H Ka /H appl ≃ 1.5 − 2.0.We examine similar calculations for different α but the boundary seems to be unchanged; see Figs. 8(d)-8(f), which are the bifurcation diagrams for α = 0.020.Because of high nonlinearity in the LLG equation, it is difficult to analytically verify these indications; however, it might be useful to design, for example, the modulation voltage for VCMA-based MRAM utilizing the parametric oscillation.We keep this question as a future work.
We note that the value of H Ka /H appl ≃ 1.5 − 2.0 is available by current technology.The typical value of H appl is on the order of 100 Oe; for example, it is 250 Oe in Ref. 23 and 720 Oe in Ref. 29 .Although the value of H appl can be further large experimentally, a large H appl might be unsuitable for practical applications; see Methods for analytical treatment of the LLG equation.On the other hand, the value of H Ka relates to the VCMA efficiency η, the saturation magnetization M, the applied voltage V appl , and the thicknesses of the free and insulating layers, d F and d I , via H Ka = (2ηV appl )/(Md F d I ).Substituting their typical values [η is about 300 fJ/(Vm) 18 , V appl is about 0.5 V at maximum, M is about 1000 emu/cm 3 , d F is about 1 nm, and d I is about 2.5 nm), H Ka can be on the order of kilo Oersted at maximum, as written in the introduction.Therefore, we believe that the value of H Ka /H appl ≃ 1.5 − 2.0 is experimentally achievable.

Conclusion
In summary, the magnetization dynamics in the voltage-controlled parametric oscillator for a large microwave voltage limit were investigated by numerical simulation of the LLG equation.As the modulation increases, the dynamical trajectory changes from a simple parametric oscillation to complex oscillations, which are still periodic but have several local maxima in the amplitude.Such dynamics will be of little preference in a switching scheme for VCMA-based MRAM.The evaluation of the bifurcation diagrams for various values of the external magnetic field and the damping constant indicated that the simple parametric oscillation is broken when the amplitude of the modulated magnetic anisotropy field becomes larger than the external magnetic field, and this bifurcation point seems to be hardly sensitive to the value of the damping constant.A further enhancement of the microwave modulation leads to chaotic and transient-chaotic behaviors, which might make the voltage-controlled parametric oscillator applicable to other usage in electric devices.These dynamics were classified from the bifurcation diagrams, Lyapunov exponent, and analyses on temporal dynamics.

Preparation of initial state
We prepare natural initial states by solving the LLG equation with thermal fluctuation 34 .For this purpose, we use Eq.(1) as the magnetic field.We add a torque, −γm × h, due to thermal fluctuation to the right-hand side of Eq. (3).Here, the components of h satisfy the fluctuation-dissipation theorem 63 , In the numerical simulation, the random field is given by where the time increment ∆t is 1 ps in this work.White noise ξ k is defined from two random numbers, ).We added Eq. ( 5) to the magnetic field and solved the LLG equation numerically using the Runge-Kutta method for the investigation of thermal equilibrium before applying microwave voltage.The value of the net perpendicular magnetic anisotropy field, , in the absence of microwave voltage is assumed to be 6.283 kOe, while the saturation magnetization is set to be 955 emu/cm 329 .The temperature T is 300 K, and the volume is V = π × 50 × 50 × 1.1 nm 3 , which is typical for VCMA experiments.The thermal fluctuation excites a small-amplitude oscillation of the magnetization around the energetically minimum state with the ferromagnetic resonance frequency.We pick the temporal directions of the oscillating magnetization and use them as the natural initial states.It should be noted that the value of H Ka at which the complex but still periodic oscillation appears weakly depends on the choice of the initial state, despite the structure of the bifurcation diagram is not changed qualitatively.It might relates to the presence of two stable phases of the parametric oscillation with respect to the microwave voltage 35 .

Analytical treatment of the LLG equation
Here, we discuss an analytical treatment of the parametric oscillation, which provides a condition to excite the oscillation.It is, however, not applicable to investigate the bifurcation from the simple parametric oscillation to the complex but still periodic oscillation shown in Figs.1(b) and 2(a).Since we are interested in the oscillation around the x axis, we introduce angles Θ and Φ as m = (cos Θ, sin Θ cos Φ, sin Θ sin Φ).For simplicity, we introduce notations k = γH Ka and h = γH appl .The LLG equation for Θ and Φ are given as dΘ/dt = − i.e., Since we are interested in an oscillation where m y and m z oscillates with the frequency f L = f /2, we introduce In the parametric oscillation states, the tilted angle Θ of the magnetization from the x axis and the phase difference Ψ are approximately constants 35 .Therefore, after averaging Eqs. ( 7) and ( 9) over a period 1/ f L , we obtain, where σ = h − 2π f L .Although the present work focuses on the case of σ = 0 only, we keep σ as finite here, for generality.The steady state solutions of Θ and Ψ after averaging are These solutions imply, for example, that for exciting the parametric oscillation, which can easily be satisfied when α ≪ 1.These solutions, however, cannot be used to discuss, for example, the bifurcation from the simple parametric oscillation to the complex oscillation because Φ above is assumed to oscillates with only a single frequency f L , while the complex oscillation is a superposition of multiple frequencies.Even for the parametric oscillation state, the above solutions may not be quantitative due to the fact that, strictly speaking for example, Θ is not constant.Although a reliability of the magnetization switching by the parametric oscillation is not the main scope of this work, let us briefly provide a comment on the relationship between the switching reliability and the value of H Ka /H appl mentioned in the introduction.First, the value of H appl should be carefully chosen.For a large H appl , the precession frequency of the magnetization becomes high.Such a fast precession makes the switching unstable because the dynamics becomes sensitive to the pulse shape, as mentioned in the introduction.When H appl is small, however, the switching is also unstable because the dynamics is highly affected by thermal fluctuation, which is also written in the introduction.Summarizing them, the value of H appl should be carefully determined, depending on the various factors, such as circuit quality controlling the pulse shape and the volume of the ferromagnetic layer.Second, a large H Ka is considered to be preferable for a reliable switching due to the 9/14 following reasons.First, as mentioned below Eqs. ( 12) and ( 13), there is a threshold value of H to excite the parametric oscillation.Second, Eq. ( 12) indicates that the cone angle Θ of the magnetization precession around H appl becomes close to 90 • when H Ka is large.In other words, the cone angle decreases as H Ka decreases.A precession with a small cone angle is unstable because such a small oscillation is easily disturbed by thermal fluctuation.Regarding these points, it will be preferable to increase H Ka /H appl , or H Ka , for a reliable switching because there is a threshold of H Ka of the parametric oscillation and it is necessary to make the oscillation robust against thermal fluctuation.However, as revealed in the main text, the precession becomes complex when H Ka /H appl becomes greatly large, which is a main message in this work.

Evaluation method of Lyapunov exponent
The evaluation method of the Lyapunov exponent is as follows 34,64 .We denote the solution of the LLG equation at a certain time t as m(t).We introduce the zenith and azimuth angles, θ and ϕ, as m = (m x , m y , m z ) = (sin θ cos ϕ, sin θ sin ϕ, cos θ ).We also introduce m (1) (t) = (sin θ (1) cos ϕ (1) , sin θ (1) sin ϕ (1) , cos θ (1) ), where θ (1) and ϕ (1) Note that ε corresponds to the distance between m(t) and m (1) (t) at t in the spherical space.Since the Lyapunov exponent characterizes the sensitivity of the dynamical system to the initial state, we study an expansion of ε with time.For this purpose, we assume a small value of ε, which is ε = 1.0 × 10 −5 in this work.For convenience, we introduce a notation, Solving the LLG equations of m(t) and m (1) (t), we obtain m(t + ∆t) and m (1) (t + ∆t), where ∆t is time increment.From them, we evaluate the distance between m(t + ∆t) and m (1) Then, a temporal Lyapunov exponent at t + ∆t is given as where D (1) = D[m(t + ∆t), m (1) (t + ∆t)].

Figure 1 .
Figure 1.(a) Schematic illustration of magnetization oscillation in a ferromagnetic/nonmagnetic trilayer.The unit vector m pointing in the magnetization direction in free layer shows a sustainable oscillation around an external magnetic field H appl in the in-plane direction when the frequency of a microwave voltage is twice the Larmor frequency f L .(b) Dynamical trajectory of the magnetization in a parametric oscillation state.Parameters are H appl = 300 Oe and H Ka = 100 Oe.The black triangle indicates the direction of the magnetization motion.(c) Time evolution of m z in the parametric oscillation state.(d) Fourier spectrum of m z .The inset shows the spectrum around the main peak in a linear scale.diagrams show broad distributions.It is revealed that the appearance of these complex, broadened structures indicated chaotic or transient-chaotic behavior, which are confirmed from evaluations of the Lyapunov exponent and analyses on temporal dynamics.

Figure 2 .
Figure 2. (a) Dynamical trajectory of the magnetization in a parametric oscillation state.Parameters are H appl = 300 Oe and H Ka = 620 Oe.The black triangle indicates the direction of the magnetization motion.(b) Time evolution of m z in the parametric oscillation state.(c) Fourier spectrum of m z .
(b).We observe a clear change of the magnetization dynamics by comparing Figs.1(b) and 2(a).The trajectory is not a simple circle in Fig.2(a).The time evolution of m z and its Fourier transformation are shown in Figs.2(b) and 2(c), respectively.Figure2(b)indicates that the magnetization dynamics is still periodic, while Fig.2(c) indicates the appearance of multipeak structure.These results

Figure 3 .
Figure 3. (a) A bifurcation diagram summarizing local maxima of m z (t) and (b) Lyapunov exponent for various H Ka , where H appl = 300 Oe.

Figure 4 .
Figure 4. (a) Dynamical trajectory of the magnetization in a parametric oscillation state.Parameters are H appl = 300 Oe and H Ka = 1200 Oe.(b) Time evolution of m z in the parametric oscillation state.(c) Fourier spectrum of m z .

Figure 5 .
Figure 5. (a) Dynamical trajectory of the magnetization in a parametric oscillation state.Parameters are H appl = 300 Oe and H Ka = 1600 Oe.The black triangle indicates the direction of the magnetization motion.(b) Time evolution of m z in the parametric oscillation state.(c) Fourier spectrum of m z .

Figure 6 .
Figure 6.Time evolution of m z near t = 0 for (a) H Ka = 620 Oe and (b) H Ka = 1600 Oe.
ζ a and ζ b , in the range of 0 < ζ a , ζ b < 1 by the Box-Muller transformation as ξ a = −2 ln ζ a sin(2πζ b ) and ξ b = −2 ln ζ a cos(2πζ b