# Andreev reflection of quantum Hall states through a quantum point contact

Mehdi Hatefipour1 Joseph J. Cuozzo2 Ido Levy1 William M. Strickland1 Dylan
Langone1 Enrico Rossi3 Javad Shabani1<EMAIL_ADDRESS>1Center for Quantum
Information Physics, Department of Physics, New York University, NY 10003, USA
2Materials Physics Department, Sandia National Laboratories, Livermore, CA
94551, USA
3 Department of Physics, William & Mary, Williamsburg, Virginia 23187, USA

###### Abstract

We investigate the interplay between the quantum Hall (QH) effect and
superconductivity in InAs surface quantum well (SQW)/NbTiN heterostructures
using a quantum point contact (QPC). We use QPC to control the proximity of
the edge states to the superconductor. By measuring the upstream and
downstream resistances of the device, we investigate the efficiency of Andreev
conversion at the InAs/NbTiN interface. Our experimental data is analyzed
using the Landauer-Büttiker formalism, generalized to allow for Andreev
reflection processes. We show that by varying the voltage of the QPC,
$V_{QPC}$, the average Andreev reflection, $A$, at the QH-SC interface can be
tuned from 50% to $\sim 10$%. The evolution of $A$ with $V_{QPC}$ extracted
from the measurements exhibits plateaus separated by regions for which $A$
varies continuously with $V_{QPC}$. The presence of plateaus suggests that for
some ranges of $V_{QPC}$ the QPC might be pinching off almost completely from
the QH-SC interface some of the edge modes. Our work shows a new experimental
setup to control and advance the understanding of the complex interplay
between superconductivity and QH effect in two-dimensional electron gas
systems.

##  I I. Introduction

In recent years, the realization of topological phases of matter has been a
focus of intense research in the field of condensed matter physics Kitaev
(2001); Beenakker (2013); Alicea (2012). This is partly motivated by the
potential of such systems to host exotic quasiparticles with non-Abelian
statistics, which could be used for fault-tolerant quantum computation through
braiding operations Nayak et al. (2008); Laucht et al. (2016); Alicea (2012).
While InAs surface quantum wells (SQWs) in proximity to superconductors have
emerged as a promising platform for the realization of topological
superconductivity Aghaee et al. (2022); Dartiailh et al. (2021); Lee et al.
(2017), it is proposed that combining rich quantum Hall physics with
superconductivity can allow access to higher topological states Mong et al.
(2014); Clarke et al. (2013).

Recent studies Hatefipour et al. (2022); Gül et al. (2022); Lee et al. (2017);
Zhao et al. (2020, 2022); Amet et al. (2016) have demonstrated the potential
of proximitizing quantum Hall edge states and superconductivity. While these
studies are performed on various platforms the signature of Andreev reflection
is observed through negative downstream resistance and a reduction in Hall
(upstream) resistance in the quantum Hall regime. Specifically in Ref.
Hatefipour et al. (2022) we showed that InAs/NbTiN systems can exhibit up to
60% Andreev conversion. This is attributed to the high efficiency of Andreev
conversion at the InAs/NbTiN interface, resulting from the strong
hybridization of the quantum Hall edge modes with the states in the
superconductor. While cleanliness of the interface is crucial in achieving
high Andreev conversion, the underlying microscopic understanding of this
negative resistance is still not well-understood. These recent experimental
works have galvanized theoretical efforts to study Andreev processes in
quantum Hall edge state transport involving superconductivity Tang et al.
(2022); David et al. (2022); Kurilovich et al. (2023); Michelsen et al.
(2023); Kurilovich and Glazman (2022); Manesco et al. (2022); Beconcini et al.
(2018); Alavirad et al. (2018); Schiller et al. (2022); Galambos et al.
(2022); Cuozzo and Rossi (2023). In addition to the studies mentioned above,
recent study Vignaud et al. (2023) has shown evidence for chiral supercurrent
in quantum Hall Josephson junctions, and Ref. Uday et al. (2023) has reported
the observation of crossed Andreev reflection (CAR) across a narrow
superconducting Nb electrode contacting the chiral edge state of a quantum
anomalous Hall insulator (QAHI). These studies provide concrete demonstrations
of the hybridization of superconductivity and quantum Hall effects.

In this work, we introduce a new tool that is often used in studying
mesoscopic features of quantum Hall physics Bid et al. (2009); Gabelli et al.
(2007); Sánchez and Moskalets (2020); Taktak et al. (2022); Baer et al.
(2014); Zimmermann et al. (2017); Nayak et al. (2008) to control the number of
edge modes that interact with the superconductor. To achieve this we place a
voltage-controlled constriction, quantum point contact (QPC) exactly on the
interface between InAs and NbTiN to control the interplay between the QH edge
modes and the superconductor.

##  II II. sample growth and preparation

Figure 1: (a) Schematic of the grown and fabricated layers of the sample (the
region of the device specified by the red rectangle in panel (c)). (b)
Schematic diagram of the device showing the InAs/NbTiN interface with a QPC in
the integer quantum Hall regime. The solid and dashed lines represent
electrons and holes, respectively. When the QPC is activated, edge states
allowed to pass through the QPC can be Andreev reflected off the NbTiN
interface. (c) An optical photo of the finished device, highlighting the
different materials labeled by color and the contacts used for current source
and drain. (d) Downstream resistance (RD) as a function of perpendicular
magnetic field and QPC voltages. RD is measured between contacts 5 and 4. (e)
The 2DEG magnetotransport data of the fabricated device, where Rxx and Rxy are
measured between contacts 6-5 and 5-3, respectively. The plateaus are labeled
by their corresponding filling factors. Additionally, we show RD (between
contacts 5 and 4) and RU (between 3 and 4’) plots with blue and red color,
respectively.

The semiconductor in our study is a molecular beam epitaxy (MBE)
heterostructure grown on a semi-insulating InP(100) wafer. To form a two-
dimensional electron gas (2DEG), a quantum well (QW) was grown on an
InxAl1-xAs buffer with a graded indium content. The QW was formed by growing a
4nm layer of In${.81}$Ga${.19}$As, a 4nm layer of InAs, and a 10nm layer of
In.81Ga.19As. A delta-doped Si layer was placed below the QW at a depth of
6nm, with a doping concentration of $n\sim 1\times 10^{12}$ $cm^{-2}$. A
schematic of the stack layers is shown in Fig. 1(a). These InAs quantum wells
have been extensively studied in the context of mesoscopic superconductivity
and topological superconducting states. Previous studies have mainly
investigated the InAs/Al interface for developing tunable qubits and detecting
topological superconductivityLarsen et al. (2015); Casparis et al. (2018);
Mayer et al. (2019); Ren et al. (2019); Fornieri et al. (2019). Recently the
work has been extended to NbTiN Hatefipour et al. (2022) investigated the
InAs/NbTiN interface in the context of semiconductor-superconductor
heterostructures. In this new study, we aim to further explore this interface
using QPCs to investigate the interplay between the integer quantum Hall
effect (IQHE) and superconductivity with even greater precision.

##  III III. Device fabrication and measurement setup

We fabricated a Hall bar using electron beam lithography and chemical wet
etching. We cleaned the surface of the device using Argon plasma etching in
the sputtering tool at 25 W power for 15 s followed by deposition of a 90 nm
thick layer of NbTiN as the superconducting contacts. The interface between
the InAs and NbTiN was 150 $\mu$m long, and a QPC was added with a separation
of 150 nm between the QPC arms. A metallic gate for the QPC arms and pads was
created by depositing 20 nm aluminum oxide (Al2O3) as the gate dielectric
followed by 5 nm Cr and 20 nm Au e-beam deposition as the gate electrodes.
Fig.1(b) shows the schematic of the device zoomed in around the QPC and
interface area. The schematic demonstrates the edge modes reflection mechanism
in IQH regime and when QPC is activated. Fig.1(b) shows an optical photo of
the finished device labeling different region by their corresponding
materials. We have fabricated two samples (A and B) to confirm our
observations however we mainly focus on Sample A data in the main text and
show all the data for sample B in the Supplementary Information (SI). The
experiment was performed in a dilution fridge at a temperature of $T\sim 30$
mK with a maximum magnetic field of 12 T. The magnetotransport experiment was
carried out using lock-in amplifiers and an AC four-point measurement
technique with a frequency of $<20$ Hz and a $<1\mu$A AC excitation current.

##  IV IV. Measurement results

###  IV.1 A. Magnetotransport data

In order to assess the mobility of the quantum well (QW) in our fabricated
device, we conducted magnetotransport experiments on the 2DEG. As shown in
Fig. 1(c), we utilized contacts 1 and 4 as current source and drain, and
measured Rxx(6-5) and Rxy(5-3) as a function of magnetic field for mobility
and density analysis. Our measurements revealed that the mobility of the QW
for sample A is approximately $\mu\sim 13700$ $cm^{2}/V\cdot s$ at an electron
density of $n=1.54\times 10^{12}$ $cm^{-2}$. This mobility value corresponds
to an electron mean free path of approximately $l_{e}\sim 280nm$. The data for
$R_{xx}$ and $R_{xy}$ as a function of magnetic field is shown in Fig. 1(e),
which indicates that the sample has a relatively high density. Several
oscillations and plateaus in the longitudinal and Hall transport data are
observed in the 12T window, respectively, corresponding to filling factors
$\nu=6,8$ and $10$. Fig.1(d) displays the downstream resistance $R_{D}$
(measured between contacts 5 and 4) as a function of QPC voltage and
perpendicular magnetic field.

To understand the transport properties of mesoscopic systems, the Landauer-
Büttiker formalism is widely used Büttiker (1986). When a superconductor is
present, Andreev reflection (AR) processes can occur at the interface between
the semiconductor and the superconductor Hoppe et al. (2000). To include the
possibility of AR, we have extended the Landauer-Büttiker formalism describing
edge state transport in a Hall bar geometry. The effect of the AR is captured
by the average Andreev reflection $A\equiv(1/\nu)\sum_{i=1..\nu}A_{i}$, where
$\nu$ is the total number of edge modes, and $A_{i}$ is the probability of
Andreev reflection, i.e., the probability of an electron-like chiral state
approaching the QH-SC interface to be converted into a hole-like chiral state
when leaving the QH-SC interface, Hatefipour et al. (2022). As $V_{QPC}$ is
increased the distance of the QH edge modes from the QH-SC interface
increases. As a consequence the strength of the superconducting correlations
for the QH edge modes decreases, especially for the modes further removed from
the edge, resulting in a decrease of $A$.

