Massively parallel and programmable photonic differential equation solver

Massively parallel and programmable photonic differential equation solver

PDF

Jiahao Wang¹^,², Wen Chen³, Zhou Zhou⁴, Dongyu Hu¹^,⁵, Zile Li¹^,²^,^*, Peng Chen³^,^*, Yan-qing Lu³, Shuang Zhang⁶, Cheng-Wei Qiu⁴, Shaohua Yu⁷, Guoxing Zheng¹^,²^,⁵^,^*

Opto-Electronic Advances | 2026, 9(5) : 250244

Less

Opto-Electronic Advances | 2026, 9(5): 250244

• Article •

Massively parallel and programmable photonic differential equation solver

Full

Jiahao Wang¹^,², Wen Chen³, Zhou Zhou⁴, Dongyu Hu¹^,⁵, Zile Li¹^,²^,^*, Peng Chen³^,^*, Yan-qing Lu³, Shuang Zhang⁶, Cheng-Wei Qiu⁴, Shaohua Yu⁷, Guoxing Zheng¹^,²^,⁵^,^*

Affiliations

¹Electronic Information School, Wuhan University, Wuhan 430072, China

²Peng Cheng Laboratory, Shenzhen 518055, China

³National Laboratory of Solid State Microstructures, and College of Engineering and Applied Sciences, Nanjing University, Nanjing 210093, China

⁴NUS Graduate School and Department of Electrical and Computer Engineering, National University of Singapore, Singapore 119077, Singapore

⁵Wuhan Institute of Quantum Technology, Wuhan 430206, China

⁶Department of Physics, University of Hong Kong, Hong Kong 999077, China

⁷Chinese Academy of Engineering, Beijing 100088, China

Published: 2026-05-15 doi: 10.29026/oea.2026.250244

Outline

Abstract

Less

Calculus equations are fundamental mathematical tools, whose numerical solution is crucial. Existing solvers with optical analog computing struggle to simultaneously integrate programmability and parallel processing, thus constraining computational speed and density. Herein, we propose a reconfigurable all-optical platform capable of solving variable-coefficient first-order ordinary differential equations in parallel. We utilize the electrically tunable liquid crystals (LCs) as computing kernels to address these equations. The solver's applicability to canonical scientific problems, such as heat conduction and resistor-capacitor circuit dynamics, is further showcased with simultaneous solving of 158 equations with only one single forward propagation of light. Experimental results confirm the efficacy of the platform in solving equations in an ultra-fast, reconfigurable, broadband, and parallel manner.

Key words

photonics computing / ordinary differential equations / photonic differential equation solver

Cite this Article

Jiahao Wang, Wen Chen, Zhou Zhou, Dongyu Hu, Zile Li, Peng Chen, Yan-qing Lu, Shuang Zhang, Cheng-Wei Qiu, Shaohua Yu, Guoxing Zheng. Massively parallel and programmable photonic differential equation solver[J]. Opto-Electronic Advances, 2026 , 9 (5) : 250244 - . DOI: 10.29026/oea.2026.250244

Full Text

Less

1　Introduction

Less

Differential equations^1,2 form a cornerstone of mathematics, underpinning a wide range of applications in physics^3,4, biomedicine⁵, engineering^6,7, and many other disciplines⁸⁻¹⁰. They are pivotal in modeling diverse phenomena, such as thermoelectric effects¹¹, compressive sensing¹², optical neural networks¹³, and nonlinear optics¹⁴. While solving differential equations is a central focus across many research areas, analytical solutions are feasible only for selected few simple cases. In most real-world applications, numerical methods that leverage computational techniques to obtain approximate solutions are essential¹⁵.

In high-throughput data processing applications, traditional data processing methods based on integrated circuits¹⁶ encounter significant challenges, including substantial transmission losses, slow processing speeds, and excessive heat generation. These issues create bottlenecks at the integration level, ultimately limiting the efficiency of data processing. All-optical computing¹⁷⁻²¹ has emerged as a transformative technology, offering high-speed, parallel, and low-energy processing of large-scale data. Its application to solving calculus equations in the spatial domain has been widely investigated. A variety of material systems and photonic computing frameworks have been proposed to achieve high-speed solutions for calculus equations, including microring resonators²²⁻²⁵, metasurfaces²⁶⁻³⁰, metamaterials³¹⁻³⁴, and others³⁵⁻³⁸. These developments have brought a wide range of applications in intelligent agents³⁹ such as wearable devices⁴⁰, pattern recognition⁴¹, and edge-enhanced depth perception⁴². A wide variety of optical design methodologies have consequently emerged⁴³⁻⁴⁶. However, existing optical architectures for solving calculus equations in the spatial domain^26,29−32 are constrained by their inability to modulate the computing kernel of the constructed equations. This limitation restricts their application to calculus equations with fixed coefficients.

