© The Author(s), Published by EDP Sciences, 2026

1 Introduction

In many real-world acoustic scenarios, the challenge is to reduce the sound pressure level at a specific receiver location within complex propagation domains. In such cases, directly modifying the sound source is often impractical, and receiver positions are usually fixed by functional, operational, or other constraints. Therefore, adjusting the acoustic properties of the surfaces that define the propagation domain is a natural and effective strategy for influencing the sound field at locations where performance is ultimately assessed.

The importance of surface-based acoustic control is widely recognized across various fields of acoustics. In environmental acoustics, it is relevant to outdoor sound propagation near buildings and façades [1–3], as well as urban noise mitigation and architectural environments [4, 5]. In industrial acoustics, surface treatments are commonly used to control noise inside ducts, machinery enclosures, and HVAC systems [6]. Similar source-surface-receiver configurations occur in ultrasonic applications, where wave interaction with engineered boundaries governs reflection and transmission phenomena [7, 8], and in underwater acoustics, where structured or impedance-tailored boundaries are investigated for reflection control, absorption, or scattering in marine acoustic environments [9–11].

Despite the diversity of these applications, they share a common physical configuration: a source radiating in the presence of one or more bounding surfaces, with acoustic performance evaluated at a specific observation point. In such settings, the ability to tailor the acoustic response of the bounding surface is central to controlling the resulting sound field. This requirement has gradually shifted attention from empirical surface treatments to design-oriented approaches, in which the surface response is explicitly modeled and engineered.

Various approaches have been proposed to control sound propagation in semi-open environments by modifying the properties of bounding surfaces. These approaches primarily aim to alter surface-related contributions to the acoustic field, thereby influencing interference phenomena commonly referred to as ground effect in environmental acoustics. Among the most widely studied solutions are porous absorbent materials [12], surface roughness modifications [13], and resonant treatments such as quarter-wavelength resonators [14] and cavity-based resonators [15].

Unlike empirical resonator-based treatments, acoustic metamaterials rely on a design-oriented framework that enables predictive control of the acoustic response through physically grounded models [16]. Acoustic metamaterials are artificially engineered structures whose effective acoustic properties arise primarily from their geometry rather than their constituent materials [17]. Their behavior is typically described in terms of effective parameters, which explicitly account for resonant mechanisms and associated visco-thermal losses [18]. This modeling approach enables a direct link between geometric design parameters and the resulting acoustic response, allowing the behavior of the structure to be predicted prior to fabrication. Within this framework, resonant elements such as Helmholtz resonators (HRs) [19–21], membrane-based systems [22], micro-perforated panels (MPP) [23, 24], or sonic crystals [21, 25] are not employed as isolated components but as part of a unified design strategy in which losses, coupling effects, and effective impedance are explicitly taken into account. One of the most notable outcomes of this approach is the possibility of achieving perfect absorption, defined as the complete suppression of reflected sound energy at the surface. This condition corresponds to impedance matching, where the real part of the normalized surface impedance equals unity and the imaginary part is zero [26, 27].

A key feature of acoustic metamaterials is their ability to provide effective sound control at wavelengths much larger than the characteristic size of the unit cell, enabling compact and subwavelength surface treatments [16, 28]. When these concepts are applied to structures that are thin compared to the wavelength and integrated directly into a two-dimensional (2D) boundary, the resulting systems are commonly referred to as acoustic metasurfaces [29, 30]. Acoustic metasurfaces extend the metamaterial paradigm to surface-based configurations [29, 30].

Acoustic metasurfaces are defined by numerous geometric and physical design parameters, whose combined effects determine their effective acoustic response. Consequently, systematic optimization techniques have become a widely adopted tool in the field [31]. Most optimization strategies reported in the literature focus on maximizing a specific intrinsic acoustic quantity associated with the metamaterial itself. For example, these include the absorption coefficient α [19, 21, 32], the scattering coefficient δ [33, 34], the transmission loss TL of the metamaterial over the widest possible frequency bandwidth [35], or the size of a bandgap [36].

By design, optimization strategies based on intrinsic metamaterial properties rely solely on quantities evaluated independently of any specific source-receiver configuration. Consequently, the interaction between the metamaterial and the surrounding medium is not explicitly considered, and the optimization process remains disconnected from the global acoustic configuration in which performance is ultimately assessed at a specific point in space. This limitation becomes critical when the metric of interest depends on the combined effects of the metasurface, surrounding boundaries, and source-receiver arrangement, resulting in physical phenomena that cannot be captured by intrinsic properties alone.

Several recent studies have highlighted that optimizing intrinsic properties does not necessarily result in optimal acoustic performance when the full propagation configuration is considered. In particular, reductions in expected performance have been reported once complete geometries, diffraction effects, or source-environment interactions are taken into account [37–39].

In this work, we address this limitation by formulating the optimization problem directly in terms of acoustic performance evaluated within the full propagation environment. However, in many practical situations, metasurfaces are still designed by optimizing intrinsic quantities of the structure itself, such as the absorption coefficient, which is evaluated under idealized excitation conditions and independently of the propagation geometry. Specifically, we consider a semi-open acoustic configuration in which an acoustic metasurface is integrated into a bounding surface, with the objective of maximizing insertion loss (IL) at a given receiver position. Unlike conventional optimization strategies that target intrinsic metamaterial properties evaluated under idealized conditions, the proposed approach explicitly accounts for the interaction between the metasurface, the surrounding boundaries, and the source-receiver arrangement.

First, the overall methodology of the study is described in Section 2, outlining the objectives and modeling assumptions. A description of the metasurface, the finite element numerical model, and the optimization method used in this approach is provided. Next, in Section 3, the optimization results are presented and compared with conventional approaches for optimizing the absorption coefficient of the metasurface alone. In addition, the robustness of the method and the understanding of the physical phenomena are discussed in Section 4. Finally, Section 5 concludes this work.

2 Methodology

2.1 Acoustic configuration

