A transfer matrix method for calculating the transmission and reflection coefficient of labyrinthine metamaterials Open Access

Physicists and engineers have been interested in controlling, directing, and manipulating sound waves for a long time. Traditionally, phased arrays have been used to generate these functionalities.^1–4 However, in these devices, the distribution of their sources in an active way needs to be driven individually using electrical techniques, leading to major drawbacks such as high cost and complexity. Acoustic metamaterials offer new possibilities to achieve these functionalities. They are man-made, specially fabricated materials featuring properties never found in nature. Such metamaterials usually gain their properties from their structure rather than their constitutive materials, through techniques such as the controlled fabrication of small inhomogeneities to achieve effective macroscopic behaviour.^5–8 These innovative engineered materials, emerging recently with rich physics and plenty of extraordinary capabilities in various wave manipulations, have intrigued scientists and engineers for more than ten years now, paving the way for many applications of acoustic and elastic wave engineering.^9–13 

As a subcategory of acoustic metamaterials, acoustic metasurfaces demonstrate the advantageous features of a planar profile and subwavelength thickness compared to bulky meta-structures, and, thus, possess additional capabilities and flexibilities in sound wave manipulation. Metasurfaces have been demonstrated to be capable of many forms of wave manipulation and, thus, have drawn significant attention from both the physics and engineering communities.^14–23 Many geometries have been used by researchers to manufacture acoustic metasurfaces such as labyrinthine structures,^14,17,18 helical structures,²⁰ multi-splits,^19,21 and Helmholtz resonators.^15,16,22,23

Labyrinthine structures have recently attracted extensive research attention due to their ability to exhibit high refractive indexes, multiple vibration modes, and extraordinary acoustic properties. These structures enable the effective speed of acoustic waves to be slowed down due to path elongation by means of folded narrow channels and may be suitable for achieving broadband control of acoustic waves with a single structural unit.

In this work, we develop an analytical model based on the transfer matrix method for the phase shift of transmission/reflection labyrinthine unit cells. The labyrinthine structure channel has been divided into two kinds of short pipes, as proposed by Yong Li et al.¹⁷ However, our analytical development led us to a different matrix. We present the amplitude and phase shift for two different labyrinthine structures, full reflection, and transmission, in the frequency range of 830–4630 Hz, which covers the peak sensitivity of human hearing. We validate our model with finite element method (FEM) simulations—using commercial software—and measurements on transmissive and full reflective units with 4 and 6 baffles perpendicular to the direction of propagation (“bars” in the following) using different experimental methods. In particular, we used point scans (in the frequency domain) and impulse response (in the time domain, using square pulses) for the transmissive units. For the reflective ones, we focused on the time domain, looking at impulse response using as sources either white noise, a maximum-length sequence (MLS) or a square pulse. The model presented here can serve as an efficient tool for designing an acoustic metasurface that provides the desired precision in the acoustic field.

We consider a specific labyrinthine structure for full reflection with a width of $ax$⁠, height $ay$⁠, and thickness $t$⁠, and consisting of $m$ and $n$ identical bars in the upper and lower boundary, respectively. The space of the zigzag channel can be controlled by changing these parameters $(m,n)$⁠. The width and length of these identical bars are $w$ and $l$⁠, respectively. The labyrinthine structure channel width can be expressed as $d=ax/(m+n)−w$ and the length of the bars $=ay−2t−d$⁠. Figure 1 shows a schematic of this specific labyrinthine structure unit for $m=2$ and $n=2$⁠.

The labyrinthine channel is divided into two types of regions. Those of the first type—A, C, E, and G in Fig. 1—have a geometric size along the $y$ axis that is much smaller than the working wavelength. Therefore, it is possible in these regions to simplify the analytical derivation by considering plane wave propagations only. Those of the second type—B, D, F, and H—have geometric sizes comparable to $ay$⁠. Therefore, the waves propagating in these regions can be expressed as a summation of the normal modes. In our model, we consider an incident plane wave travelling in the opposite direction to the $x$ axis and the reflected wave travelling in the direction of the $x$ axis, as shown in Fig. 1. The pressure, $p$⁠, and the velocity component, $u$⁠, are connected via the relationship $ik0ρ0c0u=∂p/∂x$⁠, where $k0$⁠, $ρ0$⁠, and $c0$ are wave number, density, and velocity in air, respectively. The harmonic factor $e−iωt$ with $ω=k0c0$ is omitted.

As mentioned previously, the pressure field and the velocity component in the regions $A (−w≤x≤0and t≤y≤t+d)$⁠, and $C (−2w−d≤x≤−w−d and t+l≤y≤t+l+d)$ can be expressed as a plane wave only,

