Generalized analytic model for rotational and anisotropic metasolids

An analytical approach is presented to model a metasolid accounting for anisotropic effects and rotational mode. The metasolid is made of either cylindrical or spherical hard inclusions embedded in a stiff matrix via soft claddings, and the analytical approach to study the composite material is a generalization of the method introduced by Liu \textit{et al.} [Phys. Rev. B, 71, 014103 (2005)]. It is shown that such a metasolid exhibits negative mass densities near the translational-mode resonances, and negative density of moment of inertia near the rotational resonances. The results obtained by this analytical and continuum approach are compared with those from discrete mass-spring model, and the validity of the later is discussed. Based on derived analytical expressions, we study how different resonance frequencies associated with different modes vary and are placed with respect to each other, in function of the mechanical properties of the coating layer. We demonstrate that the resonances associated with additional modes taken into account, that is, axial translation for cylinders, and rotations for both cylindrical and spherical systems, can occur at lower frequencies compared to the previously studied plane-translational modes.


INTRODUCTION
Elastic or acoustic metamaterials are structures with subwavelength units that exhibit unusual macroscopic parameters within certain frequencies, such as negative mass density [1-3], negative elastic bulk modulus [4,5], simultaneous negative mass density and elastic bulk modulus [6][7][8][9][10][11], or negative density and shear modulus [12,13], and negative index of refraction [14]. Metamaterials with unconventional constitutive parameters have broadly extended the ability of manipulation and control of mechanical wave propagation. Wave control can arise from the formation of band gaps, that are the frequency bands where the propagation becomes forbidden. Differently from the Bragg band-gaps in phononic crystals [15,16] that are produced based on the collective effects of the periodically-arranged scatterers in the medium, the spectral gaps manifested in metamaterials originates in localized resonances of the material building-blocks. Thus, in metamaterilals the production of stop bands depends on the internal structure of the building units. Another mechanism emerged to achieve extreme material-parameters is based on space coiling [17,18].
Over the last decade, research on anisotropic metamaterials has enhanced the ability to control sound and elastic wave propagation, leading to the proposal and fabrication of new devices [19][20][21][22][23][24][25][26]. The first acoustic metamaterial that was realized is a block of stiff material (epoxy) with uniformly (or periodically) dispersed locally resonant structural units that consist of a high-density lead sphere coated with a soft material. This material exhibits resonance-based spectral-gaps in low-frequency regime where the sonic wavelength is of two orders of magnitude larger than the lattice constants. Thus, this designed material was a major progress as sound attenuation is challenging in low frequencies. An analytical model has been proposed to capture the essence of physics dealing with this type of materials [2]. Since in practice the lead sphere and the host are much stiffer compared with the coating material, in this model the host and coated sphere have been approximated to be rigid within the long-wavelength limit. Considering the two-dimensional (2D) system with coated cylindrical inclusions, and three-dimensional (3D) with coated spheres, effective dynamics of the metasolid was derived in resonance-induced band-gap regimes related to the translational motions. Analytical expressions were given for the effective mass density which becomes negative near the resonance.
In this paper, generalizing the simple model proposed in Ref.
[2], we describe the dy-namics of the 3D three-component composite taking into account the rotational motions of the lead cylinder or sphere, and the matrix. It is shown that when the inclusions are cylindrical, because of their geometrical anisotropy, wave propagation in the direction perpendicular to the axis of cylinders gives rise to different effective density in comparison with the effective density describing the propagation perpendicular to that axis. Components of the anisotropic effective density are given through analytical expressions. Furthermore, we introduce and analytically calculate the effective density of the moment of inertia, as a parameter that describes the effective rotational mode of the material. It has been previously demonstrated that, when the cladding is very soft, the stop-band frequency can be very low [1]. This makes such a composite with a practical size interesting for various applications including sound and vibration isolations. Here, allowing for additional degrees of freedom (DOFs) related to rotational motions and translation along cylinder axis, we show that the model predicts spectral gaps characteristic for each mode of propagation. In particular, with the material components chosen in Ref.
[2], we demonstrate that in the case of spherical inclusions, the band gaps related to rotational mode occur in lower frequencies compared with those corresponding to translational modes, where the respective effective parameters become negative. Also, for the case of cylindrical inclusions, comparison between the stop bands for translational modes along and perpendicular to the cylinders shows that those that occur along the cylinder axis are formed at lower frequencies. In general, for materials made of spherical or cylindrical inclusions, we analyze the occurrence of local resonances related to each mode on the frequency axis, as function of coating-material properties. It turns out that accounting for additional DOFs that results in additional propagation modes improve notably the tunability of the material in terms of size and relevant frequencies. This is significant for applications such as sound and vibration proofing, which requires small-size materials enabling the low-frequency wave attenuation.
Rotational modes in periodic solid composites have been previously studied [10,19,[27][28][29][30][31]. Rotary resonances have been estimated by a simple model in 2D square lattice of glass cylinders in epoxy [27]. Numerical analysis was performed to obtain double negative properties, for density and bulk modulus, by utilizing simultaneously the local translational and rotational resonances of a chiral microstructure with a unit cell of three-component solid media composed of a chirally soft-coated heavy cylinder core embedded in a elastic matrix [10]. Periodic system of rotating resonators coated by an anisotropic heterogeneous elastic material embedded in a matrix, which was modeled by an asymptotic approach, has been reported to produce negative refraction [29]. Thus, by modelling the anisotropic material that is placed between the core and matrix as ligaments, and by coupling the translational motions with rotations, it has been numerically demonstrated that this structured medium could be employed for lensing and mode localization. Due to simultaneous translational and rotational local-resonances, a metasolid made of a chiral microstructure with a single phase material served as a system to experimentally achieve negative refraction [30]. Accounting for rotational modes, a mass-spring model was used to reproduce band gaps in a 2D phononic crystal made of rigid cylinders in epoxy [28]. This model was enhanced to describe dispersion relation for rotational modes of a 2D square array of rubber-coated steel cylinders in epoxy [31]. However, no analytical approach for rotational metasolids based on full elastodynamics has been reported so far.
Here, we study rotational modes and analyze their associated homogenized properties by modeling analytically the classical three-component subwavelength structures through taking general assumptions and using direct first-principle elastodynamics without any fitting parameters. This offers us not only essential physical insights for the effective dynamics of the metasolid but also a quick and efficient guide to design materials with local rotationalresonances. Additionally and in parallel, for each mode we have systematically calculated the effective parameters of the metasolid based on the discrete mass-spring modeling, and compared the results with those arising from the continuum elastodynamic model. This comparison clarifies the limits of the mass-spring based model, in particular in terms of describing the effective-material dynamics related to the second local-resonance phenomenon that manifests itself for all translational and rotational modes.
In the following, we introduce the elastodynamics equations at microscale in Sec. 2 for an arbitrary shape of the inclusions. We then study the case of cylindrical inclusions in Sec. 3, followed by the analysis for spherical inclusions in Sec. 4. In Sec. 5, the effective parameters for translational and rotational modes are calculated for the homogenized media including microscopically either cylindrical or spherical resonators. In the concluding Sec. 6, the main results of the paper are briefly summarized.

MICROSCALE EQUATIONS
The unit cell under consideration consists of a hard inclusion Ω a embedded in a hard matrix Ω b via a soft cladding Ω (Fig. 1). The interfaces Γ a and Γ b located between Ω a and Ω, and between Ω b and Ω, respectively, are assumed to be perfect. The materials constituting Ω b and Ω a are very stiff with respect to the material forming Ω. For this reason, and within the long-wavelength limit, we make the assumption that the inclusion Ω a and matrix Ω b are rigid while the cladding Ω is linearly elastic [2,32].
The elastodynamic formulation of the unit-cell system is based on the assumptions that only the cladding is deformable and that the displacements of the inclusion and matrix are small and the strains inside the soft layer are infinitesimal. Consequently, the kinematics of the inclusion and the matrix can be formulated bŷ where α = a, b, r is the position vector relative to a the origin O of a 3D space, andû represents the displacement field generated by a translational displacement vectorV α and an infinitesimal rotation characterized by the 3D vectorΘ α . Considering the time-harmonic motion of the foregoing unit system, we writeû(r, t) = u(r, ω)e iωt ,V α (t) = V α (ω)e iωt , andΘ α (t) = Θ α (ω)e iωt , where i = √ −1 and ω is the angular frequency. The motion of the cladding made of a linearly elastic isotropic homogeneous material is described by the following Navier equation together with the displacement continuity across the perfect interfaces Γ b and Γ a : where λ and µ are the Lamé constants and ρ is the mass density of the coating material.
For later use, it is convenient to introduce the Helmholtz decomposition: where Φ is a scalar field and Ψ a vector field. Using Eq. (3), the equation (2a) can be recast into the following ones where h = ω/c l and κ = ω/c t with c l = (λ + 2µ)/ρ and c t = µ/ρ being the celerities of the longitudinal and transverse elastic waves, respectively, in an infinite linearly elastic isotropic medium. In the following, we shall analytically solve the Helmholtz equations (4) with the boundary conditions (2b) for two special cases of major interest.

MEDIA WITH CYLINDRICAL INCLUSIONS
The first special case concerns a material with structural units consisting of three coaxial cylinders of height L, uniformly or periodically distributed. The inner cylinder Ω a is a long revolution one of radius a ≪ L while the outer cylinder Ω b is bounded by an inner cylinder of circular cross-section of radius b > a and an outer cylindrical surface with square cross-section (Fig. 2a).
Regarding the geometry of the unit cell, it is convenient to use the cylindrical coordinates associated orthonormal vectors e r , e θ , and e z (Fig. 2b). As it was assumed in Sec. 2, the inner and outer cylinders are rigid, whereas the cylindrical cladding is deformable. In addition, since the dimensions of the unit cell are such that L ≫ b > a, the six DOFs of the motion of Ω a relative to Ω b can be reduced to four: two translations in the transverse plane, one translation along the cylinder axis and one rotation about the latter. Thus, with no loss of generality, the boundary conditions (2b) for the cladding Ω (Fig. 3) can be written as  u| r=α = U α {cos θ α e x + sin θ α e y } + αΘ α e z × e r + T α e z , where U α represents the amplitude of the displacement in the plane direction e α resulting from the rotation of e x by the angle θ α around e z . The last term in the above expressions stands for the axial translation, with T α the amplitude of the displacement along e z . The last DOF, i.e. the rotation around the common axis e z , is proportional to Θ α e z × e r , where Θ α is the infinitesimal amplitude of the rotation. To link the general form of the boundary conditions (2b) to Eq. (5) associated with the case of cylindrical inclusions, it is obvious that V α ≡ U α e α + T α e z and Θ α ≡ Θ α e z (Fig. 3).
Considering the cylindrical nature of both the geometry and boundary conditions, the displacement field u is clearly independent of the axial coordinate z, thereby the potentials Φ and Ψ are also independent of z. Moreover, following Love [33], Ψ can be chosen to take the form Ψ = ψ θ (r, θ)e θ + Ψ z (r, θ)e z . The condition ∇.Ψ = 0 in Eq. (3) is now satisfied if and only if Ψ θ is independent of θ. Thus, the potential Ψ can be expressed as The general form of the solutions to the Helmholtz equations (4) can be written as where J n (resp. Y n ) stands for the n th -order Bessel function of the first (resp. second) kind.
In these solutions, A n , B n , C n , D n , E n , F n , G n , H n , A, and B are unknown constants to be determined.

Boundary conditions (5) and Helmholtz decomposition (3) expressed in cylindrical coordinates lead to
Now, we are able to determine all unknown constants involved in Eqs. (6).
Owing to the linearity of the problem, general solution for the displacement field inside the cladding Ω can be written as u = u (p) + u (a) + u (r) , where u (p) refers to the translation in xy plane, u (a) to the axial displacement in e z direction, and u (r) to the rotation around the axis (O, e z ). The first displacement field u (p) can be expressed as where the first four terms, and the last four terms describe the displacements generated by the translations along x-axis and y-axis, respectively. The unknown coefficients in the above expression are the solutions to the linear systems (A1) and (A2) given in Appendix A. As such, we find in different terms, the solution of the plane translation given in Ref. [2] corresponding to the two DOFs related to translational motions.
The displacement due to the translation along the cylinders' axis, i.e. e z , is expressed as Finally the displacements related to the rotational motions are provided by After applying the boundary conditions (7), we obtain explicitly 1 It is important to note that, these axial and rotational contributions of the motion involves only κ (and not h). This is due to the state of shear dictating such kinematics.
This constitutes a key point in order to achieve an anisotropic effective density, and obtain different characteristics for different directions.
We now use the equations of motion and the stresses that are known in the cladding, to deduce the relation between the displacements of the core cylinder and the embedding matrix. The equations of motion are written as where M a = ρ a (πa 2 L) is the mass of the core cylinder with the density ρ a , and I a = M a a 2 /2 is its moment of inertia. The stress tensor σ = µ(∇u + ∇u T ) + λ(∇.u)I, with u T the transpose of u and I the identity matrix, is known since the displacement u in the cladding is fully determined. After some simple calculations, we obtain the relationship between displacements of the core and the matrix, as follows where γ := b/a. The expressions of R (p) (ω) and g (p) (ω), involved in plane translational motion are given in the Appendix A [Eqs. (A4)], which can also be found also in Ref. [2].
The other components related to axial translation and rotation are expressed as   Table I  rigid-body mode insofar as the core vibrates as a whole within the elastic cladding. In this case, the dynamics of the structural unit suggests that the unit can be reasonably assim-ilated to a simple mass-spring system for the translational [35] as well as rotational mode [32], by taking the hard cylinder as the mass, the cladding as the spring, which is attached to the mass on one side and to the rigid matrix on its other side. The mass-spring system is exited by the motion of the matrix.
Although, as seen in Fig. 4, the mass-spring model describes well the dynamics of the material at lower frequencies including the first resonances, it fails to predict the second resonance behavior for all of the three modes. This is due to the fact that the displacement field at the second resonance which is mainly produced by the standing waves inside the coating layer, while the lead cylinder follows the motion of the matrix and remains practically at rest (see Ref. [2] for plane translation mode). We know that the spring, modeling the coating layer, has been assumed to be massless, and consequently this discrete model by its construction cannot describe the second resonance. That being said, if we consider a more complex discrete model by taking into account a mass for the mode-related springs a second local resonance in the elastic medium will appear. However, the outcome of latter model does not correspond to that of continuum model; particularly because in the mass-spring model the springs function in one dimension, while at the second local resonance, the displacement fields in the coating layer can be distributed in three dimensions taking, e.g. the shear effects into account. The distribution of displacement field in the xy plane for each of the three modes is shown in Fig. 5, close to the first and to the second resonance frequencies. In this figure, comparison of the field pattern at the first resonance for either of the modes with those arising from the same mode but at the second resonance, shows that at the second resonance the displacement field is distributed in space to a significantly greater extent than that related to the first resonance. In addition, the evaluation of the spring constants is performed in a quasi-static regime, which is further from the frequency-dependent nature of the second resonance occurring inside the elastic layer at higher frequencies. The resonance frequencies in Figs. 5a to 5h, correspond to the points (a), (b), ..., and (h) in Fig. 4. The second resonance of the structure almost coincides with the first cladding resonance. This first cladding frequency can be estimated by searching for the first root of the determinant, which is a classical technique for modal analysis. In fact, we can proceed to find the solutions, by imposing a zero-displacement for the lead cylinder within the continuum model, assuming that the resonance occur only inside the elastic medium. Thus, we are led to determine the first root of ∆ (p) , ∆ (a) , and ∆ (r) in order to find the second resonance frequency of the plane translation, axial translation, and rotation mode, respectively, through the following We refer to the equation (A1) in Appendix A for the derivation of the first relation and the definition of the associated block matrix D. The two last relations arise from the equations (11) and (12). The first and second local-resonance frequencies related to different modes, and calculated by different methods, are collected in Table II, in order to easily compare these values associated with their respective method of calculation. In particular, it can be noticed that the relations (17) provide good estimations of the second resonance frequencies.

MEDIA WITH SPHERICAL INCLUSIONS
Here, we study the same kind of problem as in the previous section with the only difference that the geometry of the structural unit is changed: the cylinder of radius a and height L ≫ a is replaced by a hard ball of radius a. An embedding matrix has now spherical cavities of radius b, that are filled with the hard spheres and a concentric elastic cladding such that a < R < b (Fig. 2a). All these structural components are made of the same materials as in the case of cylindrical system (see Table I), so that the same simplifications can be applied to the corresponding kinematics. By virtue of the analysis in the previous section with coaxial cylinders, it is evident to claim, without any calculations, that for an arbitrary motion of the embedding matrix as a rigid body, the same type of motion is induced for the embedded hard ball. The superposition principle, verified for the previous case, suggests us that translations and rotations are totally uncoupled. A translation in a fixed direction or a rotation around a fixed direction, will lead, respectively, to the same translation or rotation for the hard ball, but with possibly different amplitude and/or phase.
Moreover, assuming a total isotropy of the problem, it is now convenient to consider, Without loss of generality, we can choose Θ α = Θ α e z (see Fig. 2b for the coordinate system), so that the solution is of the form u = u θ (R, ϕ)e θ ≡ f (R) sin ϕe θ . The equation

with the boundary conditions
f (R = α) = αΘ α , following (18b). We can find the solution explicitly in terms of the first order spherical Bessel functions of the first -j 1 -and second -y 1 -kind : As in Sec. 3, we derive the relation between the displacement of the lead sphere and the embedding matrix by applying the equation of motion for the coated sphere, and using the fact that the stresses in the cladding can be known through the displacements [Eq. (19)]: where the moment of inertia of the hard sphere I a = 2M a a 2 /5 around e z . We finally obtain provided that Evolution of rotational displacement H (r) * , and translational displacements of the lead sphere with frequency are shown in Fig. 7, for a medium with the filling fraction 47%, the lead sphere radius a = 0.5 cm, and the thickness of the coating layer is 0.25 cm. We note that, here also, with our material and geometrical parameters, the first two local resonances of the rotational mode occur in lower frequencies compared to those of translational mode, as in the case of cylindrical inclusions. The first resonance frequencies for the translation f (p) 0 ≈ 380 Hz and for rotational f (r) 0 ≈ 261 Hz. Similarly to the case of cylindrical system, the frequency of these resonances that correspond mainly to the vibration of the core sphere can be estimated through mass-spring models, leading to the following expressions With our specific parameters, the above formulas give f  0 (γ), which is an increasing function between 1 for ν = 1/4, and 2.71 when ν = 1/2. Fig. 8 shows the evolution of f (p) 0 with γ for three values of ν: 1/8 (0 ≤ ν ≤ 1/4), 3/8 (1/4 < ν ≤ 1/2), and 0.47 that is the value corresponding to displacement behavior in Fig. 7.

EFFECTIVE-MEDIUM PROPERTIES: METASOLIDS
The knowledge of the kinematics of the cladding is of great interest to go forward in the dynamical homogenization of the whole structure, with the hard coated cylinder or sphere.
Indeed, our purpose is now to describe the homogenized behavior of the hard inclusion, the cladding and the matrix; and to obtain the key effective-parameters of the medium.  Table I). The rotational mode are independent of ν.

