Partial differential equations (PDEs), which describe how a function changes with respect to multiple variables at the same time, are used in science and engineering to model complex systems, such as electromagnetics and fluid dynamics.1 It is a common goal for researchers in academic and industrial fields to solve PDEs in parallel with high accuracy.2 In practical applications, such as numerical weather forecasting,3^,4 petroleum exploration,5 aircraft design,6 and nuclear fusion simulation,7 valuable supercomputers with a great number of cores and specially designed architectures are required to solve PDEs,8^–10 especially time-evolving PDEs. However, the development of computers based on electronic components faces the physical limit and reaches a bottleneck in the post-Moore era,11^–14 which requires new computing schemes.15^,16

The photonic chip, which employs photons as information carriers, possesses fast speed and high parallelism. Based on photonic chips, photonic computing has made great progress, stimulated by the development of big data.17^–19 Photonic computing manifests dramatic advantages in delocalized edge computing,20^,21 advanced machine vision,22^–24 convolution accelerators,25^–28 photonic neural networks,29^–32 and mathematical operations.33^–36 Benefiting from the development of photonic chips, schemes for ordinary differential equations, which are easier to solve compared with PDEs, were first proposed based on microring resonators.37^–41 Meanwhile, the Fredholm integral equations can be solved with inversely designed metamaterial42^,43 and metasurfaces,44 respectively. As for the optical solving of PDEs, a kind of special photonic circuit for solving Poisson equations has been reported.45^,46 A photonic iterative processor for matrix inversion is presented to solve the Poisson equation. However, unlike static equations, such as the Poisson equation, the time-evolving PDE problems have the variable of time and are more complex to solve. The error accumulates as the solution evolves over time, so a high mesh accuracy is needed to ensure the correctness of evolution. The coefficient matrix size increases quadratically when the mesh number increases. Using a photonic chip with a limited size to compute a large PDE coefficient matrix is still difficult because a large matrix size usually leads to the rapid expansion of device numbers, which brings challenges to integration and control. There is a lack of effective methods to compute PDE’s big coefficient matrix on a photonic chip. In this context, it is still challenging to solve time-evolving PDEs on a photonic chip.

Here, we introduce time-division multiplexing and matrix partition techniques in the photonic solving of time-evolving PDEs, which can solve PDEs with large coefficient matrices on a photonic chip with only a limited number of computing units. A reconfigurable photonic chip is designed and fabricated with standard complementary metal oxide semiconductor (CMOS) technology. The photonic chip, with a size of 3.7mm×2.5mm, consists of four inversely designed 1:3 power splitters and a 9×9 tunable microring array. An integrated nonlinear Kerr microcomb is employed to provide dense signal channels with a high signal-to-noise ratio and low cross talk. With our scheme, various kinds of time-evolving PDE problems are solved experimentally. The heat equation with the first order of time derivative is solved with an accuracy of above 95%, the wave equation with the second order of time derivative is solved with an accuracy of around 93.7%, and the nonlinear Burgers equation with the first order of time derivative is solved with an accuracy of above 94%. In addition, an example of static equations, that is, Poisson’s equation, is solved with a coefficient matrix size of 144×144 (mesh size of 12×12) and an accuracy of 93.5%. Furthermore, parallel solving of Poisson’s equation and Laplace’s equation is demonstrated based on a single photonic chip, with an accuracy of 95.9% and 95.8%, respectively. This work provides a powerful platform to solve different types of PDEs and lays a foundation for future on-chip photonic computers.

To solve a PDE numerically, the solving domain is discretized with the finite difference method, and the PDE is converted into a matrix form.47 The matrix equation is then solved based on the photonic chip with the time-division multiplexing and matrix partition methods. In the solving process, most of the computation tasks arise from matrix-vector multiplication (MVM). When the fineness of discretization increases, the computation cost of MVM grows in the form of a quadratic function. Using the finite difference method, the continuous solving domain is discretized into a 3n×3n mesh, where each grid point is an independent variable to be solved. These variables can be specified with a sequence number ranging from 1 to 9n2 in row-major order, which forms a variable vector with a size of 9n2. A coefficient matrix with a size of 9n2×9n2 is used to describe the relationship of the 9n2 variables. For example, if we name the coefficient matrix as R, the element Rmn is the interaction of variable um with un, which is determined by the PDE operator and the specific discretization method. By employing light to calculate the multiplication between the variable vector and the coefficient matrix, most computation tasks in the solving process can be completed during the propagation of light, offering an efficient solution to the challenges posed by the large-scale computations required for high-resolution discretization.33

