Acoustic metasurfaces are two-dimensional materials that impart non-trivial amplitude and phase shifts on incident acoustic waves at a predetermined frequency. While acoustic metasurfaces enable extraordinary wavefront engineering capabilities, they are not developed well enough to independently control the amplitude and phase of reflected and transmitted acoustic waves simultaneously, which are governed by their geometry. We aim to solve the inverse design problem of finding a geometry to achieve a specified set of acoustic properties. The geometry is modeled by discretizing the continuous space into a finite number of elements, where each element can either be filled with air or solid material. Full wave simulations are performed to obtain the acoustic properties for a given geometry. It is computationally infeasible to simulate all geometries. To address this challenge, we develop an experimental design-based algorithm to efficiently perform the simulations. The algorithm starts with a few geometries and adaptively adds geometries to the set, such that they fill the entire space of the desired acoustic properties using a small fraction of the possible geometries. We find that the geometry needs to have at least 7 × 7 elements to obtain any given acoustic property with a tolerance of 5.4% of its maximum range. This is achieved by simulating 24 000 geometries using the proposed algorithm, which is only $4.2\xd710\u22129%$ of the 563 × 10^{12} possible geometries. The method provides a general solution to the inverse design problem that can be extended to control more acoustic properties.

Acoustic metasurfaces are artificial sheets of material with subwavelength thickness formed by different unit cells to modulate propagating sound waves.^{1,2} The ultrathin feature of the metasurfaces is a result of designing the wave–matter interactions and wave manipulations within a deep subwavelength space. This deep subwavelength wave modulation functionality has enabled many applications of metasurfaces, including beam steering and forming,^{3,4} perfect sound absorbers,^{5,6} acoustic illusions,^{7} holography,^{8} signal multiplexing for high-speed communications,^{9} particle levitation, motion control,^{10} and ground cloaking,^{7,11} which are critical for the development of underwater exploration, noninvasive biomedical treatments, architectural designs, and noise control. In most of these applications, independent control of the acoustic amplitude and phase profile is crucial.

Despite the importance of metasurfaces, most existing designs lack the capability to control the acoustic amplitude and phase independently. A recent experimental study has suggested the introduction of controlled energy loss in the design of metasurfaces to realize a decoupled modulation of sound amplitude and phase.^{12} However, the intrinsic loss resulted from this design method limits the efficiency of acoustic wave modulation. Yet, the additional loss parameter increases the complexity of the problem when compared to previous geometry-based metasurface designs that prevent the development of a systematic analytical design methodology for arbitrary acoustic metasurfaces with independent amplitude and phase controllability. Moreover, the loss design demonstrated can only be used for the modulation of reflected acoustic waves.^{12}

Here, we formulate the problem as an inverse design problem, which, if solved, will enable independent phase and amplitude modulation. The acoustic properties can be changed by changing the geometry of the metasurface. Because of the large number of possible geometries, performing simulations and predicting the acoustic properties on all of them are practically infeasible. For example, a 5-by-5 meshed unit cell can produce more than 33 × 10^{6} possible geometries. This makes it challenging to understand the dependence of the acoustic properties on geometry using simulations. Researchers have solved inverse design problems by approximating the simulation model using a cheap-to-evaluate machine learning model.^{13,14} However, the acoustic metasurface design poses new challenges. First, the input space is huge and discrete. Second, because of the small-sized elements, the convolutional neural network-based techniques are not very effective. Third, the geometries have several redundancies, which are not easy to incorporate in statistical and machine learning models.

In this work, we develop an experimental design-based approach^{15} to address the above challenges. The proposed adaptive space-filling algorithm^{16} starts from a set of few initial geometries and then uses them to select geometries that efficiently fill the space of the desired acoustic properties. A key feature of the algorithm is that it is model-free and iterates through only a small fraction of the possible number of geometries to find those that span the entire space of the desired acoustic properties. Then, given the desired properties, we can pick the geometry that approximately achieves them.

To characterize all possible geometries for the design of an acoustic metasurface, we discretized the two-dimensional continuous space (a unit cell) into a finite number of elements. Each element can either be filled with air (where acoustic wave propagates) or solid material (no sound wave can propagate). Figure 1 shows a unit cell with $m=n=4$, resulting in 16 elements, where *m* and *n* are the number of elements perpendicular and parallel to the waveguide, respectively. The height of the unit cell is chosen to be one-third of the wavelength $(\lambda )$, as it cuts off higher order acoustic waveguide modes.

Our problem has some resemblance to topology optimization.^{17–19} However, we are not performing an optimization to get a specific target. Instead, we find a set of geometries that can span the space of acoustic properties, so that when a target is specified, we can quickly pick one from the set to achieve the target approximately within the specified tolerance. Thus, our approach would be more useful than topology optimization when several targets are desired, and the inverse solutions need to be found quickly.

