Hardware optimization for photonic time-delay reservoir computer dynamics

Reservoir computing (RC) is one kind of neuromorphic computing mainly applied to process sequential data such as time-dependent signals. In this paper, the bifurcation diagram of a photonic time-delay RC system is thoroughly studied, and a method of bifurcation dynamics guided hardware hyperparameter optimization is presented. The time-evolution equation expressed by the photonic hardware parameters is established while the intrinsic dynamics of the photonic RC system is quantitively studied. Bifurcation dynamics based hyperparameter optimization offers a simple yet effective approach in hardware setting optimization that aims to reduce the complexity and time in hardware adjustment. Three benchmark tasks, nonlinear channel equalization (NCE), nonlinear auto regressive moving average with 10th order time lag (NARMA10) and Santa Fe laser time-series prediction tasks are implemented on the photonic delay-line RC using bifurcation dynamics guided hardware optimization. The experimental results of these benchmark tasks achieved overall good agreement with the simulated bifurcation dynamics modeling results.


Introduction
Reservoir computing (RC), a subset of the recurrent neural network has attracted considerable attentions due to its simplicity in network training. A typical RC network consists of three layers: an input layer, a reservoir layer (hidden layer) and an output layer. The neuron connections in the reservoir layer are sparse and random but fixed while the training takes place only at the output layer by computing the output weight matrix. RC is a simple yet powerful machine learning approach and particularly effective processing sequential data in classification or prediction tasks such as time series forecasting [1], pattern classification [2][3][4], spoken digit recognition [5], signal recovery in communication systems [6,7].
Recently, driven by edge computing applications, RC algorithms have been explored to harness on various hardware platforms such as electronic RC by nonlinear analog electronic circuits [8,9] and FPGAs [10], reservoir by spintronics [11,12], photonic RC by optical node arrays [5,13] and photonic time-delay systems [14][15][16][17]. For many physical RCs, the inherent dynamics of physical devices are utilized to form random neuron connections of the reservoir. Photonic RC that has rapid signal processing speed and chip-scale integrated capability is one viable solution in developing small size computing unit for edge computing platform. Integrated photonic RC has been realized on an on-chip micro-ring or micro-cavity [18,19] and interconnected-node mesh reservoir [5,13]. Photonic delay line-based RC using virtual nodes in a temporal space with simple ring topology was first demonstrated experimentally in an optoelectronic (OE) oscillator configuration by Paquot et al [14] and Larger et al [20]. Other versions of photonic delay-line RC implementation have also been reported such as all-optical feedback loop [15] and laser cavity with feedback [16]. These reported photonic delay-line RCs are appealing as they can be constructed by off-the-shelf photonic components.
It is well known that the RC network inference result is strongly influenced by the hyperparameters of the RC network. Grid search [17,21], Bayesian optimization [22][23][24] and entropy-based hyperparameter optimization have been reported [25]. Grid search for hyperparameter optimization is commonly used to determine the lowest inference error in a global setting, but this approach often suffers from exceedingly long computation time particularly when there are multiple hyperparameters. Comparing to grid search, Bayesian optimization utilizing probability theory can greatly reduce the number of iterations.
Hyperparameter optimization for physical RC poses unique challenge as it requires the tunability and reconfiguration of a physical system. In this work, we report a practical and simple method that utilizes bifurcation diagrams for hyperparameter optimization. We first establish a time evolution equation of the optical oscillator that closely matches the physical setting of the delay-line reservoir. Next, we model the derived evolution equation to determine the dynamic regimes of the RC system on a bifurcation diagram. We find that the low inference errors are generally obtained near bifurcation points [26] when the system is at the edge of stability just before entering the oscillatory state, which allows us to obtain a much smaller range of physical settings to be tuned for RC hyperparameter optimization. We implemented this approach on three benchmark tasks, i.e. nonlinear channel equalization (NCE), nonlinear auto regressive moving average with the 10th order time lag (NARMA10) and Santa Fe time-series prediction tasks.
The paper will be organized as follows. In section 2, general background on time-delay RC network and modeling of our photonic RC implementation are introduced. Detailed derivation of hardware-based time-evolution equation is also provided in this section. Section 3 introduces the intrinsic dynamics of delay-line RC system while section 4 applies the bifurcation-dynamic guided hyperparameter optimization on three RC benchmark tasks.

