Broadband impedance modulation via non-local acoustic metamaterials

Abstract Causality of linear time-invariant systems inherently defines the wave-matter interaction process in wave physics. This principle imposes strict constraints on the interfacial response of materials on various physical platforms. A typical consequence is that a delicate balance has to be struck between the conflicting bandwidth and geometric thickness when constructing a medium with desired impedance, which makes it challenging to realize broadband impedance modulation with compact structures. In pursuit of improvement, the over-damped recipe and the reduced excessive response recipe are creatively presented in this work. As a proof-of-concept demonstration, we construct a metamaterial with intensive mode density that supports strong non-locality over a frequency band from 320 Hz to 6400 Hz. Under the guidelines of the over-damped recipe and the reduced excessive response recipe, the metamaterial realizes impedance matching to air and exhibits broadband near-perfect absorption without evident impedance oscillation and absorption dips in the working frequency band. We further present a dual-functional design capable of frequency-selective absorption and reflection by concentrating the resonance modes in three frequency bands. Our research reveals the significance of over-damped recipe and the strong non-local effect in broadband impedance modulation, which may open up avenues for constructing efficient artificial impedance boundaries for energy absorption and other wave manipulation.

where ∆f is the half-maximum relative bandwidth, f r and λ r are the resonance frequency and wavelength, R a is the normalized acoustic resistance, and L denotes the thickness of HR. For perfect absorption when R a = 1, L can be solved and written as Make ∆f = f 2 − f 1 , then where λ 1 and λ 2 are the corresponding wavelengths, and F(λ) = 1/λ is the primitive function.
The two terms in parenthesis have a dimension of length, which reflects and quantifies the frequency response of the HR. The integral over the interval from λ 2 to λ 1 is obtained via the Newton-Leibniz formula. Obviously, Eq. (S3) possesses an analogous mathematical form to Eq. (2) in the main text. This equation implies that to keep the same L, the interval ∆λ = (λ 2 − λ 1 ) decreases (increases) with λ r increasing (decreasing), i.e. the bandwidth of the resonance peak is proportional to the resonance frequency. In this sense, the quality factor Q is equivalent to the causality principle. However, for broadband scenario, the absorption spectrum can not be quantified simply by a half-maximum bandwidth and thus the quality factor is invalid.
In addition, some conclusions can be drawn from the quantitative relation contained in the causality principle. First, for a certain thickness L, the absorption coefficient decreases (increases) with an increasing (decreasing) working bandwidth to keep the same integral value. Second, for a certain frequency band, the minimal thickness L min increases (decreases) with an increasing (decreasing) absorption coefficient. However, the definition of optimality is not restricted by certain bandwidths or absorption coefficients. For instance, a material exhibiting low absorption or narrow bandwidth can still be optimal as long as the condition L = L min is satisfied.

Supplementary Note 2
The response functions in Fig. 1 In Fig. 1b, we use auxiliary function as a mathematical fitting of the absorption spectrum.
Generally, it can be summarized by a form of function expressed as A = [tan −1 (af + b) + π/2]/c+ d . By carefully reducing the coefficients a, b, and c, and choosing an appropriate intercept d, three functions increasing monotonically with descending derivative are obtained. The functions corresponding to the green, red, and blue curves are where S EN is the cross-sectional area of the embedded neck, ξ = d/ min(a, b), and a and b are the length and width of the square tube. As ξ approaches 1, the end correction decreases and converts to a small positive value, while Eq. (S5) gives a negative one. This violation would cause obvious deviation in high frequency regime. To amend this violation, we revise the end correction for ξ ≥ 0.4 as where a=2/3, b=-26/15, c=16/15 are the fitting coefficients. Fig. S1 shows that after the revision, the end correction agrees well with the result presented by Ingard [42] when ξ ≥ 0.4.