The full wave simulations of acoustic wave modulation by a single metasurface unit cell are performed in the Pressure Acoustic Module of COMSOL Multiphysics. Because we assume that the coupling between the unit cells is zero, each individual unit cell can be simulated separately.^{3,4,20–26} To calculate the amplitudes and phases of propagating acoustic waves transmitted and reflected by each unit cell, we apply a one-dimensional acoustic waveguide, as shown in Fig. 1. The top and bottom boundaries are sound hard boundaries, which confine the acoustic energy inside the one-dimensional waveguide. An acoustic plane wave is incident from the left end of the waveguide, whose phase will be set such that the incident acoustic phase is zero at the left boundary of the metasurface unit cell. The transmission and reflection are calculated by the pressure fields along two vertical lines with a distance to be an integer multiple of the targeted acoustic wavelength $(n\lambda )$ from the right and left boundaries of the metasurface unit cell, respectively, to avoid near-field effects. Note that the incident acoustic wave is known, and the reflected acoustic wave equals the difference of the calculated pressure field and the incident acoustic wave at the same line. Because the distance between the lines used for the calculation of transmission and reflection and the boundaries of the unit cell is an integer multiple of wavelength, the acoustic phases retrieved are exactly the same as the phases at the boundaries of the unit cell, and therefore, identical to the transmitted and reflected phases. The amplitudes and phases of acoustic waves transmitted and reflected from the unit cell are the acoustic properties that we intend to control.

We need to choose an optimal grid size, i.e., optimal values for *m* and *n*, which will suffice to modulate the amplitude to any value in $[0,1]$ normalized by the incident wave, and the phase to any value in [−180°, 180°]. If the grid size is too small, there may not be enough number of distinct geometrical structures that can span the space of the desired acoustic properties. If the grid size is too large, it may lead to an unnecessarily large number of geometries and makes it very expensive to physically construct the acoustic metasurface.

Once the grid size is selected, our objective is to find the geometry that can produce the desired acoustic properties. Note that the reflection and transmission amplitudes are related by the conservation of energy. For lossless acoustic elements

where *T*, *R*, and *I* are the amplitudes of transmission, reflection, and incidence, respectively. Because of this, we can independently control only three of the acoustic properties—transmission phase, reflection phase, and either the transmission amplitude or the reflection amplitude. Thus, the problem can be stated more precisely as follows. Let ** g** be the geometry of the acoustic metasurface, and

**denote the corresponding three-dimensional acoustic property. Their functional relationship is given by $x=f(g)$, where $f(\xb7)$ in our case represents the full-wave simulation code.**

*x*The inverse design problem is to find $g*$ given an $x*$ by solving

where $||\xb7||$ is a norm such as the $L\u221e$ norm and $Gm\xd7n$ is the set of $2mn$ possible geometries of grid size *m* ×* n*. The solution to this optimization problem need not be unique, but this non-uniqueness is not an issue for us. We just need one $g*$ that will produce the desired property $x*$. The main issue with this optimization problem is twofold: $f(\xb7)$ is computationally expensive and $Gm\xd7n$ is huge for large grids, and therefore, solving (2) is very hard. Furthermore, we also need to find the smallest *m* and *n*, such that

where Δ is a specified tolerance.

We performed a set of preliminary simulations with grid sizes of up to 16 elements, which corresponds to a maximum of $216=65\u2009536$ distinct geometries. This included grid sizes such as $4\xd72,8\xd72,2\xd78$, and 4 × 4. Figure 2 shows the output from the simulation of a 4 × 4 grid as a scatterplot matrix (for example, the sub-plot in the first row and the third column shows the reflection amplitude on the vertical axis and the transmission amplitude on the horizontal axis). We can see an almost perfect relationship between the amplitudes of reflected and transmitted waves, which is expected because of (1). Clearly, 16 elements are not sufficient to cover the property space and satisfy (3). Thus, we need to try larger grid sizes.

An important observation from the foregoing preliminary study is that 4 × 4 grid produced the densest set of points in the space of acoustic properties (the simulation results for all grid sizes can be found in Sec. A of the supplementary material). This indicates that the optimal grid size is likely to be square-shaped, that is, *m* =* n*. Thus, we can now focus on simulating 5 × 5, 6 × 6, etc. However, the number of simulations needed is computationally infeasible for larger grid sizes. For example, even for 5 × 5, $225\u224833.55$ × 10^{6} geometries are possible, which requires several months of simulation time on an ordinary laptop.