Photonic true-time-delay RC network
The photonic time-delay RC network studied in this work uses a cyclic topology similar to the setup in Paquot et al [14]. The input layer signals are sequentially injected to the reservoir layer via a single physical nonlinear node while the signals in the reservoir circulate on a closed delay-line loop. The internal reservoir states are often referred to as the virtual nodes traveling along the time-delay line in a series of optical pulse trains with varying intensity to form a spatial-temporal dynamic system. A good reservoir should have high dimension, nonlinearity and fading memory [27,28]. The RC network of this work was implemented on a bench-top photonic time-delay system constructed with commercially available optical and electrical components. A schematic of the photonic delay-line RC is sketched in figure 1 and the components used in the RC construction are listed in supplementary materials.
The reservoir in this work is an optical fiber oscillator that operates at the edge of stability. The capacity of the reservoir is denoted by the virtual node number N = τ /θ, where τ is the round-trip delay time and θ is the node separation in unit of time. In this work, the optical signal is coupled to the reservoir by an arbitrary waveform generator (AWG) at a sampling rate of N/τ . A reservoir connectivity matrix, W res is used to characterize the node connections within the reservoir. A node mismatch k is defined as k = (τ − τ ′ ) /θ where τ ′ denotes the sample hold time for each input signal. For synchronized operation where τ = τ ′ , W res is simply an identity matrix. To enrich the reservoir dynamics with increased interactions among virtual nodes, RC can be tuned to asynchronized regime [14] in which W res is a shifted identity matrix. In this work, we utilized a small node mismatch k = 2 for NCE and Santa Fe tasks and the W res matrix non-zero diagonal element '1' is shifted by 2 columns to the right as shown below: The input signal u (t) is acquired externally and then coupled to the physical RC via a sample-hold operation. A piecewise mask function m (t) is multiplied with u (t) to increase the diversity of the projected input data in the reservoir [14]. The discrete levels of m (t) are randomly generated in the range of [−1, 1] with a time duration of θ for each level and raw input data are normalized to u (t) ∈ [−1, 1]. The period of mask function is set equal to the duration time of each u (t) level τ ′ . In this work, the input layer is  implemented digitally in an external computer and the operation of m (t) · u (t) is carried on in MATLAB program. Duport et al [29] reported using cascade Mach-Zehnder modulators (MZMs) to perform the multiplication in optical domain so a fully hardware implementation of the input layer is feasible. The masked input data sequence m (t) · u (t) is fed into the reservoir layer through the AWG. The amplitude of m (t) · u (t) at the AWG output is to be adjusted as a key parameter closely related to β.
In the reservoir layer, a LiNbO 3 MZM serves as the single physical node to provide a nonlinear response to the encoded input optical signals. The MZM has a marked half-wave voltage (V π ) of 3.35 V (DC) and 3.7 V (RF). The DC operating point of the MZM is set by an external bias V DC . The measured transfer function of the MZM is plotted in figure 2(a) and a built-in phase offset φ os = 0.24π is measured at V DC = 0. An optical splitter is used at the MZM output to tap out 10% of the optical power to a readout photodetector (PD) in the output layer. An electrical amplifier is used to compensate the optical loss in the delay line as well as provide an easy way to adjust the input gain β and the feedback gain α in addition to the optical attenuator in the delay line loop. The electrical amplifier can be eliminated from the RC hardware setup when the laser power is sufficiently large or the loss in the delay line loop is low. The elimination of the electrical amplifier helps to reduce the power consumption of the RC network as well as reduce thermal noise in the system. The output intensity I (t) of the MZM with RF port is fed with V (t) can be expressed as where I 0 is the MZM maximum output power with insertion loss and φ 0 is defined as The RF signal V (t) to the MZM comprises the summation of feedback signals after a round trip in the delay line loop and the newly added input u (t) at the combiner and it can be expressed as where η O = I fb (t) /I (t) is the feedback PD optical power retaining ratio, η E is the total gain/loss of the electrical path (green part in figure 1), R (unit: V/W) is the multiplication of the feedback PD responsivity with the PD's effective resistance, V AWG is the maximum amplitude of AWG output and I * (t − τ ) denotes the feedback optical intensity as a function of time.
The output photocurrent (I ph ) of both the feedback PD and the readout PD is proportional to the PD input optical intensity while I ph of the feedback PD multiplied with the input impedance of the oscilloscope (50 ohm) yields a voltage that corresponds to the reservoir state value and is displayed at the oscilloscope. For future work of an integrated optical delay-line RC, a transimpedance amplifier can be employed to produce the reservoir states. In Paquot's work [14], the output signal intensity I (t) was normalized following 1] under the assumption that the reference of the output intensity (x = 0) is set at I DC = I 0 /2 which is different from how the reservoir states are read out from the oscilloscope in our work. Here, we utilized an AC coupled approach to obtain both positive and negative reservoir state. The feedback PD is connected to a high pass filter with a cut-off frequency of 30 kHz so that the DC component of the PD output is blocked. Impedance matched transmission line ensures the AC signals are efficiently delivered to record reservoir states, where a positive/negative AC component of PD output corresponds to positive/negative reservoir state reading.
For illustration, an example is shown in figure 2(b) where V DC is set at 4.9 V, corresponding to a 0.3π phase difference from the peak transmission. This gives a reference of output intensity set at I DC = I 0 cos 2 (φ 0 /2) = 0.8I 0 without any signal fed to the RF port. The instant MZM output I (t) determined by the V DC + V (t) can then give either a positive (above red line) or a negative (below red line) reservoir state value using I DC as the reference. Normalizing to I 0 , the reservoir state can be expressed as Here  (4), the time evolution equation can be expressed using hardware parameters as Note the coefficient of x (t − τ ) and u (t) of equation (5) marks the delay-line loop feedback strength α and input signal strength β respectively with the following definition: Combining the phase terms yields A transfer function of the MZM with respect to the hyperparameter φ defined by equation (8) is plotted in figure 2(c). The general form of a sinusoidal time evolution equation can then be expressed as where α, β, φ are three hyperparameters of this RC network to be optimized based on the computation tasks. Replacing the time variable t by index n, a discretized time evolution equation can be derived for convenience of theoretical studies: where n represents the n th input data and i denotes the i th virtual node on the delay line.
In the output layer, a one-step training is performed by calculating the output matrix W out via ridge regression where Y is the targeted output matrix consisting of targeted outputs y (n), and X is the reservoir states matrix concatenating the projected reservoir state vectors of each input data point. The superscript T represents the matrix transpose operation and I denotes the identity matrix. λ is a regularization parameter to prevent overfitting and is chosen according to the virtual node number N ′ and the computation task.

