Broadband acoustic double-zero-index cloaking with coupled Helmholtz resonators

Acoustic double-zero-index metamaterials (DZIM) characterized by extremely large phase velocity along with no phase changes of the wave propagation inside the materials have received tremendous attention due to the fascinating physics and potential applications. However, due to the requirement of the degeneracy of dipolar and monopolar resonances and the available resonance-induced losses, the realization of highly efficient and broadband near-zero index metamaterials is still facing challenges. Here we report that by coupling two identical Helmholtz resonators with a connecting channel, acoustic DZIM can be realized. Owing to the presence of a connecting tube, the system can generate the dipolar mode that is independently tunable and the monopolar mode that is virtually unchanged. It thereby makes the mass density (ρ) and the reciprocal of the bulk modulus (1/B) simultaneously crossing zero possible. We numerically obtain the transmission and phase, and then calculate the effective mass density and bulk modulus, which agree remarkably well with the experimental results. Finally, we successfully cloak a rectangle block inside a two-dimensional waveguide grafted by the designed acoustic DZIM array of unit cells. A broadband cloaking is experimentally demonstrated at 1865–1925 Hz, which can offer potential possibilities for vast practical applications.

Here, we report the experimental realization of acoustic DZIM structure by designing a cylindrical tube to connect two identical HRs. The double-zero (DZ)-index can be generated by fine-tuning the size of the connecting tube. We successfully cloak a rectangle block in a two-dimensional waveguide grafted by the designed 5 × 9 array of the DZIM unit cells. In the numerical and experimental validations, the phase distribution of the propagation waves maintains virtually unchanged and the impedance matches well with air. A relatively broadband acoustic cloaking is experimentally demonstrated in the frequency band of 1865-1925 Hz. Interestingly, the designed DZIM cloaking system has no intrinsic limit to the number of unit cells, and the corresponding array can be designed according to real-world applications with the advantages of simple fabrication and assembly.

Design of the acoustic DZIM structure
We present a numerical and experimental realization of an acoustic DZIM structure with 1/B and ρ being simultaneously crossing zero, as displayed in figure 1(a). It is well known that a single HR has the negative bulk modulus [18] induced by the monopolar mode [42]. Therefore, 1/B moves from the negativity branch to the positivity branch, in which exists a zero-crossing point at a certain frequency. The dynamic negative mass density is induced by the dipolar mode [42] and can be generated by the coupling between two or more HRs [43]. However, such coupling is generally difficult to tune. Here, we propose a simple yet effective way to adjust this coupling by introducing a connecting tube between these two identical HRs. By doing so, the first two eigenmodes can be obtained by paralleling two identical HRs, and we then tune the dimension of the connecting tube, by which ρ ≈ 0 can be readily adjusted, to realize DZIM at the same frequency. As listed in figure 1(a), the dimensions of the unit cell are: L = 20.3 mm, W = 18 mm, H = 4 mm, D = 8 mm, t = 2 mm, h = 6 mm, and d = 4 mm. Inside the unit cell, the spacing between the necks of HRs is a = 24 mm. The DZIM unit cell is measured by the Brüel & Kjaer type-4206 impedance tubes as shown in figure 1(b). To form the DZIM cloaking structure, the unit cells are arranged into a 5 × 9 array, with the period of the unit cell being 2a and a in x and y directions, respectively. The array is grafted to a two-dimensional waveguide, in which a 48 mm × 72 mm rectangle block is centrally embedded for the demonstration of cloaking. The inner dimension of the waveguide is 588 mm (length), 198 mm (width), and 29 mm (height) with its wall thickness of 2 mm. The resonant frequency of the HR is 1600 Hz, which is below the cutoff frequency of the waveguide, so as to ensure a plane wave propagation inside the waveguide.

