Semiconductor lasers subjected to different forms of external perturbations produce rich nonlinear dynamics [1–9]. The chaotic dynamic properties of semiconductor lasers have been extensively studied in the past decades for their potential applications in secure communication, chaotic radar, and random bit generators [10–16]. To unveil the chaotic dynamics of semiconductor lasers, classical methods usually rely on the nonlinear dynamical equations derived from physical principles [17–20]. There are also some mathematical models that reconstruct the dynamical state from chaotic time series directly [21–23]. However, the application of these models may cause substantial errors due to the imperfect modeling of crucial dynamics and limitation of system dimension.

Recent advancements in the field of machine learning have made great progress in the prediction of dynamical systems purely from past measured data of chaotic time series. In the past decades, different machine learning approaches have been adopted to realize chaotic dynamics prediction [24]. Among these machine learning approaches being used so far, the reservoir computing has proven to perform superior prediction with less computational demand [25–28]. As a specific neuromorphic computing model that deviated from the von-Neumann architecture, reservoir computing imitates the neural information processing systems and consists of three layers, namely, the input layer (sensing), the reservoir layer (processing), and the output layer (control). The input and reservoir layers of reservoir computing are fixed. Only the readout coefficients of the output layer need to be trained in an optimized way, which leads to faster and more stable training and prediction [29]. Thus, reservoir computing has been widely employed to forecast the spatiotemporal chaotic dynamics, such as the chaotic time series of semiconductor lasers. Prediction of chaotic dynamics of an optically injected laser is discussed systematically with several machine learning algorithms [deep learning, support vector machine (SVM), nearest neighbors, and reservoir computing (RC)] [30]. Reservoir computing can also be seen as a state observer to infer the dynamical variables of semiconductor lasers. After training, only one of the dynamical variables is required to reconstruct the other two variables, if the algorithm is properly trained with all three variables [31]. Additionally, it is numerically demonstrated that reservoir computing can successfully predict the continuous intensity time series and reproduce the underlying chaotic dynamical behaviors of semiconductor lasers, such as power spectrum, PDF, and attractor [32]. The automated genetic algorithm optimization approach is also introduced to efficiently select the optimal reservoir parameters and achieve accurate forecasting [33]. Furthermore, recent advances in time-delayed reservoir computing (TDRC) have made it possible to physically implement the forecasting for the chaotic dynamics of distributed feedback laser (DFB) with optical feedback, mutually coupled VCSELs, and other spatiotemporal chaotic dynamics [34–46]. On the other hand, the traditional reservoir computing is a data-driven method and makes little use of the dynamical equations that generate the training data. Therefore, it may require an additional size of training data and computational resources to achieve accurate prediction, especially when the spatiotemporally chaotic system is large and complex. While in many realistic applications scenario, some underlying knowledge about the system is available, it could facilitate the performance of reservoir computing effectively [47–49].

In this paper, we numerically and experimentally investigate the heterogeneous forecasting scheme for chaotic dynamics generated by VCSELs with knowledge-based photonic reservoir computing. With the combination of the spin-flip model and photonic TDRC, we elucidate that the heterogeneous forecasting scheme could exhibit better prediction with a reduced training dataset compared to the situation when these two methods are conducted individually. Moreover, the prediction is robust to the deficiencies in the physical model. The influence of operating parameters and the complexity of chaotic dynamics in VCSELs on the performance of knowledge-based reservoir computing is investigated systematically. Finally, the generality of our results is validated experimentally using a TDRC based on semiconductor lasers.

Then the knowledge-based dynamical model of the training data is established based on the SMF with parameter and initial condition mismatch, as follows: I′(t+Δt)=F[(1+η)I(t)],(7)where I(t) and its corresponding dynamical state [Evx(t),Evy(t),Nv(t),nv(t)] at time t with feedback strength γf=25ns−1 are derived from Eqs. (1)–(4). And I′(t+Δt) is the one-step-ahead prediction data evolving from the imperfect dynamical-based model with deviated feedback strength γf1=(1+η1)γf or dynamical state mismatch [Evx′(t),Evy′(t),Nv′(t),nv′(t)]=[(1+η2)Evx(t),(1+η2)Evy(t),(1+η2)Nv(t),(1+η2)nv(t)] [53,54]. The error parameters η1 and η2 are used to quantify the deviation between the knowledge-based dynamical model and the “true” dynamical system. However, it is noteworthy that the accuracy of this one-step-ahead prediction I′(t+Δt) by the knowledge-based dynamical model is limited and may induce substantial errors compared to the “true” data I(t+Δt).