We consider a two-dimensional semi-open acoustic propagation problem above a surface incorporating a finite-length surface treatment. Figure 1 shows the acoustic configuration, consisting of a sound source and a receiver located at coordinates (0, y_s) and (x₁ + x₂ + D, y_r), respectively, above a composite surface made up of three successive impedance regions. From left to right, the surface comprises: (i) a first segment, adjacent to the source and defined over the interval −∞≤x <  x₁; (ii) a central region of length D, corresponding to the metasurface treatment under study and characterized by an impedance Z, defined over x₁ ≤ x <  x₁ + D; and (iii) a final segment, adjacent to the receiver and defined over x₁ + D ≤ x <  +∞. Assuming a cylindrical source and invariance in the out-of-plane direction, the problem is formulated in two dimensions. This represents an infinitely extended line source in a geometry that is uniform along the transverse axis, where the essential propagation and surface interaction mechanisms are fully captured. The extension to three dimensions follows the same physical principles. To quantify the acoustic performance of the surface treatment in this configuration, the insertion loss (IL) at the receiver location is used as the evaluation metric. It is defined as follows:

Figure 1. Semi-open propagation geometry comprising a point source at distance y s , a receiver at distance y r , and a metasurface segment of length D embedded in the ground surface over a horizontal distance L.

$$\begin{array}{c}\hfill \mathrm{IL}=20{log}_{10}\left(\left|\frac{p_{}\text{ref}}{p}\right|\right)\end{array}$$(1)

where p and p _ref denote the sound pressure evaluated at the receiver location with and without the acoustic surface treatment, respectively. In the present configuration, IL is defined by comparing the sound pressure level at the receiver when the surface treatment is present to that obtained when it is replaced with a rigid boundary condition.

2.2 Helmholtz resonator metasurface

In this work, the surface treatment D is implemented using an acoustic metasurface, which enables tailored control of absorption [19, 20, 33], scattering [40], and wave manipulation [28, 41]. Such metasurfaces are typically composed of resonant subwavelength elements that provide a tunable effective acoustic impedance. Among these designs, metasurfaces based on HRs have proved particularly effective for sound control applications.

As shown in Figure 2, the metasurface used in this work consists of an array of N parallel slit-type HRs. Each resonator i is defined by its neck height h_n, i, neck length l_n, i, cavity length l_c, i, cavity height h_c, i, and spacing d_i. They are distributed along a segment D, where $D=(h_{n,1}+h_{n,N})/2+\sum _{i=1}^N{d}_i$.

Figure 2. Metasurface formed by an array of slit-type HRs. Each resonator i is defined by its neck length l n, i and width h n, i , as well as its cavity length l c, i and width h c, i . Two successive resonators, i and i + 1, are separated by a center-to-center distance d i .

A Helmholtz resonator (HR) is a structure composed of two rigid parts: a neck and a cavity connected together. In this configuration, the air in the neck acts as an oscillating mass when excited by the incident wave, while the air in the rear cavity acts as a spring due to compression. The resulting dynamic system can be described as an equivalent mass-spring system [20]. The first resonance frequency f_r, i of the i-th two-dimensional slit HR can be approximated using the classical analytical formulation [16, 20, 42]:

$$\begin{array}{c}\hfill f_{r,i}=\frac{c_0}{2\pi }\sqrt{\frac{h_{n,i}}{h_{c,i}{l}_{c,i}(l_{n,i}+{\delta }_i)}},\end{array}$$(2)

where c ₀ and δ _i are, respectively, the speed of sound and an end correction accounting for the inertial effects caused by the cross-sectional discontinuity between the neck and the cavity.

Assuming plane-wave propagation in the resonator neck and cavity, modeled in two dimensions as a slit of constant cross-section, and considering characteristic dimensions smaller than the acoustic wavelength but larger than the viscous and thermal boundary-layer thicknesses, the fluid inside the neck can be described by a complex effective density and bulk modulus according to Stinson’s model [43]:

$$\begin{array}{cc}\hfill {\rho }_{\text{eff},i}(\omega )& ={\rho }_0{[1-\frac{tanh\left(\frac{h_{p,i}}{2}{G}_{\rho }(\omega )\right)}{\frac{h_{p,i}}{2}{G}_{\rho }(\omega )}]}^{-1},\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill K_{\text{eff},i}(\omega )& =K_0{[1+(\gamma -1)\frac{tanh\left(\frac{h_{p,i}}{2}{G}_K(\omega )\right)}{\frac{h_{p,i}}{2}{G}_K(\omega )}]}^{-1},\hfill \end{array}$$(4)

where ρ _eff, i and K _eff, i are the complex equivalent density and the complex equivalent bulk modulus of the ith resonator, respectively, ω is the angular frequency; γ = C _P /C _V is the specific heat ratio of air and h _p, i being the neck and cavity slit height, h _n, i and h _c, i respectively. Similarly, the functions G _ρ (ω) and G _K (ω) are calculated as follows:

$$\begin{array}{cc}\hfill G_{\rho }(\omega )& =\sqrt{\frac{i\omega {\rho }_0}{\eta }},\hfill \end{array}$$(5)

$$\begin{array}{cc}\hfill G_K(\omega )& =\sqrt{\frac{i\omega Pr{\rho }_0}{\eta }},\hfill \end{array}$$(6)

where K ₀ = γ P ₀, and P ₀, η, and ρ ₀ are the isostatic modulus of elasticity, static pressure, dynamic viscosity, and density of air, respectively.

Then the impedance and the equivalent wavenumber of the medium (Z _eff,i , k _eff, i ) can be recalculated as follows:

$$\begin{array}{cc}\hfill Z_{\mathrm{eff},\mathrm{i}}(\omega )& =\frac{1}{\varphi }{\left({\rho }_{\text{eff},i}{K}_{\text{eff},i}\right)}^{1/2}\hfill \end{array}$$(7)

$$\begin{array}{cc}\hfill k_{\mathrm{eff},i}& =\omega {\left(\frac{{\rho }_{\text{eff},i}}{K_{\text{eff},i}}\right)}^{1/2}\hfill \end{array}$$(8)

with ϕ being the width of the slit (h _n, i , h _c, i ) relative to its periodicity d _i . The acoustic metasurface can be modeled using an equivalent model based on the transfer matrix method (TMM) [44]. In this approach, each part of the unit cell corresponding to a Helmholtz resonator is represented by equivalent series impedances that form a global matrix T _i described as:

$$\begin{array}{cc}\hfill {\mathbf{T}}_i& ={\mathbf{T}}_{\mathrm{\Delta }{l}_{n,i}^{(\text{out})}}{\mathbf{T}}_{n,i}{\mathbf{T}}_{\mathrm{\Delta }{l}_{n,i}^{(\text{in})}}{\mathbf{T}}_{c,i},\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill {\mathbf{T}}_{\mathrm{\Delta }{l}_{n,i}^{(\text{out})}}& =\left[\begin{array}{cc}1& \mathrm{j}{k}_{n,i}\mathrm{\Delta }{l}_{n,i}^{(\text{out})}{Z}_{n,i}\\ 0& 1\end{array}\right],\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill {\mathbf{T}}_{n,i}& =\left[\begin{array}{cc}cos(k_{n,i}{l}_{n,i})& \mathrm{j}{Z}_{n,i}sin(k_{n,i}{l}_{n,i})\\ \mathrm{j}\frac{sin(k_{n,i}{l}_{n,i})}{Z_{n,i}}& cos(k_{n,i}{l}_{n,i})\end{array}\right],\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill {\mathbf{T}}_{\mathrm{\Delta }{l}_{n,i}^{(\text{in})}}& =\left[\begin{array}{cc}1& \mathrm{j}{k}_n\mathrm{\Delta }{l}_n^{(\text{in})}{Z}_{n,i}\\ 0& 1\end{array}\right],\hfill \end{array}$$(12)

$$\begin{array}{cc}\hfill {\mathbf{T}}_{c,i}& =\left[\begin{array}{cc}cos(k_{c,i}{l}_{c,i})& \mathrm{j}{Z}_{c,i}sin(k_{c,i}{l}_{c,i})\\ \mathrm{j}\frac{sin(k_{c,i}{l}_{c,i})}{Z_{c,i}}& cos(k_{c,i}{l}_{c,i})\end{array}\right],\hfill \end{array}$$(13)

where ${\mathbf{T}}_{\mathrm{\Delta }{l}_{n,i}^{(\text{out})}}$, T _{n , i} , ${\mathbf{T}}_{\mathrm{\Delta }{l}_{n,i}^{(\text{in})}}$, and T _{c , i} correspond, respectively, to the correction at the outer end of the neck, the neck itself, the correction at the inner end of the neck, and the cavity. Z _n, i , k _n, i , Z _c, i , and k _c, i correspond to the equivalent impedance and wavenumber of the neck and the cavity obtained through equations (7) and (8). Each correction expression can be found in [44, 45].

Finally, the surface impedance of the metasurface is calculated using:

$$\begin{array}{c}\hfill Z_s=\frac{{\mathbf{T}}_{\mathbf{11}}}{{\mathbf{T}}_{\mathbf{21}}}·\end{array}$$(14)

From this expression, the reflection coefficient and the absorption coefficient can be derived from the surface impedance and the incidence angle θ as follows:

$$\begin{array}{cc}\hfill R(\theta )=|R|e^{j{\varphi }_r}& =\frac{Z_s-{\rho }_0{c}_0/\mathrm{c}os(\theta )}{Z_s+{\rho }_0{c}_0/\mathrm{c}os(\theta )},\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill \alpha (\theta )& =1-{|R(\theta )|}^2.\hfill \end{array}$$(16)

Similarly, the average absorption coefficient in the angular band is described as:

$$\begin{array}{c}\hfill \overline{\alpha }(\theta )=\frac{1}{N_{\theta }}\sum _{n=1}^{N_{\theta }}\alpha ({\theta }_n)\end{array}$$(17)

with N_θ being the number of incidence angles.

2.3 Numerical models

The acoustic configuration is modeled using a frequency-domain finite element method (FEM). The numerical models used in this study are shown in Figures 3 and 4. The semi-open propagation configuration shown in Figure 3 is used to evaluate IL, the central performance metric of this work. This numerical model accounts for the impedance discontinuity caused by the metasurface and interference effects inherent to the semi-open propagation domain. It includes the monopole source, the receiver, and the composite surface with the metasurface treatment, as defined in Figure 1. For simplicity, identical resonators are considered, so h_n, i = h_n, h_c, i = h_c, l_n, i = l_n, and l_c, i = l_c for i = 1, 2, 3, while the spacing between resonators is assumed to be uniform, i.e., d_i = d.

Figure 3. 2D finite element model used to evaluate the IL at the receiver location in a semi-open propagation domain.

Figure 4. 2D finite element model used to evaluate the omnidirectional absorption coefficient of a Helmholtz resonator unit cell with periodic boundary conditions (PBC), representing an infinite periodic array under plane-wave excitation.

The acoustic excitation is generated by a monopole source with a pressure amplitude set at p = 1 Pa. To replicate the semi-open propagation conditions, the computational domain is truncated using Perfectly Matched Layers (PMLs) to eliminate artificial reflections at the outer boundaries [46]. Thermoviscous losses within the metasurface are included by applying Stinson’s model in the three resonator necks and cavities as described in Section 2.2.

To provide a reference framework corresponding to the conventional design strategy based on intrinsic metasurface properties and to compare with the analytical model, the omnidirectional incidence absorption coefficient of the HR (α) is evaluated numerically using a periodic unit cell model. This configuration represents the idealized framework commonly used to characterize the intrinsic acoustic properties of metasurfaces. The model consists of a single resonator with periodic boundary conditions that reproduce the collective behavior of an infinite periodic array. This approach estimates the effective acoustic response of the metasurface formed by multiple resonators, consistent with the periodic arrangement considered in the semi-open propagation model.

The configuration is illustrated in Figure 4. Periodic boundary conditions are applied to the lateral boundaries with a spacing equal to d, corresponding to the resonator spacing defined in the semi-open configuration (Fig. 3). The excitation is a plane wave at an incidence angle from 0 to 85°. The incident wave is modeled using the Background Pressure Field condition.