Delay-line RC system dynamics
The dynamics of OE time-delay systems have been investigated in several works [30][31][32]. Passive optical systems subject to delayed feedback can be described by Ikeda delay differential equation (DDE) or slightly modified forms [33] to model the nonlinearity [34]. The dynamics of the delay-line RC used in this work can be described as where the dimensionless dynamical variable x represents the normalized MZM output power and α, φ follows the same notation as discussed in section 2. Though the detailed form of the DDEs may vary due to system configuration variations, two parameters are common for OE delayed systems in characterizing the dynamics: relaxation time τ r and delay time τ . The parameter τ r in equation (12) characterizes how fast the system reacts to the input changes and is limited by the system cutoff frequency f H . In our RC network, the slowest component is the feedback PD (and the readout PD) with a marked bandwidth of 1 GHz. The estimated τ r of the RC is τ r = 1/2πf H ≈ 0.16 ns [16] which is several orders of magnitude smaller than the oscillator cavity delay time τ . For the case of τ r /τ ≪ 1, the derivative term x ′ (t) in equation (12) can be neglected as suggested in [33]. In some delayed dynamical systems [8,10], τ r is comparable to the node separation θ to incite interactions between the neighboring nodes to further boost the reservoir dynamics. In this work, we set τ r ≪ θ and τ r ≪ τ so neighbor nodes interaction is negligible. Consider a discrete-time system with index n, equation (12) can then be simplified to Equation (13) shows an explicit mapping of the current dynamical state x (n) with the past state x (n − 1), which agrees with the time evolution equation (10) when β = 0. The stability of the system can be revealed from a Hopf bifurcation diagram. In general, the system is stable for small α and oscillatory or chaotic when α value is large.
We define α H as the α value at the onset of oscillatory state on the Hopf bifurcation diagram. When the system has no input feed, i.e. β = 0, varying φ from 0 to 2π, α H versus φ is plotted in figure 3 given the initial state x (0) = 1. Above the α H curve, the system is oscillatory or chaotic, and below it, the system is stable. The reservoir layer, utilizing the system dynamics, projects the input signals to a higher dimensional space to increase separation ability of the inputs signals. Dimensionality expansion is the underline physics of RC inference. Richer dynamics of the reservoir gives better RC inference results. Hence, the hyperparameters should be chosen so that the RC operates near the bifurcation point where the system has the richest dynamics while still in the stable regime. Below the α H (φ) curve, we marked a region in cyan bounded by α H (φ) − ∆ (red dash line) where ∆ is a small value to visualize a set of ideal hyperparameter combinations on the φ − α plane.
The α H (φ) curve in figure 3 reveals the intrinsic dynamics of the system without input signals, i.e. β = 0. As shown in equations (9) and (10), the input signal m (t) · u (t) multiplied with the input gain β produces an extra phase to the MZM cosine transfer function. When the input data is normalized to [−1, 1], the input layer signals, i.e. m (t) · u (t) projected to the reservoir layer at the combiner is proportional to [-β, β]. If β is too small, the extra phase induced by β will cause the reservoir states to only vary in a small portion of the sinusoidally varying curve of the nonlinear node. For example, when the MZM is biased near its quadrature point, we can only obtain a near linear response of the reservoir with small β, so the reservoir layer acts more like a linear classifier rather than a complex dynamic system. On the other hand, if β is large, the feedback signals produced at the PD become insignificant compared to the input at the combiner, which diminishes the memory effect of the RC. β is mainly tuned by adjusting the amplitude of the AWG output and in general has minor impact on the shape of optimized regions on φ − α plane. However, proper choice of β can give us lower error rates for different RC tasks and more relaxed restriction to α and φ conditions.

