New metamaterials with macroscopic behavior outside that of continuum elastodynamics

Metamaterials are constructed such that, for a narrow range of frequencies, the momentum density depends on the local displacement gradient, and the stress depends on the local velocity. In these models the momentum density generally depends not only on the strain, but also on the local rotation, and the stress is generally not symmetric. A variant is constructed for which, at a fixed frequency, the momentum density is independent of the local rotation (but still depends on the strain) and the stress is symmetric (but still depends on the velocity). Generalizations of these metamaterials may be useful in the design of elastic cloaking devices.


Introduction
The possibility of cloaking objects to make them invisible has attracted substantial recent attention. Alú and Engheta ( [1], [2]) following earlier work of Kerker [3], found that the scattering from an object could be substantially reduced by surrounding it by a plasmonic or metamaterial coating. Milton and Nicorovici [4], expanding on the work of [5] and [6], proved that the dipole moments of a polarizable point dipole or clusters of polarizable line dipoles would vanish when placed within a cloaking region surrounding a superlens, and this has been confirmed numerically [7]. Bruno And Lintner [8] show that only partial cloaking occurs when the polarizable object is not discrete. Ramm [9] and Miller [10] found that cloaking could be achieved by sensing and manipulating the fields near the boundary of the object.
Perhaps the greatest interest has been generated by transformation based cloaking. This type of cloaking was first discovered by Greenleaf Lassas and Uhlmann ( [11], [12]) in the context of electrical conductivity. Their idea was to apply the well known fact (see, for example, [13]) that the equations of electrical conductivity retain their form under coordinate transformations, using a singular transformation which mapped a point to a sphere (the cloaking region), and which was the identity outside a larger sphere. Perturbing the conductivity in the cloaking region corresponds to changing the conductivity at a point in the equivalent problem, and this has no effect on the fields outside the larger sphere. This type of cloaking was extended to the full time harmonic Maxwell equations by Pendry, Schurig and Smith [14], and at the same time Leonhardt [15] discovered a different transformation based, two-dimensional, geometric optics or high frequency acoustic cloaking scheme which only utilized isotropic materials. The former cloaking has been verified numerically ( [16]; [17]; [18]), placed on a firm theoretical foundation ( [19]), and experimentally substantiated in the microwave regime in an approximate way in two-dimensions by Schurig, Mock, Justice, Cummer, Pendry, Starr, and Smith [20] using Schelkunoff and Friis's idea ( [21]) of utilizing metamaterials with split ring resonators to achieve artificial magnetism. Cai, Chettiar, Kildishev, and Shalaev [22] proposed another design which may experimentally achieve approximate two-dimensional cloaking in the visible regime.
Milton, Briane and Willis [23] found that applying transformation based cloaking to continuum elastodynamics would require new materials with very unusual properties. This is because the continuum elastodynamic equations do not retain their form under coordinate transformations, but rather take the form of the equations that Willis ([24], [25], [26]) found described the ensemble averaged behavior of composite materials. One interesting feature is that the density at a given frequency is required to be matrix valued. Willis [27] had proved that his density operator had this property. It also follows from the work of Movchan and Guenneau [28] andÁvila, Griso, and Miara [29], that there exist high contrast materials with a local anisotropic effective density. Simple models with anisotropic and possibly complex density were obtained by Milton and Willis [30]. These microstructures generalized a construction of Sheng, Zhang, Liu, and Chan [31] and Liu, Chan, and Sheng [32], who established that the density can be negative over a range of frequencies. The same conclusion was reached and rigorously proved byÁvila, Griso, and Miara [29]. A comparison of Maxwell's equations, which can be written in the form [in which E is the electric field, ε the electric permittivity tensor, µ the magnetic permeability tensor, and e ijm = 1 (-1) if ijm is an even (odd) permutation of 123 and is zero otherwise] with the equations of continuum elastodynamics [in which u is the displacement field, ρ is the density, and C is now the elasticity tensor] suggests that ρ as a function of the frequency ω could share many of the same properties as ε(ω). This has been established, and in fact for every function ρ(ω) satisfying these properties one can construct a model which has approximately that function as its effective density as a function of frequency ( [30]). Cummer and Schurig [33] investigated transformation based acoustic cloaking in two-dimensions and, due to a mathematical equivalence with the electromagnetic problem, found that it could be achieved provided one could construct the necessary materials with anisotropic density. Anisotropic density is not the only new property needed to achieve transformation based elasticity cloaking. One needs materials where the constitutive law, like that in the Willis equations, couples the stress with not only with the strain but also with the velocity, and couples the momentum density not only with the velocity but also with the displacement gradient through the strain. Also, unlike in the Willis equations, one wants this dependence to be local in space and to apply to a single microgeometry rather than to an ensemble average of microgeometries. Here we show that such unusual behavior can be realized in a model such that, for a narrow range of frequencies, the associated waves have wavelength much larger than the microstructure, and the momentum density depends on the displacement gradient, and the stress depends on the velocity. Curiously the momentum density depends not only on the strain, but also on the local rotation, and the stress is not symmetric. We also construct a variant of the model for which (at a fixed frequency) the momentum density is independent of the local rotation (but still depends on the strain) and the stress is symmetric (but still depends on the velocity). Both models rely heavily on the discovery Sheng, Zhang, Liu, and Chan [31] that one can design structures which, at a fixed frequency, respond as if they have negative mass. Following their ideas, and the simplified constructions of Milton and Willis [30] which incorporate these ideas, figure 1 illustrates a structure with negative mass. In the models considered in this paper both positive and negative masses are idealized as point masses.
This work on metamaterials with macroscopic behavior outside that of continuum elastodynamics, even though they are governed by continuum elastodynamics at the microscale, is preceeded by work on metamaterials with macroscopic non-Ohmic, possibly non-local, conducting behavior, even though they conform to Ohm's law at the microscale, ( [34]; [35]; [36]; [37]; [38]; [39]) by work on metamaterials with non-Maxwellian macroscopic electromagnetic behavior, even though they conform to Maxwell's equations at the microscale ( [40]), and by work on metamaterials with a macroscopic higher order gradient or non-local elastic response even though they are governed by usual linear elasticity equations at the microscale ( [41]; [42]; [43]).
For simplicity we consider a two dimensional model, although it can be generalized to three dimensions. The model is illustrated in figure 2.
2 The momentum density in the model Figure 3 shows the unit cell in the model. The masses are approximated as point masses, and the springs may be taken to have all the same spring constant hK, which scales in proportion to the size of the unit cell. This ensures that the spring network, by itself, responds as a elastic material in the limit h → 0. Figure 1: This structure responds as if it had negative mass to oscillations at a fixed frequency above resonance, such that the spherical core oscillates out of phase with the rigid, but light, surrounding shell. The springs have the same spring constant to ensure that the effective mass is isotropic.