Homogenization for translational modes
For a given translation of the embedding matrix V b , we aim at finding the equivalent density for the embedded solid (sphere or cylinder) and the cladding, subject to the induced translation V a . The knowledge of the stress distribution in the cladding permits us to evaluate the total force applied on the cladding by the matrix, and subsequently to derive the homogenized density. Indeed we can write where V e is the volume of the region occupied by the hard inclusion Ω a and the cladding Ω, i.e. the region Ω a ∪ Ω, that is bounded by the surface Γ b . The frequency-dependent effective density ρ e is defined for each point inside Ω a ∪ Ω and must coincide with the static effective value given byφ a ρ a +φρ, withφ a andφ being the filling fraction of the coated solid and the cladding, respectively, such thatφ a +φ = 1.
In order to obtain the effective density for the whole unit cell, it remains to include the density of the matrix ρ b , so that we can finally set ρ eff = (φ + φ a )ρ e + φ b ρ b , where φ a , φ and φ b stand for the filling fraction of the coated solid, the cladding and the embedded matrix, respectively, such that φ a + φ + φ b = 1.
The case of spherical inclusions has been studied in Ref. [2], through a result with an analytical expression 2 . Therefore, here we only deal with the cylindrical inclusions, where in particular an anisotropic aspect appears due to the fact that we take into account not only the plane-translation mode, but also the axial mode. (24), the decoupling between the plane and axial translations leads to two different expressions for the effective density ρ e , depending on whether we are dealing with the axial or plane translation. The anisotropy due to the geometrical configuration of the unit cell is obviously the cause for the anisotropy of the effective mass density. We must consider the effective density ρ eff as a second order tensor, and rewrite Eq.
Thus, the final expression will be ρ eff (ω) = ρ where the elements of the block matrices are defined in Appendix A. Fig. 9 shows the evolution of ρ  eff in accordance with the fact that the first resonance occurs for the axial mode (Fig. 4). We explained in Sec. 3 that the (first) resonance of the axial should emerge at lower frequency, independently of the micro-structural parameters. This can be of particular interest, given that it is more challenging to make attenuate low-frequency waves.
We note, however, that the resonance band where the mass density becomes negative, may also be less wide for the axial mode. With the present material and geometrical parameters, the band width for the axial mode ∆f (a) ≈ 118 Hz, while that of plane-translational mode ∆f (p) ≈ 280 Hz.
The effective masses related to the plane translations for media with cylindrical or spherical inclusions, and axial translation modes, can also be calculated through the discrete mass-spring model. Once the effective spring constants k  (Figs. 4 and 7), are ignored by the discrete model for both modes and medium micro-geometries (Fig. 9).

