Thermally-robust spatiotemporal parallel reservoir computing by frequency filtering in frustrated magnets

Physical reservoir computing is a framework for brain-inspired information processing that utilizes nonlinear and high-dimensional dynamics in non-von-Neumann systems. In recent years, spintronic devices have been proposed for use as physical reservoirs, but their practical application remains a major challenge, mainly because thermal noise prevents them from retaining short-term memory, the essence of neuromorphic computing. Here, we propose a framework for spintronic physical reservoirs that exploits frequency domain dynamics in interacting spins. Through the effective use of frequency filters, we demonstrate, for a model of frustrated magnets, both robustness to thermal fluctuations and feasibility of frequency division multiplexing. This scheme can be coupled with parallelization in spatial domain even down to the level of a single spin, yielding a vast number of spatiotemporal computational units. Furthermore, the nonlinearity via the exchange interaction allows information processing among different frequency threads. Our findings establish a design principle for high-performance spintronic reservoirs with the potential for highly integrated devices.


Introduction
Physical reservoir computing is attracting interdisciplinary attentions as one of the key technologies for realizing real-time information processing required in the coming Internet of Things society, while breaking free from energy-consuming silicon-based devices (1)(2)(3).
Taking advantage of nonlinear phenomena in a physical system known as a "physical reservoir", the input data is mapped onto a high-dimensional feature space in a nonlinear manner.In addition to these two characteristics, fading memory property is considered another important feature for processing time-varying input streams (4,5).In this framework, internal states of the reservoir are represented by a read-out function in the form of vectors, and a linear transformation of these vectors provides the final output for machine learning problems, such as classification and predictive analytics.In contrast to usual neural networks, training is limited to linear regression only at the read-out layer, which enables physical reservoirs to operate at low power, high speed, and high versatility, as demonstrated in various implementations, including photonic systems (6)(7)(8)(9)(10)(11)(12)(13)(14) and electrical systems (15)(16)(17).
Toward device applications of spintronic reservoirs in realistic systems, two major challenges remain to be resolved.One is robustness against thermal fluctuations, which is crucial for retaining short-term memory (STM) long enough for practical use.Most spintronic reservoirs exploit nonlinear phenomena originating from magnetization dynamics, and are therefore inherently vulnerable to external noises that disturb spin precessions.The other challenge is designing devices with highly integrated computational units, the importance of which is evident from the problems facing silicon chips.The straightforward solution for the latter is the use of multiple read-outs (19,20,28,31,32), in which multiple internal state vectors are extracted from several local measurements, such as the dynamics of particular spins, and each is used independently for different computations.Another solution is parallelization by a selective read-out from superposed signals, such as wavelength-multiplexing (10,13) and frequency-multiplexing (14) studied in photonic systems, but such schemes have been largely unexplored for spintronic reservoirs thus far.
In this article, we propose a framework of spintronic physical reservoir computing that achieves both thermal-robustness and high-integration.We demonstrate it by utilizing a simple model of frustrated magnets as a physical reservoir, with the input by AC magnetic fields and read-out by spin dynamics.We find that the memories of input information are stored in the spin dynamics at frequencies around that of the AC field and can be preserved against thermal noise by filtering out irrelevant signals at the other frequencies.We also demonstrate parallel processing on multiple inputs in both space and time, using read-outs on different spins at different frequencies.Moreover, we show that the nonlinearity arising from the exchange interaction allows calculations that capture information across various frequencies without the need for dedicated communication channels.Our results pave the way for spintronic physical reservoir computing in realistic situations.

