Ordinary Differential Equation Solver Using Microwave Photonics with Integrated Microcombs

1. Introduction

Mathematical models play a critical role in the analysis and design of systems. Among them, ordinary differential equations (ODEs) have widespread applications in describing dynamic phenomena across various domains, such as control systems, mechanical and electrical systems, and meteorology [1,2,3,4,5]. There is a growing need for faster and more energy-efficient solutions to ODEs, especially for ultrafast processes like rocket launches or data-intensive applications such as automated manufacturing, weather prediction, and epidemic modeling. As the speed and power efficiency of electronic analog circuits approach their bottleneck in the post-Moore era [6], photonic technologies have emerged as a compelling alternative for realizing ultrafast and broadband analog computing systems [7].

Recently, various all-optical ODE solvers, which are realized by designing infinite impulse response (IIR) filters with suitable transfer functions, have been implemented based on resonators and optical feedback loops to perform recursive operations [2,8,9,10,11]. Semiconductor optical amplifiers and optical filters have been used as differentiators and wavelength converters to achieve first-order ODE solving, where the feedback configuration required precise phase control and suffered from instability [8]. Second-order ODE solvers have subsequently been demonstrated by employing cascaded microring resonators (MRRs) [9]. In addition, systems of ODEs have been solved through using an add/drop MRR and introducing two external feedback waveguides to realize the coupling between these ODEs [2]. However, their coefficients are inherently coupled, determined by the coupling and feedback parameters, which limits the reconfigurability of the ODE solvers. In practical applications, to characterize general linear time-invariant (LTI) systems, it is necessary to solve ODEs with various orders and different coefficients. To overcome this limitation, tunable interferometric couplers have been introduced to enable tunable ODE coefficients, but the realization of arbitrary coefficients remained challenging due to the interdependence between coupling strength and phase of interferometric couplers [10,11].

With the rapid development of microwave photonics (MWP) technologies, photonic analogue computing has been extended to the microwave domain. In contrast to conventional microwave signal processing based on electronics, which face intrinsic bandwidth constraints, the use of photonic hardware and technologies to process high-bandwidth microwave signals, or MWP processing [12,13,14], can provide orders of magnitude faster speeds, which is critical for high-speed processing applications [15,16]. To achieve reconfigurable MWP ODE solvers based on the aforementioned IIR filter structure (Figure 1(a)), cascading multiple MRRs, precise alignment, and feedback loops are required, which makes the system highly sensitive to fabrication tolerance and thermal drift. An MWP ODE solver has been demonstrated using an add/drop MRR to solve a first-order ODE, where the ODE coefficient was determined by the internal decay loss of the MRR [17]. Although this approach eliminates the need for an external feedback loop, the tunability of coefficients and scalability are still restricted.

Transversal filter structures based on finite impulse response (FIR) filters can be employed to remove recursive operation and enhance reconfigurability [18,19,20,21]. The transversal filter structure performs convolution between the input microwave signal and a set of tap weights (Figure 1(b)). By appropriately designing the tap weights, which define the system’s impulse response, the convolution module can directly emulate the desired transfer function corresponding to a given ODE. This approach enables independent and reconfigurable coefficient design while maintaining stability and broadband performance.

Here, we employed a transversal filter structure, which enables high reconfigurability by using an integrated optical microcomb [22], to achieve MWP ODE solvers. Optical microcombs can provide a large number of wavelength channels by using compact micro-scale resonators [23,24], which is critical for improving the processing accuracy of MWP transversal signal processors. They are also with the ability to offer broad Nyquist zones, which allow for large processing bandwidths [25,26,27]. Optical microcombs enable MWP transversal systems to have small size, weight, and power consumption (SWaP) [23,24,28,29], and excellent compatibility with monolithic integration, making them powerful alternatives to conventional multiwavelength sources, for example, discrete laser arrays [30,31,32].

We propose and experimentally demonstrate order- and coefficient-tunable MWP ODE solvers based on a transversal filter structure using a soliton crystal microcomb source, where first- and second-order ODEs with different coefficients are solved. System demonstrations of ODE solvers with an operation bandwidth of ~24.5 GHz for Gaussian pulse input signals are performed. For the simplified first-order, general first-order, and second-order ODE solvers, with root-mean-square errors (RMSEs) between the measured and theoretical curves being ~0.030, ~0.068, and ~0.091, respectively. This good alignment confirms the effectiveness of our approach for realizing high-speed reconfigurable MWP ODE solvers.

2. Operation Principle