Following our previous work Hatefipour et al. (2022), we begin by considering
the six-terminal configuration shown in Fig. 1(b). The terminals 1, 2, 3, 5,
and 6 are considered as ideal metallic leads, while contact 4 is a
superconducting lead. Let $I_{i}$ and $V_{i}$ denote the currents and
voltages, respectively, at the terminals $i=(1,2,...,6)$. Without loss of
generality, we set $V_{4}=0$. The conservation of charge equation
$\sum_{i}I_{i}=0$ is employed to express $I_{4}$ in terms of the currents at
the other leads. Using these considerations, the Landauer-Büttiker equations
take the form

$\displaystyle\begin{split}\begin{pmatrix}I_{1}\\\ I_{2}\\\ I_{3}\\\ I_{5}\\\
I_{6}\end{pmatrix}=&\frac{\nu}{R_{H}}\begin{pmatrix}1&0&0&0&-1\\\
-1&1&0&0&0\\\ 0&-1&1&0&0\\\ 0&0&2A-1&1&0\\\
0&0&0&-1&1\end{pmatrix}\begin{pmatrix}V_{1}\\\ V_{2}\\\ V_{3}\\\ V_{5}\\\
V_{6}\end{pmatrix},\end{split}$ (1)

where $\nu$ represents the total number of edge states, $R_{H}$, the Hall
resistance ($h/e^{2}$), and A, the average probability of Andreev reflection
per edge mode.

By assuming that no current flows into leads 2, 3, 5, and 6, we can simplify
the equation by setting $I_{2}=I_{3}=I_{5}=I_{6}=0$. We set the voltages at
terminals 1, 2, and 3 equal to each other $V_{1}=V_{2}=V_{3}$, and the
voltages at terminals 5 and 6 equal to each other $V_{5}=V_{6}$. Let
$I=I_{1}=-I_{4}$. With these assumptions, we can solve Eq. 1 to obtain

$\displaystyle R_{U}=\frac{V_{3}-V_{4}}{I}=\frac{R_{H}}{\nu}\frac{1}{2A}$ (2)
$\displaystyle
R_{D}=\frac{V_{5}-V_{4}}{I}=\frac{R_{H}}{\nu}\left(\frac{1}{2A}-1\right)$ (3)

The set of equations outlined above provide a theoretical framework for
studying edge state transport in the IQH/SC hybrid system.

Figure 2: (a) Extracted Andreev reflection data for three distinct filling
factors: 6, 8, and 10. This data is derived from both the RD and RU datasets,
and with a QPC voltage of 0V. (b) Dependence of the difference $R_{U}-R_{D}$
on QPC voltage. Figure 3: Measured downstream resistance ($R_{D}$) as a
function of QPC voltage for different values of $\nu$: (a) 6, (b) 8 and (c)
10. The magnetic field for each filling factor is labeled on top of each
panel.

Starting with the case where the QPC voltage is zero, we an extract $A$ using
Eq. (2-3) from $R_{U}$ and $R_{D}$. In Fig. 2(a) we show the extracted $A$ at
different filling factors for sample A studied in this work. We see that for
$V_{QPC}=0$ the value of $A$ falls within the range of $40\%$–$50\%$ and that
the extracted value of $A$ is found to be consistent between upstream and
downstream resistances. Sample B shows a similar trend at higher filling
factors filling factors ($\nu=12,~{}14,~{}16$) as shown in the SI.

Figure 3 shows the downstream resistance as a function of QPC voltage for
different fillings ($\nu=6,~{}8,~{}10$). In each case, the downstream
resistance remains constant for voltages $\gtrsim-1.5$ V before increasing. In
Fig. 4, we present the upstream resistance at the same magnetic fields and QPC
range as Fig. 3. There we also observe an increase in the resistance with
decreasing QPC voltage for $V_{QPC}<-1.5$ V. We can check the consistency of
Eqs. (2-3) with the experimental measurements by considering the difference
$(R_{U}-R_{D})$. According to Eq.(2-3), $(R_{U}-R_{D})=h/(\nu e^{2})$
regardless of the transmission through the QPC. From Fig. 2(b) we see the
difference in resistances remains fairly constant with QPC voltage $>-3$ V and
corresponds to the expected filling. This suggests that Eqs. (2-3)
consistently model the experimental data between the two edge state transport
regimes for $V_{QPC}>-3$ V. When the QPC voltage falls below $-3$ V, the RU
and RD quickly diverges and the Landuaer-Büttiker model breaks down.

Figure 4: Measured upstream resistance ($R_{U}$) as a function of QPC voltage
for different values of $\nu$: (a) 6, (b) 8 and (c) 10.  Figure 5: The
extracted $A$ parameter from both upstream (U) and downstream (D) resistance
for different filling factors $\nu$ as a function of QPC voltage.

From the experimental measurement of $R_{D}$ and $R_{U}$ shown in Figs. 3, 4,
using Eqs.(2) (3) we can obtain the evolution of $A$ with $V_{QPC}$ for
different filling factors $\nu$. Figure 5 shows the obtained average Andreev
conversion probability. We see that, regardless of the bulk filling factor,
$A$ is suppressed by the QPC. We also observe that for $V_{QPC}<-1$ V, $A$
first decreases continuously and then, for some ranges of $V_{QPC}$, exhibits
plateaus, as can be seen most clearly for the case $\nu=8$ when
$-2<V_{QPC}<-1.7$ V. The continuous decrease of $A$ with $V_{QPC}$ can be
interpreted as due to the progressive reduction of the superconducting pairing
correlations induced by the superconductor into the QH edge modes as the QPC
pushes away the edge modes from the QH-SC interface. As $V_{QPC}$ reaches
threshold values some of the edge modes further removed from the edge are
almost completely pinched off from the QH-SC interface resulting in a new,
reduced, value of $A$ that approximately is not affected by a further decrease
of $V_{QPC}$ until $V_{QPC}$ is large enough, in absolute value, to
significantly affect the pairing correlations of one more edge mode. Let the
number of edge states remaining be $\nu_{s}$. Then for the ranges of $V_{QPC}$
for which $A$ is approximately constant we can assume that $\nu_{s}<\nu$ modes
have a value of $A\neq 0$, while $\nu-\nu_{s}$ modes have $A\approx 0$. One of
the effects of the QPC is to also reduce the effective length, $L_{sc}$, of
the QH-SC interface. When $V_{QPC}=0$ we can assume $L_{sc}=L=150\mu$m. For
the largest absolute value of $V_{QPC}$ $L_{sc}$ can be taken to be equal to
the distance between the two gates forming the QPC: $L_{sc}=150~{}$nm. In both
cases $L_{sc}$ is much greater than the superconducting coherence length of
NbTiN ($\sim 10$ nm) so that crossed Andreev reflection processes can be
neglected Lee et al. (2017); Gül et al. (2022). However, one would expect that
the variation of $L_{sc}$ induced by $V_{QPC}$ would induce oscillations in
$A$. The lack of oscillations of $A$ with $V_{QPC}$ suggests that effects due
to disorder and vortices in NbTiN might play an important role resulting on an
effective averaging of $A$ over $L_{sc}$, as proposed in Ref. Hatefipour et
al. (2022), giving rise to a value of $A$ averaged along the length $L_{sc}$
that is independent of $L_{sc}$ as long as $L_{sc}$ is much large than the
superconducting coherence length $\xi$ of the superconductor.

##  V V. Conclusion

Our work demonstrates the successful fabrication of a hybrid device at the
InAs/NbTiN interface, incorporating a quantum point contact (QPC). We have
shown that by tuning the QPC voltage the effective, average, Andreev
conversion probability $A$ for QH edge modes can be tuned. We find that there
are threshold values of the QPC voltage for which some of the QH edge modes
appear to be completely pinched off by the QPC from the QH-SC interface
resulting in plateaus in the scaling of $A$ with $V_{QPC}$. The results also
show that the variation induced by the QPC of the effective length $L_{sc}$ of
the QH-SC interface does not result in oscillations of $A$. This is consistent
with the findings of Ref. Hatefipour et al. (2022) and suggests that effects
due to disorder and vortices in NbTiN must play an important role in
determining the properties of chiral Andreev states in InAs/NbTiN QH-SC
heterojunctions resulting on an effective averaging of $A$ over $L_{sc}$ that
is independent of $L_{sc}$ as long as $L_{sc}$ is much larger than the
superconductor’s coherence length. These findings advance the understanding of
QH-SC interfaces and should motivate future works to further elucidate the
details of the interplay of QH and superconducting states.

##  VI Acknowledgements

This work was supported by US DOE BES Grant No DE-SC0022245 and U.S. Army
Research Office Agreement No. W911NF1810067. The work at Sandia is supported
by a LDRD project.

Sandia National Laboratories is a multi-mission laboratory managed and
operated by National Technology & Engineering Solutions of Sandia, LLC
(NTESS), a wholly owned subsidiary of Honeywell International Inc., for the
U.S. Department of Energy’s National Nuclear Security Administration
(DOE/NNSA) under contract DE-NA0003525. This written work is authored by an
employee of NTESS. The employee, not NTESS, owns the right, title and interest
in and to the written work and is responsible for its contents. Any subjective
views or opinions that might be expressed in the written work do not
necessarily represent the views of the U.S. Government. The publisher
acknowledges that the U.S. Government retains a non-exclusive, paid-up,
irrevocable, world-wide license to publish or reproduce the published form of
this written work or allow others to do so, for U.S. Government purposes. The
DOE will provide public access to results of federally sponsored research in
accordance with the DOE Public Access Plan. ER thanks KITP, supported in part
by the National Science Foundation under Grants No. NSF PHY-1748958 and
PHY-2309135, where part of this work was performed.

References

