Injection locking at fractional frequencies of magnetic tunnel junction (MTJ)-based read sensors' ferromagnetic resonance modes

TL;DR

This study investigates the nonlinear magnetization dynamics of MTJ-based read sensors, revealing fractional frequency resonances and their dependence on magnetic and spin-torque parameters through NVNA measurements and micromagnetic modeling.

Contribution

The paper demonstrates the presence of fractional frequency resonances in MTJ sensors and links these to nonlinear spin-torque dynamics, enhancing understanding of sensor behavior.

Findings

Distinct resonances at fractional frequencies of FMR mode identified

DC response explained by low-order nonlinearity and magnetodipolar feedback

Harmonic response is strongly influenced by mutual spin-torque effects

Abstract

Being nonlinear dynamic systems, magnetic read sensors should respond to an excitation signal of a frequency considerably different from their natural ferromagnetic resonance (FMR) frequencies. Because of the magnetization dynamics' inherent nonlinear nature, the sensors' response should be measured at the DC, excitation frequency, and its multiples (harmonics). In this paper, we present results of such measurements, accomplished using a one-port nonlinear vector network analyzer (NVNA), which show distinct resonances at fractional frequencies of the free layer (FL) FMR mode. Identification of these resonances, resulting from the nonlinear nature of the spin-torque (ST)-induced magnetization dynamics, was performed using micromagnetic modeling. In particular, we show that the measured DC response at the above-mentioned fractional frequencies can be explained by a low-order nonlinearity and strong magnetodipolar feedback between magnetic layers adjacent to an MgO barrier. Additionally, we determined that the simulated harmonic response is strongly enhanced by the mutual ST effect between these layers. Finally, we demonstrate that the read sensors' nonlinear magnetization dynamics and, by extension, their harmonic response are highly sensitive to various magnetic and ST parameters. Thus, this study shows that using NVNA measurements in conjunction with micromagnetic modeling can clarify the uncertainty in the definition of these parameters.

Equations

$$I=I_{\text{DC}}+I_{\text{AC}}\sin (\omega t).$$

$$V(t)=I(t)R(t)=\left[I_{\text{DC}}+I_{\text{AC}}\sin (\omega t)\right]\times \left[R_0+R_{\text{TMR}}(t)\right].$$

$$R_{\text{TMR}}(t)=R_{\text{TMR}}^{\text{av}}+\sum _{n=1}^{\infty }{a}_n\cos (n\omega t)+\sum _{n=1}^{\infty }{b}_n\sin (n\omega t),$$

$$V_{\text{DC}}=⟨V(t)⟩=\underset{\text{“constant”}}{\underbrace{I_{\text{DC}}{R}_0+I_{\text{DC}}{R}_{\text{TMR}}^{\text{av}}}}+\underset{\text{“oscillating”}}{\underbrace{I_{\text{AC}}\frac{b_1}{2}}},$$

$$\text{derivative}\mathbf{M}t=-\gamma \left[\mathbf{M}\text{crossproduct}\left({\mathbf{H}}_{\det }+{\mathbf{H}}_{\mathrm{th}}\right)\right]-\gamma \frac{\lambda }{M_s}\left[\mathbf{M}\text{crossproduct}\left[\mathbf{M}\text{crossproduct}\left({\mathbf{H}}_{\det }+{\mathbf{H}}_{\mathrm{th}}\right)\right]\right],$$

$${\mathbf{H}}_{\mathrm{ST}}=f_J(\theta )\left[\mathbf{M}\text{crossproduct}\mathbf{p}\right].$$

$$f_J(\theta )=a_J\frac{2{\Lambda }^2}{({\Lambda }^2+1)+({\Lambda }^2-1)\cos \theta },$$

$$a_J=\frac{\hslash }{2\text{absolutevalue}e}\frac{jP}{M_s^{2}d}.$$

$$\begin{array}{rl}⟨{\mathbf{H}}_{\xi ,i}^{\mathrm{th}}⟩& =0,\\ ⟨{\mathbf{H}}_{\xi ,i}^{\mathrm{th}}(0){\mathbf{H}}_{\psi ,j}^{\mathrm{th}}(t)⟩& =2D\delta (t){\delta }_{ij}{\delta }_{\xi \psi },\end{array}$$

$$D=\frac{\lambda }{1+{\lambda }^2}\frac{kT}{\gamma \mu },$$

$$R=\frac{1}{G}=\frac{R_{\perp }}{1+\frac{{\Delta }_{\text{TMR}}}{2}\cos \theta },$$

Full text

I Introduction

