Acoustic analog computing system based on labyrinthine metasurfaces

Acoustic computing devices, including switches, logic gates, differentiator and integrator, have attracted extensive attentions in both academic research and engineering. However, no scheme of acoustic computing device with more complex functionality has been proposed, such as ordinary differential equation (ODE) solver. Here, we propose an acoustic analog computing (AAC) system based on three cascaded metasurfaces to solve the nth-order ODEs. The metasurfaces are constructed with layered labyrinthine units featuring broad amplitude and phase modulation ranges. The simulated transmitted pressure of the AAC system agrees well with the theoretical solution of ODE, demonstrating the excellent functionality. Unlike the optical ODE solver based on differentiator or integrator, whose geometry becomes more complicated for solving higher order ODE, the proposed AAC system with fixed geometry can be designed for arbitrary nth-order ODE in principle. The proposal may find applications in various scenarios such as acoustic communication, analog computing and signal processing.

constructed by layered labyrinthine units, which can provide nearly complete modulations on phase and amplitude of the incident acoustic signals 23,24 . Numerical full-wave simulations show that the transmitted pressure of the AAC system agrees well with the theoretical solution of the ODE, confirming the excellent functionality of the acoustic AAC system.

Results
Prototype for solving ODEs based on the spatial Fourier transform. The  where f(y) denotes the input function, g(y) represents the corresponding output function (equation solution), g (n) (y) is the nth-order derivative of g(y) (n ≥ 0), and a n is the constant coefficient. Imposing the FT and the inverse FT (IFT), the equation solution can be expressed as where k y is the spatial frequency. Figure 1a shows a linear space-invariant system, which can realize the solving process described by Eq. (2) [25][26][27] . In this figure, the input function propagates along the x direction, the FT block performs the spatial FT on the input function, the spatial filter applies the transfer function H(k y ), and the IFT block performs the spatial IFT. Therefore, the equation solution can be obtained by using FT and IFT blocks, along with the appropriate spatial filter.
Design of the acoustic ODE solver. To solve ODEs in acoustics, we design the AAC system as shown in Fig. 1b, where P i (y) and P t (y) are incident and transmitted signals, respectively. Analogous to the system in Fig. 1a, the AAC system consists of three sub-blocks: (i) a FT block, (ii) a well-designed SFM, and (iii) a FT block. The FT block is achieved by using a FM [with the transmission coefficient T f (y)] together with two neighboring blocks of background medium 28 . The SFM with the transmission coefficient T s (y) can function as the spatial filter. The IFT block is not considered here, since it requires the negative metamaterial 21 . Regarding these facts, the transmitted signal of the AAC system can be described as It is evident that Eq. (4) can be related mathematically to Eq. (2), provided that the incident signal P i (y), the real-space coordinate y at the SFM, and the transmission coefficient T s (y) are interpreted as the input function f(y), the spatial frequency k y , and the transfer function H(k y ), respectively. The transmitted signal P t (y) is thus proportional to the mirror image of the desired solution g(y) owing to the relation of − ∝ g y gy ( ) FT{FT[ ( )]}. As a result, the solution of the ODE can be obtained by appropriate tailoring of the transmission coefficient T s (y) to mimic the transfer function H(k y ). Unlike the optical ODE solver based on differentiator or integrator, whose geometry becomes more complicated for higher order ODE 29 , the proposed AAC system with a fixed geometry is capable for arbitrary order ODE. Besides, the AAC system can perform mathematical operations by designing the transfer function, whose functionality contains and beyond the scope of the previous analog system 21,22 .
To construct the FM and the SFM, Fig. 2a shows the layered labyrinthine unit with a width Λ = 1.2 cm and a thickness Δ = 4 cm. The operating frequency is chosen as 2500 Hz, and the background medium is chosen as air (with a density ρ = . 1 21 kg/m 3 and a bulk modulus κ = . 142 36 kPa). The labyrinthine unit consists of three tapered labyrinthine components characterized by spiral radians s 1 , s 2 and s 3 , respectively. The three spiral radians can independently vary from 0.45π to 2.55π. For an incident signal (blue arrow), the complex path-coiling of the labyrinthine unit and the induced impedance mismatch enhance the multi-scattering process, which results in the transmitted signal with tunable amplitude and phase 24 . The behavior allows the labyrinthine unit to function as an acoustic signal modulator that features broad amplitude and phase modulation ranges. Due to the property of broad modulation, the layered labyrinthine units have been used to perform mathematical operations 22 .
To demonstrate such a broad tuning, we numerical calculated the transmission coefficient T s s s ( , , ) 0 1 2 3 . The calculated transmission amplitude and phase with various spiral radians s 1 , s 2 and s 3 are shown in threedimensional color Fig. 2b and c, respectively. In Fig. 2(b,c), one coordinate point denotes one labyrinthine unit, and the value of the point denotes the transmission amplitude (or phase) of the corresponding unit. For example, to the sample unit in Fig. 2a, the transmission amplitude is 0.4 as shown by the black point in Fig. 2b, while the transmission phase is 0.63π as shown by the black point in Fig. 2c. Accordingly, the transmission coefficient of the sample unit is 0.4exp(0.63πi). In Fig. 2b and c, we can obtain a lot of labyrinthine units, whose transmission amplitudes range from 0 to 1, and transmission phases cover the entire 2π range. From these units, we can choose the desired units by minimizing Λ − T m T s s s ( ) ( , , ) 0 1 2 3 . Here, −40 ≤ m ≤ 40 is an integer, T(mΛ) is the transmission coefficient of FM or SFM. By assembling the 81 chosen units in the y direction, we can construct the FM or SFM. Therefore, the FM and SFM with the complex transmission responses can be constructed with the labyrinthine units by elaborately tailoring the spiral radians. It should be noted that the spiral radian is not a periodic function, and the range of T 0 depends on the range of spiral radian ( Fig. 2b and c). The spiral radian is thus expected to have a wide range. Nevertheless, when the spiral radian is above 2.55π, the two spiral fins will touch together, prohibiting the propagation of signal. Besides, when the spiral radian is below 0.45π, the labyrinthine component can be considered as a waveguide, whose manipulate ability is very weak. Therefore, the range of spiral radian is set as 0.45π to 2.55π.
Construct the FT block. The desired transmission coefficient of the FM can be expressed as where λ is wavelength in the air at the operating frequency, L = 4λ is the focal length. It should be noted that the focal length of the FM can be designed as a desired value 30 . Figure 3a shows the discrete transmission coefficient (circles) provided by the designed FM and the desirable continuous transmission coefficient (lines) described by Eq. (5). As can be seen, the discrete transmission phase shows a hyperboloidal profile and is in good agreement with the required transmission phase. We carry out a numerical simulation to verify the focusing functionality of the FM. For a plane incident signal, Fig. 3b illustrates the simulated intensity field distribution of the fabricated FM, where the intensity field is calculating as the square of the amplitudes of the acoustic pressure field. It is observed that the transmitted energies are focused at the focus point, which demonstrates that the fabricated FM has excellent focusing effect. The constructed FM and its two neighboring blocks of air form the FT block, and we carry out a numerical simulation to confirm its functionality of the Fourier transform. Figure 3c shows the pressure distribution of the FT block with the incident signal being P i (y) = exp(−50y 2 ). It can be seen that the transmitted pressure is strong in the central region, and tis symmetric about y = 0. Figure 3d shows the corresponding normalized transmitted pressure (circles) of the FT block. The transmitted pressure presents an axial symmetric form, and a peak of pressure is occurred at y = 0, which agrees well with the analytical result (red line). Therefore, we realize the FT block with excellent functionality based on the constructed FM. It is mentioned that the spatial frequency k y and the transverse variable y are related to each other via the relation of k y = εy with ε = − λ π L To implement the transmission coefficient in Eq. (6), we construct the SFM, whose discrete transmission coefficient is shown by circles in Fig. 4a. It is found that the transmission amplitude is in a parabolic-like profile and the transmission phase fluctuates around 0°, which agrees well with the required transmission coefficient (black lines). By inserting the designed SFM into the middle of two fabricated FT blocks, Fig. 4b shows the pressure distribution of the AAC system, where the signal P i (y) is incident on the AAC system at x = −0.5488 m (red dash line), and the transmitted signal P t (y) is obtained at x = 1.7664 m (green dash line). At the output plane, the pressure is strong around y = 0 and symmetric about y = 0. Figure 4c shows the transmitted pressure, which presents   Figure 4d shows the discrete transmission coefficient (circles) provided by the designed SFM. The transmission amplitude presents an analogous parabolic profile, and its phase decreases monotonically with increasing y. The designed transmission coefficient is consistent with the required transmission coefficient (black lines). Figure 4e shows the pressure distribution of the AAC system. It can be observed that the pressure is relatively strong around y = 0 at the output plane. Figure 4f shows the transmitted pressure of the system, which presents a valley-like form and accords with the analytical solution. The simulation results confirm that the designed AAC system could solve the third-order ODE. The designed AAC system is feasible for arbitrary incident signals at the designed operating frequency 2500 Hz (see details in Supplementary information). For example, we change the incident signal in the simulations of Fig. 4 into P i (y) [=f 2 (y) = 1280y 3 exp(−80y 2 ) − 74yexp(−80y 2 )]. The SFMs are the same for solving the second-order and the third-order ODEs with their counterparts. Figures 5a,c show the pressure distributions of the AAC system to solve the second-order and the third-order ODEs, respectively. Figures 5b,d depict the corresponding transmitted pressures. For the system to solve the second-order ODE, the pressure is weak around y = 0, where a π-phase difference is observed (Fig. 5a). On the other hand, the transmitted pressure approaches zero at y = 0, and appears a valley (peak) in the region of y < 0 (y > 0), as shown in Fig. 5b, which agree with the analytical solution g(y) = yexp(−80y 2 ). For the system to solve the third-order ODE, the pressure is weak around y = 0 (Fig. 5c), and the transmitted pressure presents a peak (valley) profile in the region of y < 0 (y > 0) (Fig. 5d). The simulated results agree with the analytical solution. Hence, the designed AAC systems for solving the ODEs are feasible for arbitrary incident signal.
Although the layered labyrinthine unit is designed with lossless assumption, the performance is also confirmed in simulation by considering visco-thermal loss since it is the inherent loss of the structure. We insert the AAC systems (Fig. 5a,c) into the COMSOL Multiphyiscs and use the thermoviscous acoustics model instead of the ideal loss-free model. For the input signal P i (y) = f 2 (y), Fig. 5e,f show the transmitted signals of AAC systems for solving the second-order and the third-order ODEs, respectively. As can be seen, the simulated results in Fig. 5e,f accord with the analytical solutions, indicating that the performance of the designed AAC system will not be severely influenced when loss is considered. Obviously, the losses would impact more or less on the amplitude and phase modulation performances. However, the operating frequency here (2500 Hz) is not high and the labyrinthine units are not resonate-based 24 . Therefore, the visco-thermal losses of the units are weak, and the functionality of the AAC system is not very sensitive to the weak visco-thermal losses.
Geometrical width of incident signal. In the simulation results (Figs 4 and 5), some discrepancies can be observed. The reasons for the discrepancies are described in Supplementary information. To improve the functionality of the AAC system, we restrict the geometrical width of incident signal by two conditions. On the one hand, the ideal FT on the incident signal P i (y) is defined as For the FT block based on FM, the integration area is limited by the width W, and the FT on incident signal P i (y) is described by which should be equal to FT ideal [P i (y)] to promise the FT block with good functionality. Comparing Eq. (9) to Eq. (8), the incident signal P i (y) must be integrable in the region (−∞, +∞), and the energy should be concentrated in the region of (−W/2, W/2), which is the first condition.
On the other hand, the FMs and SFM are low-pass filters in analogy with the optical lenses, whose cutoff frequencies are |k yc | = εW/2. The Fourier spectrum of the incident signal is S k Py ( ) FT[ ( )] y i 0 = 31 , whose all frequency contents contribute to the transmitted signal. Due to the diffraction, the Fourier spectrum S 0 (k y ) is variable in air, resulting in some high frequency contents. The FMs and SFM can complete eliminate the contents beyond k yc , leading to some discrepancies in Figs 4 and 5. Though the high frequency contents are hard to avoid in the propagation, they can be reduced by defining the geometrical width of incident signal. Here, we only consider the first FM, because a majority of high frequency contents are filtered by the first FM. The Fourier spectrum at the filtering plane (x = 0) is 32 y i y 0 2 By analyzing the Fourier spectrum S(k y ), we can explain why the results for P i (y) = f 1 (y) (Fig. 4) are better than those for P i (y) = f 2 (y) (Fig. 5). Figure 6a,b show the Fourier spectrums S(k y ) for P i (y) = f 1 (y) and P i (y) = f 2 (y), respectively. As shown in Fig. 6a, all frequency contents are in the passing band (yellow region), indicating that no content is filtered by the first FM. In this case, few discrepancies between the simulation and exact results can be observed (Fig. 4c,f). On the contrary, some frequency contents in Fig. 6b beyond |k yc |, which will be filtered by the first FM, causing some discrepancies between the simulation and exact results (Fig. 5b,d). Therefore, the distribution of Fourier spectrum S(k y ) is critical to the performance of the AAC system. To improve the accuracy of the AAC system, the ratio of low frequency contents y y should be above 98%, which is the second condition.