## References

  * Kitaev (2001) A. Y. Kitaev, Physics-Uspekhi 44, 131 (2001).
  * Beenakker (2013) C. Beenakker, Annual Review of Condensed Matter Physics 4, 113 (2013), eprint https://doi.org/10.1146/annurev-conmatphys-030212-184337, URL https://doi.org/10.1146/annurev-conmatphys-030212-184337.
  * Alicea (2012) J. Alicea, Reports on Progress in Physics 75, 076501 (2012).
  * Nayak et al. (2008) C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Rev. Mod. Phys. 80, 1083 (2008).
  * Laucht et al. (2016) A. Laucht, S. Simmons, R. Kalra, G. Tosi, J. P. Dehollain, J. T. Muhonen, S. Freer, F. E. Hudson, K. M. Itoh, D. N. Jamieson, et al., Phys. Rev. B 94, 161302 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.161302.
  * Aghaee et al. (2022) M. Aghaee, A. Akkala, Z. Alam, R. Ali, A. A. Ramirez, M. Andrzejczuk, A. E. Antipov, P. Aseev, M. Astafev, B. Bauer, et al., _Inas-al hybrid devices passing the topological gap protocol_ (2022), URL https://arxiv.org/abs/2207.02472.
  * Dartiailh et al. (2021) M. C. Dartiailh, W. Mayer, J. Yuan, K. S. Wickramasinghe, A. Matos-Abiague, I. Žutić, and J. Shabani, Phys. Rev. Lett. 126, 036802 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.126.036802.
  * Lee et al. (2017) J. S. Lee, B. Shojaei, M. Pendharkar, A. P. McFadden, Y. Kim, H. J. Suominen, M. Kjaergaard, F. Nichele, C. M. Marcus, and C. J. Palmstrom, arXiv (2017).
  * Mong et al. (2014) R. S. K. Mong, D. J. Clarke, J. Alicea, N. H. Lindner, P. Fendley, C. Nayak, Y. Oreg, A. Stern, E. Berg, K. Shtengel, et al., Phys. Rev. X 4, 011036 (2014).
  * Clarke et al. (2013) D. J. Clarke, J. Alicea, and K. Shtengel, Nature Communications 4, 1348 EP (2013).
  * Hatefipour et al. (2022) M. Hatefipour, J. J. Cuozzo, J. Kanter, W. M. Strickland, C. R. Allemang, T.-M. Lu, E. Rossi, and J. Shabani, Nano Letters 22, 6173 (2022), pMID: 35867620, eprint https://doi.org/10.1021/acs.nanolett.2c01413, URL https://doi.org/10.1021/acs.nanolett.2c01413.
  * Gül et al. (2022) O. Gül, Y. Ronen, S. Y. Lee, H. Shapourian, J. Zauberman, Y. H. Lee, K. Watanabe, T. Taniguchi, A. Vishwanath, A. Yacoby, et al., Phys. Rev. X 12, 021057 (2022), URL https://link.aps.org/doi/10.1103/PhysRevX.12.021057.
  * Zhao et al. (2020) L. Zhao, E. G. Arnault, A. Bondarev, A. Seredinski, T. F. Q. Larson, A. W. Draelos, H. Li, K. Watanabe, T. Taniguchi, F. Amet, et al., Nature Physics 16, 862 (2020), URL https://doi.org/10.1038/s41567-020-0898-5.
  * Zhao et al. (2022) L. Zhao, Z. Iftikhar, T. F. Q. Larson, E. G. Arnault, K. Watanabe, T. Taniguchi, F. Amet, and G. Finkelstein, _Loss and decoherence at the quantum hall - superconductor interface_ (2022), eprint 2210.04842.
  * Amet et al. (2016) F. Amet, C. T. Ke, I. V. Borzenets, J. Wang, K. Watanabe, T. Taniguchi, R. S. Deacon, M. Yamamoto, Y. Bomze, S. Tarucha, et al., Science 352, 966 (2016), URL https://www.science.org/doi/abs/10.1126/science.aad6203.
  * Tang et al. (2022) Y. Tang, C. Knapp, and J. Alicea, Phys. Rev. B 106, 245411 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.106.245411.
  * David et al. (2022) A. David, J. S. Meyer, and M. Houzet, _Effects of the geometry on the downstream conductance in quantum hall-superconductor hybrid systems_ (2022), eprint 2210.16867.
  * Kurilovich et al. (2023) V. D. Kurilovich, Z. M. Raines, and L. I. Glazman, Nature Communications 14, 2237 (2023), ISSN 2041-1723, URL https://doi.org/10.1038/s41467-023-37794-1.
  * Michelsen et al. (2023) A. B. Michelsen, P. Recher, B. Braunecker, and T. L. Schmidt, Phys. Rev. Res. 5, 013066 (2023), URL https://link.aps.org/doi/10.1103/PhysRevResearch.5.013066.
  * Kurilovich and Glazman (2022) V. D. Kurilovich and L. I. Glazman, _Criticality in the crossed andreev reflection of a quantum hall edge_ (2022), eprint 2209.12932.
  * Manesco et al. (2022) A. L. R. Manesco, I. M. Flór, C.-X. Liu, and A. R. Akhmerov, SciPost Phys. Core 5, 045 (2022), URL https://scipost.org/10.21468/SciPostPhysCore.5.3.045.
  * Beconcini et al. (2018) M. Beconcini, M. Polini, and F. Taddei, Phys. Rev. B 97, 201403 (2018), URL https://link.aps.org/doi/10.1103/PhysRevB.97.201403.
  * Alavirad et al. (2018) Y. Alavirad, J. Lee, Z.-X. Lin, and J. D. Sau, Physical Review B 98, 214504 (2018), ISSN 2469-9950, URL https://link.aps.org/doi/10.1103/PhysRevB.98.214504.
  * Schiller et al. (2022) N. Schiller, B. A. Katzir, A. Stern, E. Berg, N. H. Lindner, and Y. Oreg, _Interplay of superconductivity and dissipation in quantum hall edges_ (2022), eprint 2202.10475.
  * Galambos et al. (2022) T. H. Galambos, F. Ronetti, B. Hetényi, D. Loss, and J. Klinovaja, Phys. Rev. B 106, 075410 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.106.075410.
  * Cuozzo and Rossi (2023) J. J. Cuozzo and E. Rossi, _Su(4) symmetry breaking and induced superconductivity in graphene quantum hall edges_ (2023), eprint arXiv:2306.12483.
  * Vignaud et al. (2023) H. Vignaud, D. Perconte, W. Yang, B. Kousar, E. Wagner, F. Gay, K. Watanabe, T. Taniguchi, H. Courtois, Z. Han, et al., _Evidence for chiral supercurrent in quantum hall josephson junctions_ (2023), eprint 2305.01766.
  * Uday et al. (2023) A. Uday, G. Lippertz, K. Moors, H. F. Legg, A. Bliesener, L. M. C. Pereira, A. A. Taskin, and Y. Ando, _Induced superconducting correlations in the quantum anomalous hall insulator_ (2023), eprint 2307.08578.
  * Bid et al. (2009) A. Bid, N. Ofek, M. Heiblum, V. Umansky, and D. Mahalu, Phys. Rev. Lett. 103, 236802 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.103.236802.
  * Gabelli et al. (2007) J. Gabelli, G. Fève, T. Kontos, J.-M. Berroir, B. Placais, D. C. Glattli, B. Etienne, Y. Jin, and M. Büttiker, Phys. Rev. Lett. 98, 166806 (2007), URL https://link.aps.org/doi/10.1103/PhysRevLett.98.166806.
  * Sánchez and Moskalets (2020) D. Sánchez and M. Moskalets, Entropy 22 (2020), ISSN 1099-4300, URL https://www.mdpi.com/1099-4300/22/9/977.
  * Taktak et al. (2022) I. Taktak, M. Kapfer, J. Nath, P. Roulleau, M. Acciai, J. Splettstoesser, I. Farrer, D. A. Ritchie, and D. C. Glattli, Nature Communications 13, 5863 (2022), ISSN 2041-1723, URL https://doi.org/10.1038/s41467-022-33603-3.
  * Baer et al. (2014) S. Baer, C. Rössler, T. Ihn, K. Ensslin, C. Reichl, and W. Wegscheider, Phys. Rev. B 90, 075403 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.075403.
  * Zimmermann et al. (2017) K. Zimmermann, A. Jordan, F. Gay, K. Watanabe, T. Taniguchi, Z. Han, V. Bouchiat, H. Sellier, and B. Sacépé, Nature Communications 8, 14983 (2017), ISSN 2041-1723, URL https://doi.org/10.1038/ncomms14983.
  * Larsen et al. (2015) T. W. Larsen, K. D. Petersson, F. Kuemmeth, T. S. Jespersen, P. Krogstrup, J. Nygård, and C. M. Marcus, Phys. Rev. Lett. 115, 127001 (2015).
  * Casparis et al. (2018) L. Casparis, M. R. Connolly, M. Kjaergaard, N. J. Pearson, A. Kringhøj, T. W. Larsen, F. Kuemmeth, T. Wang, C. Thomas, S. Gronin, et al., Nature Nanotechnology 13, 915 (2018), URL https://doi.org/10.1038/s41565-018-0207-y.
  * Mayer et al. (2019) W. Mayer, M. C. Dartiailh, J. Yuan, K. S. Wickramasinghe, A. Matos-Abiague, I. Žutić, and J. Shabani, arXiv e-prints arXiv:1906.01179 (2019), eprint 1906.01179.
  * Ren et al. (2019) H. Ren, F. Pientka, S. Hart, A. T. Pierce, M. Kosowsky, L. Lunczer, R. Schlereth, B. Scharf, E. M. Hankiewicz, L. W. Molenkamp, et al., Nature 569, 93 (2019), ISSN 1476-4687, URL https://doi.org/10.1038/s41586-019-1148-9.
  * Fornieri et al. (2019) A. Fornieri, A. M. Whiticar, F. Setiawan, E. Portolés, A. C. C. Drachmann, A. Keselman, S. Gronin, C. Thomas, T. Wang, R. Kallaher, et al., Nature 569, 89 (2019), ISSN 1476-4687, URL https://doi.org/10.1038/s41586-019-1068-8.
  * Büttiker (1986) M. Büttiker, Phys. Rev. Lett. 57, 1761 (1986), URL https://link.aps.org/doi/10.1103/PhysRevLett.57.1761.
  * Hoppe et al. (2000) H. Hoppe, U. Zülicke, and G. Schön, Phys. Rev. Lett. 84, 1804 (2000), URL https://link.aps.org/doi/10.1103/PhysRevLett.84.1804.


# Reservoir computing with the Kuramoto model

Hayato Chiba Advanced Institute for Materials Research, Tohoku University,
Japan. Koichi Taniguchi  Department of Mathematical and Systems Engineering,
Faculty of Engineering, Shizuoka University, Japan. Takuma Sumi∗

(Last update: 13 July 2024)

###### Abstract

Reservoir computing aims to achieve high-performance and low-cost machine
learning with a dynamical system as a reservoir. However, in general, there
are almost no theoretical guidelines for its high-performance or optimality.
This paper focuses on the reservoir computing with the Kuramoto model and
theoretically reveals its approximation ability. The main result provides an
explicit expression of the dynamics of the Kuramoto reservoir by using the
order parameters. Thus, the output of the reservoir computing is expressed as
a linear combination of the order parameters. As a corollary, sufficient
conditions on hyperparameters are obtained so that the set of the order
parameters gives the complete basis of the Lebesgue space. This implies that
the Kuramoto reservoir has a universal approximation property. Furthermore,
the conjecture on the edge of bifurcation, which is a generalization of the
famous criterion the edge of chaos for designing a high-performance reservoir,
is also discussed from the viewpoint of its approximation ability. It is
numerically demonstrated by a prediction task and a transformation task.

##  1 Introduction

