Sound localization through multi-scattering and gradient-based optimization

A gradient-based optimization (GBO) method is presented for acoustic lens design and sound localization. GBO uses a semi-analytical optimization combined with the principle of acoustic reciprocity. The idea differs from earlier inverse designs that use topology optimization tools and generic algorithms. We first derive a formula for the gradients of the pressure at the focal point with respect to positions of a set of cylindrical scatterers. The analytic form of the gradients enhances modeling capability when combined with optimization algorithms and parallel computing. The GBO algorithm maximizes the sound amplification at the focal point and enhances the sound localization by evaluating pressure derivatives with respect to the cylinder positions and then perturbatively optimizing the position of each cylinder in the lens while incorporating multiple scattering between the cylindrical scatterers. The results of the GBO of the uniand multi-directional broadband acoustic lens designs are presented including several performance measures for the frequency dependence and the incidence angle. A multi-directional broadband acoustic lens is designed to localize the sound and to focus acoustic incident waves received from multiple directions onto a predetermined localization region or focal point. The method is illustrated for configurations of sound hard and sound soft cylinders as well as clusters of elastic thin shells in water.


Introduction
Metamaterials are engineered materials that can control and guide the energy flow with capabilities exceeding those possible in conventional materials, enabling the control of acoustic, electromagnetic, and mechanical waves. Despite the exciting results obtained by the physics-inspired metamaterial designs, the forward design methodology relies on tuning parameters by trial and error where the low efficiency and limited exploration of the design variations lean to miss out the optimal solution. Using the inverse design approaches, one can start from the opposite end and optimize certain objective functions describing the desired functionality. The progress in inverse design has led to the discovery of metastructures with exceptional performances. Among various methods, large-scale gradient-based optimization (GBO) [1,2] has been a popular approach to design structures containing a vast number of degrees of freedom. The adjoint method [3] was employed in optimization of elastic and acoustic metamaterials [4,5], and of nanophotonic devices [6,7]. Topology optimization was utilized to design metasurfaces [8][9][10] and acoustic cloak composed of bimode on a honeycomb lattice [11], to improve the vibro-acoustic response of 3D sandwich panels [12] and predict structural parameters of 3D lattice structures [13,14], bimode materials [15], and mechanical metamaterials [16]. Popa and Cummer [17] proposed an optimization approach which produced simplified, broadband omnidirectional designs of transformation optics devices. A level set-based topology optimization using the FEM was employed in inverse design of optical hyperlens by Otomori et al. [18].
Placement optimization is commonly used to arrange scatterers or sources in acoustic/elastodynamic/thermal fields using genetic algorithms [19]. Håkansson et al. used multiple scattering theory and genetic algorithms in the inverse design of photonic crystals [20] and flat acoustic lenses [21,22]. Lu et al. [23] proposed an acoustic cloak design based on Bézier scatterers using genetic algorithms (GA) and simulated annealing. These GA used for optimization are gradient-free. The number of iterations needed before the solution converges tend to be large as the complexity of the model grows and as the number of design parameters increases. Stochastic algorithms do not require the evaluation of the objective function, but, if available, the analytical form of the gradients can enhance the quality and precision of the global optimal solution [24,25]. The function evaluations required in non-gradient optimization algorithms are orders of magnitude more than that is required for GBO methods [26]. In comparison to gradient-free based algorithms, GBO utilizes additional derivative information with respect to the design parameters that can result in a reduced number of iterations. The accuracy of solutions obtained with GBO algorithms [27] exceeds those obtained with stochastic algorithms. Erikson [19] performed a thermal design optimization with design parameters based on geometric properties such as length, width, angle, or position of the heat sources by solving PDE constrained optimization problem using fmincon solver and sequential quadratic programming (SQP) algorithms. In this model, Erikson [19] showed that the final optimal solution depended on the boundary conditions and the profile of the ingoing velocity field applied, and considered a loop over several random initial positions to find the best optimal solution which varied as the convective field and the ambient temperature changed. The gradient-based approach was developed by Amirkulova and Norris to design an acoustic cloak [28][29][30] and lens [31,32]. Andersen et al. [27] designed a multi-directional acoustic cloak by means of the boundary element method and shape optimization using the gradient-based SQP algorithm with the fmincon solver [33]. Chen et al. [34] discussed level set and eigenfunction optimization methods to represent the topography of a dielectric media and discussed gradient-based methods to solve different material design problems. McCormick and Shepherd [35] provided a comparison of performance gradient-based, fmincon [33], and evolutionary, Borg, algorithms to minimize the vibration response within a certain region of the beam at discrete frequencies, with constraints on total mass and static compliance. However, in their simulations [35] with fmincon, the analytical formulas for the gradients were not provided and the fmincon solver computed the gradients numerically using finite differences which requires additional runtime to compute the gradients. As these numerical gradients are not exact, fmincon produces poorer results as opposed to the case when gradients are evaluated analytically and supplied to fmincon.
In this paper, we explore the directionality and frequency band efficiency associated with focusing and localization effects based on discrete scatterers and show some ways to overcome limitations. We obtain a semi-analytical formula for the gradient of the pressure at the focal point with respect to positions of a set of cylindrical scatterers using two approaches namely forward and reciprocal formulations. We combine these formulae with GBO algorithms and apply them to design a uni-and multi-incidence-directional broadband lens by means of multiple scattering theory and applying the principle of reciprocity. A broad review of the literature on single and multiple scattering and of the concepts of multiple scattering is given in [36]; a survey of more recent findings on multiple scattering from obstacles in acoustic and elastic media is provided in [37]. In our model, an acoustic lens consists of a meta-cluster of rigid cylinders, cylindrical voids, or a set of elastic shells submerged in water. We design metastructures made of conventional simple isotropic materials that are available in nature and obtain useful and interesting effects. The idea differs from earlier inverse designs that use topology optimization tools and generic algorithms. The GBO algorithm maximizes the sound amplification at the focal point by evaluating pressure derivative with respect to the cylinder positions and then perturbatively optimizing the position of each cylinder in the lens while taking into account acoustic the multiple scattering between the cylinders. Computations are performed on MATLAB using advanced parallel optimization algorithms and a MultiStart optimization solver, and supplying the gradient of absolute pressure amplitude at the focal point.
The paper is organized as follows. Section 2 starts with a definition of multiple scattering problem, the position-dependent absolute acoustic pressure at the focal point, and the principle of acoustic reciprocity for this problem. Gradients of the absolute acoustic pressure at the focal point with respect to cylinder positions are then obtained by means of multiple scattering theory for forward and reciprocal formulations. Some properties of the gradient vectors are illustrated through numerical examples in Section 2.6. Application of the closed form for the gradient of absolute acoustic pressure at the focal point to acoustic lens designs are given in Section 3. Single and broadband as well as uni-and multi-directional focusing effects are illustrated using multiple reconfigurable cylinders as the lensing and localization mechanisms.

