Reservoir Computing on Spin-Torque Oscillator Array

TL;DR

This paper demonstrates that spin-torque oscillator arrays can be effectively used for reservoir computing, with performance influenced by synchronization states and surpassing traditional echo-state networks in certain configurations.

Contribution

It introduces the use of STO arrays for reservoir computing and shows their potential advantages over conventional neural network models.

Findings

Reservoir computing is feasible with synchronized STO oscillations.

Performance improves with more STOs and peaks at the boundary between synchronization and disorder.

STO arrays outperform similar-sized echo-state networks in reservoir computing tasks.

Abstract

We numerically study reservoir computing on a spin-torque oscillator (STO) array, describing the magnetization dynamics of the STO array by a nonlinear oscillator model. The STOs exhibit synchronized oscillation due to coupling by magnetic dipolar fields. We show that reservoir computing can be performed using the synchronized oscillation state. The performance can be improved by increasing the number of STOs. The performance becomes highest at the boundary between the synchronized and disordered states. Using an STO array, we can achieve higher performance than that of an echo-state network with similar number of units. This result indicates that STO arrays are promising for hardware implementation of reservoir computing.

Equations

$$\frac{\mathrm{d}{c}_i}{\mathrm{d}t}=-\mathrm{i}\delta {\omega }_i(p_i)c_i-{\Gamma }_i(p_i)c_i-\mathrm{i}\sum _{j=1,j\ne i}^N{\Omega }_{ij}{c}_j,$$

$$\delta {\omega }_i(p_i)$$

$${\Gamma }_i(p_i)$$

$$\rho$$

$$\bar{p}$$

$$\frac{\partial {\mathbf{m}}_i}{\partial t}=-\gamma {\mathbf{m}}_i\times {\mathbf{H}}_i^{\mathrm{eff}}+\alpha {\mathbf{m}}_i\times \frac{\partial {\mathbf{m}}_i}{\partial t}+{\mathbf{T}}_i^S,$$

$$\frac{\mathrm{d}{\tilde{c}}_i}{\mathrm{d}t}=-\mathrm{i}\gamma H_{iz}^{\mathrm{eff}}{\tilde{c}}_i+\mathrm{i}\gamma m_{iz}\left(H_{ix}^{\mathrm{eff}}-\mathrm{i}{H}_{iy}^{\mathrm{eff}}\right)-{\Gamma }_i{\tilde{c}}_i,$$

$$\mathbf{H}=3\frac{\left(\mathbf{\mu }\cdot \mathbf{r}\right)\mathbf{r}}{r^5}-\frac{\mathbf{\mu }}{r^3},$$

$$H_x-\mathrm{i}{H}_y=\frac{M_{s}V}{2r^5}$$

$$\times \left[\left(x^2+y^2-2z^2\right)\tilde{c}+3(x-\mathrm{i}y{)}^2{\tilde{c}}^{\ast }+6(x-\mathrm{i}y)z{m}_z\right].$$

$$H_x-\mathrm{i}{H}_y=\frac{M_{s}V}{2r^3}\left[\tilde{c}+3\frac{(x-\mathrm{i}y{)}^2}{r^2}{\tilde{c}}^{\ast }\right].$$

$$\frac{\mathrm{d}{\tilde{c}}_i}{\mathrm{d}t}=-\mathrm{i}{\tilde{\omega }}_i{\tilde{c}}_i-\mathrm{i}\sum _{j=1,j\ne i}^N\left({\Omega }_{ij}{\tilde{c}}_j+{\Xi }_{ij}{\tilde{c}}_j^{\ast }\right)-{\Gamma }_i{\tilde{c}}_i,$$

$$\frac{\mathrm{d}{c}_i}{\mathrm{d}t}=-\mathrm{i}\delta {\omega }_i{c}_i-\mathrm{i}\sum _j\left({\Omega }_{ij}{c}_j+{\Xi }_{ij}{c}_j^{\ast }{\mathrm{e}}^{2\mathrm{i}\omega t}\right)-{\Gamma }_i{c}_i.$$

$$\frac{\mathrm{d}{c}_i}{\mathrm{d}t}=-\mathrm{i}\delta {\omega }_i{c}_i-{\Gamma }_i{c}_i-\mathrm{i}\sum _{j=1,j\ne i}^N{\Omega }_{ij}{c}_j,$$

$${\bar{c}}_i(t)=\frac{\omega }{\pi }{\int }_0^{\pi /\omega }\mathrm{d}{t}^{\prime }{c}_i(t+t^{\prime }).$$

$$\frac{\partial \mathbf{m}}{\partial t}$$

$${\Gamma }_i$$

$$=$$

$${\Gamma }_p=p_{i0}\frac{\partial {\Gamma }_i}{\partial p_i}.$$

$$\frac{\partial {\Gamma }_i}{\partial I_i}\left(I_i-I_{i0}\right)=-{\xi }_i(t).$$