Benchmark tasks using bifurcation dynamics guided hyperparameter optimization
For the photonic delay-line RC of this study, the feedback gain (α), input gain (β) and bias phase offset of the modulator (φ) are the RC hyperparameters to be adjusted, though the RC inference results are not equally sensitive to all these hyperparameters. The feedback gain and the phase value are more crucial in determining the RC network performance. Based on the theoretical study of α − φ dynamic regime while choosing the proper input gain β, the hyperparameter optimization can be simplified to a 2D grid search. Furthermore, using the simulated bifurcation diagram based on hardware parameters extracted to construct the evolution equation, there is no need to run a full space 2D optimization of α and φ while a much reduced α − φ grid range would produce good results. This bifurcation dynamic guided optimization utilizing mathematic modeling of hardware RC leads to a significant time reduction in hardware setting optimization.
We tested our method of bifurcation dynamic guided hyperparameter optimization on three commonly studied benchmark tasks: NCE, NARMA10 and Santa Fe tasks. For the RC numerical study based on equation (10), we perform a fine 2D grid search for α and φ in which φ and α are scanned in a range of 0 to 2π with a step of 0.02π and 0 to 4 with a step of 0.02, respectively. This mathematic modeling is carried on in MATLAB with fine resolution to achieve the best RC inference results. β is set at different fixed values for three benchmark tasks and more discussion of β effect on RC performance will be provided in section 4(d).
To examine how closely the experimental testing results match the simulated results, we also run a 2D grid search for α and φ configurations on the photonic delay line RC. The hyperparameter φ, tuned by changing the value of MZM DC bias voltage, is tested with a step of π/6 for each setting. Hyperparameter α and β are experimentally adjusted by tuning the optical attenuator and AWG maximum output amplitude, respectively. For the photonic RC, the reservoir enters either an oscillatory or chaotic state at large α so the photonic RC hardware settings corresponding to large α are excluded from the experimental tests. For very small α value, close to zero, the photonic RC performance is limited by the system noise originating from analog-to-digital conversion (ADC) of oscilloscope, laser relaxation oscillations, and PD dark current and shot noise, etc [17,32], so α with small values (< 0.05) are also excluded from the tests. In general, we choose α with a small interval (∼ 0.05) near chaotic states and a larger interval as α decreases. The total number of tested (φ, α) settings is 80 − 100 for each task, sufficient to distinguish the error rate difference caused by varying the hyperparameter setting. Detailed experimental settings of three tasks are listed in table 1.
For our current experimental setup, optical attenuator and DC bias voltage are tuned manually, making it time consuming to adjust the hardware for all possible combinations of φ and α values, therefore a much larger grid size is utilized compared with numerical simulation. The theoretical modeling results are plotted  in 2D color contour maps where 'islands' of good inference results are identified. Based on the optimal values of φ and α obtained by theoretical modeling, we identify four combinations of (φ, α) in the photonic delay-line RC setting that is closest to simulation. An empirical factor 0.6 is applied in determining the optimal α and scaled in the 2D plots of the experimental results. The overall optimization process is shown through the flow diagram in figure 4.