Reservoir computing is one of the methods of recurrent machine learning with
dynamical systems for computations on time series data (see [5, 10]) and has
various applications such as time series prediction, system identification,
signal generation and signal transformation (see [9, 16]). The architecture of
reservoir computing consists of two components, namely, reservoir and readout.
The reservoir is a fixed recurrent nonlinear map (dynamical system) which
appropriately maps input data to a high dimensional space. The readout is a
time-independent linear map (matrix), and only the readout is learned to
minimize the error between its output and the desired result. Since the
learning of reservoir computing is done only by simple algorithms such as
linear regression, reservoir computing has a distinctive advantage to make it
possible to learn at low cost and to optimize instantly. On the other hand,
its ability to achieve high performance is strongly dependent on the choice of
a reservoir.

One of the well-known criteria for designing a high-performance reservoir is
so-called the edge of chaos (or edge of stability) [11, 8, 3]. The edge of
chaos is the phase transition from an ordered state to a chaotic state, and it
seems that reservoir computing has its greatest computational capacity near
the phase transition point. The effectiveness of the edge of chaos has been
experimentally confirmed in various settings (see e.g. [1, 13, 15] for the
details of the edge of chaos). However, to our knowledge, there are no
theoretical guarantees or mathematical proofs for the edge of chaos.

Recently, it has been reported that the performance of reservoir computing can
be improved without the edge of chaos. Specifically, Goto, Nakajima and Notsu
[4] studied reservoir computing with the Navier-Stokes equation as a reservoir
and showed that the optimal computational performance is achieved near the
critical Reynolds number for the Hopf bifurcation. In light of this recent
study, we propose the edge of bifurcation, which is a generalization of the
edge of chaos, and means that the performance gets better when the dynamical
system undergoes some bifurcations (e.g. Hopf bifurcation).

In this paper, we tackle the theoretical analysis of the edge of bifurcation
and aim to provide mathematical insights in this topic. For this purpose, we
focus on the Kuramoto model as a reservoir, which is the most typical
mathematical model for synchronization phenomena. Recently, reservoir
computing with the Kuramoto model has been studied numerically (see [14, 17]),
but its potential has not yet been fully clarified, and there is a lack of
theoretical results on its properties. The purpose of this paper is to
theoretically clarify the approximation ability of reservoir computing with
the Kuramoto model and to provide mathematical insights on the edge of
bifurcation for reservoir computing.

The basic results of the Kuramoto model, which are the fundament of this
paper, are summarized in Section 2. A summary of our results and contributions
is as follows:

  * •

(Formulation of reservoir computing and numerical method) The first
contribution is to formulate the Kuramoto reservoir computing with the aid of
a certain integral operator for mathematical analysis. More precisely, we
propose reservoir computing by using the infinite-dimensional Kuramoto model
as a reservoir, which is a continuous limit ($N\to\infty$) of the Kuramoto
model as the number $N$ of oscillators tends to infinity, and a readout is not
a matrix but an integral operator (see Section 3). Furthermore, our reservoir
computing is easily implementable and its specific numerical calculation
procedure is also given via the Fourier series (see Section 5).

  * •

(Approximation ability of the Kuramoto reservoir) The second contribution is
to theoretically clarify the approximation ability of reservoir computing with
the Kuramoto model. We can obtain a useful representation of the Kuramoto
reservoir computing via the Fourier series. More precisely, we assume that the
integral kernel of the integral operator for the readout belongs to
$L^{2}(0,2\pi)$. Then the output function can be represented by a linear
combination of the $n$-th order parameters of a solution of the infinite-
dimensional Kuramoto model (see (3.5)). Therefore, if the family of $n$-th
order parameters constructs a complete basis of $L^{2}_{\mathrm{per}}$ (the
set of periodic $L^{2}$ functions), the Kuramoto reservoir has an enough
approximation ability. Our main theorem provides an explicit expression of the
$n$-th order parameters (Theorem 4.1) . As a corollary, a sufficient condition
for this family to be a complete basis of $L^{2}_{\mathrm{per}}$ (Corollary
4.3) will be obtained. The precise statements and proofs are given in Section
4.

  * •

(Edge of bifurcation) The third contribution is to provide some mathematical
insights on the edge of bifurcation conjecture for reservoir computing with
the Kuramoto model from a viewpoint of the approximation ability. A key
parameter is the coupling strength $K$ of the Kuramoto model, which determines
the occurrence of synchronization. There is a critical value $K_{c}$ such that
a bifurcation occurs at $K=K_{c}$ under certain conditions (see [6, 7, 2]);
the de-synchronous state is stable when $K<K_{c}$ and the synchronous state
appears when $K>K_{c}$. Hence, it is expected that this $K_{c}$ has a strong
influence on the performance of the Kuramoto reservoir computing in terms of
the edge of bifurcation. In this paper, we consider the edge of bifurcation in
some settings based on our approximation theorem and the bifurcation theory of
the Kuramoto model. As a consequence, we can see that the type of bifurcation
and the shape of bifurcation diagram are important in the performance of
Kuramoto reservoir computing. For example, it turns out that a bifurcation
should be not a pitchfork bifurcation (steady to steady), but be a Hopf
bifurcation (steady to periodic). Further, we will show that a sufficient
condition for the family of the $n$-th order parameters to be the basis of
$L^{2}_{\mathrm{per}}$ is satisfied when $K$ is slightly larger than $K_{c}$,
that guarantees the edge of bifurcation. See Section 4 for the details of
these observations. This result is also numerically supported in Section 5.

##  2 Kuramoto model

The (finite-dimensional) Kuramoto model is the system of ordinary differential
equations given by

$\frac{d\theta_{i}}{dt}=\omega_{i}+\frac{K}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i}),\quad
i=1,\cdots,N,$ (2.1)

where $i\in[N]:=\\{1,2,\ldots,N\\}$, $\theta_{i}=\theta_{i}(t)\in[0,2\pi)$
denotes the phase of an $i$-th oscillator on a circle,
$\omega_{i}\in\mathbb{R}$ is its natural frequency distributed according to
some probability density function $g(\omega)$, $K>0$ is a coupling strength,
and $N$ is the number of oscillators. We normally assume that $g(\omega)$ is
an even and unimodal function about its mean frequency $\Omega\in\mathbb{R}$:
(even) $g(\Omega+\omega)=g(\Omega-\omega)$ for any $\omega\in\mathbb{R}$;
(unimodal) $g(\omega_{1})>g(\omega_{2})$ for $\Omega\leq\omega_{1}<\omega_{2}$
and $g(\omega_{1})<g(\omega_{2})$ for $\omega_{1}<\omega_{2}\leq\Omega$. This
system was proposed by Kuramoto [6] in order to investigate collective
synchronization phenomena and was derived by means of the averaging method
from coupled dynamical systems having limit cycles. The order parameter is
usually used to measure the degree of synchronization of the oscillators
$\\{\theta_{i}\\}_{i=1}^{N}$ and is defined by

$r(t):=\frac{1}{N}\sum_{j=1}^{N}e^{i\theta_{j}(t)},$ (2.2)

where $i=\sqrt{-1}$. The order parameter gives the centroid of oscillators. If
$|r|>0$, then the synchronization occurs, while if $|r|$ is nearly equal to
zero, the oscillators are uniformly distributed (de-synchronization), see Fig.
1.

######  Remark 2.1.

For the study of the Kuramoto model, we can assume without loss of generality
that $\Omega=0$. Indeed, changing variables for (2.1) as
$\theta_{i}\mapsto\theta_{i}+\Omega t$ yields

$\frac{d\theta_{i}}{dt}=\omega_{i}-\Omega+\frac{K}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i}),\quad
i=1,\cdots,N.$ (2.3)

Then, the mean value of $\\{\omega_{i}-\Omega\\}_{i}$ is zero. Nevertheless,
we will see that for an application to a reservoir computing, nonzero $\Omega$
plays an important role (Corollary 4.3).

Figure 1: Schematic representation of oscillators on a circle. The left
represents a de-synchronous state with $|r|\approx 0$. The right represents a
synchronous state with $|r|\approx 1$.

The coupling strength $K$ is the key bifurcation parameter that determines
whether the dynamics becomes synchronous or de-synchronous, and it is
important to obtain the bifurcation diagram of $r(t)$ for the study of the
edge of bifurcation. In this paper, we address the case where $N$ is large and
consider the continuous limit of (2.1) as $N\to\infty$, namely, the infinite-
dimensional Kuramoto model given by

$\begin{cases}\dfrac{\partial\rho}{\partial
t}+\dfrac{\partial}{\partial\theta}(\rho v)=0,\\\
v=\omega+K|r|\sin(\psi-\theta),\\\ \displaystyle
r(t)=|r(t)|e^{i\psi(t)}:=\int_{\mathbb{R}}\int_{0}^{2\pi}e^{i\theta}\rho(\theta,\omega,t)g(\omega)\,d\theta
d\omega,\end{cases}$ (2.4)

where $\rho=\rho(\theta,\omega,t)$ is the density of oscillators at a phase
$\theta$ parameterized by natural frequency $\omega$ and time $t$,
$v=v(\theta,\omega,t)$ is the drift velocity of the oscillators, and $r(t)$ is
the order parameter of a solution $\rho$. For any $K>0$, the trivial steady
state solution of (2.4) exists (i.e. a completely incoherent state given by
$\rho(\theta,\omega)=1/(2\pi)$ for any $\theta,\omega$) and it corresponds to
the de-synchronous state $r\equiv 0$. Kuramoto [7] conjectured that a
bifurcation diagram of $r(t)$ is given as Fig. 2 (a):

Kuramoto conjecture

Suppose that $N\to\infty$ and natural frequencies $\omega_{i}$’s are
distributed according to an even and unimodal function $g(\omega)$ about its
mean value $\Omega$ with $g^{\prime\prime}(\Omega)<0$ (this is true for Fig. 2
(a)). If $K<K_{c}:=2/(\pi g(\Omega))$, then $|r(t)|\equiv 0$ is asymptotically
stable, while if $K>K_{c}$, the synchronous state emerges and $|r(t)|\equiv
r_{c}$ is asymptotically stable for some $r_{c}>0$ given by

$r_{c}=\sqrt{\frac{-16}{\pi
K_{c}^{4}g^{\prime\prime}(0)}}\sqrt{K-K_{c}}+O(K-K_{c})$ (2.5)

(see [12] for Kuramoto’s discussion).

The critical value $K_{c}$ is called the Kuramoto transition point. Then,
Chiba [2] mathematically formulated this problem and proved the Kuramoto
conjecture. See [2] for the precise statement.

It is known that the type of bifurcation strongly depends on the assumption
for $g(\omega)$, especially for the sign of $g^{\prime\prime}(\Omega)$. For
example, when $g(\omega)$ is a uniform distribution (i.e.
$g^{\prime\prime}(\Omega)=0$), the bifurcation diagram is given as Fig. 2 (b).
In this paper, as a reservoir, we consider both cases (a) $g$ is Gaussian (b)
$g$ is a uniform distribution.