Full text

I Introduction

Reservoir computing Jaeger2004 ; Maass2002 ; Verstraeten2007 , a framework for machine learning, has attracted much attention for realizing high-performance real-time computing for applications such as time-series modeling and prediction, signal processing, and pattern classification. Reservoir computing utilizes the complex dynamics of interacting nonlinear units. In reservoir computing, only output weights are learned in an external computer. This feature allows hardware implementations of reservoir computing, as well as reduction of the computational costs of the learning process. Hardware implementations have been proposed in several physical systems, such as analog electronic circuits Appeltant2011 , optical systems Larger2012 ; Paquot2012 , quantum-bit systems Fujii2017 , and magnetic systems Prychynenko2018 ; Nakane2018 ; Furuta2018a . Although these systems have shown good potential for reservoir computing, further improvements in computational capacity, operation speed, hardware size, energy consumption, and balance are required for practical applications.

Recently, reservoir computing on spin-torque oscillators (STOs) Kiselev2003 ; Slavin2009 has been reported Torrejon2017a ; Tsunegi2018a ; Markovic2019 ; Tsunegi2019 . (Applications of STOs to other types of computing have also been proposed Macia2011 ; Shibata2012 ; Csaba2013 ; Grollier2016a ; Kudo2017 ; Vodenicarevic2017b ; Romera2018 .) An STO is a nonlinear auto-oscillator with a power-dependent frequency shift. An STO exhibits microwave oscillation induced by a spin-polarized direct current via a spin-transfer torque, yielding microwave magnetic and electromagnetic fields. Compared with the other systems above, STOs have the advantages of nanosecond-scale fast dynamics, nanometer-scale size, and low energy consumption. The microwave fields can lead STOs to interact with each other in an array. Although the use of many interacting STOs is expected to achieve higher performance due to complex dynamics with many decrees of freedom, a single STO has been assumed in previous studies of STO reservoir computing Torrejon2017a ; Tsunegi2018a ; Markovic2019 ; Tsunegi2019 .

In this paper, we investigate reservoir computing on an array of STOs coupled by magnetic dipolar fields, as shown in Fig. 1, where couplings naturally appear when the STOs are integrated closely with each other. We perform numerical experiments on reservoir computing utilizing a nonlinear oscillator model Slavin2009 corresponding to an STO array. As a result, we show that reservoir computing can be performed in the synchronized oscillation state of the array and that its performance can be improved by increasing the number of STOs. The highest performance is obtained at the boundary between the synchronized and disordered states. The performance can become higher than that of a standard neural-network model known as an echo-state network (ESN) Jaeger2004 ; Maass2002 .

II Reservoir computing on STO array

II.1 Model of STO array

We place $N$ STOs on a square lattice with a size of ${M\times M(=N)}$ as shown in Fig. 1. In this paper, we assume that the STOs exhibit out-of-plane precession Kiselev2004 . These STOs emit large microwave dipolar fields Kanao2018 and output signals due to the magnetoresistive effect Zeng2013a ; Kubota2013 , which lead to strong couplings between the STOs and a high signal-to-noise ratio in measurements, respectively. The oscillations of the STOs are synchronized, with the same phase Belanovsky2012 ; Locatelli2015a ; Flovik2016 . [We show in Appendix A an example of such synchronized oscillations by use of a numerical simulation based on the Landau-Lifshitz-Gilbert-Slonczewski (LLGS) equation Slonczewski1996 .]

For reservoir computing, we utilize the dynamics of the oscillation powers and phases (or frequencies) of the STOs Torrejon2017a ; Tsunegi2018a ; Markovic2019 ; Tsunegi2019 . To calculate this dynamics, we use a corresponding nonlinear oscillator model Slavin2009 ; Kanao2016 , which can capture the essential features of the dynamics of the STOs, including the synchronization Slavin2006 , irrespective of the detailed magnetization configuration of the STOs. The nonlinear oscillator model can be derived from the LLGS equation under certain assumptions. (See Appendix B for details.) The powers and frequencies of the STOs usually vary slowly compared with the oscillation itself Slavin2009 ; Kanao2018 , and fast oscillation components can be eliminated in the nonlinear oscillator model. Thus this model allows us to reduce the computational costs of the simulation.

The nonlinear oscillator model for $N$ STOs is written as Slavin2009

[TABLE]

