Anisotropic mass density by two-dimensional acoustic metamaterials

We show that specially designed two-dimensional arrangements of full elastic cylinders embedded in a nonviscous fluid or gas define (in the homogenization limit) a new class of acoustic metamaterials characterized by a dynamical effective mass density that is anisotropic. Here, analytic expressions for the dynamical mass density and the effective sound velocity tensors are derived in the long wavelength limit. Both show an explicit dependence on the lattice filling fraction, the elastic properties of cylinders relative to the background, their positions in the unit cell, and their multiple scattering interactions. Several examples of these metamaterials are reported and discussed.


Introduction
Recently, there has been great interest in studying the properties of sonic crystals (SC), a name that specifically defines periodic distribution of sound scatterers embedded in a fluid or a gas [1]. Their properties in the low-frequency limit (homogenization) have been studied by several groups in the last few years [2]- [15] for their potential applications as refractive devices. The refractive properties of these systems are controlled by their effective acoustic parameters; i.e. the speed of sound, c eff , and the dynamical mass density ρ eff . The problem of calculating the effective parameters of a heterogeneous medium have been studied in the past [16,17], but always from a statistical point of view. The underlying periodicity of SC makes it possible to calculate c eff from the acoustic band structure in the low frequency limit by using a plane wave expansion (PWE) [3], but this method does not provide the value of ρ eff . More recently, Mei et al [10] working in the framework of multiple scattering theory (MST) were also able to obtain from the acoustic band structures the expressions for both c eff and ρ eff , and claim their validity for any filling fraction. However, we demonstrate here that their expressions are only valid at low volume fractions and coincide with those already obtained by these authors dealing with finite systems [11,12].
In practical situations, SC are made by clusters with a finite number of cylinders, and therefore their properties in the long wavelength limit have been studied by means of MST without doing any averaging [11,12,14]. In other words, each scatterer was considered individually and their positions and mutual interactions into ideal hexagonal and square lattices (isotropic lattices) were fully taken into account. Analytical expressions were obtained for the isotropic mass density and sound velocity as a function of the SC parameters.
In the present work, we go a step further and develop an MST-based procedure to characterize the homogenization of two-dimensional (2D) SC as a function of the positions and elastic properties of a scatterer in the unit cell. We have derived rigorous expressions for the dynamic ρ eff as well for c eff that have been used to design anisotropic acoustic metamaterials 3 whose properties can be tailored by changing the positions of the scatterers in the unit cell and/or their material constituents.
It is interesting to point out that quite recently Cummer and Schurig [18] have predicted that acoustic cloaking similar to that previously proposed for electromagnetic waves [19,20] should be possible by means of acoustic materials having anisotropic mass density and sound speed. In this regard, this work demonstrates that acoustic materials with anisotropic parameters are physically realizable. So, in a future work we will look for the recipe to build acoustic materials with the properties shown in [18].
The paper is organized as follows. Section 2 gives a brief introduction to the basic equations of acoustic wave propagation in anisotropic media. In section 3, we summarize the main ingredients of MST that will be used in the following sections. Particularly, we deduce the secular equation from which the homogenization parameters are explicitly obtained.
Afterward, in section 4, the effective parameters are rigorously derived by studying the low frequency limit of the acoustic band structure given by the secular equation. Firstly, the sound speed is analytically derived and an expression is obtained that depends on the so called anisotropy factor of the lattice. Secondly, the well-known expression for effective bulk modulus, B eff , is demonstrated to be hold for the systems. Finally, the anisotropic mass density is also derived from the previously known c eff and B eff . Results for parameters c eff and ρ eff are explicitly calculated for several lattices of rigid cylinders in air. Section 5 reports results for the case in which the elastic properties of cylinders are fully taken into account. The work is summarized in section 6.

Anisotropic wave equation
An anisotropic acoustic medium can be characterized by its anisotropic mass density tensor, ρ i j , and the scalar bulk modulus, B [21]. The state equations for the particle velocity vector v and for the acoustic pressure field P are As will be seen later, it is worth working with the reciprocal density tensor, so that equations (1) can be cast as: which has been obtained by using the unitary property i ρ −1 ki ρ i j = δ k j . Now, taking the derivative with respect to x k and adding in the subindex k From equation (1b) the wave equation for the acoustic pressure is finally obtained: 4 If we assume plane wave solutions with angular frequency ω: where k = k cos θx + k sin θŷ is the wavenumber. Then, its insertion in equation (4) gives: And the speed of sound, c = ω/k, is which defines a tensor, c i j , such that this is the well-known Wood's law [22] generalized to the case of an anisotropic medium.