The pressure and velocity component in region B $(−w−d≤x≤−w and t≤y≤t+d+l)$ can be expressed in terms of the normal modes as

where $∅n(y)$ is the transverse eigenmode and can be expressed as $∅ny=2−δ0n cos kyn(y−t)$⁠, and satisfies the orthogonality $∫∅n(y)∅m(y)dy=σδmn$⁠, with $δmn$ representing the Kronecker delta and σ being the cross section of the channel. Here, the wave number along the $y$ direction can be expressed as $kyn=nπ/(l+d)$ (n = 0, 1, 2…) and $kxn=k02−kyn2$⁠. The symbols $A+, A−, C+, and C−$ denote the coefficients in the thin channel, and $Bn+ and Bn−$ denote the coefficients of the nth modes. The signs (+) and (⁠$−$⁠) refer to the propagating coefficients in the direction of and opposite direction of the x axis, respectively.

The total pressure field in the incident region $(x≥0and 0≤y≤ay)$ can be expressed as

where $ψny=2−δ0n cos(kyn′y$⁠), is the transverse eigenmode, with the wave number along the $y$ direction, $kyn′=nπ/l+d+2t$ and $kxn′=k02−kyn′2$⁠. In Eq. (7), $Rn$ represents the reflection coefficient of the nth mode.

To take into account the potential effect of viscous and thermal losses, an additional term was added to the wave number along the $x$ direction, as expressed in Ref. 24,

and

where

Here, $γ=1.4$ is the adiabatic constant for air. In standard conditions, the speed of sound in air $c≈343 m/s$⁠, the viscous characteristic length $lv=4.5×10−8m$ and the thermal characteristic length $lt=6.2×10−8m$⁠.

Now, considering the continuity of pressure and volume velocity at the interface between region A and B at $x=−w$ and between regions B and C at $x=−w−d$⁠, we obtain the following equations:

Before the next step, we multiply Eqs. (15) and (17) by $∅m(y)$ and integrate both sides of the equations along the boundaries in the $y$ direction. After utilizing the orthogonality of $∅m(y)$⁠, the coefficients $Bn+$ and $Bn−$ can be obtained as below,

Substituting the coefficients $Bn+$ and $Bn−$ into Eqs. (14) and (16), and averaging the pressure field, we obtain the following equations:

where

To obtain $C−$⁠, we subtracted Eq. (20) from Eq. (21). To obtain $C+$⁠, we rearranged Eq. (20) before subtracting it from Eq. (21). In this way, a transfer matrix could be constructed to connect the coefficients between regions A and C as follows:

where $M1$ is the transfer matrix

$Since Φn,12≈Φn,22$⁠, we assume that $α=α1$⁠. Therefore, the transfer matrix $M1$ reduces to the following:

The transfer matrix from regions C to E is the same as the matrix $M1$⁠. Therefore, the coefficients in region G can be expressed as

Considering the continuity of pressure and volume velocity at the interface between region G and H at $x=−4w−3d$⁠, and following the same procedure as before, we obtain the following equation:

Since there is a hard boundary, the velocity at $x˘=−4w−4d$ should be zero, which in the form of an equation becomes $uHx˘,y=0$⁠. Therefore, the relation between the coefficients $Jn+$ and $Jn−$ can be obtained as

Combining Eqs. (29) and (30), the coefficients $Hn+$ and $Hn−$ can be expressed as

Substituting the coefficients $Hn+$ and $Hn−$ into the pressure continuity equation between region G and H at $x=−4w−3d$⁠, and averaging the pressure field, the coefficients in region G can be expressed as

Considering the continuity of pressure and volume velocity between the incident region and region A at $x=0$⁠, and following the same procedure as above, a transfer matrix can be constructed to connect the incident region coefficient and region A coefficients as follows:

where

where

Finally, by combining Eq. (28), Eq. (33), and Eq. (34), the reflection coefficient can be expressed as

with

For a general case of a labyrinthine structure with $m$ and $n$ identical bars, the matrix M can be written as

In this section, we consider a labyrinthine structure unit for transmission and reflection, consisting of $m$ and $n$ identical baffles, perpendicular to the direction of propagation (“bars” in the following) in the upper and lower boundary, respectively. To describe the geometrical parameters $ax, ay, t, w$ in Fig. 2, we will use the definitions provided in the previous section. For this structure, however, we will use a different expression for the channel width d, which is given by $d=(ax/m+n−1)−(m+n/m+n−1)w$⁠. Figure 2 shows a schematic of a labyrinthine structure unit for transmission and reflection with $m=2$ and $n=2$⁠.