Figure 2:  Bifurcation diagrams of the order parameter (lower panel) for the
cases where (a) $g(\omega)$ is a unimodal Gaussian (b) $g(\omega)$ is a
uniform distribution with the average $\Omega$. The bifurcation diagrams are
independent of $\Omega$ as was mentioned in the previous Remark.

##  3 Reservoir computing with the Kuramoto model

We introduce the finite-dimensional Kuramoto model with an input term

$\frac{d\theta_{i}}{dt}=\alpha
u+\omega_{i}+\frac{K}{N}\sum_{j=1}^{N}\sin(\theta_{j}-\theta_{i}),$ (3.1)

where $i\in[N]$, $\theta_{i}=\theta_{i}(\alpha u;t)\in[0,2\pi)$ denotes the
phase of an $i$-th oscillator on a circle, $\omega_{i}\in\mathbb{R}$ denotes
its natural frequency, $K>0$ is a coupling strength, $N$ is the number of
oscillators, $\alpha>0$ is a constant and $u=u(t)$ is an input. The case
$\alpha=0$ is the (usual) Kuramoto model (2.1). In reservoir computing, we
take a sufficiently small $\alpha$ and treat the input-aware model. The order
parameter $r(t)$ is similarly defined by (2.2). In addition, we define the
$n$-th order parameter of $\\{\theta_{j}\\}_{j=1}^{N}$ by

$r_{n}(\alpha u;t):=\frac{1}{N}\sum_{j=1}^{N}e^{in\theta_{j}(\alpha
u;t)},\quad n\in\mathbb{Z}.$

Similarly to Section 2, we consider the continuous limit of (3.1) as
$N\to\infty$: with given natural frequency $\omega$ at time $t$, $\rho$
satisfies the continuity equation

$\begin{cases}\dfrac{\partial\rho}{\partial
t}+\dfrac{\partial}{\partial\theta}(\rho v)=0,\\\ v=\alpha
u+\omega+K|r|\sin(\psi-\theta),\\\ \displaystyle r(\alpha u;t)=|r(\alpha
u;t)|e^{i\psi(\alpha
u;t)}:=\int_{\mathbb{R}}\int_{0}^{2\pi}e^{i\theta}\rho(\alpha
u;\theta,\omega,t)g(\omega)\,d\theta d\omega,\end{cases}$ (3.2)

where $\rho=\rho(\alpha u;\theta,\omega,t)$ is the density of oscillators at a
phase $\theta$ parameterized by a natural frequency $\omega$ at time $t$,
$v=v(\alpha u;\theta,\omega,t)$ is the drift velocity of the oscillators, and
$r(\alpha u;t)$ is the order parameter of $\rho$. The $n$-th order parameter
is defined as

$r_{n}(\alpha u;t):=\int_{\mathbb{R}}\int_{0}^{2\pi}e^{in\theta}\rho(\alpha
u;\theta,\omega,t)g(\omega)\,d\theta d\omega,\quad n\in\mathbb{Z}.$ (3.3)

We propose reservoir computing with the infinite-dimensional Kuramoto model
(3.2) as follows:

Reservoir computing and its purpose.

The reservoir computing with (3.2) consists of mappings
$u(t)\mapsto\rho(\alpha u;\theta,\omega,t)\mapsto y(t)$, where
$y(t)\in\mathbb{R}$ is an output function. The map $u\mapsto\rho$ is defined
by a solution of (3.2). Our purpose is to find the mapping (readout)
$\rho(\alpha u;\theta,\omega,t)\mapsto y(t)$ so that

$\|y_{\mathrm{tar}}-y\|\approx 0,$

for a given target function $y_{\text{tar}}(t)\in\mathbb{R}$, where
$\|\cdot\|$ is a certain norm or metric. For a finite dimensional problem,
usually a readout is given by a matrix. However, since (3.2) is an infinite
dimensional dynamical system, a readout may be a linear operator on an
infinite dimensional space. Thus, we consider the following linear integral
operator

$\rho\mapsto y(t)=\mathcal{W}_{\mathrm{out}}[\rho](\alpha
u,t):=\int_{\mathbb{R}}\int_{0}^{2\pi}W_{\mathrm{out}}(\theta)\rho(\alpha
u;\theta,\omega,t)g(\omega)\,d\theta d\omega,$ (3.4)

where $\mathcal{W}_{\mathrm{out}}$ is a linear integral operator with an
integral kernel $W_{\mathrm{out}}(\theta)$. Hence, the purpose of our
reservoir computing is to find a suitable function $W_{\mathrm{out}}(\theta)$
such that $\|y_{\mathrm{tar}}-y\|\approx 0$ holds.

For mathematical analysis, we always assume that $W_{\mathrm{out}}\in
L^{2}(0,2\pi)$ and it is real-valued in this paper. It is important to note
that our reservoir computing (3.4) is convenient because it is expressed by
the $n$-th order parameters as follows. The integral kernel
$W_{\mathrm{out}}\in L^{2}(0,2\pi)$ has a Fourier expansion

$W_{\mathrm{out}}(\theta)=\sum_{n=-\infty}^{\infty}w_{n}e^{in\theta}.$

Then, by the definition (3.3) of $r_{n}(\alpha u;t)$, we rewrite
$\mathcal{W}_{\mathrm{out}}[\rho]$ as

$\mathcal{W}_{\mathrm{out}}[\rho](\alpha
u,t)=\sum_{n=-\infty}^{\infty}w_{n}\int_{\mathbb{R}}\int_{0}^{2\pi}e^{in\theta}\rho(\alpha
u;\theta,\omega,t)g(\omega)\,d\theta
d\omega=\sum_{n=-\infty}^{\infty}w_{n}r_{n}(\alpha u;t).$

Hence, the output of our reservoir computing can be represented as

$y(t)=\sum_{n=-\infty}^{\infty}w_{n}r_{n}(\alpha u;t).$ (3.5)

This is an easily implementable setting and we can compute
$W_{\text{out}}(\theta)$ and $y(t)$ numerically by using these formulae.
Before giving its numerical algorithm, we consider mathematical results in the
next section.

##  4 Approximation ability

As mentioned in the previous section, the approximation ability of reservoir
computing (3.4) is essentially determined by the property of $\\{r_{n}(\alpha
u;t)\\}_{n}$. In this section, we focus on the analysis of $\\{r_{n}(\alpha
u;t)\\}_{n}$. Our main results are the following:

######  Theorem 4.1.

Let $K>0$, $\alpha\in\mathbb{R}$ and $\Omega$ mean value of $g(\omega)$.
Suppose that $u(t)$ is bounded and continuous on $\mathbb{R}$.

  * (i)

For a steady state solution $\rho_{*}$ of the system (3.2) with
$\Omega=\alpha=0$, its order parameters $r_{n}=:c_{n}$ are real numbers for
any $n\in\mathbb{Z}$.

  * (ii)

For the system (3.2), there exists a solution $\rho$ constructed from
$\rho_{*}$ such that the $n$-th order parameters of $\rho$ is expressed as

$r_{n}(\alpha u;t)=c_{n}e^{in(\Omega t+\alpha\int_{0}^{t}u(s)\,ds)},\quad
t>0,\,\,c_{n}\in\mathbb{R}.$ (4.1)

######  Remark 4.2.

A few remarks are in order.

  1. (a)

The system (3.2) for $\Omega=\alpha=0$ always has the steady state
$\rho_{*}=1/(2\pi)$ (de-synchronous state). In this case, $c_{n}=0$ for $n\neq
0$. Thus, we are interested in a non-trivial steady state. For example, when
$g(\omega)$ is even and unimodal or uniform distribution, it has a non-trivial
steady state $\rho_{*}$ after the bifurcation $(K>K_{c})$ to the synchronous
state as in Fig. 2. In this case, $c_{n}\not=0$ in general (see Theorem 4.6
below).

  2. (b)

For the system (3.2) with $\Omega=\alpha=0$, the steady state
$\rho_{*}=1/(2\pi)$ is stable when $K<K_{c}$ and the non-trivial steady state
$\rho_{*}$ is stable when $K>K_{c}$ in some sense (see [2] for the details).
Moreover, the solution $\rho$ in (ii) is stable if $\rho_{*}$ in (i) is
stable.

######  Corollary 4.3.

Under the same assumption as in Theorem 4.1, if $c_{n}\not=0$ for any $n$,
$\Omega\not=0$ and $|\alpha|$ is sufficiently small, then $\\{r_{n}(\alpha
u;t)\\}_{n}$ given in (4.1) is a complete basis of $L^{2}_{\mathrm{per}}(0,T)$
with $T\sim 2\pi/\Omega+O(\alpha)$.

Our main results mean that the Kuramoto reservoir has almost the same
approximation ability as the Fourier basis when $c_{n}\not=0$ for any $n$,
$\Omega\not=0$ and $|\alpha|$ is sufficiently small. See Theorem 4.6 for a
sufficient condition for $c_{n}\not=0$.

To prove Theorem 4.1, we use the following lemma, which is easily proved by
integrating $(\partial/\partial\theta)(\rho v)=0$.

######  Lemma 4.4.

Let $\Omega=0$ and $\alpha=0$. A steady state solution $\rho_{*}$ of (3.2) is
given as

$\rho_{*}(\theta,\omega)=\begin{dcases}\delta\left(\theta-\mathrm{Arcsin}\,\left(\frac{\omega}{Kr_{*}}\right)\right),\quad&|\omega|<Kr_{*},\\\
\frac{1}{2\pi}\frac{\sqrt{\omega^{2}-K^{2}r_{*}^{2}}}{|\omega-
Kr_{*}\sin\theta|},&|\omega|>Kr_{*},\end{dcases}$

where $\delta$ is the Dirac delta function and $r_{*}$ is the $1$-st order
parameter of $\rho_{*}$.

###### Proof of Theorem 4.1.

In order to prove the statement (ii), it is enough to show the statement (i).
In fact, for the system (3.2), change the variables by

$\begin{dcases}\theta=\tilde{\theta}+\Omega t+\alpha\int_{0}^{t}u(s)\,ds,\\\
\omega=\tilde{\omega}+\Omega,\\\
g(\tilde{\omega}+\Omega)=\tilde{g}(\tilde{\omega}),\\\ \rho(\alpha
u;\theta,\omega,t)=\tilde{\rho}(\tilde{\theta},\tilde{\omega},t).\end{dcases}$

Then, it is easy to verify that the system is reduced to the case
$\Omega=\alpha=0$. More precisely, let $\tilde{r}_{n}$ be the order parameters
for the case $\Omega=\alpha=0$. Then we can calculate $r_{n}$ as

