{'fraction_non_alphanumeric': 0.046478408966047685, 'fraction_numerical': 0.01428414459949456, 'mean_word_length': 4.148493848112007, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 17, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Chaotic Antiferromagnetic Nano-Oscillator driven by Spin-Torque

We theoretically describe the behavior of a terahertz nano-oscillator based on an anisotropic antiferromagnetic dynamical element driven by spin torque. We consider the situation when the polarization of the spin-current is perpendicular to the external magnetic field applied along the anisotropy easy-axis. We determine the domain of the parametric space (field, current) where the oscillator demonstrates chaotic dynamics. Characteristics of the chaotic regimes are analyzed using conventional techniques such as spectra of the Lyapunov exponents. We show that the threshold current of the chaos appearance is particularly low in the vicinity of the spin-flop transition. In this regime, we consider the mechanism of the chaos appearance in detail when the field is fixed and the current density increases. We show that the appearance of chaos is preceded by a regime of quasiperiodic dynamics on the surface of a two-frequency torus arising in phase space as a result of the Neimark-Sacker bifurcation.

We theoretically describe the behavior of a terahertz nano-oscillator based on an anisotropic antiferromagnetic dynamical element driven by spin torque. We consider the situation when the polarization of the spin-current is perpendicular to the external magnetic field applied along the anisotropy easy-axis. We determine the domain of the parametric space (field, current) where the oscillator demonstrates chaotic dynamics. Characteristics of the chaotic regimes are analyzed using conventional techniques such as spectra of the Lyapunov exponents. We show that the threshold current of the chaos appearance is particularly low in the vicinity of the spin-flop transition. In this regime, we consider the mechanism of the chaos appearance in detail when the field is fixed and the current density increases. We show that the appearance of chaos is preceded by a regime of quasiperiodic dynamics on the surface of a two-frequency torus arising in phase space as a result of the Neimark-Sacker bifurcation.

I. INTRODUCTION
Neuromorphic computing is a rapidly developing field inspired by the idea of emulating biological processes in the brain [1,2]. Since neurons demonstrate a rhythmic activity, non-linear nano-oscillators of different physical nature [3][4][5] are considered as candidates for artificial neurons. Spin-torque nano-oscillators are of special interest because their dynamical regimes are easy-tunable by means of electrical current. They are used to emulate single neurons [4], as well as small neural networks [6].
Several concepts of neuromorphic computing rely on stochastic (chaotic) oscillator behavior; examples are reservoir computing [4,[7][8][9][10] or spike-based encoding [2,[11][12][13]. Thereby, a spin-torque nano-oscillator that demonstrates controllable chaotic dynamics is of high importance for these concepts. The simplest spin-torque nano-oscillator is made of a collinear antiferromagnet (AFM) and is essentially a non-linear dynamical system with a four-dimensional phase space [14,15]. This feature makes an AFM nano-oscillator a natural candidate for the realization of chaotic dynamic regimes. The possibility of chaos in such systems was recently pointed out in Ref. 16. An advantage of AFM nano-oscillators is their fast dynamics in the THz regime. This may allow for the construction of ultrafast artificial neurons [11].
Here, we provide a general analysis of possible dynamical regimes of an AFM nano-oscillator driven by a spincurrent. Two external parameters control the system, namely the applied magnetic field and the strength of the current. The dynamical regimes were studied in a large domain of the 2D parametric space, which is of interest for spintronic applications. First, we formulate the model and derive the basic equation of motion, see Sec. II supplemented with App. A. Next, we consider possible fixed points of the system and study their stability, see Sec. III and App. C. Then we analyze the regimes in phase space where stable limit-cycles are present, and we provide analytical expressions for the parameters of the limit cycles and analyze their stability, see Sec. IV and App. D. In the last Sec. V we demonstrate the possibility of chaos and hyperchaos in the system and analyze characteristics of the chaotic dynamics within parameter space. Also, we identify the mechanism of the chaos appearance in the vicinity of the spin-flop transition.

II. MODEL
We consider a two-sublattice antiferromagnetic (AFM) film. The magnetization of each of the sublattices is modeled by means of a continuous function M i = M i (r, t) of constant amplitude |M i | = M s , where i = 1, 2. In the following, it is instructive to introduce the unit magnetization vector µ i = M i /M s together with the ferrom = (µ 1 +µ 2 )/2 and antiferromagnetic n = (µ 1 −µ 2 )/2 vector order parameters. Note that m · n = 0 and m 2 + n 2 = 1.
The magnetization dynamics is described by the set of Landau-Lifshitz equations supplemented with the Gilbert damping and the spintorque term in the form proposed by Slonczewski [19][20][21]. Here, γ > 0 is the gyromagnetic ratio, and α is the damping constant. The strength of the spin torque is determined by the constant σ = γ P J/(2|e|LM s ), where denotes the Planck constant, e is the electron charge, L is the thickness of the AFM film, J > 0 and 0 < P < 1 are the density and the rate of the spin-polarization of the current which flows perpendicularly to the film. The geometry of the studied system is represented by a pillar structure that consists of three layers: ferromagnetic polarizer, nonmagnetic spacer, and the studied AFM free layer (see Fig. 1(a,b)). The unit vector p indicates the direction of the spin polarization. Two equations of the set (1) are coupled by means of the energy functional E = E[µ 1 , µ 2 ] which includes intra-and inter-lattices exchange interactions (A1), uniaxial anisotropy with the easy-axis oriented perpendicularly to the AFM film (A3) and interaction with an external magnetic field B (A4). In the following, we consider the case when the applied field and the anisotropy field B an are much smaller compared to the exchange field B ex acting between the sublattices. The fields B ex and B an are introduced in (A1) and (A3), respectively. In this limit, the dynamics of the Néel vector n is governed by the Lagrangian [14,22] supplemented with the Rayleigh dissipation function of density [14] R where the constraint |n| = 1 is presumed. Thus, the dynamical system under consideration is controlled by three parameters: the dimensionless magnetic field b = B/B sf in units of the spin-flop field B sf = √ B ex B an , the dimensionless current density j = P J/J 0 with J 0 = B an |e|LM s / , and the modified damping coefficientᾱ = α B ex /B an . In (2) we utilized the dimensionless space-time coordinates = r/ and τ = tω 0 , where = A/(B an M s ) is the exchange length and ω 0 = γB sf is the frequency of the uniform AFM resonance. Here, A is the constant of the nonunifom exchange interaction introduced in (A2). This dimensionless system of units was used throughout the following part of the main text.
Having determined the Néel vector n by means of the theory (2), one obtains the magnetization vector in the form [14,22] In the following, we focus on the case p = −ŷ and b = bẑ. The spatially uniform solutions n = n(τ ) of (2) are determined by the equation It is convenient to enforce the constraint |n| = 1 using the spherical parameterization n = sin θ (x cos φ +ŷ sin φ) + cos θẑ. In this case, Eq. (4) reduces tö The set of nonlinear equations (5) is the main subject of this paper.

III. STATIC EQUILIBRIUM STATES
Let us start with static equilibrium solutions. For the case j = 0, the considered AFM system is well studied [23]. It is known that depending on the applied magnetic field, there are two equilibrium states: one state with uniform polarization along the anisotropy axis n||ẑ for 0 < |b| < 1 (solution A), and the spin-flop state n ⊥ẑ when |b| > 1 (solution B). In the first and latter cases the equilibrium orientation of the Néel vector is doubly and continuously degenerate, respectively. Linear eigenexcitations of both solutions are characterized by two frequency modes. The high-and low-frequency modes of solution (A) are, ω a h = 1 + b and ω a l = 1 − b, respectively. They are plotted as dashed orange lines in Fig. 1d  Applying a spin-torque changes the properties of both states (A) and (B). For state (A), the magnetizations of the sublattices take a turn in the plane perpendicular to p such that the Néel order parameter n makes an angle 0 < θ A < π/4 with the anisotropy axis, see Fig. 1a). Here, with j c = 1−b 2 satisfies static form of Eqs. (5). Note that one more solution can be obtained from (6) by means of the transformation θ → π − θ, φ → φ + π (equivalently n → −n) that also satisfies (5) and has the same properties as (6). State (A) is stable within the parameter domain b < 1 and j < 1 − b 2 , see Fig. 1c). The spin torque removes the degeneracy of the high-and lowfrequency modes at zero magnetic fields, resulting in the splitting ∝ j 2 . Namely, ω a h = 1 − 1 8 j 2 and ω a l = 1 − 1 4 j 2 if b = 0. The spin torque also reduces the critical field b c1 = √ 1 − j when the state (A) becomes unstable. In the vicinity of the instability the high-and low-frequency modes demonstrate the following asymptotic behavior ω a h ≈ c 1 +c 1 Note that Eqs. (5) allow one more static solution (A'), which one obtains from (6) by means of the transformation θ a = θ a + π/2, φ a = φ a + π. It is unstable for all values of b and j (see App. C for details) and is not considered in the main text.
The spin-torque removes the continuous degeneracy of the state (B) and leads the sublattices magnetization to orient itself parallel to the direction of the spin polarization: n = ±p. In terms of the angular parameterization  Fig. 1d). The spin-torque increases the minimal magnetic field b c2 = √ 1 + j above which the state (B) is stable. So, a small spin-torque opens a gap b c1 < b < b c2 in the vicinity of the spin-flop field, where no static solitons can exist. As it will be shown latter, chaos appears mostly in this gap.
The magnetization of the static state (B) is The analysis provided in this section was made assuming vanishing damping. Details of the stability analysis for possible static states, including the damping effects, are provided in Appendix C.

IV. REGULAR DYNAMICS AND TRANSITION TO CHAOS
Let us start with the simple case b = 0. Increasing the current from j = 0 to j = 1, we obtain the static stable solution (A) shown in Fig. 1(a) whose inclination angle varies from θ a = 0 to θ a = π/4. With a further increase in the current j > 1, equations (5) do not allow stable static solutions any more. However, in the particular case b = 0, Eq. (5b) is satisfied if φ = 0. This means that the dynamics of the Néel vector n is constrained within z − x plane. The orientation of the vector n is determined by the polar angle θ(τ ) which satisfies the equation of motion of a driven nonlinear pendulum with dampingθ +ᾱθ + sin θ cos θ = j/2.
The dynamics of the antiferromagnetic oscillator for this case was analyzed in detail in Ref. 16 and will not be considered here. We only note that for large currents j 1 the dynamics of n takes the form of a quasiuniform rotation θ ≈ Ωτ with Ω = j/(2ᾱ). On a unit sphere, this dynamics is represented by a circular limit cycle in x − z-plane, see Fig. 2a.
Now we introduce a small magnetic field b 1. This limit allows us to assume thatφ θ . I.e. the system has two well distinguished time scales, namely, vector n rapidly rotates with the angular velocityθ within a vertical plane, which slowly rotates about the z-axis with the angular velocityφ. This kind of dynamics is shown in Fig. 2b. The corresponding approximate analytical solution can be derived from (5) in the following way: Assuming that the current is large we conclude from (5a) that θ ≈ Ωτ with Ω = j cos φ/(2ᾱ). Substituting this approximation into (5b) and neglecting terms involving φ andᾱφ we obtain a solution for φ(τ ) in the implicit form The solution (10) satisfies the initial condition φ(0) = 0, which corresponds to the static state (6). Thus, if one starts with the static state (A) and increases the current step-by-step such that j > j c , then the resulting dynamics can be approximated by (10). Indeed, as it can be seen in Fig. 2g, in the limit of small fields, the approximation (10) demonstrates a good agreement with the numerical solution of (5).
As it follows from (10), for τ → ∞ the orientation of the plane of rotation of n approaches the asymptotic value φ → φ 0 = arctan(1/η). The corresponding limit cycle on the unit sphere has the following characteristics The typical time of approaching the limit cycle is ∆τ ≈ 2ᾱ/(ᾱ 2 + 4b 2 ). Note that for larger fields, the plane of rotation of vector n approaches its asymptotic orientation φ 0 in an oscillatory manner, see solid lines Fig. 2g. This behavior is not described by the approximation (10) because terms involvingφ in (5b) were neglected during the derivation of (10). Despite this, the asymptotic value φ 0 is valid for relatively large fields, see Fig. 2e. Interestingly, in the particular case b = 1, the limit cycle (11) is an exact solution of Eqs. (5). Using the monodromy matrix technique (see App. D) we found that this solution demonstrates a couterintuitive stability up to relatively high currents.
Experimentally, the described dynamics can be detected by observing the oscillations of the magnetization vector m, which are measurable in contrast to n. Using (3) one easily obtains that for the limit cycle (11) the magnetization vector reads where ε =x cos φ 0 +ŷ sin φ 0 . Thus, the presence of the doubled frequency 2Ω 0 in the spectra of electromagnetic radiation of the nano-oscillator indicates the presence of the limit cycle (11). This is reflected by the Fourier spectra in Fig. 2(a,b). Although the oscillatory part of the magnetization vector (3) vanishes for b = 0, the spectral line at 2Ω 0 is still present, see Fig. 2a. This is because the actual solution θ(τ ) of the pendulum equation Eq. (9) deviates from the uniform rotation θ ≈ Ω 0 τ .
As it follows from (11), the frequency Ω 0 decreases with the field and soθ also decreases. This results in a violation of the assumption about the separation of time scalesθ φ . It leads to a significant deviation of the trajectories on the unit sphere from the approximation (10). The resulting dynamics can be still regular (Fig. 2c) or chaotic (Fig. 2d). The appearance of chaos drastically changes the spectra of the magnetization, see Fig. 2d, and therefore can be experimentally detected using spectral analysis of the radiation emitted by the oscillator.
To track the appearance of chaos in our system, we employed a method similar to the use of Poincaré sections: We consider the intersection points of the trajectories on the unit sphere with the horizontal half-plane {z = 0, y > 0} for times τ > τ * where time τ * is much larger than all typical timescales of the system. The regular dynamics is characterized by a constant number of intersection points for all times, e.g. one point for Figs. 2a,b and three points for Fig. 2c. For chaotic or quasiperiodic behavior, the number of intersection points increases in time and becomes quite large for long times of integration of (5), see Fig. 2d. Since an intersection point lies on the circle of unit radius, it is characterized by a single polar coordinate φ. For the circular limit cycles described above φ = φ 0 . So we use this notation for the intersection point coordinates even if the limit circle is more complicated (Fig. 2c) or it does not exist (Fig. 2d). Plotting all coordinates φ 0 obtained for a range of magnetic fields, we build a typical bifurcation diagram shown in Figs. 2e,f. As one can see, chaotic regimes are intercepted by windows of regular dynamics characterized by period multiplication of the limit cycles. Note that this technique does not allow one to distinguish between chaotic and quasiperiodic dynamics. However, as it will be shown in the next section, the domains of quasiperiodicity are quite narrow.

V. CHARACTERISTICS OF THE CHAOTIC REGIME
Spectra of the Lyapunov exponents (LE) are a universal tool for studying the complicated dynamics of ergodic systems [24,25]. In order to utilize this technique we consider our system (5) as a flowẋ = f (x), with x = (θ,θ, φ,φ) ∈ R 4 being a time-dependent vector of 4-dimensional phase space. The phase space is equipped with the Euclidean metric and norm ||x|| = i . The latter enables us to introduce the distance ∆(τ ) = ||x(τ ) −x(τ )|| between two trajectories x andx. The nature of the long-term evolution of the distance ∆(τ ) for the case ∆(0) → 0 characterizes the instability of the trajectories with respect to the initial conditions. The rate of divergence of the trajectories is determined by the largest LE [25] which is shown in Fig. 3. The value of Λ 1 enables one to deduce a general features of the dynamics, e.g. for the case Λ 1 > 0 the system dynamics is chaotic. However, the complete characteristics of the dynamics of a dissipative system is determined by the types of attractors present in phase space. An attractor can be identified from the full spectrum of the LEs whose number coincides with the dimension d = 4 of the phase space [26].
The key instrument here is Oseledets multiplicative ergodic theorem [27] which guaranties that all trajectories belonging to one attractor have the same LE spectrum. This means that any trajectory can be chosen for determining the LEs, in the other words, the LEs are independent on the initial conditions. Using the algorithm proposed by Benettin at al. [28,29] we calculated all four LEs for a wide range of parameters b and j, see Fig. 3. For all points of the studied parameter space the sum of all LEs is negative [30]. This reflects the dissipative nature of our system, since negative S d = ∇ · f (x) τ indicates shrinking volumes in phase space [25]. Thus, we can characterize the dynamics of the system by the type of attractors formed in phase space.
The signs of the ordered LEs compose the signature of the attractor. Using the obtained maps of LEs we distinguish the following signatures: −, −, −, − -stable fixed point corresponding to a static equilibrium state, 0, −, −, − -stable limit cycle corresponding to periodical dynamics, +, 0, −, − -chaos, +, +, 0, − -hyperchaos. The fixed points are analyzed in Sec. III, and the corresponding magnetic structures are shown in Fig. 1ab. In the map for Λ 1 (Fig. 3) the parameter space domains which correspond to the fixed points A and B are marked in blue and are labeled "FP a " and "FP b ", respectively. The white areas of the map for Λ 1 correspond to limit cycles and are labeled as "LC". Some examples of such limit cycles are studied in Sec. IV, they are illustrated in Figs. 2abc, see also insets I and II in Fig. 5. Note that the line {j = 0; b > 1} can be considered as a limit cycle of zero frequency. It is formed due to the continuous degeneration of the fixed point B with respect to the rotation of the vector n within x − y-plane.
Red areas in the map for Λ 1 indicate chaotic dynamics; they are labeled as "C". Interestingly, in the vicinity of the spin-flop b 1, the dynamics is mostly regular. This correlates with a high stability of the limit cycle (11) for the case b = 1, see App. D for details. Also, as it follows from the shape of the LC-domain for lower fields, the current supports the stability of the limit cycle. Finally, the chaotic region is bounded for b < 1. This is in contrast to the case b > 1, where the current supports the emerging chaos.
Note that at the boundaries between the periodic and chaotic regimes, there are quite narrow regions with signatures 0, 0, −, − . They correspond to quasiperiodic dynamics on the surface of a 2-frequency tori in the phase space. Due to the narrowness, these quasiperiodic regions are not marked on the diagram in Fig. 3. Their investigation requires a significant increase in resolution of the LE maps and, therefore, high computing costs. Here, we restrict ourselves to only one example explored in Fig. 5, where the transition to chaos goes through a relatively large region of quasiperiodic dynamics.
The maps of LEs were obtained for the phase trajectory with initial conditions These initial conditions represent a possible experimental realization when the AFM sample is relaxed in the given magnetic field, and then the current is switched on step-by-step. If there is a single attractor in the phase space, then the choice of the trajectory (initial conditions) does not matter because, according to Oseledets theorem, all trajectories belonging to the same attractor have identical spectra of LEs. However, in our system, we found regimes when two attractors coexist, namely fixed point of type A coexists with limit cycle and with a strange attractor in regions "LC/FP a ", and "C/FP a ", respectively, see Fig. 3. In these regimes, the long-term dynamics is sensitive to the initial conditions. We illustrate this by plotting the basins of attraction of the fixed point of type A for different regimes, see Fig. 4. As it follows from panel a), one can reach the fixed point only by starting in the vicinity of this point (gray region). White points correspond to the attraction basin of the limit cycle. Note that the complete 4D attraction basins are path-connected, in contrast to the 2D-cross-sections presented here. The complicated structure of the attraction basin of the strange attractor (white area of panel b) allows us to conjecture that the boundary of the attraction basins has a fractal structure in 4D phase space [31]. The strange attractor is shown in Fig. 4c) by means of its natural measure, which is the probability density to find the system in a given point of the phase space. Naturally, the highest probability density is in the vicinity of the fixed points.
In the map of Λ 2 one can distinguish a small region of hyperchaos where Λ 2 > 0, which is denoted as "HC". For example for the point (b, j) = (1.18, 0.54) we obtained (Λ 1 , Λ 2 , Λ 3 , Λ 4 ) = (0.017, 0.009, 0, −0.226) with an error smaller than 10 −3 . In order to better understand how chaos develops in this regime, we explore the behavior of the Lyapunov exponents along the line b = 1.18 -the vertical dashed line in Fig. 3. The development of chaos with increasing current is illustrated in Fig. 5 by analyzing the three largest LE. The smallest LE is always negative and not important for the analysis. For small currents j < j c ≈ b 2 − 1 all LEs are negative, which corresponds to the fixed point 'B' illustrated in Fig. 1b). In the point '1', the two largest LE Λ 1 and Λ 2 are equal to each other and reach the value zero. As the current increases further, we obtain Λ 1 = 0 and Λ 2 < 0. This is the scenario where a soft-born, stable limit cycle (LC) emerges via an Andronov-Hopf bifurcation. When the current increases even further, the LC continuously increases its size and chances its form and position on the surface of the unit sphere, see inset I. In point '2', the LC experiences a jump-wise change of its form and position, see inset II. This is not a period-doubling bifurcation. After this, the two LEs Λ 2 and Λ 3 are equal to each other, reaching the value zero at point '3'. As the current keeps increasing, we have Λ 1 = Λ 2 = 0 and Λ 3 < 0. This is the scenario where a soft-born, 2-frequency torus (T) emerges via a Neimark-Sacker bifurcation (NSB) [32,33]. An example of the quasiperiodic dynamics on the torus surface is demonstrated in inset III by plotting the n-vector trajectories obtained for the time interval ∆τ = 300. In point '4' the torus changes its form and in point '5' the torus is destroyed, and chaos (C) appears, featuring a single positive LE. The latter process has a fine structure: the torus experiences synchronization, resulting in the creation of a LC, from which a subsequent torus emerges via a second NSB. In point '6', chaos is replaced by a periodic window, where a LC is present. In point '7' a 2-frequency torus bifurcates from the LC via NSB. With a small current increase, the torus is destroyed, and hyperchaos (HC) with two positive LEs develops. An ex- ample for hyperchaotic behavior is demonstrated on inset IV. Remarkably, a practically identical mechanism for the appearance of hyperchaos was recently reported for a Van der Pol oscillator with feedback loops [32]. A similar mechanism was also found for a system of two coupled Rößler oscillators [34]. All the previous analysis was made for a constant dampingᾱ = 0.1. In order to analyze how damping influences the chaotic dynamics, we consider the dependence of the Lyapunov dimension D l of the strange attractor onᾱ for the cases b < 1 and b > 1, see Fig. 6. The Lyapunov dimension is used to approximate the informational dimension D 1 using the formula of Kaplan and Yorke: where k : S k ≥ 0 and S k+1 < 0. In both cases, the behavior D l (ᾱ) is the same ifᾱ is small: for finite but vanishinḡ α one has D l d and both attractors demonstrate the same rate of decrease of D l with increasingᾱ. However, for larger damping the behavior is different. In the case b < 1, chaos is suppressed already forᾱ ≈ 0.2 when the chaotic attractor tranforms into a limit cycle. In the case b > 1, chaotic dynamics -intercepted by periodic windows -exists for much larger dampings, up toᾱ ≈ 0.8. In the latter case, the strange attractor becomes almost flat.
Chaos characterization by analyzing the magnetization. The observation of the dynamics of the Néel vector is a challenging experimental problem. On the contrary, the magnetization dynamics can be easily observed by detecting the generated ac magnetic field. For this reason, it is instructive to explore the relation between the given regime of the system dynamics and the behavior of the magnetization vector (3). In Fig. 2a-c we demonstrated how the Fourier spectra of a magnetization component can be used to distinguish regular and chaotic dynamics. Here, we analyzed the behavior of the magnetization vector for a wide range of parameters space, see  (12) and (11).
the applied magnetic field (see Sec. III), wherein |m b | is noticeably larger and exactly parallel to the applied magnetic field. The time-averaged magnetization generated by the low-field periodic dynamics is determined by the first summand in (12). It is oriented within the plane perpendicular to the field. In chaotic regimes, the magntization dynamics is rather complicated, see Fig. 7a). Remarkably, the time-averaged magnetization m τ is determined only by the control parameters b, j, α, and it is a property of the attractor realized for the given dynamics, see Fig. 7b)-e). The spatial octant, in which vector m τ lies is determined by vectors b and p, namely sign( m x τ ) = sign(bp), sign( m y τ ) = sign(p), and sign( m z τ ) = sign(b). Here, vectors b and p are applied along axesẑ andŷ, respectively.