Model analysis
The constitutive parameters (ρ and B) can be deduced from the Green's function of the metamaterial surface response [43]. At the subwavelength scale, we apply a homogenization scheme based on the eigenmode analysis for this unit cell in the direction of wave propagation. If no damping is considered, the Green's function can be expressed by eigenfunction as [44]: where ω = 2πf is the angular frequency, ω n is the nth eigenfrequency, ϕ n (x) and ϕ * n (x ) are the eigenfunction of the nth eigenmode at the unit cell's surface positions x and x , respectively.
For scalar waves, the Green's function of a homogeneous one-dimensional system G αβ is related to the heterogeneous surface response G αβ [44], as G αβ = G αβ (Z, k), where Z is the effective impedance and k is the effective wave number. The combined α and β denotes the state of the heterogeneous surface response (G) of the unit cell. For simplicity, we will only consider the lowest two surface modes: the monopolar case with α(β) = 0 and the dipolar case with α(β) = 1. In parity symmetric cases, except for the two trivial ones G 10 = G 01 = 0, the unit cell has only two independent surface modes: G 00 and G 11 , representing the monopolar mode and the dipolar mode, respectively. With the accompanying boundary conditions and the wave equations [44,45], the unit cell's surface responses are described as Here, the unit cell is grafted to a cylindrical waveguide with a length of 2a, that is x ∈ (−a, a). ϕ n (a) and ϕ n (−a) are obtained by the eigenfrequency study in COMSOL Multiphysics, where n = 500 is set in the calculation for the hybridization and maintain the numerical accuracy of the dipolar and monopolar modes. G 00 and G 11 , given by equation (2), lead to the following expressions: Based on the relationship between Green's function and constitutive parameters, the calculated G 00 and G 11 of the unit cell with different diameters of the connecting tube are shown in figure 2. There are three frequencies in the Green's function plot, named f 1 , f 2 and f 3 , in which f 1 and f 3 sit in a relatively narrow range. As illustrated in figure 2 and derived in equation (3), G 11 abruptly transits to infinity at f 1 where ρ becomes divergent, and crosses zero at f 3 where ρ crosses zero. G 00 transits to infinity at f 2 , where 1/B crosses zero. As a result, the double near-zero region can be found near f 2 and f 3 . It is further found that G 11 is largely dependent upon the dimension of the connecting tube, and can be readily tuned. For instance, as in figure 2(a) where d = 4 mm, the frequency range between f 2 and f 3 is 11.8 Hz. With the diameter decreasing to d = 0.4 mm or increasing to d = 6 mm ( figure 2(b)), G 11 can be effectively shifted to the lower or higher frequencies, respectively. The inset of figure 2(a) illustrates that the dipolar mode induces a negative effective mass density at the eigenfrequency of f 1 = 1773.2 Hz. In the meantime, the monopolar mode induces a negative effective bulk modulus which is tuned to 1/B ≈ 0 around the eigenfrequency of f 2 = 1800.2 Hz. Thus, within the overlapping region of the dipolar mode and monopolar mode, there exists an optimal diameter of the connecting tube with which ρ and 1/B can be engineered to simultaneously cross zero from negativity to positivity at the same frequency, which also ensures the impedance matching with the surrounding medium. To better understand the role of the connecting tube, the effective medium theory [46] is applied to retrieve the constitutive parameters of the unit cell and particular attention is paid to the frequency regions where ρ ≈ 0 and 1/B ≈ 0. As shown in figure 3, 1/B ≈ 0 occurs in the frequency region of 1740-1802 Hz, with the bandwidth of 56 Hz. As for ρ ≈ 0, the corresponding frequency moves monotonically towards to the higher end by increasing the diameter of the connecting tube d and decreasing its height h, with the bandwidth of 1175 Hz (1455-2630 Hz). It is clear that ρ is much more sensitive than 1/B, indicating that the dimensional tuning of the connecting tube is effective in achieving DZIM. Its mechanism is as follows: (1) the monopolar mode, corresponding to the infinite point of G 00 and 1/B ≈ 0, is generated by the same interaction force between air inside the unit cell and the incident acoustic pressure field. The acting force is produced by the incident acoustic pressure field at the neck of the HRs. The reaction force produced by air is related to the volume of each HR (V single = 6570.72 mm 3 ) and the connecting tube (V max = 580.78 mm 3 ). Although the connecting tube participates the air compression-expansion process, its volume is only 1/24 of the whole volume of the unit cell, thus has little influence to the whole. Therefore, the corresponding frequency of 1/B ≈ 0 is less sensitive to the varying dimensional parameters of the connecting tube. (2) The dipolar mode, corresponding to the infinite point of the G 11 where near ρ ≈ 0, creates a negative response such that the motion of the center of mass density is out of phase with the incident pressure field. Inside the unit cell, the dipolar mode is mainly generated by the coupling between two HRs, in which the interaction path of air within two identical HRs has been changed due to the existence of the connecting tube. Therefore, the dynamic balance of mass density is greatly influenced by the diameter or height of the connecting tube. The DZIM can also be achieved by various dimensions of the connecting tube, as illustrated by the DZ line of the constitutive parameters in figure 3. Interestingly, this DZ line indicates that with proper dimensions of the connecting tube, we can also readily realize the DZIM unit cell at different frequencies.