Void Rigid
where c is a parameter between 0 and 1 which controls the inclination of the rods joining E and F to B and D: these rods have length h √ 1 + c 2 . In the limit h → 0 we assume the (infinitesimal) physical displacements at the points x A , x B , x C and x D derive from some smooth complex valued displacement field u(x), and in particular where and x 0 = 0 because the origin is at the center of the unit cell.
At the macroscopic level we choose not keep track of the displacements u E and u F . These are internal hidden variables which however can be recovered from u B and u D . They have the form where the complex vector s remains to be determined. Since there is a rigid rod connecting the points D and E we have, in the limit h → 0, the constraint that w − s must be perpendicular to the line joining D and E, i.e.
Similarly since there is a rigid rod connecting the points B and E we have, in the limit h → 0, the constraint These equations have the solution One can easily see that if c is replaced by −c then s gets replaced by −s, which justifies the formula (2.4) for u F . Now suppose the masses at the points E and F , which we call a mass pair, are respectively chosen to have the form where m is a positive real constant, and δ is a possibly complex constant with a nonnegative imaginary part. Then at any time t the physical momentum in the unit cell is, to second order in h, Since the area of the unit cell is 2h 2 we see that the complex momentum density is in which −iωu 0 is the complex velocity of the unit cell, and s depends on the deformation gradient through (2.7).