MST
A comprehensive account of MST can be found in textbooks [23] and articles [24,25]. Therefore, this section briefly reports the basic ingredients to understand the rest of the paper. Consider a cluster of N parallel cylinders with arbitrary transversal section located at R α (with α = 1, 2, . . . , N ) and embedded in an acoustic medium characterized by its sound speed c b . Let us also assume that an external field, P ext , with frequency ω and wavenumber k impinges the cluster: where k = ω/c b and (r, ϕ) are the polar coordinates of an arbitrary point in the 2D space. The total pressure field will be given by P = P sc + P ext , where P sc is the total scattered field by all the individual α cylinders: where H q is the qth order Hankel function of first kind, and (r α , ϕ α ) are the coordinates with the origin translated to the center of the α-cylinder, i.e. r α = r − R α , as shown in figure 1. (A α ) q are the coefficients to be determined. The total field impinging on the α-cylinder can be expressed by where the coefficients (B α ) q are related to the (A α ) q through the T matrix [26]: being (T α ) qs the elements of the T matrix associated with the α-cylinder. Expressions for the T matrix for a fluid-like cylinder and for an elastic cylinder are known and can be found elsewhere [25,27]. After a few manipulations, the solution of the problem is: where and

Band structure calculation
If the cylinders are ordered in a 2D lattice, their α-positions are defined by the Bravais lattice R α = n 1 a 1 + n 2 a 2 , where n 1 and n 2 are integers and a 1 and a 2 are the primitive vectors: a 2 = a 2 cos φx + a 2 sin φŷ.
Bloch theorem relates the A-coefficients at any arbitrary α-site with those corresponding to the cylinder at the origin of coordinates: where K is a Bloch wavenumber. This relationship applied to (13) allows to obtain: whereG (19) has been obtained by considering that all cylinders are equal (i.e. T α = T ; ∀α). Also, it can be cast in terms of sums in the reciprocal lattice by using the method developed in [28]: where In these expressions, G h = K + h 1 b 1 + h 2 b 2 is a vector (in the reciprocal space) obtained by the translation of primitive reciprocal lattice vectors b 1 and b 2 and adding the Bloch wavevector K. V d is the area of the 2D unit cell (V d = |a 1 × a 2 | in figure 2) and R min is the smaller distance between corners of the unit cell. With no external field (A ext = 0), the equation (18) becomes where the subindex 0 has been omitted for simplification. This is a set of coupled linear equations that in matrix form is M · A = 0, the elements of the M matrix are: The solution of the secular equation det M = 0 determines the acoustic bands K(ω).

Effective parameters for lattices of rigid cylinders
In what follows, we present a rigorous derivation of the effective parameters for the general case of lattice cylinders embedded in a nonviscous liquid or a gas. The theory developed here is valid for all kind of material cylinders: rigid, fluidlike and full elastic cylinders. Firstly, the effective speed of sound is obtained by studying the low frequency limit of acoustic bands. Then, the effective bulk modulus is obtained from T 00 , the first diagonal element of the T matrix corresponding to a cluster that is anisotropic. Finally, the anisotropic effective mass density is obtained from the relation c 2 eff = B eff /ρ eff . In practical applications the condition of rigid cylinders (i.e. ρ ≈ ∞) is the common situation encountered when working with solid cylinders (made of, for example, metals like lead, iron or aluminium) embedded in air. So, in this section, we present a comprehensive analysis of results for this particular case and results for elastic cylinders will be studied in section 5.