A magnetic tunnel junction (MTJ) is the basis of modern read sensors commonly used in hard disk drives. As shown in Fig. 1, a typical MTJ represents a nanoscale multilayered structure composed of the free and pinned layers (FL and PL, respectively). The PL is composed of two antiferromagnetically coupled layers (PL1 and PL2) to reduce the stray field from the PL onto the FL. The antiferromagnetic (AFM) layer fixes the magnetization orientation of the PL1 via the direct exchange coupling between the AFM and PL1. The side bias field keeps the FL magnetization orthogonal to the PL. The stack’s shape is tapered to increase the FL’s magnetic stability. The magnetization direction of the PL is assumed to be perfectly fixed while that of the FL rotates in response to the field from the magnetic media. The change in the relative orientation between the FL and PL magnetizations translates into variations in the sensor’s tunnel magnetoresistance (TMR), which indicate the difference between the “up” and “down” bits representing the recorded information [1, 2].

MTJs are nonlinear dynamic systems. At room temperature, thermal fluctuations result in random magnetization dynamics of both the FL and PL, which can be quantified by their ferromagnetic resonance (FMR) modes [3]. Another source of magnetization dynamics in MTJs is the spin-torque (ST) effect, which describes a direct transfer of angular momentum from the spin-polarized electrons to the local magnetization [4]. Even in the macrospin approximation, magnetization dynamics are inherently nonlinear because oscillations of magnetization components arise from the precession of the magnetization vector with a constant magnitude, leading to the following nonlinear relation between these components: $m_{x}^{2}+m_{y}^{2}+m_{z}^{2}=1$ [5]. ST-induced magnetization dynamics result in large-angle magnetization precession where the nonlinear nature of magnetization dynamics becomes especially pronounced.

The fundamental theory of nonlinear oscillations states that their nonlinearity results in new resonances such that oscillations of frequency close to $f_{0}$ (natural frequency) can be excited by an external force with a frequency considerably different from $f_{0}$. Namely, a resonant condition might occur at every excitation frequency $pf_{0}/q$, where $p$ and $q$ are positive integers. In practice, however, $p$ and $q$ should be small because the resonance strength rapidly decreases with increasing order of nonlinearity [6].

Specific cases of forced oscillations in nonlinear systems excited at $nf_{0}$ and $f_{0}/n$, where $n$ is a positive integer, are called harmonic (or super-harmonic) and sub-harmonic injection locking, respectively. Injection locking is observed in numerous types of physical systems. Most often, however, this term is associated with electronic oscillators and laser resonators. In optics, injection locking has been used to improve the frequency stability of lasers and reduce the frequency noise of laser diodes. In electronic systems, injection locking has been used to increase the pull-in (or “capture”) range and reduce the output phase jitter in phased-locked loops [7].

In ST-driven systems, e.g., spin-torque nano-oscillators (STNOs), super- and sub-harmonic injection locking can be used to generate microwave and millimeter wave signals. In Ref. [8], Keatley et al. excited DC-biased spin-torque vortex oscillators (STVOs) with an AC signal at multiples (harmonics) of their fundamental frequency. The response was measured at a fractional frequency of the excitation signal corresponding to the fundamental frequency. In Ref. [9], Lebrun et al. presented an experimental study of both super- and sub-harmonic injection locking to an AC excitation signal in DC-biased STVOs. It resulted in pure phase locking with no phase slips and an output power of $>\!\!1~{}\mu\textrm{W}$ observed at room temperature and zero magnetic field. In Ref. [10], Carpenteri et al. performed a numerical study of both super- and sub-harmonic injection locking in STNOs based on hybrid spin-valves composed of two FLs and orthogonal PLs. In Ref. [11], Quinsat et al. reported the results of a comprehensive synthesis of experimental and numerical studies of the super-harmonic injection locking in STNOs. This work was inspired by the analytical theory of fractional synchronization presented in Ref. [12], which, in turn, concluded that for large AC power levels the coupling between the oscillating magnetic layer and the source becomes nonlinear, resulting in a fractional synchronization regime [13].

In this work, we study the sub-harmonic injection locking in MgO-based MTJs. Combining the nonlinear vector network analyzer (NVNA) measurements and micromagnetic simulations, we show that this locking is determined not only by the intrinsic nonlinearity of the magnetization dynamics, but also by the magnetodipolar feedback between the FL and PL. Furthermore, we demonstrate that the orthogonality of magnetization directions of the FL and PL1 also plays an important role, facilitating the read sensors’ response at 1/2 the FL FMR frequency $f_{\textrm{FL}}$. We observe distinct peaks in the DC response at 1/2, 1/3, 1/4, and 1/6 of $f_{\textrm{FL}}$. Performing corresponding simulations, we show that these resonances are due to the combination of the low-order nonlinearity, magnetodipolar feedback, and the mutual ST effect between the FL and PL1. Importantly, strong magnetodipolar feedback permits sub-harmonic injection locking within a wide range of integer fractions of $f_{\textrm{FL}}$.

The presented work suggests that NVNA measurements in conjunction with micromagnetic modeling can greatly assist in clarifying the uncertainty in defining the system’s magnetic and ST parameters. Furthermore, such a synthesis of nonlinear measurements and modeling identified the need to consider mutual ST between the FL and PL in micromagnetic modeling.