Surface impedance analysis based on the transfer matrix method
As mentioned in the main text, the lumped-parameter impedance model is not valid for CNEHR in high frequency regime. Therefore, we developed the transfer matrix method to analyze the surface impedance of the CNEHR. Then, according to the acoustic continuity conditions (the continuity of sound pressure and volume velocity), the transfer matrix at different position inside the CNEHR can be deduced to relate the sound pressure and volume velocity at different positions (Fig. S2). This involves several equations expressed as when i takes different number (8,7,6,5,4,3,2), denotes the transfer matrix, sound pressure and volume velocity at x 8 ,x 7 ,x 6 ,x 5 ,x 4 ,x 3 ,x 2 ,x 1 .
Note that all the complex wave numbers and complex characteristic impedance in each part of the structure are calculated after considering the viscous and thermal dissipation in the acoustic boundary layer on the wall [43]. The following section gives the specific transfer matrix. The details of the deduction can be found in our previous work [26]. At the upper port of EN 1 (x 8 ), the end correction of the acoustic impedance can be added by employing the transfer matrix where k and ω are the wave number and angular frequency of air in free space, Z 0 = ρ 0 c 0 and η are the characteristic impedance and dynamic viscosity of static air, δ 1 is the end correction of EN 1 relative to incident field, kZ 0 2η/ωρ 0 = √ 2ηωρ 0 represents the surface resistance at the port of a circular tube, and S EN 1 is the cross-sectional area of EN 1 . Inside the channel of EN 1 (x 7 ), the transfer matrix is where k EN 1 and Z EN 1 = ρ EN 1 c EN 1 are the complex wave number and complex characteristic impedance of air inside EN 1 . At the lower port of EN 1 (x 6 ), the region A around EN 1 splits the wave into two flows, where the boundary conditions can be expressed as P Considering the end correction of acoustic impedance and the boundary conditions, it can be deduced that where k A and Z A = ρ A c A are the complex wave number and complex characteristic impedance of air in region A, δ 2 is the end correction of EN 1 relative to the cavity, and x 7 x 4 x 5 S A = wh − π(d 1 /2 + b) 2 is the cross-sectional area of region A. Inside the region B (x 5 ), the transfer matrix is where k B and Z B = ρ B c B are the complex wave number and complex characteristic impedance of air inside region B, and S B = wh is the cross-sectional area of region B.
Similarly, the transfer matrix at x 4 , inside EN 2 (x 3 ), at x 2 , and inside region E (x 1 ) can be obtained as follows: where k C and Z C = ρ C c C are the complex wave number and complex characteristic impedance of air in region C, δ 3 is the end correction of EN 2 relative to the cavity, and S EN 2 are the cross-sectional areas of regions C and EN 2 ; where k EN 2 and Z EN 2 = ρ EN 2 c EN 2 are the complex wave number and complex characteristic impedance of air inside EN 2 , and S EN 2 is the cross-sectional area of EN 2 ; where k D and Z D = ρ D c D are the complex wave number and complex characteristic impedance of air in region D, and S D = S C is the cross-sectional area of region C; where k E and Z E = ρ E c E are the complex wave number and complex characteristic impedance of air in region E, and S E = S B is the cross-sectional area of region E. Then the total transfer matrix T can be expressed as: Combined with the boundary condition V b = 0 (hard boundary), the surface impedance of the CNEHR is written as Z an = T 1 /T 3 .
To verify the acoustic impedance obtained above, we fabricated a 3D-printed sample of the CNEHR with the geometric parameters t 1 =5 mm, t 2 =20 mm, t 3 =15 mm, d 1 =3.9 mm, d 2 =4.8 mm, w=10 mm, h=10 mm, L 1 =60 mm, L 2 =40 mm, and b=1.5 mm, and assembled a 2 × 2 array with it (Fig. S2b). It can be observed from  and k is the wave number of air in free space. After considering the intrinsic loss of the square tube, the reflection coefficient should be rigorously 1: where H c I = e −jk s and H c R = e jk s are the transfer functions of the incident and reflected waves after considering the loss, and k is the complex wave number in the square tube.
The complex wave number k can be solved with Eq. (S17) and Eq. (S18), which yields the reflection coefficient r sp of the tested sample: where H sp 12 = p sp 2 /p sp 1 with p sp 2 and p sp 1 being the sound pressure detected by Mic. 1 and Mic. 2 when the tested sample is placed at the end of the square tube. Then the absorption coefficient and acoustic impedance of the sample can be obtained with A = 1 − |r sp | 2 and Z a = (ρ 0 c 0 /S + r sp )/(ρ 0 c 0 /S − r sp ).