LTI systems are commonly demonstrated by a constant-coefficient linear ODE, which is defined as [10]

$${\displaystyle \frac{d^{n}y\left(t\right)}{d{t}^n}}+{\displaystyle {\sum }_{i=0}^{n-1}}{a}_i{\displaystyle \frac{d^{i}y\left(t\right)}{d{t}^i}}={\displaystyle {\sum }_{k=0}^m}{b}_k{\displaystyle \frac{d^{k}x\left(t\right)}{d{t}^k}}$$

(1)

where x(t) is the input signal of an LTI system, and y(t) is the output signal representing the solution of the ODE. a_i and b_k are the constant coefficients. Based on Eq. (1), a simplified first-order ODE, a general first-order ODE, and a second-order ODE can be expressed as [10,11]

$${\displaystyle \frac{dy\left(t\right)}{dt}}+a_{0}y\left(t\right)=x(t)$$

(2)

$${\displaystyle \frac{dy\left(t\right)}{dt}}+a_{0}y\left(t\right)=b_1{\displaystyle \frac{dx\left(t\right)}{dt}}+b_{0}x(t)$$

(3)

$${\displaystyle \frac{d^{2}y\left(t\right)}{d{t}^2}}+a_1{\displaystyle \frac{dy\left(t\right)}{dt}}+a_{0}y\left(t\right)=b_2{\displaystyle \frac{d^{2}x\left(t\right)}{d{t}^2}}+b_1{\displaystyle \frac{dx\left(t\right)}{dt}}+b_{0}x(t)$$

(4)

After Fourier transformation from Eqs. (2) – (4), transfer functions of the simplified first-order ODE H_fs (ω), the general first-order ODE H_fg (ω), and the second-order ODE H_s(ω) can be given by [33]

$$H_{fs}\left(\omega \right)={\displaystyle \frac{1}{j\omega +a_0}}$$

(5)

$$H_{fg}\left(\omega \right)={\displaystyle \frac{b_{1}j\omega +b_0}{j\omega +a_0}}$$

(6)

$$H_s\left(\omega \right)={\displaystyle \frac{{b_2\left(j\omega \right)}^2+b_{1}j\omega +b_0}{{\left(j\omega \right)}^2+a_{1}j\omega +a_0}}$$

(7)

where j = $\sqrt{-1}$, and ω is the angular frequency.

To implement the ODE solvers, we use a transversal filter structure as shown in Figure 1(b) that offers high reconfigurability in terms of its spectral transfer function, which can be expressed as [25,34]

$$H(\omega )=\underset{k=0}{\text{∑}^{M-1}}{h}_k{e}^{-j\omega n\Delta T},$$

(8)

where ω is the angular frequency of the input microwave signal to be processed, M is the tap number, h_k (k = 0, 1, 2, …, M−1) is the tap weight of the k^th tap, and ΔT is the time delay between adjacent taps. By properly designing the various tap weights h_k (k = 0, 1, 2, …, M-1), different signal processing functions can be realized by using a single system without changing the hardware.

An optical microcomb is used to serve as the multi-wavelength source, which generates multiple wavelength channels that act as discrete taps for the transversal signal processor. The generated optical microcomb is spectrally shaped according to the designed tap weights h_k (k = 0, 1, 2, …, M−1). Next, all of the wavelength channels of the shaped optical microcomb are imprinted with the input microwave signal via an electro-optic modulator (EOM), leading to the generation of multiple microwave replicas. Following this, the modulated optical signals are transmitted through a dispersive medium, for example, a single-mode fibre (SMF), to introduce time delays ∆T between adjacent wavelength channels, which progressively separate the microwave replicas. Finally, the delayed replicas are summed upon photodetection via a photodetector. ODE solvers with different orders and coefficients can be realized by properly setting the tap weights h_k (k = 0, 1, 2, …, M−1). Tap weights are calculated through inverse Fourier transformation (IFT) from Eqs. (5) – (7).