II Nonlinear characterization

We use the advances in microwave measurement techniques to determine the read sensors’ DC and harmonic responses, which comprise nonlinear characterization. Due to the broadband nature of and variability in the GHz-range FMR in read sensors, the nonlinear characterization was performed over a wide frequency range (the excitation frequency ranging from 1 to 15 GHz). Additionally, due to the presence of magnetic shields and bias magnets (Fig. 1), our samples are unsuitable for characterization with an applied magnetic field. The shields would distort the external field such that one would be unable to control the relative orientation of the FL and PLs to the degree required for such studies. Therefore, all measurements presented in this work were performed without an applied magnetic field, and the excitation was a combination of AC power and DC bias current.

II.1 DC response

We start with the definition of the voltage induced in the read sensor, which is subjected to both DC and AC currents:

[TABLE]

In the presence of the TMR effect, the system’s resistance has two constituents: magnetization-independent ohmic resistance $R_{0}$ and magnetization-dependent, TMR-caused resistance $R_{\textrm{TMR}}(t)$. Hence, the resultant time-dependent voltage response is

[TABLE]

If $I_{\textrm{DC}}$ is less than the critical current above which the ST-induced steady-state oscillations arise, the only excitation present in the system is the periodic AC current. Thus, the time-dependent magnetization resulting in the time-dependent TMR term is a periodic function with a period $T=2\pi/\omega$. Due to various sources of nonlinearities in the system, the TMR-caused resistance can be non-harmonic and subsequently be expressed as a Fourier series:

[TABLE]

where the Fourier coefficients $a_{n}$ and $b_{n}$ have the dimensionality of resistance. The averaged $R_{\textrm{TMR}}^{\textrm{av}}$ is determined by the nonlinearity of the magnetization oscillations and is highly sensitive to ST asymmetry. This issue is discussed in more detail in Section IV. Substituting Eq. (3) into Eq. (2) and then time-averaging the result produces the following DC voltage response:

[TABLE]

where all other higher-order terms disappear after time-averaging due to the orthogonality properties of harmonic functions.

In Eq. (4), the “constant” term comes from the magnetization-independent (ohmic resistance) and time-averaged magnetization-dependent, TMR-caused resistances. The “oscillating” term involves mixing the excitation signal with the TMR oscillations only at the same frequency [14]. From Eq. (4), it can also be seen that the contribution from the “oscillating” term depends on the phase difference $\Delta\phi$ between these two signals: its magnitude is at its maximum at $\Delta\phi=0^{\circ}$ or $180^{\circ}$ (in- or out-of-phase) and its minimum at $\Delta\phi=90^{\circ}$ or $270^{\circ}$ (quadrature).

In this study, the typical MTJ’s dimensions are $30\times 35~{}\textrm{nm}^{2}$ in the $x$-$y$ plane and $\approx~{}\!\!25~{}\textrm{nm}$ in the $z$-direction. Its resistance-area product and TMR ratio are $0.4~{}\Omega\mu\textrm{m}^{2}$ and $90\%$, respectively. Experimentally, the sensor’s DC response can be extracted using the steady-state DC readout measurement technique [15, 14]. In the simplest measurement configuration (Fig. 2), the MTJ-based read sensor is excited with an AC signal for different DC bias currents, while the DC voltage across the sensor is measured with a sourcemeter.

In this study, six Western Digital production read sensors were characterized. The representative sample’s DC readout measurements shown in Fig. 2 were obtained using the measurement technique introduced in Fig. 2 (utilizing a Keysight PNA-X N5247A as the AC source). Our main observation is the presence of peaks emerging in the DC readout measurements with non-zero DC bias current at frequencies that are the integer fractions ($1/2$, $1/3$, $1/4$, and $1/6$) of the FL’s resonant precession frequency $f_{\textrm{FL}}$. A similar observation was made in a micromagnetic study of the current-perpendicular-to-plane (CPP) spin valve heads under ST excitation [16]. A subsequent study of ST-induced magnetization dynamics in thin magnetic nanoelements demonstrated that nanoelements of a certain size can exhibit splitting of their magnetization precession trajectory into limiting sub-cycles [17]. In the spectrum, these sub-cycles corresponded to peaks at frequencies lower than that of the complete motion cycle. In Ref. [17], this regime precedes the state of chaos characterized by strongly inhomogeneous large-angle magnetization precession, which produces chaotic trajectories and almost continuous spectrum with no distinct FMR peaks. In both micromagnetic studies, however, the state near bifurcation (characterized by the presence of peaks at the fractional frequencies of $f_{\textrm{FL}}$) was observed in a narrow range of corresponding parameters: applied DC bias currents (2.8 to 3.6 mA) in Ref. [16] or sizes of nanoelements (40 to 52 nm) in Ref. [17]. In contrast, we observe peaks at the FL’s fractional frequencies in a wide range of applied AC and DC signals, which requires an alternative explanation.