where ${t,c_{i},\delta\omega_{i}(p_{i}),\Gamma_{i}(p_{i})}$, and $\Omega_{ij}$ (${i=1,2,\cdots,N}$) are the time, complex amplitudes representing the slowly-varying magnetization dynamics, the nonlinear frequency shifts originating from magnetic anisotropy, the effective damping rates taking account of Gilbert damping and the spin-transfer torque, and the coupling strengths between the STOs caused by the dipolar fields, respectively. Here, a uniform high-frequency oscillation factor has been eliminated (see Sec. B.1). $\delta\omega_{i}(p_{i})$ and $\Gamma_{i}(p_{i})$ are supposed to depend on a dimensionless oscillation power ${p_{i}=\left|c_{i}\right|^{2}}$. Assuming a linear dependence of $\delta\omega_{i}(p_{i})$ and $\Gamma_{i}(p_{i})$ on $p_{i}$ Slavin2009 ; Kanao2016 , we express them as

[TABLE]

where $p_{i0},N_{f},\Gamma_{p}$, and ${\xi_{i}(t)}$ are the stationary oscillation power, the nonlinear frequency shift, the rate of damping of power, and the modulation by an input signal, respectively, in units of frequency Slavin2009 . $\delta\omega_{i}(p_{i})$ and $\Gamma_{i}(p_{i})$ describe the properties of a single STO, and ${\delta\omega_{i}(p_{i0})=\Gamma_{i}(p_{i0})=0}$ holds for its stationary oscillation state for a constant current. We add the input signal to the current, which, via the spin-transfer torque, leads to a modulation of $\Gamma_{i}(p_{i})$. This modulation is represented by $-\xi_{i}(t)$. It is supposed that $\delta\omega_{i}(p_{i})$ does not directly depend on the current, and instead is affected through the changes in $p_{i}$ Slavin2009 . Equations (2) and (3) can be regarded as linear approximations of $\delta\omega_{i}(p_{i})$ and $\Gamma_{i}(p_{i})$ for small deviations of $p_{i}$ from $p_{i0}$ and for an input signal in the current. (See Sec. B.2 for details.)

In the array, the dipolar fields between the nearest-neighbor STOs $\Omega$, and between the next-nearest-neighbor STOs, ${\Omega/\sqrt{8}}$, are taken into account, as the dipolar field is inversely proportional to the cube of the distance. Differences between the STOs are modeled by taking $p_{i0}$ as a Gaussian random variable with a mean of $p_{0}$. In this work, the input signal is applied to an STO at the corner of the array, namely STO1 in Fig. 1, and thus only $\xi_{1}(t)$ is modulated, and ${\xi_{i}(t)=0}$ for ${i\neq 1}$. Based on reported values Kanao2016 ; Quinsat2010 , we choose the following values: ${p_{0}=0.5}$, ${N_{f}/(2\pi)=0.1}$ GHz, ${\Gamma_{p}=0.1}$ ns*-1*, ${\left|\xi_{1}(t)\right|<0.1}$ ns*-1*, and ${\Omega=-0.1}$ ns*-1*, unless otherwise stated. The value of $\Omega$ is estimated in Sec. B.1. The standard deviation of $p_{i0}$ is set to 0.01. The denominator $p_{i0}$ in Eq. (3) is approximated by $p_{0}$. Equation (1) is numerically integrated with the initial condition ${c_{i}=\sqrt{p_{0}}}$, which is equivalent to all the magnetizations of a free layer of STOs being in the same direction. This situation can be realized in a strong uniform external field.

Reservoir computing utilizes the dynamics of the waveforms induced by the input signals. Figure 2 shows examples of waveforms of the input signal $\xi_{1}(t)$ and of nine STOs obtained from Eq. (1). The waveforms after a sufficiently long time are shown. The oscillation of $\mathrm{Re}\hskip 1.99997ptc_{i}$ in Fig. 2(a) is induced by the terms $\delta\omega_{i}(p_{i})$ and $\Omega_{ij}$ in Eq. (1). It can be seen that the oscillation phases are nearly synchronized. Fluctuations of the oscillations are induced by the input signal. Since the input signal is applied to STO1, the power of STO1 is directly modulated by the input signal, as shown in Fig. 2(b). The effects of the input signal propagate to the other STOs. When the input signal is constant, the dynamics of the powers is damped in several tens of nanoseconds. The timescale of the damping of the powers is determined by ${1/\Gamma_{p}(=10}$ ns). The waveforms of the STOs at symmetric positions (see Fig. 1), i.e., (STO2 and 4), (STO3 and 7), and (STO6 and 8), almost overlap, but not completely, because of the deviation of $p_{i0}$.

II.2 Reservoir computing