Supplementary Note 5
Reduction of excessive response via increasing the number of subunits Here, we will demonstrate that increasing the number of subunits can reduce the excessive response of the metamaterials. We investigate the absorption spectra of metamaterials consisting of different numbers of subunits. For comparison, the absorption spectra of three metamaterials with 20, 25 and 36 subunits are shown in Fig. S4. According to the causality principle, the red, black and blue spectra respectively lead to minimal thicknesses of 97.3 mm, 97.3, and 96.7 mm, which are close to the actual thickness (100 mm). Therefore, all the structures reach the minimal thickness. However, different numbers of subunit are observed to yield distinct excessive responses. Actually, the number of subunit puts a strong constraint on the resonance characteristics of the metamaterial. Fewer subunits lead to bigger crosssectional area for each subunit, which causes comparatively larger quality factor for each resonance mode. Generally, this signifies that the absorption peak contributed by each mode near 320 Hz is flatter. Hence, after coupling all the modes, the overall absorption spectra of the metamaterials with 20 and 25 subunits from 0 Hz to 320 Hz are inherently flatter than the spectrum of the one with 36 subunits, indicating larger excessive response and a bigger slope intuitively. Because the minimal thicknesses are almost the same, a higher absorption coefficient under 320 Hz inevitably leads to lower absorption above 320 Hz. Eventually, the absorption from 320 Hz to 6400 Hz is weakened and oscillates evidently. Therefore, to reduce the excessive response, we can employ a relatively larger number of subunits. It should be noted that when the thickness of the subunits' side walls is fixed, employing too many subunits will reduce the volume fraction and increase the minimal thickness or decrease the absorption efficiency. The number of subunits should be a suitably large number.
Excessive response is an intrinsic feature of the coupled-resonant system that can be weakened by altering the resonance modes to improve absorption performance in the operating spectrum. Especially in low frequency impedance-matching design, the excessive response plays a more important role in the absorption performance modulation. Boundary" due to the huge impedance mismatch between air and the resin walls. Besides, the inner surfaces of the structures are set as "Thermoviscous Boundary Layer Impedance Interface" to calculate the thermoviscous effect. An incident plane wave guided by a square duct normally impinges on the metamaterial (Fig. S5). As a consequence, the pressure and energy flux field in Fig. 3e are obtained by setting a Cut Plane 0.1 mm above the metamaterial. The energy flux is defined as I x(y) = 1 2 Re(p * v x(y) ), where p is the sound pressure, v x(y) is the particle velocity with respect to the x(y) direction, and I x(y) is the energy flux with respect to the x(y) direction.

Influence of the arrangement of subunits
To explore the influence caused by the arrangement of subunits, we simulated the absorption spectra of two sets of comparative structures (Fig. S6). The first structure was composed of 2 subunits. We changed the arrangement in terms of the distance between the two subunits to evaluate the effect of distance on absorption. The second structure was composed of 12 subunits. We changed the arrangement in terms of the order of subunits to evaluate the effect of multi-subunits' order on absorption. The resulting absorption spectra  Table III. The parameters of the 12-subunit structure are those listed in Supplementary Note 11 Table ?? unit 13-24.
in both cases have small deviation after the change in subunit order, which verifies that subunit arrangement is not the primary influence factor in the non-local coupling and can be ignored for our structures in the interest of frequency. Additionally, we note that although the coupling among the subunits is not sensitive to their arrangement in our design, the coupling is greatly affected by other important factors such as the resonance frequency, mode density and acoustic resistance.

Resonance mode of each subunit of the metamaterials in the main text
To further demonstrate the near-continuous resonance mode of the metamaterial in Fig. 3 and Fig. 4 in the main text, we calculate the absorption spectra of all the subunits using where A n is the absorption coefficient of the n-th subunit. The overall absorption spectra of the metamaterial and the individual CNEHRs are shown in (Fig. S7). The numerous absorption peaks demonstrate the near-continuous resonance mode of the metamaterial.

Design strategy for the non-local metamaterials
In our work, we appropriately harnessed the nonlocality in the near field to suppress the impedance oscillation at anti-resonance frequency. However, with the bandwidth broadening, more subunits are required to ensure broadband modulation of the impedance. Meanwhile, the coupling becomes more complicated. Thus, we employed optimization algorithm (genetic algorithm) to accomplish the non-local metamaterials in the main text. In the optimization, specific resistance and reactance values are expected in the target frequency band. As stated in the main text, the impedance matching condition requires the resistance to be 1 and the reactance to be 0. However, the inevitable small oscillation of the impedance will easily results in a resistance smaller than 1 in some frequency bands. Hence, to ensure the over-damped condition, the objective of resistance is larger than 1 in this study. The objective with respect to resistance and reactance is written as where f i is the chosen frequencies in the target band, and r a (f i ) and x a (f i ) is the normalized resistance and reactance at f i . During the iteration, the objective gradually decreases and reaches a minimum where the optimized parameters of the device are obtained. The obtained results lead to the optimal absorption performance under our concept-based design, i.e., the given type and number of subunits. Because our results demonstrate that the relatively more intensive mode density can better suppress the dispersion of impedance and a smaller volume fraction of the material (φ) allows a thinner structure, we think it is possible to improve the structure by further increasing the mode density and decreasing the volume fraction.
Note that in our design, keeping the same length L for all the CNEHRs ensures that the metamaterial possesses a bigger effective volume fraction, which enables the metamaterial to approach the minimal thickness governed by the causality principle.

