Solving ordinary and partial differential equations using an analog computing system based on ultrasonic metasurfaces

Wave-based analog computing has recently emerged as a promising computing paradigm due to its potential for high computational efficiency and minimal crosstalk. Although low-frequency acoustic analog computing systems exist, their bulky size makes it difficult to integrate them into chips that are compatible with complementary metal-oxide semiconductors (CMOS). This research paper addresses this issue by introducing a compact analog computing system (ACS) that leverages the interactions between ultrasonic waves and metasurfaces to solve ordinary and partial differential equations. The results of our wave propagation simulations, conducted using MATLAB, demonstrate the high accuracy of the ACS in solving such differential equations. Our proposed device has the potential to enhance the prospects of wave-based analog computing systems as the supercomputers of tomorrow.


Architecture and working principle of the analog computing system (ACS)
Architecture.Our proposed ACS (Fig. 1a) is made up of three key components: an Ultrasonic Fourier Transform (UFT) block, a spatial filtering metasurface (SFM), and another UFT block.The pressure fields at the input and output planes of the ACS are P I x, y and P O x, y , respectively.The side length of the entire ACS' square cross-section is L.
Referring to Fig. 1b, each UFT block consists of three main parts: a fused silica substrate layer, the ultrasonic metalens, followed by another fused silica substrate layer.If the input pressure field of the ultrasonic wave that is made to pass through the UFT block is P , the output pressure field obtained at the other end of the UFT block is P FT , which has been shown to be proportional to P 's Fourier transform 37 .A key condition for obtaining the UFT through this block is that both the thickness of the substrate layers and the focal length of the metalens must be f 37 .The metalens' thickness is t m , as shown in Fig. 1b.
Our proposed ACS is designed to operate at a frequency of f wave = 1.7 GHz, which is a high ultrasonic fre- quency.This enables greater compactness, making it easier to integrate the ACS into CMOS-compatible chips.Each substrate layer has a thickness of f = 1.0886 mm and is made of fused silica (which we chose for its isotropy as a material).In fused silica, the speed of ultrasonic waves is v wave = 5880 m s -1 , from which we can calculate the wavelength to be = v wave /f wave = 3.46 µm.
In Fig. 2a, the metalens is made up of several unit cells, each of which has a thickness of t m = 16 µm and a square cross-section of side length 3 µm (a subwavelength feature).Each unit cell (Fig. 2b) is composed of a square cuboid made of Si with a cylindrical post made of SiO 2 embedded in it.According to the theoretical working principle of the ACS, the ultrasonic metalens ought to obey a paraboloidal phase profile such that the pressure field P FT would be proportional to the FT of P .Due to the limited number of distinct unit cells available, however, discretization is required.Therefore, the cylindrical post radius of each unit cell must correspond to the interpolated phase shift at that point (after discretization).The process of interpolation transforms the ideal phase map (Fig. 2c) into the discretized phase map (Fig. 2d) that is later used for phaseto-radius mapping.
In addition, the transmission coefficient function T x, y of the SFM must correspond to the transfer function (TF) H k x , k y required to solve a particular ordinary or partial differential equation.

Working principle.
Uy and Bui 37 have previously determined that the input P and the output P FT of the UFT block are approximately (see Table 1 for list of approximations) related by where j is the imaginary unit, is the wavelength, k is the wavenumber, and the operator F denotes the FT.Now, let f x, y and g x, y be the input and output, respectively, of a particular ODE or PDE.It can be shown that f x, y and g x, y are mathematically related by the equation where the operator F −1 denotes the inverse FT and H k x , k y is the transfer function for a certain ODE or PDE.
At first glance, Eqs. ( 2) and (3) might seem to suggest that the ACS cannot compute the accurate result.For one, the UFT block yields an output P FT that is only proportional to-but not actually equal to-the FT of the input P .Moreover, it is essential to note that the correct output is obtained by taking the inverse FT of H k x , k y F f x, y , whereas the second UFT block calculates the FT (not the inverse FT).However, these concerns do not actually hinder the ACS from yielding the desired output.In fact, the mirror image of the correct output g x, y is given by the equation (2) Table 1.List of approximations required to achieve P FT ∝ F {P}. www.nature.com/scientificreports/which is the additive inverse of the magnitude of the ACS' output P O x, y .Therefore, to achieve the desired result g x, y , we simply need to take the mirror image of the ACS' output P O x, y and subsequently obtain its additive inverse.It is also important to recognize that the transmission coefficient can only have amplitudes of less than or equal to 1.If this condition is satisfied by H k x , k y , then all is well, and H k x , k y would also be the transmission coefficient function.However, this is not necessarily satisfied by all transfer functions.In such cases, we would have to normalize the transfer function such that it satisfies this condition.Consequently, we would also not obtain the exact output g x, y but rather a scaled version of it.Nevertheless, we can recover the desired result by appropriately rescaling the output.