II.2 Harmonic response

The steady-state DC readout measurement technique allows analysis only of the DC spectral component. Since nonlinear effects are involved, a thorough nonlinear harmonic analysis is highly desired. Such an analysis can be accomplished by using an NVNA. Since the characterized magnetic read sensors are one-port devices, this analysis is limited to the one-port scenario.

In contrast to the linear VNA, which measures the magnitude ratio and phase difference between the incident and reflected waves ($A_{1}$ and $B_{1}$, respectively) only at the excitation frequency, the NVNA measures the actual $A_{1}$ and $B_{1}$ waves’ magnitudes and phases at the excitation frequency as well as harmonic components to which energy may be transferred due to the device’s nonlinear characteristics (Fig. 4) [18, 19]. In this measurement, the excitation frequency $f_{i}$ of the incident wave $A_{1}$ was swept from 1 to 15 GHz. At each frequency, the fundamental of the reflected wave $B_{1}$ was measured at $f_{i}$. Naturally, the second harmonic of the reflected wave $B_{1}$ was measured at $2f_{i}$ (2 to 30 GHz) and the third harmonic was measured at $3f_{i}$ (3 to 45 GHz). We emphasize that the resultant measurements represent the frequency response of $B_{1}$, not the power spectrum. A power spectrum is the distribution with frequency of the power content of the signal [20], whereas the frequency response [21] curves in Fig. 5 represent the magnitude of the corresponding harmonic component of $B_{1}$ as a function of the excitation frequency and are linearly proportional to voltage.

In Fig. 5, the magnitudes of the fundamental, second, and third harmonics of $B_{1}$ are counterposed to the DC readout curve. In Fig. 5(a), the fundamental of $B_{1}$ exhibits distinct peaks at frequencies corresponding to $f_{\textrm{FL}}$, $f_{\textrm{FL}}/2$, $f_{\textrm{FL}}/4$, and (less pronounced) $f_{\textrm{FL}}/6$. These peaks have their counterparts on the DC response curve. In addition, strong peaks are observed in-between $f_{\textrm{FL}}$ and $f_{\textrm{FL}}/2$. They might correspond to resonances at frequencies $pf_{\textrm{FL}}/q$ (where $p$ and $q$ are positive integers) and have no analogues in the DC response. Understanding the nature of these additional peaks requires further investigation. In Figs. 5(b) and 5(c), the second and third harmonics of $B_{1}$ also reveal distinct peaks at $f_{\textrm{FL}}/2$, $f_{\textrm{FL}}/4$, and (less pronounced) $f_{\textrm{FL}}/6$. Interestingly, the FL mode and the peak located at $f_{\textrm{FL}}/3$ do not produce the second and third harmonics of $B_{1}$.

The NVNA measurements imply that the read sensor’s nonlinear characteristics give rise to nonlinear oscillations under the AC excitation signal with the frequency considerably lower than, but still an integer ratio of, the FL mode. Naturally, these dynamics also produce a measurable DC contribution corresponding to this excitation frequency.

III Micromagnetic modeling

We employ micromagnetic modeling to understand the underlying physics that enables the DC and harmonic responses at frequencies corresponding to the FL FMR mode’s fractional ratios. Due to the model’s complexity and uncertainty in the definition of some of its magnetic and ST parameters, we find it vital to discuss in detail the modeling methodology and modeling considerations.

III.1 Modeling methodology

All numerical simulations presented in this work were performed using the micromagnetic software MicroMagus, which solves the modified Landau-Lifshitz-Gilbert (LLG) equation for the magnetization M [22]

[TABLE]

using one of the optimized Runge-Kutta or Bulirsch-Stoer algorithms with an adaptive step-size control for both $T=0$ and $T>0$. In Eq. (5), the precession constant $\gamma$ is defined via the absolute value of the gyromagnetic ratio $\gamma_{0}$ as $\gamma=\gamma_{0}/(1+\lambda^{2})$. The damping constant $\lambda$ is equal to the corresponding damping $\alpha$ in the LLG form where the magnetization derivative is present on both sides of the equation.

The deterministic field $\bf{H}_{\rm det}$ contains four standard contributions: external, anisotropy, exchange, and magnetodipolar interaction fields. In our case, the ST effect is taken into account via an additional effective field term, which in the standard Slonczewski formalism  [23, 24] has the form

[TABLE]

Here, the dimensionless ST amplitude $f_{J}$ depends on the angle $\theta$ between the magnetization M and the spin-polarization direction p as follows [24, 25]:

[TABLE]

where the factor $a_{J}$ is given by

[TABLE]

In Eqs. (7) and (8), $e$ is the electron charge, $j$ is the electric current density, and $d$ is the thickness of the magnetic layer subjected to ST. The asymmetry parameter $\Lambda$ strongly depends on the sample configuration and various transport coefficients [24, 25]. When $\Lambda=1$, the ST effect is assumed to be symmetric. $P$ is the degree of spin polarization of the electrical current 333To describe the ST’s angular dependence in MgO-based junctions, we used the common formulations for metallic junctions. Recent work (Ref. [40]) shows that the effective difference in angular dependence for these two scenarios is insignificant..