Effective speed of sound
In the long wavelength limit, k → 0 (i.e. ω → 0), the dispersion relation, K(ω), for sound waves propagating in a periodic medium becomes linear and an effective speed of sound can be obtained from c eff = lim ω→0 (ω/|K(ω)|), which can also be cast as: where c eff = c eff /c b is the effective speed of sound relative to that of the background. Hereafter, an overlined variable will be used to denote the corresponding quantity normalized to that of the background. Moreover, k is the wavevector in the embedded background; i.e. k = |k| = ω/c b and K(k) = K cos θx + K sin θŷ (θ being the polar angle defined by K).
To study this limit, it is convenient to define frequency-normalized coefficientsÂ q such that Here, we are dealing with lattices of circular-shaped cylinders, which have diagonal T matrices. Then, the equation for coefficientsÂ q iŝ and we obtain for the M matrix: It was shown in [15] that the T matrix elements for both the elastic and the fluid-like cylinders have the same asymptotic form at large wavelengths. In fact, the elements can be given (in the low frequency limit) as a function of quantitiesT q that are independent of k: In this limit T q S H s−q k |s|+δ s0 −|q|−δ q0 ≈ iT q S Y l k |s|+|q|+δ s0 +δ q0 and, therefore, the dispersion relation is determined from: where the matrix elements ofM arê The secular equation (29) can be solved only numerically. The dimension ofM is determined by the truncation on the angular momenta needed to get convergence. For |q| Q max and |s| Q max , the dimension ofM is (2Q max + 1) × (2Q max + 1). Therefore, the calculation of the dispersion relation is not an easy task and to obtain an analytical expression for c eff (k) is impossible if one follows this solving procedure. However, appendix B shows that there exists a block matrix of dimension 3 × 3 containing only the dependence on c eff . Also, it is demonstrated The secular equation associated with this matrix can be analytically solved for c 2 eff : It is important to remember that and contains lattice sums on the reciprocal space. These quantities as well as the rest of the variables in the expressions above are fully described in appendix B.
Note that the expression for c 2 eff takes the form of an angle-dependent speed of sound in an anisotropic medium, where the components of the velocity tensor are An alternative expression is where  (37), for several values of the ratio, a 2 /a 1 , as a function of the angle, φ, between lattice vectors. Anisotropy disappears for the 2D isotropic lattices (square and hexagonal) corresponding to the parameters described in table 1. Table 1. Effective properties of several 2D topologies studied here. a denotes the length of the smaller primitive translation vectors of the corresponding lattice (see figure 2).
Primitive basis Topology Symmetric We see that anisotropy in (32) and (33) comes from factor . In factor , which is given in (B.39), the main contribution to the anisotropy comes from (0) (see (B.20)). The value of (0) given in (B.20) allows to introduce the so called parameter of anisotropic strength A , which is defined as:  The strength of anisotropic effects predicted by A is analyzed in figure 4, where the diagonal elements of the sound speed tensor, c 2 x x and c 2 yy , are plotted for the case of rigid cylinders in three different anisotropic lattices and compared with the corresponding results for the hexagonal lattice, which is isotropic. Note that the maximum possible value for R a , which is determined by the touching condition between neighboring cylinders, depends on the lattice geometry. It is seen in figure 4 that the more anisotropic behavior corresponds to the case a 2 = 2a 1 and φ = 75 • , which has the larger value of A (see figure 3).

Isotropic lattices.
The case of isotropic lattices deserves special attention because results has been published by two different research teams [3,10]. The solution for this case is easily obtained from matrix (31) by introducing the isotropy condition, that is = 0. After straightforward manipulations, the secular equation can be cast as: The analytical solution of this equation is: This solution contains relevant terms of multiple scattering interaction that has been forgotten in the solution given by Mei et al [10], who were also working in the framework of MST. Particularly, our results reduced to those in [10] when we impose in (39) the condition = 1.
In other words, when it is assumed that multiple scattering interactions are neglected. In fact, it has been shown by us [11,12] that this condition is only valid at low filling fractions. The parameter is responsible for the abrupt decrease of speed of sound when the filling fraction approaches the condition of close-packing as shown in figure 4. This behavior is not shown in figure 1 of [10]. Our results fully agree with those found in Krokhin et al, which used a PWE (see figure 1 in [3]).