Physical reservoir and input/output
Our spintronic reservoir receives a random sequence of input binary digits { !} in the form of local magnetic fields, and transforms them nonlinearly through the spin dynamics (Fig. 1a).As a prototypical model of frustrated magnets, we consider an antiferromagnetic Heisenberg model with classical spins on an  ×  ( = 128) triangular lattice with open boundary conditions.The Hamiltonian is given by where  " = 5 " * ,  " + ,  " ' 7 represents the classical spin at site  (| " | = 1) and the sum of ⟨, ⟩ is taken for nearest-neighbor sites;  # ' () is the time-dependent magnetic field for the input (Fig. 1b) and Λ is a set of lattice sites used as the input terminals for the physical reservoir (Fig. 1c).Hereafter we take  = 1 as the energy unit.In the absence of magnetic fields, the spin configuration of this system has a 120 ∘ structure with three sublattices at low temperature (33,34), which might be realized experimentally, for instance, in ACrO2 (A = Li, Na) (35)(36)(37)(38), Ba3CoSb2O9 (39,40), and κ-(BEDT-TTF)2Cu2(CN)3 (41,42).
The input is a time series information of the binary bit  ! with an interval of time, which we take  -. = 12, and converted into the magnetic field as where  -. and  -. is the norm and frequency of the magnetic field, respectively (Fig. 1b).
We take  -. = /(2 -. ) with an integer  so that  # ' () varies continuously in time.The magnetization dynamics under the input field is simulated by the stochastic Landau-Lifshitz-Gilbert (LLG) equation where  is the Gilbert damping constant and  " 011 is the effective magnetic field at site  consisting of the input magnetic field, the exchange magnetic field, and the thermal field at temperature  (43, 44) (see the Methods section).We take  = 0.1 in the following calculations.The time profiles of the  components of the read-out spins  # ' are observed every Δ -. =  -. / 2 (Fig. 1d) where we take  2 = 120 (45).Then, the internal states of the reservoir,  !, are defined by an ( 2 + 1) -dimensional vector including an additional component with the constant value 1 as The final output is obtained by linear transformation of the internal state vector  !=  !⋅  where  is a weight vector.With sufficiently large training data, the weight  is trained so that each  !becomes close to the desired output  Z ! for a given problem.The reservoir performance is evaluated by the determination coefficient  / between the sequence of outputs  = { !} and targets  Z = { Z !}:  / = cov / (,  Z)/[ / () / ( Z)] where cov and σ 2 denote covariance and variance, respectively. / is close to 1 when  and  Z are well matched, and approaches 0 otherwise (see the Methods section).