As shown in Fig. 1, the knowledge-based photonic reservoir computing is constructed by an optically injected semiconductor laser under optical feedback which consists of three layers.

In the input layer: The time series of training data I(t) and the output of the imperfect knowledge-based dynamical model I′(t+Δt) are concatenated to reconstruct the input signal u(t), thereby mitigating the issue of insufficient RC data driven solely by pure data. Then, u(t) is multiplied by the modified [0.1, 1] binary random mask matrix M(t) with a sampling period T and an input scaling factor γ to obtain the electrically injected masked signal S(t). It will be modulated onto the carrier wave of the light source: u(t)=[I′(t+Δt)I(t)],(8)S(t)=γ×u(t)×M(t).(9)

In the reservoir layer: The modulated carrier waveform from the light source is injected into a semiconductor laser with optical feedback to excite transient dynamics. The nonlinear dynamics of the reservoir laser is described by the Lang–Kobayashi equations [18,34]: $$\begin{gathered} Er(t)=1+iαr2[G(t)−1τp]Er(t)+σEs(t)exp(iΔwst) \\ +kfEr(t−τf)exp(−iwrτf), \end{gathered}$$(10)Nr(t)=pfJth−Nr(t)τe−G(t)|Er(t)|2,(11)G(t)=g[Nr(t)−N0]/1+ε|Er(t)|2),(12)where Er(t) denotes the complex, slowly varying electric field amplitude of the reservoir laser; Nr(t) represents the carrier density; αr is the linewidth enhancement factor; G(t) is the optical gain; N0 is the transparent carrier density; Es(t)=Is2{1+exp[iS(t)]} is the output of the light source, which is modulated by S(t) with a Mach–Zehnder intensity modulator; σ is the photonic injection strength; kf stands for the self-feedback strength of reservoir laser; Is is the optical intensity of the light source; and Δws=ws−wr=2πΔfs is the frequency detuning between the modulated light source and the reservoir laser, where ws and wr are the angular frequencies of the light source and reservoir laser, respectively. The values of other typical parameters are shown in Table 2.Table 2.

Values of Parameters for the Reservoir Lasers

Then NR neuron states of reservoir Xi(t) are extracted from the intensity of the reservoir laser, Ir(t)=|Er(t)|2, with time interval θ. In this work, the feedback delay time τf=(N+1)×θ. Additionally, the reservoir state X^i(t) of knowledge-based photonic reservoir is further extended by introducing the one-step prediction I′(t+Δt) of the imperfect dynamical model, as shown in Eq. (13). The former R1×Kv represents the connection to I′(t+Δt) from the imperfect knowledge-based model, while the latter RN×Kv connects the reservoir state matrix to the output. This approach not only integrates data information but also enhances the nonlinearity of the reservoir layer: X^i(t)=(I′(t+Δt)Xi(t)).(13)

In the output layer: The ridge regression algorithm is conducted to train and optimize the output connection weights vector W^out. The predicted results of knowledge-based photonic reservoir computing can be obtained by calculating y¯(t+Δt)=∑i=1NRW^outX^i(t). To evaluate the prediction accuracy quantitatively, the normalized mean-square error (NMSE) is calculated between the predicted y¯(t) and the “true” time series y(t) as follows: NMSE=1m∑t=1m[y¯(t)−y(t)]2σ[y(t)],(14)where m denotes the number of prediction samples. The value of the NMSE equals zero when the predicted time series closely matches the “true” time series. We set the threshold value for accurate prediction when the NMSE is less than 0.1.

