Multiple scattering by scatterers suspended in composites is a fundamental topic with the extensive literature. The estimation of the effective properties of composites is very significant for the design and preparation of composites. So far, many related approaches and predictive models have been proposed. As early as 1945, Foldy et al. proposed the MST for waves that are isotropically scattered by randomly distributed scatterers. Since the specific configuration of the amorphously distributed scatterers is not important, the configuration of the obstacles can be considered as the average of one state in an ensemble. Foldy considered the first order of the scattering coefficient for a single scattering and obtained Eq. 1 to calculate the frequency-dependent complex wavenumber   K , where K 0 is the wavenumber of the matrix phase, and N is the number of the scatterers (Foldy, 1945). Lax generalized Foldy’s results to anisotropic scatterers by inducing a quantum-mechanical viewpoint and using the quasi-crystalline approximation (Lax, 1951; 1952). However, in both Foldy’s and Lax’s studies, the effective wavenumber is expressed in terms N and the cumulative effect of the only forward far-field scattering amplitude, without considering the important role of the backscattered amplitude. K 2 = K 0 2 − i 4 N f ( θ ) . (1)

Waterman and Truell further improved Foldy’s method by introducing quadratic terms for the number of scatterers, laying a foundation for other multiple scattering theories. Their most important result is that the wave propagation in the scattering medium can be described by the far-field amplitude f ( θ ) of a single scatterer, which is an essential parameter of the scattering medium. As observed from Eq. 2, the Waterman and Truell framework has taken both forward ( θ = 0 ) and backscattered ( θ = π ) far-field scattering amplitudes into consideration (Waterman and Truell, 1961). ( K / K 0 ) 2 = 1 − 4 i N f ( 0 ) / K 0 2 − 4 N 2 [ f ( 0 ) 2 − f ( π ) 2 ] / K 0 4 , (2) f ( θ ) = ∑ n = 0 ∞ ε n T n cos ( n θ ) , (3) where ε n = 1 for   n = 0 and ε n = 2 for n ≥ 1 . The expansion coefficients T n are determined by the boundary conditions of the particular problem considered. By comparing Eq. 1 and Eq. 2, it is noted that when the far-field amplitude is equal to the forward far-field amplitude in the Waterman and Truell framework, Eq. 2 is simplified to Eq. 1 in Foldy’s framework. Twersky considered different incidence field angles and expanded the normal incidence formula in the Waterman and Truell framework into a random incidence angle α i n for two dimensions (Twersky, 1962). ( K / K 0 ) 2 = 1 − 4 i N f ( 0 ) / K 0 2 − 4 N 2 c o s 2 α i n [ f ( 0 ) 2 − f ( π − 2 α i n ) 2 ] / K 0 4 , (4) where α i n is the incidence angle, and f ( 0 ) and f ( π − 2 α i n ) can be obtained from Eq. 3. If the normal incident to the slab ( α i n = 0 ) is considered, Eq. 4 in the Twersky framework is equal to Eq. 2 in the Waterman and Truell framework. Lloyd and Berry (Lloyd and Berry, 1967) further extended and corrected the work by Waterman and Truell, and they showed that the formula for calculating the wavenumber should be ( K / K 0 ) 2 = 1 − 4 i N f ( 0 ) / K 0 2 + 4 N 2 { f ( 0 ) 2 − f ( π ) 2 + ∫ 0 π 1 / sin ( θ / 2 ) d d θ [ f ( θ ) 2 ] d θ } / K 0 4 . (5)

Eq. 5 is different from Eq. 2 by the Waterman and Truell framework, both in the signs of the forward and backscattered terms and in the additional modified term of N 2 . Linton and Martin also presented a two-dimensional model for circular cylinders inspired and based on the hole correction, assuming that the exclusion distance is small compared to the wavelength (Linton and Martin, 2005): ( K / K 0 ) 2 = 1 − 4 i N f ( 0 ) / K 0 2 + ( 8 N 2 / π K 0 4 ) ∫ 0 π cot ( θ / 2 ) d d θ [ f ( θ ) 2 ] d θ . (6)