In the LLG equation (5), thermal effects are accounted for by the thermal field term $\bf{H}_{\rm th}$ describing random fluctuations induced by the interaction of the ferromagnet with the thermal bath. Components of this fluctuation field have the following statistical properties:

[TABLE]

meaning that these fluctuations are assumed to be uncorrelated in space and time ($i$, $j$ are the discretization cell indices; $\xi,{}\psi=x,{}y,{}z$). The noise power $D$ is proportional to the system temperature $T$:

[TABLE]

Here, $\mu$ is the magnetic moment magnitude of a discretization cell. Unless stated otherwise, all micromagnetic simulations in subsequent sections were performed at $T=0$ K.

The TMR response introduced in Section II.1 is the quantitative description of the TMR’s dependence on the angle $\theta$ defined between the magnetizations of adjacent layers [27]:

[TABLE]

where $R_{\perp}$ is the sensor resistance at the orthogonal state and $\Delta_{\textrm{TMR}}$ is the TMR ratio.

Equation (11) is different from the generally used formulations describing a linear variation of the MTJ’s resistance with $\cos\theta$, which appears to be valid only for small values of TMR [27].

III.2 Modeling considerations

In this study, the MgO-based MTJ devices’ PL1, PL2, and FL are made of a similar CoFe/CoFeB alloy. The layers’ saturation magnetization and exchange stiffness constant were assumed to be $M_{s}=1\times 10^{3}~{}\textrm{G}$ and $A=1\times 10^{-6}~{}\textrm{erg/cm}$, respectively. The intrinsic Gilbert dampings of the FL and PL2 were set to $\lambda_{\rm FL}=\lambda_{\rm PL2}=0.01$. The IrMn-exchange-pinned damping of the PL1 was initially assumed to be an order of magnitude larger, $\lambda_{\rm PL2}=0.1$ [28]. The interlayer coupling strengths were set to values similar to those used in Ref. [29]: $J_{1}=0.04~{}\textrm{erg/cm}^{2}$ between the FL and PL2 (“orange-peel” coupling) and $J_{2}=-1.6~{}\textrm{erg/cm}^{2}$ between the PL1 and PL2 (strong antiferromagnetic coupling via the Ru interlayer).

For each experimental frequency point at a given AC power level and DC bias, we had to perform an independent time-domain micromagnetic simulation. For linear systems that respond only at the excitation frequency, fast broadband excitation can be accomplished via sinc pulses or multifrequency signals (e.g., Schroeder-phased harmonic signals) that are optimal for uniform excitation of all system modes [30, 31]. This is done so that, e.g., the AC-susceptibility in the whole frequency range of interest can be sampled in a single simulation run. For nonlinear systems, such an approach is not feasible because 1) nonlinear systems respond not only at the excitation frequency but also at its multiples (i.e., harmonics) which strongly interact with each other and 2) by exciting the system with a strong ST pulse, reaching a dynamic equilibrium state requires significant time, which may be much longer than the duration of the optimal sinc pulse.

To enable faster frequency sweeps around the ranges of interest (i.e., FL FMR frequency and its fractional ratios), the following simplifications to the original read sensor design were adopted: 1) The AFM layer was excluded from simulations, but its effect on the PL1 was taken into account by the corresponding exchange bias field of 1000 Oe [29]. 2) The left and right side bias magnets were also excluded from simulations. Instead, the side bias demagnetizing field was calculated in the quasi-static solver and then included as an external field in all subsequent dynamic simulations. In this model, each magnetic layer is discretized in-plane into $N_{x}\crossproduct N_{y}=32\times 20$ cells. No discretization was performed in the out-of-plane direction. We verified that introducing such a discretization even for the thickest layer (i.e., FL) did not lead to any significant changes in final results. With these simplifications, MicroMagus simulations required $\approx\!\!1$ hour to collect 60 nsec of magnetization dynamics at $T=0$ K and $\approx\!\!15$ hours to simulate 200 nsec at $T=300$ K.

To the best of the authors’ knowledge, there have not been any comprehensive micromagnetic studies on the nonlinear magnetization dynamics used to predict the MTJ’s harmonic response. In general, these dynamics and, by extension, its harmonic response should be sensitive to the various magnetic and ST parameters described above. To test this, we excited the read sensor model via a combination of AC and DC signals (as in the experiment in Fig. 5) at 1/4 the frequency of the FL FMR mode.