A sketch map of a photonic computing system is shown in Fig. 1(a), where red lines denote the light path and gold lines denote the electric path. The system uses an optical frequency comb as the light source, and the generated optical signals are separated into nine distinct wavelength channels by a dense wavelength demultiplexer (DWDM). A set of variable optical attenuators (VOAs) is used to modulate the intensity of each wavelength channel, representing the values of the vector elements. These modulated multiwavelength signals are then recombined by a wavelength multiplexer and coupled into the photonic chip through grating couplers for computation.

The chip consists of an inverse-designed 1:9 power splitter and a 9×9 microring array. The power splitter evenly divides the input optical signals into the microring array. The microring array adjusts the resonant wavelengths of the microrings by tuning the electrode voltages, enabling the loading of the input matrix. To accommodate enough available wavelength channels, the radius of each microring is designed as 5.5μm, providing a free spectral range of ∼16nm. MVM is performed with the microring array, and the computation results are read out by optical power meters through the output gratings and stored in the cache. In addition, an external controller coordinates the modulation of optical signals, manages the solving process, and processes the data generated by photonic computing. Figure 1(b) shows a microscopic optical image of the photonic chip mounted on a printed circuit board (PCB), with each PCB containing two photonic chips. The electrode voltages are connected to the chip via gold wires to precisely adjust the optical properties of the microrings. The microscope image above in Fig. 1(c) shows horizontal silicon waveguides with a width of 450 nm and vertical metal wires. As shown on the left side of the top panel, four 1:3 power splitters are cascaded to allocate optical signal evenly to nine rows of microrings. An optimization algorithm is employed to design the structure, which is capable of considering both footprint and performance.48^–53 The optimization target of the 1:3 power splitter is to send a signal of a wavelength band to three ports as evenly as possible. The design process and simulated performance of the inversely designed power splitter can be found in Note 2 in the Supplementary Material. The detailed structure of the microrings is shown below. The horizontal distance between the centers of two microrings is 160μm. In the right bottom image, white color is used to denote silicon waveguides. Each microring is covered by a TiN electrode denoted with a dark green color. The yellow color denotes the layer of AlCu, which is used to conduct current.

The coefficient matrix of a PDE grows quadratically as the mesh accuracy increases. To handle MVM for such large-scale matrices, we propose to introduce time-division multiplexing and matrix partition techniques. As shown in Fig. 2(a), the large coefficient matrix is divided into smaller patches of size 3×3, with each patch representing the local interaction between adjacent grid points in the discretized solving domain. These 3×3 matrix patches are then loaded onto different parts of the microring array, referred to as MVM cores, where each core independently computes one matrix patch and operates in parallel with other cores. To further enhance computing efficiency and flexibility, time-division multiplexing is employed, where the input vector is divided into small slices and loaded across different time slots. The vector slices are loaded onto adjacent wavelength channels through VOAs and a DWDM. In each time slot, signals of the unwanted wavelength channels are attenuated to zero, and the external controller collects the sum through power meters of the needed rows. It is worth pointing out that the matrix partition technique will not greatly slow down the computing speed. On the one hand, if the matrix is not partitioned, the chip scale must be expanded to fit the big matrix, which will increase the time of light propagation. On the other hand, the vector to be multiplied with the unpartitioned matrix will also be large, which will increase the time of the interaction between the photonic chip and external controller each time. Thus, the matrix partition technique will not bring a significant slowdown to the solution.

The Kerr microcomb possesses the ability to generate multiple equidistant wavelength signals, and its unique narrow peaks and low cross talk characteristics provide significant advantages in optical communication and signal processing fields.54^–58 These features enable the optical frequency comb to achieve efficient spectral control while significantly reducing interference in multichannel signal processing, thus improving signal purity and system stability. Such advantages have driven the rapid development of optical computing architectures, with applications spanning high-speed data transmission and optical neural networks.25^,27^,59^,60 Here, a Kerr microcomb with comb spacing of around 0.8 nm at the telecommunication C-band is used as a signal source; the spectra in different wavelength ranges are shown in Note 5 in the Supplementary Material. Using DWDM, nine comb lines are selected, with power individually controlled by VOAs to load vector data. The measured output signals of the nine signal channels are shown with a blue line in Fig. 2(b). The distance between the two channels is 0.8 nm. Each signal can be tuned independently with a VOA to load the vector element. The red line in Fig. 2(b) shows the measured spectrum of the test ring of the photonic chip. The free spectrum range is from 1539.35 to 1555.13 nm, with a width of 15.78 nm, which covers the selected wavelength channels.