Next, we explain reservoir computing on the STO array. We evaluate its performance by following the setting used in Ref. Fujii2017 . The $k$th input signal, $s_{k}$, is introduced into STO1 at a rate $1/T_{s}$ during the time interval between ${T_{s}(k-1)}$ and $T_{s}k$. The powers are measured at the end of each interval, in order to fully incorporate the effects of the latest signal. The measured values are denoted by ${x_{ki}=p_{i}(T_{s}k)}$. The output data are ${y_{k}=\sum^{N}_{i=0}x_{ki}w_{i}}$, where $w_{i}$ is the $i$th weight for reservoir computing, and ${x_{k0}=1}$ is introduced to treat a constant term. The weights are determined in a training process by minimizing $\sum_{k}(y_{k}-\bar{y}_{k})^{2}$, so that a target series $\bar{y}_{k}$ is reproduced. These $w_{i}$ are given by ${{\bf w}=X^{+}\bar{{\bf y}}}$, where $X$ is a matrix and ${\bf w}$ and $\bar{{\bf y}}$ are column vectors, with components $x_{ki}$, $w_{i}$, and $\bar{y}_{k}$, respectively, and where $X^{+}$ is the Moore-Penrose pseudoinverse matrix of $X$ Fujii2017 . The dynamics of the STO array is calculated for a time of $5000T_{s}$. The first $1000T_{s}$ is discarded, and the subsequent $3000T_{s}$ and the last $1000T_{s}$ are used for training and evaluation, respectively.

We evaluate the performance by using short-term memory (STM) and parity-check (PC) tasks Fujii2017 ; Furuta2018a ; Tsunegi2018a ; Tsunegi2019 , because these tasks can estimate capacities for short-term memory and nonlinearity, which are necessary for applications to real-time computing. The STM task measures the STM capacity by outputting delayed inputs. The PC task can estimate nonlinearity because the required function is a nonlinear mapping. [We also evaluate a nonlinear autoregressive moving average (NARMA) task Fujii2017 ; Atiya2000 and show the result in Appendix C.] For the STM and PC tasks, $s_{k}$ is a random sequence of [math]s and $1$s. We convert $s_{k}$ to $\xi_{1}(t)$ by use of a square wave with value ${\xi_{1}(t)=0.1(2s_{k}-1)}$ ns*-1*. The targets of the STM and PC tasks are ${\bar{y}_{k}=s_{k-\tau}}$ and ${\bar{y}_{k}=Q\left(\sum^{\tau}_{l=0}s_{k-l}\right)}$, respectively, for ${\tau=0,1,2,\cdots}$. The parity function $Q(S)$ is 0 for even $S$ and 1 for odd $S$. The capacities are measured by ${C=\sum^{\tau_{\rm max}-1}_{\tau=0}C_{\tau}}$, where $C_{\tau}$ is the square of the correlation coefficient of $\bar{y}_{k}$ and $y_{k}$ called a $\tau$-delay capacity ${\left(0\leq C_{\tau}\leq 1\right)}$. Here, ${\tau_{\rm max}=500}$. A larger $C_{\tau}$ means higher performance, and $C_{\tau}$ reaches its maximum value of 1 when ${\bar{y}_{k}=y_{k}}$. The performances are evaluated for 20 sets of $\left(p_{i0},s_{k}\right)$, and the average values of the performances and standard deviations are calculated.

III Results and discussion

III.1 Performance improvement with increasing $N$

Figures 3(a) and 3(b) show the STM and PC capacities as a function of the number of STOs $N$ and the sampling time $T_{s}$. The STM capacity increases with $N$ in the range shown (${N\leq 400}$). The PC capacity saturates at large $N$. The behaviors of the capacities at larger $N$ are discussed below. The STM and PC capacities depend on $T_{s}$ in opposite ways; namely, the STM capacity decreases with increasing $T_{s}$, while the PC capacity increases. The STM capacity decreases with increasing $T_{s}$ because the effect of the input signal decays when $T_{s}$ is long. On the other hand, the PC capacity increases with increasing $T_{s}$ because a longer $T_{s}$ allows larger variations of the powers such that nonlinear effects appear. These behaviors occur in the range of $2$ ns ${\lesssim T_{s}\lesssim 50}$ ns, which can be understood from the timescale of the damping of power, ${1/\Gamma_{p}=10}$ ns. Outside this range of $T_{s}$, the powers are no longer susceptible to influence from the input signals.

The data in Figs. 3(a) and 3(b) are summarized in Fig. 3(c). The capacities are distributed below a curve on which the product of the STM and PC capacities is a constant (equal to $85$), implying a trade-off relation between the STM and PC capacities. This kind of relation has been reported as a memory-nonlinearity trade-off Dambre2012 .

Figure 4 shows the capacities at each $\tau$ (the $\tau$-delay capacities) for ${N=400}$. The capacities decrease from 1 to 0 with increasing $\tau$ for ${T_{s}\geq 10}$ ns ${=1/\Gamma_{p}}$. The capacities for ${T_{s}<1/\Gamma_{p}}$ are small compared with 1 even at small $\tau$, which indicates that $T_{s}$ is too short for the powers to fully respond to the input signals. The STM capacities for ${T_{s}\leq 2}$ ns show finite values at large $\tau$, indicating that the effects of the input signals remain.