Results
A function that is space-limited has non-zero magnitude values only within in a finite region in the space domain, whereas a function that is bandlimited is one that has a finite spectral width that includes all spatial frequency components with non-zero magnitude values.It should be noted that it is not possible for a function to be both space-limited and bandlimited 38 .Furthermore, functions that are space-limited but not bandlimited, such as the rect function, are not of particular interest in the context of solving differential equations.In this paper, we thus consider two main kinds of functions: (1) bandlimited but not space-limited and (2) neither space-limited nor bandlimited.
In the wave propagation simulations, we used the Sinc and Gaussian functions as archetypal examples for each kind.The Sinc functions follow the general form f (x) = sinc(x/w) , where the parameter w ∈ R + is an indication of the Sinc function's geometric spread.The Gaussian functions, on the other hand, take the form f (x) = exp −πx 2 /γ 2 , where the parameter γ ∈ R + serves as a measure of the Gaussian function's geometric spread.Additionally, we chose the root-mean-squared error (RMSE) after normalization as the error metric used to assess the accuracy of the ACS' output, relative to the analytical solution.

Ordinary differential equations (ODE). Mathematical basis.
Generalizing the derivative property of the FT 39 , Consider a general n th -order inhomogeneous ODE Taking the FT of both sides of Eq. ( 6), we obtain from which we can deduce that Note that Eq. ( 8) can be re-expressed as By definition, the spatial frequencies are k x = −2πx/ f and k y = −2πy/ f , so we have dk x = − 2π/ f dx and dk y = − 2π/ f dy .Therefore, using these substitutions and then replacing x with −x and y with −y, It is apparent from Eq. (10) that f (x) is the input and is the transfer function needed to solve nth-order inhomogeneous ODEs of the form given in Eq. ( 6) 9 .
Simulation results.In the simulations, we used the ODE www.nature.com/scientificreports/For the first type of function, the solution g(x) is a Sinc function with parameter w = 18 .The RMSE after normalization, based on our simulation, was 0.00737.From Fig. 3a, we can observe that the analytical solution and the ACS' output are in excellent agreement, especially near the center.There are some minor discrepancies towards the edges, but these are actually expected since we know that the FT is only obtained in the paraxial region (see approximations listed in Table 1).Furthermore, metalens aberration-as a result of discretization of the metalens' phase profile-also contributes, albeit very minimally, to the observed deviations.It should also be noted that there is some undersampling due to truncation of the input function f (x) as well as aliasing arising from bandlimiting of the output function g(x).
The simulation for the second type of function was conducted with a Gaussian function with parameter γ = 48 as the solution g(x) .In this case, the RMSE after normalization was found to be 0.00975.Figure 3b shows that the output of the ACS and the analytical solution agree very well with each other.Granted, there are small lobes towards the edges, largely due to metalens aberration as well as the paraxial approximation requirement not being met.There are also minor errors associated with undersampling of the input and bandlimiting of the output.Be that as it may, the output nonetheless matches the analytical solution very well.
Optimization of accuracy.Figures 3c and d show the relationship between the RMSE and the geometric spread parameters w and γ , respectively, for solving ODEs.We can observe that the RMSE initially decreases then increases as w or γ is increased.This two-part trend in the RMSE is mainly caused by two factors: (1) the level of aliasing and undersampling and (2) the validity of the paraxial approximations.
Firstly, as w or γ is initially increased, there is reduced aliasing in the output plane of the first UFT block and reduced undersampling in the input plane of the second UFT block, resulting in a fall in the RMSE.However, as w or γ continues to increase past a certain threshold, there is increased undersampling in the input plane of the first UFT block and increased aliasing in the output plane of the second UFT block that drive the observed rise in the RMSE.
Secondly, the paraxial approximations initially become more valid as w or γ increases, then it becomes less valid as w or γ increases further.As w or γ initially increases, the energy becomes less highly concentrated near the center of the first UFT block's input plane and less extensively spread out in the first UFT block's output plane.The energy also becomes less extensively spread out in the second UFT block's input plane and less highly concentrated near the center of the second UFT block's output plane.Then, as w or γ increases further, the energy becomes more extensively spread out in the first UFT block's input plane and more highly concentrated near the center of the first UFT block's output plane.Moreover, the energy becomes more highly concentrated near the center of the second UFT block's input plane and more extensively spread out in the second UFT block's output plane.
Therefore, to optimize accuracy, the input function f (x) can be scaled parallel to x such that the geometric spread parameter w or γ takes on a moderate value.The output P O x, y of the ACS can then be appropriately rescaled back to obtain the desired solution g(x).
Refer to Supplementary Figs.S2 to S5 for additional diagrams in support of the above explanation.