NCE task
NCE task is widely studied in the RC community [14,15,35,36]. The input d (n) is randomly generated from the set {-3, −1, 1, 3} following uniform distribution. The nonlinear channel provides a transformation following where υ (n) is the i.i.d. Gaussian noise term with zero mean. In this work, the signal-to-noise ratio is set to be 28 dB and a dataset length of 3000, a training length of 1200, and a testing length of 1600 is used. The symbol error rate (SER) of the 2D RC simulated result is plotted in a color contour map with varying α and φ, shown in figure 5(a). Each contour line corresponds to an SER reading while the filled-in colored area represents the level of SER. The SER contour plot is close to symmetric about φ = π. SER lower than 0.01 is widely achieved over a large φ range near two quadrature points for α ∈ [0.25, 1]. SER lower than 0.001 is found to form a small 'island' near φ = π and a larger 'island' near φ = 2π (0), both marked in dark blue. Numerical modeling of the RC gives SER < 6 × 10 −4 at φ = 1.82π, α = 0.64 marked by red 'x' mark in figure 5(a).
The α − φ grid search was performed on the photonic delay-line RC and the SER results are shown in figure 5(b). Two areas of SER < 6 × 10 −4 are identified in the experimental grid search: one is between φ = 0 and φ = 0.5π and the other one between φ = 1.5π and φ = 2π. The achieved SER is on the same order as Paquot et al [14] that reported SER of 1.3 × 10 −4 and we expect the upper bound of our best SER to be smaller if more data points are utilized for testing. These low SER areas in the photonic RC testing are closely matched with the RC numerical modeling of equation (10). Due to φ being continuous with a period of 2π, these two areas are indeed connected and correspond to the quadrature point φ = 0 (2π). The general trend of the SER by photonic RC matches well with the RC by numerical modeling. Four closest grid points to the optimal (φ, α) were identified, marked by green 'x' marks in figure 5(b) on the photonic RC 2D SER map and summarized in table 2. Red 'x' mark represents the optimal (φ, α) condition suggested by numerical study and four closest photonic RC hardware settings are marked by green 'x' marks. To examine how the numerical RC matches with the photonic RC, we performed 2D α-φ grid search by both numerical and experimental approach. In practice, the photonic RC hyperparameter optimization can be reduced to examining those four (φ, α) settings closely matching the bifurcation dynamics analysis to set the hardware parameters of the photonic RC. This would allow significant time reduction in tedious hardware adjusting.

NARMA10 task
The emulation of the nonlinear auto regressive moving average model, i.e. NARMA10 task, is a commonly used benchmark prediction task for nonlinear dynamic systems. It was first introduced in [1] and studied as a benchmark prediction task in [5,8]. The discrete time input u (n) and targeted output y (n) follows the dynamic y (n + 1) = 0.3y (n) + 0.05y (n) · 9 ∑ i=0 y (n − i) + 1.5u (n) u (n − 9) + 0.1 (15) where u (n) follows a uniform distribution over the range of [0, 0.5] and is randomly generated. The goal of this task is to predict the output y (n) given the input u (n). A dataset with a length of 3000 is used to run the task, with a training length of 1000, testing length of 1800 and first 200 datapoints discarded. The RC performance of NARMA10 task is characterized by the normalized mean square error (NMSE) [29]. The color map of NMSE with respect to φ, α by numerically modeled RC is shown in figure 6(a). There are five 'islands' that correspond to low NMSE on the bifurcation diagram. A large 'island' of ideal (φ, α) combinations was found near φ = π, just below the system bifurcation, on which the lowest NMSE by the RC numerical simulator is found to be 0.0256 at φ = 0.88π, α = 0.92, marked by red 'x' mark as shown in figure 6(a). Those two isolated dark blue islands near α = 2 are not experimentally measurable due to chaotic response of the photonic RC system. Those two low NMSE areas at bottom left and right of the bifurcation diagram (α ∼ 1), near symmetric and connected at φ = 0 and 2π considering an unfolded diagram, are also ideal operation regimes for NARMA10 task. If we consider NMSE ⩽ 0.2 acceptable, both lower 'islands' near two quadrature points with α ∼ 1 meet the criterion while the 2D photonic grid search matches well with the numerical RC modeling.  The color map trends of experimental results shown in figure 6(b) overall matches the numerical RC simulator results, where the low NMSE [dark blue in figure 6(b)] area is also found near two quadrature points and close to bifurcation. Those 4 closet grid points on the photonic RC to the optimal (φ, α) by numerical RC simulator were marked by green 'x' marks on Fib. 5(b) the calculated NMSE results are summarized in table 3. Using photonic RC, the lowest NMSE is 0.0930 corresponding to setting 2 on table 3, better than Paquot's [14] result NMSE = 0.168 obtained from similar experimental setup as ours. Photonic RC testing overall yields higher prediction error rate than the RC simulator mainly because of the quantization noise in the ADC process at the output layer [17]. Both simulated and experimental results show that peak and valley of MZM transmission should be avoided.