Homogenization for rotational modes
In order to obtain the effective parameter related to the rotations, we apply first the angular equation of motion to the sub-structure Γ b∪Γ formed by the hard inclusion and the cladding. It reads where I e denotes the resulting moment of inertia around e z for the embedded rigid body and its surrounded cladding that move in phase. Following a similar procedure regarding translation modes, we can establish the expression of the density of the moment of inertia i e (ω) such that dI e = i e dV e . Since this unusual quantity is obviously intensive in thermodynamic sense, its homogenization in the medium should be possible. Tensorial behavior is involved to ensure the anisotropic aspects of the problem, hence we write : Similar to the case of translations, the formal expression of the effective density of the moment of inertia can be written as i eff (ω) = (φ a + φ)i e (ω) + φ b i b I.
The spherical system gives rise to an effective density of moment of inertia i eff ≡ i eff I, provided that where, h * 2 (ω) = 5 γκb j 2 (κb)y 1 (κb) − j 1 (κb)y 2 (κb) j 1 (κa)y 1 (κb) − j 1 (κb)y 1 (κa) The isotropy of i eff is clear and it is evident that i eff (ω) coincides with the static value i 0 := φ a i a + φi + φ b i b in the quasistatic limit ω → 0; with i a = 2 5 ρ a a 2 , i = 2 5 ρ b 5 − a 5 b 3 − a 3 , and Employing the same method that was used in the case of spherical system, for media with cylindrical inclusions, we obtain with where i a = 1 2 ρ a a 2 , i = 1 2 ρ(b 2 + a 2 ), and  frequencies is similar to that we observed earlier for the effective mass densities. Here also, after calculating the spring constants for analogous mass-spring system (Appendix B), and then obtaining the displacement Θ a /Θ b (Figs. 4 and 7), the effective mass densities can be expressed according to the discrete model. Again, from Fig. 10 we observe that negative mass densities related to the lower resonance frequency band are appropriately described by the discrete model, while those associated with the higher frequency band, that is largely due to the rotational resonance inside the coating, are disregarded.

