Data Loading...
Large-amplitude coherent spin waves excited by spin ... Flipbook PDF
Large-amplitude coherent spin waves excited by spin-polarized current in nanoscale spin valves I. N. Krivorotov Departme
109 Views
102 Downloads
FLIP PDF 738.97KB
PHYSICAL REVIEW B 76, 024418 共2007兲
Large-amplitude coherent spin waves excited by spin-polarized current in nanoscale spin valves I. N. Krivorotov Department of Physics and Astronomy, University of California, Irvine, California 92697-4575, USA
D. V. Berkov and N. L. Gorn Innovent Technology Development, Pruessingstraße 27B, D-07745 Jena, Germany
N. C. Emley, J. C. Sankey, D. C. Ralph, and R. A. Buhrman Cornell University, Ithaca, New York 14853-2501, USA 共Received 17 March 2007; revised manuscript received 27 May 2007; published 13 July 2007兲 We present spectral measurements of spin-wave excitations driven by direct spin-polarized current in a free layer of nanoscale Ir20Mn80 / Ni80Fe20 / Cu/ Ni80Fe20 spin valves. The measurements reveal that large-amplitude coherent spin-wave modes are excited over a wide range of bias current. The frequency of these excitations exhibits a series of jumps as a function of current due to transitions between different localized nonlinear spin-wave modes of the Ni80Fe20 nanomagnet. We find that micromagnetic simulations employing the LandauLifshitz-Gilbert equation of motion augmented by the Slonczewski spin-torque term 共LLGS兲 accurately describe the frequency of the current-driven excitations including the mode transition behavior. However, LLGS simulations give qualitatively incorrect predictions for the amplitude of excited spin waves as a function of current. DOI: 10.1103/PhysRevB.76.024418
PACS number共s兲: 72.25.Ba, 75.75.⫹a, 75.20.⫺g
I. INTRODUCTION
The recent discovery of persistent current-driven excitations of magnetization in magnetic nanostructures1–13 has created new opportunities for studies of magnetization dynamics in extremely nonlinear regimes inaccessible with conventional techniques such as ferromagnetic resonance 共FMR兲. It was recently demonstrated14 that a spin-polarized current can excite motion of magnetization in metallic nanomagnets with precession cone angles over 30°, values far exceeding those achievable in typical FMR experiments performed on bulk and thin-film samples. There are two reasons why it is possible to have such large-amplitude currentdriven motions of magnetization in nanomagnets: 共i兲 suppression of Suhl instability processes15,16 due to quantization of the magnon spectrum in the nanomagnet17–28 and 共ii兲 efficient amplification of spin waves by spin-transfer torque that can act approximately as negative magnetic damping.1,2 The possibility of exciting large-amplitude oscillations of magnetization in magnetic nanostructures by spin-polarized current provides a unique testing ground for theories of nonlinear magnetization dynamics in ferromagnetic metals.29–32 Most importantly, it gives an opportunity to test the validity of the Landau-Lifshitz-Gilbert 共LLG兲 equation for the description of large-amplitude motion of magnetization. The LLG equation is phenomenological in nature, and thus its applicability must be tested in every new type of experimental situation. This equation has proved to be largely successful in the description of persistent small-angle magnetic excitations33 and transient large-angle magnetization dynamics34,35 in thin films of ferromagnetic metals 共with some notable exceptions36兲. However, it is not known a priori that the LLG equation is suitable for the quantitative description of a persistent magnetization precession with very large amplitude. For example, the phenomenological Gilbert damping term parametrized by a single constant in 1098-0121/2007/76共2兲/024418共14兲
the LLG equation may prove to be an approximation suitable for description of small-angle dynamics but not valid in general. Recently, large-angle persistent motion of magnetization was studied in thin films of Ni80Fe20 by time-resolved measurements, and a large increase of apparent damping was observed in the nonlinear regime.37 However, measurements of intrinsic damping in continuous ferromagnetic films are obscured by generation of parametrically excited spin waves that give rise to at least a large portion of the increased damping found in Ref. 37. This generation of parametrically pumped spin waves is expected to be suppressed in nanoscale ferromagnets,23 and thus information on the amplitude dependence of intrinsic damping can, in principle, be accessed. A number of recent models predict nontrivial angular dependence of damping38,39 and suggest how it may depend on the rate and amplitude of magnetization precession40 in metallic magnetic nanostructures. These predictions remain largely untested primarily due to the difficulty of exciting persistent large-amplitude magnetization dynamics in nanomagnets. In this work we report a detailed comparison of experimentally measured spectra of current-driven magnetization oscillations in elliptical Py 共Py⬅ Ni80Fe20兲 nanoelements to the results of full-scale micromagnetic simulations for these structures and thus test the validity of the micromagnetic LLG approach for the description of strongly nonlinear oscillations of magnetization in magnetic nanostructures. We find that although simulations based on the LLG equation augmented by the Slonczewski spin torque term1 共the LLGS equation兲 can successfully mimic many properties of currentdriven magnetization dynamics such as the current dependence of the excitation frequency and abrupt frequency jumps with increasing current, they qualitatively fail to reproduce the dependence of the amplitude of current-driven spin waves as a function of current. Our results demonstrate the deficiencies of the current LLGS implementation for the
024418-1
©2007 The American Physical Society
PHYSICAL REVIEW B 76, 024418 共2007兲
KRIVOROTOV et al.
description of spin-torque-driven magnetization dynamics and suggest the need for modification of this implementation for a quantitative description of large-amplitude magnetization motion. We suggest that it may be necessary to introduce a nonlinear dissipation or to consider effects of spin transfer from lateral spin diffusion that are not contained in our calculation.
(a)
Py
Cu
MP
Py
IrMn
θp
Cu
(c)
1
(d)
0.8
D R/D Rmax
R (Ohm)
5.95
5.9
0.6 0.4
5.85
0.2 0
5.8 -700
-350
0
350
700
-700
-350
0
350
700
Magnetic Field (Oe)
Magnetic Field (Oe) 6.15
(e)
6.3
(f) 0 Oe
6.05
R (Ohm)
dV/dI (Ohm)
The current-perpendicular-to-plane 共CPP兲 nanopillar spin valves for our experiments are prepared by magnetron sputtering of continuous magnetic multilayers onto an oxidized Si wafer followed by a multistep nanofabrication process.7 As a first step of the sample preparation process, a multilayer of Cu共80 nm兲 / Ir20Mn80 共8 nm兲Py共4 nm兲 / Cu 共8 nm兲 / Py 共4nm兲 / Cu共20 nm兲 / Pt共30 nm兲 is deposited onto a thermally oxidized Si 共100兲 wafer by magnetron sputtering in a highvacuum system with a base pressure of 2 10−8 Torr. The 80-nm Cu layer is used as the bottom electrode of the CPP spin valve. The Pt capping layer is employed for protection of the multilayer from oxidation during the nanopillar fabrication process. The multilayer is deposited at room temperature in a magnetic field of approximately 500 Oe applied in the plane of the sample and post-annealed at T = 250 ° C for 80 min in the same field. We use a subtractive process employing e-beam lithography, photolithography, and etching of the multilayer in order to define nanoscale spin valves of approximately elliptical shape with major and minor axes of 130 nm and 60 nm, respectively, and with Cu electrodes making contact with the top and bottom of the spin valve as shown in Fig. 1共a兲. The role of the antiferromagnetic Ir20Mn80 layer in the spin-valve structure is twofold: 共i兲 to pin the direction of magnetization of the fixed Py nanomagnet at a nonzero angle with respect to the easy axis of the free nanomagnet using the exchange bias effect as shown in Fig. 1共b兲 and 共ii兲 to suppress current-driven excitations of magnetization in the fixed nanomagnet due to the giant enhancement of Gilbert damping observed in exchange-biased ferromagnets.41,42. The nominal direction of the exchange bias field set during the multilayer deposition and subsequent annealing is in the plane of the sample at 45° with respect to the major axis of the ellipse. However, within a set of 40 samples we found significant 共±35° 兲 sample-to-sample variations of the exchange bias direction, as determined from the StonerWohlfarth 共SW兲 fitting procedure described below. These sample-to-sample variations of the exchange bias direction are not surprising in a magnetic nanostructure and may be attributed to finite-size effects43–45 as well as to resetting of the exchange bias direction due to sample heating that occurs during lithography and ion milling process employed to define the nanopillar structure. In this paper we report experimental results for the most extensively studied sample although qualitatively similar results were obtained for the other samples from the set of 40. The quantitative differences between the samples can be correlated with differences of the shapes of the hysteresis loop of resistance versus field,
MF
H
6
II. EXPERIMENT A. Sample preparation and characterization
HEB
(b)
Cu
6.1
0 Oe
680 Oe
5.95
680 Oe
5.9
5.85
5.7
5.75 -10
-5
0
Current (mA)
5
10
-10
-5
0
5
10
Current (mA)
FIG. 1. 共Color online兲 共a兲 Schematic side view of the nanopillar spin valve used for studies of magnetization dynamics. 共b兲 Schematic top view of the spin valve with approximate directions of magnetizations of the pinned, M P, and the free, M F, nanomagnets as well as the direction of positive external magnetic field H and exchange bias field HEB. 共c兲 Experimentally measured resistance of the nanopillar as a function of the external magnetic field 共circles兲 and a macrospin Stoner-Wohlfarth fit to the data 共solid line兲 with the parameters described in text. 共d兲 Resistance versus field obtained from micromagnetic simulations using the giant magnetoresistance 共GMR兲 asymmetry parameter = 0.5, the exchange bias field magnitude HEB = 1600 Oe and its direction EB = 30° obtained from the macrospin fit shown in 共c兲. 共e兲 Differential resistance of the sample as a function of bias current measured at H = 0 Oe 共red兲 and H = 680 Oe 共blue兲. 共f兲 dc resistance of the sample as a function of bias current obtained from the data in 共e兲 by numerical integration.
such as that shown in Fig. 1共c兲, and ultimately to variations of the direction of the exchange bias field. Samples with similar resistance versus field hysteresis loops exhibit similar spectral properties of the current driven magnetization oscillations. All measurements reported in this paper were made at T = 4.2 K. We determine the direction and magnitude of the exchange bias field for each nanopillar sample by fitting the Stoner-Wohlfarth model to the experimental resistance– versus–magnetic-field hysteresis loop, such as that shown in Fig. 1共c兲. For all measurements and simulations reported in this paper, we apply the external magnetic field in the plane of the sample at 45° with respect to the ellipse major axis and approximately perpendicular to the exchange bias direction as shown in Fig. 1共b兲. Stoner-Wohlfarth simulations
024418-2
PHYSICAL REVIEW B 76, 024418 共2007兲
LARGE-AMPLITUDE COHERENT SPIN WAVES EXCITED…
show that this choice of the bias field direction results in a weak dependence on the magnitude of external magnetic field for the equilibrium angle between magnetic moments of the free and pinned layers. According to Stoner-Wohlfarth simulations, the equilibrium angle between magnetic moments of the free and the pinned layer varies between 34° and 36° in the field range from 300 Oe to 1100 Oe. The solid line in Fig. 1共c兲 is a four-parameter Stoner-Wohlfarth fit to the data with the following fitting parameters: the exchange bias field magnitude HEB, its direction EB, the magnetoresistance 共MR兲 asymmetry , and the MR magnitude ⌬R. The MR asymmetry parameter 共Refs. 46 and 47兲 describes a deviation of the angular dependence of the giant magnetoresistance 共GMR兲 from a simple cosine form: R共兲 = R0 + ⌬R
1 − cos2共/2兲 . 1 + cos2共/2兲
共1兲
Here is the angle between magnetic moments of the pinned and the free layers. The Stoner-Wohlfarth fit shown in Fig. 1共c兲 yields HEB = 1.6± 0.5 kG, EB = 共30± 6兲°, = 0.5± 0.3, and ⌬R = 0.161± 0.007 ⍀. Two other parameters used in the Stoner-Wohlfarth simulations are the uniaxial shape anisotropy field HK of the elliptical Py nanomagnets and the average dipolar coupling field between the fixed and the pinned layers, Hdip. The value of HK = 600 Oe was obtained as the saturation field along the in-plane hard axis of the nanomagnet by employing micromagnetic simulations 共OOMMF兲.48 The value of Hdip = 80 Oe was obtained by numerical integration of the dipolar coupling energy of the two uniformly magnetized Py nanomagnets. The value of obtained from our fitting procedure is significantly less than that reported for a similar structure in Ref. 47 共 ⬇ 2兲. The difference is probably due to the different values of the effective Py/Cu interfacial and Py bulk resistances in our spin valves, possibly due to interdiffusion of metallic layers of the spin valve during the annealing process.49 To test the validity of the Stoner-Wohlfarth approach for fitting the quasistatic MR hysteresis loop, we calculate the MR loop for this sample by employing full micromagnetic simulations50 with the values of HEB, EB, and obtained from the SW fit. Other input parameters for micromagnetic simulations were obtained by direct measurements. The saturation magnetization M S of a 4-nm-thick Py film sandwiched between two Cu films and subjected to the same heat treatment as the spin valves under study was measured by superconducting quantum interference device 共SQUID兲 magnetometry and was found to be M S = 650 emu/ cm3 at T = 5 K. The Gilbert damping parameter = 0.025 共needed for the dynamic simulations described in Sec. III below兲 for these samples was measured at T = 7 K by a pump-probe technique described in Ref. 14. This value significantly exceeds that measured for Py nanomagnets at room temperature by spintorque FMR spectroscopy 共 = 0.01兲.51 A similar increase of the Gilbert damping parameter in nanoscale Py elements at low temperature was previously observed in Ref. 52 and was attributed to exchange coupling of the nanomagnet to an antiferromagnetic oxide layer formed along the perimeter of the nanomagnet.
The result of the full-scale micromagnetic simulation is shown in Fig. 1共d兲. We find that the SW model is a reasonable approximation for the quasistatic hysteresis loop in that the coercivity predicted by micromagnetic simulation is ⬃80% of that given by the SW model and the shapes of the Stoner-Wohlfarth and micromagnetic hysteresis loops are similar. However, we could not obtain a quantitatively correct fit of the measured GMR loop using full-scale micromagnetic simulations. The main difference was that the simulated loops were always narrower than the loop measured experimentally 关Fig. 1共c兲兴: the difference between the left and right coercive fields for the loop shown in Fig. 1共c兲 is ⌬Hc = 680 Oe while the maximum ⌬Hc = 600 Oe was obtained for various directions of the exchange bias field that we have used in our simulations. Regarding this discrepancy we note that in full-scale simulations we do not have at our disposal the anisotropy field HK and the dipolar coupling Hdip as adjustable parameters—the corresponding effective field contributions are calculated from the material saturation magnetization and the sample geometry. Taking into account that the width of the GMR hysteresis loops varied substantially among different samples as discussed in Sec. IV, we do not consider this discrepancy as being qualitatively significant. The difference between the SW simulations and fullscale micromagnetic modeling means, first, that the fit parameters obtained from the SW approximation should be considered not as exact values, but rather as reasonable guesses, and second, that some magnetic properties of the system under study 共e.g., surface anisotropy, sample shape imperfections, and the possible presence of antiferromagnetic oxides along the perimeter of the free layer nanomagnet52兲 are still not included in our model. Because the main goal of this paper is the study of dynamic system properties, we postpone the discussion of this quite interesting problem to future publications. Figure 1共e兲 shows the measured differential resistance of the sample as a function of direct current flowing through the sample for H = 0 Oe and for H = 680 Oe. Positive current in this and subsequent figures corresponds to the flow of electrons from the free to the pinned layer. Figure 1共f兲 shows the dc sample resistance R = V / I for H = 0 Oe and H = 680 Oe as a function of direct current, obtained by numerical integration of the differential resistance data in Fig. 1共e兲. The quasiparabolic increase of the resistance with increasing current 共most clearly seen for negative currents兲 is due to a combination of Ohmic heating53 and the Peltier effect in the nanopillar junction.54 The other features in the plots of dV / dI versus I and R共I兲 such as hysteretic switching of resistance at H = 0 Oe or peaks in the differential resistance at H = 680 Oe are due to changes of magnetic state of the nanopillar. At fields below the coercive field of the free layer, we observe current-induced hysteretic switching between the low- and high-resistance states. For fields exceeding the coercive field, the time-averaged resistance of the sample R共I兲 undergoes a transition from the low-resistance state to an intermediate-resistance state under the action of direct current as shown in Fig. 1共f兲 关e.g., R共680 Oe兲 = R共0 Oe兲 − 0.27⌬R for I = 10 mA兴. As we demonstrate below, this intermediate-resistance state is a state of persistent currentdriven magnetization dynamics for the free nanomagnet.
024418-3
PHYSICAL REVIEW B 76, 024418 共2007兲
KRIVOROTOV et al. (a)
spectral density V共f兲 of the GMR signal due to oscillations of magnetization at the nanopillar as V共f兲 = 共RS + R0兲冑P共f兲 / R0.55 In this expression, R0 = 50 ⍀ is the characteristic impedance of all components of the microwave circuit shown in Fig. 2共a兲 except for the nanopillar itself, and RS = 26 ⍀ is the resistance of the nanopillar junction and leads. We define the normalized rms amplitude spectral density S共f兲 as the rms voltage spectral density V共f兲 divided by the maximum rms GMR voltage signal amplitude 具Vmax典 共where = 2 f兲 = 冑具关共I⌬R / 2兲sin共t兲兴2典 / 2 = I⌬R / 2冑2 achievable due to 360° uniform rotation of magnetization in the sample plane at a given current bias:
(b)
(c)
(d)
S共f兲 =
V共f兲 V共f兲 = 冑8 . 具Vmax典 I⌬R
共2兲
The dimensionless integrated signal amplitude Sint, Sint =
冑冕
⬁
S共f兲2df ,
共3兲
0
FIG. 2. 共Color online兲 共a兲 Circuit schematic for measurements of magnetization dynamics driven by a direct current. 共b兲 Normalized rms amplitude spectral density S共f兲 共defined in text兲 generated by the spin valve under a dc bias of 6.15 mA. 共c兲 S共f兲 as a function of current for the nanopillar spin valve measured at H = 680 Oe. 共d兲 S共f兲 at 3.7 mA and H = 680 Oe. B. Measurements of current-driven oscillations of magnetization
To measure the current-driven excitations of magnetization directly, we employ a spectroscopic technique developed in Ref. 10. Figure 2共a兲 schematically shows the measurement setup employed for detection of the current-driven excitations of magnetization. In this setup, direct current flowing perpendicular to the layers of the spin valve sample excites magnetization oscillations in the free Py nanomagnet, which give rise to a temporal variation in the resistance of the spin valve due to the GMR effect, R共t兲. Since the sample is current biased 共Idc兲, the temporal variation of the resistance generates on ac voltage, V共t兲 = IdcR共t兲. This ac voltage is amplified with a microwave signal amplifier, and its spectral content is recorded with a spectrum analyzer. A spectrum measured at zero-dc-bias current is subtracted from all spectra in order to eliminate a small background due to thermal and electronics noise. Figures 2共b兲 and 2共d兲 show representative examples of typical spectra generated by the spin valve under direct current bias. The signals shown in Figs. 2共b兲 and 2共d兲 are the normalized rms amplitude spectral density S共f兲 defined below. This quantity characterizes the amplitude, frequency, and coherence of oscillations of magnetization. To calculate S共f兲, we start with the power spectral density measured with the spectrum analyzer, Pan共f兲. This quantity is corrected for frequency-dependent amplification and attenuation in the circuit between the spectrum analyzer and the nanopillar sample in order to obtain the power spectral density P共f兲 of the signal emitted by the sample into a 50-⍀ transmission line. This latter quantity is used to calculate the rms voltage
reaches its maximum value Sint = 1 / 冑2 for the maximum possible GMR voltage signal due to 360° uniform rotation of magnetization in the sample plane: Vmax共t兲 =
I⌬R sin共2 ft兲. 2
共4兲
The integrated signal amplitude Sint is a convenient dimensionless scalar quantity that characterizes the amplitude of magnetization precession. Its square is directly proportional to the integrated power emitted by the device. This dimensionless quantity is also convenient for comparison of experimental data to the results of micromagnetic simulations. A typical experimentally measured spectrum S共f兲 for our samples is characterized by a single frequency 关the fundamental peak and higher harmonics such as that shown in Fig. 2共b兲兴. However, for some values of the bias current, two peaks that are not harmonically related to each other are observed 关Fig. 2共d兲兴. Figure 2共c兲 shows a summary of spectra generated by the sample as a function of the direct current bias Idc measured at a fixed value of the applied magnetic field H = 680 Oe. The most important features of these data are the following. 共i兲 The frequency of the current-driven excitations decreases with increasing current. This decrease of frequency with increasing current can be explained as a nonlinear effect arising from the dependence of the frequency of precessing magnetization on the precession amplitude.10,29,30,56,57 共ii兲 The frequency of the current-driven excitations exhibits downward jumps at I ⬇ 3.7 mA and 4.85 mA. The current values at which the frequency jumps occur coincide with the positions of the peaks in the plot of differential resistance versus current 关Fig. 1共e兲兴. A double-peak structure in the spectrum such as that shown in Fig. 2共d兲 is observed only for currents near frequency jumps, indicating that the apparent jumps are in fact nonhysteretic crossovers between two excitations with different frequencies. As the current is increased across the transition region, the emitted power is
024418-4
PHYSICAL REVIEW B 76, 024418 共2007兲
LARGE-AMPLITUDE COHERENT SPIN WAVES EXCITED…
gradually transferred from the excitation with the higher frequency to the excitation with the lower frequency. We also observe that the linewidths of the current-driven excitations increase in the current intervals where two excitations coexist 关e.g., compare Figs. 2共b兲 and 2共d兲兴. In the current intervals where a single large-amplitude mode is excited, spectral lines as narrow as 10 MHz are observed while for currents in the mode transition regions spectral lines as wide as 250 MHz are found. The increase of the linewidth of the excitation indicates the decrease of its temporal coherence.57–60 The linewidth increase is observed in all transition regions suggesting that the decrease of coherence of the current-driven spin waves is induced by interaction between the two excited spin-wave modes. 共iii兲 Modes with very low power visible only on the logarithmic amplitude scale of Fig. 2共c兲 are observed for currents above 3.7 mA. These modes are not harmonically related to the dominant modes. Although these modes emit low integrated power, they may play an important role in determining the coherence of the dominant spin-wave excitations.61
共p = 130° 兲 values of the equilibrium angle between magnetization and current polarization to study the effect of the spin polarization direction on the current-driven dynamics. The decisive advantage of numerical simulations is the possibility to study and understand the influence of all relevant physical factors separately. For this reason we start from the “minimal model,” where the influence of the Oersted field and thermal fluctuations is neglected and the spin torque is assumed to be symmetric 关⌳ = 1 in Eq. 共5兲兴, and then switch on in succession all the factors listed above to analyze their influence on the magnetization dynamics. Minimal model. Results for this model, for which the Oersted field and thermal fluctuations are not included and the asymmetry parameter ⌳ = 1 关f J共兲 = const兴, are presented in Fig. 3共a兲. For the spin orientation angle p = 150° used for the first simulation series, the critical spin-torque value for the oscillation onset was found to be acr J ⬇ 0.308共2兲. Using the simplest expression for the spin torque given, e.g., in Ref. 62, one can easily derive the relation between the reduced spin-torque amplitude aJ and other device parameters as aJ =
III. NUMERICAL SIMULATIONS
Full-scale micromagnetic simulations of the currentinduced magnetization dynamics in the nanopillar described above were performed by solving the stochastic LLG equation of motion for the magnetization M共r兲 using our commercially available simulation package MICROMAGUS 共Ref. 50兲 supplemented by a spin injection module. Details of the simulation methodology are given in the Appendix. In these simulations, the spin torque in the LLGS simulations has the form ⌫ = −f J共兲关M ⫻ 关M ⫻ p兴兴 where the dimensionless spintorque amplitude f J depends on the angle between the magnetization M and the unit vector p of the polarization direction of the electron magnetic moments 共in the spin-polarized current兲. We use, in general, the asymmetric angular dependence of the spin torque amplitude f J共兲, f J共 兲 = a J ·
2⌳2 , 共⌳2 + 1兲 + 共⌳2 − 1兲cos
共5兲
given in Refs. 46 and 62, where ⌳ is the asymmetry parameter related to the GMR asymmetry parameter from Eq. 共1兲 via ⌳2 = 1 + . As will be demonstrated below, variations in the direction of the spin current polarization p 共opposite to the magnetization of the pinned layer MP兲 in the spin torque term ⌫ = −f J共兲关M ⫻ 关M ⫻ p兴兴 can result in qualitative changes of magnetization dynamics and thus the orientation of p is a very important parameter of the problem under study. As explained in Sec. II, the direction of p could not be determined quantitatively from the available experimental data. To understand the dependence of the magnetization dynamics on the orientation of p, we first study magnetization dynamics for p directed opposite to the exchange bias field extracted from the GMR hysteresis fit as described in Sec. II 共i.e., the angle between p and the positive direction of the x axis was set to p = 150°兲. Then we perform two additional simulation sets for larger 共p = 170° 兲 and smaller
ប j 1 1 P , 2 兩e兩 d M 2S
共6兲
where e is the electron charge, j is the electric current density, d is the thickness of a magnetic layer subject to a spin torque, and P is the degree of spin polarization of the electrical current. Using the definition of the current density, j = I / Selem 共I is the total current and Selem is the area of the nanopillar cross section兲, and substituting the values for the experimentally measured critical current Icr ⬇ 2.7 mA and the threshold for the oscillation onset acr J ⬇ 0.3 found in simulations, we obtain that the polarization degree of the electron magnetic moments is P ⬇ 0.32. From the relation between the critical current Icr and the critical spin torque amplitude acr J , the proportionality factor between the spin-torque amplitude aJ used in the simulation and the experimental current strength I 共in mA兲 is ⬇ 0.11 共mA−1兲 共whereby aJ = I兲 for this spin polarization direction 共 p = 150° 兲. The simulated spectral lines are very narrow 共mostly ⬍100 MHz兲 for all values of the spin-torque amplitude acr J 艋 aJ 艋 2.0, which means that for this simplest model a transition to a quasichaotic regime similar to that found in Ref. 63 does not occur in the interval of currents studied. For this reason we show in Fig. 3共a兲 only the positions of the spectral maxima of the M z component as a function of aJ 共red circles兲. In addition to the narrow lines, this minimal model also reproduces two other important qualitative features of the experimental results 关see Fig. 2共c兲兴: a rapid decrease of the oscillation frequency with increasing current immediately after the oscillation onset and 共ii兲 downward frequency jumps at higher current values. The first feature, the rapid decrease of the oscillation frequency immediately after the oscillation onset, is a nonlinear effect due to the rapid growth of the oscillation amplitude with increasing current. In the nonlinear regime, the frequency decreases with increasing amplitude because the length of the precession orbit grows faster than the magnetization velocity. The corresponding effect was obtained ana-
024418-5
PHYSICAL REVIEW B 76, 024418 共2007兲
KRIVOROTOV et al.
FIG. 3. 共Color online兲 共a兲 Dependence of the frequency of magnetization oscillations, f, on spin-torque amplitude aJ calculated in the “minimal model” 共see text for details兲. Gray-scale maps represent spatial distributions of the oscillation power for chosen aJ values 共bright corresponds to maximal oscillation power兲. 共b兲 Oersted field effect: f共aJ兲 without 共open circles兲 and with 共open triangles兲 the Oersted field included in the simulations; arrows indicate the positions of frequency jumps, and straight lines are guides to the eye. 共c兲 Effect of the torque asymmetry: f共aJ兲 for the symmetric 共⌳ = 1.0, open triangles兲 and asymmetric 共⌳2 = 1.5, solid triangles兲 spin torque.
lytically in Refs. 29, 30, and 57 and observed numerically in our full-scale micromagnetic simulations63 of the experiments published in Ref. 10 The second important observation, the existence of downward frequency jumps with increasing current, cannot yet be explained using analytical theories, and such jumps are absent in the macrospin description of current-driven magnetization dynamics. Spatially resolved spectral analysis of our simulation data reveals that these jumps correspond to transitions between strongly nonlinear oscillation modes 关see spatial maps of the magnetization oscillation power in Fig. 3共a兲兴. With each frequency jump, the mode becomes more localized but the oscillation power is still concentrated in one single-connected spatial region which has no node lines. We
discuss these modes in more detail below when analyzing results for different current polarization directions. An analytical theory of the nonlinear eigenmodes of a resonator having the correct shape would be required to achieve a thorough understanding of this phenomenon. In the minimal model we also observe for the current strength aJ ⬎ 1.5 the so-called “out-of-plane” coherent precession regime for which the magnetization acquires a nonzero time-average component perpendicular to the sample plane. This regime is characterized by frequency increasing with current and is well known from analytical consideration and numerical simulations.64,65 The out-of-plane regime was experimentally observed for a nanopillar sample with a 2-nm-thick free Py layer and low critical current, Icr = 1.4 mA, and thus relatively small Oersted field.66 However, for the devices with a 4-nm-thick free layer that we study, this type of mode is an artifact of the minimal model 共due to the absence of the Oersted field兲 and it was not observed experimentally. Effect of the Oersted field. The effect of the Oersted field is demonstrated in Fig. 3共b兲 where the dependences of the oscillation frequency on the spin-torque magnitude aJ are shown without 关red circles, identical to Fig. 3共a兲兴 and with the Oersted field 共green triangles兲. To compute the Oersted field, we have used the proportionality constant between the spin-torque magnitude aJ and experimental current value I 共in mA兲 assuming that the simulated threshold value acr J ⬇ 0.31 corresponds to the experimentally measured critical current Icr ⬇ 2.7 mA. The results shown in Fig. 3共b兲 demonstrate that the Oersted field has two major effects on magnetization dynamics. First, this field eliminates the out-of-plane precession: inspection of magnetization trajectories shows that for all aJ values they correspond to “in-plane” steady-state oscillations. As a consequence, the Oersted field eliminates the upward frequency jump in the f共aJ兲 dependence. The suppression of the out-of-plane mode occurs because the Oersted field is a strongly inhomogeneous in-plane field that keeps magnetization close to the plane of the sample. Second, the Oersted field shifts the frequency jumps to lower current values. This can be explained as follows: the Oersted field is highly inhomogeneous, with its maximal values at the edges of the elliptical element. For this reason it should suppress magnetization oscillations at the element edges, thus favoring spatial oscillation modes localized near the element center such as those shown in Fig. 3共a兲. Hence the transition from the homogeneous mode to more localized ones should occur for lower currents when the Oersted field is taken into account. Effect of the spin-torque asymmetry. It can be seen directly from Eq. 共5兲 that for / 2 ⬍ p ⬍ 3 / 2 the spin-torque magnitude in the case of positive GMR asymmetry 共 ⬎ 0, ⌳ ⬎ 1兲 is larger than for the symmetric 共 = 0, ⌳ = 1兲 case. This difference is expected to result in a decrease of the steady-state precession frequency at a given current value because larger spin torque results in larger amplitude of magnetization oscillations. For the system studied in this paper, the GMR asymmetry is relatively low 共the Stoner-Wohlfarth fit of the quasistatic GMR curve gave the value = 0.5 and ⌳2 = 1.5兲, so that the expected frequency decrease is quite
024418-6
PHYSICAL REVIEW B 76, 024418 共2007兲
LARGE-AMPLITUDE COHERENT SPIN WAVES EXCITED…
weak, but it is nevertheless clearly visible when comparing the f共aJ兲 dependences for symmetric and asymmetric cases in Fig. 3共c兲. The frequency decrease is greatest in the lowcurrent region where the dependence of the oscillation frequency on current is very steep. The asymmetric torque also leads to a small decrease of the threshold current for the oscillation onset, which becomes acr J ⬇ 0.27. A more important effect of introducing the spin torque asymmetry is a shift of the frequency values where the transitions between nonlinear eigenmodes of the system 共accompanied by the frequency jumps as explained above兲 take place. Even for the relatively low value ⌳2 = 1.5 used in our simulations, this shift is significant 关see Fig. 3共c兲兴. The reason, again, is that the asymmetric torque form gives a larger spin torque magnitude for a given current value, so that the transitions to more localized modes occur earlier. Influence of thermal fluctuations. To take into account the influence of thermal fluctuations, we have first to estimate the real temperature of the sample. Although the experiments were performed at liquid helium temperature T = 4.2 K, Joule heating of our multilayer nanoelement due to the direct current through the device was unavoidable.53 To estimate the maximal temperature of the nanoelement, we have measured 共i兲 the temperature dependence of its resistance R0共T兲 in the absence of any dc current by heating the whole setup and 共ii兲 the dependence of the resistance on current R0共I兲 measured at positive current and H = 680 Oe as shown in Fig. 1共f兲. The increase of the resistance with current is due both to the excitation of coherent magnetization oscillations and Ohmic heating; therefore, an estimate of the sample temperature from R0共I兲 gives an upper bound on the temperature of the sample. Comparison of R0共T兲 and R0共I兲 shows that the nanoelement temperature does not exceed T ⬇ 60 K for the highest current I = 10 mA used in the measurements. Taking into account that this maximal temperature is relatively low, we have simply adopted a linear interpolation between the lowest temperature T = 4 K for I = 0 mA and T ⬇ 60 K for I ⬇ 10 mA 共with I converted into the reduced spin torque amplitude aJ兲 for our simulations. The dependences of the excitation frequency on current f共aJ兲 for T = 0 and T ⬀ aJ with the proportionality factor calculated as explained above are compared in Fig. 4共a兲 共for both simulation sets the effect of the Oersted field and the spin torque asymmetry with ⌳2 = 1.5 are included兲. It is clear that due to the relatively low temperatures, thermal fluctuations have a minor effect both on the oscillation frequency and on the positions of the frequency jumps. Influence of the spin current polarization direction. The polarization direction p of the electron magnetic moments in the dc current is expected to be one of the most important parameters of the problem. First, the onset threshold for oscillations should depend strongly on this polarization direction.56,67 Second, the relative strength of the Oersted field 共with respect to the spin-torque magnitude aJ兲 also should depend on p, because the Oersted field for different p is computed assuming that the threshold value acr J always corresponds to the experimentally measured critical current Icr ⬇ 2.7 mA. Since p could not be accurately determined from the fit of the quasistatic MR hysteresis loop 共see dis-
(a)
(b)
(c)
(d)
FIG. 4. 共Color online兲 共a兲 Effect of thermal fluctuations on the frequency of magnetization oscillations: f共aJ兲 for T = 0 共triangles兲 and for T ⬀ I ⬀ aJ 共crosses兲. Simulated rms-amplitude spectral densities S共f兲 for oscillations 共b兲 before the first frequency jump, 共c兲 between the first and second frequency jumps, and 共d兲 after the second frequency jump. The linewidth of the peaks marked with ␦ is below the resolution limit of our numerical simulations. See the text for the detailed analysis of these spectra.
cussion in Sec. II A兲, we have carried out additional series of simulation runs to study the effect of the spin current polarization direction on the magnetization dynamics. The results of these simulations are summarized in Fig. 5, where we show the dependences of the oscillation frequency on the spin-torque magnitude f共aJ兲 关Fig. 5共a兲兴 and on the spin-torque magnitude normalized by the threshold value acr J for the corresponding angle f共aJ / acr J 兲 关Fig. 5共b兲兴. The dependence f共aJ兲 for p = 150°—i.e., for the case which a detailed analysis has been presented above—is shown in this figure with open circles. For the increased polarization orientation angle p = 170° 共open triangles兲 the onset threshold for the magnetization dynamics decreases cr from acr J 共 p = 150° 兲 ⬇ 0.27 to aJ 共 p = 170° 兲 ⬇ 0.15 in a qualitative agreement with Slonczewski’s prediction for the macrospin model 共Icr ⬃ 1 / 兩cos p兩兲 关Ref. 1兴 experimentally confirmed in Ref. 67. The first frequency jump with increasing current is still present, but instead of the second jump we observe a kink in the f共aJ兲 curve 关see Fig. 5共a兲兴. We note that
024418-7
PHYSICAL REVIEW B 76, 024418 共2007兲
KRIVOROTOV et al.
FIG. 5. 共Color online兲 共a兲 Dependence of the frequency of oscillations on spin torque amplitude aJ for different polarization angles of the spin current, p, as defined in Fig. 1共b兲. 共b兲 The same frequencies plotted as functions of the normalized spin torque amplitude aJ / acr J.
the importance of the Oersted field relative to the spin-torque effect in this case is much larger than for p = 150° for the following reason. The Oersted field is always computed assuming that the critical value of the spin-torque magnitude acr J corresponds to one and the same physical current value Icr ⬇ 2.7 mA. This means that if acr J decreases, the same Oersted field corresponds to smaller aJ values so that the importance of the Oersted field effect increases relative to the spintorque action. For a smaller polarization angle 共 p = 135°, results are shown on Fig. 5 with crosses兲 the critical value of aJ increases 关acr J 共 p = 135° 兲 ⬇ 0.66兴, so that the influence of the Oersted field is weaker than for p = 150°. This leads, in particular, to the reappearance of the out-of-plane oscillation regime, which manifests itself in the increase of the oscillation frequency with increasing aJ. Recall that the out-ofplane precession regime was found for p = 150° in the absence of the Oersted field, but was suppressed by this field as explained above 关see Fig. 2共b兲兴. For the angle p = 135° the Oersted field is not strong enough to eliminate this regime when the spin-torque magnitude increases. To compare magnetization dynamics for various spin polarization angles we plot the frequency of oscillations for all three values of p studied as a function of spin-torque mag-
nitude normalized to its critical value, aJ / acr J 关Fig. 5共b兲兴. The most striking feature of the f共aJ / acr J 兲 curves for various angles p is that they all nearly collapse onto the universal cr cr f共aJ / acr J 兲 dependence for aJ / aJ values up to aJ / aJ ⬇ 2. This region includes, in particular, the fast frequency decrease after the oscillation onset 共see the discussion of this nonlinear effect above兲 and the first frequency jump arising for all spin polarization directions at almost the same value of aJ / acr J ⬇ 1.5. These results clearly demonstrate that the initial nonlinear rapid frequency decrease and the first frequency jump are universal for the system under study, whereas further behavior of the magnetization dynamics 共in particular, the existence of the second frequency jump兲 are much more subtle features and thus may vary from sample to sample. The first frequency jump is always present because it marks the transition from the homogeneous to a localized oscillation mode 关see Fig. 3共a兲兴 which is always accompanied by an abrupt change of the oscillation frequency. The next frequency jump for the situation when the Oersted field is neglected corresponds to the transition between the modes with different 共but symmetric兲 localization patterns—before the second jump the mode is localized in the direction along the major ellipse axis only, whereas after this jump the new mode is confined in both directions 关compare second and third maps in Fig. 3共a兲兴. This latter transition is strongly disturbed by the Oersted field, which leads, in particular, to strongly asymmetric spatial mode patterns for localized modes. This may eliminate the qualitative differences between modes with different localization patterns that give rise to the frequency jumps. To understand the degree of agreement that may be expected between the simulation results and the experimental data for our samples, it is instructive to examine sample-tosample variations of experimentally measured magnetization dynamics for samples with different directions of the exchange bias field. Figure 6 shows resistance as a function of field 共left column兲 and the corresponding dependence of the frequency of the current-driven excitations on current 共right column兲 for three representative samples from the set of forty samples studied. These samples have different directions of the exchange bias field as determined from the Stoner-Wohlfarth fit of the hysteresis loop of resistance versus field 关共a兲,共b兲 EB = 22°, 共c兲,共d兲 EB = 38°, 共e兲,共f兲 EB = 48°兴. Figure 6 along with Figs. 1 and 2 illustrates typical sampleto-sample variations of the quasistatic hysteresis loop and the frequency of current-driven excitations among nominally identical samples. While it is clear that there are significant sample-to-sample variations of the current dependence of the frequency of excitations, the initial decrease of frequency with current as well as downward frequency jumps is always present in these data. We note that there are correlations between theoretically predicted trends in magnetization dynamics as a function of the exchange bias direction shown in Fig. 5 and experimental data such as those shown in Fig. 6. In particular, we usually observe two frequency jumps for samples with relatively large exchange bias angle 关EB ⲏ 30°, e.g., Figs. 6共e兲 and 6共f兲兴 and one jump for samples relatively small exchange bias angle 关EB ⱗ 30°, e.g., Figs. 6共a兲 and 6共b兲兴. However, it is not always the case that
024418-8
PHYSICAL REVIEW B 76, 024418 共2007兲
LARGE-AMPLITUDE COHERENT SPIN WAVES EXCITED…
(a)
(b)
(c)
(d)
(e)
(f)
FIG. 6. 共Color online兲 共a兲, 共c兲, 共e兲 Examples of hysteresis loops of resistance versus field of different nominally identical samples. 共b兲, 共d兲, 共f兲 Frequency of persistent magnetization dynamics as a function of current for the samples in 共a兲, 共c兲, 共e兲 measured at H = 600 Oe.
samples with similar direction of exchange bias as determined by the Stoner-Wohlfarth fitting procedure 关e.g., Figs. 1共c兲 and 6共c兲兴 show identical dependence of frequency on current 关Figs. 2共c兲 and 6共d兲兴. We attribute these differences to sample shape imperfections. IV. DISCUSSION
Having developed an understanding about how the various parameters influence the simulated dynamics, we can proceed with the analysis of magnetization oscillation spectra and a comparison with the experimentally observed magnetoresistance power spectra. In Fig. 4共a兲 we display the dependence of the simulated spectral maximum frequencies on the spin-torque amplitude aJ, in the presence of thermal fluctuations. These simulations take into account all the physical factors that are generally included in a state-of-the-art micromagnetic model. Spectral amplitudes of magnetoresistance oscillations are displayed in Figs. 4共b兲–4共d兲. The spectra can be divided into the following three groups: 共i兲 from the oscillation onset to the first frequency jump 关group A, Fig. 4共b兲兴, 共ii兲 from the first to the second frequency jump 关group B, Fig. 4共c兲兴, and after the second frequency jump 关group C, Fig. 4共d兲兴. For all groups, the frequency of the spectral maximum decreases monotonically with the spin-torque magnitude aJ. The dependencies of the linewidth and the integrated spectral power 关see Fig. 7共b兲兴 on the spin-torque magnitude require special discussion. For the first group—spectra from the oscillation onset to the first frequency jump—the linewidth for small aJ is relatively large 共⬃100 MHz兲 due to a relatively large influence
FIG. 7. 共Color online兲 共a兲 Comparison of the experimentally measured 共solid triangles兲 and simulated 共open circles兲 dependence of the frequency of oscillations on current. 共b兲 Experimentally measured 共solid triangles兲 and simulated 共open circles兲 integrated rmsamplitude spectral density Sint as a function of current.
of thermal fluctuations on small-amplitude motion of the magnetization.61 The oscillation amplitude grows rapidly with increasing current 共compare spectra for aJ = 0.30, 0.31, and 0.32兲 and the linewidth strongly decreases 共to ⬃20 MHz for aJ = 0.32兲, which is due to an increasing contribution from the spin-torque-driven dynamics resulting in the effective suppression of the influence of thermal fluctuations and thus in the decrease of the linewidth.61,68 When the current is increased further and approaches its value for the first jump, the contribution of the second nonlinear oscillation mode 共which will dominate the spectrum after the first frequency jump兲 becomes visible, leading to line broadening and a decrease of the maximal spectral amplitude 共see spectra for aJ = 0.34, 0.36兲. After the first jump, the amplitude of magnetization precession becomes large 共Sint ⬎ 0.3兲 and the relative influence of thermal fluctuations on the motion of magnetization becomes small. For this reason, the linewidth for most spectra of the second group 共except for those close to the second frequency jump兲 is extremely small. In fact, it is below the resolution limit of our simulations 共⌬f min ⬇ 10 MHz兲, thus being in a good agreement with experimental observations. When approaching the current value of the second frequency jump, the line width starts to increase again 共and the maximal spectral power decreases兲 due to the influence of the next nonlinear mode.
024418-9
PHYSICAL REVIEW B 76, 024418 共2007兲
KRIVOROTOV et al.
For the last spectral group, the linewidth and the maximal value of the spectral power exhibit the same nonmonotonic behavior. However, the line broadening for the large current values 共aJ ⬎ 0.90兲 in this region is not due to the next incipient frequency jump but due to the onset of spatially incoherent magnetization dynamics 共note that the maximal experimentally used current value Imax = 10 mA corresponds to aJ ⬇ 1.0兲. The line broadening for I ⬎ 9.0 mA is also observed experimentally and is clearly visible in Fig. 2共c兲. However, for values of aJ greater than the second frequency jump, the width of the simulated spectral peaks is substantially larger than the width of spectral lines measured for corresponding current values 共computed as I = aJ / 兲. Another important difference between experiment and the simulation results is that the narrowest spectral lines found in the simulations exist between the first and the second frequency jumps, while the narrowest lines observed experimentally occur after the second frequency jump. Before proceeding to a direct comparison with the experimental oscillation frequencies and amplitudes, we note an important difference between the magnetization dynamics of the Py elliptical nanomagnet simulated in this paper and that of the Co elliptical nanoelement studied in detail previously in Ref. 63. For the Co element in Ref. 63 the coherence of the magnetization oscillations was lost already for currents very close to the onset of the steady-state oscillations, followed by a transition to a completely chaotic regime.69 In contrast to this behavior, magnetization dynamics of the Py element studied here remains nearly coherent up to current values several times larger than the critical current. This difference cannot be attributed to much lower temperatures for which the experiment discussed here has been performed 共compared to room temperatures used in Ref. 10兲, because the transition to the chaotic regime slightly above acr was observed in Ref. 63 already for simulations performed at T = 0. The difference can also not be due to a slightly higher element thickness used here 共hPy = 4 nm compared to hCo = 3 nm in Ref. 63兲, because the much higher exchange constant of Co 共ACo = 3 ⫻ 10−6 erg/ cm; see Ref. 63 for details兲 when compared with the Py exchange 共APy = 1.3⫻ 10−6 erg/ cm兲 should at least compensate this slightly larger thickness of the Py nanoelement. We argue that this important discrepancy in the behavior of the two quite similar systems studied here and in Ref. 63 is due to the very different character of the nonlinear magnetization oscillation modes of these nanoelements. Whereas in Ref. 63 several oscillation modes with a quite complicated localization patterns arose and coexisted when the oscillation amplitude increased 共see spatial maps in Figs. 1, 3, and 4 in Ref. 63兲, in this work we have found that for each given current value there is a single nonlinear eigenmode where the oscillating spins are confined in a localized area of the nanomagnet without any node lines between these oscillating spins. It seems plausible that the transition to a quasichaotic behavior from a single mode would be inhibited compared to the case of several coexisting modes with different spatial profiles. We believe that the physical reason for excitation of a single mode in the case of substantially noncollinear magnetizations of the pinned and the free layers is that spin torque is nearly spatially uniform. Indeed, in the case of
nominally collinear magnetizations, the direction and magnitude of spin torque exerted on the free layer exhibits strong spatial variations due to spatially nonuniform magnetization direction predicted by micromagnetics. This results in local magnetizations of the free and pinned layers making small negative angles with respect to each other in some parts of the sample and small positive angles in other parts of the sample. Since spin torque is proportional to the small angle between magnetizations of the free and pinned layers, the case of nominally collinear magnetizations gives rise to strongly spatially nonuniform spin torque. In the case of noncollinear magnetizations, small variations of the magnetization direction over the sample area result in small deviations of spin-torque direction and magnitude from their average values. A spatially nonuniform pattern of spin torque is more likely to couple to multiple oscillation modes of the nanomagnet. In the case of nearly constant uniform torque, the coupling to the longest-wavelength mode is expected to be the strongest. Figure 7 presents a direct comparison between experimental data and results of LLGS simulations. First, we show in Fig. 7共a兲 the current dependence of the magnetization oscillation frequency as measured experimentally 共solid triangles兲 and as obtained from micromagnetic simulations 共open circles兲. For plotting the simulation data as f共I兲 we have used the conversion from the spin torque amplitude aJ to the current strength I in the form I = aJ / with the conversion factor ⬇ 0.1 computed as explained above. The simulations reproduce the current dependence of the oscillation frequency fairly well, except for the position of the second frequency jump, which occurs in simulations at a current about 20% higher than in the experiment. However, taking into account that a nanomagnet with perfect edges was simulated and that the simulations did not contain any adjustable parameters 共except the conversion factor 兲 the agreement between simulations and experiment can be considered as very satisfactory, as far as the oscillation frequency is concerned. The fact that micromagnetic LLGS simulations successfully reproduce the frequency jumps in the current dependence of the frequency of oscillations is a significant success of the LLGS micromagnetic approach. Since these jumps result from transitions between modes with different degrees of spatial localization, LLGS simulations in the macrospin approximation would not be able to describe such transitions. This shows that while LLGS simulations in the macrospin approximation10,56 are useful for an initial qualitative understanding of many properties of current-driven magnetization dynamics such as the nonlinear shift of the oscillation frequency with current and the existence of different types of nonlinear oscillation modes 共e.g., in-plane and out-of-plane precession modes10,66兲, a quantitative description of persistent current-driven dynamics requires a micromagnetic approach. Despite the good agreement between experiment and simulations for the dependence of the oscillation frequency on current, our micromagnetic simulations could not closely reproduce the corresponding dependence of the oscillation amplitude on current. Figure 7共b兲 shows the experimentally measured 共solid triangles兲 and simulated 共open circles兲 integrated signal amplitudes Sint as functions of the bias current.
024418-10
PHYSICAL REVIEW B 76, 024418 共2007兲
LARGE-AMPLITUDE COHERENT SPIN WAVES EXCITED…
The general trend of the measured oscillation amplitude is to increase gradually with increasing current, together with a series of dips at currents corresponding to the frequency jumps shown in Fig. 2共c兲. In contrast, the LLGS simulations predict a rapid increase of the oscillation amplitude just above the critical current for the onset of oscillations, followed by a slow gradual decrease. Some minor anomalies on the simulated Sint共I兲 around the frequency jumps can be seen; however, they are far less pronounced than the corresponding experimentally observed dips. Taking into account the good agreement between simulations and experiment for the oscillation frequency, the discrepancy for the amplitude of oscillations is very surprising and requires a detailed analysis. The failure of our LLGS simulations to predict the correct dependence of the oscillation amplitude on current indicates that the standard micromagnetic LLGS approach for spin-torque-driven excitations in nano-magnets requires modifications. Below we propose some possible routes towards improvement of the theoretical description of spin-torque-driven excitations in nanomagnets. One possible way of solving this problem would be to introduce a nonlinear dissipation 关a dependence of the dissipation parameter on the rate of the magnetization change dm / dt in the form = 0共1 + q1共dm / dt兲2 + ¯ 兲兴 as suggested in Ref. 40. In making such an attempt, one should keep in mind that a too strong nonlinearity 共large values of the nonlinear coefficient q1兲 would destroy the good agreement between simulated and measured oscillation frequencies, especially for the initial part of the f共aJ兲 dependence where the transition between linear and nonlinear oscillation regimes is observed. However, a moderate nonlinearity could weakly affect the oscillation frequency for small to moderate oscillation amplitudes 共small aJ兲, while improving the coherence of the magnetization oscillations for large currents. 共If the dissipation coefficient increases with increasing dm / dt, then it should strongly suppress the short-wavelength excitations that lead to incoherent magnetization oscillations.兲 In this way, one would obtain higher oscillation powers and narrower linewidths for larger currents, thus improving the agreement between theory and experiment. Clearly, this subject requires further investigation. Another possible way of reconciling theory and experiment for the current dependence of both frequency and amplitude of the excited modes would be the generation of spinwave modes that are more spatially nonuniform than those shown in Fig. 3共a兲. Indeed, if only a part of magnetization of the nanomagnet moves with large amplitude 共e.g., edge modes兲, both a significant nonlinear shift of frequency and a relatively small average measured amplitude will be observed. Furthermore, the growth of the average amplitude of such nonuniform spin-wave modes is likely to proceed via a gradual spatial growth of the oscillating domain, which should give rise to a gradual increase of the measured amplitude and result in a dependence of the amplitude on current similar to the experimentally observed dependence shown in Fig. 7共b兲. A possible mechanism leading to excitation of strongly spatially nonuniform modes is the instability of magnetization arising from lateral spin transport in spinvalve structures.70–72 A theoretical test of this scenario re-
quires the development of a micromagnetic code that explicitly treats magnetization dynamics coupled to spatially nonuniform spin-dependent electrical transport, which is beyond the scope of this work. Softening of the spin-wave spectrum by spin-polarized current73,74 could also be an important factor to be taken into account for reconciling the theory of current-driven excitations with the experimental results presented in this paper. V. CONCLUSIONS
In conclusion, we have measured the spectral properties of current-driven excitations in nanoscale spin valves with noncollinear magnetizations of the free and pinned ferromagnetic layers. We find that spin-polarized current in these devices excites a few coherent large-amplitude nonlinear modes of magnetization oscillation in the free layer. Different modes are excited in different current intervals. We find that the amplitude and the coherence of the current-driven excitations decrease in the current intervals where transitions between these modes take place. We simulate the response of magnetization to spin-polarized current in our samples by employing LLG micromagnetic simulations with a Slonczewski spin-torque term.46 These LLGS simulations capture a number of features of the experimental data: 共i兲 the decrease of frequency of the excited oscillation modes with increasing current, 共ii兲 downward jumps of the frequency of excitations with increasing current resulting from transitions between different oscillation modes, and 共iii兲 the high degree of coherence 共narrow spectral linewidth兲 of the excited modes. However, the LLGS simulations give qualitatively incorrect predictions for the amplitude of the excited modes as a function of current. Simulations predict rapid growth of the oscillation amplitude above the threshold current for the onset of spin-wave excitations, followed by a slow decrease of the amplitude. This is in sharp contrast to the more gradual increase of the oscillation amplitude with current observed in our experiment. Our results demonstrate that additional factors possibly including nonlinear damping and/or lateral spin transport need to be taken into account for a quantitative description of large-amplitude magnetization dynamics driven by spin-polarized current in magnetic nanostructures. ACKNOWLEDGMENTS
The authors thank J. Miltat, D. Mills, and A. Slavin for many useful discussions. This research was supported by the Deutsche Forschungsgemeinschaft 共DFG Grant No. BE 2464/4-1兲, the Office of Naval Research, and the National Science Foundation’s Nanoscale Science and Engineering Centers program through the Cornell Center for Nanoscale Systems. We also acknowledge NSF support through use of the Cornell Nanoscale Facility node of the National Nanofabrication Infrastructure Network and the use of the facilities of the Cornell Center for Materials Research. APPENDIX: NUMERICAL SIMULATION METHODOLOGY
In our full-scale micromagnetic simulations magnetization dynamics are simulated by solving the stochastic LLG
024418-11
PHYSICAL REVIEW B 76, 024418 共2007兲
KRIVOROTOV et al.
equation of motion for the magnetization Mi of each discretization cell in the following form: dMi fl = − ␥关Mi ⫻ 共Hdet i + Hi 兲兴 dt −␥
fl 关Mi ⫻ 关Mi ⫻ 共Hdet i + Hi 兲兴兴. MS
共A1兲
Here ␥ = ␥0 / 共1 + 2兲, where ␥0 共⬎0兲 is the absolute value of the gyromagnetic ratio and is the reduced dissipation constant 共it is equal to the constant ␣ in the LLG equation writ˙ = −␥ 关M ⫻ H兴 + 共␣ / M 兲关M ⫻ M兴兲. The deten in the form M 0 S terministic effective field Hdet acting on the magnetization of i the ith cell includes all standard micromagnetic contributions 共external, anisotropy, exchange, and magnetodipolar interaction fields兲. In addition, this deterministic field includes the spintorque effect in the following way. The spin torque in the Slonczewski formalism is taken into account adding the term G = −f J共兲关M ⫻ 关M ⫻ p兴兴 to the equation of motion in the Gilbert form 共see, e.g., Ref. 56兲
␣ dM ˙兴 = − ␥0关M ⫻ Heff兴 + 关M ⫻ M dt MS − ␥0 f J共兲关M ⫻ 关M ⫻ p兴兴,
共A2兲
where the dimensionless spin-torque amplitude f J depends on the angle between the magnetization M and the unit vector p of the polarization direction of the electron magnetic moments 共in the spin-polarized current兲. From the computational point of view, this additional torque can be put eff = Heff + f J关M ⫻ p兴, after which into the effective field as HST Eq. 共A2兲 can be converted to the numerically more convenient form 共A1兲 in a standard way. We use the following asymmetric angular dependence of the amplitude f J共兲 共Refs. 46 and 62兲: f J共 兲 = a J
2⌳2 . 共⌳2 + 1兲 + 共⌳2 − 1兲cos
共A3兲
Here aJ gives the 共constant兲 value of the spin-torque amplitude for the symmetric torque 共⌳ = 1兲; the asymmetry parameter ⌳ can in principle be computed when the device configuration and various transport coefficients are known and is related to the GMR asymmetry parameter in Eq. 共1兲 via ⌳2 = 1 + . Expression 共A3兲 is strictly valid only for symmetrical spin valves46,62 with identical ferromagnetic layers and identical top and bottom leads. The expression for f J共兲 in an asymmetric device is more complex62 and involves effective resistances of the ferromagnetic layers and leads. However, we use a simplified expression 共A3兲 for f J共兲 in our simulations for three reasons. First, we do not expect the spin-torque asymmetry of our device 共with respect to the above-mentioned effective resistances兲 to be large because the thicknesses of two ferromagnetic layers of the spin valve are identical and the thicknesses of nonmagnetic leads are not very different. Second, the spin diffusion length of Py 关⬃5 nm 共Ref. 75兲兴 is similar to the thickness of Py layers in our spin-valve structure, which substantially decreases the
influence of the transport properties of the leads on spin torque.62 Third, Eq. 共A3兲 can be considered as the simplest form 共apart from the form with f J = const兲 for studying the effect of the asymmetry of the spin-torque angular dependence on the magnetization dynamics. The more complex expression derived in Ref. 62 can be investigated after the effect of the simplest form of spin-torque asymmetry given by Eq. 共A3兲 is understood. The random fluctuation field Hfli in Eq. 共A1兲 represents the influence of thermal fluctuations and has standard ␦-functional spatial and temporal correlation properties: 具Hfl,i典 = 0,
具Hfl,i共0兲Hfl,j共t兲典 = 2D ␦共t兲␦ij␦
共A4兲
共i , j are the discretization cell indices; , = x , y , z兲, with the noise power D proportional to the system temperature D = / 共1 + 2兲共kT / ␥兲; here, denotes the magnetic moment magnitude for a single discretization cell. The justification for using ␦-correlated random noise for a finite-element version of an initially continuous system with interactions can be found, e.g., in Ref. 76. The remaining simulation methodology is similar to that described in Ref. 63. We simulate spin-torque-driven excitations in the free ferromagnetic layer only. We neglect magnetostatic and RKKY interactions between the free and pinned 关antiferromagnetic- 共AF-兲 coupled兴 Py layers. This approximation is justified because, first, the RKKY exchange coupling via the thick 共hsp = 8 nm兲 Cu spacer is negligibly small and, second, the dipolar field acting on the free layer from the fixed one is on average 共⬃80 Oe兲 much smaller than the external field. The free layer 共130⫻ 60⫻ 4 nm3 ellipse兲 is discretized into 50⫻ 24⫻ 1 cells; we checked that further refinement of the grid did not lead to any significant changes in the results. Magnetic parameters of the free Py layer used in simulations are: saturation magnetization M S = 650 emu/ cm3 共measured by SQUID magnetometry as explained in Sec. II兲; exchange constant A = 1.3⫻ 10−6 erg/ cm 共standard value for Py兲; the random magnetocrystalline anisotropy was neglected due to its low value 共Kcub = 5 ⫻ 103 erg/ cm3兲 for Py. The dissipation parameter is set to = 0.025 共see also Sec. II兲. As was shown above, the influence of the Oersted field HOe induced by the current flowing through the spin valve may be very important. The calculation of this field is also a highly nontrivial issue. In principle, its precise evaluation requires the exact knowledge of the three-dimensional current distribution in the device itself and especially in adjacent electrical contact layers, which is normally not available. For this reason the Oersted field is usually computed assuming that the current is distributed homogeneously across the nanopillar cross section. Further, one of the following approximations is used: 共i兲 one assumes that HOe is created by the infinitely long wire with the cross section corresponding to that of the nanopillar 共in our case the ellipse with la ⫻ lb = 130⫻ 60 nm2兲 or 共ii兲 the contribution to the Oersted field from the current inside the nanopillar itself only is included 共i.e., HOe is created by the piece of the wire with the length equal to the nanopillar height htot兲. Both approximations deliver the same result for nanopillars with the height much
024418-12
PHYSICAL REVIEW B 76, 024418 共2007兲
LARGE-AMPLITUDE COHERENT SPIN WAVES EXCITED…
larger than their characteristic cross-section size 关htot Ⰷ max共la , lb兲兴, which is, however, very rarely the case for experiments performed up to now in the nanopillar geometry. In particular, in our situation the opposite inequality is true 共Ia ⬎ htot兲. Taking into account that the first approximation is also reasonably accurate for the system where the distribution of currents in the nanopillar and adjacent leads is axially symmetric, and that the geometry of the electric contacts in our device is also highly symmetric, we have chosen the first method to calculate HOe. However, we point out once more that due to the importance of the influence of the Oersted field more precise methods for its calculation are highly desirable. Magnetization dynamics was simulated by integrating Eq. 共A1兲 with the spin-torque term included using the addition-
J. C. Slonczewski, J. Magn. Magn. Mater. 159, L1 共1996兲. Berger, Phys. Rev. B 54, 9353 共1996兲. 3 M. Tsoi, A. G. M. Jansen, J. Bass, W. C. Chiang, M. Seck, V. Tsoi, and P. Wyder, Phys. Rev. Lett. 80, 4281 共1998兲; 81, 493共E兲 共1998兲. 4 E. B. Myers, D. C. Ralph, J. A. Katine, R. N. Louie, and R. A. Buhrman, Science 285, 867 共1999兲. 5 J.-E. Wegrowe, D. Kelly, Y. Jaccard, P. Guittienne, and J.-P. Ansermet, Europhys. Lett. 45, 626 共1999兲. 6 J. Z. Sun, J. Magn. Magn. Mater. 202, 157 共1999兲. 7 J. A. Katine, F. J. Albert, R. A. Buhrman, E. B. Myers, and D. C. Ralph, Phys. Rev. Lett. 84, 3149 共2000兲. 8 J. Grollier, V. Cros, A. Hamzic, J. M. George, H. Jaffres, A. Fert, G. Faini, J. Ben Youssef, and H. Legall, Appl. Phys. Lett. 78, 3663 共2001兲. 9 M. Tsoi, A. G. M. Jansen, J. Bass, W. C. Chiang, V. Tsoi, and P. Wyder, Nature 共London兲 406, 46 共2000兲. 10 S. I. Kiselev, J. C. Sankey, I. N. Krivorotov, N. C. Emley, R. J. Schoelkopf, R. A. Buhrman, and D. C. Ralph, Nature 共London兲 425, 380 共2003兲. 11 W. H. Rippard, M. R. Pufall, S. Kaka, S. E. Russek, and T. J. Silva, Phys. Rev. Lett. 92, 027201 共2004兲. 12 M. Covington, M. Al Haj Darwish, Y. Ding, N. J. Gokemeijer, and M. A. Seigler, Phys. Rev. B 69, 184406 共2004兲. 13 B. Oezyilmaz, A. D. Kent, D. Monsma, J. Z. Sun, M. J. Rooks, and R. H. Koch, Phys. Rev. Lett. 91, 067203 共2003兲. 14 I. N. Krivorotov, N. C. Emley, J. C. Sankey, S. I. Kiselev, D. C. Ralph, and R. A. Buhrman, Science 307, 228 共2005兲. 15 P. W. Anderson and H. Suhl, Phys. Rev. 100, 1788 共1955兲. 16 H. Suhl, J. Phys. Chem. Solids 1, 209 共1957兲. 17 K. Y. Guslienko, R. W. Chantrell, and A. N. Slavin, Phys. Rev. B 68, 024422 共2003兲. 18 J. Jorzick, S. O. Demokritov, B. Hillebrands, M. Bailleul, C. Fermon, K. Y. Guslienko, A. N. Slavin, D. V. Berkov, and N. L. Gorn, Phys. Rev. Lett. 88, 047204 共2002兲. 19 R. Arias, P. Chu, and D. L. Mills, Phys. Rev. B 71, 224410 共2005兲. 20 R. D. McMichael and M. D. Stiles, J. Appl. Phys. 97, 10J901 共2005兲. 21 J. P. Park, P. Eames, D. M. Engebretson, J. Berezovsky, and P. A. 1
2 L.
ally optimized Bulirsch-Stoer algorithm77 with the adaptive step-size control. 共The adaptive step-size control is especially important when the magnetization state significantly deviates from a homogeneous one.兲 For each current value 共each value of aJ in our formalism兲 the dependence of the magnetization on time for every discretization cell was saved for the physical time interval ⌬t = 400 ns. The spectral analysis of these magnetization “trajectories” was performed using either 共i兲 the Lomb algorithm 共as described in Ref. 63兲 especially designed for nonevenly spaced sequences of time moments as provided by the adaptive integration method or 共ii兲 interpolation of the “raw” results onto an evenly spaced temporal grid and usage of the standard fast Fourier transform routines. Results of both methods turned out to be equivalent within the statistical errors.
Crowell, Phys. Rev. B 67, 020403共R兲 共2003兲. B. Zhu, Z. G. Liu, V. Metlushko, P. Grutter, and M. R. Freeman, Phys. Rev. B 71, 180408共R兲 共2005兲. 23 M. Kostylev, J. G. Hu, and R. L. Stamps, Appl. Phys. Lett. 90, 012507 共2007兲. 24 G. De Loubens, V. V. Naletov, Olivier Klein, J. Ben Youssef, F. Boust, and N. Vukadinovic, Phys. Rev. Lett. 98, 127601 共2007兲. 25 P. E. Wigen, M. L. Roukes, and P. C. Hammel, Top. Appl. Phys. 101, 105 共2006兲. 26 T. Mewes, J. Kim, D. V. Pelekhov, G. N. Kakazei, P. E. Wigen, S. Batra, and P. C. Hammel, Phys. Rev. B 74, 144424 共2006兲. 27 M. Grimsditch, F. Y. Fradin, Y. Ji, A. Hoffmann, R. E. Camley, V. Metlushko, and V. Novosad, Phys. Rev. Lett. 96, 047401 共2006兲. 28 T. M. Crawford, M. Covington, and G. J. Parker, Phys. Rev. B 67, 024411 共2003兲. 29 S. M. Rezende, F. M. de Aguiar, and A. Azevedo, Phys. Rev. Lett. 94, 037202 共2005兲. 30 A. N. Slavin and P. Kabos, IEEE Trans. Magn. 41, 1264 共2005兲. 31 A. Yu. Dobin and R. H. Victora, Phys. Rev. Lett. 90, 167203 共2003兲. 32 A. Misra and R. H. Victora, Phys. Rev. B 73, 172414 共2006兲. 33 S. Y. An, P. Krivosik, M. A. Kraemer, H. M. Olson, A. V. Nazarov, and C. E. Patton, J. Appl. Phys. 96, 1572 共2004兲. 34 T. J. Silva, C. S. Lee, T. M. Crawford, and C. T. Rogers, J. Appl. Phys. 85, 7849 共1999兲. 35 S. G. Reidy, L. Cheng, and W. E. Bailey, Appl. Phys. Lett. 82, 1254 共2003兲. 36 J. Lindner, K. Lenz, E. Kosubek, K. Baberschke, D. Spoddig, R. Meckenstock, J. Pelzl, Z. Frait, and D. L. Mills, Phys. Rev. B 68, 060102共R兲 共2003兲. 37 T. Gerrits, M. L. Schneider, A. B. Kos, and T. J. Silva, Phys. Rev. B 73, 094454 共2006兲. 38 Y. Tserkovnyak, A. Brataas, and G. E. W. Bauer, Phys. Rev. Lett. 88, 117601 共2002兲. 39 B. Heinrich, Y. Tserkovnyak, G. Woltersdorf, A. Brataas, R. Urban, and G. E. W. Bauer, Phys. Rev. Lett. 90, 187601 共2003兲. 40 V. Tiberkevich and A. Slavin, Phys. Rev. B 75, 014440 共2007兲. 41 P. Lubitz, M. Rubinstein, J. J. Krebs, and S. F. Cheng, J. Appl. Phys. 89, 6901 共2001兲. 22 X.
024418-13
PHYSICAL REVIEW B 76, 024418 共2007兲
KRIVOROTOV et al. 42 R.
L. Compton, M. J. Pechan, S. Maat, and E. E. Fullerton, Phys. Rev. B 66, 054411 共2002兲. 43 M. Fraune, U. Rudiger, G. Guntherodt, S. Cardoso, and P. Freitas, Appl. Phys. Lett. 77, 3815 共2000兲. 44 J. Eisenmenger, Z. P. Li, W. A. A. Macedo, and I. K. Schuller, Phys. Rev. Lett. 94, 057203 共2005兲. 45 I. V. Roshchin, O. Petracic, R. Morales, Z. P. Li, X. Batlle, and I. K. Schuller, Europhys. Lett. 71, 297 共2005兲. 46 J. C. Slonczewski, J. Magn. Magn. Mater. 247, 324 共2002兲. 47 S. Urazhdin, R. Loloee, and W. P. Pratt, Phys. Rev. B 71, 100401共R兲 共2005兲. 48 M. J. Donahue and D. G. Porter 共unpublished兲. 49 C. G. Lee, J. G. Jung, V. S. Gornakov, R. D. McMichael, A. Chen, and W. F. Egelhoff, J. Magn. Magn. Mater. 272, 1887 共2004兲. 50 D. V. Berkov and N. Gorn, MICROMAGUS, package for micromagnetic simulations, http://www.micromagus.de 51 G. D. Fuchs, J. C. Sankey, V. S. Pribiag, L. Qian, P. M. Braganca, A. G. F. Garcia, E. M. Ryan, Zhi-Pan Li, O. Ozatay, D. C. Ralph, and R. A. Buhrman, arXiv:cond-mat/0703577 共unpublished兲. 52 N. C. Emley, I. N. Krivorotov, O. Ozatay, A. G. F. Garcia, J. C. Sankey, D. C. Ralph, and R. A. Buhrman, Phys. Rev. Lett. 96, 247204 共2006兲. 53 I. N. Krivorotov, N. C. Emley, A. G. F. Garcia, J. C. Sankey, S. I. Kiselev, D. C. Ralph, and R. A. Buhrman, Phys. Rev. Lett. 93, 166603 共2004兲. 54 A. Fukushima, H. Kubota, A. Yamamoto, Y. Suzuki, and S. Yuasa, J. Appl. Phys. 99, 08H706 共2006兲. 55 D. M. Pozar, Microwave Engineering 共Wiley, New York, 1998兲, p. 89. 56 J. Z. Sun, Phys. Rev. B 62, 570 共2000兲. 57 G. Bertotti, C. Serpico, I. D. Mayergoyz, A. Magni, M. d’Aquino, and R. Bonin, Phys. Rev. Lett. 94, 127206 共2005兲. 58 J. V. Kim, Phys. Rev. B 73, 174412 共2006兲. 59 J. Grollier, V. Cros, and A. Fert, Phys. Rev. B 73, 060409共R兲 共2006兲. 60 K. Mizushima, K. Kudo, and R. Sato, J. Magn. Magn. Mater. 316, E960 共2007兲.
61 J.
C. Sankey, I. N. Krivorotov, S. I. Kiselev, P. M. Braganca, N. C. Emley, R. A. Buhrman, and D. C. Ralph, Phys. Rev. B 72, 224427 共2005兲. 62 J. Xiao, A. Zangwill, and M. D. Stiles, Phys. Rev. B 70, 172405 共2004兲. 63 D. V. Berkov and N. L. Gorn, Phys. Rev. B 72, 094401 共2005兲. 64 M. D. Stiles and J. Miltat, in Spin Dynamics in Confined Magnetic Structures III, Springer Series Topics in Applied Physics, Vol. 101 共Springer-Verlag, Berlin, 2006兲. 65 A. N. Slavin and V. S. Tiberkevich, Phys. Rev. B 72, 094428 共2005兲. 66 S. I. Kiselev, J. C. Sankey, I. N. Krivorotov, N. C. Emley, A. G. F. Garcia, R. A. Buhrman, and D. C. Ralph, Phys. Rev. B 72, 064430 共2005兲. 67 F. B. Mancoff, R. W. Dave, N. D. Rizzo, T. C. Eschrich, B. N. Engel, and S. Tehrani, Appl. Phys. Lett. 83, 1596 共2003兲. 68 Q. Mistral, J. V. Kim, T. Devolder, P. Crozat, C. Chappert, J. A. Katine, M. J. Carey, and K. Ito, Appl. Phys. Lett. 88, 192507 共2006兲. 69 K. J. Lee, A. Deac, O. Redon, J. P. Nozières, and B. Dieny, Nat. Mater. 3, 877 共2004兲. 70 M. L. Polianski and P. W. Brouwer, Phys. Rev. Lett. 92, 026602 共2004兲. 71 M. D. Stiles, J. Xiao, and A. Zangwill, Phys. Rev. B 69, 054408 共2004兲. 72 B. Ozyilmaz, A. D. Kent, J. Z. Sun, M. J. Rooks, and R. H. Koch, Phys. Rev. Lett. 93, 176604 共2004兲. 73 Y. B. Bazaliy, B. A. Jones, and S. C. Zhang, Phys. Rev. B 57, R3213 共1998兲. 74 J. Fernandez-Rossier, M. Braun, A. S. Nunez, and A. H. MacDonald, Phys. Rev. B 69, 174412 共2004兲. 75 J. Bass and W. P. Pratt, J. Magn. Magn. Mater. 200, 274 共1999兲. 76 C. Bayer, J. Jorzick, B. Hillebrands, S. O. Demokritov, R. Kouba, R. Bozinoski, A. N. Slavin, K. Guslienko, D. V. Berkov, N. L. Gorn, and M. P. Kostylev, Phys. Rev. B 72, 064427 共2005兲. 77 W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in Fortran: the Art of Scientific Computing 共Cambridge University Press, Cambridge, England, 1992兲.
024418-14