The corresponding radiofrequency (RF) amplitude response of the simplified first-order, general first-order, and second-order ODE solvers with different coefficients as a function of tap numbers (M) is shown in Figure 2(a) − (c), respectively. For simplified first-order ODEs, the coefficients are (1) a₀ = 0.2×10¹⁰, (2) a₀ = 2×10¹⁰, and (3) a₀ = 20×10¹⁰. For general first-order ODEs, the coefficients are (1) a₀ = 0.5×10¹⁰, b₀ = 1×10¹⁰, b₁ = 1×10¹⁰, (2) a₀ = 2×10¹⁰, b₀ = −8.5×10⁹, b₁ = 1, and (3) a₀ = 10×10¹⁰, b₀ = 1×10¹⁰, b₁ = 1. For second-order ODEs, the coefficients are (1) a₀ = 5×10²⁰, a₁ = 5×10¹⁰, b₀ = 1.5×10²¹, b₁ = 1.5×10¹⁰, b₂ = 1, (2) a₀ = 1×10²², a₁ = 1×10⁸, b₀ = −1×10²⁰, b₁ = −1×10²⁰, b₂ = 1, and (3) a₀ = 10×10³⁰, a₁ = 5×10¹⁰, b₀ = 4×10¹⁴, b₁ = 5×10⁹, b₂ = 1. ODEs with different orders and coefficients can be solved without requiring optical feedback or cascaded MRRs to perform recursive operations, verifying the reconfigurability of the proposed system. The simulated results obtained with different tap numbers (M = 5 − 45) are shown together with the ideal transfer functions. As can be seen, the discrepancies between the ideal and simulated amplitude RF response are improved with tap numbers for all the three ODE solvers. For a small M, the limited spectral sampling introduces errors, leading to ripples and distortions in the amplitude response.

We use the root-mean-square error (RMSE) to quantitatively analyze these discrepancies, which is defined as [24]

$$\mathrm{RMSE}=\sqrt{\underset{i=1}{\text{∑}^n}{\displaystyle \frac{(Y_i\text{–}{y}_i)2}{n}}}$$

(9)

where n is the number of sampled points, Y₁, Y₂, …, Y_n are the values of the ideal results, and y₁, y₂, …, y_n are the values of the simulated results. Figure 3(a) − (c) shows the RMSEs between the ideal and simulated RF response of the simplified first-order, general first-order, and second-order ODE solvers presented in Figure 2 as a function of tap number M. As expected, for all cases, the RMSEs are inversely proportional to tap numbers, confirming that the processing accuracy is influenced by the number of taps. As M increases, the discrepancies are gradually suppressed, and RMSEs reach a minimal value (typically < 0.01), indicating agreement between the simulated RF response of the proposed system and the ideal transfer function. This trend is consistent across different coefficient sets, verifying the robustness of the system. We note that when tap numbers increase, RMSEs initially decrease dramatically and reduce more gradually as tap numbers become larger, which indicates that the need for large tap numbers is not strong. For the second-order ODE solver (2), both b₀ and b₁ are negative, which shifts the system zeros toward the imaginary axis [35], enhancing the high-frequency gain of the transfer function as shown in Figure 2(c-ii), so the RMSE shows small improvement as the tap number increases.

3. Results

We used the experimental setup shown in Figure 4 to implement MWP ODE solvers, which included a comb generation (Figure 4(a)) and a signal processing module (Figure 4(c)). In our experiments, we employed a soliton crystal microcomb as the multi-wavelength source [19]. Soliton crystal microcombs represent a unique class of optical frequency combs, characterized by self-organized ensembles of multiple co-propagating solitons that form a crystal-like structure in the angular domain [22,36]. As shown in Figure 4a, the MRR used for generating soliton crystal microcombs was fabricated on a complementary metal-oxide-semiconductor (CMOS)-compatible doped silica glass platform [37,38], with a quality (Q) factor of ~1.9 million and a radius of ~592 μm, yielding a free spectral range (FSR) of ~49 GHz (i.e., ~0.4 nm).

The initiation of the soliton crystal microcomb was achieved via a simple pump wavelength sweeping method [39], where a continuous-wave (CW) pump laser (Yenista Optics), amplified to ~32.1 dBm by an erbium-doped fibre amplifier (EDFA, IdealPhotonics), was swept in wavelength from blue to red until modulation instability oscillations emerged. As the detuning between pump wavelength and MRR’s resonance was progressively adjusted, a stable soliton crystal oscillation state was achieved, generating over 90 wavelength channels across the C band at a pump wavelength of ~1551.3 nm, as shown in Figure 4(b). As can be seen, the soliton crystal microcomb exhibits a palm-shaped spectrum which, although often considered a potential limitation for practical use, has in fact been shown not to present a significant drawback. This is partly attributed to the higher power conversion efficiency of soliton crystal microcombs compared with single soliton states [22]. In addition, the soliton crystal microcomb demonstrated remarkable long-term power stability, with a measured relative standard deviation of −14 dB over 66 hours [40], highlighting its stability for use in practical MWP systems.