The HR geometry and loss modeling are identical to those used in the insertion-loss model. The difference lies not only in the number of resonators represented but, more fundamentally, in the physical system being modeled. The IL model describes a finite array of three resonators embedded in a bounded surrounding environment, whereas the present setup uses periodic boundary conditions (PBC) with lateral spacing equal to d to represent an infinite periodic arrangement of identical resonators, thereby suppressing edge effects. This model thus provides a reference configuration for comparison with the approach based on IL introduced previously.

Thermoviscous losses in the metasurface and at all fluid-solid interfaces are accounted for in all models. To achieve this, a low reduced frequency (LRF) model is applied to the necks and cavities of the metasurface [47], while a boundary layer impedance (BLI) boundary condition is applied to all fluid-solid interfaces [48].

Since accurate prediction of IL and absorption peaks requires sufficient spatial resolution, the mesh size in both numerical models is chosen according to the criterion proposed in [49]:

$$\begin{array}{c}\hfill h_{\max }=\frac{c_0}{10f_{\max }},\end{array}$$(18)

where h _max denotes the maximum element size, and f _max is the maximum frequency considered in the study. The mesh is further refined locally in critical regions, such as the resonator necks and near the source, to ensure convergence of the solution. Several convergence studies have been conducted on these models and have shown that this criterion is sufficient to obtain adequate results.

2.4 Optimization procedure

An optimization procedure is applied to the semi-open propagation model used to evaluate IL presented on Section 2.3 as well as to analytical models based on the TMM. The aim is to identify the resonator parameters that yield optimal acoustic performance within each modeling framework. In the semi-open model, the objective is to maximize IL at a specified receiver location, directly targeting acoustic performance in the full propagation environment.

The optimization procedure is structured as a numerical design loop in which the variables h _n , h _c , l _n , l _c and d are iteratively updated to optimize a prescribed acoustic performance metric in each numerical model. For each candidate set of variables, the acoustic response of the system is evaluated using a frequency domain FEM formulation, which computes the pressure field in the corresponding numerical model. The objective function, either the omnidirectional absorption coefficient, reflexion coefficient or the IL at a specified receiver location, is then calculated from the simulated acoustic field. Based on this evaluation, the variables are updated and the process is repeated until convergence is achieved. In this way, the FEM model establishes a direct link between the design variables and the resulting acoustic performance within each modeling framework.

The optimization framework is implemented in a Python 3.11 environment integrated with COMSOL Multiphysics v6.3, the FEM software used in this study. Interaction between the two platforms is enabled by the MPh 1.3 module, which provides programmatic access to and control of the FEM model. This setup allows automatic updating of design variables at each iteration, regeneration of the geometry, remeshing, and execution of frequency-domain analysis without manual intervention.

The algorithm used is a differential evolution (DE) algorithm from the family of metaheuristic algorithms [50]. Developed in the 1990s based on the work of Storn and Price [51], this algorithm is inspired by Darwin’s theory of evolution and genetic principles such as crossover, mutation, and selection to find the global optimum of a function [50, 52]. The advantages of this algorithm are that it does not require gradient computation, allowing it to work with non-differentiable functions, its simplicity of implementation, and its ability to efficiently explore a large portion of the solution space [50].

Figure 5 summarizes the optimization workflow. Starting with an initial population of candidate geometries, the acoustic response is evaluated using the corresponding numerical model, and the population is iteratively updated according to the DE algorithm until the stopping criterion is met. The main steps of the differential evolution (DE) algorithm are described below [50]. For each generation G, each individual solution, denoted by j, is represented by a vector of dimension N, corresponding to the number of parameters, such that

Figure 5. Workflow of the DE optimization loop.

$$\begin{array}{c}\hfill {\mathbf{x}}_j^G=({\mathbf{x}}_{j,1},{\mathbf{x}}_{j,2},\dots ,{\mathbf{x}}_{j,N}),\end{array}$$(19)

with j = (1, 2, …, N P), where N P is the total population size. This size is calculated by multiplying the number of parameters N by P, the population size parameter. In our study, the parameter vector for individual x _j is written as follows:

$$\begin{array}{c}\hfill {\mathbf{x}}_j=(h_{n,j},l_{n,j},h_{c,j},l_{c,j},d_j),\end{array}$$(20)

corresponding to the geometric parameters of the metasurface to be optimized, as described in Figure 2.

The first phase of the algorithm is initialization, where the initial population is randomly generated in the parameter space as follows:

$$\begin{array}{c}\hfill {\mathbf{x}}_j^0={\mathbf{x}}_{\text{min}}+\text{random}(0,1)\times ({\mathbf{x}}_{\text{max}}-{\mathbf{x}}_{\text{min}}),\end{array}$$(21)

where x _min and x _max are the minimum and maximum values of the parameter x, respectively, and rand(0, 1) denotes a random number uniformly distributed between 0 and 1.

From this initialization, for each target individual ${\mathbf{x}}_j^G$ in generation G, a mutant vector ${\mathbf{v}}_j^{G+1}$ is generated by combining the differences between other individuals in the population:

$$\begin{array}{c}\hfill {\mathbf{v}}_j^{G+1}={\mathbf{x}}_{r_1}^G+F\times ({\mathbf{x}}_{r_2}^G-{\mathbf{x}}_{r_3}^G),\end{array}$$(22)

where F is the scale factor, and r ₁, r ₂, and r ₃ are distinct random indices chosen from the population. From this mutant, a cross is made between the target vector and the mutant to produce a trial vector denoted $u_j^{G+1}$ such that