Partial differential equations (PDE). Mathematical basis. Consider the partial differential equation
Taking the FT of both sides of Eq. ( 13), which can be re-arranged to get In its integral form, Eq. ( 15) can be re-written as Substituting dk x = − 2π/ f dx and dk y = − 2π/ f dy and replacing x with −x and y with −y then yield Therefore, f x, y is the input and is the transfer function needed to solve PDEs of the specified form 28 .
Simulation results.In the simulations, we used the PDE (13) κ 1 ∇ 2 g x, y + κ 2 g x, y = f x, y .The simulation for the first type of function was carried out with a two-dimensional Sinc function with parameter w = 18 as the solution g x, y .The RMSE after normalization, based on our simulation, was 0.00643.We can observe from Fig. 4a that the ACS' output and the analytical solution are in excellent agreement, particularly near the center.There are some small deviations towards the edges, which can be attributed to the paraxial approximation not being met.In addition, metalens aberration also contributes to the discrepancies.There is also some undersampling due to truncation of the input function f x, y as well as aliasing arising from bandlimiting of the output function g x, y .
For the second type of function, the solution g x, y is a two-dimensional Gaussian function with parameter γ = 18 .In this case, the RMSE after normalization was found to be 0.00266.Figure 4b shows that the output of the ACS and the analytical solution agree very well with each other, especially in the paraxial region for the same reason cited above.Due to metalens aberration as well as the paraxial approximation requirement not being met, there are some small lobes towards the edges.What is more, undersampling of the input f x, y and aliasing of the output function g x, y also contribute to the observed discrepancies.Nonetheless, the output for the significant magnitude values matches the analytical solution very well.
Optimization of accuracy.Figure 4c and d show the relationship between the RMSE and the geometric spread parameters w and γ , respectively, for solving PDEs.We can observe that as w or γ is increased, the RMSE initially falls and subsequently rises.This trend in the RMSE can be largely ascribed to two key factors: (1) the level of undersampling and aliasing and (2) the paraxial approximations' validity.
Firstly, as w or γ is initially increased, there is less aliasing in the first UFT block's output plane and less undersampling in the second UFT block's input plane, so the RMSE decreases.However, as w or γ continues to increase beyond a certain value, there is increased undersampling in the first UFT block's input plane and increased aliasing in the second UFT block's output plane, resulting in the observed increase in the RMSE.
Secondly, the paraxial approximations initially become more valid as w or γ initially increases, then it becomes less valid as w or γ increases further.Initially, as w or γ increases, the energy becomes less highly concentrated near the center of the input plane of the first UFT block and less extensively spread out in the output plane of the first UFT block.The energy also becomes less extensively spread out in the input plane of the second UFT block and less highly concentrated near the center of the output plane of the second UFT block.Then, as w or γ increases further, the energy becomes more extensively spread out in the input plane of the first UFT block and more highly concentrated near the center of the output plane of the first UFT block.Additionally, the energy becomes more highly concentrated near the center of the input plane of the second UFT block and more extensively spread out in the output plane of the second UFT block.
Thus, accuracy can be optimized by scaling the input f (x) parallel to x so that the geometric spread param- eter takes on a moderate value.We can then appropriately rescale back the ACS' output P O x, y to obtain the desired solution.
Refer to Supplementary Figs.S6 to S9 for additional diagrams supporting the above explanation.

Conclusion
This paper introduces an analog computing system (ACS) that uses ultrasonic waves and metasurfaces to solve ordinary and partial differential equations.Through our simulations, we have clearly demonstrated the ability of the proposed ACS to yield highly accurate results when solving both types of differential equations involving both types of functions.In contrast to other studies in existing literature, a key contribution of our paper is the exploration of how the accuracy of the ACS' output may be optimized through the selection of an appropriate (moderate) value of the geometric spread parameter w or γ.Our study's findings are anticipated to advance the development of wave-based analog computing systems, potentially surpassing the constraints of digital computers.This has far-reaching implications for the fields of computing and signal processing, hopefully laying the foundation for technological breakthroughs in the future.