Figure 4(c) shows the schematic of the signal processing module. The first Waveshaper (Finisar) was employed to flatten the spectrum of the initially generated soliton crystal microcomb to reduce the power difference between different comb lines. The flattened microcomb was then amplified and fed into an EOM (iXblue), where the optical signal was modulated by an input microwave signal. We used a Gaussian pulse with a pulse width of ~0.1 ns as the input signal, which was produced by an arbitrary waveform generator (AWG, Keysight). The input signal was broadcast to each wavelength channel, generating multiple replicas of the microwave signal. Then, the replicas propagated in an SMF with a length of L = ~2.153 km, a dispersion of D₂ = ~17.4 ps/nm/km, which introduced a time delay of ΔT = L × D₂ × Δλ = ~15.0 ps, corresponding to an operation bandwidth of ~33.4 GHz. We note that the maximum operational bandwidth of the proposed system is limited by the repetition rate of the soliton crystal. Large crosstalk between adjacent taps occurs for microwave operation beyond 24.5 GHz, which is half of the soliton crystal’s repetition rate. Next, according to the designed tap weights, the comb lines were shaped by using the second Waveshaper. Finally, the delayed and weighted replicas were combined through photodetection via the balanced photodetector BPD (Finisar).

A two-stage feedback control was implemented to mitigate static and slow-varying errors [41]. The first stage reduced intensity deviations introduced during the comb flattening process. The shaped optical spectrum at the output of the first Waveshaper was measured using an optical spectrum analyzer (OSA, Anritsu) and transmitted to a computer, where it was compared with the desired tap weights. The second stage was to further compensate for errors caused by imperfections of the impulse response. It was performed channel by channel with the same input microwave signal modulated onto each comb line. The tap weights (i.e., the peak intensities of the impulse response) were measured from the BPD recorded by the oscilloscope (OSC, Keysight). These measured results were then subtracted from the designed tap weights to generate error signals, which were used to calibrate the attenuation of comb line intensities in the Waveshaper. All components of the system were included in the two feedback loops, enabling compensation for errors resulting from static and slowly varying error sources introduced by different components.

According to the results in Figure 3, we selected specific tap numbers (i.e., number of wavelength channels) for different ODE solvers. We used 15, 11, and 15 comb lines as discrete taps for the three simplified first-order ODE solvers. We employed 25, 17, and 11 comb lines for the three general first-order ODE solvers (1) – (3), respectively. For second-order ODE solvers, we used 17 comb lines for (1) and (2), and 15 comb lines for (3).

Figure 5 shows the temporal processing results of aforementioned three simplified first-order ODE solvers (1) – (3). Figure 5(a) shows the measured tap weights for the three ODE solvers with different coefficients compared with the desired tap weights. The measured comb line magnitude (blue cross) exhibits agreement with the desired weights (red circle), which confirms the high spectral shaping accuracy resulting from the two-stage feedback control. Figure 5(b) demonstrates the corresponding temporal output of the ODE solvers. The measured outputs (yellow) show good agreement with the theoretical outputs (red), and the waveform of the input Gaussian signal (blue) is also shown. The calculated RMSEs between the theoretical and measured results are ~0.030, ~0.026, and ~0.026, respectively, which indicates the high accuracy of the ODE solvers. Minor temporal broadening and slight amplitude differences are mainly attributed to the fast-varying noise from the microcomb and BPD, which cannot be compensated for by the two-stage feedback control [29].

Figure 6 shows the corresponding processing output of three general first-order ODE solvers (1) – (3) mentioned in Figure 2(b). Figure 6(a) shows the tap weights corresponding to three different coefficient sets, where the measured tap weights (blue cross) closely match the ideal targets (red circle). The corresponding temporal results are presented in Figure 6(b), where the measured outputs (yellow) exhibit good agreement with the theoretical outputs (red). The calculated RMSEs between the theoretical and measured outputs for (1) – (3) are ~0.020, ~0.054, and ~0.068, respectively. The results confirm that the proposed system can accurately realize general first-order ODE solvers.

We also note that, in our system, both positive and negative ODE coefficients have been realized by simply designing specific tap weights. Implementing negative ODE coefficients using an IIR-based photonic solving system is inherently challenging, where the sign and magnitude of the ode coefficient are determined by the amplitude and phase of the feedback paths. Hence, introducing a π-phase shift in the optical feedback signal is necessary to realize negative coefficients, which demands precise control of optical phase and temperature, as slight drifts can lead to phase mismatching and instability.