The largest values of the STM and PC capacities obtained in the present study are $46.4$ at ${N=400}$ and ${T_{s}=2}$ ns, and $6.4$ at ${N=196}$ and ${T_{s}=50}$ ns, respectively. These values are larger than the values for an ESN with 500 units reported in Ref. Fujii2017 , which are about $17.6$ and $6.0$, respectively. Remarkably, the largest value of the STM capacity here is larger than twice the value for the ESN.

A possible reason for this large STM capacity is that the STOs have a phase degree of freedom, where the phases are time integrations of the oscillation frequencies in the past. In the STO array, the phases affect the time evolution of the powers, and a memory for longer time can be obtained from the powers. In an ESN, on the other hand, each unit has only one variable and shows damped motion in reservoir computing, which suppresses the past memory. Studies of reservoir computing on general nonlinear oscillator networks are an important topic for future work. Another difference which can lead to higher performance of an STO array than of an ESN is in the network structure. The STOs here are coupled by short-range interaction. This kind of weak connection allows various cooperative dynamics suitable for reservoir computing Jaeger2004 , which is discussed in the next section (Sec. III.2). On the other hand, the ESN used for comparison is coupled by all-to-all connections, which tend to suppress the diversity of the dynamics.

III.2 Performance enhancement near transition

Finally, we investigate the dependence of the STM and PC capacities on the coupling, and find that the capacities depend strongly on the coupling strength. Figures 5(a) and 5(b) show the STM and PC capacities as a function of the coupling strength $|\Omega|$ and of $N$ for ${T_{s}=5}$ ns. The largest capacities appear around ${|\Omega|=0.14}$ ns*-1*.

To clarify the reason for these dependences of the capacities, we examine the oscillation states of the STO array under the corresponding conditions. The oscillation states are characterized by an order parameter $\rho$ Kuramoto1987 ; Acebron2005 and an average power ${\bar{p}}$,

[TABLE]

which are shown in Figs. 5(c) and 5(d), respectively. The order parameter $\rho$ characterizes the uniformity of the oscillation phase, and reaches its largest value of $1$ when all the STOs are perfectly synchronized.

Figure 5(c) shows that $\rho$ is nearly $1$ in most of the region, indicating that all $c_{i}$ have nearly the same phases and thus the oscillations are synchronized. $\rho$ decreases as $N$ is increased for the same $|\Omega|$. The regions where the STM and PC capacities are finite in Figs. 5(a) and 5(b) are included in the region where $\rho$ is rather large, ${\left(\rho\gtrsim 0.5\right)}$, in Fig. 5(c).

To see more clearly the relation between the performance and the oscillation state, we plot the capacities and the order parameter as a function of $|\Omega|$ in Fig. 6, where ${N=400(=20^{2}),784(=28^{2})}$, and ${1600(=40^{2})}$. As shown in Figs. 6(a) and 6(b), with increasing $|\Omega|$ the STM and PC capacities increase almost monotonically for ${|\Omega|\lesssim 0.15}$ ns*-1* and then suddenly decrease. As shown in Fig. 6(c), the order parameter $\rho$ exhibits similar behavior, which indicates a transition from the synchronized oscillation state to a disordered state around ${|\Omega|=0.15}$ ns*-1* (see below for details). From the results in Figs. 5 and 6, we conclude that reservoir computing can be performed in the synchronized oscillation state of the STO array and the performance is highest near the transition to the disordered state, which is consistent with the reported behavior that the capacities of reservoir computing can be enhanced at the edge of this kind of transition Verstraeten2007 ; Appeltant2011 .

The dependence of the oscillation state on the coupling strength observed in Fig. 6(c) can be understood as follows. The explanation is divided into two parts, for small and large $|\Omega|$. For small $|\Omega|$, the power is almost constant, ${\bar{p}\simeq p_{0}=0.5}$, as seen in Fig. 5(d), and thus the phase freedom is important. Then, the model of Eq. (1) can be reduced to a two-dimensional Kuramoto model with a short-range interaction Slavin2009 ; Kuramoto1987 ; Acebron2005 . In this Kuramoto model, it is known that stronger couplings lead to more ordered synchronization, while a larger number of oscillators causes a more disordered oscillation state Lee2010 ; Flovik2016 . These results for the Kuramoto model can explain the dependence of $\rho$ on $|\Omega|$ and $N$ for ${|\Omega|\lesssim 0.15}$ ns*-1* in Fig. 6(c), where $\rho$ increases with $|\Omega|$ and decreases with $N$. For large $|\Omega|$, on the other hand, the average power deviates from ${p_{0}=0.5}$, as seen in Fig. 5(d), and thus nonlinear effects of $\delta\omega_{i}(p_{i})$ and $\Gamma_{i}(p_{i})$ in Eq. (1) can appear. In this situation, it is known that the power degree of freedom causes disordered unsteady states, including chaotic states Matthews1990 ; Acebron2005 . In these states the synchronization is suppressed, which corresponds to low $\rho$ for ${|\Omega|\gtrsim 0.15}$ ns*-1* in Fig. 6(c). These disordered unsteady states Matthews1990 ; Acebron2005 can also explain the large standard deviations observed for the same region in Fig. 6(c), which indicates that the oscillation states are sensitive to small changes in the conditions. Between the synchronized state and the disordered state, for small and large $|\Omega|$, respectively, a transition occurs, causing a sharp change in $\rho$.