The parameters of the rings in the microring array of the photonic chip are the same. By tuning the voltages added to each ring, the resonant wavelength of a ring can be matched with one of the signal channels. By finely tuning the deviation between the resonant wavelength of the ring and the center wavelength of the matched signal channel, the transmission of the drop port can be set from 0 to 1. In this way, the coefficient matrix can be loaded onto the photonic chip. The relationship between the resonant wavelength of the ring and the added voltage determines the way to tune the photonic chip. The measured movement of the resonant wavelength when the different voltages are added is shown in Fig. 2(c).

As an example of the matrix partition method, the loading of the Laplacian operator is demonstrated here, which is used in solving PDE problems below. When the grid mesh size is set as 3m×3m and m is an integer larger than 1, the coefficient matrix consists of four nonzero patches. When m is more than 2, the matrix patches remain the same. By tuning the voltages on the electrodes of the microring array, the patches can be loaded. The top left part of the coefficient matrix of the Laplacian operator with a mesh size of 6×6 is shown in Fig. 2(d), where the duplicate negative elements are removed. Matrix patches are marked with red squares. The measured normalized spectra of nine rows of the microring array with nine DWDM wavelength channels opened, are shown in Fig. 2(e). The central wavelengths of the nine channels are marked with arrows on the top of the figure. The spectra are measured using a telecom optical spectrum analyzer (model: Yokogawa AQ6370D), with zero-attenuating in all wavelength channels. By switching the wavelength of the input light and measuring the output power with power meters, the power matrix can be obtained. The switching of the input wavelength is achieved by setting the attenuating value of one channel as zero and the one of the unwanted channels as maximum. The measured normalized power matrix of the microring array is shown in Fig. 2(f). Although there are differences in light power between different channels and different rows, what matters is to make the relative values in each 3×3 matrix patch as close to the target matrix as possible. Considering practical experimental conditions, the matrix patches of the coefficient matrix are loaded successfully.

Time-evolving PDEs are important in describing the evolution process of the state. The calculation of the solution at a time point is based on the solutions at the previous time points, which calls for frequent MVM operation. The accumulated computation error in the iteration process can be fatal to the final solution, and the noise in analog optical computing is usually unavoidable, which makes the optical solving of the time-evolving PDEs elusive. Here, we show that time-evolving PDEs can be solved well with our solving architecture. The demonstrated cases include the heat equation, which is of the first derivative of time; the wave equation, which is of the second derivative of time; and the Burgers equation, which has a nonlinear term.

The first case of the time-evolving PDE problem is the heat equation, which has the first order of time derivative and is used in the description of thermal conductivity. Here, we consider the heat conduction in a 2D rectangular area with isotropic thermal diffusivity constant. The heat equation is written as follows: $$\begin{gathered} ∂u(x,y,t)∂t=a∇2u(x,y,t)+q(t),−3≤x≤3,−3≤y≤3, \\ u(x,−3)=u(x,3)=u(−3,y)=u(3,y)=0. \end{gathered}$$(1)where a is the thermal diffusivity constant, which is set as 1, and t is time. The solving time is set as 5 with a time step of 0.025. The source term, q(t), is defined as 1 in the period of 0 to 1.25 and 0 the rest of time, which corresponds to heating and cooling, respectively. The initial temperature distribution is 0. To solve this equation, the finite difference method is applied to discretize the solving domain, and the PDE is converted into the matrix form.61 Time-division multiplexing and matrix partition techniques are used to solve the matrix equation with our photonic chip. The deduction of the iteration format of the heat equation can be found in Note 7 in the Supplementary Material. The dominant operator in the heat equation is the Laplacian operator, ∇2, which is the origin of the MVM operation when solving the heat equation. The coefficient matrix of the Laplacian operator is loaded in the way described in the former section. The spectra are shown in Fig. 2(e). Figures 3(a)–3(d) give the results of the heat equation under a mesh size of 6×6. The evolution process of the solution at grid (3,3) is shown in Fig. 3(a). The blue solid line represents the simulation results, and the red dotted line represents the experimental results, which shows that the experimental results agree well with the simulation. The accuracy of the solution at different times in the evolution process is shown in Fig. 3(b), which is defined as follows: $$\begin{gathered} Ak=1−1N∑i,j|ue,i,j,k−us,i,j,k|max(us)−min(us), \\ max(us)=max1⩽i⩽o(max1⩽j⩽p(max1⩽k⩽q(us,i,j,k))), \\ min(us)=min1⩽i⩽o(min1⩽j⩽p(min1⩽k⩽q(us,i,j,k))). \end{gathered}$$(2)

