A Frequency-Phase Potential for a Forced STNO Network: an Example of Evoked Memory
AA Frequency-Phase Potential for a ForcedSTNO Network: an Example of EvokedMemory ∗ Frank HoppensteadtCourant Institute of Mathematical SciencesNew York University
Abstract
The network studied here is based on a standard model in physics,but it appears in various applications ranging from spintronics to neu-roscience. When the network is forced by an external signal commonto all its elements, there are shown to be two potential (gradient) func-tions: One for amplitudes and one for phases. But the phase potentialdisappears when the forcing is removed. The phase potential describesthe distribution of in-phase/anti-phase oscillations in the network, aswell as resonances in the form of phase locking. A valley in a potentialsurface corresponds to memory that may be accessed by associativerecall. The two potentials derived here exhibit two different forms ofmemory: structural memory (time domain memory) that is sustainedin the free problem, and evoked memory (frequency domain memory)that is sustained by the phase potential, only appearing when the sys-tem is illuminated by common external forcing. The common forcingorganizes the network into those elements that are locked to forcingfrequencies and other elements that may form secluded sub-networks.The secluded networks may perform independent operations such aspattern recognition and logic computations. Various control methodsfor shaping the network’s outputs are demonstrated. ∗ KEYWORDS: Spin Torque Nano Oscillators, neuromorphic engineering, gradient po-tential, spintronics, phase-locked loops, neural energy. Supported in part by NeurocircLLC. I thank Eugene Izhikevich for commenting on an earlier version of this paper. a r X i v : . [ n li n . AO ] A ug n the terminology used here, an oscillation is a signal that is character-ized by amplitude, frequency, and phase deviation (timing), and an oscillatoris any device that may produce such signals. An oscillatory neural network(ONN) is a network of oscillators designed to achieve some objective, such aslearning, storing and retrieving memory, computing, controlling subsystems,or switching. The oscillator network studied here is based on an electromag-netic oscillator called a spin torque nano oscillator (STNO). The space scaleof an STNO is 100 or fewer nanometers and one operates at GHz frequencies.Standard methods are used to describe the response of the oscillators whenthey are connected to others and when the array is exposed to a commonoscillating background. The responses are described in terms of potential(gradient) functions for amplitudes and phases, where valleys in the poten-tial manifolds correspond to stable patterns of electromagnetic oscillations.Neuromorphic aspects of this work are partly guided by the following ob-servations from neuroscience: Neurons and networks of them may have oscil-latory behavior [1]; neural tissue has stability properties, including phase andfrequency locking [2]; EEGs suggest that particular background frequencybands may be important in brain functions [3]; resonance phenomena arepresent in neural activity [4, 5]; and phase coding has been identified in thebrain [4, 6, 7, 8]. There are also hardware considerations: A single oscillatormay efficiently replace a larger network of off-on devices [9, 10]; and, there arewell-developed technologies for frequency-based computation, control, signalprocessing, phase detection, and information flow.An interesting question in brain science that partly motivates this workis: How might interactions between background fields, whatever their source,and networks of electrically active cells function? The specific network ofspin-torque nano-oscillators studied here gives some insight to the function-ality of such a structure. A common background signal is shown to enable anetwork to process and classify information, create and sustain memory, andcontrol the flow of information using frequency based operations. This is anextension of work in [11, 12, 13]. This is not a study of the brain.The novelty here is deriving a frequency-phase landscape for an STNOnetwork and simplifying its analysis by introducing a phase factor that un-folds singularities in the STNO model.The STNO model is described first, then it is converted to phase-amplitudecoordinates. The amplitude potential is derived, and then the system isfurther reduced to a flow on a torus where the phase potential is derived.Approximations to the system’s states are derived, and these are used to2stimate the power spectrum density and to describe various controls forbursting, self-organization, and hysteresis. Simulations of the system illus-trate its phase potential. Next, a uniform forcing is applied to a weaklyconnected network of STNOs, and various patterns of phase locking andin-phase/anti-phase responses are described. When the external forcing isreplaced by the combined output of the oscillators, referred to here as globalfeedback, then the phase deviations between oscillators are shown to be syn-chronized and intervals of phase locking appear. Memory, information flow,self-organization, computation, perturbations by random noise, and contin-uum models are discussed. The STNO model studied here is derived from the Landau-Lifshitz-Gilbert-Slonczewski model (LLGS) [14, 15]. The model describes activity in a thinlayer of material that contains magnetic elements. A magnetic element hasa magnetization vector in three dimensions, and a complex number z is usedto represent the projection of the magnetization vector onto the layer. Theoscillator is at rest when z = 0, meaning that the vector is orthogonal to thelayer. According to Maxwell’s equations an electric current in the directionof the magnetization may apply torque to the vector causing it to precess.Oscillations in magnetization are useful in a number of applications includingcomputation and memory [15].The model is [13, 14, 15, 16]˙ z = b z | z | + λz − | z | z + g ( | z | )) z + εf ( t ) (1)where ˙ z = dz/dt , = √−
1, the parameters b and λ represent parametricforcing, the frequency modulation term g is a smooth, real-valued function,and external forcing is represented by f , that has effect when the system isotherwise at rest ( z = 0 here). The phase factor, bz/ | z | . Including the phase factor bz/ | z | , whichdoes not appear in [13, 14, 15, 16], is based on the following reasoning: Theequilibrium amplitudes ( R ) with b = 0 , ε = 0, are determined by solving theequation 0 = λR − R λ = 0 , R = 0. There are three real solutions ofthis equation (counting multiplicity), R = 0 , (cid:112) λ + , e π (cid:112) λ + , where λ + =( λ + | λ | ) / λ . These solutions may be designated aszero, spin up, and spin down, respectively, to acknowledge their signs. Theuniversal unfolding of this singularity is0 = b + λR − R where b is the unfolding parameter. When b >
0, there is a unique positiveroot of this equation, denoted here by R = R ∗ ( b, λ ) . It follows that R ∗ > (cid:112) λ + and lim b → +0 R ∗ ( b, λ ) = (cid:112) λ + . There are two neg-ative real roots that appear through a fold bifurcation when λ is sufficientlylarge; those will be considered elsewhere. The phase factor simplifies calcu-lations since 1 /R often appears in the course of analysis, and R = 0 cannothappen if b (cid:54) = 0. However, the phase factor does introduce into equation (1)an essential singularity at z = 0, which is repelling.The phase factor bz/ | z | preserves rotational symmetry in the free problem( f ≡ z is a solution, then so is ze ν for any real number ν . In particular,both z and its anti-phase twin ze π = − z are solutions of the free problem.A calculation (not presented here) shows that z ( t, b ) − z ( t,
0) = O ( b / ) + O ( ε )uniformly in t ≥ t ∗ ( b ) >
0, where t ∗ ( b ) goes to zero as b → +0, accounts foran initial transient [17]. No physical argument is presented here for includingthe phase factor. External forcing . The external (or background) forcing f comprises M modes of oscillation: f ( t ) = M (cid:88) m =1 C m e (Ω m t +Ψ m ) (2)having (real) amplitudes C m , frequencies Ω m , and phase deviations Ψ m .Without loss of generality, the frequencies may be relabelled so Ω < Ω < · · · < Ω M . These frequencies may be incommensurable, in which case f isquasi-periodic. Also, C m ≥ , for m ∈ , M , since their signs may be set The notation 1 , M = { , . . . , M } is used.
4y adjusting Ψ m by adding π or not. The background field is assumed to beweak relative to the dynamics of the oscillator, as represented by the smallpositive parameter ε in (1). Phase-amplitude coordinates.
Changing z to polar coordinates, z ( t ) = R ( t ) exp( θ ( t )), takes equation (1) into˙ R + R ˙ θ = (cid:2) b + ( λ + g ( R )) R − R (cid:3) + ε M (cid:88) m =1 C m e (Ω m t +Ψ m − θ ) , after dividing both sides by exp( θ ). Equating real and imaginary parts givesa system of two equations, one for the amplitude R ˙ R = ( b + λR − R ) + ε M (cid:88) m =1 C m cos(Ω m t + Ψ m − θ ) , (3)and one for the phase θR ˙ θ = g ( R ) R + ε M (cid:88) m =1 C m sin(Ω m t + Ψ m − θ ) . (4) Let z be determined from equation (1) where f is given in (2). Following aresome facts about this system. The amplitude potential function used here is defined by V ( R ; b, λ ) = (cid:18) R − λ R − bR (cid:19) . (5)With this, equation (3) may be rewritten as˙ R = − ∂V∂R ( R ; b, λ ) + ε M (cid:88) m =1 C m cos(Ω m t + Ψ m − θ ) , ε → ≤ t < ∞ . Thepotential function V ( R ; b, λ ) describes a landscape of amplitudes indexed by b and λ .The isolated minima of V define equilibria for the system when ε = 0,and they are stable under persistent disturbances. Stable under persis-tent disturbances here means that small perturbations of a gradient system˙ x = −∇ F ( x ), say ˙ y = −∇ F ( y ) + εg ( t, y ), preserve stable behavior nearminima of F . That is, if x ∗ is an isolated minimum of F and y (0) is near x ∗ , then, under natural conditions on g , y ( t ) = x ∗ + error, where the errorgoes to 0 as ε → ≤ t < ∞ . Stability under persistentdisturbances is justified by Lyapunov’s stability theory since the potentialdefines a Lyapunov function at each of its minima [18, 19]. The convergence of R to R ∗ occurs first among the hierarchy of time scalesin this problem as ε → b >
0, and in most of the calculationsdone here we assume that R is at its asymptotic value R = R ∗ ( b, λ ). Theremainder of this paper is based on (1) with g ( | z | ) replaced by the constant ω ∗ = g ( R ∗ ( b, λ )), since the difference R ( t ) − R ∗ = O (exp(( λ − R ∗ ) t/ε ).The surface R = R ∗ defines a torus that is indexed by θ and (cid:126) Ω t where (cid:126) Ωis the vector of forcing frequencies, and the model studied in the remainderof this section is ˙ θ = ω ∗ + ε M (cid:88) m =1 ˆ C m sin(Ω m t + Ψ m − θ ) (6)where ˆ C m = C m /R ∗ . This equation represents the reduction of (1) to adifferential equation on an M + 1-dimensional torus.The differences in forcing frequencies are assumed to satisfy | Ω k − Ω l | (cid:29) ε for all k (cid:54) = l . Setting θ = ω ∗ t + φ in equation (6) gives an equivalent equation,but now one for the phase deviation φ . Equation (6) becomes˙ φ = ε M (cid:88) m =1 ˆ C m sin((Ω m − ω ∗ ) t + Ψ m − φ ) . (7) For example, the perturbation g ( t, y ) is bounded, L as a function of t , and differen-tiable as a function of y . ≤ t ≤ T /ε for some fixedtime T ; that is, for the slow time variable s = εt ranging over 0 ≤ s ≤ T . Themethod of averaging may be applied to (7) as follows. Define the averagingoperator acting on an integrable function F ( t ) to be (cid:68) F (cid:69) ≡ εT (cid:90) T/ε F ( t (cid:48) ) dt (cid:48) . Applying this operator to (7) while viewing φ = φ ( s ) as being an independentfunction of the slow time s , gives (cid:68) dφds (cid:69) = M (cid:88) m =1 ˆ C m (cid:68) sin((Ω m − ω ∗ ) t + Ψ m − φ ( s )) (cid:69) . This average evaluates to d ¯ φds = M (cid:88) m =1 ˆ C m (cid:18) cos( ¯ φ − Ψ m ) − cos( ¯ φ − Ψ m − ν m ) ν m (cid:19) (8)where ¯ φ replaces φ in the average and ν m = (Ω m − ω ∗ ) T /ε .This equation may be rewritten as d ¯ φds = − ∂W∂ ¯ φ ( ¯ φ, ω ∗ , ε )where W ( ¯ φ, ω ∗ , ε ) = − M (cid:88) m =1 ˆ C m sin( ¯ φ − Ψ m ) − sin( ¯ φ − Ψ m − ν m ) ν m . (9)An equivalent form of W that is useful for analysis and simulation is W = − M (cid:88) m =1 ˆ C m ( cosc( ν m ) sin( ¯ φ − Ψ m ) + sinc( ν m ) cos( ¯ φ − Ψ m )) . (10) The unnormalized cardinal functions used here are sinc x ≡ sin x/x , and cosc x ≡ (1 − cos x ) /x . Note that sinc(0) = 1 , cosc(0) = 0, and cosc( x ) = ( x/
2) sinc( x/ . | φ ( t ) − ¯ φ ( εt ) | = O ( ε ) uniformly for 0 ≤ t ≤ T /ε [19]. The isolated minima of W are equilibria for ¯ φ , and these arestable under persistent disturbances. Note for later use that ∂W∂ω (cid:12)(cid:12)(cid:12) ω = ω ∗ = Tε M (cid:88) m =1 ˆ C m ( cosc ∗ ( ν m ) sin( ¯ φ − Ψ m )+ sinc ∗ ( ν m ) cos( ¯ φ − Ψ m )) (11)where sinc ∗ ( x ) = d sinc( x ) /dx and cosc ∗ ( x ) = d cosc( x ) /dx . How wide are the valleys in W near its minima? Define X K ( ω ∗ ) = (Ω K − ω ∗ ) R ∗ εC K . It is shown in Appendix 1 that if the phase locking condition | X k | < θ ≈ Ω K t + φ ∗ K where φ ∗ K = Ψ K − arcsin X K , and the output frequency is Ω K . The error inthis approximation is shown in Appendix 1.In summary, if (12) holds, then the amplitude approaches R ∗ , the fre-quency is locked to Ω K , and the phase deviation φ ∗ K = Ψ K − arcsin X K isstable. In this situation z = R ∗ e (Ω K t + φ ∗ K ) + remainderwhere R ∗ and φ ∗ K are at or near stable minima of V and W , respectively, andthe remainder is small for large T R ∗ /ε and for small (Ω K − ω ∗ ).Estimates of the width of wells in W may be derived from (12). They arethe frequency intervals Ω m − εC m R ∗ < ω < Ω m + εC m R ∗ for m ∈ , M ; this interval lies within a well of W located near Ω m provided | εC m /R ∗ | is sufficiently small. This shows to what extent the power in theforcing signal ( | C m | ) controls the size of phase locking intervals.8ather than a detailed analysis of the extrema of W , their nature isdemonstrated in the limit T /ε → ∞ , where W in (10) becomes W ∞ ( ¯ φ, ω ) = − M (cid:88) m =1 ˆ C m δ ( ω − Ω m ) cos( ¯ φ − Ψ m ) . (13)where δ is Kronecker’s function: δ (0) = 1 but δ ( x ) = 0 for x (cid:54) = 0. Thefunction W ∞ has minima at ω = Ω m , ¯ φ = Ψ m for m ∈ , M , and their depthis proportional to ˆ C m ≡ C m /R ∗ . By continuity, W will have minima nearthese values, as shown by direct calculation that is not done here. Power spectrum . The power spectrum density (PSD) of z is derived inAppendix 2. The result is that P SD z ( ω ∗ ) ≈ R ∗ sinc ((Ω K − ω ∗ ) T /ε ) + O ( ε/T ) . (14)for ω ∗ ∼ Ω K . The notation ω ∗ ∼ Ω K means that the phase locking condition(12) is satisfied.It is useful to consider also the rotation number for this system ρ ( ω ∗ ) = lim t →∞ θ ( t, ω ∗ ) t since θ is accessible in this simple model. This represents the output fre-quency of the oscillator, and it gives an alternate depiction of the PSD, seeFigure 1. All of the data in (1) are candidates for being control variables. The datamay be set by external or internal control laws. Three examples of parametriccontrol are given:BURSTING: The amplitude approximation described above is | z | ≈ R ∗ ( b, λ ).If the parameter λ is slowly varying, say λ = λ ( µt ) where µ (cid:28) ε , and it ex-ecutes a cycle increasing through the value λ = 0, and back again, then R ∗ ( b, λ ( µt )) varies from low values near R = b / when λ < R = (cid:112) λ ( µt ) when λ >
0, and then back to near 0. In thatcase, the form of z shown in the previous result for ω near Ω K is z ≈ (cid:112) λ ( µt ) + e (Ω K t + φ ∗ K ) . (cid:112) λ ( µt ) + and the frequency within a burst is controlledby ω . In biological settings, λ would typically have fill-and-flush dynamics,for example, exhaustible resources are needed to sustain its activity and oncedepleted must be replenished, and in electronic circuits by a Schmitt trigger[20] for example.SELF-ORGANIZATION: If ω ( µs ) is controlled to be slowly varying throughits range, then the solution locks onto various forcing frequencies as ω ( µs )passes through intervals of phase locking. Furthermore, if W is used as afeedback control rule, say d ¯ ωds = − µ ∂W∂ω ( ¯ φ, ¯ ω ) , (15)where the right hand side is given in (11), then for ¯ φ near Ψ K and for ¯ ω nearΩ K , d ¯ ωds = − µ ˆ C K (cid:16) sin( ¯ φ − Ψ K ) / − K − ¯ ω ) cos( ¯ φ − Ψ K ) / (cid:17) where h.o.t. denotes higher order terms in ¯ ω − Ω K . This equation showsthat ¯ ω = Ω K is stable under persistent disturbances, and the frequencies tryto organize at the forcing frequencies, depending on their initialization.HYSTERESIS: If b is controlled by negative feedback from the amplitude,say ˙ b = − µR + g ( t ) , ˙ R = b + λR − R , then differentiating the second equation gives van der Pol’s equation¨ R + (3 R − λ ) ˙ R + µR = g ( t ) . Therefore, with b being controlled by negative feedback from R (now b maytake on negative values), the forced system may exhibit a wide range ofbehaviors, including negative differential resistance, hysteresis, and chaoticbehavior [17, 21]. Consider equation (1) with M = 4:˙ z = (cid:18) b z | z | + λz − | z | z (cid:19) + ωz + ε (cid:88) m =1 C m e (Ω m t +Ψ m ) . (16)10he gradient potentials for amplitudes and phases in this case are, respec-tively, V ( R, b, λ ) = R / − λR / − bR and W = − (cid:88) m =1 ˆ C m ( sinc( ν m ) cos( φ − Ψ m ) + cosc( ν m ) sin( φ − Ψ m )) (17)where ν m = (Ω m − ω ) T /ε . The output frequencies are plotted in Figure 1and W is plotted in Figure 2.Figure 1: Frequency locking: Output frequencies as a function of ω and λ . Projecting plateau regions onto the λω -plane define regions where theoutput frequency is locked, as read off the vertical axis. In this simulation C m = 10 , Ψ m = 0 , Ω m = 2 m , m ∈ , b = 0 . ε = 0 .
1, and the outputfrequency ρ = θ (10) /
10. Note that for λ <
0, subharmonics appear, asindicated by the finer stair-step structure at λ = − . (cid:126) Ω = (2 , , ,
16) for various values of ω, λ . These locking regionscorrespond to evoked memory in the system when it is illuminated by energyfrom the background field f [12]. Figure 1 shows the changes in the system’s11esponse as the parameter ω varies through an interval containing all of theforcing frequencies. The output frequencies that result from ω being in vari-ous places are shown in terms of the rotation number ρ ( ω, λ ) = lim t →∞ θ ( t ) /t .Figure 2: The potential function W in (17) is plotted using the same data asin Figure 1 and T = 1. If ω is near Ω k , then φ is captured in an energy wellcorresponding to frequency Ω k . The depth of a well is related to the powerin the corresponding forcing mode ( | C m | ). The width of a well is controlledby the parameter εC m /R ∗ . In this simulation, Ψ = [ π, π, , φ near π for ω near 2 and 4. Wells for φ near 0 mod 2 π are seen for ω near 8 and 16.The phase locking condition (12) gives an estimate of the domain of at-traction of φ ∗ K (in terms of ω ). Condition (12) is satisfied when (Ω K − ω ) R ∗ / ( C K ε ) is small and when C K is large. However, as noted after equa-tion (24), the error in that approximation holds for each R ∗ >
0, but notuniformly for R ∗ >
0. While making R ∗ near zero may make X K smaller, theerror in averaging in (24) may be large; that notwithstanding, Figure 1 indi-cates that phase locking strengthens when λ < R ∗ ≈
0. The analysisjust completed explains the structure in Figure 2 for λ >
0. Additional wrin-kles in W , those not located at the driving frequencies Ω m , are not spuriousminima, since they may correspond to various harmonic responses [21] thatare not explored further here. The simulation figures may be viewed as beinglike Bode plots for the phase response ( φ ∗ ). The results of the analysis maybe plotted in various other ways depending on the application.12 Weakly Connected STNO Networks
A weakly connected network of N oscillators based on (1) is described by theequations ˙ z n = b n z n | z n | + ( λ n + ω n ) z n − | z n | z n (18)+ εf ( t ) + ε N (cid:88) k =1 A n,k e ψ n,k z k for n ∈ , N , where f defined in (2) represents a forcing that is common to allelements in the array, each b n >
0, and, per the earlier discussion, the centerfrequencies ω n are constants. The connection matrix is A = ( A k,l exp( ψ k,l ))for k, l ∈ , N , where A k,l ≥ A k,k = 0. Again, a connection’s sign may beadjusted by adding to ψ k,l either π or 0. The connection matrix A describesinteractions between network elements, and it may represent electrical oroptical connections from one oscillator to another, or torque transfer. Thebackground signal is stronger than the connections between oscillators asreflected by the relative sizes of the small parameters ε and ε . The analysisis simplified by the fact that to order ε these oscillators are independent soeach component in the array may exhibit self-sustained oscillations, and eachhas its own potential function.Converting to polar coordinates ( z n = R n e θ n ) and taking R n = R ∗ ( b n , λ n )leads to equations for the phases˙ θ n = ω n + ε M (cid:88) m =1 ˆ C n,m sin(Ω m t + Ψ m − θ n ) (19)+ ε N (cid:88) k =1 ˆ A n,k sin( θ k − θ n + ψ n,k )where ˆ C n,m = C m /R ∗ n , and ˆ A n,k = A n,k R ∗ k /R ∗ n . The free problem ( f ≡
0) isa weakly connected Kuramoto network˙ θ n = ω n + ε N (cid:88) k =1 ˆ A n,k sin( θ k − θ n + ψ n,k )that has been extensively studied (e.g., [16]). This not a gradient systemunless A satisfies some restrictive conditions.13s before, we focus on the phase deviations φ n = θ n − ω n t :˙ φ n = ε M (cid:88) m =1 ˆ C n,m sin((Ω m − ω n ) t + Ψ m − φ n ) (20)+ ε N (cid:88) k =1 ˆ A n,k sin(( ω k − ω n ) t + ( φ k − φ n ) + ψ n,k )Averaging system (20) over T , gives˙¯ φ n = − ε ∂W∂ ¯ φ n ( ¯ φ n , ω n , ε ) + O ( ε )where W is defined in (10).Denote by S K the set of center frequencies ω n with ω n ∼ Ω K for each K ∈ , M . Then the calculation in Appendix 1 shows that θ n ≈ Ω K t + φ ∗ K . For n such that ω n / ∈ S ≡ ∪ k ∈ ,M S k d ¯ φ n ds = ε N (cid:88) k =1 ,ω k ∼ ω n ˆ A n,k sin( ¯ φ k − ¯ φ n + ψ n,k ) . This is referred to as being a network secluded from f . It is further structuredby those frequencies that match others and those that don’t. The latteroscillators wander, the former create coherent subnetworks called affinitygroups. This shows that the forcing divides the network into those oscillatorsthat are locked to a forcing mode and those that are not, but the lattermay form affinity groups of oscillators that are Kuramoto networks havingidentical center frequencies.Note that the function˜ W ( (cid:126)φ, (cid:126)ω ) = N (cid:88) n =1 W ( φ n , ω n , ε )where W is given in (10) defines a potential function for the averaged networkto leading order in ε . 14 .1 Global Feedback. Electroencephalograms (EEGs) describe oscillatory electromagnetic activityin a brain. The signals being measured are believed to be generated bysources within the neocortex [3]. One view is that EEGs represent the com-posite activity of a population of cells, which in turn creates an oscillatoryenvironment in which the network operates.To study how the STNO network operates in this case, we replace f + ε A z by N (cid:80) z n : ˙ z n = b z n | z n | + ( λ n + ω n ) z n − | z n | z n + εN N (cid:88) k =1 z k . (21)No assumption is made about the physical positioning of the oscillators; theymay be contiguous or widely dispersed. The only thing they share is the inputfrom their aggregate output.Converting to polar coordinates, setting φ n = θ − ω n t , and averaging asfor (7) gives ˙¯ φ n = ε N N (cid:88) n =1 R ∗ k R ∗ n c n,k sin( ¯ φ k − ¯ φ n ) + O ( ε ) (22)where c n,k = sinc (cid:18) T ( ω k − ω n ) ε (cid:19) . The coefficients may be symmetric since c k,l = c l,k for all k, l ∈ , N , but ingeneral R k /R l (cid:54) = R l /R k .If all R ∗ n = R ∗ are identical, a standard argument, e.g. Theorem 11.4 [16],shows that for each n ∈ , N the phase deviation ¯ φ n has a limit, say ¯ φ n → ¯ φ ∗ n .Suppose the center frequencies have a common random distribution, say foreach n ∈ , N , ω n ∈ N (Ω , σ ), which is a normal gaussian random variablehaving mean Ω and variance σ . Then the coefficients c n,k = sinc (cid:18) T ( ω k − ω n ) ε (cid:19) form a symmetric matrix, and the arguments of sinc are identically dis-tributed random variables having a normal distribution. Averaging over N gives the averaged coefficients¯ c = (cid:90) ∞−∞ sinc( η ) dw ( η )15here w is the probability measure for N [24]. Setting ξ n = (Ω − ω n ) t − ¯ φ n in (22) and averaging over N and over T /ε gives˙ ξ n = ε N ¯ c N (cid:88) k =1 sin( ξ k − ξ n ) . A standard argument shows that the variables ξ n have limitslim t →∞ ξ n ( t ) = χ ∗ n , and θ n ≈ ¯ ω t + χ ∗ n , where ¯ ω is the average of the set { ω n } . The approxima-tion is in a probabilistic sense not made precise here, except to say that theapproximation degrades with increasing variance T σ /ε . In this case theglobal feedback synchronizes the oscillators with the phase deviation distri-bution χ ∗ and the output frequencies are ¯ ω . A simulation of this is shownin Figure 3 where both the center frequencies and the amplitudes R ∗ n arerandomly chosen.Figure 3 shows that global feedback within the network may stabilizeit at a common frequency, and small random perturbations of the data donot significantly impact this result. If there are additional oscillators thathave their center frequencies similarly clustered but sufficiently separated,the feedback will stabilize the cluster to a single frequency and the PSD of z will have spikes at the locked frequencies. No assumption is made aboutthe spatial location of oscillators; but, the results here show coherence in thefrequency domain when the composite signal feedbacks to the network.There are many more interesting problems here for stochastic analysis,but we consider only this case to illustrate frequency and phase-locking fea-tures of the global feedback network. In particular, global feedback maysynchronize the population of phase deviations and lock the oscillator fre-quencies. The aims of this paper were to demonstrate how the network (18) can createand sustain memory, process and classify information, and perform compu-tations.
Memory:
The valleys in a potential may be used as memory. Thesemay be accessed by associative recall since valleys correspond to basins of16igure 3: This simulation plots the input and output frequencies of 200oscillators when the center frequencies form a gradient ω n = πn/ η n for η n ∈ N (0 , . n ∈ , , and the amplitudes are randomly chosen: R ∗ n ∈ U , which is the uniform distribution on (0 , ω n ,triangle markers) and output frequencies ( ρ n , square markers) are plotted for200 oscillators. Equation (22) is solved for 0 ≤ t ≤ ε = 0 . , λ = 1 . , and T = 1 .
0. This demonstrates that global feedback stabilizes the outputfrequency of the population. In most cases, the output frequencies convergeto the average of the center frequencies, which is ¯ ω = 0 . ρ = ¯ ω . The deviationsof ρ n from ¯ ω are largely due to random variations in the ratios R ∗ n /R ∗ k .attraction of their minima. That is, if the system is initialized in a valley,it will reliably converge toward a corresponding minimum of the potential.The analysis here identified the phase potential function, located its basinsof attraction, and estimated their scope.There are two kinds of memory in system (1), one called here structural,or time domain, memory, and the other evoked, or frequency domain, mem-ory. The function V with b > , λ >
0, has a single positive minimum at R ∗ ≈ √ λ , although the analysis here supports V having more wells. Evokedmemory is described by the phase potential. Phase states are stabilized byexternal forcing, somewhat similar to the up position of a pendulum beingstabilized by vibration of its support point, but which becomes unstable when17he vibration is removed [12]. The phase potential W represents M evokedmemory units per oscillator. In the network there are possibly M × N evokedmemory units, depending on the distribution of center frequencies. Information flow:
The forcing frequencies in f define channels for trans-ferring information using phase deviations. The oscillators act as phase-locked loops [22] in that their phases converge to forcing phases. If infor-mation is encoded in the forcing phases, say Ψ m ( t ), as in phase modulationsignaling, the network will track these phases in a stable way, and the resultsmay be transmitted by the outputs { z n } . The power spectrum density (PSD)of { z n } for n ∈ , N can be used to describe the channels for communicationand their robustness. Computation:
In-phase and anti-phase oscillations in ONN may be usedfor a variety of computations. Phase deviations may carry digital informationin oscillations. For example, the bit ‘0’ corresponds to Ψ = 0 and the bit‘1’ to the anti-phase oscillation where Ψ = π . z carries this informationforward in the various channels established by the forcing frequencies, andthe information may be retrieved from z by correlations [13]. Free system:
If the forcing is not present ( f ≡ Random noise perturbations:
Stability under persistent disturbancesas used here also applies to systems perturbed by small random noise. Largedeviation random perturbations may be studied using the theory developedin [23, 24] for random dynamics in gradient systems. However, with randomnoise comes the possibility that R = 0 even if b >
0, so the solution may bedriven through the essential singularity at z = 0, and so it may experienceunknown phase shifts, or cycle slips, which might be near π or not. Continuum model:
For a planar array of STNO each magnetic ele-ment is described by its magnetization, taken here to be the complex-valuedfunction z ( x, t ) of time t and a spatial variable x ∈ E . The LandauLif-shitzGilbertSlonczewski equation (see [14, 15]) describes precessional motionof magnetization in a layer, and it accounts for spin-transfer torque fromone oscillator to another. The analog of this in our context is the nonlinear18chr¨odinger equation [26] ˙ z = (1 + β ) D ∇ z − F ( z ) + G ( | z | ) z + εf ( t ) . (23)where D accounts for torque transfer between elements, F ( z ) = bz/ | z | + λz − | z | z , G is a smooth function satisfying G (0) >
0, and β is a dampingcoefficient. Direct calculation shows that this problem, to leading order in ε ,has spin waves of the form z = R ∗ e ( (cid:126)k · x +(1+ β ) D | (cid:126)k | t ) e ω ∗ t having wave vector (cid:126)k and frequency ω ∗ = G ( R ∗ ). Spin waves and vortexsolutions have been found for this equation when b = 0 , ε = 0, and excitationsin such media may propagate as spin waves that may carry information forcomputations [13, 15, 25]. A second centered-difference discretization schemefor (23) will have the form of (18). The response to external forcing by a single STNO with phase factor wasstudied first. It was shown that there are potential functions for both theamplitudes and phases. The phase potential was illustrated in the simulationsin Figures 1 and 2. The averaged potential function W ∞ gives a usefulapproximation to W , and it provides a convenient description of the locationof energy wells and of the distribution of in-phase/anti-phase behavior whenthe forcing phase deviations are 0 or π , respectively. If one defines the energyof an oscillator to be proportional to its frequency, the potential function W in (9) represents a landscape of energy for the forced oscillator. Theanalysis here shows how, given the Fourier decomposition of the forcing f ,there are potential functions, V for amplitude and W for phase deviationsfor weakly forced networks, as well. The data in V , W and f shape theresponses of the system, and these data may be used as control variables toplace equilibria and to shape their domains of stability, both in amplitudeand phase. These results facilitate using ONN for training and learning, forexample by showing how to manipulate outputs to produce patterns thatagree with target outcomes. 19 eferences [1] G. Buzsaki (2006) Rhythms of the Brain, Oxford U Press.[2] R. Guttman, et al. (1980) Frequency Entrainment of Squid Axon Mem-brane, J Membr Physiol 56:9-18.[3] P. L. Nunez, R. Srinivasan (2005) Electric Fields of the Brain: TheNeurophysics of EEG, 2nd Ed., Oxford U Press.[4] E. M. Izhikevich (2001) Resonate-and-Fire Neurons, Neural Networks,14:883-894.[5] E. M. Izhikevich, et al., (2003) Bursts as a unit of neural information:selective communication via resonance. Trends in Neurosci. 26 (3):161-167.[6] J. O’Keefe, M. L. Recce (1993) Phase Relationship Between Hippocam-pal Place Units and the EEG Theta Rhythm. Hippocampus, 3 (3): 317-330.[7] P. Fries, D. Nikoli, W. Singer (2007) The Gamma Cycle, Trends inneuroscience, 30 (7): 309-16.[8] A. Winfree, (1980), The Geometry of Biological Time, Springer-Verlag,New York.[9] J. Torrejon, et al. (2017) Neuromorphic computing with nanoscale spin-tronic oscillators, Nature, 547:428-431 (27 July 2017)[10] M. Romero, et al., (2018) Vowel recognition with four cou-pled spin-torque nano-oscillators, Nature, v 563 (8 Nov, 2018).(doi.org/10.1038/s41586-018-0632-y).[11] F. Hoppensteadt, E. Izhikevich (1999) Oscillatory Neurocomputers withDynamic Connectivity, PRL, vol. 82, Consider (6) with ω ∗ being near to Ω K for some K ∈ , M . The equationmay be rewritten as˙ φ = ε ˆ C K sin((Ω K − ω ∗ ) t + Ψ K − φ ) + ε M (cid:88) m =1 ,m (cid:54) = K ˆ C m sin((Ω m − ω ∗ ) t + Ψ m − φ ) . Let a new variable ξ be defined by ξ = (Ω K − ω ∗ ) t + Ψ K − φ . This gives onthe slow time scale dξds = Ω K − ω ∗ ε − ˆ C K sin ξ − M (cid:88) m =1 ,m (cid:54) = K ˆ C m sin (cid:16) (Ω m − Ω K ) sε + Ψ m − Ψ K + ξ (cid:17) . (24)The last terms are higher frequency terms since | Ω m − Ω K | (cid:29) ε for m (cid:54) = K .The method of averaging [19] is applied to this equation over 0 ≤ s ≤ T .The higher frequency terms average to O ( ε/ ( T R ∗ )); this estimate holds foreach R ∗ >
0, but not uniformly in R ∗ . This gives to leading order in εd ¯ ξds = Ω K − ω ∗ ε − ˆ C K sin ¯ ξ, (25)and the method shows that ξ ( t ) − ¯ ξ ( εt ) = O ( ε ) uniformly for 0 ≤ s ≤ T . Equation (25) may also be written as a gradient system, and it has anequilibrium ¯ ξ = arcsin X K where X K ( ω ∗ ) = (Ω K − ω ∗ ) R ∗ εC K , if the phase locking condition (12) holds. The equilibrium is stable underpersistent disturbances.As a result, for all data satisfying (12), θ = ω ∗ t + φ = ωt − ξ + (Ω K − ω ∗ ) t + Ψ K ≈ Ω K t + φ ∗ K φ ∗ K = Ψ K − arcsin X K . It is shown in [17] that ξ = φ ∗ K + O (exp( − ˆ C K s )),and so after an initial transient, ¯ ξ equilibrates to φ ∗ .The phase deviation φ is not constant, but has frequency Ω K − ω ∗ . Therelation between ˙ φ and ω ∗ − Ω K is essentially linear for ω ∗ near Ω K ; this factfacilitates the design of phase locking devices [27]. z . The power spectrum density (PSD) of f given in (2) is P SD f ( ρ ) = M (cid:88) m =1 | C m | δ ( ρ − Ω m ) . The approximation to z derived earlier, may be used to approximate the PSDof z using the formula P SD z ( ρ ) = |F z ( ρ ) | where F z is the Fourier transform of z given by F z ( ρ ) = εT (cid:90) T/ε e − ρt z ( t ) dt ≈ εT (cid:90) T/ε e − ρt R ∗ e θ ( t ) dt. The notation ω ∗ ∼ Ω K means that the phase locking condition (12) is satis-fied, so θ ≈ Ω K t + φ ∗ K ( ω ∗ ); and if ω ∗ ∼ Ω K , then F z ( ρ ) ≈ εT (cid:90) T/ε e − ρt R ∗ e (Ω K t + φ ∗ K ) dt = R ∗ e φ ∗ K sinc((Ω K − ρ ) T /ε ) dt + O ( ε/T ) . The result is that
P SD z ( ω ∗ ) ≈ R ∗ sinc ((Ω K − ω ∗ ) T /ε ) + O ( ε/T ) ..