Problem Definition
We consider multiple acoustic scattering in two dimensions. Starting from the linear wave equation for an acoustic medium with sound speed c, and assuming time harmonic dependence e − i ωt of frequency ω results in the the Helmholtz equation for the acoustic pressure p(x), where k = ω/c is the wavenumber, c is the acoustic speed, and x is the position vector of point P with respect to origin O (see Figure 1). The factor e − i ωt is omitted in the following. The total field p is the sum of incident p inc and scattered p sc pressure fields: The incident field is the plane wave of unit amplitude in direction ψ Our objective is to maximize the acoustic pressure at the predefined focal point x f by rearranging the positions of M scatterers denoted by the set of vectors {r m } as shown in Figure 1. For simplicity, we take these to be circularly cylindrical scatterers, which may be either rigid cylinders or thin elastic shells illustrated in Figure 1. We define p f or equivalently p f ({r m }) to be the total acoustic pressure at the focal point for a given set {r m }, i.e., As a measure of the scattering we use the sound amplification, SA f , at the focal point defined as where the absolute value is |p| = (pp * ) 1/2 , and * denotes the complex conjugate. The objective is to maximize SA f by rearranging the cylinders to find an optimal set {r m }. The critical quantity that we use in the process is the gradient of the acoustic pressure at the focal point, defined next.