Here, Ak is the overall accuracy of the solution at the kth time point and N is the grid number. ue,i,j,k and us,i,j,k mean the experimental and simulation solution of grid point (i,j) at the kth time point. To describe the accuracy of the evolution tendency of the solution, we normalized the absolute error to the range of the simulation solutions. The accuracy at t=0 is 100% because it is the initial condition. Although the accuracy fluctuates because of the random noise in practical measurement, the accuracy is above 95% during the evolution process. The accuracy at the final time point is ∼96.9%. Figures 3(c) and 3(d) show the experimental and simulation results at the time points of 0.25, 1.25, and 3.75, which agree well with each other. To describe the deviation of the experimental results from the simulation ones, the error is defined as Ei,j,k=|ue,i,j,k−us,i,j,k|max(us)−min(us).(3)

Ei,j,k is the deviation of the grid point (i,j) at the kth time point. The average error is 0.2%, 2.1%, and 3.8% at three time points. Equation (2) defines the global accuracy of the solution, whereas Eq. (3) defines the local error at each grid point, which provides insights into the correctness of the experimental results from different angles.

The time derivative of the heat equation is of the first order. Our scheme can also be applied to time-evolving PDE problems with the second order of time derivative. The case used is a wave equation describing the oscillation of a wave with constant boundary conditions, which is widely used in electromagnetics and acoustics. The solved wave equation is given as follows: $$\begin{gathered} ∂2u(x,y,t)∂t2=a2∇2u(x,y,t),−3≤x≤3,−3≤y≤3; \\ u(x,−3)=u(x,3)=u(−3,y)=u(3,y)=0. \end{gathered}$$(4)

Here, a is the wave speed, which is set as 2, and t is time. The solving time is from 0 to 1, which corresponds to 1/4 of the oscillation period. The time step is set as 0.005. The initial state is set as cos(πx/6)cos(πy/6). The finite difference method is applied to discretize the solving domain, and the PDE is converted into matrix form. The mesh accuracy here is 6×6. The deduction of the iteration format of the wave equation can be found in Note 8 in the Supplementary Material. The results of the wave equation are shown in Figs. 3(e)–3(h), with the solved equation shown at the top. The evolution process of the simulation and experimental solution at grid (3,3) is plotted in Fig. 3(e). The tendency of the experimental solution is generally correct, although the accumulated error is higher than the heat equation because of the higher time order of the wave equation. The overall accuracy of the experimental solution in the evolution process is given in Fig. 3(f), which is above 85%. The middle value of accuracy in the evolution process is 93.7%. The comparison of the experimental solution and simulation one at time points of 0.05, 0.40, and 0.9 is shown in Figs. 3(h) and 3(h). The average error is 0.08%, 4.8%, and 13% at these three time points. The accumulated error can be reduced by improving the modulation and measuring accuracy of experimental equipment in the future.

Our scheme can also be applied to nonlinear PDEs if the equation can be converted into matrix form and loaded on the photonic chip. The equation solved here is the Burgers equation, which has a first-order derivative of time and a second-order derivative of space. The Burgers equation is usually used in fluid dynamics and aerodynamics to describe the formation process of shock waves. Burgers equation is also a classic nonlinear PDE because it has a nonlinear term. Nonlinear PDEs usually have no analytic solutions, and the only way to solve them is through numerical simulation.62 As the mesh size increases, the computation costs for direct simulation of nonlinear PDEs become tremendous, which is a big burden on computing energy consumption and calculation time. Below, we show that the Burgers equation can be solved using our photonic chip.