The decrease of $\rho$ with increasing $N$ is due to the short-range interaction between the STOs. By designing the STOs in the array to couple through a long-range interaction, e.g., by using a common electrode Kudo2017 , one can make $\rho$ finite for large $N$ Kuramoto1987 ; Acebron2005 , allowing the synchronization of a large number of STOs.

IV Summary and outlook

We investigate reservoir computing on an STO array coupled by dipolar fields. Numerical experiments on reservoir computing are performed by using a corresponding nonlinear oscillator model, which exhibits synchronized oscillation. As a result, we show that the performance of reservoir computing is improved by increasing the number of STOs in the synchronized oscillation state. The values of the STM and PC capacities for ${N\leq 400}$ can be larger than those for an ESN with 500 units. In particular, the maximum STM capacity is more than twice that of an ESN. The largest capacities are obtained at the boundary between the synchronized and disordered states. We believe that the results can be explained by an enhancement of the capacities at the edge of chaos. The present results indicate that STO arrays are promising for the hardware implementation of reservoir computing. The performance can be further improved by introducing methods in reservoir computing not included here, such as time multiplexing. In this paper, the effects of thermal fluctuations Quinsat2010 , which are large in nanometer-sized magnetizations, are not included. This issue is left for future work.

To realize the reservoir computing on an STO array, further considerations of several factors are necessary. We expect that an array of pillar-shaped magnetic tunnel junctions (MTJs) like those in magnetoresistive random access memories Bhatti2017 or an array of giant magnetoresistive elements can be utilized as the STO array. MgO-based MTJs can exhibit large output signals due to the tunneling magnetoresistive effect. Typical ferromagnetic materials based on elements such as Fe and Co will yield enough coupling by dipolar fields, as estimated in Appendixes A and B. At the current time, a few STOs have been coupled by dipolar fields Locatelli2015a . To increase the number of STOs, a one-dimensional array is the first candidate because the design of the electrodes which apply currents and extract signals will be simpler. With an increasing number of lines, the array approaches a two-dimensional one. Three-dimensional design of the electrodes will be necessary for integrating a large number of STOs into a two-dimensional array. Instead, reservoir computing may be possible when the number of measured signals is smaller than the number of STOs, for example by measuring signals from a fraction of the STOs or measuring sums of signals from several STOs, which reduces the complexity of the electrode configuration. In this case, estimation of the performance degradation caused by the reduced number of measured signals is necessary. Although we have assumed that the input signals are applied to one STO, reservoir computing can be performed by applying input signals to several STOs, which increases the freedom of design of the array. Input signals with waveforms other than a square wave, such as a distorted square wave or a triangle wave, will introduce similar fluctuations of the powers when the rates and amplitudes are suitable. In the measurement of the powers as the outputs, signal-processing technologies for amplitude modulation can be used Nakamura2018 .

Acknowledgements.

The authors would like to thank S. Tsunegi, M. Goto, K. Kudo, and R. Sato for valuable discussions.

Appendix A SYNCHRONIZATION OF STOS IN LLGS EQUATION

We show an example of the magnetization dynamics of STOs by solving the LLGS equation. The LLGS equation for the magnetizations of a free layer of STOs ${\bf m}_{i}$ ($i=1,2,\cdots,N$) normalized by the saturation magnetization $M_{s}$ is given by

[TABLE]

where $N$, $t$, $\gamma$, and $\alpha$ are the number of STOs, the time, the gyromagnetic ratio, and the Gilbert damping, respectively. ${\bf H}^{\rm{eff}}_{i}$ and ${{\bf T}^{S}_{i}=\sigma I_{i}{\bf m}_{i}\times({\bf m}_{i}\times{\bf e}_{p})}$ are an effective field and a spin-torque term with efficiency $\sigma$, current $I_{i}$, and unit vector in the direction of the spin polarization of the current ${\bf e}_{p}$.