Numerical and experimental validation of DZIM unit cells
The unit cell and array are designed and solved by the frequency-domain study in the commercial finite element method package (COMSOL Multiphysics). Considering the unavoidable viscous losses in experiments, each composite part of the DZIM unit cell is modelled with reasonable losses in the simulations. The eigenstates of the DZIM unit cell are calculated by the eigenfrequency study in COMSOL Multiphysics. Furthermore, the DZIM unit cells are fabricated by the standard 3D printing using the photosensitive resin and measured by the Brüel & Kjaer type 4206 impedance tube apparatus, as illustrated in the figure 1(b).
The comparisons between the theoretical, simulated, measured physical properties of the unit cell are shown in figure 4 with the key dimensions of the connecting tube: d = 4 mm and h = 6 mm. Considering the unavoidable viscous losses and the 3D printing error (±0.2 mm) in our experiments, the measurement results agree remarkably well with the theoretical predictions and simulations in terms of ρ and 1/B. Comparing figures 4(a) and (b) with figure 2(b), ρ jumps from positivity to negative infinity at f 1 , where G 11 becomes divergent; then ρ crosses zero at f 3 , and turns to positivity at high frequencies. Similarly, 1/B crosses zero at f 2 where G 00 becomes divergent. Therefore, there is a frequency region in which ρ and 1/B are both negative, and this makes it possible for them crossing zero simultaneously.
The near-zero region of ρ and 1/B retrieved from simulations and experiments is shown in the inset of figure 4(b). The values of ρ and 1/B that are closest to zero in the simulations occur at 1797 Hz, where ρ ≈ −0.008, 1/B ≈ −0.005 (green arrow). This crossover point shifts to 1803 Hz in the experiment, where ρ ≈ 1/B ≈ −0.12 (blue arrow). At 1810 Hz, ρ ≈ −0.002, and 1/B ≈ −0.102 (cyan arrow). It is interesting to see that, ρ is more sensitive to frequency, indicating the tuning of the connecting tube is vital in designing this type of DZIM. The widened frequency range in the measurement is mainly due to the viscous losses and the printing errors, which may be in favor in designing broadband DZIM. However, the deteriorated transmission is unavoidable and must be carefully balanced.
From figure 4(c), we can see the theoretical transmission coefficient matches very well with the simulation and measured results. The transmission has two dips (f , f 1 ) and one peak (f 3 ), corresponding to the two divergent points and one zero point of the phase, respectively ( figure 4(d)). The transmission properties are determined by three factors: (i) the local resonance (f ) of a single HR which is precisely where G 00 equals G 11 ( figure 2(a)), (ii) the dipolar resonance (f 1 ) inducing the divergent ρ which can be effectively tuned by the coupling between the two HRs via the connecting tube, and (iii) the effective impedance (Z = √ ρB) that matches well with the background air around the transmission peak (f 3 ) where both ρ and 1/B are near zero (figure 5). Note that when the viscous losses are considered in the simulations, the transmission can no longer reach unity but 0.72 at f 3 and the measured transmission is 0.602 as losses may be much larger than those set in the simulations. By inverting the transmission and reflection coefficients [46], we calculate the corresponding impedance (Z) and the refractive index (n). Then the effective mass density (ρ) and the reciprocal of the bulk modulus (1/B) are obtain fromZ = √ ρB and n = ρ/B. Here, these effective parameters of the lossless case are displayed in figure 5. When the real parts of the constitutive parameters (ρ and 1/B) are in the white region (1780-1806.8 Hz), corresponding to the negative-index behavior (i.e., real(n) < 0, imag(n) ≈ 0), the phase velocity and group velocity in the unit cell are in opposite directions. In the light-yellow region (1806.8-1814 Hz), the real parts of ρ and 1/B are opposite in sign, and they are both near zero, i.e., −0.2846 < real(ρ) < −0.0038, and 0.0023 < real(1/B) < 0.0253. While the real(n) is very close to zero (−0.0096 < real(n) < −0.0016) in this region, the refractive index n is purely imaginary. Meanwhile, the real(Z) is much closer to zero with the non-zero imag(Z), indicating that the acoustic waves in this unit cell are evanescent. It is interesting to see that the unit cell achieves DZ-index at 1815 Hz, where the real parts of the constitutive parameters (ρ and 1/B) are both very close to zero, i.e., real(ρ) = 0.027, real (1/B) = 0.0287. The corresponding real(n) is 0.0279 and the real(Z) is approaching unity (0.9696), which means the unit cell matches well with the background air. In terms of the transmission coefficient (T ), it is important to note that our Helmholtz-type DZIM unit cell exhibits a rather high transmission (∼1) in the defined frequency range. In addition, the double-near-zero band is demonstrated in the light-green region (1814.2-1819.4 Hz) where the optimum point occurs at 1815 Hz. After this frequency band, the real parts of the constitutive parameters depart from zero, resulting in that the real(n) is no longer near zero and the impedance matching is deteriorating.