The Burgers equation solved here is given as follows: $$\begin{gathered} ∂u(x,y,t)∂t+u(x,y,t)∇u(x,y,t)=ν∇2u(x,y,t),−1<x,y<1; \\ u(−1,y)=u(1,y)=u(x,−1)=u(x,1)=0, \end{gathered}$$(5)where u represents the speed of fluid media at a given position and time. ν is the coefficient of viscosity. To solve the PDE with our architecture, the finite difference method is applied to discretize the equation. The matrix forms of the Del operator, ∇, and the Laplacian operator, ∇2, are loaded on the photonic chip after removing the diagonal elements. The detailed deduction of the iteration format and full matrices of these two operators can be found in Note 9 in the Supplementary Material. It is worth noting that the matrices of Del and Laplacian operators have common matrix patches, which makes it easy to load them both on the microring array. The measured normalized spectra of nine rows of the microring array after tuning the voltages added on the electrodes and the normalized power matrix are given in Note 10 in the Supplementary Material, which demonstrates that the target matrix patches are loaded on the chip successfully.

For the Burgers equation solved here, without losing generality, the boundary condition is set as 0, and the initial condition is defined as (x2−1)2×(y2−1)2. The evolution time is from 0 to 3.0 s. The time step is 0.1. The mesh is set as 9×9 for a fine mesh demonstration. The solving results of the Burgers equation are given in Figs. 3(i)–3(l), with the solved equation at the top. The evolution process of the simulation and experimental solutions at grid (2,8) is presented in Fig. 3(i). The overall accuracy of the experimental solution at different times in the evolution process is given in Fig. 3(j), which is above 94% in the solving process. The comparisons of the experimental and simulation results at the time points of 0.20, 0.80, and 1.40 are given in Figs. 3(k) and 3(l). The change of the wave shape shows the process of the formation and dying out of shock waves. The average error is 4.4%, 5.2%, and 4.4% at these three time points.

The above results show the ability of our scheme to solve time-evolving PDE problems. Thanks to the high signal-to-noise ratio and low cross talk of the optical frequency comb, our scheme can handle various types of time-evolving PDEs, especially the wave equation, which has the second order of time derivative. To better illustrate the advantage of our scheme with the optical comb, we give the solving results performed with a supercontinuous laser source in Note 11 in the Supplementary Material.

The solving of PDEs is usually limited to coarse mesh size, which leads to low accuracy. In principle, our scheme can solve PDEs with unlimited mesh size without increasing the device number. Here, our scheme is applied to solving the Poisson equation with a mesh size of 12×12 experimentally, which is 5 times higher than previous reports. The general form of the Poisson equation solved here is given as ∂2u(x,y)∂2x+∂2u(x,y)∂2y=−2sin(x)cos(y),−3≤x≤3,−3≤y≤3.(6)

The equation can be written in finite difference form and converted into a matrix format. Applying the scaling operation given in Note 1 in the Supplementary Material, the iteration equation can be solved using photonic computing. The detailed deduction of the iteration format and the full matrix can be found in Note 12 in the Supplementary Material. The boundary condition is $$\begin{gathered} u(−3,y)=−u(3,y)=−sin(3)cos(y),−3≤y≤3; \\ u(x,−3)=u(x,3)=cos(3)sin(x),−3≤x≤3. \end{gathered}$$(7)

The solving results are shown in Fig. 4. Figures 4(a) and 4(b) compare the simulation and experimental results, which agree well with each other. The error distribution of the experimental solution is given in Fig. 4(c), which is defined as $$\begin{gathered} Ei,j=|ue,i,j−us,i,j|max(us)−min(us), \\ max(us)=max1⩽i⩽o(max1⩽j⩽p(us,i,j)), \\ min(us)=min1⩽i⩽o(min1⩽j⩽p(us,i,j)). \end{gathered}$$(8)

The legend shows the range of the color, which is from 0 to 1. From the color distribution, the error of the solution is low. The error value of each variable is also written in Fig. 4(c). The average error of the experimental solution is 6.5%. Figure 4(d) counts the number of grid points with different accuracy, which is defined as 1−Ei,j. The average accuracy of the experimental results is 93.5%, with a standard deviation of 0.0507.