Wave propagation.
The sound propagation through anisotropic lattices has its own interest and will be discussed in a separate paper, which will be published elsewhere. However, it is possible to advance the behavior expected by simply looking at the so called refractive index ellipsoid, which we have introduced here in acoustics in a manner similar to that in optics: Two index ellipsoids have been plotted in polar coordinates in figures 5 and 6 for two different anisotropic lattices and for several values of cylinder radius R a . In figure 5 it should be noted that the principal axes are rotated 22.5 • with respect to the x-axis of the lattice. However, it is remarkable how the principal axes are slightly rotated with respect to the x y-axis for the lattice studied in figure 5, and more important is the tilted angle which depends on the filling fraction of the lattice. The wave propagation will respond to the index ellipsoid, in such a way that slow propagation is expected along the direction defined by the longer side of the ellipsoid and faster propagation will take place along the direction defined by the smaller side of the ellipsoid.

Effective bulk modulus
The effective parameters of a cluster of cylinders in which the underlying lattice is isotropic (hexagonal or square) were obtained in [12] from the scattered field by the given cluster. In brief, the method associates to the cluster an effective T matrix, T eff , and assumes that (in the long wavelength limit) this matrix must be equal to that of a homogeneous isotropic cylinder. So, it was demonstrated that for a cluster of N cylinders the coefficient of the lower order term in the series expansion of element (T eff ) 00 is where (T α ) 00 is the corresponding coefficient of the α-cylinder in the cluster. For the case of equal cylinders, the matrix elements (T α ) 00 are all identical and where R a and B a are the radius and bulk modulus of cylinders, respectively. Here, we are dealing with a cluster based on anisotropic lattices and, therefore, it is expected that such cluster behaves (in the regime of large wavelengths) as an effective anisotropic fluid-like cylinder with some effective radius R eff that also has to be determined. Appendix A shows that the coefficientT 00 of an anisotropic fluid-like cylinder is equal to that of an isotropic fluid-like cylinder. Therefore, following the method developed in [12] i π R 2 where the left side in this expression is the coefficient (T eff ) 00 of the homogenized cluster with effective bulk modulus B eff and R eff being its effective radius. R eff can be obtained by a simple approach: it is assumed that the fraction of volume occupied by the cylinders is equal to the filling fraction of the underlying lattice, f , that is Now, with the value of R eff it is possible to obtain B eff from (43): This is the standard averaging of bulk moduli that has been previously shown to be valid for isotropic 2D lattices of composites made of elastic cylinders [14,15]. With B eff and the velocity tensor, we are ready to determine the reciprocal density tensor and, then, to fully characterize the anisotropic medium.

Effective mass density
With the expressions derived above, the reciprocal density tensor are derived from ρ −1 eff (θ ) = c 2 eff (θ)/B eff : It is important to note that the reciprocal density tensor (and the effective density) does not depend on the bulk modulus of background and cylinder. In other words, the effective density only depends on the lattice structure, its filling fraction and the density of cylinders relative to the background. The elastic nature of cylinders will be only present in the effective density for high filling fractions, where the higher orders of the T matrix will be present in both the and factors. Figure 7 plots the behavior of the reciprocal density tensor as a function of the cylinder radius (in units of a) for three different anisotropic lattices. For the sake of comparison, the results for the isotropic lattices, hexagonal and square, are depicted in figure 7(a). Note that the stronger anisotropy is achieved for the case (c), which corresponds to the lattice having the larger value of anisotropy strength A (see figure 3).

Isotropic lattices.
For the case of isotropic lattices ( = 0) the effective velocity (39) can be also cast as: in which: The second factor in the right hand side of (47) is the effective bulk modulus B eff given in (45). It can be demonstrated that for low enough f (i.e. =1) the expression (48) reduces to that obtained by Berryman [17] for the dimensionality parameter d = 2 (see also equation (2) in [10]).