Interestingly, we observed that due to certain geometrical symmetries, some distinct geometries have the same acoustic properties. The description of each redundancy and algorithms to discard redundant geometries are provided in Sec. B of the supplementary material. However, eliminating redundancies is still insufficient to simplify higher grid-sizes. For example, for the 5 × 5 grid-size, even with a 99% reduction in the number of geometries under consideration, it will be infeasible to simulate the remaining 1%, or 335 500 geometries.

As we cannot simulate all possible non-redundant geometries, we propose the following approach to solve the inverse design problem. We will find a set of *N* geometries $G={g1,\u2026,gN}$, such that

where $xi=f(gi)$ for $i=1,\u2026,N$ and $X=[\u2212180\xb0,180\xb0]2\xd7[0,1]$. Then, for any $x*$, we will find the geometry $gi*$, where $i*=arg\u2009mini||x*\u2212xi||$. In other words, we will find the geometry that will produce the desired acoustic properties with a worst-case error of Δ.

In the experimental design literature,^{15} finding $X={x1,\u2026,xN}$ that minimizes the left side of (4) is called a minimax design.^{16,27} However, the problem here is different and much more complex. The minimization should be done over the geometries $G$ and not $X$. Hence, the existing algorithms for finding minimax designs cannot be used.^{28} Therefore, we propose the following algorithm to solve the problem.

Our algorithm begins by simulating the acoustic properties of a few randomly generated geometries. Let $G0={g1,\u2026,gI}$ denote the set of *I* initial geometries, and $X0={x1,\u2026,xI}$ denote the set of their corresponding acoustic properties. Then, the geometry corresponding to the most sparsely located acoustic property (say $xsparse$) is “perturbed,” i.e., changed slightly. See Sec. C in the supplementary material for details on perturbation. This will generate a new geometry (say $gnew$) whose acoustic property (say $xnew$) is likely to be in the sparsely populated region of the properties space, thereby filling the “gaps” in the space. Note that $gnew$ is discarded if it is redundant with a geometry that has already been simulated. The most sparsely located property is the one being the farthest from its nearest-neighbor on either side of the property along each dimension. Let $R={r1,\u2026,rI}$ be the vector of the distances to the farthest nearest-neighbor for each acoustic property and $sparse=arg\u2009maxiri$. Then, the geometry corresponding to $xsparse$ is perturbed. See Sec. D in the supplementary material for details on finding the most sparsely located property, which is motivated by an algorithm used for crystal structure prediction.^{29}

This new acoustic property $xnew$ corresponding to $gnew$ is accepted only if it is farther than a threshold distance *t* (defined later) from its nearest neighboring property. Otherwise, $xsparse$ is perturbed again. The purpose of the threshold distance *t* is to avoid selecting geometries with very similar properties, which reduces the computational cost of the algorithm. The threshold distance is updated in each iteration, irrespective of acceptance or rejection of $gnew$. If *t _{i}* is the threshold distance in the

*i*th iteration, and $dmin,i$ is the distance of $xnew$ to its nearest neighbor in this iteration, then the threshold distance for the next iteration is given by $ti+1=0.5(ti+dmin,i)\u2009for\u2009all\u2009i>1$. We stop the algorithm when the largest gap in the desired acoustic properties space is less than Δ. We call our algorithm “minimax perturbation algorithm,” which is summarized as Algorithm 1.

1: Input $G,X,tol$ |

2: $i\u2190I$ |

3: Compute $R={r1,\u2026,rI}$ |

4: Compute $D={d1,\u2026,dI}$ {Distance to the nearest neighbor of each acoustic property} |

5: $t\u2190$ mean() D |

6: $Gapmax\u2190max(D)$ |

7: while $Gapmax>\Delta $ do |

8: sparse $\u2190arg\u2009maxiri$ |

9: Perturb $gsparse$ to generate $gnew$ |

10: if $gnew$ is redundant, goto (9) |

11: $xnew\u2190COMSOL_simulation(gnew)$ |

12: Compute $dmin,\u2009Distance\u2009to\u2009the\u2009nearest\u2009neighbor\u2009of\u2009xnew$ |

13: if $dmin>t$ then |

14: $G\u2190$ append($G,gnew$) |

15: $X\u2190$ append($X,xnew$) |

16: $i\u2190i+1$ |

17: Update $R,D$ |

18: $Gapmax\u2190max(D)$ |

19: end if |

20: $t\u21900.5(t+dmin)$ |

21: end while |

22: Output $G,X$ |

1: Input $G,X,tol$ |

2: $i\u2190I$ |

3: Compute $R={r1,\u2026,rI}$ |

4: Compute $D={d1,\u2026,dI}$ {Distance to the nearest neighbor of each acoustic property} |

5: $t\u2190$ mean() D |