Real-world problems often involve variations in parameters, necessitating adjustments to the coefficients of the calculus equations. This creates an urgent demand for solvers capable of handling variable-coefficient calculus equations. While utilizing microring resonators or pixel metamaterials can solve time-domain variable-coefficient ordinary differential equations (ODEs), they lack the inherent parallel processing capability of spatial domain signal processing. Furthermore, the need for input signal multiplexing processing further restricts the speed of equation solving. As an important branch of calculus equations, first-order ODEs constitute a fundamental component. Many methods for solving higher-order ODEs rely on transforming them into first-order ODEs. Therefore, developing solutions for variable-coefficient first-order ODEs is of great importance.

In this study, we propose and experimentally demonstrate a spatially all-optical variable-coefficient first-order ODE solver (SAO-VODE Solver) based on electrically tunable liquid crystals (LCs). With their high transparency and electrically tunable properties, LCs have become one of the most attractive reconfigurable material platforms for optical computing. By carefully arranging the orientation distribution of the LCs' director, we align the LC transfer function to the Taylor series expansion of the first-order ODE transfer function. By precisely adjusting the voltages applied to the LCs, we can modulate the coefficients of the first-order ODEs and solve multiple equations simultaneously across multiple wavelengths. We further showcase the solver’s applicability to canonical scientific problems, such as heat conduction under varying temperatures and resistor-capacitor circuit dynamics with varying input voltages. Our method combines the parallel computing capability of spatial-domain processing with the reconfigurability of time-domain processing, delivering an efficient all-optical computing solution with high speed, full parallelism, and multitasking capability.

2　Principles

Less

2.1　SAO-VODE Solver concept

The concept of SAO-VODE Solver is illustrated in Fig. 1(a). Two Fourier lenses are employed to construct a classic 4f system, and the tunable LC solver is placed at the confocal plane of the 4f system as the computing kernel for ODEs. By adjusting the voltage applied to the LC solver to impart different phases to the input signals, it implements different transfer functions corresponding to ODEs with varied coefficients. Leveraging the inherent parallel processing capability of spatial domain optical signal processing, multiple input signals g(x) with varying waveforms can be processed simultaneously. The solutions f(x) are obtained independently and simultaneously on the image plane of the 4f system. Our approach combines both spatial signal processing capability, multitasking ability, and reconfigurability, as illustrated in Fig. 1(b).

Figure 1(c) illustrates the simulation results of the solutions of trigonometric function g₁(x) = sin(120πx) and quadratic polynomial function g₂(x) = 500x² + x + 0.01 under different voltages U_n (n=1, 2, 3, 4). When a constant voltage is applied to the LCs, their dispersive properties facilitate the solution of equations with variable coefficients for incident light of different wavelengths. Furthermore, by utilizing the electrically tunable characteristics of LCs, we can achieve solutions to differential equations with variable coefficients without the need to change the wavelength. These two properties of LCs significantly enhance the capability for high-throughput signal processing.

2.2　Taylor expansion and Fourier transform methods

We employed the Fourier transform method to solve the first-order ODEs with variable coefficients. Without loss of generality, we consider one-dimensional scenarios in which input signals and the orientation distribution of the LCs’ director depend solely on the horizontal coordinate x. As shown in Fig. 2(a), a general first-order ordinary differential equation can be expressed as follows

