{'fraction_non_alphanumeric': 0.0389426385459256, 'fraction_numerical': 0.011820426668825648, 'mean_word_length': 4.524150268336315, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Analytical investigation of modulated spin torque oscillators in the framework of coupled differential equations with variable coefficients

Modulation of Spin Torque Oscillators (STOs) is investigated by analytically solving the time-dependent coupled equations of an auto-oscillator. A Fourier series solution is proposed, leading to the coefficients being determined with a linear set of equations, from which a Nonlinear Amplitude and Frequency Modulation (NFAM) scheme is obtained. In this framework, the NFAM features are related to the intrinsic STO parameters, revealing a frequency-dependence of the harmonic-dependent modulation index that allows a modulation bandwidth to be defined for these devices. The presented results expose a rich parameter space, where the modulation and the STO's operation conditions define the observed modulation features. The Fourier-series representation of the time signal is suitable for studying periodic perturbations on the auto-oscillator equation.


I. INTRODUCTION
Spin torque oscillators 1 (STOs) are spin torque [2][3][4] driven devices whose high-frequency tunability and nanosized dimensions are promising for diverse technological applications. The oscillators' functionality can be very broad, depending on the specific application and its operation frequency. Two phenomena are of particular use, namely (i) the ability to phase lock or synchronize to another periodic signal, and (ii) the possibility of modulating a base-band signal. 5 In the case of STOs, both phenomena have been experimentally observed.
Phase locking of STOs has been extensively studied in the past few years. Injection locking of an STO to an external source is already well-established in both experimental 6-8 and numerical [9][10][11][12][13][14] works. Phase locking can also be achieved by mutually coupling several STOs. In this case, a distinction must be made between nanocontact 15 and nanopillar 16 geometries. In nanocontacts, synchronization mediated by propagating spin waves 17,18 has been satisfactorily achieved both experimentally [19][20][21] and numerically 22 for GMR-based STOs. In contrast, the mutual locking of nanopillars has been mainly successful in numerical studies. [23][24][25] Modulation of STOs has also been shown to be experimentally straightforward giving rise to several applications. STOs nanopillars have been proposed as hard-drive read-heads where the stray field from the perpendicular media acts a modulating source 26 both in GMR 27 and TMR 28 spin valves. From the perspective of communication applcations, Frequency Modulation (FM) has been observed in GMR nanocontacts [29][30][31][32][33][34] and more recently in nano-oxide layer (NOL) nanocontact STOs. 35 Digital communication schemes, in particular Frequency Shift Keying (FSK), has been demonstrated in vortex-based nanocontact STOs up to the limit of analog FM. 36,37 Similar modulation effects are most likely to be observed in other STO geometries that take advantage of Perpendicular Magnetic Anisotropy (PMA) materials 38 and nanopillars located in a microstrip resonator. 39 In the pioneering frequency modulation experiment, 29 the authors took advantage of the strong phase-power coupling of GMR-based STOs to obtain a frequencymodulated voltage output from a base-band current tone. Although the gross features of frequency modulation were present, significant discrepancies from the model were observed: (i) an unexpected carrier frequency shift as a function of modulation strength, and (ii) asymmetric sideband amplitudes in contrast to the expected n-th order Bessel functions. A mathematical model was proposed by Consolo et al. 40 to fit these features. Here, the temporal signal s(t) was defined to be both amplitudemodulated and frequency-modulated where ω c is the carrier frequency, and the amplitude A and instantaneous frequency ω i are polynomial expansions of a base-band message x(t). This combined action of amplitude and frequency modulation was named "Nonlinear Amplitude and Frequency Modulation" (NFAM). Using this approximation, it was possible with very good accuracy to fit the GMR-STO modulation data with a polynomial expansion up to the third order. 30 However, the information acquired from these coefficients does not reflect the intrinsic mechanism of NFAM in STOs, nor its consequences for a specific choice of experimental conditions and device characteristics. To investigate the origin of NFAM in STOs, we here derive the modulation spectrum as a function of the intrinsic STO parameters. The modulation spectrum of an STO is obtained by perturbing the auto-oscillator general equation 41 with a slow time-varying tone. Such perturbation creates coupled phase and amplitude variations, which lead to NFAM. It is shown that a Fourier series gives an adequate representation of s(t), in which case the problem can be reduced to determining the coefficients from a linear set of equations. Moreover, the Fourier coefficients obtained by this method give quantitative information, such as the carrier frequency shift, the modulation index, FIG. 1. Nanocontact pseudo spin valve consisting of a thick (fixed) and a thin (free) magnetic layer (gray) decoupled by a thin metallic spacer (black). A nanocontact is patterned on top of the free layer, allowing the device to be driven with a spatially confined current. It is assumed that the autooscillator Eq. (2) is a good approximation for such a structure. The free-layer magnetization m is free to precess about the applied field Happ. The fixed layer is assumed to be magnetized predominantly in plane, and its magnetization PS subtends an angle γ0 relative to the free layer magnetization. A steady precession is achieved when a spin polarized current exerts sufficient spin torque on the free layer to compensate for the magnetic damping. and the modulation bandwidth of the STO. This paper is divided as follows: in Section II, the general formulation of the problem and its solution is given. It is shown that a carrier frequency shift appears as a consequence of power and phase coupling. Moreover, the modulation index is obtained by calculating the power spectral density, and shows an unexpected modulationfrequency dependence. In Section III, the analytical solution is evaluated for different modulation and operation conditions, and shows a qualitative agreement with experimental observations. Concluding remarks are given in Section IV.

II. MODULATION OF A NONLINEAR AUTO-OSCILLATOR
In this section, the power spectral density (PSD) of a frequency-modulated nonlinear auto-oscillator is analytically calculated, using the general model proposed by Slavin and Tiberkevich: 41 where p is the oscillation's power, c = √ pe −iφ is the oscillation's amplitude, φ is its time-dependent phase, and ω(p), Γ + (p), and Γ − (p) are respectively the powerdependent oscillation frequency, damping, and negative damping. The power dependencies of these quantities are treated in Ref. 41 as polynomial expansions. We assume that the model of Eq. (2) is a good description of the nanocontact geometry schematically shown in Fig. 1.
Here, a metallic nanocontact is patterned on top of an extended GMR spin valve. A strong magnetic field H app is applied to ensure the saturation of the magnetically active or "free" layer material. When a bias current is driven through the nanocontact, spin torque is exerted on the underlying free-layer area, counteracting the magnetic damping action, and leading to a small-angle precession about the applied field. Depending on the applied field angle, the free layer might exhibit propagating spin waves or a so-called localized bullet. [42][43][44][45] In this paper, we restrict our attention to a perpendicular applied field without loss of generality.
In order to solve the modulation problem, we recall the experiment of Ref. 29 where an ac current is added to the dc bias current. Consequently, we define the slow time-varying current as where µ is the modulation strength, defined as the ratio between the ac amplitude and the dc current, and ω m is the modulation frequency. We assume that ω m ≪ ω ST O (where ω ST O is the STO's free-running frequency) in order to be considered a modulating frequency. On the other hand, if ω m ≈ ω ST O , the STO might instead become injection-locked. The current is thus included in the negative damping parameter, Γ − (p), which can be generally expanded as a polynomial in power where σ = ǫ γ/2eM S V , ǫ the spin-polarization efficiency, γ the gyromagnetic ratio, e the electron charge, M S the free layer's saturation magnetization, V the free layer's current-carrying volume, and the coefficients q i , i = 0, 1, 2, ... are assumed to describe the sample's characteristics. It is further assumed that I(t) creates a power perturbation of the form p = p 0 (1 + 2δp), where p 0 is the STO's free-running power. This approximation is valid as long as the modulation strength is small, which we consider a common scenario in STO modulation experiments (for instance, µ < 0.1 in both Ref. 29 and 30). Separating Eq. (2) into real and imaginary parts, and expanding in power to first order, we obtain a set of differential equations with variable (time-dependent) coefficients for the STO's power and phase: where we define the constants Here, is the restoration rate, α is the Gilbert damping parameter, ξ = I dc /I th is the supercriticality parameter, and I th is the threshold current for spin-torque-driven oscillations. The restoration rate is assumed to be constant with respect to the free-running oscillations. Qualitatively, it gives a measure of how fast the STO reacts to a perturbation, and so plays a fundamental role determining the synchronization speed for these devices. 13,14,25 We propose a Fourier series as a solution of Eq. (5)a, based on the periodicity of the time-dependent term: Introducing this solution into Eq. (5)a gives an infinite system of equations for the Fourier coefficients. It is possible to further reduce the problem to the determination of A n , leading to a system of linear equations which can be solved numerically (see Appendix A). The solution, expressed as a single sinusoid, is The notation of Eq. (8) explicitly shows a harmonicdependent phase ψ n = arctan(A n /B n ). This phase plays a fundamental role in the form of the PSD, as discussed below. The first indication of NFAM is visible in the first right-hand term of Eq. (8b), where the oscillation frequency is shifted by 2πf s = 2νΓ p A 0 . From the exact solution given in Eq. (A3), we obtain Consequently, the shift is directly proportional to the modulation strength and to the STO nonlinearities, via the constant C 2 .
In order to identify the next feature of NFAM-namely the asymmetry of the sideband-the PSD must be calculated. By expressing the sinusoidal functions as exponentials and performing Taylor expansion (Appendix B), the PSD can be expressed as a series of convolutions where we define the complex variable X n = B n + iA n (hence ψ n = arg(X n )), the notation(·) represents the complex conjugate, and the shifted carrier frequency is written as The product symbol represents in this case a series of convolutions. The real harmonic-dependent modulation index is defined as Comparison with the FM modulation index, h = ∆f /f m , suggests that the peak frequency deviation for STOs can be defined as 2π∆f n = 2νΓ p |X n | for given harmonic frequency and modulation conditions. Due to the form of |X n |, it is expected to be linearly dependent on the modulation strength µ and proportional to [(nω m ) 2 + (2Γ p ) 2 ] −1/2 . The latter dependence is of particular interest in terms of the modulation bandwidth (MBW) discussed in Section III.
Several features of the spectrum can be readily identified from the right-hand side of Eq. (10). First we can see the shifted carrier frequency which, by virtue of the properties of convolution, shifts the base-band spectrum in frequency. Second, a nonlinear amplitudemodulation (NAM) spectrum arises from the power fluctuations. Third, a series of FM spectra whose harmonics expand with n can be seen to arise from the combined contributions of the phase modulation and the nonlinearly enhanced power fluctuations. In the case of a weakly nonlinear oscillator (ν ≤ 1), the FM series of convolutions can be neglected, and the spectrum reduces to a NAM spectrum. For a particular applied field angle, the condition ν = 0 is satisfied, and the spectrum reduces to pure AM. 46 On the other hand, if the time-dependent term is neglected (C 2 ≈ 0), only the first harmonic coefficients of Eq. (A3) will be non-zero, and the PSD will take the form of pure FM. 41 The convoluting terms in Eq. (10) define the power of each harmonic. Although it is tedious to obtain a meaningful analytical expression, it can be inferred that the power of the harmonics is sensitive to the phase ψ n . We can illustrate this by approximating each convoluting term to the second harmonic, and solving for the first upper and lower sideband power. By expanding the absolute values to first order, we obtain the power difference of the sideband ∆ = 2J 1 (β 1 )(1 + A 0 )p 0 , which is in general non-vanishing. This crude approximation does not represent the asymmetry of Eq. (10), which is affected up to the fifth harmonic, as shown in Appendix A.
Summarizing this section, we have obtained a fairly complex spectrum by solving the set of Eq. (5). This complexity, along with the interdependence of several key parameters such as the modulation strength and frequency, limits the analytical insight that can be gained. However, the frequency shift of the carrier can be explicitly obtained by considering a Fourier-series approach, and the sideband asymmetry arises as a consequence of the solution's harmonic-dependent phases. The solution of the spectrum leads to the definition of the harmonic- dependent modulation index, from which the STO peak frequency deviation can be defined.

III. NUMERICAL RESULTS
In this section, we consider a first order expansion of the STOs parameters in order to evaluate the PSD of Eq. (10). The coefficients are calculated up to the tenth harmonic by solving the linear system of equations Eq. (A3) in matrix form. We consider the nanocontact geometry of Fig. 1 with the parameters µ 0 H app = 1 T, µ 0 M S = 0.8 T, γ = 28 GHz/T, α = 0.01, and ν = 100. Evaluating the STO parameters Γ + (p), Γ − (p), and ω(p) according to auto-oscillator theory, 41 we obtain the frequency vs supercriticality dispersion shown in Fig. 2a. As expected, the frequency is equal to the FMR frequency at the threshold of the oscillation (ξ = 1).
Three operation points (OP), OP1, OP2, and OP3, were selected in Fig. 2a, such that their restoration rates were Γ p /2π = 11.2, 44.8, and 156.8 MHz, respectively. Each OP was modulated, and their PSDs, calculated as the FFT of Eq. (B1), are shown in Fig. 2(b-d) as a function of β 1 and a fixed modulation frequency, f m = 100 MHz. The first-harmonic modulation index is used as a reference, since it is assumed to be the main contribution to the sideband power. In order to evaluate Eq. (A3), the modulation strength is calculated from β 1 , returning different ranges for each OP. The frequency axis is normalized to the modulation frequency and centered on the carrier frequency. On this scale, all operation points share the same number of visible sidebands, as well as the position of the carrier. It is observed that the frequency shift is more pronounced near the threshold (OP1). This is expected from Eq. (9), shown as the white lines in Fig. 2(b-d), since the modulation strength is higher near both the threshold and the curvature of Fig 2(a). Such curvature-dependence has been observed in experiments 30 away from the oscillation threshold, presumably due to higher order nonlinearities. These effects can be included in the present framework through Eq. (6b), as the derivative of the negative damping may be substantially different. A numerical fit of experimental data is not within the scope of the present paper, but we argue that such a fit could be achieved by measuring the auto-oscillator parameters, as has recently been shown in Ref. 7.
The power of the normalized carrier and of the first upper and lower sidebands is shown in Fig. 3(a-c) for each OP and with modulation frequency f m = 100 MHz. The combined action of the frequency shift and the harmonicdependent phase introduced by the Fourier series results in power asymmetry between the sidebands. This asym-metry is notably higher for OP1 in correlation with its similarly enhanced frequency shift. In the color plots of Fig. 3(d-f), the power difference of the first upper and lower sidebands ∆ is shown as a function of β 1 and f m . The white horizontal line denotes the slice shown in the panels (a-c), respectively. Being a common feature, the sideband asymmetry increases together with both the modulation strength and the frequency. These dependencies can be qualitatively understood from Eq. (10) and Eq. (A3). A stronger µ increases the power of the higher harmonic coefficients, so that they become non-negligible in the PSD. Thus, an increase in the asymmetry is expected. On the other hand, the dependence on the modulation frequency arises from the relation between A n and B n . Roughly, one can approximate A 1 = ω m /2Γ p B 1 , so that the complex variables in the PSD can be expressed as X 1 = B 1 (1 + iω m /2Γ p ). From here, it is clear that as ω m increases past the parameter 2Γ p , the imaginary term becomes dominant or, in other words, the harmonic phase increases towards ±π/2. Consequently, the power contribution of each harmonic to the sideband becomes heavily weighted by the phase, leading to enhanced asymmetry.
Finally, we discuss the features of the peak frequency deviation ∆f n in relation to the modulation bandwidth. The frequency dependence of ∆f for each OP and fixed µ = 0.05 is shown in Fig. 4, both by evaluating the firstharmonic modulation index, Eq. (11) (solid lines), and by calculating the maximum of the instantaneous frequency from Eq. (8b) (circles). The magnitude of the deviation is shown relative to the OP oscillator frequency, while the logarithmic scale is used to enhance the frequency dependence. Both approaches agree very well, so that the first harmonic modulation index suffices to perform the following calculations approximately. A practical consequence of this feature is that choosing a maximum modulation strength does not guarantee the frequency excursion estimated from the STO's frequency vs current characteristics. This is closely related to the modulation bandwidth of oscillators.
The modulation bandwidth (MBW) gives a measure of the frequency range in which an oscillator has optimal modulation properties. A common criterion is the 3 dB power attenuation that is likewise used to characterize filters. Since the power in a FM scheme depends on the Bessel functions, the methods used to measure the MBW rely on indirectly estimating the degradation of the linearity of the modulation index. In the present framework, the modulation index is defined from the calculation of the power spectrum (Appendix B), and one can directly estimate the MBW from Eq. (11). This can be accomplished by considering a vanishing modulation frequency and an arbitrary µ, so that the first-harmonic modulation index is β i . The task is then to keep the ratio µ/ω m constant while looking for β f = β i / √ 2. Performing this calculation gives the result that the MBW for STOs is 2Γ p . Since the lack of linearity in our framework is given by 2π∆f , one can now understand Fig. 4 as a STO transfer function under modulation having a char-

IV. CONCLUSIONS
The NFAM spectrum is obtained from the autooscillator general model by considering the power perturbation created by the modulating signal. These power fluctuations are enhanced by the STO's nonlinearity, resulting in a spectrum consisting of a series of convolutions between NAM and FM spectra. Its implicit dependence on the total damping parameter and nonlinearities suggests that the NFAM characteristics will exhibit sample-to-sample variation. In fact, such characteristics have been observed in experiments, 30 and their dependence on the modulation parameters qualitatively agrees with the results shown here. In contrast to the originally proposed model, 40 we found a frequency dependence of the modulation index which leads to the estimation of the modulation bandwidth for STOs. Moreover, the STO under modulation has a low-pass filter behavior with a cutoff frequency given by 2Γ p , which usually depends on the operation point of a specific sample. Although the modulation bandwidth for STOs has yet to be measured experimentally, we expect a correlated degradation of STO modulation characteristics in frequency-dependent studies. On the other hand, the Fourier series solution proposed in this paper can, in principle, reduce to a linear set of coupled equations any STO geometry described by the auto-oscillator general equation Eq. (2) perturbed by slow time-varying signal. Support from the Swedish Research Council (VR) is gratefully acknowledged. JohanÅkerman is a Royal Swedish Academy of Sciences Research Fellow supported by a grant from the Knut and Alice Wallenberg Founda-tion.

Appendix A
The Fourier series solution proposed in Eq. (7) is introduced into Eq. (5a). The variable term can be expanded using trigonometric identities: This series can be then rearranged by expanding the summation and identifying equal harmonics. Changing the index accordingly, we obtain This expression treats the harmonics separately, so that Eq. (5a) can be solved by collecting harmonic terms. The resulting system of equations for the coefficients is There are two possible solutions to the set of Eq. (A3). One can, for instance, express a system of equations depending only on the As, so that A n ∝ A n+2 , A n+1 , A n−1 , A n−2 . This new system can easily be solved numerically using matrix algebra. The error will depend on the size of the matrix or, in other words, on the maximum harmonic where the series is truncated. The total error is shown in Fig. 5(a) for two frequencies of different orders of magnitude. The case for N = 20 harmonics is used as a reference. It is observed that in both cases, N ≥ 5 returns a total error smaller than 10 −5 %. Another solution is grounded on the fact that the harmonic coefficients decay with n, so that we can assume in general that A n ∝ A n−1 , A n−2 . This approach can be implemented numerically as a recurrent series, and this is shown for different values of N and f m in Fig. 5(b-c). Here we observe that N = 5 also converges to a minimum error in this approximation. However, as β 1 is swept out, it is clear the the error increases towards the minimum of the first sideband (indicated by a dashed line). It is noteworthy that the maximum error in Fig. 5(c) is slightly shifted, which corresponds to the impact of higher harmonic terms.

Appendix B
The power spectral density (PSD) of the proposed solution is obtained from Eq. (8), which defines both timedependent power and phase variations. The PSD is defined as P SD = p 0 * f (1 + δp) * f (cos(φ)) , wheref is the Fourier transform. In the following we perform each Fourier transform separately, and express the result as a convolution. The first term on the right hand side has the form of a NAM. The phase introduced by the sine function leads us to define the complex variable X n = B n + iA n , so that f (1 + δp) = (1 + A 0 )δ(0) (B2) + ∞ n=1X n 2 δ(nω m ) + X n 2 δ(−nω m ), where the notation δ(x 0 ) = δ(x − x 0 ) is used and δ is the Kronecker delta. The second term on the right hand side of Eq. (B1) is calculated by using Euler's formulae, and subsequently expanding in Taylor series. We obtain where a n = 2νΓ p A n /(nω m ), b n = 2νΓ p B n /(nω m ) and c.c is the complex conjugate. The two summations can be multiplied by expanding and rearranging the terms. Defining 2x n = b n + ia n , it becomes possible to rewrite the product terms as (−1) k x nxn k!(k + j)! (x n e) j + (−x n e −1 ) j .
From this expression, we can identify the summations on k as Bessel functions of order j and argument β n = 2 √ x nxn (upon renormalization of the coefficients by β n /2). The product takes the form ∞ n=1 (·) = J 0 (β n ) The Fourier transform of Eq. (B5) is easily performed by expressing the result as a product of convolutions. By further evaluating x n in terms of the complex variable X n , and convoluting each harmonic contribution, we obtain Eq. (10).