Thermal robustness by frequency domain filtering
First, we examine the STM task to evaluate how long the input information is retained in the reservoir.The target output for this task is  Z !=  !34 where  is the delay step after the input.Here, we apply the input magnetic field at frequency  -. = 8/ -. to one specific terminal spin at site , namely, Λ = {  }, and use the spin dynamics at the same site,  # ' (), as a read-out to construct the internal state vectors.We place  around the center of the system to avoid uninteresting behavior near the edges.
Figure 2a shows temperature dependence of the reservoir performance for the STM task.
At absolute zero temperature, the memories are recovered with high accuracy of  / > 0.95 up to delay  = 7, and decay gradually for larger , exhibiting the fading memory property required for a physical reservoir.In contrast, at finite temperature, the reservoir rapidly loses the STM older than  = 2 even at extremely low temperature  = 10 35 , and furthermore, remaining memories of  = 0 and 1 are also lost as the temperature increases.This is an indication that the STM stored in the spin dynamics is extremely fragile against thermal noise, as commonly seen in various spintronic physical reservoirs (18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28).
To clarify how the thermal noise disturbs the read-out signals, we analyze the Fourier spectrum of the dynamics of  # ' () .Figure 2d displays the spectra at three different temperatures.At  = 0, the spectrum shows strong peaks at around the frequency of the input magnetic field, ± -. , as the terminal spin linearly follows the time-dependent magnetic field via the Zeeman coupling.The wavy structure repeating every 1/(2 -. ) also appears, which originates from the random switching of the input magnetic field at every  -.
according to the input bit.Here, the finite peak width at ± -. reflects the probabilistic nature of the bit sequence { !} in the input magnetic field, otherwise there would be only a deltafunctional peak at exactly ± -. .In contrast, at finite temperature, thermal noise disturbs the spin dynamics, obscuring these characteristic peaks.At  = 0.001, the wavy structure disappears due to the thermal noise, but the peaks at around ± -. are still present.The latter peaks are also mostly drowned out at a higher temperature,  = 0.1.
From this observation, we introduce a frequency filter in order to mitigate the disturbances caused by thermal agitations.Specifically, we only retain signals within the frequency windows of  -. ≤ || <  -. + 1/(2 -. ) , and use the frequency-filtered spin dynamics  #,1-670809 ' () for computation instead of the whole spin dynamics,  # ' (); see Fig. 2c and the Methods section.Figure 2b shows the reservoir performance for the STM task with the frequency filter.In clear contrast to the case without frequency filtering in Fig. 2a, the STM persists up to  = 5 with almost no information loss even at finite temperature, and the older memories are lost gradually.This indicates, surprisingly, that the input information can be recovered only by the signals within the narrow frequency bandwidth of 1/(2 -. ) including  -. , and that the memories within this range are resistant to thermal fluctuations.
Figure 2e shows the distributions of the STM in frequency domain.We introduce frequency filters that pass signals within the frequency window of ( − 1)/(2 -. ) ≤ || < /(2 -. ), and evaluate the performance for the STM task while changing a positive integer .At  = 0, the memories are encoded within a broad frequency range, but  / gradually diminishes as the frequency window deviates from  -. .Since the information embedding process from the input magnetic field is mainly governed by the linear Zeeman interaction, the input information is prominently restored in the spin dynamics at around  -. .At  = 0.001, the system performance at around  -. remains largely unchanged from that at  = 0, although it deteriorates significantly at other frequencies.A comparison with Fig. 2d demonstrates that the input information is retained only at frequencies where the spin dynamics is preserved against thermal noise.As the signals at the other frequencies contribute as noises, the reservoir without the frequency filter is vulnerable to thermal agitations.At  = 0.1, the reservoir is almost non-functional even at around  -. because of the increased noise level as shown in Fig. 2b.In our configuration, we employ a single spin as the input and read-out, however, we verify that the utilization of multiple spins or the entire system does not alter the fundamental behavior.We note that increasing the amplitude of the input magnetic field broadens the temperature range within which the system can operate effectively (see Supplementary Note 1 and Supplementary Fig. 1).
The influence of the input information in frequency domain is closely tied to the typical time scale of the system.In our input protocol, the input magnetic field is decomposed into the product of an AC magnetic field at frequency  -. and a pulse square field that randomly changes its sign every  -. depending on the supplied bit.This switching duration is the only relevant time scale for the randomness of the input, and thus the information is stored in units of 1/ -. in frequency domain.In other words, the memories retained at around ,  ± / -. and − + / -. with a positive interger  are almost identical (see Supplementary Note 2 and Supplementary Fig. 2).This observation suggests that filters of ( − 1)/ (2 -. ) ≤ || < /(2 -. ) are the narrowest frequency bandwidth that enables maximum recovery of STM without overlapping information.