$\begin{split}r_{n}(\alpha
u;t)&=\int_{\mathbb{R}}\int_{0}^{2\pi}e^{in\theta}\rho(\alpha
u;\theta,\omega,t)g(\omega)\,d\theta d\omega\\\
&=\int_{\mathbb{R}}\int_{0}^{2\pi}e^{in\tilde{\theta}}e^{in(\Omega
t+\alpha\int_{0}^{t}u(s)\,ds)}\tilde{\rho}(\tilde{\theta},\tilde{\omega},t)\tilde{g}(\tilde{\omega})\,d\tilde{\theta}d\tilde{\omega}\\\
&=\tilde{r}_{n}(t)e^{in(\Omega t+\alpha\int_{0}^{t}u(s)\,ds)},\end{split}$

Therefore, once we prove $\tilde{r}_{n}(t)=:c_{n}\in\mathbb{R}$ when
$\tilde{\rho}=\rho_{*}$ is a steady state, we conclude the statement (ii).

Suppose $\Omega=0$ and $\alpha=0$, and let $\rho_{*}$ be a steady state. By
Lemma 4.4, its $n$-th order parameters are expressed as

$\begin{split}r_{n}&=\int_{\mathbb{R}}\int_{0}^{2\pi}e^{in\theta}\rho_{*}(\theta,\omega)g(\omega)\,d\theta
d\omega\\\
&=\int_{|\omega|<Kr_{*}}\int_{0}^{2\pi}e^{in\theta}\delta\left(\theta-\mathrm{Arcsin}\,\left(\frac{\omega}{Kr_{*}}\right)\right)g(\omega)\,d\theta
d\omega\\\
&\qquad\qquad\qquad+\frac{1}{2\pi}\int_{|\omega|>Kr_{*}}\int_{0}^{2\pi}e^{in\theta}\frac{\sqrt{\omega^{2}-K^{2}r_{*}^{2}}}{|\omega-
Kr_{*}\sin\theta|}g(\omega)\,d\theta d\omega\\\
&=:a_{n}+b_{n}=:c_{n}.\end{split}$

Since $g(\omega)$ is an even function about $0$, we write

$\begin{split}a_{n}=\int_{|\omega|<Kr_{*}}e^{in\mathrm{Arcsin}\,\left(\frac{\omega}{Kr_{*}}\right)}g(\omega)\,d\omega=\int_{|\omega|<Kr_{*}}e^{-in\mathrm{Arcsin}\,\left(\frac{\omega}{Kr_{*}}\right)}g(\omega)\,d\omega.\end{split}$

By Euler’s formula, this implies that $a_{n}$ is a real number for any $n$ as

$a_{n}=\int_{|\omega|<Kr_{*}}\cos\left(n\cdot\mathrm{Arcsin}\,\left(\frac{\omega}{Kr_{*}}\right)\right)g(\omega)\,d\omega\in\mathbb{R}.$
(4.2)

As for $b_{n}$, by making the change
$(\theta,\omega)\mapsto(\theta+\pi,-\omega)$, we have

$b_{n}=\frac{1}{2\pi}\int_{|\omega|>Kr_{*}}\int_{-\pi}^{\pi}e^{in\theta}e^{in\pi}A(\omega,\theta)g(\omega)\,d\theta
d\omega=\begin{cases}-b_{n}\quad&\text{if $n$ is odd},\\\ b_{n}\quad&\text{if
$n$ is even},\end{cases}$

where

$A(\omega,\theta):=\frac{\sqrt{\omega^{2}-K^{2}r_{*}^{2}}}{|\omega-
Kr_{*}\sin\theta|}\in\mathbb{R}.$

This implies $b_{n}=0$ if $n$ is odd (this fact will be used later). For any
$n$, we can also write $b_{n}$ as

$b_{n}=\frac{1}{2\pi}\int_{|\omega|>Kr_{*}}\left(\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}+\int_{\frac{\pi}{2}}^{\frac{3\pi}{2}}\right)e^{in\theta}A(\omega,\theta)g(\omega)\,d\theta
d\omega,$

and by making the change $(\theta,\omega)\mapsto(\theta+\pi,-\omega)$ for the
second integral, we have

$\displaystyle b_{n}$
$\displaystyle=\frac{1}{\pi}\int_{|\omega|>Kr_{*}}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos(n\theta)A(\omega,\theta)g(\omega)\,d\theta
d\omega$
$\displaystyle=\frac{2}{\pi}\int^{\infty}_{Kr_{*}}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos(n\theta)A(\omega,\theta)g(\omega)\,d\theta
d\omega\in\mathbb{R}.$ (4.3)

Thus, $c_{n}\in\mathbb{R}$ and the proof of Theorem 4.1 is completed. ∎

To prove Corollary 4.3, we use the following lemma.

######  Lemma 4.5.

Let $\\{g_{k}\\}_{k}$ be a family of elements in a Hilbert space $H$ with
inner product $(\cdot,\cdot)_{H}$. Suppose that $f=0$ if $f\in H$ and
$(f,g_{k})_{H}=0$ for any $k$. Then $\\{g_{k}\\}_{k}$ is complete in $H$.

###### Proof of Corollary 4.3.

Define $T>0$ by

$2\pi=\Omega T+\alpha\int_{0}^{T}u(s)\,ds.$

For any $f\in L^{2}(0,T)$, we suppose

$\int_{0}^{T}f(t)e^{in(\Omega t+\alpha\int_{0}^{t}u(s)\,ds)}\,dt=0$

for any $n$. Put $\tilde{t}=\Omega t+\alpha\int_{0}^{t}u(s)\,ds$. If $u$ is
bounded and continuous, $\Omega\not=0$, and $|\alpha|$ is sufficiently small,
then the map $\eta:t\mapsto\tilde{t}$ is one-to-one and of $C^{1}$. Hence, the
inverse function theorem gives

$\int_{0}^{2\pi}\frac{f(\eta^{-1}(\tilde{t}))}{\Omega+\alpha
u(\eta^{-1}(\tilde{t}))}e^{in\tilde{t}}\,d\tilde{t}=0$

for any $n$. Since $\\{e^{in\tilde{t}}\\}_{n}$ is an orthogonal basis of
$L^{2}(0,2\pi)$, this implies that

$\frac{f(\eta^{-1}(\tilde{t}))}{\Omega+\alpha u(\eta^{-1}(\tilde{t}))}=0.$

Since $\Omega+\alpha u(\eta^{-1}(\tilde{t}))\neq 0$ if $|\alpha|$ is
sufficiently small, we have $f=0$. Therefore, $\\{e^{in(\Omega
t+\alpha\int_{0}^{t}u(s)\,ds)}\\}_{n}$ is a complete basis of $L^{2}(0,T)$ by
Lemma 4.5. Thus, we conclude Corollary 4.3. ∎

Now we give a sufficient condition for $c_{n}\not=0$.

######  Theorem 4.6.

Assume that $g(\omega)$ is either uniform distribution or even and unimodal
$C^{2}$ function satisfying $g^{\prime\prime}(0)<0$ with the mean value
$\Omega=0$. Then, there exists $\varepsilon>0$ such that when
$K_{c}<K<K_{c}+\varepsilon$, the constants $c_{n}$ given in Theorem 4.1 are
nonzero for $n=0,1,\cdots$.

######  Corollary 4.7 (edge of bifurcation).

Under the same assumptions with Theorems 4.1 and 4.6, if $\Omega\neq 0$ and
$|\alpha|$ is sufficiently small, then there exists $\varepsilon>0$ such that
when $K_{c}<K<K_{c}+\varepsilon$, the order parameters $\\{r_{n}(\alpha
u,t)\\}_{n}$ for the system (3.2) gives a complete basis of
$L^{2}_{\mathrm{per}}(0,T)$ with $T\sim 2\pi/\Omega+O(\alpha)$.

Corollary 4.7 immediately follows from Corollary 4.3 and Theorem 4.6. Thus, we
can express any function in $L^{2}_{\mathrm{per}}(0,T)$ as a linear
combination of the order parameters $r_{n}(\alpha u,t)$ of the system (3.2).
According to (3.5), we can find numbers $\\{w_{n}\\}_{n}$ (i.e. readout
$W_{\text{out}}(\theta)$) to obtain any function $y(t)\in
L^{2}_{\mathrm{per}}(0,T)$.

In the proof of Theorem 4.6, we use the second mean value theorem for definite
integrals.

######  Lemma 4.8.

Let $a<b$. Suppose $f:[a,b]\to\mathbb{R}$ is integrable and
$\varphi:[a,b]\to\mathbb{R}$ is monotone and bounded. Then there exists
$\xi\in(a,b)$ such that

$\int_{a}^{b}f(x)\varphi(x)\,dx=\varphi(a)\int_{a}^{\xi}f(x)\,dx+\varphi(b)\int_{\xi}^{b}f(x)\,dx.$

###### Proof of Theorem 4.6.

It is trivial that $c_{0}(=r_{0})=1$ by the definition. Hence, we assume
$n\geq 1$. We prove the theorem only for the even and unimodal case, since the
proof for the uniform distribution is similar. Since we assume $K_{c}<K$, the
constant $c_{1}(=r_{*})$ is positive because of (2.5) ($r_{*}\sim
O(\sqrt{K-K_{c}})$). Thus, we assume $n\geq 2$. By (4.2) and (4.3), $c_{n}$ is
given by

$\displaystyle c_{n}$ $\displaystyle=a_{n}+b_{n}$
$\displaystyle=\int_{|\omega|<Kr_{*}}\cos\left(n\cdot\mathrm{Arcsin}\,\left(\frac{\omega}{Kr_{*}}\right)\right)g(\omega)\,d\omega$
$\displaystyle\qquad+\frac{2}{\pi}\int^{\infty}_{Kr_{*}}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos(n\theta)A(\omega,\theta)g(\omega)\,d\theta
d\omega.$

(for the uniform distribution $g(\omega)$, we can calculate this integral and
a proof is more easy).

First, $a_{n}$ is estimated as follows: By $\omega\mapsto Kr_{*}x$, we have

$a_{n}=2Kr_{*}\int^{1}_{0}\cos(n\cdot\mathrm{Arcsin}\,x)g(Kr_{*}x)dx.$

Due to Lemma 4.8, there exists $0<\xi<1$ such that

$a_{n}=2Kr_{*}g(0)\int^{\xi}_{0}\cos(n\cdot\mathrm{Arcsin}\,x)dx+2Kr_{*}g(Kr_{*})\int^{1}_{\xi}\cos(n\cdot\mathrm{Arcsin}\,x)dx.$

We estimate it around the bifurcation point ($K\sim K_{c}$, $Kr_{*}\sim 0$ and
$\xi\sim 1$) as