Effect of the side bias and magnetic shielding. In Fig. 1, the read sensor model is shown without the magnetic shields. In a real system, the top and bottom magnetic shields isolate the sensor from adjacent bits and large writer fields. Material parameters of both shields are typical of permalloy. The shields are sufficiently larger than the read sensor. Being unable to micromagnetically model the full-size shields, we reduced their dimensions to $L_{x}\crossproduct L_{y}\crossproduct L_{z}=300\times 36\times 100~{}\textrm{nm}^{3}$ and applied the periodic boundary conditions along the $x$-axis to avoid the influence of artificial magnetic “surface charges” from the vertical (in the $y$-$z$ plane) shield surfaces. By doing this, we estimated the quasi-quantitative effect of these shields on the side bias field (Fig. 6).

Figure 6 shows that in the presence of magnetic shields, the simulated side bias field, both onto the FL and especially onto the pinned layers PL2/PL1, sufficiently decreased. This means that the total side bias field (i.e., taking into account the presence of magnetic shields) strongly depends on the shield configuration. Additionally, variations in the spacing between the side bias and the MTJ stack within the fabrication tolerances may affect the side bias field strength. The side bias magnetization is also known to be accurate within approximately $\pm 10\%$. Hence, without having more accurate information concerning the geometry and magnetic parameters of the shields and side bias, we have to adjust the side bias field to match the most reliable experimental results.

The most direct experimental observation is the measured FL FMR frequency. This value can be used to adjust the total external field onto the FL. Based on the obtained value of $f_{\textrm{FL}}\approx 7~{}{\rm GHz}$, we have determined that the increase in the simulated side bias field onto the FL by a factor of $\approx\!\!1.7$ is necessary to reproduce this frequency. Moreover, we have found that the proposed adjustment of the side bias field onto the FL has shifted the system’s dynamic regime from quasi-chaotic towards pure phase-locking when excited at fractional frequencies. On one hand, the quasi-chaotic regime results in a stronger DC response. On the other hand, this regime was observed only within a narrow range of magnetic and ST parameters and thus is most probably absent in our system. Hence, the quasi-chaotic regime should be avoided, which also justifies the increase in the FL side bias field.

Adjustment of the side bias field onto the pinned layers PL2/PL1 is a more subtle issue. Here, our main criterion was that the side bias field onto the pinned layers should be such that it allows strongly nonlinear magnetization dynamics, leading to the generation of the higher-order harmonics and the sought after DC response (Fig. 7). Performing test simulations, we discovered that a factor of $\approx\!\!0.2$ reduction in the simulated side bias field onto the PL2/PL1 is necessary to obtain the measurable DC response and higher-order harmonics.

Effect of damping. Another ambiguity is the value of the effective damping $\lambda$ included in the LLG equation used for simulations. A lower damping constant corresponds to lower power absorption by the magnetic system, consequently leading to nonlinear oscillations with a larger amplitude and, by extension, a stronger harmonic response. While 0.01–0.02 is a typical damping constant for the FL and PL2 [32], the PL1 damping was initially set to a much higher value, $\lambda_{\rm PL1}=0.1$ [28]. In Ref. [28], Smith et al. suggested the following mechanisms that can explain a factor of 10 increase in the PL’ “standard” damping: 1) PL-FL spin-pumping and 2) strong interfacial exchange coupling at the IrMn/PL interface.

In test simulations, we have observed that the pinned layer PL1 damping strongly affects the read sensor’s dynamics (Fig. 8) even though this layer is assumed to be pinned. This effect is due to the large magnetodipolar interaction not only between the PL1 and PL2, but also between the PL1 and FL. The latter interaction is strong because the FL is relatively thick ($\approx\!\!7~{}{\rm nm}$). These couplings lead to a significant energy transfer to the PL1 and substantial energy dissipation resulting from the overdamped dynamics of the PL1. As a consequence, neither a noticeable DC response nor a significant harmonic response is observed at $\lambda_{\rm PL1}=0.1$ (Fig. 8(b)).

Based on these observations, we compared the system studied in Ref. [28] with our sensor composition. In contrast to our MTJ-based stack, Smith et al. studied current-perpendicular-to-plane giant-magnetoresistive (CPP-GMR) spin-valve stacks. Their samples, being just a pinned layer coupled to an AFM layer, did not have a PL2/PL1 with Ru in between. This can lead to significant differences between the effective dampings of the pinned layer adjacent to the AFM layer in these two systems. In Ref. [33], Mohammadi et. al. reported an inverse-thickness-squared dependence of damping for exchange-biased CoFe layers and increased damping (but still lower than 0.1 for a 3-nm-thick CoFe layer) via spin-pumping. Thus, based on the arguments presented above and partially on the results from Ref. [33], we have set the PL1 damping constant to the same value $\lambda_{\rm PL1}=0.01$ as for other layers. Still, further studies of this question are highly desired.

Mutual ST effect between the FL and PL2. Within a multilayered stack, the ST term (6) can be taken into account only on the FL if and only if the following statements are true: 1) The FL is significantly thinner than all other layers. Being a surface effect, the ST is more efficient for thinner layers. 2) The PL is usually pinned to the AFM layer by an exchange bias coupling, which is much stronger than the external field onto the FL.