6: $Gapmax\u2190max(D)$ |

7: while $Gapmax>\Delta $ do |

8: sparse $\u2190arg\u2009maxiri$ |

9: Perturb $gsparse$ to generate $gnew$ |

10: if $gnew$ is redundant, goto (9) |

11: $xnew\u2190COMSOL_simulation(gnew)$ |

12: Compute $dmin,\u2009Distance\u2009to\u2009the\u2009nearest\u2009neighbor\u2009of\u2009xnew$ |

13: if $dmin>t$ then |

14: $G\u2190$ append($G,gnew$) |

15: $X\u2190$ append($X,xnew$) |

16: $i\u2190i+1$ |

17: Update $R,D$ |

18: $Gapmax\u2190max(D)$ |

19: end if |

20: $t\u21900.5(t+dmin)$ |

21: end while |

22: Output $G,X$ |

We used the minimax perturbation algorithm to simulate geometries for 5 × 5 and monitored the maximum gap $maxx\u2208Xmini=1:N||x\u2212xi||$ in the property space. We chose a Δ around 5% of the individual ranges of the properties. After simulating about 10 000 geometries, we found that the maximum gap is consistently much higher than the tolerance Δ. The same thing happened with 6 × 6 grids. So, we moved onto 7 × 7 grids. After simulating about 24 000 geometries, we found that the maximum gap reduced to about 5.4%, which is a satisfactory result. Thus, we conclude that 7 × 7 grid-size is the optimal grid size. Figure 3 shows the scatterplot of the acoustic properties for the 7 × 7 grid size. Note that 24 000 geometries are only $4.2\xd710\u22129%$ of the $249\u2248563$ × 10^{12} possible geometries. This shows the efficiency of the proposed algorithm for solving inverse design problems.

We performed a simulation exercise to demonstrate the utility of our methodology. We aim to identify the geometries to create a metasurface that produces two unique focused acoustic waves on the transmission/reflection sides of the metasurface. These unique focal points require a custom reflection/transmission phase profile that maintains a constant reflection/transmission amplitude while satisfying energy conservation. The required phases and amplitudes are shown by solid red points in Fig. 3. Since the 24 000 geometries identified by our algorithm fully span the space of the acoustic reflection and transmission coefficients, we quickly selected ten geometries with the desired reflection/transmission coefficients to produce the two distinct focal points on either side of the metasurface. The phase profile for the unique focal points is calculated from Li *et al.*^{4} Figure 4 shows the focal points generated with the waves reflected (left) and transmitted (right) from the acoustic metasurface (center).

In conclusion, we have developed an experimental design-based algorithm to find the acoustic metasurface geometries that efficiently fill the space of the acoustic properties. These geometries can be stored and quickly queried for a given acoustic property, thus approximately solving the inverse design problem. The proposed algorithm seems to be suitable for the inverse design of acoustic metasurface because (i) it is capable of exploring the huge space of geometries using relatively few simulations, which is an advantage over the traditional optimization-based approaches and (ii) it is a model-free method and can easily incorporate symmetries and redundancies of geometries without additional feature engineering, which is an advantage over the modern machine learning-based approaches. For this study, we have set the mesh-size (7 × 7) to fully span the space of the transmission and reflection coefficients within a 5% margin of error. Future studies could further decrease the mesh size to cover this space with a higher density of points and create resonant structures that achieve negative effective material properties.

Acoustic metasurfaces with independent phase and amplitude control significantly improve the wavefront engineering capabilities of metasurfaces for beam steering, acoustic holography, and particle manipulation for biomedical, architectural, and underwater acoustic applications.

See the supplementary material for (a) preliminary simulation results with the 4 × 4 grid, (b) algorithms to discard redundant geometries, (c) details on perturbing a geometry, and (d) method to find $xsparse$.

This research was supported by the U.S. National Science Foundation under Grant Nos. ECCS-2102129 and CMMI-1921646.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**Arvind Krishna:** Conceptualization (equal); Formal analysis (equal); Investigation (equal); Methodology (equal); Project administration (equal); Software (equal); Validation (equal); Visualization (equal); Writing – original draft (equal); Writing – review and editing (equal). **Steven R Craig:** Formal analysis (equal); Methodology (equal); Software (equal); Visualization (equal); Writing – review and editing (equal). **Chengzhi Shi:** Conceptualization (equal); Methodology (supporting); Visualization (supporting); Writing – original draft (supporting); Project administration (equal); Resources (equal); Supervision (equal). **V. Roshan Joseph:** Conceptualization (supporting); Formal analysis (supporting); Investigation (supporting); Project administration (equal); Resources (equal); Supervision (equal); Visualization (supporting); Writing – original draft (supporting); Writing – review and editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.