In the tracks of photonic crystals, phononic crystals (Dowling, 2008) have provided a fillip for research in acoustic stop band structures (Kushwaha et al., 1993; Kafesaki and Economou, 1999) within which light or sound is prohibited to propagate due to multiple scattering between periodically spaced inclusions. In 2000, Liu et al. provided the first numerical and experimental evidence of frequency dispersive elastic parameters of locally resonant structures for elastic waves in three-dimensional arrays of thin coated spheres (Liu et al., 2000): The effective parameters were shown to turn negative where low frequency stop bands occur. This important work paved the way toward acoustic analogs of electromagnetic meta-materials, such as fluid-solid composites for the control of pressure waves (Auriault and Bonnet, 1985; Auriault, 1994; Liu et al., 2000; Goffaux et al., 2002; Ho et al., 2003; Elford et al., 2004; Hirsekorn et al., 2004; Li and Chan, 2004; Movchan and Guenneau, 2004; Fang et al., 2006; Mei et al., 2006; Guenneau et al., 2007b; Hou et al., 2007; Wang et al., 2008; Chalmers et al., 2009; Norris, 2009; Lemoult et al., 2011; Auriault and Boutin, 2012; Boutin, 2013; Boutin and Becot, 2015). Using asymptotic methods for fields in multi-structures (Kozlov et al., 1999; Movchan et al., 2002), it has been proposed to use arrays of cylinders with a split ring cross section as building blocks for two-dimensional localized resonant acoustic structures displaying negative refraction (Movchan and Guenneau, 2004; Guenneau et al., 2007a), based on analogies with split ring resonators (SRRs), introduced by Pendry in the context of electromagnetic waves almost 20 years ago (Pendry et al., 1999). Magnetic activity of metamaterials occurs near resonances of SRRs, which was essential for instance in the design of the first electromagnetic cloak (Schurig et al., 2006), and similar designs have been proposed for acoustic cloaks (Craster et al., 2013).

In a series of articles, the research group of Auriault developed asymptotic models of locally resonant structures (Auriault and Bonnet, 1985; Auriault, 1994; Auriault and Boutin, 2012), some of which predate the birth of acoustic metamaterials with the seminal work of the research group of Ping Shen (Liu et al., 2000). Low frequency stop bands in arrays of 2D Helmholtz resonators with elongated necks have been studied in (Movchan and Guenneau, 2004; Guenneau et al., 2007a) as mentioned above, with further predictions and measurements of sound transmission in elastic shells in air in Krynkin et al. (2010). However, 3D Helmholtz resonators with elongated neck position inside the cavity and periodically distributed in a fluid have been introduced in (Boutin, 2013; Boutin and Becot, 2015). Therein, Boutin and coauthors have shown that the macroscopic fields (pressure P and velocity V) satisfy the following equation of mass conservation and generalized Darcy Law:

where K is the effective bulk modulus (including thermal effects) of the matrix, ω the resonant frequency, i² = −1, Φ and T are the porosity and dynamic permeability tensor (including viscosity effects) when considering the resonators as perfectly rigid spheres, η is the viscosity of the fluid-matrix, ∣Y∣ is the volume of the 3D unit cell Y and Q is the flux emitted by the resonator into the fluid matrix in response to the pressure acting on it. The case of the Helmholtz resonator with a single elongated neck has been not only derived, but also experimentally tested, in Boutin (2013) and Boutin and Becot (2015). And the fluid need not be inviscid.