Position-Dependent Acoustic Pressure Field at the Focal Point
We employ multiple scattering theory to solve for the pressure at the focal point. Each scatterer, denoted by index m, radiates as a localized source, with multipole strengths B (m) n where n denotes the dependence on angle, e i nθ . Thus, according to (A8), the scattered field in the neighborhood of cylinder S m is given by where B (m) n are the unknown coefficients, and x m is the position vector of point P with respect to the cylinder center at O m : The functions V ± n (x), defined in Equation (A9), are where H (1) n is the Hankel function of the first kind of order n, arg x denotes the argument of vector x, i.e., the angle that the vector x makes with the positive x-axis, and it is measured as arg x ∈ [0, 2π) and arg(−x) = (arg(x) ± π) mod 2π [38]. Introducing Equations (3) and (6) into Equation (2), and evaluating it at the focal point P f , defines the total pressure at the fixed focal point: In order to apply the solution Equation (9) in practice, the infinite sum must be truncate. Let N be the truncation value of the infinite sum in Equation (9) chosen so that the sum converges. In practice, the value of N depends upon frequency and typically taking the truncation number N ≈ 2.5 ka is adequate [29]. Let us introduce the vectors b, v ∈ C M×(2N+1) . The components of the scattering coefficient vector n = {B n (r j )} found from Equations (A21) and (A22), and written in column vector form: (2) . . .
The elements of the dual vector v = {V (j) where is the position of focal point with respect to local coordinate at O m and r m is the position of scatterer with respect to origin at O. Then, the total pressure field at the focal point can be written in vector form as The Scattering Strengths The scattering problem is solved once we have the scattering strengths, represented by the column vector b which solves the linear system of equations where (2) . . .
The elements of the column vector a ∈ C M×(2N+1) are defined such that a = {A i n e i kx m , for a normal incidence ψ = 0, where x m = e 1 · r m , i n e ik(x m cos ψ+y m sin ψ) e − i nψ , for a general incident angle ψ, for n ∈ [−N, N] and N ∈ Z. Here, the matrix X is the interaction matrix that defines the coupling between each scatterer of the configuration (see Appendix A.1 for details) The matrix T (j) is the transition or T-matrix for scatterer j, and P j,m = P j,m ql is a Toeplitz matrix that depends on the position vector r jm depicted in Figure 1. The matrix P j,m takes into account the interaction between the scatterers, whereas the transition matrix T (j) depends on the shape and the physical properties of the material of cylinder, as well as the boundary conditions on the interfaces. The T-matrix is diagonal for scatterers with rotational symmetry. In general, for an obstacle with no rotational symmetry, e.g., an elastic thin cylindrical shell with internal spring-mass attachments [39], the T-matrix is nondiagonal. The off-diagonal terms couple different azimuthal orders of the incident and scattered field. At low frequency, the T-matrix may be truncated as a 3 × 3 matrix, in which case off-diagonal elements can be interpreted in terms of Willis coupling [40]. In this work, we consider 2-dimensional configurations of circularly cylindrical scatterers, for which the T-matrices become diagonal, see in [37] for specific details. In particular, P j,m = P j,m ql for j, m ∈ (1, M), j = m, l, q ∈ (−N, N), P j,m ql = P ql (r jm ) = H (1) l−q (kr jm )e i(l−q) arg r jm (18) where r jm = r j − r m , arg x ∈ [0, 2π). Alternatively, the acoustic pressure field at fixed focal point can be determined in a simpler form using the principle of acoustic reciprocity which is described next.

Acoustic Reciprocity
Here, we will illustrate the use of acoustic reciprocity [41] to define the pressure at the focal point P f due to a source located in a far-field. The key result is that the acoustic pressure at position x i due to a monopole source at position x f is identical to the pressure at x f resulting from a monopole of the same strength located at x i . This reciprocity enables us to relate the response for an incident plane wave in terms of the far-field Green's function.
We apply the reciprocity result to this problem by first considering two points x f and x i illustrated in Figure 2. In the forward problem, FP, the source is at x i and the pressure is computed at x f : where p(x; y) is the total acoustic pressure at x for a point source at y, which corresponds to a non-zero delta term on the right hand side of Equation (1). In the reciprocal problem, RP, the source is at x f and the pressure is computed at x i : By reciprocity we have We take x i in the far-field (see Figure 2), i.e., such that k|x i | 1, and choose the monopole amplitude such that the source at x i yields a plane wave of unit amplitude in the vicinity of x = 0, which is the incident wave for FP. Therefore, the source, or incident wave, for RP is Note that in this case in the local coordinates of multipole centered at O m : where x m and r f m are correspondingly the position vectors of an arbitrary point P and a source S with respect to O m . Thus, using the Graf's addition theorem (A10) in the neighborhood of cylinder S m (see Figure 1), we have for m = 1, M, n ∈ Z, where U ± n and V ± n are defined by (A11). Ultimately we consider |x m | < |r f m | to compute pressure and apply boundary conditions at the interface |x m | = a m . Thus, for |x m | < |r f m |, near cylinder S m we have where D We now consider the reciprocal problem by examining a monopole source located at position x f at the focal point P f = (r = r f , θ = 0). Then, according to Equation (A23), the far-field radiated by source at P f is where the far-field amplitude function in the direction of x i is The scattering coefficientsB are components of the dual vectorb that satisfies the multiple scattering condition where X is the interaction matrix defined by Equation (17). The elements of the dual vector (1) d (2) . . .
The far-field amplitude function (28) can then be expressed as and where γ = 2 i kπ|x i | e i k|x i | . By reciprocity principle, Equation (21), the resultant total pressure at the fixed focal point P f at θ = 0 by a plane wave incidence from θ = π can be written as This is to be compared with the solution (13) obtained by solving the direct multiple scattering problem. Of course, by reciprocity both are equal.

The Gradient Vectors q j and p j
The real valued vector q j is defined as gradient of the absolute value of total pressure field |p f ({r m })| with respect to positions r j . It can can be evaluated as where the associated complex valued gradient vectors are

The Pressure Gradient by Solving the Forward Problem
We now determine the explicit form of the gradient vector using the direct and reciprocal formulations, from Equations (13) and (33), respectively.
The pressure at the focal point follows from Equations (13) and (14) as The complex valued gradient vectors p j of Equation (35) can be evaluated by noting In component form, the last Equation (37) can be written as where Define the derivative function V ± n (x) as Then based on the fact that the T-matrices are independent of position, the gradient of the components of the matrix X are [29] where O nl are components of the zero matrix. The gradients in Equation (41) have the form [29] ∂P nl ∂r j The gradients of the components of v with respect to the positions r j of the j-th scatterer are where The gradient of the components of a with respect to the position r j of the j-th scatterer is where for a plane wave incidence in e 1 direction: and, hence, The plane wave incidence in e ψ direction is given by Equation (A4):

The Pressure Gradient by Solving the Reciprocal Problem
We now determine the explicit form of the gradient vector by solving the reciprocal problem. The pressure gradient at the focal point follows from Equations (29) and (33) as where theˆindicates the use of the principle of reciprocity. Of course, p j of (48) should equal p j , and their equivalence is verified numerically below. In component form, (48) can be written as where l } satisfies Equation (29). Other gradient components and variables are defined the same as in Section 2.4.

Numerical Examples
Before considering a sound localization problem, numerical examples will be presented for the gradient vector q j which has some interesting characteristics. Computations are performed on MATLAB for configurations of uniform rigid cylinders of radius a. We consider two problem formulations: forward, FR, and reciprocal, RP, and make comparison of simulation times for both formulations. The reciprocal formulation of inverse design of acoustic lens takes lesser time to compute the absolute pressure and it gradients at the focal point.

Illustration of Principle of Reciprocity
In order to illustrate the application of reciprocity in acoustic lens design, we present numerical examples using two configurations. Figure 3a (top) shows the first configuration of a vertical setup of three cylinders of radius a and of distance 2.2a between each two adjacent cylinders. Figure 3 illustrates the variation of absolute pressure p f and phase angle θ f at a focal point, and the gradient vector components q xj q xj , q yj and q yj as functions of normalized wavenumber ka for FR and RP formulations. The absolute pressure at a focal point that is at the point (10a, 0) was computed at values of normalized wavenumber ka ∈ (0, 15] using both the FP and RP methods and they coincide perfectly as shown in Figure 3a. Figure 3b presents θ f , the phase angle of the total pressure at the focal point versus the non-dimensional wavenumber ka. The gradients in the x direction are shown in Figure 3c and were computed using the FP method. The corresponding gradients that were computed using the RP method are exactly the same and they are not shown just for the clarity of the figure. The gradients in the y direction are shown in Figure 3d. They were computed using both the FR and the RP methods and they perfectly overlap. Due to the symmetry of the given configuration with respect to the x-axis, the gradients of cylinders 1 and 3 in the x direction coincide while in the y direction they have an opposite sign. Figure 4 shows the non-symmetric configuration of two rigid cylinders each of radius a, where the first is centered at the origin and the second is centered at the point (1.5a, 3a). The figure illustrates the dependency of absolute pressure p f and phase angle θ f at a focal point (10a, 0) and the gradient vector components q xj q xj , q yj and q yj as functions of normalized wavenumber ka for FR and RP formulations for a non-symmetric set of two rigid cylinders. This is given to emphasize that both FP and RP methods give the same results for the nonsymmetric configurations. The results of computing the absolute pressure p f , the phase angle θ f , and the absolute pressure gradients in the x and y directions are the same for both the FR and the RP methods. The same symmetric and non-symmetric configurations were used to illustrate the non-linearity of the problem in hand. The absolute pressure was computed for the values of non-dimensional wavenumber ka ∈ (0, 15] and for the incident angle ψ ∈ [−π/4, π/4]. Figure 5c,f show the absolute pressure versus both the non-dimensional wavenumber and the incident angle. The 3D mesh plots illustrate the existence of a lot of local maxima and minima, so a small deviation in the starting optimization point and/or the incident angle and the non-dimensional wavenumber will lead to a different optimized configuration. Figure 5a,d shows the absolute pressure versus the non-dimensional wavenumber and they can be seen as the cross sectional views of the 3D plots at normal incidence: ψ = 0. Figure 5b,e shows the absolute pressure versus the incident angle and can be seen as the cross-sectional views at ka = 7.5.

CPU Timing
In order to compare between the CPU time of the two computing methods, two configurations were used as depicted in Figure 6a,b. Figure 6 shows the measured time. The elapsed times were used in computing the absolute pressure only, and the absolute pressure and the gradients sequentially varying the wavenumber ka and using both the FP and the RP methods were measured. This was done at different discrete values of ka ∈ (0, 15]. It can be seen from Figure 6 that the computation time generally increases for all computations as the non-dimensional frequency increases. Furthermore, generally the computation time in the configuration with higher number of cylinders is higher. This is due to the increase in the rank of the matrices to be evaluated during the computation process. In addition, the computation time increases for both FP and RP formulations when gradients are evaluated. However, providing gradients improves the performance of GBO algorithms allowing the faster convergence and more accurate and more reliable computations [33]. For 1 × 3 configuration, FP with gradients takes the most compute time and FP without gradients is the least computation time; the compute times of RP with gradients and FP without gradients are of similar order. For 2 × 3 configuration, the computation times of different methods can be put in a descending order as follows: FP with gradients, RP with gradients, FP without gradients, and RP without gradients. This depicts that using the RP method is computationally less expensive than the FP method.

Sound Localization by Maximizing the Absolute Pressure at the Focal Point
The sound localization was implemented by maximizing the absolute pressure at the focal point and by employing a gradient based optimization. Computations are performed on MATLAB using advanced parallel optimization algorithms, the MultiStart optimization solver combined with fmincon solver and SQP algorithms, and supplying the gradients of absolute pressure amplitude at the focal point |p f |. Figure 7 illustrates initial random configurations of cylinders used in simulation where the locations of the control cylinder positions are denoted by pink dots. For given scatterer location points, the gradients q j are in the direction of greatest increase of absolute pressure amplitude at the focal point |p f |. As seen in the example of Figures 3 and 4, this provides the optimal directions to increase |p f | by incremental displacement. Thus, in order to maximize the |p f | produced by acoustic lens, we will move the cylinders inside the circular region in the direction of the gradient vectors, i.e., q j using gradients evaluated by FP and RP formulations. This is achieved in MATLAB with the use of Global Optimization and Parallel Computing Toolboxes, and by supplying the derived analytic form of the gradients vectors from Equation (34).

Acoustic Lens Design
In order to study the focusing and localization efficiency, we consider three cases: In the first case, we maximize the absolute pressure at a fixed value of frequency and normal incidence while supplying the gradients q j . In the second case, we maximize the root mean square of a set of the absolute pressure over a range of frequencies. In the third case, we maximize the root mean square of a set of the absolute pressure over a range of frequencies and a set of the incident angles of the impinging plane wave. In each case, we supply the analytical formulas for the gradients of objective functions. The cost functions are nonconvex with many local maxima. The global maximum is unbounded in each considered cases; we solve the non-convex optimization problems with nonlinear constraints. For simplicity and demonstration of implementation of the proposed approach, we consider a configuration of uniform cylinders/shells of radius a.

Constrained Optimization Problem
The idea is to maximize the absolute pressure at the focal point, |p f |, by rearranging the position of each cylinders and starting with random configuration of scatterers as depicted in Figure 7. The design vector, x M , consist of all scatterer position and can be represented in a column vector form as where M is the total number of scatterers used in lens design. The nonlinear constraints for the design vector components are as follows: 1. The cylindrical scatterers are constrained to move inside a fixed circular region with the radius r = R out , see Figure 7: where r j = (x j , y j ) and j = 1, 2, ..., M.
2. In order to avoid overlapping the distance between the centers of cylinders/shells are constrained by [29]: As was shown in our earlier work [29] on acoustic cloak design, the closest proximity according to the constraint (53) should be far greater than the viscous skin depth δ v = (2ν/ω) 1/2 [42]. In water, at temperature T = 5 • C with kinematic viscosity ν = 1.5182 mm 2 /s [43], this is δ v = 0.000695m/ √ s √ f for frequency f (Hz); therefore, the undamped acoustic model is accurate for frequencies in kHz and above, which is the range of practical interest. We use direct optimization methods such as SQP algorithms to get the best performance by taking advantage of the analytical form for the gradients of objective functions combined with parallel computing. It is suggested in [33] to include a gradient evaluation in the objective function for faster or more reliable computations. In our computations, we used MATLAB built fmincon solver which is a nonlinear programming solver that finds a local minimum of constrained nonlinear multivariable function. Therefore, to maximize the absolute pressure, we changed the sign of objective function and its gradients.

Broadband Multi-Directional Lens Design
The procedure for a broadband multi-directional GBO is as follows. First, we define the cost function as the root mean square (RMS) of a set of absolute pressure at the focal point |p f | s over some range of normalized wavenumbers k i a (i = 1, 2, ..., N k ) and over a range of incident angles ψ l (l = 1, 2, ..., N l ): Then, we define the broadband gradient vectors with respect to positions r j : which can be found in terms of the individual single frequency gradients as where |p f | RMS is defined by (54) and q j (k i a) are defined by (34) and evaluated at normalized wavenumbers k i a (i = 1, 2, ..., N k ) and angles of incidence ψ l (l = 1, 2, ..., N l ). For a single angle of incidence ψ 1 , Equations (54) and (56) reduce to the form and Both FP and RP formulations allow designing broadband unidirectional lens. Note that a formulation of RP reciprocal problem is given at normal incidence. Therefore, a reciprocal formulation can only be used to design broadband unidirectional lens, whereas a forward formulation can be used in all the cases: a single frequency at single incidence, i.e., unidirectional lens optimized at single frequency; broadband at single incidence, i.e., broadband unidirectional lens; and broadband multi-directional lens designs.

Numerical Examples of Sound Localization
In this section, the results of the gradient-based optimization of the uni-and multidirectional acoustic lens designs are presented including several performance measures for the frequency and the incidence angle dependencies. The optimized model is an improvement over the initial random configuration of cylinders. Numerical results are demonstrated for configurations of rigid cylinders, cylindrical voids, and empty thin elastic nickel cylindrical shells of thickness h = 0.1a. We consider cylinders and shells submerged in a medium with the acoustic properties of water: ρ 0 = 1000 kg/m 3 , c 0 = 1480 m/s. All computations are performed for an incident plane wave propagating from left to right. We illustrate the effect of acoustic lens device on plane wave incidence localized at focal point which resembles the focusing effect of acoustic Luneburg lens.
Various values are taken for the non-dimensional wavenumber, ka; the total number of scatterers in the focusing lens, M; the number of discrete frequency points, N k ; and the number of discrete incidence angle points N l for a broadband localization and focusing as defined in Equation (54). Greater accuracy is observed, as expected, with increased number of scatterers M. However, large values of ka and M require longer computation times, and some numerical experimentation is necessary to find the smallest values for which the pressure field can be localized to the desired degree in the acoustic lens.
Numerical computations are performed on MATLAB using parallel optimization algorithms with a MultiStart optimization solver, and supplying the gradients of absolute pressure |p f | with respect to position vectors. We start with an initial random planar configuration of cylinders or thin shells of uniform size and radius a = 0.0075 m, e.g., as shown in Figure 7a-e, and consider a focal point located at x f = (0.0865, 0) m. We use MATLAB built function fmincon combined with SQP algorithm and As mentioned in fmincon documentation [33]: "a solver can reach a point x such that x is feasible, but finite differences around x always lead to an infeasible point. In this case, a solver can fail or halt prematurely. Providing a gradient allows a solver to proceed." Therefore, supplying the gradients accelerates the optimization process by searching a wide variety of start points in parallel and allows the solver to reach the optimum values which may otherwise be impossible [29]. Here, we present results for a ring configuration of M cylinders running 200 different scenarios (initial random configurations) for a single frequency optimization and running 100 scenarios for broadband and multi-incidence-angle optimizations concurrently on 48 CPUs of Dell workstation, described next.
We evaluated the absolute pressure |p f | rel relative to the incident background pressure at the focal point and the magnification factor (MF) defined as where |p inc (x f )| is the incident plane wave of unit amplitude at the focal point, |p f | and |p f | initial are correspondingly the absolute pressure at the focal point for the optimized final and initial configurations.

Lens Design at Single Frequency and Normal Incidence
First, we present the localization effect for a plane wave incident on configuration of rigid cylinders depicted in Figures 8 and 9 where various discrete values of normalized wavenumber, ka, and the total number of scatterers in acoustic lens, M, were considered. The optimized lens device in each example is a configuration of identical rigid cylinders located inside of circular region of radius R out = 0.0715 m. Figure 8 shows the real part of total acoustic pressure field p at normal plane wave incidence for a configuration of M = 50 rigid (or sound hard) cylinder at discrete values of wavenumbers ka = 1, 1.5, and 2. Figure 9 illustrates the absolute total acoustic pressure |p| for the optimized configurations of M = 13, 31, and 50 cylinders ka = 1, 1.5, and 2. The value of |p f | and magnification factor, MF, are depicted in Figure 9 for considered values of M and ka under each subfigure. Next, we investigate the sound localization effect for a plane wave incident on configuration of void cylinders considering discrete values of normalized wavenumber, ka, and the total number of scatterers in acoustic lens, M. Figure 10 illustrates the absolute total acoustic pressure |p| computed at normal plane wave incidence for configuration of M = 13, 31, and 50 cylindrical voids (or sound soft cylinders) at normalized wavenumbers ka = 1, 1.5, and 2. The absolute total acoustic pressure |p| is computed for the optimized configurations of M = 13, 31, and 50 cylinders ka = 1, 1.5, and 2. The value of |p f | and magnification factor, MF, are depicted in Figure 10  Finally, we present the sound localization and focusing mechanism for empty elastic thin cylindrical shells submerged in water. We consider empty thin elastic nickel cylindrical shells of thickness h = 0.1a with mechanical properties: ρ = 8850 kg/m 3 , c p = 5480 m/s where a is the outer radius of scatterer. Figures 11 and 12 illustrate the localization effect for a normal plane wave incident on configurations of elastic thin cylindrical shells at various discrete values of normalized wavenumber, ka, and the total number of scatterers M. Figure 11 displays the distribution of absolute total pressure filed |p| for the optimized configurations of 50 nickel thin cylindrical shells (a) at ka = 0.75, (b) ka = 1.5, and (c) ka = 2. The elastic thin cylindrical shells can display large scattering at low frequency due to internal resonances of longitudinal and flexural waves. This is confirmed by results in Figure 11 Figure 11a. Figure 11 also shows that sound focusing becomes more noticeable with increase of frequencies. Figure  It is noticed that, in general, the sound focusing and localization effects are enhanced with the increase of frequency and number of scatterers. The total field distributions for void cylinders and elastic thin shells differ both qualitatively and quantitatively from the total field produced configuration of rigid cylinders. Both void cylinders and elastic thin shells produced high amplification of sound at focal point, i.e., MF values are high, compared to initial configuration. The design can be further improved by considering the cylinder radii as an additional design parameter which will require deriving analytical formulas for the gradients of absolute pressure with respect to cylinder radii and will be studied elsewhere.

Performance of Different Solvers for Lens Design at Single Frequency and Normal Incidence
In this section, we investigate the performance of GBO and GA, compare the best final optimal results produced, and measure the time taken by the optimization functions in all cases. We consider configurations of M = 13 rigid and elastic scatterers at discreet values of wavenumbers ka = 0.75, 1.5, and 2 for convenience. For GBO, we first run the fmincon solver with SQP algorithms and then we implement MultiStart optimization combined with fmincon and SQP algorithms. In GBO, we start each time with the feasible design vector x M 0 , i.e., the initial random configuration depicted in Figure 7c, and use the central difference for finite differences.
For GA, we solve the optimization problem by running ga solver from Global Optimization Toolbox on MATLAB. The ga solver starts with its own initial random population and we use "rng default" option for reproducibility, analysis and testing of numerical results. We analyze the results for running the ga solver by evaluating fitness functions in parallel using 48 workers (CPUs) concurrently. We keep the maximum number of iterations before the algorithm halts at default option [44], i.e., MaxGenerations = 100 * 2M = 2600 and vary MaxStallGenerations option to obtain finer solutions. The MathWorks documentation [44] indicates that the MaxGenerations option determines the maximum number of generations that the ga takes and MaxStallGenerations option controls the number of steps the ga looks over to see if it is making progress and can enable the solver to obtain a better solution by allowing more function evaluations. We observed that running the ga solver with above mentioned options, the nonlinear constraints (53) were not always satisfied, the (52) constraints were slightly violated in each of considered cases.
Numerical simulations presented in Table 1 through Table 6 are performed on MAT-LAB using Dell Workstation with 48 CPUs. In each case, the first column represents the best values of absolute pressure |p f | obtained by a solver and the second column shows the corresponding elapsed time to compute the optimal values. The elapsed compute time is given in seconds. Table 1 through Table 3 compares the results of fmincon, MultiStart, and ga solver runs for configurations of M = 13 rigid cylinders. Table 1 depicts the results produced by fmincon with three different options, namely, running fmincon without (w/o) gradients sequentially and in parallel using concurrently 48 workers (CPUs), and running fmincon with gradients sequentially. The results presented show that running fmincon without (w/o) gradients sequentially take 5 to 16 times longer time compared to running fmincon with gradients sequentially. Parallel computing allows to reduce compute time approximately to 3.5 to 7 times. Produced optimal values vary with change of wavenumber. For ka = 0.75 fmincon with gradients produced better results and for ka = 1.5 and 2 fmincon without gradients yield higher optimal values. Table 1. The comparison of results for configuration of M = 13 rigid cylinders at ka = 0.75, 1.5, and 2 using fmincon solver with three different options: running fmincon without (w/o) gradients sequentially and in parallel using concurrently 48 workers, and running fmincon with gradients sequentially. The total elapsed compute time is given in seconds. When analytically evaluated gradients are supplied to fmincon, we can run MultiStart in parallel considering different scenarios concurrently. Table 2 illustrates results for GBO using MultiStart combined with fmincon and SQP algorithms. We run 20, 48, and 200 parallel scenarios with MultiStart using 48 workers concurrently, and providing gradients to fmincon which runs sequentially for each scenario. Considering 20 scenarios, MultiStart already outperforms single fmincon run results shown in Table 1 producing higher optimal values in similar compute times. Increasing the number of scenarios improves the optimal value produced by MultiStart.  Table 3 shows the results for running ga solver in parallel using concurrently 48 workers with three different options: MaxStallGenerations = 50, 500, and 1300. The ga produced poorer results and took longer to find the best results. The ga solver violated nonlinear constraints (52), although it is still acceptable as predicted optimal scatterer positions are not more one cylinder radii away from constrained circular region R out . However, for ka = 0.75, running ga solver with MaxStallGenerations = 50 resulted in infeasible solution violating both (52) and (53) constraints and allowing cylinders to overlap. Figure 13 depicts the final scatterer configurations obtained by ga solver for ka = 0.75 with three different options: MaxStallGenerations = 50 (a), 500 (b), and 1300 (c). For MaxStallGenerations = 500, both constraints (52) and (53) are violated, and for MaxStallGenerations = 50 and 1300 the constraints (52) are slightly violated. The results show that ga struggles to solve the nonconvex optimization problem with nonlinear constraints. Table 3. The comparison of results for configuration of M = 13 rigid cylinders at ka = 0.75, 1.5, 2 running ga solver in parallel using concurrently 48 workers with three different options: MaxStallGenerations = 50, 500, and 1300. In each case, the constraints (52) are slightly violated and in some cases predicted final cylinder positions are up to one cylinder radii away from a circular region of radius R out . At ka = 0.75, ga with option MaxStallGenerations = 500 produces an infeasible solution violating both (52) and (53) constraints. Tables 4-6 analyze the results of fmincon, MultiStart, and ga solver runs for configurations of M = 13 thin elastic shells with the same options used for the configurations of rigid scatterers presented above in Table 1 through Table 3. Analyzing parallel runs of fmincon w/o gradients and MultiStart runs starting with 20 scenarios indicate that MultiStart produced much more enhanced results within a comparable compute time. In addition, the MultiStart solver starting from 200 various random configurations yields the best results and outperforms fmincon solver with and without providing gradients as well as the ga solver. The ga solver running with parallel option performs much poorer producing much lower values of |p f | and violating constraints (52). Increasing MaxStallGenerations option values with ga solver overall has tendency to obtain more improved results but also required longer compute times from 3 to 6 times. The GBO results illustrated in this section show the advantages of using MultiStart with fmincon which produced the highest values of absolute pressure |p f | when combined with parallel computing. The difference becomes more noticeable with the increase of number of scatterers and wavenumber. We observed that MultiStart integrated with fmincon supplying gradients was able to converge and obtain optimal solutions for configurations with larger number of scatterers and high frequencies for which single fmincon runs without supplying gradients were not able to converge and/or to obtain optimal values in reasonable time. Table 4. The comparison of results for configuration of M = 13 thin elastic shells at ka = 0.75, 1.5, 2 using fmincon solver with three different options: running fmincon without (w/o) gradients sequentially and in parallel using concurrently 48 workers, and running fmincon with gradients sequentially.   | on ka and ψ simultaneously via 3D parametric surface plots. Here, the lower blue mesh denoted by "1" represents the results for a non-optimized configuration and the higher mesh indicated by "2" exhibits the results for the configuration optimized. The results in Figure 14b-d for the non-optimized configurations "1" show overall low values of |p f | RMS in comparison to the mesh surface plots for optimized configurations "2" as illustrated; the results for "2" show clearly the amplification of |p RMS f |. For the fixed value ka = 0.45, the surface plots in Figure 14d show that the optimized configuration "2" produces the peak performance ka = 0.45 and near the incidence angles −π/4 and π/4; the surface plot "2" produced lower values of |p RMS f | at wavenumbers near ka = 0.55 for all ψ values considered which also coincides with results in Figure 14a, denoted by the red curve "1". The trade-off between the increase of the frequency band and focusing efficiency leads to shift of peak performance for broadband focusing achieving the maximum performance at ka = 0.55 which is not intuitive. Figure 15 illustrates the performance of broadband unidirectional lens at different values wavenumbers ka and displays the distribution of absolute total acoustic pressure field |p| for the optimized unidirectional lens, the configuration of M = 50 rigid cylindrical scatterers. The top row of Figure 15 shows results for the absolute pressure |p| for values of wavenumber ka at which the lens was optimized, i.e., (a) ka = 0.35, (b) ka = 0.45, and (c) ka = 0.55. The bottom row of Figure 15 displays the absolute total field for values of ka at which the lens was not optimized, i.e., (d) ka = 0.225, (e) ka = 0.4, and (f) ka = 0.585. Figure 15a-c clearly shows the focusing effect. Figure 15d,f illustrates that focusing amplitudes reduces, and beyond these values, i.e., 0 < ka ≤ 0.225 and ka > 0.585, the focusing effect diminishes. Importantly, Figure 15e exhibits clearly the sound localization at the desired focal region at wavenumber ka = 0.4 at which the configuration was not optimized.  Figure 16 displays the distribution of absolute total acoustic pressure field |p| for the optimized unidirectional lens, the configuration of M = 50 rigid cylindrical scatterers at ka = 0.55 varying values of incident angle ψ. The unidirectional lens is optimized for ka = 0.35, 0.45, and 0.55 at normal incidence, i.e., ψ, and the performance is shown at angles of incidence ψ = π/18, π/8, π/4, π/18, −π/8, −π/4. The unidirectional lens is able to localize the sound well at normal incidence and at small incident angles ψ = ±π/18 and with increase of ψ angle values focusing effect deteriorates. Figure 17 presents |p f | RMS , as a function of incidence angle and the normalized wavenumber ka considering a simultaneous multi-frequency and multi-angle optimization of |p f | RMS at a broad range of the incidence angle ψ ∈ [−π/4, π/4] and 3 values of wavenumber ka for configurations of M = 50 rigid cylinders. Figure 17a-c displays the variation of |p f | RMS with the angle of incidence at three fixed values of wavenumber ka = 0.35 (a), ka = 0.45 (b), and ka = 0.55 (c) for the optimized and initial unoptimized configurations cylinders. Figure 17d illustrates the 3D surface plot of |p RMS f | versus ka and ψ simultaneously. Here, the lower blue mesh denoted by "1" represents the results for a non-optimized configuration and the higher mesh indicated by "2" exhibits the results for the configuration optimized at three values of wavenumbers and five incident angles. Figure 17d shows that the non-optimized surface mesh "1" depicts the low values of |p RMS f | in all considered values of ka and ψ, and the surface mesh "2" is slightly lifted up with its highest amplification around ψ = −π/4 and ψ = π/4 at the largest values of ka considered, i.e., ka = 0.55. As can be seen from these figures, optimizing |p f | RMS for angle of incidence ψ ∈ [−π/4, π/4] using five equidistant values of ψ, the overall |p RMS f | show low values in the incidence range considered which suggests that this range is too ambitious for the small number of discrete incident angles and frequencies considered. Therefore, we narrow it down to ψ ∈ [−π/8, π/8].  Here, the red-dashed curve corresponds to configuration optimized at single value of ψ = 0 corresponding to normal incidence and the solid blue curve mean results for an unoptimized initial configuration at normal incidence. Figure 18a-c shows that adding more angles of incidence as a parameter in optimization broadens the range of application of multi-directional acoustic lens and improves its performance. Although optimizing at three values of ψ = −π/8, 0, π/8] shifted up |p f | RMS at edge values ψ = −π/8 and ψ = π/8], it reduced the peak at ψ = 0. The peak performance at ψ = 0 was further enhanced by optimizing |p f | RMS at five equidistant values of ψ ∈ [−π/8, π/8] across all three values of ka as can be seen from these figures. Figure 18d illustrates the dependency of the objective function |p RMS f | on ka and ψ simultaneously via 3D parametric surface plots. Here, the lower blue mesh denoted by "1" represents the results for a non-optimized configuration and the higher mesh indicated by "2" exhibits the results for the configuration optimized at three values of wavenumber ka and at five equidistant values of incident angles ψ ∈ [−π/8, π/8]. In comparison to Figure 17d, the results in Figure 18 show better performance of GBO algorithm to produce a broadband multi-incident sound localization effect. The optimized design illustrates broadband performance as compared to initial configuration response although it is not included as an optimization criterion. Figure 16. The absolute value of the total acoustic pressure field |p| for the unidirectional lens optimized at normal incidence, i.e., the configuration of M = 50 rigid cylinders at wavenumber ka = 0.55 and different angles of incidence ψ. The performance of optimized unidirectional lens at normal incidence, i.e., ψ = 0, is depicted in Figure 15c. We notice that the trade-off between the increase of range of incidence and focusing efficiency leads to reduced peak performance for simultaneous broad-range incidence and broadband focusing. Results in Figure 17 for broader incidence angle range ψ ∈ [−π/4, π/4] are poorer compared to range ψ ∈ [−π/8, π/8] given in Figure 18. The results for a broader range of ψ, i.e., ψ ∈ [−π/4, π/4], can be improved increasing the number of incidence angles in RMS computation, considering non-equidistant values of ψ, and choosing more points (angles) around region where |p f | RMS values are the lowest.