SUMMARY AND CONCLUSIONS
A metasolid composed of a distribution of either cylindrical or spherical hard inclusions coated by a soft material has been studied. Generalizing an existing continuum model [2] that describes resonance-induced band gaps related to plain-translational modes, we provide analytical analysis to investigate anisotropic properties in the case of cylindrical inclusions, and also rotational modes for both cylindrical and spherical systems.
We showed that the translational mode along the cylinders' axis and the rotational mode in both cylindrical and spherical systems could generate band gaps at low frequencies. The band gaps that originate in local resonance phenomena at long-wavelength regime, are tunable in terms of frequency by changing the coating's elastic properties and geometrical parameters of the structural unit. Furthermore, for all translational and rotational modes we found that the discrete mass-spring model could provide appropriate description of the material's kinematics including first resonances, but neglects the second resonances that occur at higher frequencies and mainly arise from resonances inside the coating. The mass-spring model also provides simple analytical expressions allowing us to study easily the kinematics of first resonances in terms of the material parameters. We demonstrated that, while the macroscopic material dynamics is dictated by the effective mass densities for translational modes, it is characterized by its effective moment of inertia regarding the rotational modes. Through analytical expressions, it was shown that components of frequency-dependent anisotropic effective densities for translational modes, and effective density of moment of inertia for rotational mode, become negative near corresponding resonance frequencies. These effective parameters have been also calculated based on the mass-spring model, and their limits have been clarified by comparison to the continuum model.
These results establish a step forward for metasolid characterization via introducing the density of moment of inertia as an effective-medium parameter. They give insight to improving the tunability of the metasolid for various applications, such as sound and vibration isolation, and can serve as a physical-based quick guide for optimal design of the material to exhibit desired properties.
Here, we will give the explicit form of the solutions to the Helmholtz equations (4), combined with the boundary conditions for the case of cylindrical inclusions in Section 3.
The boundary conditions (7) and the solutions (6) lead to the following first conclusions                          C n = 0 = D n = G n = H n , ∀n ≥ 2; A 0 C 0 = 0 = B 0 C 0 E 0 G 0 and F 0 G 0 to be determined; A 1 C 1 , B 1 C 1 , E 1 H 1 , and F 1 H 1 to be determined; A 1 D 1 , B 1 D 1 , E 1 G 1 , and F 1 G 1 to be determined.
The expressions of the unknowns E 0 G 0 and F 0 G 0 are established in the article. The coefficients A 1 C 1 , B 1 C 1 , E 1 H 1 and F 1 H 1 involved in (8) where the matrices used in the bloc-matrix formalism are Similarly, the coefficients A 1 D 1 , B 1 D 1 , E 1 G 1 and F 1 G 1 in (8)  The solutions of the above linear systems can be obtained with bloc-matrix formalism, which leads to the final expression for the displacement field involved in (8).
It is easy now to compute the ratio between the relative plane translations U a and U b .
Using (13), after some formal calculations, we obtain We note that the relation (A3) is not only obtained for the plane translation in e x direction (given after solving (A1)), that is U a cos θ a = H (p) (ω)U b cos θ b , but also for the plane translation in e y direction (given after solving (A2)), i.e., U a sin θ a = H (p) (ω)U b sin θ b .
Although this seems clear in the physical sense, the mathematical treatment is relatively tedious.
for translational and rotational motions, respectively. In the above equations k eff is the spring constant, equivalent of the coating-layer that is represented by a spring in 1D tension/compression, and related to the (plane and axial) translational motion of the layer.
Likewise, C eff is the spring constant describing the torsional stiffness of the coating layer that is taken to be analogous to a spring in 1D tension/compression, related to the rotational motion of the layer. For cylindrical inclusions, the expressions of these constants are obtained as where k The expressions for spring constants k eff for media with cylindrical and spherical inclusions can also be found in Ref. [32].