VI. CONCLUSIONS
The simultaneous action of magnetic field and spintorque on AFM nano-oscillator can result in complicated nonlinear dynamics of the AFM order parameter. Depending on the control parameters (b, j), the system demonstrates different types of nonlinear behavior: regular dynamics along limit circles, quasiperiodic dynamics along 2-frequency tori, as well as chaotic dynamics. The latter can be of interest for neuromorphic computing relying on stochastic elements. We show that the threshold current of the appearance of chaos is especially low in the vicinity of the spin-flop transition. This Appendix aims to introduce the number of material parameters e.g. exchange B ex and anisotropy B an fields, and also to provide the derivation of the equations of motion in exchange approximation B ex B an . The latter derivation was previously made in Ref. 14, and we provide it here for the sake of text coherence. For the derivation without spin torques, see also Refs. 35-37. We assume that the energy density of the system can be presented in the form E = E ex + E an + E zee . Here, the first summand represents the exchange contribution with density [23,38] which includes uniform exchange between the sublattices (B ex ) and isotropic nonuniform intra-(A ) and inter-(A ) lattices exchange interactions. In terms of the vectors n and m the exchange energy density (A1) has the form The energy of uniaxial anisotropy E an has the density [23,38] where we assumed that the anisotropy easy-axis is oriented alongẑ, and also introduced the anisotropy fields B an = (K − K )/M s andB an = (K + K )/M s . Finally, the interaction with an external magnetic field B is In terms of the vectors n and m the set of Landau-Lifshitz equations takes the following forṁ Based on the form of the energy functional E, we conclude that in the static case |m| ∼ B/B ex 1 for a large exchange field B ex B,B an . In this limit, we neglect in Eq. (A5a) all terms proportional to m except the leading linear term ∝ B ex m. From this we find thaṫ n ≈ γn × [mB ex − B], where it was also assumed that B ex B an . Using this in the considered limit |n| ≈ 1, we solve the latter equation with respect to m: The equations of motion (5) can be rewritten as a first-order differential problemẋ = f (x) with x = (θ,θ, φ,φ) t . There are six different fixed points x * : f (x * ) = 0, which correspond to three distinct physical states, since n → −n. The three static, uniform solutions are visualized in Fig. 8: for finite field b and current j solutions (A) and (A ) are tilted off the z-and x-axis, respectively, by an angle θ A . At the critical current j c , both fixed points collide and annihilate in a saddle-nodebifurcation at θ A = π/4. Solution (B) persists for all parameters and is always oriented parallel/anti-parallel to the direction of spin polarization p.
The stability of a fixed point x * can be determined by solving the linearized systemẋ =Lx withL = ∂f ∂x x=x * being the constant matrix andx = x − x * are small perturbations. Solutions have the formx(τ ) =x 0 e iωτ , where ω = −iλ and λ is an eigenvalue ofL.
For solution A the eigenvalues are cumbersome expressions, but for solution B, they can be determined more where the eigenfrequencies ω b h and ω b l are defined in Sec. III. The evolution of ω with b for finite dampinḡ α is shown in Fig. 8.
For solution (A)ᾱ = 0.5 was chosen in order to feature the splitting in the imaginary part. For b < b c1 , the two branches Im(ω) are both positive and they separate linearly with b. For b > b c2 , the splitting is still linear in b, but now one of the branches is negative, indicating the instability of solution (A). This branch is negative only for finite dampingᾱ, a phenomenon which is referred to as "dynamic instability" and which is known also e.g. from structural mechanics [39].
For solution (A ) there is always a negative branch of the imaginary part of ω, thus solution (A ) is never stable.
For b < b c2 there is always at least one negative branch of Im(ω) of solution (B) and thus it is unstable. At b c1 pairs of eigenvalues collide on the imaginary axis in a reversible Hopf-bifurcation. This is analogous to the upward Ziegler double pendulum, which displays also a twoby-two reversible Hopf-bifurcation, depending on the follower load [39,40]. In the field range b c1 < b < b c2 solution (B) experiences the flutter instability. Solution (B) gains stability at b c2 (ᾱ). For b > b c2 there is at first a single, positive branch of Im(ω), which is then splitting at b c3 = 1 +ᾱ 4 +4j 2 4ᾱ 2 into three branches. Here, b c3 indicates the threshold, where the in-plane, linear oscillations become overdamped, as the corresponding branch of the real part becomes zero at b c3 .
Appendix C: Stability analysis of static fixed points.
The equations of motion (5) can be rewritten as a first-order differential problemẋ = f (x) with x = (θ,θ, φ,φ) t . There are six different fixed points x * : f (x * ) = 0, which correspond to three distinct physical states, since n → −n is a symmetry. The three static, uniform solutions are visualized in Fig. 8: for finite field b and current j solutions (A) and (A ) are tilted off the z-and x-axis, respectively, by an angle θ A . At the critical current j c , both fixed points collide and annihilate in a saddle-node-bifurcation at θ A = π/4. The solution (B) persists for all parameters and is always oriented parallel/anti-parallel to the direction of spin polarization p.
The stability of a fixed point x * can be determined by solving the linearized systemẋ =Lx withL = ∂f ∂x x=x * being the constant matrix andx = x − x * are small perturbations. The solutions have the formx(τ ) =x 0 e iωτ , where ω = −iλ and λ is an eigenvalue ofL.
For solution A the eigenvalues are cumbersome expressions, but for solution B, they can be determined more easily fromL where the eigenfrequencies ω b h and ω b l are defined in Sec. III. The evolution of ω with b for finite dampinḡ α is shown in Fig. 8.
For solution (A)ᾱ = 0.5 was chosen in order to feature the splitting in the imaginary part. For b < b c1 , the two branches Im(ω) are both positive and they separate linearly with b. For b > b c2 , the splitting is still linear in b, but now one of the branches is negative, indicating the instability of solution (A). This branch is negative only for finite dampingᾱ, a phenomenon which is referred to as "dynamic instability" and which is known also e.g. from structural mechanics [39].
For solution (A ) there is always a negative branch of the imaginary part of ω, thus, solution (A ) is never stable.
For b < b c2 there is always at least one negative branch of Im(ω) of solution (B) and thus it is unstable. At b c1 pairs of eigenvalues collide on the imaginary axis in a reversible Hopf-bifurcation. This is analogous to the upward Ziegler double pendulum, which displays also a twoby-two reversible Hopf-bifurcation, depending on the follower load [39,40]. In the field range b c1 < b < b c2 solution (B) experiences a flutter instability. It gains stability at b c2 (ᾱ). For b > b c2 there is at first a single, positive branch of Im(ω), which is then splitting at b c3 = 1 +ᾱ 4 +4j 2 4ᾱ 2 into three branches. Here, b c3 indicates the threshold, where the in-plane, linear oscillations become overdamped, as the corresponding branch of the real part becomes zero at b c3 .
FIG. 9. Evolution of the Floquet multipliers (absolute value |λi| of the monodromy matrix) for the limit cycle at b = 1.0 and α = 0.1. The limit cycle is remarkable stable with just one island of instability 1.615 < j < 1.824, where the greatest Floquet multiplier is larger than one. There are only two more areas of instability for larger current: 2.366 < j < 3.073 and 4.265 < j < 6.953. sinθ sinφx + cosφẑ + cosθŷ is the Néel vector in the rotated frame of reference. The exact solution forθ(τ ) andφ(τ ) can be obtained from the relationñ = n, where n = sin(Ω 0 τ )(cos φ 0x + sin φ 0ŷ ) + cos(Ω 0 τ )ẑ.
The monodromy matrixM =Φ(T ) corresponds to the fundamental matrixΦ of the linearized dynamical systeṁ x =Lx, evaluated at the orbit period T = 2π/Ω 0 of the limit cycle. In order to determine it, the matrixL = ∂f ∂x x=x cycle is evaluated at the exactly know solution of the limit cycle, in the rotated frame of reference. Next, the linearized system is solved numerically, starting from Φ(0) =1, for one orbit period T .
The absolute value of the eigenvalues of the resulting monodromy matrixM are the Floquet multipliers; their evolution is shown in Fig. 9 depending on the current j.
The limit cycle is remarkable stable, having only three instability islands, where the greatest Floquet multiplier is larger than one, see Fig. 9. The first instability island 1.615 < j < 1.824 is consistent with a small island of chaos, where the greatest Lyapunov exponent is greater than zero in Fig. 3. There are only two more areas of instability for larger current: 2.366 < j < 3.073 and 4.265 < j < 6.953.
The presence of this stable limit cycle explains why chaotic dynamics is almost absent along the line b = 1 in the Lyapunov spectra, see Fig. 3.