Santa Fe task
The Santa Fe laser time-series prediction task is another well studied RC benchmark task [9,18]. We performed one-step ahead prediction using a dataset with 4100 data points where the first 100 data points are discarded, the next 3000 are used for training, and the remaining 1000 data points are used for testing.
The 2D NMSE plots by numerical RC simulator and photonic RC are plotted in figure 7, both indicating higher prediction error near 0.5π and 1.5π at low α values. Compared to the NCE and NARMA10 tasks, Santa Fe task is less sensitive to hyperparameter settings as suggested by both numerical RC simulator and the photonic RC. The simulated contour map in figure 7(a) shows a band of all φ values near α ∼ 0.8 that gives NMSE < 0.01. Simulation reveals a lowest NMSE of 0.0054 at φ = 1.36π, α = 1.1. Four closely matched hyperparameter settings on the photonic RC are identified, two of which correspond to the modulator transmission valley so not considered. The other two matched settings are summarized in table 4. Test setting 1 gives a NMSE of 0.0341, which is 6.3 times higher than the RC simulator result mainly due to the quantization noise similar as NARMA10 task. Both settings result in lower NMSE than Larger's result of NMSE = 0.124 that used similar photonic RC setup.

Hyperparameter β effect
In this work, we performed a simple analysis of the effect of β to the RC system on the numerical RC simulator. We scanned β from 0.5 to 2.5 with a relatively large step of 0.1. The choice of lower end of β is due to the concern of photonic RC noise, while the AWG output amplitude level determines the highest β. A 3D grid hyperparameter scan with lower resolution of α and φ is performed and the lowest SER or NMSE varying with β for those three benchmark tasks are plotted in figure 8. In the tested range of β, the NMSE   only changes for a factor 2.4 and 3 for the NARMA 10 task and Santa Fe task, respectively. As a comparison, NMSE can vary several orders of magnitude with (φ, α) choices. For the NCE task, SER < 10 −3 is obtained for β < 1. In our photonic RC test, β = 0.7 was set that gives SER = 0. For the NARMA10 task, a valley of low NMSE values is obtained for β in the range of 0.7-1.1. In the photonic RC test, β = 1 was selected as larger β is preferred for relatively higher signal-to-noise ratio. For the Santa Fe task, the NMSE decreases monotonically with β that shows a different trend as the NCE and NARMA 10 task. During the Santa Fe task test, due to an equipment failure, we changed to a different AWG (Tektronix 70 000) that has a maximum output voltage V AWG = 0.25V, corresponding to β = 1.8. We expect reduced NMSE for the Santa Fe task if a higher AWG output voltage is available.

Conclusion
In this work, we explored a practical engineering approach in selecting proper hyperparameter values of a photonic delay-line RC through bifurcation dynamic analysis. The time-evolution equation of the photonic RC is derived and modeled numerically based on the hardware settings. The intrinsic dynamics and bifurcation behaviors of the time-delay system are analyzed to give insight on how each hyperparameter affects the RC dynamics, so a reduced range of grid search is possible in determining the optimal operation conditions of a RC. It shows that a good RC performance can be obtained by tuning the delay-line system close to the bifurcation point while operating in the steady-state regime. The low-error-rate region shows good alignment with the α H (φ) curve near quadrature points in both the numerical study and the experimental hardware testing for all three benchmark tasks. A general good agreement in the 2D bifurcation diagram is obtained between the numerical model and the photonic RC experimental test. Among the three hyperparameters of the photonic RC, delay-line loop feedback strength α, modulator phase φ, and input signal strength β, the photonic RC appears to be more sensitive to the choice of α and φ so a 2D grid optimization for α and φ was performed at a fixed β value chosen from both empirical as well as a simple numerical analysis. Using the numerical RC modeling results, 4 experimental settings of (α, φ) are exanimated on the photonic RC and assessed by NMSE or SER. Good agreement is obtained between the photonic RC test results and the numerical simulation.

Data availability statement
The data cannot be made publicly available upon publication because they are not available in a format that is sufficiently accessible or reusable by other researchers. The data that support the findings of this study are available upon reasonable request from the authors.