The corresponding processing results of the second-order ODE solvers in Figure 2(c) are presented in Figure 7. In Figure 7(a), the measured tap weights are shown and the ideal tap weights are also presented for comparison. The measured tap weights aligned well with the ideal tap weights. Figure 7(b) shows the waveform of the input signal and the theoretical as well as measured processing output. The RMSEs between the theoretical and measured output for different ODE coefficients (1) – (3) are ~0.027, ~0.065, and ~0.091, respectively. The measured processing results are consistent with the theoretical output, which indicates good agreement between experimental results and theory.

4. Discussion

Our results of MWP ODE solvers based on an integrated soliton crystal microcomb highlight a step toward fully reconfigurable and broadband analog computing systems. The transversal filter structure intrinsically avoids the stability and tunability limitations of ODE solvers based on IIR structure. The convolutional operation between the tunable tap weights and the signal to be processed makes our system a recursion-free platform with straightforward programmability.

In addition, the use of a two-stage feedback control method enables accurate calibration of tap coefficients, leading to high processing accuracy with RMSEs < 0.1. After calibration, there are still residual errors induced by components in the system that cannot be compensated by using the two-stage feedback control method. We infer that these errors are mainly induced by the noise floor of microcombs as well as fast-varying fluctuations from the microcomb and BPD. Self-calibrating photonic integrated circuits have been demonstrated [42,43], where the impulse response of the system was characterized by an optical reference path, enabling the establishment of a Kramers-Kronig relationship. Amplitude and phase distortions can be retrieved through Fourier-domain processing. This approach opens new opportunities for implementing a more precise feedback control in the MWP ODE solvers.

The operation bandwidth of the ODE solvers is jointly determined by the time delay between adjacent wavelength channels and half of the soliton crystal’s repetition rate. Large operation bandwidth can be achieved through employing a shorter spool of SMF and/or a microcomb source with a larger comb spacing, although at the expense of providing fewer taps.

Recent advances have driven rapid progress toward fully integrated microcomb-based MWP transversal filter systems [17,44]. Employing integrated microcomb sources to replace bulky multi-wavelength sources has already provided many advantages in terms of SWaP, cost, and complexity, further gains can be achieved by advancing the integration level of the entire system. In principle, all the components in the microcomb-based MWP transversal filter system can be monolithically integrated on a single chip, and on-chip CW lasers [45], optical amplifiers [46], EO modulators [47], dispersive elements [48], optical spectral shapers [10,49], and PDs [50] have all been demonstrated. Based on these integrated components, some sophisticated subsystems have also been demonstrated, including microcomb generation module consisting of heterogeneously integrated pump lasers and microresonators [51] and spectral shaping arrays [52,53].

Complex-valued tap weights have been achieved in MWP transversal filter systems recently [54], which provides a solution to extend the proposed ODE solver to support complex-valued taps, broadening its applications in modeling wave propagation, quantum-inspired dynamics, and coupled field systems. In addition, time-division multiplexing and matrix partitioning were used to load operators on a reconfigurable MRR array and thereby solve partial differential equations (PDEs) with > 90% accuracy on chip [55]. This describes a broader scope in which photonic hardware is no longer limited to first- or second-order ODE solvers. While the proposed microcomb-based transversal filter system primarily solves ODEs, it has the potential for solving PDEs, as linear PDEs can be converted into coupled or independent ODEs that can, in principle, be implemented using parallel ODE solver architectures. This work has broad implications for microcombs [56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86] and their applications to microwave photonics, neuromorphic processors and communications. [87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145] The addition and use of 2D materials [146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181] will add extra functionality to microcomb chips for potential applications to quantum photonics [182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198] and other areas. [199,200,201,202,203,204,205,206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224].

Conclusions

In summary, we propose and experimentally demonstrate order- and coefficient-tunable MWP ODE solvers based on an integrated microcomb source, capable of solving constant-coefficient linear first- and second-order ODEs with different coefficients. System demonstrations of the ODE solver with an operation bandwidth of ~24.5 GHz are performed, with RMSEs of ~0.030, ~0.068, and ~0.091 for the simplified first-order, general first-order, and second-order ODE solvers, respectively. The good agreement between measured and theoretical results confirms the effectiveness of our approach in realizing high-speed MWP ODE solvers with high reconfigurability. Our approach presents an effective way to implement versatile MWP ODE solver, offering a promising solution for applications in control systems, mechanical and electrical systems, and meteorology.