Parallel reservoir computing in frequency and spatial domains
The above concept of frequency filtering can be extended to implement frequency division multiplexing, which allows for the simultaneous transmission of multiple signals encoded at different frequencies.As an illustration, we prepare three sequences of random input bits { !: }, { !/ }, { ! ; } and examine the STM task of each input bit: the target output is with delay step  and  = 1, 2, 3. We choose different frequencies  -. : = 8/ -. ,  -. / = 45/(2 -. ),  -. ; = 36/ -. for each binary sequence, and transform the bits  !< to the magnetic field at the assigned frequency in the same manner as Eq.2: superposition of these three magnetic fields is given to one input terminal spin at site  (Fig. 3a).Herein the input bits  !< are nonlinearly transformed through the dynamics of the terminal spin  # ' () and the internal state vectors  ! are then extracted for computations.Figure 3b shows the Fourier spectrum of  # ' () at zero temperature, where three peaks corresponding to ± -. : , ± -. / , ± -.
; appear in addition to the wavy structure present over a wide frequency range.
Figure 3c represents the reservoir performance for the STM tasks of three input bit sequences in frequency domain.Each information is embedded into the spin dynamics at around its input frequency and not affected by the other input bits.Consequently, the memories can be selectively recovered using a frequency filter that passes the input frequency where the desired bits are provided.Indeed, for all ,  / at around each  -. < is greater than 0.9 up to delay  = 5 and gradually decreases for larger  due to the fading memory property, which indicates that the performance is almost the same for different parallel threads.This clearly demonstrates that the spin dynamics at different frequencies can be regarded as mutually independent parallel threads, allowing for frequency division multiplexed computational units.This parallelization scheme works at finite temperature in the same manner as Fig. 2e by filtering out irrelevant signals.Worth noting that the degree of parallelization can be further increased by allocating more frequencies for computation; we confirm that a spacing of 1/ -. or more is sufficient for such parallelization.
Our reservoir can also be parallelized in spatial domain with multiple inputs and readouts.To demonstrate this, we perform the STM task by selecting two input terminal spins at sites  : and  / separated by 17 spins, namely Λ = {  : ,  / } in Eq. 1, denoted as terminal 1 and terminal 2, respectively (Fig. 4a).A bit sequence of m !# ! ,: n 5m !# ! ,/ n7 is given to terminal 1 and m !# " ,: n 5m !# " ,/ n7 is given to terminal 2 at frequency  -. : ( -. / ); in total, four input bit sequences are simultaneously supplied to the system.Figure 4b shows the performance for the STM tasks of different input sequences at zero temperature.For computation, the signals received from terminal 1 and terminal 2 are independently utilized with frequency filters.
When the spin dynamics at terminal 1 is used as the read-out, the memories of m !# ! ,: n and m !# ! ,/ n can be reproduced at around each input frequency with  / > 0.9 up to delay  = 5, whereas the information of m !# " ,: n and m !# " ,/ n cannot be recovered at all.The opposite is true when utilizing signals from terminal 2. This illustrates that the spin dynamics at different sites can be treated as independent parallel computational threads, which realizes spatial multiplexing of computational units while preserving the frequency multiplexing.
Due to the availability of single spins, the minimal components of magnetic materials, as read-outs, our proposed scheme offers the potential for high-density integration at the nanoscale on one magnetic chip; we confirm that a lattice spacing more than 10 spins is sufficient for such parallel computation in our model.Moreover, simultaneous parallelization in spatial and frequency domains would allow for extreme integration of computational units in spatiotemporal domain by further increasing the number of frequency bands and I/O terminal spins.

Logic gate operation of bits in different threads
Applying ), respectively.Here, we observe the spin dynamics on the nearest neighbor site from the input terminal site, denoted as  == , for the read-out from the reservoir; the spin dynamics at the terminal site  is also traced for a comparison (Fig. 5a).
The target output  Z ! is the bit where two bits  !: and  !/ are processed by a logic gate, the truth table of which is summarized in Fig. 5b. Figure 5c represents the Fourier spectra of the spin dynamics on the terminal site  and the neighboring site  == at zero temperature.While the overall behaviors, especially strong peaks at ± -.
: and ± -. / , are common between the two spectra, a significant difference is found in additional peaks in the spectrum from the site  == at around the second harmonics frequencies: ±2 -. ).With read-out on the nearest neighbor site, the reservoir operates very well with  / > 0.9 for both linearly separable and inseparable tasks.These findings clearly indicate that the nonlinear processes via the exchange interaction enable two bits in different threads to exert an influence on one another.The importance of this coupling is evident from the very low  / when utilizing the input spin as the read-out in Fig. 5e, where the dynamics is dominated by the linear Zeeman coupling.Multi-bit tasks such as gating require mixing different information within the same reservoir dynamics, otherwise different bits cannot establish a nonlinear relationship.This is thus an operation that takes full advantage of the principle of superposition in the frequency domain, and is almost impossible with other parallelization schemes such as spatial domain multiplexing.We note, however, that the second harmonics in the dynamics of the nearest neighbor spin are vulnerable to thermal noise, making linearly inseparable gating almost infeasible at finite temperature.This issue would be resolved by enhancing nonlinearity in models with different interactions, as we shall discuss later.