In the example, we choose an array of typical STOs with in-plane magnetized free and fixed layers, having a rectangular shape. A small (${50^{2}\times 10}$ nm3) and thus uniform ${\bf m}_{i}$ is assumed. A perpendicular field $H_{z}$ is applied, yielding an out-of-plane precession (OPP) for large $I_{i}$ Kiselev2004 . The STOs are placed on a square lattice, and all the magnetizations ${\bf m}_{i}$ are coupled by dipolar fields Kudo2012 ; Kudo2015 . $I_{i}$ and the initial directions of ${\bf m}_{i}$ are randomized. Figure 7 shows an example of the waveforms. Owing to the randomness, the waveforms are initially disordered. After several oscillation periods, the oscillations are synchronized. This state is obtained for other initial conditions, indicating the robustness of the synchronized state. This kind of synchronization has been reported in vortex-type STOs coupled by dipolar fields Belanovsky2012 ; Locatelli2015a ; Flovik2016 .

Appendix B FROM THE LLGS EQUATION TO THE NONLINEAR OSCILLATOR MODEL

B.1 Coupling between STOs

Here, we show the correspondence between the LLGS equation and the nonlinear oscillator model. In the OPP, the $x$ and $y$ components of ${\bf m}_{i}$ are mainly oscillating, and we introduce complex amplitudes ${\tilde{c}_{i}=m_{ix}-\mathrm{i}m_{iy}}$. From Eq. (6), the equation of motion of $\tilde{c}_{i}$ is

[TABLE]

where the term including ${\bf H}^{\rm{eff}}_{i}$ is transformed exactly, while the damping and spin-torque terms are represented by using an effective damping rate $\Gamma_{i}$, which is discussed in Sec. B.2. Here, we assume uniform magnetizations in each STO. In the OPP, $H^{\rm eff}_{iz}$ mainly determines the precession axis and frequency ${\tilde{\omega}_{i}=\gamma H^{\rm eff}_{iz}}$. Besides, ${H^{\rm eff}_{ix}-\mathrm{i}H^{\rm eff}_{iy}}$ is due to the dipolar fields from the other STOs. Although the prefactor includes $m_{iz}$, this prefactor depends little on $\tilde{c}_{i}$, because $m_{iz}$ is almost constant for the OPP. The term including ${H^{\rm eff}_{ix}-\mathrm{i}H^{\rm eff}_{iy}}$ reduces to the coupling term as follows. We approximate, for simplicity, the magnetization of a free layer of volume $V$ by a point magnetic moment ${{\bm{\mu}}=M_{s}V{\bf m}}$. The dipolar field at the position ${{\bf r}=(x,y,z)}$ arising from this moment at the origin is given by

[TABLE]

Then,

[TABLE]

On the same plane, i.e., ${{\bf r}=(x,y,0)}$,

[TABLE]

Thus Eq. (7) can be written as

[TABLE]

where ${\Omega_{ij}=-\gamma m_{iz}M_{s}V/\left(2|{\bf r}_{i}-{\bf r}_{j}|^{3}\right)}$ and $\Xi_{ij}=-3\gamma m_{iz}M_{s}V\left[x_{i}-x_{j}-\mathrm{i}\left(y_{i}-y_{j}\right)\right]^{2}/\left(2|{\bf r}_{i}-{\bf r}_{j}|^{5}\right)$. Here ${\bf r}_{i}$ represents the position of STO $i$.

We assume that the STOs generate a common frequency $\omega$ when isolated in the steady state, and that the variation of the frequencies is caused by the nonlinear frequency shifts $\delta\omega_{i}$. Thus, the oscillation frequencies are represented by ${\tilde{\omega}_{i}=\omega+\delta\omega_{i}}$. Then, the uniform oscillation with frequency $\omega$ is decoupled by ${\tilde{c}_{i}=c_{i}\mathrm{e}^{-\mathrm{i}\omega t}}$ as

[TABLE]

Note that an oscillating factor with a frequency $2\omega$ appears in the coefficient of $c^{*}_{j}$. This term does not affect the dynamics of $c_{i}$ when $c_{i}$ varies slowly, as shown in the following. We assume that the timescale of the change of $c_{i}$ is much longer than the oscillation period ${2\pi/\omega}$ and that $c_{i}$ can be approximated by a constant during $2\pi/\omega$. Then, by averaging Eq. (12) for a time $\pi/\omega$, the term including $c^{*}_{j}$ vanishes owing to the factor $\mathrm{e}^{2\mathrm{i}\omega t}$, and hence the nonlinear oscillator model,

[TABLE]

is obtained. This averaging is equivalent to using a variable

[TABLE]

The numerical order of magnitude of $\left|\Omega_{ij}\right|$ can be estimated as follows. In the typical case of ${m^{2}_{iz}=0.5}$, ${M_{s}=1.6\times 10^{3}}$ emu/cm3, and ${V=50^{2}\times 10}$ nm3, $\left|\Omega_{ij}\right|$ is $0.031$ ns*-1* and $0.25$ ns*-1* for center-to-center distances of $200$ nm and $100$ nm, respectively, indicating that $\left|\Omega_{ij}\right|$ is of the order of $0.01$ ns*-1* to $0.1$ ns*-1*.

B.2 Effective damping rate and nonlinear frequency shift