Based on the effective wavenumber in the Waterman and Truell framework, Angel and Aristégui further described the effective medium by the reflection coefficient and the transmission coefficient (Angel and Aristégui, 2005; Aristégui and Angel, 2007). R = − Q e x p ( − i K 0 d ) × [ 1 − e x p ( 2 i K 0 d ) ] /   [ 1 − Q 2 e x p ( 2 i K d ) ] , (7) T = ( 1 − Q 2 ) × exp [ i ( K − K 0 ) d ] /   [ 1 − Q 2 e x p ( 2 i K d ) ] , (8) where Q = ( Z 0 − Z ) / ( Z 0 + Z ) , Z , indicates the effective impedance of the composite, and Z 0 indicates the impedance of the matrix phase; from the formula Z = ρ ω / K , we can obtain Z / Z 0 = ( K 0 / K ) × ( ρ / ρ 0 ) , (9) where ρ / ρ 0 is the effective mass density (EMD) of the composite, and ρ 0 indicates the mass density of the matrix phase. Angel and Aristégui also obtained the EMD in Eq. 10 and effective bulk modulus (EBM) as M / M 0 in Eq.11. ρ / ρ 0 = 1 − 2 i N [ f ( 0 ) − f ( π ) ] / K 0 2 , (10) M / M 0 = [ 1 − 2 i N [ f ( 0 ) + f ( π ) ] / K 0 2 ] − 1 , (11) where M 0 = ρ 0 c 0 2 , c 0 , indicates the sound velocity of the matrix phase. Obtaining Eq. 12 by taking Eq. 10 into Eq. 9, then the transmission and reflection coefficients can be calculated. Z / Z 0 = ( K 0 / K ) × { 1 − 2 i N [ f ( 0 ) − f ( π ) ] / K 0 2 } . (12)