Discussion
In this paper, we have proposed a framework for a physical reservoir computing in magnetic materials utilizing frequency domain dynamics.The input through an AC magnetic field and the selective read-out with frequency filters enable both thermal noise reduction and frequency division multiplexing, which are of great practical significance in the design of spintronic physical reservoir computing devices.In addition, by utilizing multiple read-outs, a substantial number of computational units can be integrated through parallelization in spatiotemporal domain.The inherent nonlinearity in the interacting spin systems allows mixing of information stored in the spin dynamics at different frequencies.
In magnetic materials, the linearly introduced information via the Zeeman coupling is nonlinearly processed by the exchange interaction.This is where the differences in materials are most pronounced.For instance, an increase in connectivity and complexity of interactions would enhance the nonlinearity of the spin dynamics, and longer-ranged interactions would propagate the input information over a wider region.When electrons in the system exhibit itinerant nature rather than localized one, the synergy between charge and spin degrees of freedom might alter the information transfer and processing.Moreover, in quantum spin systems, the reservoir takes advantage of the large dimensionality of Hilbert space, where quantum correlations and quantum entanglement would strongly influence the performance of the reservoir (31,32,46,47).The impact of specific interaction on information processing will be a subject of future research.
The complexity of the magnetic structure contributes to the high dimensionality of the feature space where information is projected by a physical reservoir, providing ample resources for machine learning tasks.Here, the dimension corresponds to the number of spins whose dynamics is linearly independent, which determines the number of parameters in the feature space.In general, the spin dynamics in a uniform environment, such as a ferromagnetic state, can be described with a few parameters, while those in a complex environment require many parameters.Our reservoir with a frustrated magnet, consisting of only three sublattices, yields low dimensional internal states, yet proves sufficient for the tasks demonstrated; we confirm magnets with solely nearest-neighbor ferromagnetic interactions are too simple to be employed for computation.Physical reservoirs with more complex magnetic textures would project information onto higher dimensional feature space, which would enhance computational performance for complicated tasks, as presented in skyrmion systems (27,28).
In addition to its thermal robustness and parallel processing capabilities, our scheme exhibits versatility in its broad applicability to a variety of spintronic physical reservoirs with input by AC signals and read-out by time-varying signals.In magnetoelectric systems, for instance, the information might be embedded and extracted via either magnetic or electrical AC signals, while taking advantage of the nonlinearity of the spin-charge coupled dynamics for information processing (48,49).The input information should be stored within the spin and charge dynamics at the characteristic frequency bands, where frequency filtering would play an important role.Consequently, our implementation-independent protocol would be beneficial for realizing thermal robustness and frequency division multiplexing in various types of spintronic physical reservoirs.Our discovery opens up the way for hardware implementation of large-scale intelligent systems that exploits magnetization dynamics for real-time computation.

LLG calculation
To simulate the magnetization dynamics, we employ the stochastic LLG equation summarized in Eq. 3. The effective magnetic field includes the exchange field  " 0> and the thermal field  " 7? .The former is given by the  " derivative of the Hamiltonian (Eq. 1) as  " 0> = −ℋ/ " , and the latter is introduced as a white noise with the following properties: where Δ denotes the time step (we take Δ=0.005),() the Dirac delta function, and  "# the Kronecker delta function.The time evolution is traced by the fourth-order Runge-Kutta method.At each temperature, equilibrium spin configurations are obtained by relaxation for a sufficient time starting from a complete in-plane 120 ∘ structure, and they are used as an initial state of the reservoir.The temperature dependences of the energy and specific heat are in agreement with the previous results using Monte Carlo calculations (33,50).

Scheme for training and testing
We take sequences of binary digits { !}, and each data  ! is fed to the system by the out-of-plane input magnetic field to the input terminal spin.The internal state is characterized by the induced out-of-plane dynamics at the read-out spin as the 1 × ( 2 + 1)dimensional vector  !(Eq.4) that is then linearly transformed to the output  ! by the and an  7C × ( 2 + 1) matrix , and correspondingly the target output vectors are given by an and an . The weight should be trained to satisfy  Z 78 =  78 ⋅ , which is impossible for overdetermined systems with  78 ≫  2 + 1.
Alternatively, we optimize  to minimize the least squares error between vectors on both sides, which is analytically given by  =  78 ' ⋅  Z 78 where  F is the Moore-Penrose pseudoinverse-matrix.The output in the testing phase is calculated by  =  7C ⋅  and the similarity to the desired output  Z 7C is evaluated using the determination coefficient