In our case, both statements are false. In the read sensor design, the FL is much thicker than the PL1 (7 nm versus 3 nm). Furthermore, the external (side bias) field onto the FL is approximately the same as the exchange bias field onto the PL1, which is antiferromagnetically coupled to the PL2 (both fields are $\approx\!\!1000$ Oe). Thus, the FL-PL2 coupling must be considered in any adequate treatment of magnetization dynamics in the read head system. To account for this interaction, we have included the ST terms (6) on both PL2 and FL. The direction of the electron polarization used to compute the ST effect on the FL has been adjusted based on the magnetization direction of the PL2, and vice versa. The ST parameters $\Lambda$ and $P$ were set to be the same for both layers.

Test simulations have confirmed the importance of the mutual ST effect between the FL and PL2 (Fig. 9(a)). In the presence of the ST effect on the PL2, the peak at $f_{\textrm{FL}}/2$ is $\approx\!\!5$ times larger (red dots) than when it is absent (blue dots). The influence of the mutual ST effect is especially pronounced in the the DC response’s “constant” term (Fig. 16(a) in Appendix A).

Stability of the dynamic regime with respect to thermal fluctuations. Nonlinear dynamic systems are prone to chaotic behavior [34]. In our model, e.g., we achieved a quasi-chaotic regime when the side bias field onto the FL was too low. Hence, it was necessary to verify the stability of the phase-locking regime with respect to thermal fluctuations. We accomplished this by comparing the simulated DC responses obtained at zero temperature and $T=300~{}\textrm{K}$. Figure 9(b) shows that room-temperature fluctuations did not disturb the phase-locking regime and allowed the peak at $f_{\textrm{FL}}/2$ to remain intact.

IV Results and discussion

The micromagnetic modeling considerations proposed in Section III.2 helped to calibrate the model. The next step is to relate the simulated DC and harmonic responses to the experimental ones.

To aid the reader, in Fig. 10 we show the representative sample’s experimental DC response (extracted from Fig. 2) limited to the frequency range from $f_{\textrm{FL}}/2$ to $f_{\textrm{FL}}$ ($f_{\textrm{FL}}$ being the FL’s natural FMR frequency).

IV.1 DC response

Effect of the ST parameters. The spin-torque parameters appear naturally in the derivation of Slonczewski’s approximation for asymmetric ferromagnetic/non-ferromagnetic/ferromagnetic multilayers [35]. Experiments and theory have provided estimated values for both $P$ and $\Lambda$ [25, 36, 37]. Because of the uncertainty in both parameters, $P=0.35$ and $\Lambda=1.5$ are only first estimates.

The uncertainty in the ST parameters provided the main degree of freedom in achieving sufficient qualitative agreement between the experimental and simulated DC responses. The characteristic “foldover” FL FMR profiles evolve with increasing spin polarization factors $P$ and asymmetry parameters $\Lambda$ (Fig. 11). Moreover, higher $P$ and $\Lambda$ significantly broaden the FL mode’s linewidth. We estimated our samples’ ST parameters by matching the FL mode’s linewidths in the experimental and simulated DC responses. The ${\sim}0.5~{}\textrm{GHz}$ experimental linewidth corresponds to $P=0.4$ and $\Lambda=4$. These values are higher than our first estimates provided in the published works: $P=0.35$ and $\Lambda=1.5$. Such an increase, however, is necessary to match the experimental FL mode’s linewidth along with achieving the measurable DC response at $f_{\textrm{FL}}/2$.

When excited around a fractional frequency of the FL mode, increasing ST parameters facilitate the strongly nonlinear magnetization dynamics accompanied by the generation of the higher-order harmonics and measurable DC (e.g., as in Fig. 8(a)). The two contributions to the total DC response contain additional physical insights. Figures 14 and 15 in Appendix A show that the DC response at the FL FMR frequency is primarily determined by mixing of the TMR oscillations with the excitation signal [“oscillating” term in Eq. (4)]. On the other hand, the DC response at 1/2 the frequency of the FL FMR mode is defined by the “constant” term.

Response at the FL mode’s “sub-harmonics.” The magnetodipolar interaction between the FL and PL2 causes positive feedback, which qualitatively affects the system dynamics. This feedback can be thought of as follows. The deviation of the FL magnetization from its preferred orientation generates a stray field. This field causes the PL2 magnetization to deviate in the direction opposite to that of the FL. This, in turn, results in an even larger deviation of the FL magnetization due to the influence of the PL2 stray field, thus producing positive feedback between the magnetization dynamics of the FL and PL2 [38]