Recently, Luppé and Conoir further obtained the EMD, EBM, and effective impedance as ρ / ρ 0 = 1 − 2 i N ( f ( 0 ) − f ( π ) ) / K 0 2 + ( 4 N 2 / K 0 4 ) { f ( 0 ) 2 / 2 − f ( π ) 2 / 2       +   J ( 0 ) [ f ( π ) − T 0 − f ( 0 ) ] / f ( 0 ) + T 0 I ( π ) [ f ( 0 ) − T 0 ] / f ( 0 ) +   T 0 J ( π ) / f ( 0 )   −   T 0 I ( 0 ) [ 2 f ( π ) − f ( 0 ) − T 0 ] / f ( 0 ) } , (13) M / M 0 = ( ρ / ρ 0 ) ( K 0 2 / K 2 ) , (14) Z / Z 0 = K 0 / K × { 1 − 2 i N ( f ( 0 ) − f ( π ) ) / K 0 2 + ( 4 N 2 / K 0 4 ) { f ( 0 ) 2 / 2 − f ( π ) 2 / 2   +   J ( 0 ) [ f ( π ) − T 0 − f ( 0 ) ] / f ( 0 ) + T 0 I ( π ) [ f ( 0 ) − T 0 ] / f ( 0 )     +   T 0 J ( π ) / f ( 0 )   −   T 0 I ( 0 ) [ 2 f ( π ) − f ( 0 ) − T 0 ] / f ( 0 ) } , (15) with J ( α ) = − 1 / ( 2 π ) × ∫ 0 π cot ( θ / 2 ) { d d θ [ f ( θ ) f ( α − θ ) + f ( − θ ) f ( α + θ ) ] } d θ , (16) I ( α ) = − 1 / ( 2 π ) × ∫ 0 π cot ( θ / 2 ) { d d θ [ f ( α − θ ) + f ( α + θ ) ] } d θ . (17)

The effective mass density and effective impedance in the Angel and Aristégui framework and Luppé and Conoir framework are equal when only the first order of N 2 is considered. The reasons for the differences of the aforementioned formulas are worth discussing, and errors will occur when different boundary regions are selected during integration. Table 1 summarizes the assumptions, abbreviations, and calculation formulas of the previously mentioned five MST models.

Multiple scattering signals are composed of coherent parts and incoherent parts. Since enough scatterer position configurations are averaged, the incoherent part vanishes in MST theoretical calculations. In practice, it is not easy to obtain coherent wave characteristics by using only one sample of random composite because it requires a large number of local measurements along the composite to restore the averaging procedure. Therefore, numerical simulations are a good choice.

The numerical simulations were finished by COMSOL Multiphysics software in the pressure acoustic frequency domain. The designed numerical simulation geometric model is shown in Figure 1D. A slab of screen S contains N identical scatterers with same radius a , so the fraction is ϕ = N π a 2 / S . Their positions are selected by a random draw function. The periodicity of the model in the direction perpendicular to sound propagation is simulated by applying periodic boundary conditions on both sides of the screen. Using perfectly matched layer (PML) boundary conditions to simulate the outer boundary of the fluid domain, the incident plane wave is insonified on the left boundary of the model. The reflection coefficient and transmission coefficient of an incident plane wave passing through the screen of scatterers are calculated by using acoustic boundary conditions. We have made 100 simulations to achieve the convergence of the averaged field.

We considered the propagation of coherent waves through steel scatterers randomly immersed in water. Figure 2A and Figure 2B present the modulus of complex transmission and reflection coefficients corresponding to steel scatterers in water, according to MST frameworks (Table 1) and FEM simulations. Both in Figure 2A and Figure 2B, the transmission coefficients obtained by four MST frameworks (WT, LB, LM, and LC) agree with each other, except the Foldy framework (F) as expected. It indicates that Foldy only considered the forward far-field scattering amplitude without considering the important role of the backscattered amplitude. Therefore the scattering amplitude of | R | is significantly lower, and the scattering amplitude of | T | is higher than the other four MST models. We observed that the calculation results of the four frameworks qualitatively correspond with the FEM simulation results, especially at a low frequency (0–1.0 MHz). At a low frequency, we also observed that the FEM simulation value is between the theoretical value of LC and WT frameworks and more close to the LC framework. In the LC framework, the quadratic coefficient of the scatterer number is modified, and the results show the difference when compared with the WT framework, indicating the correction to N 2 is necessary. Depending on the long-wavelength limit, only considering the lowest-order truncation coefficient of T n in the expansion of the scattering fields can be sufficient. Simpler formulas for calculating the effective wavenumber in the LB framework Eq. 18 and LM framework Eq. 19 can be obtained by using the zero and first orders of T n   ( n ≤ 2 ) ( K / K 0 ) 2 = 1 − 4 i N ( T 0 + 2 T 1 + 2 T 2 ) / K 0 2 − 32 N 2 ( 4 T 1 2 / 3 +   304 T 2 2 / 105 + T 0 T 1 + 8 T 0 T 2 / 3 + 78 T 1 T 2 / 15 ) / K 0 4 , (18) ( K / K 0 ) 2 = 1 − 4 i N ( T 0 + 2 T 1 + 2 T 2 ) / K 0 2 − 32 N 2 ( T 1 2 +   2 T 2 2 + T 0 T 1 + 2 T 0 T 2 + 4 T 1 T 2 ) / K 0 4 . (19)

It is worth noting that the approximation formulas for LB and LM frameworks are very similar, which is also evidenced by the fact that the LM framework curves are particularly close to the LB framework curves, as shown in Figure 2A and Figure 2B.

AMMs with negative EMD and/or EBM have been studied extensively, which can cause a lot of exotic acoustic properties. To enhance our analysis, we have investigated the EMD and EBM for the microstructure composites. Figure 2C and Figure 2D show the real parts of the EMD and EBM, according to two different MST frameworks, the Angel and Aristégui framework (Eq. 10 and Eq. 11) and LC framework (Eq. 13 and Eq. 14). As observed in Figure 2A and Figure 2B, the transmission and reflection coefficients of steel scatterers in water display no resonance peaks. Therefore, the real parts of the EMD and EBM are decreasing from the static limit to unit one with increasing frequency. At a high frequency, the real parts of the EMD and EBM approach unity, and the imaginary parts of the EMD and EBM tend to zero. These results demonstrate that the EMD and EBM tend to reach that of water at a high frequency; in other words, there is no contribution in the density and modulus from the steel scatterers.

For the second case, rubber scatterers in water, it is notable to see resonances at specific frequency regions. Figure 2E and Figure 2F show the moduli of transmission and reflection coefficients corresponding to rubber scatterers in water. For transmission coefficients (Figure 2E), the frequency regions of the resonances peaks in the MST models are consistent with those in FEM simulation, except the Foldy framework misses three main resonance peaks in 0.4–3 MHz. Once again, it indicates that Foldy only considered the forward far-field scattering amplitude. However, out of resonant frequency regions (around 1.02 and 1.87 MHz), MST models provide results significantly different from FEM simulations. It can be explained that the scatterers reach a steady state at the specific narrow resonant frequency, and a large amount of energy is stored at this frequency region. As less energy is transmitted in the fluid, this reduces the intensity of interactions between scatterers meanwhile. Once out of the steady state, MST models only calculate the coherent part but the incoherent part superposition decreases the amplitude of the transmitted wave. We also noticed that the oscillation of the WT framework at resonance frequency regions appeared visibly smoother than that of LB, LM, and LC frameworks, which can be regarded as the influence of the microstructure type. LB, LM, and LC frameworks have an additional term of N 2 compared to the WT framework. When using rubber as scatterers, the velocity within the scatterers is very low in comparison with that of the matrix phase, exhibiting strong Mie-type resonance (Leighton, 1997). In other words, the additional term of N 2 causes more oscillations at resonant frequency regions due to the type of microstructure. As for the reflection coefficient (Figure 2F), the calculation results of MST models have the same trend line of change with FEM simulations but provide a lower value. The possible reason is that the incoherent part superposition decreases the amplitude of the transmitted wave and increases the amplitude of the reflected wave, but the multiple scattering theory only calculates the coherent part. For Mie-type scatterers, the results based on MST models and FEM simulations are inconsistent in the non-resonant frequency region, compared with hard-type scatterers. It indicated that the MST models are only suitable for observing the trend change and resonant frequency regions of Mie-type scatterers.

Figure 2G and Figure 2H show the real parts of the EMD and EBM corresponding to rubber scatterers in water. We noted that in six particular frequency regions, the moduli of transmission and reflection coefficients show obvious resonance peaks, which are around 0.57, 1.02, 1.44, 1.87, 2.28, and 2.75 MHz, respectively, within the ultrasound frequency range of 0.4–3 MHz. At six resonance peak frequencies, most of the incident energy is transmitted. As shown in Figure 3, rubber scatterers exhibit a variety of resonance modes at these specific frequencies, whereas at the non-resonant frequencies (Figure 3A), most of the energy is reflected by the composites. In the aforementioned narrow frequency regions, the EMD or EBM also show a negative value. There are negative resonant peaks in the EMD near 0.57 and 1.44 MHz (Figure 2G). For the Angel and Aristégui framework (AA framework), negative resonant peaks are formed in the EBM near 1.02 and 1.87 MHz. For the LC framework, negative resonant peaks in the EBM are observed near 0.57, 1.02, 1.44, and 1.87 MHz (Figure 2H). In the LC framework, the quadratic coefficient of the scatterer number is modified, which results in the difference compared with the AA framework. We can see that EBM results exhibit three more peaks compared to the AA framework (Figure 2H), but these are not shown in Figure 2D. Once again, it can be considered to be influenced by the type of the microstructure. When using ultra-slow Mie-type microstructure scatterers, the correction to N 2 term of the effective wavenumber should be considered. In summary, Figure 2 illustrates that both the EMD and EBM can be made negative near the resonance frequency by choosing appropriate soft resonant scatterers, opening the possibility for searching ultrasonic metamaterials with desired negative properties. These narrow resonance frequency regions depend on the scatterer mechanism, size, and fraction. In this work, we mainly discussed the influence of the microstructure type (scatterer mechanism) and illustrated that soft rubber scatterers can show negative acoustic parameters compared with hard steel scatterers.

The aforementioned results illustrate that the MST theoretical calculation models are capable of analyzing ultrasonic metamaterials with negative effective properties. We chose soft rubber scatterers with low velocity and density suspended in water and a high fraction of scatterers to exhibit strong Mie-type resonances to recognize negative effect properties (Ba et al., 2017). In addition, the application of this ultrasonic metamaterial is further explored, depending on the research on complementary acoustic metamaterials (Shen et al., 2014); the imaging aberration skull layer is used to demonstrate whether the ultrasonic metamaterial in this study can enhance the ultrasound transmission through the skull. A numerical model was further achieved by using FEM simulations by COMSOL Multiphysics software. As shown in Figure 4, after the ultrasonic metamaterial layer is added in front of the skull layer ( ρ = 2000 k g / m 3 ,     c = 2500 m / s ) , the ultrasound energy through the skull was significantly increased at 1.02 MHz. By analyzing the ultrasound intensity magnitude along the axis of ultrasound, it is found that the transmitted ultrasound intensity increases by up to 200% when the ultrasonic metamaterial layer is added. It is worth noting that the ultrasonic metamaterial exhibits negative EBM at 1.02 MHz, indicating that the EBM is the main factor affecting ultrasound wave penetration through the skull layer.