There is a vast amount of literature on the role of the neck in such Helmholtz resonators (Groeneweg, 1969; Alster, 1972; Orris et al., 1974; Gaunaurd and Uberall, 1982; Lim et al., 1990; Hinders et al., 1995; Selamet and Dickey, 1995; Mead, 1996; Aberg et al., 1997; Selamet et al., 1997; Chen et al., 1998; Baird et al., 1999; Seo et al., 2005; Ivansson, 2006; Duan et al., 2007; Zhao et al., 2007; Zhou et al., 2010; Wang and Mak, 2012; Li et al., 2013). Recently, periodic structures consisting of SRRs have been revisited thanks to refined homogenization techniques, which allow to rigorously link the geometrical parameters to frequency dependent effective parameters (Mercier et al., 2017). Analogs of SRRs have been proposed for in-plane elastic waves (Guenneau et al., 2007b), that allows for bending and rotational vibrations having lower resonant frequencies than compressional vibrations, and thus make good candidates for building blocks of mechanical metamaterials, which might make possible some mechanical cloak, see Kadic et al. (2013) for a review of metamaterials beyond electromagnetics. Milton, Briane and Willis provided a thorough mathematical frame for cloaking for certain types of elastodynamic waves in structural mechanics, in the framework of modified Willis equations (Milton et al., 2006). On the other hand, coupled in-plane pressure and shear waves were numerically shown to be detoured around a finite size obstacle by a specially designed cloak with an anisotropic heterogeneous elasticity tensor without the minor symmetries (Brun et al., 2009) without resorting to Willis media. Nonetheless, SRRs like in Guenneau et al. (2007b) make an interesting candidate to achieve such Cosserat-type metamaterials. Actually, a similar type of resonant elastic structure already led to a negatively refracting medium (Zhu et al., 2014), where the elastic chirality (which is a hallmark of a Cosserat medium) was put forward in conclusive experiments. The need for mechanical metamaterials with elasticity tensor without the minor symmetry motivates the extension of homgenization results of Mercier et al. (2017) to the Navier system. Importantly, the acoustic wave equation is invariant under coordinate changes (Norris, 2008, 2009), so acoustic metamaterials via geometric transform can thus achieve enhanced control of pressure waves via a simpler route. A simple example of a periodic system is a periodic lattice, whose dynamic response is well defined by Green's functions. Dynamic Green's functions in periodic lattices and their asymptotics were analyzed in Movchan and Slepyan (2014) and Vanel et al. (2016). Localization of electromagnetic waves within a two-dimensional grating of spheres embedded in the three-dimensional space and the high-frequency homogenisation approximations were studied in Maling et al. (2017).

A downfall of SRRs is that they do not allow for doubly negative acoustic parameters, so one would need a second type of inclusion within the periodic cell to achieve a goal of a negative refractive index in order to design a super lens like Pendry proposed in his seminal paper (Pendry, 2000). In 2004, Li and Chan proposed a design of negative acoustic metamaterial based on a multiple scattering theory approach (Li and Chan, 2004), which has been revisited recently by applied mathematicians using homogenization theory for bubbly fluids (Ammari et al., 2017a,b) with hybridization of Minnaert resonances (Minnaert, 1933). This is somewhat related to homogenization of high-contrast periodic structures (Figotin and Kuchment, 1998). Potential applications of bubbly media in super lensing and total absorption have been experimentally shown (Lanoy et al., 2015; Leroy et al., 2015). Such soft metamaterials (Brunet et al., 2015) are currently investigated by many groupings in the world, but they might not always be compatible with industrial processes.

Although SRRs might not achieve double negative acoustic parameters, they seem to be quite straightforward to engineer for ultrasonic waves. Fang et al. experimentally demonstrated a dynamic effective negative stiffness in a chain of air filled Helmholtz's resonators for ultrasonic waves (Fang et al., 2006). Moreover, it has been also shown using asymptotic techniques that surface water waves propagating within an array of fluid filled SRRs display a negative effective density (Farhat et al., 2009). A focussing effect through a finite array of such resonators was numerically achieved, with a resolution of a third of the wavelength. The negative effective density, reinterpreted as a negative effective gravity, has been further confirmed theoretically and experimentally (Hu et al., 2011, 2017). Fluid filled SRRs have been experimentally tested in a 17 meter long water channel for filtering effects with a dike using grooved vertical cylinders (Dupont et al., 2017). Similar filtering effects have been experimentally demonstrated for sound waves interacting with a doubly periodic array of Helmholtz's resonators shaped as soda cans (Lemoult et al., 2011). However, we note that in that case, there is no neck, just a hole, as soda cans are thin shells. In the present paper, we would like to revisit the concepts of SRRs in the case of pressure waves propagating in a three-dimensional array of 3D Helmholtz resonators with elongated necks. In order to simplify the mathematical setup, we shall consider airborne acoustic waves, in which case rigid (Neumann type) boundary conditions can be considered on resonant elements. Compared to problems of linear elasticity, the present study does not deal with dynamic degeneracies at low frequencies, which may occur for certain type of geometries of elastic systems, resulting in a group of very small eigenvalues being separated from the remaining spectrum.