Effective parameters for lattices of elastic cylinders
When the ratio between acoustic impedances of cylinders and background Z a /Z b is not large enough the condition of rigid cylinders (i.e. ρ = ∞) is not valid and the sound propagation inside the cylinders has to be taken into account. This is the usual case when working with solid cylinders embedded in water. Therefore, the full elastic properties of cylinders must be considered in the corresponding T matrix. As a consequence, the effective parameters of metamaterials based on solid cylinders embedded in a fluid, like water, present a rich variety of behavior depending of the ratio Z a /Z b , the lattice topology and the fraction of    lattices (square and hexagonal) and two anisotropic lattices. As anisotropic lattice we have studied one (a 1 = a 2 and φ = 75 • ) characterized by a very small anisotropic strength parameter (A = 0.001) and another (a 1 = 3a 2 and φ = 75 • ) in which this parameter is more than one order of magnitude larger (A = 0.042). Results for the slightly anisotropic lattice (A = 0.001) in the left panels of figures 8-10 show that the values of their effective parameters are in between of those calculated for the hexagonal (φ = 60 • ) and square (φ = 90 • ) lattices and the difference between diagonal elements is very small. Results for the stronger anisotropic lattice are depicted in the right panels of the same figures. They show that diagonal elements show appreciable differences that increase with cylinder radius and should be observable in acoustic experiments. Also note that the difference between diagonal elements decreases with decreasing density. Therefore, we can conclude that in order to observe strong anisotropic effects in lattices of solid cylinders embedded in a fluid, we have to select a lattice with a large value of A made of cylinders with a density as large as possible in comparison with that of the fluid background.

Summary
To summarize, this work has introduced a method to create acoustic metamaterials having anisotropic dynamical mass density and sound speed. The method is based on the properties of SC in the homogenization limit. It was shown that 2D arrays of cylinders ordered in lattices with symmetries other than the square and hexagonal symmetry behave in the range of large wavelengths as effective acoustic metamaterials having acoustic parameters (mass density and velocity) that are anisotropic. Analytical results for both parameters have been rigorously derived and numerical calculations have been presented for relevant examples like the case of rigid cylinders in air and some elastic cylinders embedded in water. This work should be considered as a demonstration that anisotropic fluidlike acoustic materials can be physically realizable. The possible applications of these structures strongly depend on the properties of wave propagation through them and will be the topic of our next work. However, we can foresee that these structures could be the basis of designing the class of anisotropic materials needed to demonstrate the acoustic cloaking recently predicted by Cummer and Schurig [18].
where J q (·) and H q (·) are the derivatives of Bessel and Hankel function of order q with respect to its argument and The expression ρ −1 : ∇ defines the following operator (in polar coordinates): The first diagonal element T 00 can be obtained from (A.1), which in terms of matrix elements is r T sr H rq = −J sq . (A.8) When both q and s are equal to zero, this expression reduces to It can be demonstrated that the second term in the right-hand side is of higher order in k than J 00 . Now, in the limit of k → 0, the k's dependent exponential in (1) sq and (2) sq can be expanded in powers of k as follows: By including this expansion in (1) 00 and (2) 00 and after integrating in θ it is found that H 00 ≈ −1 + O(k ln k), (A.14) which can be used to obtain T 00 as T 00 = −J 00 /H 00 .
Finally, the coefficient of lower order in k is easily determined Note that this expression is the same as that obtained for the case of a homogeneous isotropic cylinder [12].

B.2. Matrix elements with s > q
To calculate these nondiagonal elements we split the lattice sum into two separate terms: For k → 0 we use the asymptotic expressions of Bessel functions to get Only those terms independent of k will survive to the limit. So, the condition to have the factor S v s−q = f (k) is |q| + |s| + δ q0 + δ s0 = 2 → (q, s) = (0, 1), (−1, 0), (−1, 1). (B.12) In the same manner, the factor S G Therefore, the contributions to the matrix are

B.3. Matrix elements with q > s
In this case, we know that S Y s−q = S Y −(q−s) = S Y * (q−s) , and, by remembering that iT q is a real number, we havê

B.4. Final expression ofM
By defining Q ≡ 2(Q max − 1), we can reorder theM matrix aŝ where η = (ρ a − ρ b )/(ρ a + ρ b ) and ζ = B b /B a − 1, B a (ρ a ) and B b (ρ b ) being the elastic moduli (densities) of cylinders and background, respectively. Finally, θ is the angle of Bloch wavevector with x-axis and (0) is the anisotropy factor defined as The matrix D has the form where η e = (ρ a − ρ b )/(ρ a + ρ b ) for fluid like cylinders and η e = 1 for full elastic cylinders. For the B and C matrices we have

B.5. The 3 × 3 secular equation
The secular equation (29) is impossible to solve analytically since it depends on the number of angular momenta employed to ensure convergence. Instead, we should look for a more suitable secular equation so that its analytical solution will not depend on such convergence condition. To achieve this goal let us define the matrix X : The expression for the diagonal terms are: and therefore