$\displaystyle a_{n}$
$\displaystyle=2Kr_{*}g(0)\int^{1}_{0}\cos(n\cdot\mathrm{Arcsin}\,x)dx$
$\displaystyle\qquad+(Kr_{*})^{3}(g^{\prime\prime}(0)+o(Kr_{*}))\cdot\int^{1}_{\xi}\cos(n\cdot\mathrm{Arcsin}\,x)dx,$

where we used $g^{\prime}(0)=0$ because $g$ is an even function. The first
integral is calculated as

$\int^{1}_{0}\cos(n\cdot\mathrm{Arcsin}\,x)dx=-\frac{\cos(n\pi/2)}{n^{2}-1},$

(I) Suppose that $n$ is an odd number. Then, the above value is zero. Further,
we have proved in the proof of Theorem 4.1 that $b_{n}=0$. Thus, we have

$c_{n}=a_{n}+b_{n}=(Kr_{*})^{3}(g^{\prime\prime}(0)+o(Kr_{*}))\cdot\int^{1}_{\xi}\cos(n\cdot\mathrm{Arcsin}\,x)dx.$

To estimate it, define

$F_{n}(\eta)=\int^{1}_{\xi}\cos(n\cdot\mathrm{Arcsin}\,x)dx=\int^{\pi/2}_{\eta}\cos(ny)\,\cos
y\,dy,$

where we put $x=\sin y,\,\xi=\sin\eta$. Note that $\xi=1$ corresponds to
$\eta=\pi/2$. Then, we can verify that
$F_{n}(\pi/2)=F_{n}^{\prime}(\pi/2)=F_{n}^{\prime\prime}(\pi/2)=0$ and

$F_{n}^{\prime\prime\prime}(\pi/2)=-2n\sin(\pi n/2)\neq 0\,\,(n:\text{odd}).$

This proves that $F_{n}(\eta)\neq 0$ for small $Kr_{*}>0$ and there exists
$\varepsilon$ such that $c_{n}\neq 0$ for $K_{c}<K<K_{c}+\varepsilon$. Note
that $\varepsilon$ is independent of $n$ because so is the factor
$g^{\prime\prime}(0)+o(Kr_{*})$.

(II) Suppose $n$ is an even number. In this case, $b_{n}\neq 0$ and

$\displaystyle c_{n}$
$\displaystyle=-2Kr_{*}g(0)\frac{\cos(n\pi/2)}{n^{2}-1}+(Kr_{*})^{3}(g^{\prime\prime}(0)+o(Kr_{*}))\cdot\int^{1}_{\xi}\cos(n\cdot\mathrm{Arcsin}\,x)dx$
$\displaystyle\qquad+\frac{2}{\pi}\int^{\infty}_{Kr_{*}}\int_{-\frac{\pi}{2}}^{\frac{\pi}{2}}\cos(n\theta)A(\omega,\theta)g(\omega)\,d\theta
d\omega.$ (4.4)

Let us estimate the third term $b_{n}$. We apply Lemma 4.8 with

$f(\theta):=\cos(n\theta)\quad\text{and}\quad\varphi(\theta):=\frac{2}{\pi}\int^{\infty}_{Kr_{*}}A(\omega,\theta)g(\omega)\,d\omega.$

Let us show that $\varphi(\theta)$ is monotonically increasing. It is easy to
see that the integral exists and $\partial A/\partial\theta>0$ for
$\omega>Kr_{*}$ and $-\pi/2<\theta<\pi/2$. This yields

$\varphi(\theta_{1})-\varphi(\theta_{2})=\frac{2}{\pi}\int^{\infty}_{Kr_{*}}(A(\omega,\theta_{1})-A(\omega,\theta_{2}))g(\omega)\,d\omega>0$

for any $-\pi/2\leq\theta_{2}<\theta_{1}\leq\pi/2$. Hence, Lemma 4.8 is
applicable to obtain

$\displaystyle b_{n}$
$\displaystyle=\varphi\left(-\frac{\pi}{2}\right)\int_{-\frac{\pi}{2}}^{\xi}\cos(n\theta)\,d\theta+\varphi\left(\frac{\pi}{2}\right)\int_{\xi}^{\frac{\pi}{2}}\cos(n\theta)\,d\theta$
$\displaystyle=\frac{1}{n}\sin(n\xi)\left(\varphi\left(-\frac{\pi}{2}\right)-\varphi\left(\frac{\pi}{2}\right)\right).$

To calculate it, note that

$A\left(\omega,-\frac{\pi}{2}\right)-A\left(\omega,\frac{\pi}{2}\right)=\frac{\sqrt{\omega^{2}-K^{2}r_{*}^{2}}}{\omega+Kr_{*}}-\frac{\sqrt{\omega^{2}-K^{2}r_{*}^{2}}}{\omega-
Kr_{*}}=-\frac{2Kr_{*}}{\sqrt{\omega^{2}-K^{2}r_{*}^{2}}}$

when $\omega>Kr_{*}$. This provides

$\displaystyle\varphi\left(-\frac{\pi}{2}\right)-\varphi\left(\frac{\pi}{2}\right)=-\frac{4Kr_{*}}{\pi}\int_{Kr_{*}}^{\infty}\frac{g(\omega)}{\sqrt{\omega^{2}-K^{2}r_{*}^{2}}}d\omega$
$\displaystyle=-\frac{4Kr_{*}}{\pi}\int_{Kr_{*}}^{Kr_{*}+\delta}\frac{g(\omega)}{\sqrt{\omega^{2}-K^{2}r_{*}^{2}}}d\omega-\frac{4Kr_{*}}{\pi}\int_{Kr_{*}+\delta}^{\infty}\frac{g(\omega)}{\sqrt{\omega^{2}-K^{2}r_{*}^{2}}}d\omega,$

where $\delta$ is a small positive number. The second term is of order
$O(Kr_{*})$ as $Kr_{*}\to 0$. For the first term, we can verify that it is
nonzero and $O(Kr_{*}\log(Kr_{*}))$ as $Kr_{*}\to 0$. Putting

$\varphi\left(-\frac{\pi}{2}\right)-\varphi\left(\frac{\pi}{2}\right)=B,$

we obtain

$b_{n}=B\frac{\sin(n\xi)}{n},\quad B\sim O(Kr_{*}\log(Kr_{*})),\,B\neq 0.$
(4.5)

Now we have estimations of all three terms of (4.4) for small $Kr_{*}>0$. Case
(II-1) If the third term is the leading term as $Kr_{*}\to 0$, $c_{n}\neq 0$
because $B\neq 0$. Case (II-2) If the first term is the leading term as
$Kr_{*}\to 0$ (it may happen, for instance, accidentally $\sin(n\xi)=0$), then
obviously $c_{n}\neq 0$. Case (II-3) If the second term is the leading term
(it happens when the first and the third terms cancel each other). Then, the
proof is the same as the case (I). For even $n$, we can verify that
$F_{n}(\pi/2)=F^{\prime}(\pi/2)=0$ and

$F_{n}^{\prime\prime}(\pi/2)=\cos(n\pi/2)\neq 0\,\,(n:\text{even}).$

Thus, we obtain $c_{n}\neq 0$.

This completes the proof of Theorem 4.6. ∎

In summary, we have obtained the following results to explain the edge of
bifurcation:

  * •

We need $K_{c}<K$ (after bifurcation) because when $K_{c}>K$, the trivial
steady state $r_{n}\equiv 0\,(n\neq 0)$ is stable.

  * •

In particular, just after the bifurcation $K_{c}<K<K_{c}+\varepsilon$, we
could prove $r_{n}\neq 0$.

  * •

We need $\Omega\neq 0$ because when $\Omega=0$, Corollary 4.3 does not hold.
This means that we need a bifurcation from the steady state to periodic state
(if $\Omega=0$, the bifurcation is from the steady state to the steady state).

  * •

Under the situation above, $\\{r_{n}(\alpha u,t)\\}_{n}$ consists of the
$L^{2}_{\mathrm{per}}$ basis and we can approximate any $L^{2}_{\mathrm{per}}$
function by using (3.5).

##  5 Experiments

Based on the above argument, we demonstrated the performance around the edge
of bifurcation in the Kuramoto reservoir computing model using numerical
simulation. Two tasks were used to evaluate the performance: a prediction task
and a transformation task.

Since the $n$-th order parameter $r_{n}(\alpha u;t)$ can be numerically
computed, we can obtain an approximate vector $\hat{w}_{\mathrm{out}}$ of
$\\{w_{n}\\}_{n=0}^{\infty}$ in the following steps: it is an implementable
model by the following procedure based on (3.5).

Algorithm. (Kuramoto reservoir computing)

Require: Data set (input/target)
$\\{u(t_{i}),y_{\mathrm{tar}}(t_{i})\\}_{i=1}^{T}$ and hyperparameter
$N\in\mathbb{N}$.

  1. 1.

For training time $t=t_{1},t_{2},\ldots,t_{T}$, solve

$\begin{split}(y_{\mathrm{tar}}(t_{1}),\cdots,y_{\mathrm{tar}}(t_{T}))&=(w_{0},w_{1},\cdots,w_{2N})\times\\\
&\begin{pmatrix}r_{0}(\alpha u;t_{1})&r_{0}(\alpha
u;t_{2})&\cdots&r_{0}(\alpha u;t_{T})\\\ \mathrm{Re}(r_{1}(\alpha
u;t_{1}))&\mathrm{Re}(r_{1}(\alpha u;t_{2}))&\cdots&\mathrm{Re}(r_{1}(\alpha
u;t_{T}))\\\ \mathrm{Im}(r_{1}(\alpha u;t_{1}))&\mathrm{Im}(r_{1}(\alpha
u;t_{2}))&\cdots&\mathrm{Im}(r_{1}(\alpha u;t_{T}))\\\ \vdots&&&\vdots\\\
\mathrm{Re}(r_{N}(\alpha u;t_{1}))&\mathrm{Re}(r_{N}(\alpha
u;t_{2}))&\cdots&\mathrm{Re}(r_{N}(\alpha u;t_{T}))\\\
\mathrm{Im}(r_{N}(\alpha u;t_{1}))&\mathrm{Im}(r_{N}(\alpha
u;t_{2}))&\cdots&\mathrm{Im}(r_{N}(\alpha u;t_{T}))\\\ \end{pmatrix}\\\
&=:\hat{w}_{\mathrm{out}}R.\end{split}$

(note that $r_{0}=1$ by the definition).

  2. 2.

Calculate
$\hat{w}_{\mathrm{out}}=(y_{\mathrm{tar}}(t_{1}),\cdots,y_{\mathrm{tar}}(t_{T}))R^{-1}$
(Moore-Penrose inverse).

  3. 3.

Return $y(t)=\hat{w}_{\mathrm{out}}(r_{0}(\alpha u;t),\mathrm{Re}(r_{1}(\alpha
u;t)),\mathrm{Im}(r_{1}(\alpha u;t)),\ldots,\mathrm{Im}(r_{N}(\alpha
u;t)))^{\top}$.