For TDRC to achieve good performance, it must satisfy two key characteristics of coherence and short-term memory, which are mainly affected by the parameters of reservoir lasers. Moreover, the input gain coefficient γ can adjust the amplitude of the input signal that is modulated into reservoir lasers and regulate the prediction performance. Therefore, we investigate the performance of isolated TDRC in the parameter spaces of the input gain coefficient γ, self-feedback strength kf, and injection strength from the light source σ. First, Fig. 2(a) depicts the predictive capability NMSE versus the input gain coefficient γ, and the minimum NMSE = 0.0954 is achieved when the value of γ is around 1. Then we choose the best input gain coefficient γ=1 for the mask signal and alter the parameters of the reservoir laser. Figures 2(b) and 2(c) show the performance of chaotic time series forecasting with different self-feedback strength kf and injection strength from the light source σ. It can be seen that the smallest NMSE = 0.0683 is obtained with kf=8ns−1 and σ=15ns−1, where the dynamics of the reservoir laser is close to neutral stability [54].

In the knowledge-based model, the deviation of ηi is introduced to quantify the discrepancy between the constructed “imperfect” knowledge-based model and the “true” spin-flip model of a VCSEL that might be encountered in a real-life situation. ηi∈[−1,1] means that the value of parameter mismatch is in the range of −100%to100%. In this work, the discrepancy of the “imperfect” knowledge-based model stems from deviated self-feedback strength γf1=(1+η1)γf and dynamical state mismatch [Evx′(t),Evy′(t),Nv′(t),nv′(t)]=[(1+η2) ·Evx(t),(1+η2)Evy(t),(1+η2)Nv(t),(1+η2)nv(t)] in each iteration step. Figure 3 explores the influence of error coefficients η1 and η2 on the performance of heterogeneous forecasting of chaotic dynamics in a VCSEL. The value of γf is set to be 25ns−1 for η1=0. The dynamical state [Evx′(t),Evy′(t),Nv′(t),nv′(t)] is derived from the simulation of the “true” spin-flip model of the VCSEL with η2=0. We emphasized that the values of the NMSE in Figs. 3(a) and 3(b) remain unchanged as the error coefficients are not applied to the isolated TDRC. While it can be clearly seen from Figs. 3(c) and 3(d) that even though the smallest NMSE = 0.000252 is achieved at η1=0 and η2=0, our proposed heterogeneous forecasting scheme still behaves better compared to isolated TDRC with relatively large value of ηi, which means that the predictive performance of chaotic dynamics in the VCSEL could be enhanced significantly with the introduction of the “imperfect” physical model. For the scenario that η2=−1, as the one-step prediction I′(t+Δt) of the “imperfect” knowledge-based model tends to zero infinitely, the heterogeneous forecasting scheme has no effect on the prediction. At this point, the values of the NMSE for the heterogeneous TDRC and the isolated TDRC are both equal to 0.0656. We elucidate that the error coefficient of the proposed heterogeneous forecasting scheme is not limited to deviated feedback strength and dynamical state mismatch; it can also be applied to other parameters of the “true” physical models.

Then we use the above optimized parameters of isolated TDRC and induced error coefficients to evaluate the heterogeneous forecasting performance of chaotic dynamics generated by a VCSEL, and in order to eliminate the fluctuation caused by random masks, the values of NMSE are averaged over different sets of mask signals. The red line in Figs. 4(a) and 4(b) depict the chaotic time series of the VCSEL generated by numerical simulation of the spin-flip model, while the blue line plots the predicted VCSEL chaotic time series by isolated TDRC and heterogeneous TDRC, respectively. And the green line in Figs. 4(c) and 4(d) shows the difference between the predicted values and the original data. It can be observed that, in both cases, the predictive time series are similar to the original data. However, a smaller value of NMSE = 0.000634 is observed in the case of heterogeneous TDRC, which is 2 orders of magnitude lower.

We elucidate that the performance of heterogeneous forecasting for chaotic dynamics generated by a VCSEL is robust against different original signals as well as parameter changes in the spin-flip model. Figure 5 plots the values of NMSE as a function of self-feedback strength γf and current coefficient μ in the spin-flip model of a VCSEL as shown in Eqs. (1)–(4). Different chaotic dynamics of a VCSEL are generated with the self-feedback strength γf∈[10,30]ns−1 and current coefficient μ∈[2,3]. In the “imperfect” knowledge-based model of heterogeneous TDRC, the deviated self-feedback strength γf1 will be γf∈[20,60]ns−1 with η1=1. As shown in Figs. 5(a) and 5(b), the proposed heterogeneous forecasting scheme can successfully predict the chaotic time series of VCSELs over a relatively broader range of self-feedback strengths γf and current coefficients μ. Additionally, the values of NMSE for heterogeneous TDRC [Figs. 5(c) and 5(d)] are an order of magnitude lower than that of isolated TDRC, which indicates that our proposed scheme has better performance and stronger robustness.