Frequency filtering
The spin dynamics in the training phase  78 () and in the testing phase  7C () are

Input
Terminal node

Physical reservoir Output
Read-out node every Δ -. =  -. / 2 and transformed to an ( 2 + 1)-dimensional internal state vector  ! ; see Eq. 4. The final output  ! is calculated as  !=  !⋅  with the weight vector .Input    The same plot as d with frequency filters centered on ±( -. / −  -. : ) and ±( -. : +  -. / ).The performance of the nearest neighbor spin becomes much better for all gating operations including the inseparable gates, owing to the nonlinear processes via the exchange interaction.

Supplementary
2 + 1) × 1-dimensional weight .In total  7A7 = 6000 inputs are used:  B = 2000 to warm up the system,  78 = 2000 for training of , and  7C = 2000 for testing the reservoir performance.Since they are sufficiently large numbers compared to the number of tuning parameters  2 + 1 ( 2 = 120) , the generalization performance is preserved without overfitting.The vectors { !} at training phase and test phase are summarized to an

Fig. 2 :a 6
Fig. 2: Reservoir computing in the Heisenberg antiferromagnet at finite temperature.a Reservoir performance for the STM task at several temperatures  when using the whole spin dynamics (without frequency filtering).Delay step  denotes discretized time length after input and  / is the determination coefficient.The STM is extremely fragile against thermal noise, especially for  ≥ 2. b The same plot as a with the frequency filter of  -. ≤

Fig. 3 :Fig. 4 :
Fig. 3: Parallel reservoir computing by frequency division multiplexing with the superposed input magnetic fields.a Schematic of the preparation of input magnetic field.Multiple bit sequences are converted to magnetic fields at different frequencies, and their superposition is given as the input magnetic field.b Fourier spectrum of the spin dynamics the STM of m !# ! ,: n and m !# ! ,/ n, and the bottom two figures show the STM of m !# " ,: n and m !# " ,/ n.The gray lines correspond to  -.: and  -. / .The STM of m !# ! ,: n and m !# ! ,/ n can be recovered only by the signals from the terminal 1, and the STM of m !# " ,: n and m !# " ,/ n requires the signals from the terminal 2 for reproduction.Each memory is retained at around the corresponding input frequency.A straightforward extension of this scheme would allow for high-density integration of computational units in both spatial and frequency domains.

Fig. 5 :
Fig. 5: Logic gate operations of two input bits stored in different frequencies.a Schematic of the input terminal spin (the blue arrow) and its nearest neighbor spin (the red arrow).The former is directly affected by the input magnetic field through a probe represented by the yellow needle, while the latter is affected via the exchange interaction between the spins.b Truth table for logic gates with input bits  !: and  !/ .c Fourier spectra gates, XOR and XNOR, are nearly impossible with these frequency filters.e
(9,29)aracteristics of frequency division multiplexing, different information retained in parallel at different frequencies can be accessed simultaneously by selecting signals at appropriate frequencies.To demonstrate this, we examine logic gate tasks between input bits in different parallel threads(9,29).Two binary sequences { !: } and { !/ } are simultaneously given to one input terminal spin at site  encoded by the frequency  -.
: = 8/ -. and  -. / = 45/(2 -. , which originate from the square of magnetic fields at the same frequency, are delta-functional since the sign information of binary bits is 2 -. ) and  -./ − 1/(2 -. ) ≤ || <  -. / + 1/(2 -. ) .Figure5dplots the reservoir performance for each task, where the blue and red bars exhibit  / with read-outs on the input site  and the nearest neighbor site  == .By utilizing signals from the input site, the reservoir shows moderate performance with  / ∼ 0.6 for linearly separable tasks, namely AND, OR, NAND, and NOR; conversely, the performance for linearly inseparable tasks, XOR and XNOR, is extremely poor with  / ∼ 0. When utilizing the nearest neighbor spin for the read-out, the performance becomes slightly better, but  / for the inseparable tasks are still low.Thus, although the signals containing input frequencies of ± -.
in both spectra, while only that of the nearest neighbor spin shows peaks at ±2 -.Reservoir performance for the logic gate tasks with frequency filters centered on ± -.