The labyrinthine structure channel is divided into two regions, as described for the full reflection structure, which means that only the plane wave is considered for the thin channels and a summation of normal mode is employed for the long channels. Therefore, the pressure field and velocity component in these regions, along with the incident region, can be expressed as defined above. For the transmitted region, (⁠$x≤−4w−3d$ and $0≤y≤ay$⁠), the pressure field can be expressed as the following:

where $Tn$ represents the transmission coefficient of the nth mode.

Following the same procedure for the full reflection case, the transfer matrix connecting the coefficients between region K and the incident region can be expressed as

where

Considering the continuity of pressure and volume velocity between region G and the transmitted region at $x=−4w−3d$⁠, the coefficients in region G can be expressed as the following:

and

where

Combining Eq. (41) with Eq. (44), the reflection coefficient can be written as

Substituting the reflection coefficient R into Eq. (41) and combining with Eq. (43), the transmission coefficient can be written as

Numerical simulations were carried out by the finite element solver in the commercial software COMSOL Multiphysics, using its Acoustics module. Simulations featured a three-dimensional (3D) model of each single metamaterial structure inserted in a cuboid waveguide of the same base size, terminating with a perfectly matching layer. A hard-wall boundary condition was associated with the solid's boundaries in these simulations, while periodic conditions were used at the sides of the waveguide. The incident wave entered the metamaterial structure from the side opposite to the perfectly matching layer with unitary amplitude and zero phase. A study in the frequency domain was made in the selected frequency range of 830–4630 Hz, and the total acoustic pressure field and sound pressure level resulting in the output at a specific checkpoint were recorded. Then, the transmitted and reflected pressure field were obtained by subtracting two pressure fields with and without the metamaterial.

Metasurface samples consisting of 10 × 10 labyrinthine structure units with 4 bars (m = 2 and n = 2) and 6 bars (m = 3 and n = 3), were fabricated from thermoplastics (polylactic acid, also known as PLA), using a 3D printer (Stratasys F170). The PLA surface finish can be defined by the layer resolution, which is $250 μm$⁠. The labyrinthine structure units have the same geometrical parameters of the structure unit used in the model and simulation. Therefore, the metasurface samples have a height of 25 mm and a base of 125 × 145 mm. Figure 3 shows a photograph of the metasurface sample for transmission. This procedure was repeated both for the reflective and the transmissive units.

A commercial loudspeaker (JVC CS-J520X) for cars was used as a source. Its emission has previously been characterized outdoors in the range of interest and found to be equivalent to a standard reference source. The measurements were conducted in a room lined with absorbing material, but in order to further minimize spurious reflections, the loudspeaker was inserted into a circular plastic waveguide (diameter, 161 mm) to channel the waves on the sample. The length of the tube (270 mm) ensured that fully developed plane waves reached the metasurface sample.

For the transmissive case [Fig. 4(a)], the metasurface sample was fixed to one of its ends, with the loudspeaker attached to the other end (see Fig. 4). Two experimental measurements were carried out for measuring the transmission coefficient, a point scan (in the frequency domain) from 800 to 5000 Hz with a step of 200 Hz²⁵ and an impulse-response method (in the time domain), using a 100 μs square pulse as an excitation signal for generating an impulse response. For the point scan, measurements were taken using an omni-directional B&K 4138 microphone, designed for measuring in confined spaces with a very high frequency response of 6 Hz–140 kHz, and connected to a conditional amplifier (B&K NEXUS 2690–0S2). The microphone was positioned in front of the metasurface sample at 100 mm. For the time domain measurements, we used one omni-directional microphone (1/2 in. free field, model Norsonic, type 1201/30323) placed in front of the metasurface sample at 45 mm. The measurements were carried out with and without the metasurface in three different environments: a small room, whose walls were covered with absorbing material, on a flat surface outdoors, and in a large university room. Received signals were acquired using a 4-channel Picoscope of the 2000 series and acquisition software developed for the purpose. The arbitrary wavefront generator (AWG) of the Picoscope was used to generate the source signals.

For the reflective case [Fig. 4(b)], a two-microphone transfer function method inspired by impedance tubes was used to measure the reflection coefficient.²⁶ The metasurface was fixed at 157 mm from the end of the circular waveguide. Two nominally identical microphones (“Mic1” and “Mic2,” ½ in. free-field, model Norsonic, type 1201/30323) were mounted between the waveguide and the metasurface, with the first positioned from the metasurface (L) and the second from the first microphone (S). Data were acquired as in the transmissive case, and a program compiled in matlab was developed to compute the transfer function between the two microphone positions. The AWG of the Picoscope 2200 b was used to generate the three different source signals for these measurements: a square pulse (100 s), white noise, and a 10-bit MLS.