The FL precession exhibits an elliptical trajectory. Thus, the $m_{x}$ component moves back and forth twice during one oscillation cycle, hence its oscillation frequency is approximately twice that of the $m_{y}$ component [39]. Thus, if $f_{\textrm{FL}}$ is the FL’s resonant precession frequency, then $f(m_{x})=2f(m_{y})=2f_{\textrm{FL}}$. In this physical picture, the magnetodipolar interaction induces coupling between the FL and PL2. Due to the large $y$-component of the FL’s stray field induced by strong $m_{y}$ oscillations with the frequency $f_{\textrm{FL}}$, the above-mentioned coupling induces the $m_{y}$ oscillations of the PL2 with the same frequency $f_{\textrm{FL}}$. These oscillations, in turn, result in the $m_{x}$ oscillations of the PL2 with the frequency $f_{\textrm{FL}}/2$ (Fig. 12). The feedback via the PL2’s magnetodipolar field leads to the oscillation of the FL with the frequency $f_{\textrm{FL}}/2$ as well as all its harmonics. Thus, the FL power spectrum would contain the FL mode, its “sub-harmonic” at $f_{\textrm{FL}}/2$, and their higher-order modes, all of which are produced due to the nonlinear nature of the FL’s $m_{y}$ oscillations and the contribution of the PL2’s $m_{x}$ and $m_{y}$ oscillations excited by the presence of the magnetodipolar interaction.

Therefore, to facilitate the DC response at 1/2, 1/3, 1/4, and 1/6 the frequency of the FL FMR mode, it is enough to combine the magnetodipolar feedback between the FL and PL2 with a low-order nonlinearity. The DC response at 1/5 the frequency of the FL FMR mode is not present in the measurement: it requires a fifth-order nonlinearity, but the fifth harmonic is usually too small to evoke phase locking.

IV.2 Harmonic response

Figure 13 shows the fundamental and second harmonic of TMR counterposed to the total DC response. Whereas the TMR response at 1/2 the FL FMR frequency $f_{\textrm{FL}}$ has a weak fundamental and a strong second harmonic, the TMR response at the FL FMR frequency has a strong fundamental and a weak second harmonic.

The former observation implies that the MTJ’s nonlinear characteristics give rise to nonlinear oscillations under the AC excitation signal with the frequency corresponding to 1/2 the FL natural precession frequency. As the TMR response at the second harmonic is stronger than at the fundamental, the FL FMR mode seems to phase lock to this harmonic. As for the latter observation, the FL FMR mode phase locks to the TMR oscillations at the excitation frequency, thus producing a strong fundamental at the FL FMR frequency.

The major discrepancy between the measurements in Figs. 5(a) and 5(b) and simulations in Fig. 13 is the absence of the strong fundamental of $B_{1}$ at the FL FMR frequency $f_{\textrm{FL}}$. We attribute this difference to weak coupling (due to convergence to steady-state phase locking) between the NVNA’s incident wave and magnetization precession. We suggest that the experimental magnetization precession has a less stable phase relation with the NVNA’s incident signal than in the model.

In Fig. 5(a), the ripply response in the fundamental of $B_{1}$ suggests that there are more peaks at fractional frequencies of the FL mode than what we have identified. In particular, the enlarged DC responses in Fig. 10 indicate a resonant feature between 1/2 the FL FMR mode and the FL FMR mode. We qualitatively replicated this feature in the simulated TMR response [Figs. 13(a) (inset) and 13(b)]. Contrary to the experiment, however, the FL FMR mode seems to phase lock to the second harmonic, not to the TMR oscillations at the excitation frequency corresponding to 2/3 the FL FMR frequency $f_{\textrm{FL}}$. This is clearly seen in the time- and frequency-domain TMR responses for two selected excitation frequencies: 2/3 the frequency of the FL FMR mode and 5 GHz (Fig. 17 in Appendix A). The magnetic system’s harmonic response (specifically, of the second harmonic) is stronger if the excitation frequency corresponds to a fractional frequency of the FL FMR mode.

V Conclusion

We discovered that all characterized magnetic sensors’ DC responses reveal peaks at frequencies that are the integer fractions (1/2, 1/3, 1/4, and 1/6) of the devices’ natural FL FMR frequency $f_{\textrm{FL}}$. These peaks, in turn, generate the corresponding second and third harmonics of $B_{1}$. To understand the underlying physics that enabled the DC response at “sub-harmonics” of the FL mode, we employed micromagnetic modeling.

A comprehensive micromagnetic study suggested that the experimentally observed DC response at 1/2, 1/3, 1/4, and 1/6 of $f_{\textrm{FL}}$ can be defined by a low-order nonlinearity and strong magnetodipolar feedback between the FL and PL. As the PL is significantly thinner than the FL, additionally accounting for the ST effect on this layer notably enhanced the ST-driven harmonic response. Interestingly, the orthogonality of the FL and PL also facilitated the magnetic sensors’ response at 1/2 the FL FMR frequency $f_{\textrm{FL}}$.

Most importantly, strong magnetodipolar feedback permitted sub-harmonic injection locking within a wide range of integer fractions, which can be used in the development of a new generation of frequency multipliers.

Appendix A