Supplementary Note 10
Effect of non-local coupling on impedance manipulation The coupled mode theory is provided to describe the effect of non-local coupling on impedance matching. For simplicity but without loss of generality, we investigate a system coupled with three modes, where the mode amplitudes are expressed as where S + i represents the incident waves, ω 1(2,3) is the resonance frequency of mode 1(2,3), γ 1(2,3) is the corresponding radiative decay rate, and Γ 1(2,3) represents the dissipative decay rate associated with the intrinsic loss. The third group of terms in Eqs. S22-S24 represents interactions among the modes, i.e., the non-local coupling effect. In the absence of these coupling terms, the resonances become local and we denote the local resonances as a 1,0 , a 2,0 , and a 3,0 .
Figures. S8-S10 illustrate the influences of non-local coupling with increasing mode density, where we can find that the real and imaginary parts of the resonances' amplitudes are altered due to the interactions among the resonances. Here, the individual non-local factor (F ni ) and the overall non-local factor (F no ) are introduced to evaluate the influence of the non-local coupling on the individual modes and the overall system, respectively, which are expressed as  S9. Mode amplitude of the local and non-local acoustic systems under the relatively moderate mode density. ω 1 /(2π), ω 2 /(2π) and ω 3 /(2π) are 480 Hz, 600 Hz and 745 Hz, respectively. γ 1 , γ 2 , and γ 3 are 0.04ω 1 , 0.04ω 2 , 0.04ω 3 , respectively. Γ 1 , Γ 2 and Γ 3 are 0.08ω 1 , 0.08ω 2 , 0.08ω 3 , respectively.
above. Further, the reflection coefficient of the coupled system, r, is given as [46] Accordingly, the acoustic impedance (Z) of the system can be written as As shown in Fig. S11, it can be observed that with increasing mode density, the individual non-local factors of mode 1 (F ni1 ), mode 2 (F ni2 ) and mode 3 (F ni3 ) and the overall non-local factor (F no ) of the corresponding systems are enhanced. These results indicate that the relatively more intensive mode density can support stronger effect of non-local coupling on each mode. Furthermore, it can be found evidently from Figs. S11(b)-(c) that the dispersion of impedance occurs at the frequencies with the small F no , which will reduce the efficiency of impedance modulation and result in absorption dips for sound-absorbing metamaterials.
However, this problem can be well solved when we construct the system with denser modes to support strong non-local coupling. In this manner, the dispersion of impedance can be suppressed, leading to improved modulation for impedance matching.  S11. (a)-(c) The individual non-local factor, overall non-local factor, and acoustic impedance of the system with relatively sparse mode density. (d)-(f) The individual non-local factor, overall non-local factor, and acoustic impedance of the system with relatively moderate mode density. (g)-(i) The individual non-local factor, overall non-local factor and acoustic impedance of the system with relatively intensive mode density.

Geometric parameters of several structures in the main text
To support our theory, several typical CNEHRs were designed to demonstrate the overdamped nature of optimal sound-absorbing system. In addition, we presented two non-local metamaterials to realize broadband absorption. The wall thicknesses of the front, middle, back panel (b 2 ) and the embedded neck (b 1 ) of the metamaterials are set as different values to guarantee maximal use of the space as well as sufficiently strong stiffness. The 3D printer can technically reach a precision of 0.1 mm, which is at the same level as that of the parameters of the presented acoustic metamaterials. The resonance frequency may slightly shift from the theoretical value due to fabrication imperfections, but it would not significantly affect the overall impedance modulation owing to the intensive mode density. The geometric parameters are listed in the following tables.  Fig. 2d-e. Note that since the acoustic impedance of the subunits in the metamaterials will be enlarged due to the change of surface area, their neck lengths t 1 are relatively shorter than the single structure presented here.