Let us start with an illustrative numerical result for a spectral problem for the Helmholtz operator within a periodic cubic array of rigid spheres: the unknown is a pressure wave field, here sound in air (wave speed c = 340m.s⁻¹). Neumann boundary conditions are prescribed on the contour of each defect and standard Floquet-Bloch conditions are set on an elementary cell of the periodic structure. The finite element formulation was implemented in the COMSOL Multiphysics Package to compute the eigenvalues and to generate the corresponding eigenfields. We present in Figure 1 the periodic structure we want to study and in Figure 2 the corresponding dispersion diagram for eigenfrequencies ω (in unit of rad.s⁻¹) as a function of the Floquet-Bloch parameter k (in unit of m⁻¹): along the horizontal axis we have the values of modulus of k, where k stands for the position vector of a point on the contour ΓXMU within the irreducible Brillouin zone. We note the absence of bandgaps with the presence of rigid spheres. This lack of intervals of forbidden frequencies motivates the present study: how can one create a stop band without further increasing the size of the rigid spheres?

Let us first recall the finite element set-up. We consider the Helmholtz equation:

where ρ (kg.m⁻³), K (Pa) are the density and bulk modulus of the medium and ω (rad.s⁻¹) is the angular wave frequency of the pressure field p.

Due to the periodicity of the lattice, we look for solutions of (2) in terms of Floquet-Bloch waves (Floquet, 1883; Bloch, 1928; Bensoussan et al., 1978; Brillouin, 1978; Wilcox, 1978; Kittel, 1986). So, for a cubic array of unit cells Y,

where k_x, k_y and k_z are components of the Bloch vector k within the Brillouin zone Y^* = [0, π]³ (Joannopoulos et al., 1995; Gazalet et al., 2013).

The implementation in the finite element package is fairly straightforward (Hladky-Hennion et al., 1991; Nicolet et al., 2004). We first multiply equation (2) by a smooth test function v and using integration by parts, we obtain the so-called weak form of the time-harmonic acoustic equation

where ∂f/∂n = ∇f·n with n the unit outward normal to the boundary ∂Y of Y, and ds the infinitesimal surface element on ∂Y.

We note that the weak formulation holds for heterogeneous fluids as ρ and K can be spatially varying. In particular, this model works for domains such as a homogeneous fluid filled with a periodic array of rigid obstacles. For the finite element implementation, (4) is discretised using test functions taking values on nodes of a tetrahedral mesh of the basic cell (first order tetrahedral elements), see e.g., (Nicolet et al., 2004) for further details. From (4), we note that setting rigid boundary conditions on an inclusion amounts to assuming Neumann (natural) homogeneous data, whereas transmission conditions at the interface between various fluid phases mean that the quantity ρ⁻¹∂p/∂n is preserved across the interface. We note that in the case of airborne pressure waves, the contrast in density between air and inclusions made of metal or even polymer is sufficiently large to assume Neumann data, but this simplification does not hold if we replace air by water.