Discussion
In this work, we have developed an acoustic AAC system to solve the nth-order ODEs. The proposed system consists of two FMs and one SFM, which are constructed by layered labyrinthine metamaterials. The FM forms the acoustic FT block, while the SFM exhibits a complex transmission coefficient to mimic the desired transfer function of the ODE. We have demonstrated that the AAC system is capable to solve the nth-order ODEs for arbitrary incident functions. Due to the inherent merits of metasurfaces, such as low power consumption and small size, the ODE solver features fixed geometry, high efficiency and broad adjustability. The proposed system may provide an effect way to solve ODEs in acoustic domain, which could be a functional component for acoustic communication, analog computing and signal processing.

Methods
Simulations. Throughout the paper, Finite Element Method based on commercial software COMSOL Multiphysics TM 5.2a is employed for the simulations. The labyrinthine unit applied in the simulations is stiff curvature (sound hard boundaries). To calculate the transmission coefficient of the labyrinthine unit with different spiral radians (Fig. 2b and c), the plane wave radiation boundary condition is imposed on the incident and transmitted boundaries, and the periodic boundary condition are employed in the y direction. For Figs. 3-5, plane wave radiation boundaries are imposed on the outer boundaries of simulated domain to eliminate the interference from reflected wave. The designed FM and SFM with realistic structures will reflect some signals, which results in the transmitted signals containing the desired diffraction signals P t (y) and some noise signals. To decrease the reflection between the FM and SFM, their gap should be a large value and hence we choose the focal length as L = 4λ in this paper. The large L-value promises a balance between the system size and the reduced scattered signals (see Supplementary information).
Data availability. The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.