$$\begin{array}{c}\hfill u_{j,k}^{G+1}=\{\begin{array}{cc}{v}_{j,k}^G\hfill & \text{if random}(0,1)\le \mathrm{CR}\text{or}k=k_{\text{random}}\hfill \\ x_{j,k}^G\hfill & \text{otherwise},\hfill \end{array}\end{array}$$(23)

where CR and k _random correspond to the crossover rate and a random index, respectively, to ensure that at least one parameter is taken from the mutant vector.

Finally, a selection is made by comparing the trial vector with the target vector. If the trial vector yields a better result, it replaces the target vector in the population, so

$$\begin{array}{c}\hfill {\mathbf{x}}_j^{G+1}=\{\begin{array}{cc}{\mathbf{u}}_j^{G+1}\hfill & \text{if}{J}_{\mathrm{obj}}({\mathbf{u}}_j^{G+1})\le J_{\mathrm{obj}}({\mathbf{x}}_j^G)\hfill \\ {\mathbf{x}}_j^G\hfill & \text{otherwise},\hfill \end{array}\end{array}$$(24)

where J _obj(⋅) is the objective function to be minimized.

For the IL optimization calculated in the model shown in Figure 3, the objective function is defined at a prescribed target frequency f₀ as

$$\begin{array}{c}\hfill J_{\mathrm{obj}}^{\mathrm{IL}}=-\mathrm{IL}(f_0),\end{array}$$(25)

so that maximizing the IL is reformulated as minimizing the objective function. Convergence of the optimization is assessed using a relative tolerance criterion [53]. A tolerance of 0.01 is adopted. To ensure that the solution does not correspond to a local optimum, the optimization is repeated with increasing population sizes of 15 and 20, starting from an initial population of 10. When all runs converge to the same objective value, the solution is considered to represent the optimal configuration.

Similarly, for the omnidirectional absorption coefficient evaluated in the periodic numerical model described in Figure 4 and calculated in the TMM model in Section 2.2, the objective function is defined as

$$\begin{array}{c}\hfill J_{\mathrm{obj}}^{\alpha }=\frac{1}{N_{\theta }}\sum _{n=1}^{N_{\theta }}{(1-\alpha ({\theta }_n,f_0))}^2,\end{array}$$(26)

so that minimizing the objective function promotes absorption over a broad range of incident angles [44, 54]. A standard approach is to maximize the angular average of α from 0° to 88°. However, under the locally reactive impedance assumption, the absorption coefficient tends to drop significantly at near-grazing incidence. To mitigate this effect, we also consider maximizing the minimum value of α over the same angular range, thereby enforcing performance at grazing angles.

Another approach aimed at maximizing insertion loss in this type of problem can be formulated analytically by considering the metasurface within the full source-surface-receiver configuration. The goal is to produce destructive interference at the receiver between the direct wave and the wave reflected (scattered) by the metasurface. After mathematical derivation, the equivalent criterion is to ensure that |R|≈1 and that the phase of the wave reflected by the metasurface is in phase opposition with the direct wave, i.e., ϕ _r  = π + kΔr, with $\mathrm{\Delta }r=\sqrt{L^2+{(y_r-y_s)}^2}-\sqrt{{(x_r-x_1)}^2+y_r^2}$ denoting the path difference between the direct source-receiver path and the path involving reflection from the metasurface. The objective function can thus be written as:

$$\begin{array}{c}\hfill J_{\mathrm{obj}}^R=\frac{1}{N_{\theta }}\sum _{n=1}^{N_{\theta }}|R({\theta }_n,f_0)-e^{j{\varphi }_r}|.\end{array}$$(27)

The advantage of this formulation is that it accounts for the spatial configuration of the source, metasurface, and receiver in the optimization process. The same criterion is evaluated over the same angular range as in the absorption-based optimization.

The hyperparameters used in the DE algorithm for each optimization problems (IL_opt, α _opt and R _opt), including the maximum number of iterations, population size P, scaling factor F, and crossover rate CR, are listed in Table 1 and follow commonly adopted reference values in the DE literature [50, 52]. These settings are applied consistently to each models. Finally, to comply with the geometric constraints, the penalty method is used. This involves adding a penalty to the objective function to prevent solutions that do not satisfy the constraints from being considered.

3 Application and results

3.1 Configuration and parameters

To demonstrate the proposed optimization methodology, an environmental acoustics scenario involving traffic noise propagation near a road is used. In this configuration, the sound source represents road traffic flow, and the acoustic impact is evaluated at a receiver location. The metasurface is embedded in the ground surface between the source and the receiver. From a practical perspective, the proposed optimization framework enables the design of metasurfaces to reduce road traffic noise and, more generally, to mitigate noise pollution along roadsides.

The geometric parameters of the source and receiver positions follow the standard values used for vehicle pass-by noise measurements [55]. These parameters are summarized in Table 2, where y_s, y_r, L, and x₁ denote the source height, receiver height, source-receiver distance (L = x₁ + D + x₂), and source-metasurface distance, respectively, as shown in Figure 1.

The position of the metasurface relative to the source, x₁, is chosen based on previous studies of impedance discontinuities in road traffic noise scenarios with similar source-surface distances [56, 57].

The bounds defining the admissible range of the optimization variables are summarized in Table 3. For each parameter, a minimum value B_min and a maximum value B_max are specified to delimit the search space explored by the optimization algorithm. These limits are chosen to ensure geometrically feasible resonator configurations while remaining consistent with the target frequency range considered in this study. Constraints are also applied and are expressed below:

$$\begin{array}{c}\hfill d-h_c\ge 1[\mathrm{mm}],\end{array}$$(28)

which aims to prevent geometric configurations in which neighboring resonators overlap when the periodicity becomes comparable to or smaller than the cavity width. Another constraint is introduced to ensure that the resonator maintains Helmholtz-type behavior during the optimization process and does not degenerate into a quarter-wavelength resonator as the geometric parameters are varied. This constraint is expressed as:

$$\begin{array}{c}\hfill h_c-h_n\ge 20[\mathrm{mm}].\end{array}$$(29)

In this analysis, four optimized metasurface configurations are considered: (1) maximizing the IL at the target frequency f ₀ (IL_opt); (2) maximizing the mean absorption coefficient over a broad range of θ at f ₀, corresponding to omnidirectional absorption (${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$); (3) maximizing the minimum absorption over a broad range of θ, corresponding to omnidirectional absorption optimized for grazing angles (${\alpha }_{\mathrm{opt}}^{\min }$) and (4) minimizing the average difference between the reflection coefficient and a target reflection coefficient over a broad range of θ to produce destructive interference (R _opt).

While IL_opt is computed using the numerical model described in Figure 3, ${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$, ${\alpha }_{\mathrm{opt}}^{\min }$, and R _opt are optimized using the TMM. The resulting geometric configurations are then transferred, and the respective insertion losses are calculated using the model shown in Figure 3. The resulting geometric configurations are also verified using the numerical model shown in Figure 4 to calculate the absorption coefficient as a function of wave incidence angles.

The study is conducted for a target frequency f ₀ = 500 Hz, which is representative of the low-frequency range typically associated with road traffic noise. Focusing on a single frequency allows a clear comparison between the two design strategies. The influence of other target frequencies is discussed later in the paper.

3.2 Optimization convergence

Each optimization performed by TMM takes between 2 and 10 min, with durations proportional to the population size parameter. A calculation of the insertion loss using the model described in Figure 3 takes approximately 8 min to cover the range from 100 to 1000 Hz. Meanwhile, a calculation of insertion loss and absorption coefficient as a function of the angle of incidence at the target frequency f₀ took 2 s. For all configurations tested, the numerical model shown in Figure 3 consists of approximately 60 000 elements. In contrast, the model shown in Figure 4 requires approximately 1000 to 2000 elements. Finite element model calculations are performed on a HP EliteOne 840 Desktop PC running Windows 11 Pro. The system is equipped with an Intel(R) Core(TM) i7-14700 processor (2.10 GHz), 64 GB of RAM (5600 MHz), and an NVIDIA GeForce RTX 3050 Ti Laptop GPU (4 GB VRAM), while the TMM calculation is performed on a HP workstation equipped with an AMD Ryzen 7 PRO 8840U processor (8 cores, 3.30 GHz, Radeon 780M Graphics), 64 GB RAM (5600 MHz), and 2 GB of dedicated AMD Radeon(TM) Graphics memory.

The optimization process exhibited stable convergence of the objective function, with no significant oscillations observed in the final iterations. For completeness, the evolution of representative design parameters IL_opt and α_opt is shown in Figure 6, as all parameters displayed similar convergence behavior and are therefore not reported individually.

Figure 6. Convergence of the DE optimization for the ILopt (top) and ${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$ (bottom) configurations with 200 iterations each for f = f 0.

Top of Figure 6 displays the convergence behavior for IL_opt using three population sizes (P = 10, 15, and 20). In all cases, convergence is achieved after approximately 110 iterations, yielding an IL of about 7.4 dB at the target frequency. Bottom of Figure 6 shows the convergence behavior of the objective function associated with the mean absorption coefficient for different population sizes. Convergence occurs after approximately 20 iterations for a fitness value of $\overline{\alpha }(f_0)\approx 0.905$. In the case of α _opt ^min, the resulting fitness value indicates that a minimum absorption coefficient of 0.531 is expected over all incidence angles. Finally, the value obtained for R _opt is not discussed, since it is nonphysical in this case.

3.3 Optimized geometries

The optimized geometrical parameters and the corresponding values of IL(f ₀) obtained for the different population sizes are reported in Table 4. For IL_opt, the solutions obtained for the different population sizes are very similar. In particular, the parameters h_n, l_n, and d converge towards the upper or lower limits of the search bounds (see Tab. 3). The resulting IL is consistently close to 7.4 dB for every population size. This phenomenon is also observed for the R_opt case, with very similar geometric configurations for all variables and insertion loss values ranging from 5 to 5.9 dB. In contrast, the geometries obtained by ${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$ and ${\alpha }_{\mathrm{opt}}^{\min }$ differ more significantly. Nevertheless, ${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$ provides significantly higher insertion loss values ranging from 2.4 to 3.2 dB, compared to approximately 1.2 dB. Since the optimization results are largely insensitive to the population size once convergence is reached, only the best-performing insertion loss value configuration for each objective function is retained for subsequent analysis. These selected configurations are highlighted in bold in Table 4.

3.4 Comparison of optimized configurations

The acoustic responses of the four optimized configurations are now compared. To assess the differences between the four design strategies, each configuration is evaluated using the two performance metrics considered in this work, namely IL and α(θ). This cross-comparison allows examination of how each optimization performs when assessed with its defining metric, as well as when evaluated using the alternative metric.

Figure 7 presents the IL obtained for each optimized configurations. The IL is evaluated as a function of frequency for the metasurface optimized for each optimized configurations, allowing direct comparison between each design strategies.

Figure 7. Comparison of the insertion loss obtained for the four optimized configurations: ILopt, ${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$, ${\alpha }_{\mathrm{opt}}^{\min }$ and R opt.

Each configuration results in positive IL values, indicating a reduction in sound pressure level at the receiver compared to the rigid reference surface. However, a clear difference emerges between the various designs. At the target frequency f₀, the configuration obtained using IL_opt achieves a peak value of approximately 7.4 dB, whereas ${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$ and ${\alpha }_{\mathrm{opt}}^{\min }$ yield an IL values of about 4 dB and 6 dB, respectively. Similarly, R _opt yields a peak insertion loss of approximately 5.9 dB at f ₀, although its maximum is slightly shifted in frequency to about 6.6 dB. This demonstrates that maximizing absorption does not necessarily lead to optimal performance in the full propagation environment. Furthermore, taking the problem’s geometry into account yields results much closer to the optimal solution.

At frequencies slightly below the target frequency, each curve shows an increase in sound level (IL <  0), which can be attributed to constructive interference effects and the possible formation of a surface wave [58]. This behavior is consistent with the presence of a wave component propagating approximately parallel to the boundary, whose amplitude decreases exponentially with height and by $1/\sqrt{D+x_2}$, where D + x ₂ is the distance between the beginning of the surface discontinuity and the receiver [59]. Under these conditions, particularly at large distances and close to the surface, the contribution of this surface wave component may become significant [58]. Nevertheless, it is important to note that this traditional interpretation of the surface wave has recently been challenged by new research based on precise measurements [60].

At the target frequency, this amplification is followed by a peak in IL caused by destructive interference between the direct wave and the wave scattered by the metasurface. Beyond this frequency, the IL gradually decreases and approaches zero. Similar behavior has been reported in previous studies [14, 15].

Given the geometrical arrangement of the source, metasurface, and receiver, a reference IL of approximately 6 dB at f ₀ can be considered a baseline value. This level corresponds to the ideal case in which the specular reflected wave from the surface is completely suppressed, so only the direct field from the source contributes to the pressure at the receiver. However, due to the resonant behavior of the metasurface and the resulting interference between the direct and scattered fields, IL values exceeding this level may occur.

Furthermore, due to the distance between the source and the receiver, interference between the direct wave and the wave reflected by the ground creates a classical interference pattern in which the first destructive interference occurs at a much higher frequency than the target frequency. Consequently, in the frequency range considered here, the direct and reflected waves do not produce destructive interference, and the metasurface response remains essentially uncoupled from this interference mechanism.

The comparison between the TMM and FEM results is shown in Figure 8, where the absorption coefficient is plotted as a function of the angle of incidence at f = 500 Hz for each optimization strategies. For completeness, both the mean-optimized (${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$) and minimum-optimized (${\alpha }_{\mathrm{opt}}^{\min }$) responses are reported. Overall, good agreement is observed between the two approaches over most of the angular range, validating the numerical FEM model against the analytical TMM predictions. Small discrepancies appear at large angles of incidence, where the locally reactive assumption underlying the TMM model becomes less accurate, particularly for the minimum-based optimization. Nevertheless, the main trends are consistently captured, including the degradation of the mean-optimized response at near-grazing incidence and the improved robustness of the minimum-based criterion.

Figure 8. Comparison of the absorption coefficient obtained by TMM and FEM for the two optimized configurations: α mean and α min.

To better understand the mechanisms responsible for the differences observed in IL, the spatial distribution of the pressure and velocity field at the target frequency is analyzed. Figure 9 shows the comparison of the real part of the acoustic pressure field, Re{p}, and the acoustic particle velocity field v at f = f₀ for the two optimized configurations IL_opt and ${\alpha }_{\mathrm{opt}}^{\min }$. For the sake of simplicity and clarity, only these two configurations are presented, as IL_opt and R_opt produce nearly identical field distributions, while the mean-based and minimum-based α_opt configurations also show similar spatial patterns. Therefore, the selected cases represent distinct design strategies. The displayed fluid domain spans from 0 to 0.4 m in the y-direction and from 1.5 m to 2.5 m in the x-direction.

Figure 9. Comparison of the pressure field Re{p} and the acoustic velocity v at f = f 0 in the semi-open propagation model for the configuration: ILopt (top) and ${\alpha }_{\mathrm{opt}}^{\min }$ (bottom).

In Figure 9 the configuration optimized for insertion loss IL_opt (shown on top), pronounced interference patterns appear in the propagation domain above the metasurface. Regions of constructive and destructive interference form in the propagation domain, indicating significant redistribution of the scattered acoustic field. The pressure amplitudes above the metasurface, particularly near the resonator neck exits, are relatively high, indicating that the resonators act as secondary acoustic sources, re-radiating energy into the surrounding field. In addition, coupling occurs between the resonators, showing a phase shift between each resonator.

The redistribution of the acoustic field is also particularly evident in Figure 9, where the velocity fields, represented by black arrows, highlight the distinct spatial behavior of the two configurations. Indeed, a comparison of the configurations shows, in particular, that the metasurface resulting from configuration ${\alpha }_{\mathrm{opt}}^{\min }$ has a very localized effect on the velocity field at the outlet and within the neck, whereas in the case of configuration IL_opt, the field is significantly affected well beyond the neck of the metasurface, with notable changes in the field’s direction.

These observations indicate that the attenuation mechanism associated with the IL_opt configuration involves not only energy dissipation within the resonators but also the spatial redistribution of the acoustic field. In this configuration, the metasurface can be interpreted as a structure that redirects acoustic energy away from the receiver through scattering and interference effects, rather than serving solely as a dissipative absorber. Such effects are not captured by the omnidirectional absorption coefficient, which is defined under locally reacting and plane-wave assumptions. Therefore, optimization using R _opt offers a more physically consistent approach for this configuration, as it explicitly accounts for interference effects based on the known positions of the source and receiver. However, since this criterion does not fully account for all loss mechanisms present in the numerical model, which includes fluid-solid interaction and minor discrepancies between the TMM and FEM, the agreement with the IL_opt results remains limited.

3.5 Spatial robustness of the optimized configuration

To assess the sensitivity of IL_opt to receiver position, the zone of interest defined by (x, y)∈[0.5, 2]×[6, 8] m² in Figure 3 is analyzed (see Fig. 10). The insertion loss at the target frequency ranges from 5 to 8.5 dB within this region, with higher values observed for receiver positions closer to the metasurface, corresponding to smaller x values. This spatial distribution shows that the attenuation produced by the metasurface is not confined to a single point but extends over a significant region downstream of the surface treatment.

Figure 10. Insertion loss at f 0 as a function of receiver position within the zone of interest.

This spatial behavior can be attributed to the formation of a near-field shadow region immediately downstream of the metasurface, caused by the abrupt impedance discontinuity at x = x ₁. The resulting phase and amplitude redistribution of the scattered field leads to locally enhanced destructive interference and reduced sound pressure levels near the surface.

In the current configuration, the source-metasurface distance is fixed at x ₁ = 1.5 m. Smaller values of x ₁ would increase the insertion loss, as a larger portion of the acoustic field would interact with the impedance discontinuity and contribute to the formation of the shadow region [56]. These results indicate that, although the optimization is performed at a single receiver location, similar IL levels are achieved over a wide surrounding region.

4 Discussion

The methodology was further evaluated at different target frequencies f ₀ (300, 700, and 1000 Hz). In all cases, the same qualitative behavior is observed: the configuration obtained through IL_opt systematically yields higher IL values than those obtained by optimizing the absorption coefficient.

Nevertheless, when the same optimization bounds are maintained, the gap between the two optimizations tends to decrease as frequency increases. In fact, the amplitude of the IL peak also appears to depend on the selected frequency. For example, at 300 Hz, an IL amplitude of more than 10 dB is achieved, while at 1000 Hz, the maximum reaches 4.5 dB. Conversely, optimizing the absorption coefficient yields results typically in the range of 3 to 3.5 dB maximum IL for these different frequencies. This indicates that, although IL_opt remains more effective, the relative improvement becomes smaller at higher frequencies.

This trend can be attributed to the frequency dependence of the radiation impedance Z _r , which is influenced by both frequency and the geometrical parameters of the slits, including their width and periodicity [61]. In the presence of a source, this type of metasurface can be modeled as an array of baffled pistons [14]. In this interpretation, each slit acts as an individual radiating element that scatters part of the incident field, consistent with Figure 9, where the pressure scattered by the metasurface IL_opt can be interpreted as secondary acoustic sources. As frequency increases, the interaction between the resonators and the incident field changes, altering the balance among scattering, interference, and dissipation mechanisms, and ultimately reducing the achievable IL.

Although the present analysis provides clear insight into the influence of the acoustic environment on metasurface optimization, some limitations of the current framework should be noted. The study is based on a two-dimensional numerical model and considers a monopole source as well as optimization at a single target frequency. While these assumptions allow a clear comparison between intrinsic and environment-based optimization strategies, more complex acoustic scenarios involving broadband excitation, three-dimensional propagation, or extended source distributions may introduce additional mechanisms that could influence the optimal design.

5 Conclusion

This work aimed to identify the most appropriate strategy for designing the characteristics of an acoustic metasurface to maximize its performance in a specific propagation scenario, rather than optimizing an isolated metamaterial parameter such as the omnidirectional absorption coefficient. To this end, insertion loss optimization was performed at a target frequency f ₀ within the framework of a semi-open configuration comprising a source, a receiver, and a finite-sized metasurface. The results show that optimizing the acoustic metasurface within its full propagation configuration yields an IL that is 4.2 dB and 6.2 dB higher than that achieved by absorption-based optimization. The results of an optimization based on the reflection coefficient R, which incorporates the spatial configuration of the source, receiver, and metasurface, support this approach by capturing the interference mechanisms governing insertion loss. Notably, configurations that maximize IL exhibit a lower absorption coefficient than that obtained through absorption-based optimization, highlighting that high insertion loss is primarily achieved through interference and scattering rather than energy dissipation. Furthermore, the resonance frequency of this metasurface is higher than the target frequency, indicating a shift of the resonant behavior relative to the target frequency.

This behavior can be interpreted as a trade-off between thermoviscous dissipation, partly reflected in the absorption coefficient, and the spatial distribution of the scattered pressure field induced by the acoustic metasurface. In this configuration, the apertures of the resonators direct incident energy away from the receiver. These findings indicate that the absorption coefficient alone does not fully characterize the acoustic performance, particularly in semi-open configurations. Therefore, the metasurface must be analyzed within its complete acoustic environment when the objective is to maximize a given performance metric. This conclusion is consistent with recent works [38, 39], which show that perfect absorption, although often the primary objective, is not necessarily optimal in configurations involving non-normal or complex wave incidence.

The robustness of the approach was further evaluated at additional target frequencies selected below the first destructive interference of the geometrical configuration. In all cases, the same qualitative trends were observed. For spatial robustness, the increase in IL is not limited to the specific receiver position used during optimization but extends over the surrounding area. This behavior can be attributed to the formation of a near-field shadow zone downstream of the metasurface, resulting from the impedance discontinuity at the surface, which redistributes the acoustic energy in its vicinity.

Over recent decades, acoustic metamaterials have shown strong potential for sound control due to their unconventional acoustic properties. However, their design is often conducted independently of the full propagation environment in which they are intended to operate. This work takes an initial step towards integrating acoustic metamaterials within their actual application environment, demonstrating that performance metrics defined at the system level can result in substantially different optimal configurations. This perspective opens the possibility of extending metamaterial-based sound control strategies to more realistic and complex acoustic environments.

Acknowledgments

We would also like to thank the reviewers for their comments, which have greatly improved the content of this work.

Funding

This research was partially developed under grant PID2022-138321NBC22 funded by MICIU/AEI/10.13039/501100011033 and ERDF/EU.

Conflicts of interest

The authors declare no conflict of interest.

Data availability statement

Data are available on request from the authors.

Author contribution statement

Robin Mafféïs: Conceptualization, Methodology, Software, Validation, Formal analyzis, Investigation, Writing – Original Draft, Writing – revised draft. Judicaël Picaut: Conceptualization, Methodology, Funding acquisition, Supervision, Writing – review and editing. Ruben Picó: Conceptualization, Methodology, Funding acquisition, Supervision, Writing – review and editing.

References

All Tables

All Figures

	Figure 1. Semi-open propagation geometry comprising a point source at distance y s , a receiver at distance y r , and a metasurface segment of length D embedded in the ground surface over a horizontal distance L.
In the text

	Figure 2. Metasurface formed by an array of slit-type HRs. Each resonator i is defined by its neck length l n, i and width h n, i , as well as its cavity length l c, i and width h c, i . Two successive resonators, i and i + 1, are separated by a center-to-center distance d i .
In the text

	Figure 3. 2D finite element model used to evaluate the IL at the receiver location in a semi-open propagation domain.
In the text

	Figure 4. 2D finite element model used to evaluate the omnidirectional absorption coefficient of a Helmholtz resonator unit cell with periodic boundary conditions (PBC), representing an infinite periodic array under plane-wave excitation.
In the text

	Figure 5. Workflow of the DE optimization loop.
In the text

	Figure 6. Convergence of the DE optimization for the ILopt (top) and ${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$ (bottom) configurations with 200 iterations each for f = f 0.
In the text

	Figure 7. Comparison of the insertion loss obtained for the four optimized configurations: ILopt, ${\alpha }_{\mathrm{opt}}^{\mathrm{mean}}$, ${\alpha }_{\mathrm{opt}}^{\min }$ and R opt.
In the text

	Figure 8. Comparison of the absorption coefficient obtained by TMM and FEM for the two optimized configurations: α mean and α min.
In the text

	Figure 9. Comparison of the pressure field Re{p} and the acoustic velocity v at f = f 0 in the semi-open propagation model for the configuration: ILopt (top) and ${\alpha }_{\mathrm{opt}}^{\min }$ (bottom).
In the text

	Figure 10. Insertion loss at f 0 as a function of receiver position within the zone of interest.
In the text