Let us now consider a periodic array of defects Ω₁, …, Ω_N embedded in an elementary cell Y = ]0;1[³. Let p(x, y, z) satisfy the Helmholtz equation in Y\⋃j=1NΩj¯. We also assume that p satisfies Neumann boundary conditions on the contours of defects, where n denotes the unit outward normal to the boundary ∂Ω_j of a defect Ω_j:

We would like to consider a particular case when the defects Ω₁, …, Ω_N are spherical shells with thin air channels connecting an air-filled interior cavity to the exterior surrounding air. These defects can be modeled as multistructures (Kozlov et al., 1999) in the following way,

where a_(N) and b_(N) are given constants and Πε(N)(j) is the thin channel.

In this section, we derive an asymptotic approximation of the field in thin channels Πε(j), see Figure 3,

where l_j is the length of the j^th bridge, εh_j(t) the radius of its varying cross-section D_ε (parametrized by a positive real t). Here, ε is a small positive non-dimensionnal parameter. To derive the asymptotic expansions, we introduce the scaled variables ξ = (y/ε, z/ε).

Without loss of generality, and for the sake of simplicity, we drop the superscript j. In Π_ε, the time-harmonic wave equation takes the rescaled form

with the Neumann boundary conditions

The field p is approximated in the form

To leading order, we obtain

Hence, p⁽⁰⁾ = p⁽⁰⁾(x) (it is ξ−independent). Assuming that p⁽⁰⁾ is given, we derive that the function p⁽¹⁾ satisfies the following model problem on the scaled cross-section of Π_ε

The condition of solvability for the problem has the form:

Hence, we have shown that to the leading order we can approximate the field p in the thin channel Π_ε by the function p⁽⁰⁾ which satisfies the Helmholtz's equation in one-space dimension. We now assume that the field is periodic over the cell since it is localized. This shows that the average of the eigenfield over the macro-cell vanishes. Indeed, let χ₁ denote the value of the field in the large body Σ of the multi-structure Ω and let χ₂ (which we normalize to 1) denote the value of the field within the complementary area of the macro-cell Y \ Ω excluding the thin channels. Taking v = 1 in (4), we deduce that

This shows that the average of the field p over Y vanishes, hence by neglecting the small volume of the thin channels, we obtain

where meas_Σ and meas_{Y \ Ω} denote, respectively the volumes of Σ and Y \ Ω.

We now consider two cases. The first one is the study of an array of simple spherical shells with either one or six thin channels, and the other one is the study of an array of double spherical shells with one thin channel in each shell. Since we have q thin channels, we have q separate eigensolutions V_j, (j = 1, …, q), corresponding to the vibrations of thin domains Πε(j)

We note that V_j(0) is equal to a non-zero constant and also that V_j(l_j) = χ₁. Next, we need to take into account that some boundary layers occur at the end regions of thin ligaments. These boundary layers are characterized by exponential decay when the boundary conditions for the functions V₀ (the leading term in the asymptotic expansion of the thin bridge solutions) are chosen in a specific way. In our case, integrating (2) over Σ∪Πε(1)∪⋯∪Πε(q) and applying the divergence (or Gauss) theorem, we obtain to order O(ε)

where

All the channels are connected to Ξ, hence, V₁(l₁) = … = V_q(l_q) = V. We note that the boundary layer condition (18) can be interpreted as Newton's second law.

The solution of the problem (16)−(18) (that one can interpret in terms of a spring-mass model like in Figure 3) has the form

where c=K/ρ and the frequency ω is given as the solution of the following equation:

where we invoked (18). Looking at a first low frequency, we deduce an explicit asymptotic approximation

This estimate actually holds for the frequency ω₂ of the upper edge of the phononic band gap. We note that if we take V(0) = 0 instead of V(0) = meas(Ξ)/meas(Y \ Ω), we estimate the frequency of the lower-edge of the phononic band gap. We also notice that the boundary layer condition (18) is only valid at order 0(ε), so it needs to be refined if one wants to improve the frequency estimate in (22).

Thus far, we have only analyzed dispersion curves, but it has been known for over 80 years that one needs to be extra careful regarding analysis of say electronic band structures (Bouckaert et al., 1936), as there are examples showing that stop band edges might be reached strictly within the Brillouin zone, and so it is not enough to describe its edges to characterize the band spectrum (Harrison et al., 2007). We would like to investigate dynamic effective anisotropy effects near the edge of the low frequency stop band, and localization effects induced by local resonators such as in Movchan et al. (2006),Bigoni et al. (2013),Craster et al. (2013), and Llewellyn Smith and Davis (2010) for a frequency on the acoustic band.

Here, we provide representation for isofrequencies near the lower edge of the first stop band for the case of a Helmholtz resonator with one thin channel as in Figure 4. However, we observed very similar features in all other cases. We first note in Figure 13 that the distortion of the isofrequency surfaces increases when we move toward the lower edge of the stop band. A small change in the frequency leads to a dramatic change in the isofrequency surface. Besides, the observation of isofrequency (or slowness) contours in Figure 13 confirms that anisotropy increases near the edge of the stop band as contours become more and more elongated. Researchers in photonics used vanishing group velocity near stop band edges to achieve self-guiding of light in 2D photonic crystals in the past 20 years (Kosaka et al., 1999; Witzens et al., 2002; Chigrin et al., 2003). Such an effect is known in the applied mathematics and waves literature as a dynamic effective anisotropy (Slepyan et al., 1987; Ayzenberg-Stepanenko and Slepyan, 2008; Craster et al., 2009, 2010; Colquitt et al., 2012; Schnitzer, 2017) : for an isotropic medium, the isofrequency surface is spherical, and the corresponding isofrequency contours are circular. On the edge of the stop band, the extreme elongation of isofrequency contours reveals an extreme dynamic anisotropy. We note that some dynamic effective anisotropy of 3D dynamic lattices displays similar isofrquency surfaces to ours near stop band edges (Vanel et al., 2016). This could be used to achieve as aforementioned highly directive phenomena, in a way similar to what was demonstrated experimentally for transverse electromagnetic waves (Ceresoli et al., 2015). Another striking effect is that of wave localization, when one creates one or more defects in a periodic structure. Usually, this is achieved in high frequency stop bands, when the wave wavelength is on the same order as the periodicity. Here, the wave localization is achieved either at resonances of resonators, see upper panel in Figure 14, or as shown in the lower panel of the same figure, at a frequency ω1*=559.2042 rad.s⁻¹ located in a part of the acoustic band sandwiched between the low frequency stop band and another region of vanishing group velocity near the X symmetry point in Figure 4, so this has the additional feature of being below the Bragg frequency regime. We note that from Figure 14, this localized mode can be approximated as a monopole given by the first derivative of the spherical Bessel function of the first kind j_n(ωr/c), where r is the radial position. The lowest root corresponding to the lowest frequency of the cavity is the frequency given by the root of the transcendental equation (simply derived by separation of variables in Helmholtz's equation (2) written in spherical coordinates and assuming that ∂p/∂n = 0 on the boundary of the defect) :

Bearing in mind that the first root of j1′ is 2.0816 and taking a = 1.25m as the radius of a dashed circle within which the defect mode seems to be mostly confined in Figure 14, (38) leads to the frequency estimate ω₁ = 566.19 rad.s⁻¹, which is in good agreement with the finite element computation ω1*=559.2042 rad.s⁻¹. This method of the effective sphere is the analogous concept to effective disk in doubly periodic structures with defects (Movchan et al., 2007).

In this paper, we have seen that it is possible to sculpt the Bloch spectrum of three-dimensional phononic crystals almost ad libitum simply by digging some holes and adding cavities in rigid spheres periodically arranged along a cubic lattice. One of the main achievements of our numerical study is the appearance of ultra-low frequency stop bands at frequencies predicted quantitalively by an asymptotic model, that can be viewed as a 3D counterpart of (Movchan and Guenneau, 2004). We also conducted some elementary shape optimization (by varying the size, diameter and number of channels in a rigid sphere of constant radius) in order to enhance the control of the location and the number of low frequency stop bands, thanks to our asymptotic estimates. Our asymptotic results could be further used for a homogenization study in the spirit of the 2D analysis conducted in Mercier et al. (2017) to address the effective parameters of our acoustic system. However, the acoustic metamaterial we have analyzed would allow only frequency dependent effective bulk modulus, and in order to achieve a frequency dependent effective density and thus double negative acoustic parameters near resonances, one needs to introduce a second type of resonator in the periodic cell. Therefore, achieving a design of an acoustic superlens through negative refraction requires a more complex design. However, a cubic array of rigid spheres does not support any complete stop band, even in the densely packed configuration, so the acoustic metamaterial which we studied has markedly different filtering properties compared with phononic crystals with rigid spheres: multistructures make possible tunable ultra-low frequency stop bands (associated with very flat dispersion curves i.e., localized eigenmodes). Finally, we illustrate in Figure 14 a possible application of the ultra-low frequency stop band in order to localize a mode of a wavelength larger than the pitch of the array of resonators. One can also envisage sculpting a line defect thanks to the removal of a few resonators in a larger macrocell, that would make a low frequency waveguide. We further note that considering shells of same diameter with holes of varying diameters, or spheres with varying diameters and identical or varying holes would make possible graded phononic crystals with a larger low frequency stop band, thanks to the overlap of numerous low frequency stop bands. This could lead to a three dimensional counterpart to the acoustic rainbow put forward in Zhu et al. (2013). Indeed, many embedded spheres would make this possible, since the band diagram 11 suggests N embedded spheres would lead to N low frequency stop bands that might hybridize. Similarly, one could envisage to reflect, detour, or focus, pressure waves using the low frequency stop band within which effective parameters are expected to take negative values as it is now well-established for the two-dimensional counterpart of such kind of acoustic metamaterials. We hope our work will foster experimental efforts toward 3D acoustic metamaterials for airborne and underwater sound filtering effects. At first sight, the implementation of the former seems more straightforward, since the latter would require further theoretical and numerical investigation of conversion of pressure waves propagating in the fluid into pressure and shear waves propagating in the shells. For instance, in the derivation of the asymptotic estimates, the boundary layer analysis will be more involved. Nevertheless, it might be easier to perform underwater experimental characterization of the sonic metamaterial.

GD performed all numerics, AM derived the asymptotic formulae, SG wrote the manuscript and SE participated in the physical interpretation of results. All authors reviewed the manuscript and contributed to the study.