For both cases—full reflection and transmission—the labyrinthine structure unit dimension was fixed to be $ax=25 mm$ and $ay=12.5 mm$ to ensure the subwavelength property of the structure. The bar width, $w$⁠, and wall structure thickness, $t$⁠, were fixed to be 1 mm to facilitate the labyrinthine structure unit fabrication. A program compiled in matlab was developed to calculate the theoretical transmission and reflection coefficient (amplitude and phase shift) predicted by Eq. (37), Eq. (46), and Eq. (47). The transmission and reflection coefficient with losses, obtained by replacing $kxn$⁠, $k′xn$ and $k0$ with $kxn_Th$⁠, $k′xn_Th$ and $kTh$⁠, respectively in Eq. (37), Eq. (46), and Eq. (47). To validate our theoretical prediction, commercial software based on the FEM—COMSOL Multiphysics version 5.3 plus the Acoustic module—was employed for the simulations, using the same geometrical parameters of the analytical model.

First, we present the comparison between our analytical model with and without the effect of viscous and thermal losses and COMSOL predictions. Figure 5 shows the amplitude and phase shift of the analytical reflection coefficient, with (green cicle marker) and without losses (blue square marker) and simulated (red cross-marker), at the frequency range of 830–4630 Hz, for the closed labyrinthine cells, with 4 bars (⁠$m=2$ and $n=2$⁠) and 6 bars (⁠$m=3$ and $n=3$⁠). It can be seen from Figs. 5(c) and 5(d) that there is no appreciable difference between the phase shift of the reflection coefficient of the analytical model with and without losses. A slight difference, lower than 0.1 dB, can be observed in the amplitude [Figs. 5(a) and 5(b)] from the frequency 2000 to 4630 Hz for 4 bars and in the frequency range [830–1430 Hz] and [3830–4630 Hz] for 6 bars. As the width, d, of the channel is large compared to the thermal and viscous wavelengths (d = 5.25 mm for 4 bars and 3.16 mm for 6 bars in the reflective case), it is not surprising that thermoviscous effects are neglegible. In terms of amplitude [Figs. 5(a) and 5(b)], the difference between the reflection coefficients predicted by the analytical model without losses and by the COMSOL simulation barely exceeds 0.001 dB. In terms of phase shift [Figs. 5(c) and 5(d)], there is a maximum difference of 1° between the phase shift predicted by the simulations and the analytical results, for both the 4- and 6-bar structures.

Negligible effects of the thermo visous losses on the transmission coefficient were observed also for the transmissive (“open”) unit cells, both in amplitude and phase (see Fig. 6). Even in this case, the width of the channel is in fact much larger than the scale of these effects (d = 7 mm for the 4-bars structure and 3.8 mm for the 6-bars metamaterial). For the transmissive (“open”) unit cells, a good agreement on the predicted phase shift was also observed between the theory and simulation, as shown in Figs. 6(c) (4 bars) and 6(d) (6 bars). The blue square and red cross-markers, shown in this figure, represent the theoretical without losses and simulated results, respectively. For 4 bars, the maximum value of the difference in the explored range of frequencies (830–4630 Hz) is 2.6 degrees. Above 3030 Hz, the difference between the analytical and simulation results is even smaller, as it does not exceed 0.7 degrees. For 6 bars, the maximum value of the difference is 4.7 degrees and is observed at frequency 2830 Hz.

The amplitude of the simulation results for both structures—4 and 6 bars—shows the same behaviour as the analytical results. A small difference between the simulation and the theoretical results is observed, and the maximum value of the difference is 0.3 dB and 0.6 dB for 4 and 6 bars, respectively; see Figs. 6(a) and 6(b).

After benchmarking it with FEM simulations, a measurement campaign was run to validate our model. Metasurface samples consisting of 10 × 10 labyrinthine structure units that have the same geometrical parameters of the model and simulation, with 4 and 6 bars, for full reflection and transmission, were designed, manufactured, and tested over a frequency range of 830–4630 Hz.

First, we tested the transmittive case. Experimental measurements of the amplitude and phase shift of the transmission coefficient were carried out using point scan (in the frequency domain) and impulse response (in the time domain) for structures made of 4 and 6 bars. In the absence of a fully anechoic environment, the potential influence of room reflections on the results was averaged out by repeating the measurements in different rooms and outdoors. The experimental results (averaged over all the experimental rooms and with one standard deviation uncertainties) are compared in Fig. 7 with the analytical model.