Methods
Ultrasonic metalens designing process.The process of designing the ultrasonic metalens is detailed below.Refer to Supplementary Fig. S1 for a flowchart summarizing the method.
To start, it is crucial to conduct unit cell simulations in order to establish the relationship between the radius of a cylindrical post and the associated phase shift for that particular unit cell.This correlation (Phase-to-radius Mapping) serves as a reference point for subsequent steps.Next, an array of the Ideal Phase Map can be generated.It consists of phase values at sampled points following the paraboloidal phase profile required to achieve P FT ∝ F {P} theoretically 37 .The next step is to use the MATLAB function interp1 to inter- polate the nearest available phase value from the unit cell simulations, and this results in the Discretized Phase Map, consisting of phase values that have a corresponding radius from the unit cell simulations.Subsequently, the phase-to-radius mapping can be used to create the Radius Map, an array of radius values at each sampled point.Finally, the MATLAB function viscircles can be used to generate a figure of the metalens, comprising unit cells with cylindrical posts whose radii correspond to the radius at that point (according to the Radius Map).     1 were applied), we conducted semianalytical simulations of the propagation of ultrasonic waves through our proposed ACS.The output of each simulation was then compared with the analytical solution.
To conduct these simulations, we used MATLAB to implement the code because it is far less computationally costly as opposed to Finite Element Method (FEM) software like COMSOL Multiphysics.This is to enable us to carry out simulations involving arrays of significantly larger size-key to understanding the proposed ACS' true capabilities.
Propagation within each UFT block.Through an FFT-based convolution approach, we can obtain the pressure field P M− x, y at the plane before the ultrasonic metalens.This approach involves convolving the zero-padded input pressure field array P S (ξ , η) with the convolution kernel To avoid circular convolution errors, the N × N array P S (ξ , η) has to be padded by at least N − 1 zeros 37,43 .According to convention, the convolution kernel array and the zero-padded pressure field array are of the same size 37,43 .The N × N subarray at the center of the larger array generated as the output of FFT-based convolution is P M− x, y .
Subsequently, we apply the phase shift due to the discretized metalens to obtain N × N array P M+ x, y , which represents the pressure field at the plane after the metalens.This can be done by performing an elementwise multiplication of the N × N array P M− x, y and the N × N array exp iφ discretized .Note that φ discretized is the metalens' discretized phase profile after the interpolation step in the ultrasonic metalens designing process.
Following this, we use FFT-based convolution to convolve the zero-padded array P M+ x, y with the con- volution kernel in order to obtain the N × N output pressure field array P O (u, v).
The convolution kernels in Eqs. ( 21) and ( 22) were derived from exact solutions to the Kirchhoff-Helmholtz Integral 37 .
Propagation through the SFM.The SFM theoretically applies the transfer function (TF) needed to solve a particular ODE or PDE.To simulate ultrasonic wave propagation through the SFM, we perform element-wise multiplication of the output array of the first UFT block and the transmission coefficient array T x, y of the SFM.The result is then the input array for the second UFT block.We subsequently repeat the same process described above for the simulation of wave propagation through a UFT block, which ultimately produces the output P O x, y of the proposed ACS.
Simulation parameters.The process of selecting the most appropriate simulation parameters involves three key considerations.Firstly, convolution requires that the spacing between adjacent metalens unit cells and that between the sampled points of the pressure fields must be the same 37,43 .Furthermore, an appropriate length L for the cross-section of the UFT block should be chosen, keeping in mind that the geometric spread parameter w or γ must take on a moderate value so that the significant space and spatial frequency components are within the sampled array bounds (appropriate truncation and bandlimiting).Finally, the focal length f ought to satisfy the condition which can be derived by considering the sampling requirements of the convolution kernels' exponential phase term 36,37,[40][41][42][43][44] .
With these considerations in mind, the values of the parameters used in the simulations are presented in Table 2.

Figure 1 .
Figure 1.Schematic of the ACS and the UFT block.(a) The figure features a schematic of the proposed ACS, which consists of three main parts: a UFT block, an SFM, and another UFT block.(b) The figure features a schematic of the UFT block, which has three key components: a substrate layer, the ultrasonic metalens, and another substrate layer.

Figure 3 .
Figure 3. Ordinary Differential Equation (ODE) Simulations.(a) Sinc function: Input (left) and output (right) magnitude profiles for w = 18 .The blue line with circled data points indicates the simulated output of the ACS, whereas the orange line shows the analytical solution.(b) Gaussian function: Input (left) and output (right) magnitude profiles for γ = 48 .The blue line with circled data points indicates the simulated output of the ACS, whereas the orange line shows the analytical solution.(c) Relationship between the RMSE and the parameter w for the Sinc function.(d) Relationship between the RMSE and the parameter γ for the Gaussian function.

Figure 4 .
Figure 4. Partial Differential Equation (PDE) Simulations.(a) Sinc function: Input (left) and output (right) magnitude profiles for w = 18 .The blue line with circled data points indicates the simulated output of the ACS, whereas the orange line shows the analytical solution.(b) Gaussian function: Input (left) and output (right) magnitude profiles for γ = 18 .The blue line with circled data points indicates the simulated output of the ACS, whereas the orange line shows the analytical solution.(c) Relationship between the RMSE and the parameter w for the Sinc function.(d) Relationship between the RMSE and the parameter γ for the Gaussian function.

Table 2 .
Values of simulation parameters used.