Acoustic cloaking of the DZIM array
To demonstrate the cloaking effect, we conduct the array experiment to acoustically hide a rectangle block at the center of a two-dimensional waveguide, as shown in figure 6(a). The DZIM array is also fabricated by the standard 3D printing and measured by an acoustic sweep device. Six speakers connected to the power amplifier (Brüel & Kjaer type-2734) are placed on the left side of the rectangular waveguide to mimic a  plane wave source figure 1(a). The measured frequency range is 1200-2200 Hz. The free field 1/8-inch microphone (Brüel & Kjaer type-2670) with the diameter of 1 cm is inserted into the waveguide to measure the outgoing field amplitude and phase in the test area in figure 6(a) by stepping with 0.18 mm × 0.18 mm increments in the x-y plane. The results are recorded by the data-acquisition hardware (NI PXI-4461 in NI PXIe-1071) and postprocessed by the NI LabVIEW software. As displayed in figure 6(a), there is a 18 × 7 array of the interspaced cylindrical holes located atop the waveguide on the left side of the block. The compound sound pressure fields of the incident wave and the reflected wave are measured in these holes. The other holes at the corners of the waveguide are designed to calibrate the position of the microphone inserted into the waveguide. The 5 × 9 array of the DZIM unit cells is grafted to the bottom of the waveguide as shown in figure 6(b).
Figures 6(c) and (d) display the pressure field distribution inside the waveguide with and without the DZIM array under a plane wave excitation, respectively. We observe that only a minor phase change on either side of the DZIM array in figure 6(c), as if the rectangle block is absent. Within the test area, our measured pressure fields at 1900 Hz (figure 6(e)) maintain a perfect plane waveform despite the edge diffraction in the waveguide, which agree remarkably well with the simulation results. The frequency shifts to a higher frequency, which is due to the viscous losses and the 3D printing errors. As the higher the frequency corresponds to the shorter the wavelength, the whole waveform in the test area slightly shifts to the left while the phase characteristics keep consistent. For comparison, we removed the DZIM array and place the block inside the waveguide, the compound pressure fields on the left side of the block and the outgoing pressure fields lose the uniform phase distribution and display a classic diffractive pattern, which can be observed in figures 6(d) and (f), indicating that the acoustic cloaking effect no longer exists.
Finally, the measured pressure fields with the DZIM array are recorded at different frequencies (figure 6(g)) and mapped onto the dashed area in figure 6(c). Visually, the uniform phase distribution is gradually established from 1865 Hz, reaching the optimum at 1900 Hz, and then deteriorates after 1925 Hz. A relatively wideband DZIM array with the bandwidth of 60 Hz is experimentally demonstrated.

Conclusions
In this study, we have designed a very simple acoustic DZIM structure consisting of two coupled HRs via a connecting channel, in which DZ in both the effective mass density and the reciprocal of bulk modulus can be readily designed and characterized. The dipolar mode that can be independently adjusted by fine-tuning the size of the connecting tube, while the monopolar mode that is virtually unchanged. Besides, we prove that various dimensions of the connecting tube can be utilized for our DZIM design. Acoustic cloaking is numerically and experimentally demonstrated by our designed DZIM array, and more importantly, the cloaking bandwidth of 60 Hz is successfully demonstrated in the experiments. It is observed that although the transmission efficiency is inevitably hampered by the viscous losses within the unit cells, the viscous losses within the connecting channel between adjacent HRs would flatten the near-zero ρ curve, which as a result can effectively broaden the near-zero-index bandwidth of the DZIM unit cell as well as its acoustic cloaking. This is rooted from the local resonance mechanism of our coupled HRs. Since the DZ-index properties depend on the local resonance of the unit cell, the DZIM array with different numbers or arrangement of the unit cells can be designed for DZ-index applications. The simplicity in the unit cell design makes it a suitable building block for the potential acoustic metamaterial devices with complex functionalities, such as wavefront manipulation, and phase matching.