As above, we can indeed obtain the approximation $\hat{w}_{\mathrm{out}}$ of
$\\{w_{n}\\}_{n=0}^{\infty}$ (i.e. the real Fourier series of an integral
kernel $W_{\mathrm{out}}(\theta)$ in (3.4)) and use our reservoir computing
for tasks such as time series prediction.

###  5.1 Prediction task

In this task, reservoir computing predicted the waveform one step after the
input waveform. Fig. 3 (a) shows the composition of the Kuramoto reservoir
computing system given in (3.4). $N$ oscillators receive the input $\alpha
u(t)$(3.1), and $y(t)$ is the output by the product of
$\hat{w}_{\mathrm{out}}$ and $\\{r_{n}(\alpha u;t)\\}$(3.5). Note that
$\hat{w}_{\mathrm{out}}$ was trained using $y_{\mathrm{tar}}(t)$ during the
training period. The training period to obtain $\hat{w}_{\mathrm{out}}$ was
set to be 30 (s), and the testing period to the output $y(t)$ was set to be 30
(s). These training and tests were performed after 240 (s) to stabilize the
order parameters. In the prediction task, $u(t)$ was $\sin(\hat{\omega}t)$,
and $y_{\mathrm{tar}}(t)$ was $\sin(\hat{\omega}(t+\Delta t))$. The values of
these parameters are summarized in Table 1 below.

The $\omega$ of each oscillator in the Kuramoto model was obtained from the
$g(\omega)$ distribution. We simulated two cases; the Gaussian distribution
with mean $\Omega=0.5$ and standard deviation 0.4, or uniform distribution on
the interval $[0,1]$ (i.e. the mean value was also $\Omega=0.5$). When the
parameters were set in this way, the bifurcation points were
$K_{\mathrm{c}}\sim 0.638$ (Gaussian distribution), $0.637$ (uniform
distribution), respectively (Fig. 3 (b)). The reservoir computing performance
was quantified by the mean squared error represented as

$\frac{1}{T^{\prime}}\sum_{i=1}^{T^{\prime}}\|y_{\mathrm{tar}}(t_{i})-y(t_{i})\|^{2},$

where $T^{\prime}$ is the number of the test step. As a result, the Kuramoto
reservoir computing could predict the sine-wave one step after the input when
$K>K_{c}$ (Fig. 4 (a)). The mean squared error was large when $K<K_{c}$, and
the sine-wave was not well predicted, however, it became smaller when $K$
exceeded $K_{c}$ (Fig. 4 (b)). On the other hand, if $\Omega$=0, the mean
squared error remained large even when $K>K_{c}$. In comparison with the case
where $g(\omega)$ was uniformly distributed, for the Gaussian case, the mean
squared error still existed when $K$ was around $K_{c}$ (Fig. 4 (c)). This is
because, for the Gaussian, the order parameter continuously changes at
$K=K_{c}$, while for the uniform distribution, it suddenly jumps to the value
$r_{1}=\pi/4$ (see Fig. 2 and Fig. 3 (b)).

Figure 3:  The Kuramoto reservoir computing system. (a) Schematic diagram of the Kuramoto reservoir computing model and two types of tasks. (b) The relationship of the order parameter and the bifurcation with coupling strength $K$ when $g(\omega)$ is Gaussian distribution or uniform distribution in the numerical simulation. The solid line indicates the mean of 10 trials; the shade indicates the 95% confidence interval. Parameter |  Description |  Value
---|---|---
$N$ |  the number of oscillators |  500
$\Delta t$ |  time step |  0.01 (s)
$\alpha$ |  coefficient of the input |  0.01
$\hat{\omega}$ |  angular velocity of the input signal |  0.5 (rad/s)
$u(t)$ |  input signal |  $\sin(\hat{\omega}t)$
$y_{\mathrm{tar}}(t)$ |  target signal |  $\sin(\hat{\omega}(t+\Delta t))$
Table 1: Parameter settings. Figure 4:  The performance evaluation of the
Kuramoto reservoir computing in the prediction task. (a) The output $y(t)$
(blue line) and the target function $y_{\mathrm{tar}}(t)$ (orange line) when
the $\Omega=0.5$ with $K=0.60$ ($<K_{c}$), $K=0.63$ ($<K_{c}$), $K=0.64$
($>K_{c}$) and $K=0.70$ ($>K_{c}$). The bifurcation point is $K_{c}\sim
0.638$. The training was for 30 (s), shown in shading, followed by 240 (s) for
order parameter stabilization. And then a test period of 30 (s) during which
performance was evaluated. (b) The relationship between the mean squared error
and the coupling strength $K$ when $\Omega$=0 and 0.5. The solid line
indicates the mean of 10 trials; the shade indicates the 95% confidence
interval. (c) Comparison of Gaussian and uniform distribution of $g(\omega)$,
with setting $\Omega$ to 0.5.

###  5.2 Transformation task

In this task, the reservoir computing was evaluated for performance through
the conversion of the input waveform to another waveform. Here, the input
$u(t)$ was set to a sine-wave and $y_{\mathrm{tar}}(t)$ to a triangle wave
with the same period as the input signal (Fig. 5 (a)). Other parameters were
the same as those set in the prediction task.

When $\Omega=0.5$, the Kuramoto reservoir computing could convert the sine-
wave and generate the triangle wave when $K>K_{c}$ (Fig. 5 (a)). As with the
prediction task, the mean squared error was large when $K<K_{c}$, and the
triangle wave was not well generated, but it became smaller when $K$ exceeded
$K_{c}$ (Fig. 5 (b)). On the other hand, when $\Omega$=0, the mean squared
error remained large even when $K>K_{c}$. For the case where $g(\omega)$ was
uniformly distributed, as in the previous task, a difference in the mean
squared error appeared when $K$ was around $K_{c}$ (Fig. 5 (c)). The common
results obtained from these two tasks suggest that the Kuramoto reservoir can
approximate various types of target function $y_{\mathrm{tar}}(t)\in
L^{2}_{\mathrm{per}}(0,T)$ just after the bifurcation, which supports our
mathematical results.

Figure 5:  The performance evaluation of the Kuramoto reservoir computing in
the transformation task. (a) The output $y(t)$ (blue line) and the target
function $y_{\mathrm{tar}}(t)$ (orange line) with $\Omega=0.5$ for coupling
strengths: $K=0.60$ ($<K_{c}$), $K=0.63$ ($<K_{c}$), $K=0.64$ ($>K_{c}$) and
$K=0.70$ ($>K_{c}$). The training period is shown in shading. (b) The
relationship between the mean squared error and the coupling strength $K$ for
$\Omega$ values of 0 and 0.5. The solid line indicates the mean of 10 trials;
the shade indicates the 95% confidence interval. (c) Compares the impact of
Gaussian and uniform distributions of $g(\omega)$, with $\Omega$ set to 0.5.

## Acknowledgment

This work was supported by JST, CREST Grant Number JPMJCR2014, Japan (to HC
and KT), JST Moonshot R&D Program Grant Number JPMJMS2023, Japan (to TS).

## References

  * [1] Thomas L Carroll.  Do reservoir computers work best at the edge of chaos?  Chaos: An Interdisciplinary Journal of Nonlinear Science, 30(12), 2020.
  * [2] Hayato Chiba.  A proof of the Kuramoto conjecture for a bifurcation structure of the infinite-dimensional Kuramoto model.  Ergodic Theory and Dynamical Systems, 35(3):762–834, 2015.
  * [3] James P Crutchfield and Karl Young.  Computation at the onset of chaos.  Entropy, complexity, and the physics of information, 8:223, 1990\.
  * [4] Ken Goto, Kohei Nakajima, and Hirofumi Notsu.  Twin vortex computer in fluid flow.  New Journal of Physics, 23(6):063051, 2021.
  * [5] Herbert Jaeger.  The “echo state” approach to analysing and training recurrent neural networks-with an erratum note.  Bonn, Germany: German National Research Center for Information Technology GMD Technical Report, 148(34):13, 2001.
  * [6] Y. Kuramoto.  Self-entrainment of a population of coupled non-linear oscillators.  In International Symposium on Mathematical Problems in Theoretical Physics, volume 39, pages 420–422. Springer, Berlin, 1975.
  * [7] Y. Kuramoto.  Chemical oscillations, waves, and turbulence.  Springer Series in Synergetics, 19. Springer-Verlag, Berlin, 1984.
  * [8] Chris G Langton.  Computation at the edge of chaos: Phase transitions and emergent computation.  Physica D: nonlinear phenomena, 42(1-3):12–37, 1990.
  * [9] Mantas Lukoševičius and Herbert Jaeger.  Reservoir computing approaches to recurrent neural network training.  Computer science review, 3(3):127–149, 2009.
  * [10] Thomas Natschläger, Wolfgang Maass, and Henry Markram.  The “liquid computer”: A novel strategy for real-time computing on time series.  Telematik, 8(ARTICLE):39–43, 2002.
  * [11] Norman H Packard.  Adaptation toward the edge of chaos.  Dynamic patterns in complex systems, 212:293–301, 1988.
  * [12] Steven H Strogatz.  From Kuramoto to Crawford: exploring the onset of synchronization in populations of coupled oscillators.  Physica D: Nonlinear Phenomena, 143(1-4):1–20, 2000.
  * [13] Gouhei Tanaka, Toshiyuki Yamane, Jean Benoit Héroux, Ryosho Nakane, Naoki Kanazawa, Seiji Takeda, Hidetoshi Numata, Daiju Nakano, and Akira Hirose.  Recent advances in physical reservoir computing: A review.  Neural Networks, 115:100–123, 2019.
  * [14] Liang Wang, Huawei Fan, Jinghua Xiao, Yueheng Lan, and Xingang Wang.  Criticality in reservoir computer of coupled phase oscillators.  Phys. Rev. E, 105:L052201, May 2022.
  * [15] Toshiyuki Yamane, Seiji Takeda, Daiju Nakano, Gouhei Tanaka, Ryosho Nakane, Shigeru Nakagawa, and Akira Hirose.  Dynamics of reservoir computing at the edge of stability.  In Neural Information Processing, pages 205–212, Cham, 2016. Springer International Publishing.
  * [16] Min Yan, Can Huang, Peter Bienstman, Peter Tino, Wei Lin, and Jie Sun.  Emerging opportunities and challenges for the future of reservoir computing.  Nature Communications, 15:2056, 03 2024.
  * [17] Zhihao Zuo, Zhongxue Gan, Yuchuan Fan, Vjaceslavs Bobrovs, Xiaodan Pang, and Oskars Ozolins.  Self-evolutionary reservoir computer based on Kuramoto model.  arXiv preprint arXiv:2301.10654, 2023.