The time-division multiplexing and matrix partition method of our scheme makes our chip flexible for solving PDEs. The microring array size of our photonic chip is 9×9. Assuming a PDE problem takes four 3×3 patches to solve, our chip can support the solving of at least two PDEs at the same time. Here, parallel solving of Laplace’s and Poisson equations is demonstrated experimentally, which shows the computing power and flexibility for PDE-solving of our scheme. The solved Laplace’s equation has the following form: ∂2u(x,y)∂2x+∂2u(x,y)∂2y=0,−3≤x≤3,−3≤y≤3.(9)$$\begin{gathered} u(−3,y)=−u(3,y)=−sin(3)cos(y),−3≤y≤3; \\ u(x,−3)=u(x,3)=cos(3)sin(x),−3≤x≤3. \end{gathered}$$(10)

The solved Poisson equation has the following form: $$\begin{gathered} ∂2u(x,y)∂2x+∂2u(x,y)∂2y=−2sin(x)cos(y),−3≤x≤3,−3≤y≤3, \\ u(−3,y)=y+3,u(3,y)=y−3,−3≤y≤3, \\ u(x,−3)=−3−x,u(x,3)=3−x,−3≤x≤3. \end{gathered}$$(11)

When the source term on the right side of the equation is set as 0, the Poisson equation transforms into Laplace’s equation. Thus, the coefficient matrix of Laplace’s equation is the same as that of the Poisson equation.

Figure 5(a) shows the top left part of the coefficient matrices of two PDEs. The matrix patches are loaded on the microring array. The normalized spectra of nine rows of the microring array are shown in Fig. 5(b). The normalized power matrix is shown in Fig. 5(c). There are two sets of identical matrix patches, which correspond to Laplace’s and Poisson equations. The MVM operation of a matrix patch is performed at one time of loading and detection. The 3D visualization of Laplace’s equation is shown in Fig. 5(d). A comparison between simulation and experimental results is shown in Fig. 5(e). The error distribution is shown in Fig. 5(f), with the error value written for each variable. The average error of Laplace’s equation is 4.4%. The definition of the error here is given in the equation. The solving results of the Poisson equation are shown in Figs. 5(g)–5(i). The average error of the Poisson equation is 4.2%. The similarity of the average error of the two equations can be understood easily. During the solving process, two PDEs share the same photonic chip and experience the same environment noise; thus, their errors are similar. Although the PDEs demonstrated here have the same derivative order and the same coefficient matrix, PDEs that can be solved in parallel are not limited to these. Only if the coefficient matrices can be partitioned into matrix patches, which are no more than nine and loaded on one chip, two PDEs can be solved in parallel by addressing the matrix patches through the time-division multiplexing method.

There are several methods to mitigate the experimental error and error accumulation in PDE solving. First, temperature is a key factor in ensuring the correctness of MVM in the experiment because each microring is tuned by its heater, controlled by electrodes. To reduce the temperature cross talk among rings, thermal isolation can be introduced in the chip design by etching the blank area between adjacent rings. Second, to reduce the cross talk among the controlling electrical currents, independent controlling wires and ground wires can be used. Third, the temperature in the laboratory can be kept low and stable to reduce environmental temperature noise. In addition, optimization algorithms can also be developed to tune the voltages of the heaters to compensate for environmental noise automatically.63 Moreover, the self-calibration technique of the microring can be applied to improve the accuracy and stability of the microring array.64

In our scheme, the matrix tuning is finished before the formal solving of the PDEs. The loading of the matrix patches does not bring any delay in the iterative solving process. Frequent loading of the input vector slices is the main source of the time cost in photonic computing. The solving speed is currently limited by the modulation speed of the modulators. The response time of VOAs used here is milliseconds, and the possible way to mitigate this tuning delay is to adopt advanced LiNbO₃ electro-optic modulators with high modulation speed of up to 100 GHz in the future, which will reduce the modulation time to the picosecond level.65^–67 Meanwhile, the high-speed InGaAs photon detectors over 220 GHz and analog-to-digital converters can also be used to match the modulation speed.68^,69 With these advanced components, the photonic computing time cost for one time of MVM is ∼18ps, which is higher than the state-of-the-art GPU (NVIDIA GeForce RTX 4090 D). A detailed discussion and comparison of the solving time between our photonic chip and GPU is provided in Note 14 in the Supplementary Material. Because the vector-matrix multiplication is finished as light propagates on the chip, the energy consumption of the photonic computing is mainly attributed to the input light energy. The output power of each microcomb tooth is ∼1mW in our experiment. The light power in the photonic system is ∼9mW for nine microcomb teeth, which is much lower than the GPU’s energy consumption of the watt level. If equipped with a high-speed modulator and photon detectors, the computing speed of our photonic system will be up to 15.3 TOPS (tera operations per second). The computing energy efficiency is 15.3 TOPS/9 mW = 1700 TOPS/W. The description of the time cost and energy consumption of the experimental system is given in Note 15 in the Supplementary Material.