Conclusions
We have demonstrated that the GBO technique combined with the principle of acoustic reciprocity leads to an effective method to optimally arrange scatterers so as to produce acoustic focusing. The approach uses semi-analytic expressions for the gradient of the focal pressure with respect to the scatterer positions. Two problem formulations have been considered for the inverse design of the acoustic lens: forward and reciprocal. A comparison of simulation times for both formulations shows that the reciprocal formulation takes less time to compute the absolute pressure and its gradients at the focal point. The specific type of acoustic lens considered consists of a meta-cluster of rigid cylinders or cylindrical voids or thin elastic shells submerged in water. The computations were performed using advanced parallel optimization algorithms on MATLAB, a MultiStart optimization solver, and semi-analytic expressions for the gradient of pressure at the focal point. Numerical results for the uni-and multi-directional broadband acoustic lens designs indicate that significant sound pressure amplification can be achieved by optimal scatterer rearrangement using this analytical-numerical technique. The greatest magnification in amplitude is realized using void-like scatterers.
The presented gradient assisted inverse design model provides a means to establish realistic strategies for an acoustic lens device and for its practical underwater applications. The proposed method can be extended to design 3D lenses including air-born sound and underwater acoustic, elastodynamic, or electromagnetic wave localization and focusing effects. The gradient assisted inverse design of acoustic lens can be further enhanced and implemented by integrating the gradient information, i.e., Equation (37) with deep reinforcement learning algorithms [45] and generative networks [46] which have potential to search for the globally optimized devices over a broad range of parameters and can provide better solutions than ones produced by the-state-of-the-art optimization algorithms. Consider acoustic multiple scattering by M obstacles. For simplicity, the obstacles are assumed to be cylinders S m , (m = 1, M) centered at r m . A schematic configuration of cylindrical elastic shells is given in Figure A1. The incident field in the neighborhood of cylinder S m is where x m is a position vector of point P with respect to the centers of multipoles at O m (see Figure A1): the functions U ± n (x) are U ± n (x) = J n (k|x|)e ± i n arg x .
Here, arg x ∈ [0, 2π) and arg (− x) = arg x ± π mod 2π, and J n (x) is the Bessel function of the first kind of order n. In Equation (A1), the unknown coefficients A (m) n for a plane wave incidence are derived assuming no source term: p S = 0, and that the incident wave is the plane wave of unit amplitude in direction ψ: p inc = e ike ψ ·x ⇒ A (m) n = e ikê(ψ)·r r r m e in( π 2 −ψ) = e ik(x m cos ψ+y m sin ψ) e in( π 2 −ψ) .
A point source normalized by its amplitude at the origin is defined by and r = |x − x |. For a source at point S depicted in Figure A1, the coefficients A (m) n follow from Equation (A6), noting that in the local coordinates of multipole centered at O m : for m = 1, M, n ∈ Z, where the Graf's addition theorem (A10) is used for |x m | < |x m | (see Figure 1).
The total scattered field p sc , considered as a superposition of the scattered fields by all cylinders, may be expanded as a sum of multipoles: where p In order to use boundary conditions on the surface of each cylinder S m , we will express the total field using Graf's theorem ([47] Equation (9.1.79)): The functions U ± n (x) and V ± n (x) possess the properties U ± n (−x) = (−1) n U ± n (x), V ± n (−x) = (−1) n V ± n (x) = V ± −n (x) = (−1) n V ± n (x). (A11) Let r mj = r m − r j (= −r jm ) be a position vector of multipole O m with respect to multipole O j . As x = r m + x m = r j + x j → x m = x j + (r j − r m ), the total field p in the neighborhood of cylinder S j can be written as Then, noting the properties of V + n (x), Equation (A11), and using Graf's theorem, we obtain for |x j | < l j , where l j = min |r jm |: where P nl (x) ≡ V + l−n (x).
Here, the matrix P = [P nl ] is equal to the transpose of Martin's S = [S nl ] matrix [36], P = S T . The total incident field impinging on the cylinder S j is a sum of the last two terms on the right hand side of Equation (A13), i.e.,