The effective damping rate $\Gamma_{i}$ is assumed to be a function of the current $I_{i}$, while the nonlinear frequency shift $\delta\omega_{i}$ is not Slavin2009 . The reason for this can be illustrated as follows. Since $\Gamma_{i}$ and $\delta\omega_{i}$ describe the properties of a single STO, we focus on the single-STO case, in which the LLGS equation (6) is written as

[TABLE]

Here, the Gilbert damping term is transformed to this form [the second term in the right-hand side of Eq. (15)] by assuming ${\alpha^{2}\ll 1}$ and ${|\alpha\sigma I|\ll 1}$, which is valid for typical STOs. We consider the OPP and assume ${\bf H}^{\rm{eff}}$ in the $z$ direction. By supposing ${{\bf e}_{p}=e_{px}\hat{{\bf x}}+e_{pz}\hat{{\bf z}}}$ (where $\hat{{\bf x}}$ and $\hat{{\bf z}}$ are unit vectors in the $x$ and $z$ directions), the spin-torque term can be divided into two terms. The term including $e_{pz}\hat{{\bf z}}$ is collinear with the Gilbert damping term and can balance and cancel that term, allowing a stationary oscillation state. From these terms, the effective damping rate can be expressed as ${\Gamma=\left(\alpha\gamma H^{\rm{eff}}_{z}-\sigma Ie_{pz}\right)m_{z}}$, which depends explicitly on $I$ and is almost constant for the OPP. On the other hand, the term including $e_{px}\hat{{\bf x}}$ can have a component parallel to the first torque term ${-\gamma{\bf m}\times{\bf H}^{\rm{eff}}}$, which determines the oscillation frequency. The total contribution of the $e_{px}\hat{{\bf x}}$ term during one precession, however, is small because the component in the direction of the first torque term changes sign during the precession and cancels out Taniguchi2013 . The component in the direction of the first torque is proportional to ${-\left({\bf m}\times{\bf H}^{\rm{eff}}\right)\cdot\left[{\bf m}\times\left({\bf m}\times e_{px}\hat{{\bf x}}\right)\right]=H^{\rm{eff}}_{z}e_{px}m_{y}}$, and the sign of $m_{y}$ changes for the OPP during one precession. Thus the dependence of the oscillation frequency on the current can be neglected, and it is assumed that $\delta\omega$ is not a explicit function of $I$. Note that $\delta\omega$ is indirectly affected by $I$ as a result of changes in the power $p$.

The forms of $\Gamma_{i}$ and $\delta\omega_{i}$ used in Eqs. (2) and (3) are obtained as follows. In a stationary oscillation state of an individual STO with power $p_{i0}$ and current $I_{i0}$, $\Gamma_{i}$ is zero because the spin-transfer torque and the Gilbert damping are balanced and cancel each other. When a small modulation of $I_{i}$ is introduced as an input signal, $p_{i}$ deviates from $p_{i0}$ as a result of time evolution. Then $\Gamma_{i}$ can be linearly approximated as

[TABLE]

where the first term is expressed by using a damping rate of the power Slavin2009

[TABLE]

We express the second term yielded by the input signal as

[TABLE]

The minus sign is because the spin-transfer torque works as a negative damping, which is proportional to the current. As we assume that ${\delta\omega_{i}=0}$ at ${p_{i}=p_{i0}}$ and $\delta\omega_{i}$ does not depend explicitly on $I_{i}$, $\delta\omega_{i}$ is approximated as ${\delta\omega_{i}\simeq N_{f}\left(p_{i}-p_{i0}\right)}$. This kind of description of the current modulation has been developed in Ref. Slavin2009 .

Appendix C NARMA TASK

Here, we show the performance in the NARMA task, where series defined by nonlinear recurrence relations are emulated, and the ability to predict nonlinear time series is measured. We perform a NARMA task defined by the following equation: $y_{k+1}=0.4y_{k}+0.4y_{k}y_{k-1}+0.6s^{3}_{k}+0.1$ for $s_{k}$ consisting of uniform random numbers in ${(0,0.2)}$, called the NARMA2 task Fujii2017 ; Atiya2000 . We convert $s_{k}$ to $\xi_{1}(t)$ by use of a square wave with value ${\xi_{1}(t)=s_{k}-0.1}$ ns*-1*. The performance is quantified by the normalized mean-squared error ${\epsilon=\left[\sum_{k}\left(\bar{y}_{k}-y_{k}\right)^{2}\right]/\left(\sum_{k}\bar{y}^{2}_{k}\right)}$, which is smaller when $\bar{y}_{k}$ and $y_{k}$ are closer, indicating higher performance.

The results for the NARMA2 task are shown in Fig. 8. The error for ${T_{s}=20}$ ns decreases with increasing $N$ by three orders of magnitude, and the value of the error is comparable to that in previous studies Fujii2017 .