The stress in the model
The acceleration of the material will also generate stress. To see this, note that to leading order in h, the masses at E and F must be accelerated by complex forces F = −ω 2 mhu and −F acting on the respective masses (the physical forces are obtained by multiplying these by e −iωt and taking the real part). Let F ED and F F D be the forces which the rods ED and F D exert on the vertex D and let F EB and F F B be the forces which the rods EB and F B exert on the vertex B. Balance of forces at E and F requires that Since these forces are aligned with their respective rods we have where the constants α and β are such that (3.1) is satisfied: Forgetting for a moment the springs and surrounding material, to maintain the motions a force needs to act on the vertex D to balance the forces which the rods ED and F D exert on this vertex. Similarly a force −ht needs to act on the vertex B to balance the forces which the rods EB and F B exert on this vertex. On the other hand, no additional forces need to be exerted on the vertices A and C. Now consider a rectangular sample (with dimensions, say of order √ h, small compared to the wavelength) as in figure 2 with (complex) forces acting on the boundary vertices to maintain the motions. These forces will have two components: an elastic (or possibly viscoelastic) component to compensate for the (complex) stress σ E caused by the (complex) strain in the spring network, and an inertial component to compensate for the inertial stress σ I caused by the acceleration −ω 2 u of the material. The compensating inertial component will be zero on the left and right sides of the network. At the top 5 vertices in figure 2 the compensating inertial component will be approximately 2ht at each vertex. The extra factor of 2 arises because of the extra load that the springs at the top boundary carry to support the inertial component of the forces at the 4 vertices immediately below the top. (At the other interior nodes these forces are balanced). Thus the compensating inertial component per unit length, is t on the top, −t on the bottom, and zero on the sides. By the definition of stress this should be equated with σ I n, where n is the unit normal to the boundary. Thus we deduce that Of course, in the limit h → 0, the stress in the spring network (excluding the masses and connecting rods) will be governed by a normal relation where C is the elasticity tensor of the network. Our choice to define the macroscopic stress through the forces at the vertices of the unit cell (rather than as a volume average of the microscopic stress field, which would result in a symmetric stress) makes good physical sense and is consistent with our choice to define the macroscopic displacement field through the values of the displacement at the vertices of the unit cell. These choices bear some resemblance to the way Pendry, Holden, Robbins and Willis [44] define macroscopic electromagnetic fields through field values at the boundary of the unit cell. In defining the stress in this way it is important that the unit cell be chosen so that each mass pair (of positive and negative masses) is not split by the boundary of the unit cell: each mass pair and their four connecting rods is regarded as an indivisible unit in the same way that in defining the electric polarization field one would not choose to place the boundary in a dielectric material so as to split the positive and negative charges of the constituent atoms.

Summary
In the limit h → 0 the constitutive law takes the form where σ = σ E + σ I is the total stress, v = −iωu is the velocity, and ρ = (δ/2)I is the density. In the basis where σ, p, ∇u and v are represented by the vectors the third-order tensors S and D, from (2.7), (2.10) and (3.5), are represented by the matrices with their solutions will have wavelength much larger that the microstructure. In this circumstance the usual continuum elastodynamic equations normally apply. However this model shows that this is not always the case. The model has some unusual features. Although the net amount of mass in a region of unit area is independent of h, the total amount of positive (or negative) mass in a region of unit area scales like 1/h. However in any physical model, as usual, one never in fact takes the limit h → 0, but instead one sets the microstructure as small as required to get a reasonable approximation to the limiting behavior for a desired set of applied fields. Also the amount of positive (or negative) mass in the unit cell only represents the amount of apparent mass, not gravitational mass, which could be quite different if the apparent mass arises from substructures which are close to resonance. Of course the analysis given here needs some rigorous justification to check, for example, that the displacements can be approximated by a smooth field u(x) in the limit h → 0. Also from a practical point of view it needs to be checked that the behavior is stable to sufficiently small variations in the masses or spring constants from cell to cell.
The model can easily be generalized. One simple generalization has additional masses and connecting rods as illustrated in figure 4, with masses at the points G and H which, when the material is at rest, are located at where c ′ and m ′ are positive constants. By superposition the constitutive law will take the form (4.1) with a density ρ = δI and with the third-order tensors S and D represented by the matrices In particular if, at a given frequency, m ′ c ′ = mc we obtain a model in which the momentum does not depend on the local rotation and the stress is symmetric. With this condition satisfied one has the relations where  The general form of these tensors S ′ and D ′ corresponds with that required for elasticity cloaking: they are real and satisfy the identities S ′ ijk = S ′ jik = D ′ kij . They may serve as a good starting point in the quest to find materials with a suitable combination of properties C, S ′ , D ′ and ρ needed for elasticity cloaking. However, quite apart from this, the novel properties of this new class of metamaterials may have other important applications.
In the design of these models perhaps the key feature is that in the narrow frequency range under consideration the amount of negative mass in the unit cell (due to hidden internal masses moving out of phase with the motion) almost balances the amount of positive mass in the unit cell. This ensures that the momentum density approaches a constant as the cell size approaches zero, while keeping non-trivial the stresses caused by the moving masses in the unit cell. It will be interesting to see if other models with this key feature also have a local constitutive law of the form (4.1), with non-zero coupling terms S and D.