$ A{f}\,{'}\left(x\right)+Bf\left(x\right)=g\left(x\right) \;, $

(1)

where A and B are the coefficients, g(x) is a known signal and f(x) is the function to be solved. By employing the Fourier transform method and first-order Taylor expansion when k_x approaches 0 m⁻¹, the transfer function of Eq. (1) can be expressed as follows:

$ \widehat{H}\left({k}_{x}\right)=\frac{F\left({k}_{x}\right)}{G\left({k}_{x}\right)}=\frac{1}{{\mathrm{i}}A{k}_{x}+B}\approx \frac{1}{B}\left(-{\mathrm{i}}\frac{A}{B}{k}_{x}+1\right) \;, $

(2)

where F(k_x) and G(k_x) are the Fourier transforms of f(x) and g(x), respectively. k_x denotes the spatial frequency along the x-axis direction.

An electrically tunable planar liquid-crystal can be used to implement the transfer functions in Eq. (2). This implementation remains effective even when the parameters A and B vary. The relevant transmission coefficient (see Supplementary Section 1 for its derivation), which governs this behavior, is given by

$ H\left(\theta,\varphi\right)=\mathit{\mathit{\mathrm{cos}}}\frac{\varphi}{2}-{\mathrm{i}}{\mathrm{sin}}\frac{\varphi}{2}{\mathrm{sin}}2\theta\; , $

(3)

where θ represents the in-plane orientation angle of the LCs director. By placing the LCs at the confocal plane of the 4f system, and setting θ = [arcsink_xλf/(L/2)]/2 and k_x = x/(λf), where L denotes the size of the LCs, Eq. (3) can be reformulated as

$ H\left(k_x,\varphi\right)={\mathrm{cos}}\frac{\varphi}{2}\left(-{\mathrm{i}}{\mathrm{tan}}\frac{\varphi}{2}.\frac{\lambda f}{L}k_x+1\right)\; . $

(4)

By comparing Eq. (2) and Eq. (4), we can see that Eq. (4) exhibits an expression identical to Eq. (2) upon the removal of an overall constant coefficient. Furthermore, the phase delay φ associated with the LCs is correlated with the coefficients A and B in the ODEs, by the relation as follows

$ \varphi=2\mathrm{arctan}\frac{AL}{B\lambda f}\; . $

(5)

According to Eq. (5), we can construct the transfer functions of the first-order ODEs with variable coefficients by adjusting the working wavelength, or applying different voltages on the LCs to attach different phase delays φ. This forms the basic principle of the ODE solving via liquid-crystal photonic solver.

Figure 2(b) shows the designed orientation distributions of θ based on the above analysis, where the azimuth angle of the LCs varies along the x-coordinate. The inset illustrates the 3D direction of the LC molecules, including the in-plane orientation θ and the out-of-plane angle α. The right panel of Fig. 2(b) illustrates the relationship between the coefficient A/B, the phase delay φ, and the out-of-plane angle α. To demonstrate the effectiveness of our proposed method for solving the first-order ODEs, we present simulation results obtained under different ODE coefficients while maintaining a fixed input signal g(x). As illustrated in Fig. 2(c), the input signal expression is g(x) = cos(40x) + 0.8sin(30x + π/4) + 0.5sin(20x + π/7) + 3. Figure 2(d–f) illustrates the theoretical transfer function H(k_x), alongside the transfer function Ĥ(k_x), and the first-order Taylor expansion approximation of H(k_x), implemented by the LCs. We can observe that as the spatial frequency k_x approaches 0 m⁻¹, the designed transfer function Ĥ_m(k_x) = −iAk_x/B² + 1/B (m=1, 2, 3) exhibits high consistency with the theoretical transfer function H_m(k_x) = 1/(iAk_x + B) (m=1, 2, 3). As the results, Fig. 2(g–i) presents the output signals f(x) after processing through both the theoretical transfer function and the transfer function derived through Taylor expansion. It is noted that the simulation results in Fig. 2(g–i) have been appropriately translated and scaled for easy comparison. It is observed that the theoretical results closely align with the simulation results while discrepancies between the two transfer functions lead to a slight deviation between the theoretical and simulation outcomes. As the coefficient A/B increases, the range of k_x where H(k_x) and Ĥ(k_x) match well becomes progressively narrow. To ensure the LC solver processes only the aligned components, the spatial frequency range of the kernel can be limited by increasing the focal length f. In our experiments, constrained by the focal length of the lens, we conducted demonstrations using a relatively small coefficient A/B. Nonetheless, the proposed method is fully scalable to scenarios involving large coefficient A/B.

3　Results and discussions

Less

3.1　Characterization of SAO-VODE Solver

The experimental setup is schematically illustrated in Fig. 3(a). A 50× magnification objective lens and lens 1 are arranged in a confocal configuration to collimate and expand the incident beam. An aperture is employed to control the beam size. A signal mask based on a 1951 USAF resolution test chart (see Supplementary Section 3, for the photograph of the signal mask), is used to generate various input signals g(x) by selectively illuminating different regions of the resolution chart, as illustrated in Fig. 3(e). The 1951 USAF resolution test chart is positioned at the front focal plane of the Fourier lens 2, while the solver is placed at the confocal plane of the 4f system (f = 150 mm). A linear polarizer and an analyzer are positioned before and after the LC device, respectively, to control the polarization states of the input and transmitted light. The output signals are captured by a charge-coupled-device (CCD) camera located at the back focal plane of the Fourier lens 3.

Figure 3(b) shows the microscopic polarized image of the designed LC solver (see Methods for the fabrication process of the LC solver). The working area (2.12 mm × 2.12 mm) is marked within an orange dashed box. Figure 3(c) illustrates the polarization conversion efficiency (PCE) and phase curves of LCs at the wavelengths of 590 nm and 630 nm. The solid lines correspond to the PCE as a function of voltage U (see Supplementary Section 1 for more details about PCE measurement and phase curves calculation). Figure 3(d) shows the variable coefficient A/B varies with respect to voltage U. The zoomed-in photography framed by the orange box in Fig. 3(e) highlights the actual two-dimensional signal region, with each column corresponding to different input signals.

To validate the broadband working characteristics of the proposed method, two distinct operating wavelengths, 590 nm and 630 nm were selected. Using Eq. (5), the phase values φ corresponding to different coefficients A/B can be calculated. By comparing these phase values with the curve shown in Fig. 3(d), the corresponding voltage values can be determined (see Supplementary Section 4 for more details about the relationship between the voltage U and different coefficients A and B at various wavelengths λ). The proposed SAO-VODE Solver is able to process input signals directly, without requiring the analog-to-digital signal conversion. The input signals are obtained directly from the spatial domain capture. In the first column of Fig. 3(f), the input signals are generated by illuminating different areas in the 1951 USAF resolution chart with a laser operating at 590 nm and 630 nm, corresponding to the intensity distribution marked at various locations on the signal mask in Fig. 3(e).

Firstly, to demonstrate the capability of the proposed SAO-VODE Solver in processing multiple input signals, we obtained ODE solving results for three different input signals g_m(x) (m=1,2,3), as illustrated in the second column of Fig. 3(f). The third column of Fig. 3(f) demonstrates that by adjusting the voltage U, the coefficient A/B can be modified, enabling the corresponding solution results to be obtained. The fourth column illustrates that expected solutions for the ODEs can also be achieved at different wavelengths. To avoid variations in focal length for incident light of different wavelengths caused by lens dispersion. In this work, we employ an achromatic lens (LBTEK, MAD410-A) designed for the visible spectrum, the focal length tolerance of the lens is ±1%. During the experiment, the signal mask, Fourier lens, liquid crystal solver and CCD (charge-coupled device) are all mounted on a five-axis translation stage (PDV, SZ500), which provides an adjustment precision of 5 μm in each of the x, y, and z directions. The slight focal length variation introduced by the incident wavelength can be compensated by adjusting the five-axis translation stage to maintain the conjugate planes of the 4f system. All input and output signals have been normalized for visualization. To quantitative evaluation of the deviation between the experimental results f_d(x) and theoretical results f_t(x), we define the solving error as E = (∫ (f_t(x) − f_d(x))²dx) / (∫(f_t(x))²dx × ∫(f_d(x))²dx). According to quantitative evaluation results, the deviation between the experimental and theoretical results is mostly within 0.05. The results show that the experimental results agree well with the theoretical predictions across multiple wavelengths and parallel input signals, discontinuous signals are directly processed in the spatial domain, better reflecting the input signals distribution typically encountered in practical scenarios, though minor discrepancies exist due to the transfer function approximation. Additionally, positional offsets between the zero frequency of the input signals and the center position of the LC solver, as well as the limited design focal length, can introduce errors in the results (see Supplementary Sections 5, 6, and 7 for more details about the error analysis). It is noted that the experimental results in Fig. 3(f) have been appropriately translated and scaled for easy comparison (see Supplementary Section 8, for the details of the selection of the scaling and translation coefficients). To the best of our knowledge, this represents the first demonstration of a method to solve first-order variable-coefficient ODEs in the spatial domain.

3.2　Optical solutions to heat conduction and resistor-capacitor circuit problems

Solving differential equations has significant applications in thermodynamics, electromagnetics, and other scientific domains. Here, we present two practical examples demonstrating the application of the proposed SAO-VODE Solver for solving differential equations. The first example addresses the transient heat transfer problem of small-sized spheres in a gradually varying low-temperature environment, as illustrated in Fig. 4(a). According to Newton’s law of cooling, the temperature T_o(t) of the high-temperature object and the environment temperature T_e(t) are governed by the following ODE:

$ \frac{1}{k}\frac{\mathrm{d}T_{\mathrm{o}}\left(t\right)}{\mathrm{d}t}+T_{\mathrm{o}}\left(t\right)=T_{\mathrm{e}}\left(t\right)\; , $

(6)

where k denotes the thermal conductivity. Another example is the resistor-capacitor (RC) circuit, as illustrated in Fig. 4(b). For a simple RC series circuit with a time-varying voltage V(t), the charge q(t) between the capacitor satisfies the following ODE:

$ R\frac{\mathrm{d}q\left(t\right)}{\mathrm{d}t}+\frac{1}{C}q\left(t\right)=V\left(t\right)\; , $

(7)

here, R denotes the resistance and C denotes the capacitance. To generate various input signals T_e(t) and V(t), a custom-designed signal board is employed, as indicated by the white dashed region in Fig. 4(c). The signal board consists of 768 × 768 pixels, each with dimensions of 41 μm × 41 μm. Different columns of the signal board encode signals with variable coefficients. In the experiment implementation, temporal variations of the input signals are mapped onto spatial variations along the x-axis. The left half of the signal board represents a quadratic polynomial, expressed as T_e(x) = 500x² + x + (100 + 2m) × 10⁻⁴, where m gradually increases from 0 to 384 toward the center of the signal board. The right half corresponds to a sinusoidal wave, expressed as V(x) = sin [(60 + n) πx] + 1, where n gradually decreases from 384 to 0 when moving away from the center. Two representative regions from the signal board are selected as input signals for the heat conduction scenario and RC circuit scenario, respectively, as indicated by input function sets (IFS) 1 and 2 in Fig. 4(c). The intensity distributions of the input signals along the x-axis direction of the input signals are shown for four distinct lines, denoted as g_m (m=1, 2, 3, 4).

Figure 4(d) illustrates the corresponding output function sets (OFS) obtained experimentally for the heat conduction scenario with 1/k = 7.5 × 10⁻³ and the RC circuit with R = 1.5 × 10⁴ and C = 0.5 × 10⁻⁶. Similarly, Fig. 4(e) illustrates the experimental results for cases with 1/k = 2.5 × 10⁻³ and the RC circuit with R = 5 × 10³, while keeping C as a constant. The resulting output signals, denoted as f_m (m=1, 2, 3, 4), correspond to the solutions of input signals g_m (m=1, 2, 3, 4) under the respective scenarios.

Additionally, parallel solution results are illustrated. In Fig. 4(c), the input signal window size is 6.5 mm × 7 mm, corresponding to a simultaneous solution of 158 equations (6.5 mm/41 μm = 158). Among these signals, four representative signals are selected to facilitate a direct comparison between the experimental and theoretical results. The excellent agreement observed between them underscores the reliability and accuracy of the proposed approach. These two practical examples validate the capability and applicability of the SAO-VODE Solver in solving real-world differential equations, thereby demonstrating its potential for broader application scenarios.

At last, as shown in Table 1, a comparison is presented between this work and other existing methods for solving calculus equations. Our method solves differential equations with a single forward propagation of light waves, eliminating the need for multiple iterations of light signals to achieve stability, as required by methods based on metastructures (MS), metagratings (MG), and pixel metamaterials (PMM). Moreover, with a single operation, we experimentally achieve the parallel solving of 158 equations, which represents an improvement of 10² magnitude in the number of equations solved per computation compared to existing methods for solving calculus equations. Existing spatial-domain computing methods such as acoustic metasurface (AMS), MS and MG based on inverse-designed approach, are difficult to be integrated with reconfigurable or programmable devices.

4　Conclusions

Less

We have proposed and experimentally demonstrated an electrically tunable liquid crystal-based optical solver capable of simultaneously solving multiple first-order ODEs with variable coefficients in the spatial domain. The approach begins with a Taylor expansion-based analysis of the transfer function corresponding to first-order ODEs, which is matched to the complex amplitude modulation characteristics of the liquid crystal. By electrically tuning the applied voltage to attach different phases, the coefficients of the differential equations can be flexibly adjusted. The system’s capability is validated experimentally through a series of tests involving both synthesized and physically meaningful input signals. Furthermore, we demonstrate its application to practical scenarios in thermal conduction and electrical circuits. The experimental results exhibit excellent agreement with theoretical predictions, thereby confirming the accuracy and robustness of the proposed method. Notably, parallel computation of up to 158 differential equations is achieved within a single measurement, and this capacity can be further scaled by increasing the sensor area.

The parallel processing characteristics inherent to spatial domain signals, combined with the high-speed nature of optical solutions, offer substantial potential for high-throughput data processing applications. Additionally, the system supports signal processing across multiple wavelengths, further enhancing its versatility. It is anticipated that inputting time-varying signals will enable even greater processing throughput, extending the applicability of the proposed approach.

Compared to metasurface- and metamaterial-based approaches, the LC solver stands out as a lithography-free solution with greater flexibility and lower fabrication complexity, but achieves equation-solving capability as good as those lithography-made devices. Metasurface-based devices exhibit static computing capabilities after fabrication, whereas the LC solver can achieve reconfigurable and adaptable functionalities, benefiting from its soft-matter building blocks. Overall, the proposed approach provides a reconfigurable, ultra-fast, and highly compact solution for multi-channel broadband parallel signal processing, which holds promising application potential in fields such as all-optical diffractive neural networks, feature recognition, and signal processing.

5　Materials and methods

Less

5.1　Fabrication of the SAO-VODE Solver

The fabrication of the SAO-VODE Solver primarily involves two steps: the preparation of the LC cell and ultraviolet (UV) photopatterning (Supplementary Fig. S2). First, the ITO glass substrates undergo ultrasonic cleaning and UV-Ozone treatment. Next, a 0.35% solution of the sulphonic azo-dye SD1 in dimethylformamide is spin-coated onto the substrates. After curing at 100 °C for 10 minutes, the 3.23 μm-thick cell is then sealed using spacer-doped epoxy glue between two SD1-coated substrates. To transfer the desired azimuth pattern to the SD1 layer, we use a digital-micromirror-device (DMD)-based UV microlithography system, which includes a light source, components for dynamic pattern generation and focusing, and a monitor. Specifically, a UV beam carrying the designed pattern is reflected onto the DMD (Discovery 3000, Texas Instruments), focused by a tunable lens, polarized by a motorized polarizer, and projected onto the LC cell. The empty cell coated with SD1 was exposed to linearly polarized UV light through a multistep partly overlapping exposure process with synchronized polarizer rotation, reorienting the SD1 molecules perpendicular to the UV polarization. Upon injecting the E7 LCs, the patterned SD1 effectively guides the orientation of the LC directors through intermolecular interactions, resulting in a fully assembled LC cell with the desired azimuth orientations.

5.2　PCE measurement

To calibrate the PCE of the LC solver at three specific wavelengths, we utilize a super-continuum light source (YSL SC-pro) and a commercial optical power meter (Thorlabs, PM100D), as shown in Supplementary Fig. S1. An aperture is placed in front of the LC solver to minimize the influence of stray light. In addition, a left-hand circular polarizer (LCP polarizer) and a right-handed circular analyzer (RCP analyzer) are utilized to eliminate the unwanted co-polarized light. By applying square wave signals with a frequency of 1 kHz and varying peak-to-peak amplitude generated by a signal generator (Tektronix, AFG31052) to the LC solver, we obtained a series of intensities of selected wavelengths at different applied voltages. Consequently, the PCE of the SAO-VODE Solver can be measured at different wavelengths.

References

Less

Arutyunov AV, Zhukovskaya ZT, Zhukovskiy SE. On nonlinear boundary value problems for differential inclusions. Differ Equations 59, 1443–1450 (2023).

Guo B, Phong D H. On L^∞ estimates for fully non-linear partial differential equations. Ann Math 200, 365–398 (2024).

Nurisso M, Morandini M, Lucas M et al. Higher-order Laplacian renormalization. Nat Phys 21, 661–668 (2025).

Kira M, Koch SW, Smith RP et al. Quantum spectroscopy with Schrödinger-cat states. Nat Phys 7, 799–804 (2011).

Ouyang XY, Sutradhar S, Trottier O et al. Neurons exploit stochastic growth to rapidly and economically build dense dendritic arbors. Nat Commun 16, 5903 (2025).

Paolino A, Nava G, Di Natale FD et al. Learning aerodynamics for the control of flying humanoid robots. Commun Eng 4, 111 (2025).

He WW, Li JZ, Kong X et al. Multi-level physics informed deep learning for solving partial differential equations in computational structural mechanics. Commun Eng 3, 151 (2024).

Zhou Z, Zhang YH, Xie YX et al. Electrically tunable planar liquid-crystal singlets for simultaneous spectrometry and imaging. Light Sci Appl 13, 242 (2024).

Liang X, Zhu D, Dai Q et al. All-optical multi-order multiplexing differentiation based on dynamic liquid crystals. Laser Photonics Rev 18, 2400032 (2024).

Yang HW, Xie WC, Chen HF et al. Spin–orbit optical broadband achromatic spatial differentiation imaging. Optica 11, 1008–1016 (2024).

Takahagi A, Hirai T, Alasli A et al. Observation of the transverse Thomson effect. Nat Phys 21, 1283–1289 (2025).

Giorgetta FR, Potvin S, Deschênes JD et al. Free-form dual-comb spectroscopy for compressive sensing and imaging. Nat Photonics 18, 1312–1319 (2024).

Xue ZW, Zhou TK, Xu ZH et al. Fully forward mode training for optical neural networks. Nature 632, 280–286 (2024).

Huang CY, Luo Y, Zhao YL et al. A conformal mapping approach to broadband nonlinear optics on chip. Nat Photonics 18, 471–477 (2024).

Soetaert K, Petzoldt T, Setzer RW. Solving differential equations in R: package deSolve. J Stat Softw 33, 1–25 (2010).

Li WS, Gong XS, Yu ZH et al. Approaching the quantum limit in two-dimensional semiconductor contacts. Nature 613, 274–279 (2023).

Deng M, Cotrufo M, Wang J et al. Broadband angular spectrum differentiation using dielectric metasurfaces. Nat Commun 15, 2237 (2024).

Ding XM, Zhao ZH, Xie P et al. Metasurface-based optical logic operators driven by diffractive neural networks. Adv Mater 36, 2308993 (2024).

Huo PC, Tan L, Jin YQ et al. Broadband and parallel multiple-order optical spatial differentiation enabled by Bessel vortex modulated metalens. Nat Commun 15, 9045 (2024).

Wang XW, Hao HJ, He XY et al. Advances in information processing and biological imaging using flat optics. Nat Rev Electr Eng 1, 391–411 (2024).

Wang XW, Wang H, Wang JL et al. Single-shot isotropic differential interference contrast microscopy. Nat Commun 14, 2063 (2023).

Abro KA, Das B. A scientific report of non-singular techniques on microring resonators: an application to optical technology. Optik 224, 165696 (2020).

Wu JY, Cao P, Hu XF et al. Compact tunable silicon photonic differential-equation solver for general linear time-invariant systems. Opt Express 22, 26254–26264 (2014).

Yang T, Dong JJ, Lu LJ et al. All-optical differential equation solver with constant-coefficient tunable based on a single microring resonator. Sci Rep 4, 5581 (2014).

Yuan HY, Du ZC, Qi HX et al. Microcomb-driven photonic chip for solving partial differential equations. Adv Photonics 7, 016007 (2025).

Abdollahramezani S, Chizari A, Dorche AE et al. Dielectric metasurfaces solve differential and integro-differential equations. Opt Lett 42, 1197–1200 (2017).

Chizari A, Abdollahramezani S, Jamali MV et al. Analog optical computing based on a dielectric meta-reflect array. Opt Lett 41, 3451–3454 (2016).

Cordaro A, Edwards B, Nikkhah V et al. Solving integral equations in free space with inverse-designed ultrathin optical metagratings. Nat Nanotechnol 18, 365–372 (2023).

Zhang WX, Qu C, Zhang XD. Solving constant-coefficient differential equations with dielectric metamaterials. J Opt 18, 075102 (2016).

Cotrufo M, Sulejman SB, Wesemann L et al. Reconfigurable image processing metasurfaces with phase-change materials. Nat Commun 15, 4483 (2024).

Silva A, Monticone F, Castaldi G et al. Performing mathematical operations with metamaterials. Science 343, 160–163 (2014).

Estakhri NM, Edwards B, Engheta N. Inverse-designed metastructures that solve equations. Science 363, 1333–1338 (2019).

Fu PY, Xu ZM, Zhou TK et al. Reconfigurable metamaterial processing units that solve arbitrary linear calculus equations. Nat Commun 15, 6258 (2024).

Zangeneh-Nejad F, Sounas DL, Alù A et al. Analogue computing with metamaterials. Nat Rev Mater 6, 207–225 (2020).

Zangeneh-Nejad F, Fleury R. Performing mathematical operations using high-index acoustic metamaterials. New J Phys 20, 073001 (2018).

Zangeneh-Nejad F, Fleury R. Topological analog signal processing. Nat Commun 10, 2058 (2019).

Zuo SY, Wei Q, Tian Y et al. Acoustic analog computing system based on labyrinthine metasurfaces. Sci Rep 8, 10103 (2018).

Tang YH, Chen RY, Lou MH et al. Optical neural engine for solving scientific partial differential equations. Nat Commun 16, 4603 (2025).

Luo XG. Digital optics and optical intelligent agent. Appl Phys Lett 127, 010501 (2025).

Zhao BH, Cheng JW, Wu B et al. Integrated photonic convolution acceleration core for wearable devices. Opto-Electron Sci 2, 230017 (2023).

Han YN, Xiang SY, Song ZW et al. Pattern recognition in multi-synaptic photonic spiking neural networks based on a DFB-SA chip. Opto-Electron Sci 2, 230021 (2023).

Liu XY, Zhang JC, Leng BR et al. Edge enhanced depth perception with binocular meta-lens. Opto-Electron Sci 3, 230033 (2024).

Badloe T, Lee S, Rho J. Computation at the speed of light: metamaterials for all-optical calculations and neural networks. Adv Photonics 4, 064002 (2022).

Lee S, Park C, Rho J. Mapping information and light: trends of AI-enabled metaphotonics. Curr Opin Solid State Mater Sci 29, 101144 (2024).

So S, Mun J, Park J et al. Revisiting the design strategies for metasurfaces: fundamental physics, optimization, and beyond. Adv Mater 35, 2206399 (2023).

Yang Y, Kang H, Jung C et al. Revisiting optical material platforms for efficient linear and nonlinear dielectric metasurfaces in the ultraviolet, visible, and infrared. ACS Photonics 10, 307–321 (2023).

Cavicchioli G, Melloni A, Miller D A B et al. Programmable photonic architecture solving systems of ordinary differential equations. In International Conference on Photonics in Switching and Computing 1–3 (IEEE, 2023). http://doi.org/10.1109/PSC57974.2023.10297288.

Appendix

Less

Year 2026 volume 9 Issue 5

PDF

101

Cite this Article

BibTeX

Article Info

doi: 10.29026/oea.2026.250244

Receive Date：2025-09-15
Online Date：2026-07-02
Published：2026-05-15

Article Data

Affiliations

History

Received：2025-09-15
Accepted：2026-01-20

Affiliations

¹Electronic Information School, Wuhan University, Wuhan 430072, China

²Peng Cheng Laboratory, Shenzhen 518055, China

³National Laboratory of Solid State Microstructures, and College of Engineering and Applied Sciences, Nanjing University, Nanjing 210093, China

⁴NUS Graduate School and Department of Electrical and Computer Engineering, National University of Singapore, Singapore 119077, Singapore

⁵Wuhan Institute of Quantum Technology, Wuhan 430206, China

⁶Department of Physics, University of Hong Kong, Hong Kong 999077, China

⁷Chinese Academy of Engineering, Beijing 100088, China

Corresponding:

ZL Li, E-mail: lizile@whu.edu.cn

P Chen, E-mail: chenpeng@nju.edu.cn

GX Zheng, E-mail: gxzheng@whu.edu.cn

References

Share

https://castjournals.cast.org.cn/joweb/oea/EN/10.29026/oea.2026.250244

Share to

Scan QR to access full text

Cite this article

BibTeX

Citations

表12种不同金属材料的力学参数

科 Family	属数 Number of genus	种数 Number of species	占总种数比例 Percentage of total species (%)	属 Genus	种数 Number of species	占总种数比例 Percentage of total species (%)
鹅膏菌科Amanitaceae	2	11	5.26	鹅膏菌属 Amanita	10	4.78
小菇科 Mycenaceae	2	12	5.74	丝盖伞属 Inocybe	5	2.39
多孔菌科 Polyporaceae	8	14	6.70	蜡蘑属 Laccaria	5	2.39
红菇科 Russulaceae	3	23	11.00	小皮伞属 Marasmius	6	2.87
				小菇属 Mycena	11	5.26
				光柄菇属 Pluteus	5	2.39
				红菇属 Russula	17	8.13
				栓菌属 Trametes	5	2.39

关闭全屏

BibTeX
EndNote
RefWorks
TxT

Ref.	Methods	PT (cycles)	N_eqs (single operation)	Reconfigurable	Operation domain
ref.²³	MRR	NG	1	√	TD
ref.³⁷	AMS	NG	1	×	SD
ref.²⁹	DMM	NG	1	×	SD
ref.³²	MS	50	1	×	SD
ref.²⁸	MG	60	1	×	SD
ref.³³	PMM	30	1	√	TD
ref.⁴⁷	MZI	1	3	√	TD
This work	LCs	1	158	√	SD

图1 Schematic representation of the spatially all-optical variable-coefficient first-order ordinary differential equation solver. (a) The architecture of SAO-VODE Solver, featuring a 4f system, a liquid crystal ODE solver to construct the computing kernel, multiple input signals g_m(x) (m=1, 2, 3) and corresponding solution output signals f_m(x) (m=1, 2, 3) under varied voltages (for different coefficients A/B). (b) The characteristic classification of existing calculus equation processors. This work combines both spatial signal processing capability, multitasking ability, and reconfigurability. (c) Illustration of the solutions of trigonometric function sin(120πx) and quadratic polynomial function 500x² + x + 0.01 under different voltages U_n (n=1, 2, 3, 4). The abbreviation 'Norm. Amp.' has the full name of 'Normalized Amplitude'.

图2 The optical computing flowchart and the simulation results of differential equations solutions with variable coefficients. (a) The optical computing flowchart. The input and output signals are denoted as g(x) and f(x), respectively. The spatial Fourier transform and spatial inverse Fourier transform are implemented by using two Fourier lenses. (b) The two-dimensional orientation distributions of the LC solver and the relationship between the coefficient A/B, the phase φ, and the out-of-plane angle α. The orientation distribution in the blue dashed region is further illustrated with the orange ellipse indicating the orientation distributions along the x-axis direction. The diagram in the corner illustrates the direction of the LC molecules, characterized by the in-plane orientation ${\theta} $ and the out-of-plane angle α. The right panel illustrates the relationship between the coefficient A/B, the phase φ, and the out-of-plane angle α. (c) The illustration of the input signal g(x). The abbreviation 'Norm. Amp.' has the full name of 'Normalized Amplitude'. (d-f) Illustrations of the theoretical transfer function H_m(k_x) = 1/(iAk_x + B) (m=1, 2, 3) and the designed transfer function Ĥ_m(k_x) = −iAk_x/B² + 1/B (m=1, 2, 3) in the complex plane, with varied coefficients A and B. Re.: real part, Im.: imaginary part. (g-i) Results of the theoretical output signal f_m(x)-theor. (m=1,2,3) and the simulated output signal f_m(x)-sim. (m=1,2,3) for variable coefficients A and B. In the simulation, A/B is set as 0 (d, g), 1 (e, h) and 10 (f, i) respectively. The abbreviation 'theor.' has the full name of 'theoretical', the abbreviation 'sim.' has the full name of 'simulated'.

图3 Experimental setup and measurement of the LCs properties. (a) Experimental setup for ODE solving. The signal mask is the 1951 USAF resolution chart. (b) Photograph of the fabricated LC solver sample under orthogonally polarized light. The scale bar is 1 mm. (c) LC polarization conversion efficiency (PCE) and phase delay φ curves versus applied voltages at the wavelengths of 590 nm and 630 nm. (d) The coefficients curves of A/B of ODEs versus applied voltages on LC at the wavelengths of 590 nm (blue) and 630 nm (red). The image inside the red and blue box further showcase the details of the curves. (e) The illustration of the input signals g_m(x) (m=1,2,3). The image inside the orange box represents the input signal, and the yellow dashed lines indicate the input signal used in (f). (f) Comparison of the experimental results and theoretical results of ODE solution f_m(x) (m=1,2,3) under different wavelengths λ, voltages U, coefficients A/B, and input signals g_m(x) (m=1,2,3). The abbreviation ‘Norm. Amp.’ has the full name of ‘Normalized Amplitude’. E represents the deviations between the experimental results and theoretical results.

图4 Diagrams of practical examples for solving differential equations. (a) Schematic diagram of the heat conduction, a small high-temperature sphere with temperature T_o(t) and thermal conductivity k, is placed in a low-temperature environment with varying temperature T_e(t). (b) A schematic diagram of the RC series circuit, a varying voltage V(t) is applied. (c) The bright-field microscopic images of the signal board under parallel polarization states. The signal board consists of sinusoidal signals and second-order polynomial signals, with signal intensity varying along the x-axis direction and signal coefficient varying along the y-axis direction. IFS 1 and 2 present the experimental input signals of the second-order polynomial type and sinusoidal type, respectively. The right side shows the signal intensity distributions along the x-direction corresponding to the white dashed lines g_m (m=1,2,3,4). The red line represents the design signal, and the purple line represents the actual input signal. (d, e) The comparison between experimental and theoretical results for parallel solving under different equation coefficients. f_m (m=1,2,3,4) corresponds to the solution results of input signals g_m (m=1,2,3,4). The experimental results of OFS 1 and 2 correspond to the solution results of IFS 1 and 2. Both the experimental and theoretical results are subjected to simple sliding window smoothing for easy comparison. Both the input and output signals are normalized.

Articles: Latest Articles; Most Read; Collections

Updates: Events; News; Multimedia

About: About Us

Contact

No. 86 Xueyuan South Road, Haidian District, Beijing

100081

010-62199257

qkjq@cast.org.cn

Copyright © 2025 China Association for Science and Technology. All rights reserved. For all open access content, the relevant licensing terms apply.
Sponsored by the Office of the Leading Group for Cybersecurity and Informatization of CAST, and supported by Science and Technology Review Publishing House