In order to evaluate how the complexity of chaotic dynamics generated by a VCSEL affects the predictive performance, we introduce a rotating polarizer into the extended cavity of a VCSEL and variable polarization optical feedback is realized by rotating the polarizer with angle θp. As a result, the polarization properties and the complexity of chaotic dynamics generated by the VCSEL are manipulated. Figures 6(a) and 6(b) illustrate the Poincaré sphere to characterize the polarization state of chaotic dynamics generated by a VCSEL with optical feedback angles θp=10° and θp=50°. For θp=10°, the points on the Poincaré sphere are mainly concentrated near the positive X axis, indicating the dominance of the XP mode. While for θp=50°, the points are spread across the entire Poincaré sphere, and no dominant mode exists, indicating that the chaotic dynamics of the VCSEL are more complex. Figure 6(c) calculates the predictive performance NMSE of isolated TDRC and heterogeneous TDRC as a function of the polarizer angle θp. In both cases, the values of the NMSE achieve its maximum value around θp=50°, which corresponds to the most complex chaotic dynamics of a VCSEL. Additionally, Fig. 6(d) displays the zoomed-in trajectory of the NMSE for heterogeneous TDRC. It can be clearly seen that the performance of prediction is significantly improved with our proposed mechanism, and superior predictive performance could still be preserved with the most complex original signal.

The conventional form of isolated TDRC is a purely data-driven method and there is no knowledge about the physical model that generates the original data. However, in heterogeneous TDRC, an additional “imperfect” physical model of a VCSEL is introduced into photonic reservoir computing to mitigate the deficiency of the data-based approach, which yields improved processing speed, increased accuracy, and simplified parameter optimization. In order to demonstrate the pronounced benefits of our proposed scheme, the predictive performances of isolated TDRC and heterogeneous TDRC are investigated systematically in the parameter spaces of pf×Δfs and σ×kf, which are the operating parameters of reservoir lasers. The sampling interval of virtual nodes is fixed at θ=10ps, and the number of virtual nodes is NR=50. Thus, the sampling period of the mask signal is T=NR×θ=500ps, corresponding to a data processing speed of about 1/T=2GSa/s. Figure 7(a) first presents the NMSE evolution of isolated TDRC in the parameter space of pf×Δfs. To map the input signal to a higher-dimensional space and achieve better prediction performance, a balance between coherence and short-term memory is necessary. As shown in Fig. 7(a), the minimum value of NMSE=0.0340 for isolated TDRC is obtained with Δfs=−6GHz and pf=1.16. Predictive performance can be further optimized by fine-tuning the reservoir laser parameters. Figure 7(b) calculates the NMSE as a function of kf and σ. It is shown that the injection strength from the light source needs to be strong enough to realize injection locking, while the self-feedback strength kf needs to be finely tuned to manipulate the dynamics of reservoir lasers close to the edge of chaos. On the other hand, as displayed in Figs. 7(c) and 7(d), improved predictive performance of reservoir computing can be achieved in a wide range of parameter spaces of reservoir lasers, which can mitigate tedious parameter optimization effectively of the reservoir computing.