If the number of the used photonic chips can be increased, our scheme can also be extended for parallel solving of more PDEs, which is discussed in Note 17 in the Supplementary Material. To increase the robustness to fabrication error of the photonic devices, topological photonic crystals can be applied in the structure design.70 It is also worth noticing that machine learning has been proven an effective tool for PDE solving, which can discover governing PDEs from data, learn effective coordinate systems, and represent solution operators.71 Developing photonic neural networks for PDE solving could be a possible direction in the future.

A scheme based on a microcomb-driven photonic chip for solving time-evolving PDEs is realized in this work. We introduce time-division multiplexing and matrix partition techniques into the photonic solving of PDEs, which can solve PDEs with fine mesh accuracy without increasing the number of photonic devices. An integrated nonlinear Kerr microcomb is employed to provide dense signal channels with a high signal-to-noise ratio, which enables our scheme to solve different kinds of time-evolving PDE problems. Various kinds of time-evolving PDE problems are solved experimentally. The heat equation with the first order of time derivative is solved with an accuracy of above 95%, the wave equation with the second order of time derivative is solved with an accuracy of around 92.5%, and the nonlinear Burgers equation with the first order of time derivative is solved with an accuracy of above 94%. An example of static PDEs, that is, the Poisson equation, is solved with a mesh size of 12×12 and an accuracy of 93.5%. Parallel solving of Laplace’s and Poisson equations is realized on a single photonic chip with an accuracy of ∼96%, which shows the computing power of our scheme for PDE solving. Our photonic chip provides a powerful platform for overcoming the challenges of solving different types of PDEs on a single photonic chip with high accuracy and in parallel. This work opens avenues to solve mathematical problems based on photonic chips and will promote the development of on-chip photonic computing.

Hongyi Yuan received his BS in 2021 and MS in 2024 degree from Beijing Institute of Technology, Beijing, China. He studied nanophotonic devices under the guidance of Professor Cuicui Lu at Beijing Institute of Technology. He is currently a PhD candidate in the School of Material Science and Engineering of Shanghai Jiao Tong University. His research interests include optical computing, nanophotonic devices, 2D materials, and polaritons.

Zhuochen Du received his BS degree in 2022 from the School of Physics at Peking University. Since 2022, he has been pursuing a PhD at the Institute of Modern Optics, School of Physics, Peking University, under the supervision of professor Xiaoyong Hu. His research interests are silicon photonics and topological photonics.

Huixin Qi received her BS degree from the School of Science, Beijing Jiaotong University, in 2020. She is currently a PhD student in the School of Physics at Peking University under the guidance of professor Xiaoyong Hu. Her research interests focus on high-speed and low-energy all-optical integrated chips.

Guoxiang Si received his BS degree from the School of Physics and Electronic Science at Shandong Normal University in 2023. He has been currently pursuing a master's degree at Beijing Institute of Technology under the guidance of Professor Cuicui Lu. His research interests focus on the field of photonic chips and intelligent algorithms.

Cuicui Lu is a professor in the School of Physics at Beijing Institute of Technology. She received her PhD from Peking University in 2015. Her research interests include topological photonics and nanophotonics.

Yan Yang is a professor at the Institute of Microelectronics of the Chinese Academy of Sciences. She received the PhD from Nanyang Technological University in 2016. Her research interests include silicon photonics and electronic-photonic integration.

Xiaoyong Hu is a professor in the School of Physics at Peking University. He received his PhD in physics from the Institute of Physics, Chinese Academy of Science, in 2004. He worked as a postdoctoral fellow with Prof. Qihuang Gong at Peking University from 2004 to 2006. Then he joined Prof. Gong’s research group. Prof. His current research interests include photonic crystals and nonlinear optics.

Biographies of the other authors are not available.