Figure 7 shows that, for both 4 and 6 bars, the analytical theory predicts the amplitude results within $±2$ dB for both point scans and the pulse scans. The difference between the two measurement methods, however, becomes apparent in the phase shift results, where the pulse scans are closer to the analytical model. We attributed the difference to potential local effects, a hypothesis that was further investigated in the reflective case, using two microphones.

The two-microphone transfer function was used to determine the properties of the reflective metamaterial structures with either 4 or 6 bars. These measurements proved more challenging since the results could easily be affected by unwanted reflections (from the room) and by diffraction effects from the border of the sample (which was $145×125$ mm in size, with the wavelength of the sound ranging from 68 to 480 mm in the range of frequencies explored) but the method gave greater insight on the sources of uncertainties. The average results with statistical uncertainties (averaged for each method over all the experimental rooms, with one standard deviation) are shown in Fig. 8, in comparison with the predictions for the analytical model.

It can be seen from Figs. 8(a) and 8(b) that the theory overestimates the measurements taken with the pulse method by an average of 3 dB. Also, apparent in the results with this method are two peaks, one at approximately 880 Hz (⁠$≈390 mm$⁠) and one at approximately 2180 Hz (⁠$≈157 mm$⁠), which are due to reflections in the measurement set-up that were not windowed out. In our opinion, these uncontrolled reflections influenced the point scans in Fig. 7, since these covered a much larger area in front of the metamaterial. These two peaks disappear when the signal used is continuous: measurements were within 2 dB from the theory when using a MLS or white noise as a source [Figs. 8(c) and 8(d)].

Figures 9–11 show the phase shift for 4 and 6 bars for the three type of measuremts, pulse, noise, and MLS, respectively.

In the case of 4 bars, the pulse measurement [Fig. 9(a)] already show good agreement between theory and experimental results whatever the room, i.e., in the large room (four repeats), outdoor (one repeat), and in the small room with absorbing material on the walls (four repeats). Also, good agreement can be observed for 6 bars [Fig. 9(b)], where in most of the range the theory falls within the statistical uncertainty on the measurements. The only exception is outdoors, where the measurements deviate from the theory below 1400 Hz (i.e., where the amplitude peak due to spurious reflections was observed) and in the range of frequency 3490–4000 Hz. The deviation in the latter range was attributed to diffraction effects from the border of the sample, which was less than two wavelengths wide already at these frequencies. Figure 9(a) supports this hypothesis by showing the response (in phase) of a fully reflective plate: the impinging wave sees the metasurfaces as fully reflecting up to 2000 Hz (probably due to the external baffle/bar covering most of the surface at these wavelenghts), with the labyrinth-controleld delay picking up only at higher frequencies. Incidentally, at 2000 Hz the metasurface is approximately $λ2000/2π$ thick.

Like in the case of amplitude, the other signals led to a better agreement between theory and experiments. In the case of white noise, results of 4 bars [Fig. 10(a)], show a good agreement with the theory except for the frequency below 1500 Hz. Experimental results of MLS show a trend similar to the one predicted by the analitical method, but results for 4 bars appear shifted by $π/2$ [Fig. 11(a)].

For 6 bars, the analytical results, while white noise in the large room shows some agreement with the values predicted by the theory [Fig. 10(b)], the results obtained in the small room or using MLS are very far from the theory [Fig. 11(b)]. We could not determine the reasons behind this difference in this study, as they stayed in all the experimental conditions we used: they did not depend on the room or on the sampling or on the processing.

An analytical model based on the transfer matrix method for the phase shift of a full reflection and transmission/reflection labyrinthine structure unit was developed in this study without the effect of viscous and thermal losses. The results show a slight effect on the amplitude of the reflection coefficient and no effect on the transmission coefficent, which is due to the large channel width that we used in this study. A numerical simulation based on the FEM using COMSOL Multiphysics was also developed, and the phase shift agrees well with the analytical prediction. Measurements on transmissive and full reflective units with 4 and 6 bars using different experimental methods were carried out to compare with the theoretical and numerical results. In this case, the agreement with the analytical predictions is within 10% for most frequencies. We found that some of the measurement methods lead to an aagreement within 2 dB, while others are completely out of range, thus, pointing out the challenges in characterizing this type of acoustic metamaterial.

The model presented here can serve as an efficient tool for designing an acoustic metasurface that provides a desired precision in the acoustic field and opens the door for more automated wavefront manipulation and engineering of acoustic waves.

The authors acknowledge funding from the Engineering and Physical Sciences Research Council (EPSRC-UKRI) through project AURORA-EP/S001832/1.