Moreover, we demonstrate that the heterogeneous forecasting scheme can realize an improved predictive performance with accelerated information processing speed and reduced training dataset. To quantify the advantages of the heterogeneous forecasting scheme quantitatively, the relationship of the prediction accuracy, information processing speed (characterized by the reciprocal of θ×NR), and size of the training dataset Kv is investigated in Fig. 8. The number of virtual nodes NR determines whether the reservoir can fully learn the features of the training data. To investigate the improvement in prediction performance under different processing speeds, Figs. 8(a1) and 8(a2) calculate the values of NMSE as a function of different numbers of virtual nodes NR and virtual node intervals θ. The number of virtual nodes NR ranges from 20 to 400, and the virtual node interval θ ranges from 4 to 30 ps, which corresponds to a data processing speed of about 0.0833–12.5 GSa/s. As shown in Figs. 8(b1) and 8(b2), compared with the NMSE between the isolated and heterogeneous TDRC, smaller NMSEs are obtained with reduced virtual nodes and shortened virtual node intervals, resulting in an accelerated information processing speed. Additionally, the size of the training data Kv is crucial for balancing prediction performance NMSE and training efficiency. Figure 8(a3) shows that the value of NMSE=0.00279 for heterogeneous TDRC at Kv=4000 is comparable to that of isolated TDRC at Kv=10,000, which improves the training efficiency effectively.

Figure 9 illustrates the all-fiber experimental setup for heterogeneous TDRC based on semiconductor lasers. A DFB (KG-BF-1550-F-S-FA) is adopted as the master laser. The output of the master laser is first amplified by an erbium-doped fiber amplifier (EDFA) and sent into the Mach–Zehnder intensity modulator after passing through a polarization controller. Then the signal is modulated by the masked signal S(t), which is generated by an arbitrary waveform generator (AWG, Keysight 8195A, 64 GSa/s). The reservoir laser is composed of a DFB with delayed optical feedback, and 50% of the reservoir laser’s output will be combined with the modulated signal from the master laser and fed back to the reservoir laser by a 50:50 coupler. The polarization controller is adjusted to realign the modulated optical signal with the polarization of the reservoir laser. Both the master and reservoir lasers are driven and temperature stabilized by the laser current and temperature controllers (Thorlabs, PROG8000). The neuron states of reservoir lasers are detected by a photodetector (PD, HP 11982A) and recorded by a high-speed digital oscilloscope (OSC, WaveMaster, 813Zi-B, 40 GSa/s). In the experiment, the virtual node interval for the reservoir is set to be θ=100ps, and the number of virtual nodes is set to NR=50. Consequently, the TDRC system’s clock period is T=θ×NR=5ns.

To demonstrate the effectiveness of the heterogeneous forecasting scheme for chaotic dynamics generated by a VCSEL with knowledge-based photonic reservoir computing, Fig. 10 experimentally investigates the predictive performance NMSE of isolated TDRC and heterogeneous TDRC. As shown in Figs. 10(a) and 10(b), the predictive results could mimic the original signals successfully in both cases. However, the predictive performance NMSE of heterogeneous TDRC is an order of magnitude lower than that of isolated TDRC. During the experiment, the optical delay of the reservoir laser is fixed at 16.7 ns, the temperature of the master and reservoir lasers is maintained at 22.7°C and 22°C, and the bias current is set to be 70 mA and 10 mA, respectively. The optical spectrum is measured by an optical spectrum analyzer (Yokogawa), and the center wavelengths of the master and reservoir lasers are 1549.006 nm and 1548.784 nm. Figures 10(c) and 10(d) plot the values of the NMSE as a function of self-feedback strength kf and optical injection strength σ. It can be seen that the experimental results agree well with the previous numerical calculations, and superior performance of the heterogeneous forecasting scheme is validated again in the whole range of parameter spaces, indicating that the heterogeneous forecasting mechanism based on the imperfect physical model helps improve the performance of traditional isolated TDRC.

In summary, we have proposed a heterogeneous scheme for forecasting the chaotic dynamics in VCSELs with knowledge-based photonic reservoir computing. An additional “imperfect” physical model is introduced to the traditional isolated TDRC which is solely based on data. The theoretical and experimental results indicate that the proposed heterogeneous forecasting scheme could achieve improved performance compared with both methods applied individually. We elucidate that the inclusion of an “imperfect” physical model could mitigate the tedious parameter optimization of traditional isolated TDRC and yield increased accuracy, accelerated information processing speed, and reduced training data size. Furthermore, it is demonstrated that even though the prediction accuracy will be manipulated by the complexity of chaotic dynamics, the performance of heterogeneous TDRC could still be preserved with different original signals. Our results suggest new insights into the prediction of chaotic dynamics of semiconductor lasers and potentially lead to solutions for optical secure communication.