On modifications of Newton's second law and linear continuum elastodynamics

Graeme W Milton, John R Willis


In this paper, we suggest a new perspective, where Newton's second law of motion is replaced by a more general law which is a better approximation for describing the motion of seemingly rigid macroscopic bodies. We confirm a finding of Willis that the density of a body at a given frequency of oscillation can be anisotropic. The relation between the force and the acceleration is non-local (but causal) in time. Conversely, for every response function satisfying these properties, and having the appropriate high-frequency limit, there is a model which realizes that response function. In many circumstances, the differences between Newton's second law and the new law are small, but there are circumstances, such as in specially designed composite materials, where the difference is enormous. For bodies which are not seemingly rigid, the continuum equations of elastodynamics govern behaviour and also need to be modified. The modified versions of these equations presented here are a generalization of the equations proposed by Willis to describe elastodynamics in composite materials. It is argued that these new sets of equations may apply to all physical materials, not just composites. The Willis equations govern the behaviour of the average displacement field whereas one set of new equations governs the behaviour of the average-weighted displacement field, where the weighted displacement field may attach zero weight to ‘hidden’ areas in the material where the displacement may be unobservable or not defined. From knowledge of the average-weighted displacement field, one obtains an approximate formula for the ensemble averaged energy density. Two other sets of new equations govern the behaviour when the microstructure has microinertia, i.e. where there are internal spinning masses below the chosen scale of continuum modelling. In the first set, the average displacement field is assumed to be observable, while in the second set an average-weighted displacement field is assumed to be observable.


1. Introduction

The recent recognition of possibilities for designing materials for ‘cloaking’ objects to electric fields and more generally to electromagnetic radiation (Greenleaf et al. 2003; Alú & Engheta 2005; Leonhardt 2006; Milton & Nicorovici 2006; Pendry et al. 2006) and for realizing phenomena such as ‘superresolution’, where line or point sources were found to have almost singular images both in the quasi-static limit (Nicorovici et al. 1994) and beyond it (Pendry 2000), has stimulated research into similar possibilities for acoustic or elastodynamic waves (Fang et al. 2006; Milton et al. 2006). The present work is concerned with the application of Newtonian mechanics to composite objects and materials containing such objects as material elements. It is argued that, at the level of observation, some generalization of Newton's laws will appear to hold, even when the underlying microstructure conforms strictly to Newtonian mechanics. Similar generalizations are required in the context of a continuous medium, if it is inhomogeneous on a length-scale smaller than the scale of practical observation. It has been understood for many years that the elastostatic response of a composite is expressible in terms of a stress–strain relation which is non-local in space (Beran (1968), Willis (1983), Diener et al. (1984) for random media, Bakhvalov & Panasenko (1984) for media with periodic structure). In the context of inhomogeneous elastodynamics, effective constitutive relations are non-local in space and time. Furthermore, stress is coupled not only to strain but also to material velocity, and momentum density is coupled not only to velocity but also to strain. This follows directly from the formulation of inhomogeneous elastodynamics developed by Willis (1981a,b) and was stated explicitly by Willis (1997). Perhaps most surprising is the fact that the ‘effective density’ operator is a second-order tensor: this was demonstrated rigorously by Willis (1985).

Even earlier, Berryman (1980) using approximation schemes had suggested that for wave propagation in fluid–solid composites, the effective density is not simply the average density. This was confirmed beyond doubt recently by Mei et al. (2006), following the discovery, by Sheng et al. (2003), of the breakdown in the conventional mass law of sound transmission in solid composites near resonance. In addition, Ávila et al. (2005) studied an elastic matrix containing very compliant inclusions and found that the effective density is frequency dependent and anisotropic. Thus, there is already convincing theoretical and experimental evidence that the dynamic density is not the same as the static density.

The present work contributes towards the objective of designing materials with unusual response to elastic waves by first building upon the work of Sheng et al. (2003), Liu et al. (2005) and Ávila et al. (2005), and considering rigid materials containing cavities in which reside masses connected to the cavity wall by springs. (All these models have the characteristic feature that their macroscopic properties depend on the resonant properties of substructures and it is known that such resonant substructures give rise to interesting effects, such as the colours of Venetian glasses associated with suspensions of nearly resonant small metal particles (Maxwell Garnett 1904) or enhanced or negative magnetic permeability associated with split ring resonators (Schelkunoff & Friis 1952; Pendry et al. 1999)). In our models, the motion of the rigid material apparently violates Newton's law owing to the vibration of the internal masses: ‘force equals mass times acceleration’ only applies if mass is replaced by ‘effective mass’, which is a non-local operator in time (equivalently, a function of frequency under any purely harmonic forcing). Furthermore, unless the microstructure is specially constructed so as to be isotropic, the ‘effective mass’ is a second-order tensor, not a scalar. The observation of Willis (1985) should therefore not be a total surprise. However, in the spring models, the macroscopic velocity is taken to be the velocity of the rigid matrix, whereas in Willis' equations the macroscopic velocity is the average velocity which is not well defined in the spring models, or more generally in composites with voids because it is unclear what one should take for the velocity in the void phase. Clearly, there is a need for some generalized continuum elastodynamic equations which govern the effective response of bodies with or without voids. These new, more general, effective relations are developed here; they use some weighted average of the local velocity as the macroscopic velocity, where one is free to take the weighting function to be zero in the void phase or in ‘hidden regions’ which are not accessible to measurement.

In one of our models of seemingly rigid bodies, we allow a spinning top to be in the hidden region. Consequently, it turns out that the body has an effective moment of inertia which is frequency dependent. This type of construction motivates considering a hierarchical composite, in which each ‘material element’ itself has microinertia, of the type just described. This requires some further generalization of the formulation of Willis (1981a,b). We first take the average displacement field as the observable and then we extend the analysis to allow an average-weighted displacement field as the observable.

2. Mass density matrices and modifications to Newton's second law of motion for seemingly rigid bodies

The concept of a mass density matrix may be foreign but was actually introduced by Willis while deriving ‘effective’ constitutive relations for elastic wave propagation in random media, although his definition of the mass density matrix differs from the one that we adopt. In fact, it is quite easy to obtain materials with an anisotropic mass density matrix as we will show in this section. Incidentally, it is well known that electrons and holes in semiconductors can have anisotropic effective masses, due to their interaction with the periodic potential (e.g. eqn (12.29) in Ashcroft & Mermin (1976)). In particular, the effective mass is strongly anisotropic for silicon, germanium and bismuth.

The models we investigate here are based on the models of Sheng et al. (2003) and Liu et al. (2005) which show that the effective mass density at a given frequency is not simply the average mass density and can even be negative. We modify the model so that one gets a complex anisotropic mass density, and for simplicity we take the matrix to be rigid so that one does not have to worry about the effect of the elasticity of the matrix.

First, consider the one-dimensional model of figure 1, where n cylindrical cavities of length d have been carved out from a bar of rigid material. In the centre of each cavity is a lead sphere of mass m and radius r which is attached to the ends of the cavity with two, possibly viscoelastic, springs each having the same complex spring constant K. For simplicity, we neglect gravity and assume the springs are unstretched in their equilibrium configuration, i.e. when they have length d/2−r. We assume that everything varies harmonically with time with frequency ω. The spring constant K may depend on the frequency ω. The external force acting on the whole bar F(t), the force f1(t) the spring on the left of each mass exerts on the bar and the force f2(t) the spring on the right of each mass exerts on the bar are given byEmbedded Image(2.1)where Embedded Image, Embedded Image and Embedded Image are complex. We haveEmbedded Image(2.2)in which P(t) is the momentum of the whole bar.

Figure 1

A one-dimensional material where the mass depends on the frequency ω and can be negative.

In the first cavity, the left wall is at X(t) and the centre of the lead sphere is at x(t), whereEmbedded Image(2.3)where U and u are the complex displacements of the rigid bar and each lead ball, respectively. The velocity of the rigid bar isEmbedded Image(2.4)while the velocity of each lead sphere isEmbedded Image(2.5)Assuming the rigid bar has a mass M0, and applying Newton's law of motion to each lead sphere givesEmbedded Image(2.6)or equivalentlyEmbedded Image(2.7)However, v is unobservable since it is in the hidden part of the bar and we need to relate Embedded Image to the observable velocity Embedded Image. Hooke's law for each spring impliesEmbedded Image(2.8)Substituting back in equation (2.6) givesEmbedded Image(2.9)which has the solutionEmbedded Image(2.10)and this equals U in the limit of infinitely stiff springs (K→∞). ThusEmbedded Image(2.11)where we call M the effective momentum mass or p-mass. Clearly, the mass M can be complex and huge near resonances (when ω2≈2K/m). Moreover, the real part of the p-mass can be negative, as recently discovered by Ávila et al. (2005), Liu et al. (2005) and Mei et al. (2006). A similar effect occurs in electromagnetism where Pendry et al. (1999) claimed that composites with near resonant conducting substructures, such as arrays of split ring resonators (suggested by Schelkunoff & Friis (1952) for enhancing the magnetic permeability due to resonant effects), could have a negative magnetic permeability. This was subsequently justified rigorously in two dimensions by Felbacq & Bouchitté (2005). Asymptotic analysis of the elastic properties of arrays of split ring structures for antiplane elasticity (Movchan & Guenneau 2004) and planar elasticity (Guenneau et al. in press) shows that they exhibit low-frequency bandgaps, perhaps due to a negative effective density.

If each spring behaves according to a Maxwell model of a dashpot in series with a purely elastic spring, thenEmbedded Image(2.12)where k and η are real and the modulus of the purely elastic spring component and the coefficient of viscosity of the dashpot component, respectively. Substituting back in equation (2.11) givesEmbedded Image(2.13)which is exactly of the same form as the single oscillator model for the dielectric constant ϵ(ω) as a function of ω: e.g. eqn (7.129) in Jackson (1975).

The material illustrated in figure 2 is the two-dimensional extension of this model. It has an anisotropic complex p-massEmbedded Image(2.14)where K and L are the complex spring constants in the horizontal and the vertical directions. Newton's law for this material at the fixed frequency ω now takes the formEmbedded Image(2.15)where Embedded Image is the complex velocity. (The physical force F(t), momentum P(t) and velocity V(t) are obtained by multiplying Embedded Image, Embedded Image and Embedded Image by e−iωt and taking the real part.) Instead of using springs, one could of course fill each cavity with an elastically anisotropic (and possibly viscoelastic) material, with the lead sphere being inserted in the centre.

Figure 2

Schematic of a material with an anisotropic mass density ρ at a given frequency ω.

A further generalization is sketched in figure 3. It has an anisotropic complex effective p-massEmbedded Image(2.16)where the rotation matrices Rj vary from inclusion to inclusion according to their orientations. If each of the springs behave according to a Maxwell model of a dashpot in series with a purely elastic spring, thenEmbedded Image(2.17)where kj and j are real and the moduli of the purely elastic spring component while ηj and νj are real and the viscosities of the dashpot component for the respective springs. Substitution back in equation (2.16) givesEmbedded Image(2.18)Now suppose, instead of harmonically varying velocities, we just applied some time varying velocity V(t) and observed what the momentum P(t) was. There would be some linear constitutive relationEmbedded Image(2.19)where causality implies that the possibly singular integral kernel H(s) is zero when s<0, since the inertial force cannot depend on the accelerations at future times. Fourier transformation of (2.19) gives equation (2.15) whereEmbedded Image(2.20)is the Fourier transform of H(s). Causality then implies that M(ω) is an analytic function of ω in the upper half-plane Im(ω)>0, since the above integral converges there. Taking the complex conjugate of equation (2.20) givesEmbedded Image(2.21)where Embedded Image denotes the complex conjugate of z.

Figure 3

Schematic of a material which can approximate any mass density response function ρ(ω) that is consistent with causality.

Let us assume that H(s) is symmetric for all s which implies that M(ω) is symmetric for all ω, and assume thatEmbedded Image(2.22)

Now at a fixed real frequency the force F(t) and velocity v(t) are given byEmbedded Image(2.23)The average work done on the system in a cycle of oscillation will beEmbedded Image(2.24)This will be non-negative for all choices of Embedded Image if and only if ImM(ω) is positive semi-definite for all real ω>0.

Following Milton et al. (1997), it is convenient to introduce the variable z=ω2 and to study the functionEmbedded Image(2.25)The properties of M(ω) imply that (i) Embedded Image, (ii) G(z) is an analytic function of z in the plane except possibly along the positive real axis, (iii) G″(z0+iδ) is positive semidefinite for positive infinitesimal values of δ when z0 is real and positive, and (iv) G(∞)=0. Thus, G(z) is a symmetric matrix valued Stieltjes function of −z, and as such has the integral representation,Embedded Image(2.26)where the distribution B(u) is symmetric positive semi-definite and possibly singular. This function can be approximated to an arbitrarily high degree (for Im(ω)>δ0 and for any fixed δ0>0) by replacing the distribution by a sum of a large number say n of delta functions givingEmbedded Image(2.27)where without loss of generality (by having some of the uj equal to one another) we can assume the positive semidefinite matrices Bj all to have rank one. This impliesEmbedded Image(2.28)Now for comparison in the model (2.18), we can let all j's approach zero and all ηj's approach infinity (so the springs are almost purely elastic). Then equation (2.18) agrees with equation (2.28) withEmbedded Image(2.29)In other words, the models of figure 3 are sufficiently rich to mimic, in the continuum limit, the response kernel H(s) of any material assuming H(s) is symmetric and equation (2.22) holds.

Of course one could argue that Newton's second law does not need to be changed and holds for ordinary bodies, and in fact the analysis here did assume that Newton's law held at the level of the microstructure. But was this correct or just a good approximation? After all the lead spheres in our models have defects and microstructure at one level or another, and so we should have really used a model like equation (2.18) for the response of the lead sphere if we wanted to be more accurate. The analysis of Willis (1985) predicts that any random composite will have a frequency-dependent mass term. It could be argued that the failure of Newton's law in our models is due to our failure to take into account the moving ‘hidden mass’. But there is always some hidden mass, and in fact an enormous amount of it if some modern cosmological theories are to be believed. Rather than look to deeper and deeper levels of microstructure (perhaps far beyond the reach of experiments) to explain everything using Newton's law, it seems much more sensible to allow modifications of Newton's second law to explain the response of bodies.

There does not currently appear to be any need for modifying Newton's second law of motion at the atomic level. In our linear models, the dependence of mass on frequency is due to the fact that there is dissipation of energy of the springs into heat through damping in the dashpots. This causes ImM(ω) to be non-zero and by analogy with the Kramers–Kronig relations, it follows that ReM(ω) must depend on frequency. At an atomic scale, it is hard to see where the analogous damping would come from, if at all.

In our models, we could have incorporated springs with a nonlinear response and accordingly the resulting Newton's law would be nonlinear. Interestingly, there is increasing evidence that a semi-empirical nonlinear modified Newton's law (MOND, or modified Newtonian dynamics) first introduced by Milgrom (1983a,c), is needed to account for the motions of distant stars in galaxies and other astronomical observations. In a nutshell, the gravitational force is F=GMm/r2 where M is the mass at the centre and m is the mass of the star. The centripetal acceleration of the star assuming it moves in a circle of radius r at speed v is a=v2/r. It is observed that v is approximately constant for distant stars in galaxies. In its simplest form, Milgrom's idea to account for this is to propose that F is proportional to a2 for sufficiently small accelerations. We are not suggesting that a simple spring model with nonlinear springs would explain MOND, but we are suggesting that one could allow for a nonlinear Newton's second law of motion; however, only linear laws are considered in the present work. Recently, Bekenstein (2004) has developed a plausible relativistic version of MOND.

So far, we have had just a finite number of inclusions in the body. But by periodically repeating the microstructure, appropriately scaling quantities and taking the homogenization limit as the cell size shrinks to zero, one obtains ‘n-oscillator’ materials with a p-mass density matrix as a function of frequency taking the formEmbedded Image(2.30)which, in the special case when Embedded Image for all j, has the same basic structure as that found by Ávila et al. (2005) for an elastic matrix containing very compliant inclusions. As follows from our analysis of the analytic properties of M(ω), the function ρ(ω) shares the same basic properties as the permittivity tensor Embedded Image(ω) as a function of frequency. By allowing the microstructure to vary smoothly with position, one gets models where all the real positive parameters ρ0, Embedded Image, Embedded Image and Embedded Image, and the real rotation matrices Rj (satisfying Embedded Image) all depend smoothly on x. Although these materials are all rigid, one could insert spherical or other shaped inclusions of them in an elastic matrix and in this way obtain elastic bodies with anisotropic, possibly complex, effective mass densities. Elastic wave propagation in such bodies will be governed by the p-mass of the inclusions not by their gravitational mass.

The extension of these models to three dimensions is obvious and will not be given here.

3. The momentum and angular momentum of a seemingly rigid body containing a spinning top

Here, as illustrated in figure 4, we consider a rigid body with a spherical cavity of radius r1 containing a spinning top which we model as a ring of mass m and radius r0 spinning without friction on its axis at frequency ω0, but subject to small torques of frequency ω, not necessarily equal to ω0. The struts connecting the ring to its supporting cylindrical axis, and the axis itself are all assumed to be of negligible mass. The rigid shell outside the cavity has mass M0 and (tensor valued) moment of inertia K0.

Figure 4

Schematic of a body containing a spinning ring which appears to have a frequency-dependent moment of inertia.

Let us introduce a very small parameter ϵ, which roughly speaking is the amplitude of the oscillations in the system. The centre of the ring (which coincides with the centre of the cavity) has position and velocityEmbedded Image(3.1)where the vector u is complex with its real and imaginary parts not necessarily being collinear. If y denotes the position of a point in the shell when ϵ=0 then its position when ϵ is small isEmbedded Image(3.2)and the associated velocity of this point isEmbedded Image(3.3)in which ϵω(t) is the physical angular velocity derived from the complex valued angular velocity vector Embedded Image (not to be confused with the frequency ω), whose real and imaginary parts are not necessarily collinear. The centre of mass of the shell is assumed to be at ys when ϵ=0 and atEmbedded Image(3.4)when ϵ is small. By making the substitution y=y1≡(0, 0, r1) in equation (3.2), we obtain the formulaEmbedded Image(3.5)for the point where the top of the axis meets the cavity wall.

Points on the spinning ring are located atEmbedded Image(3.6)whereEmbedded Image(3.7)measures the position of the point on the ring relative to the ring centre in the rotating reference frame where the ring is at rest. We assume that, apart from the small oscillations, the ring is spinning around the x3-axis, so that the rotation matrix R(t) takes the formEmbedded Image(3.8)in whichEmbedded ImageandEmbedded Imageis a complex antisymmetric matrix, which reflects the small wobbles in the spin axis and perturbations in the spin speed. (Observe that [I+ϵRe(Be−iωt)] is an infinitesimal rotation). The right-hand side of equation (3.2) also gives the point x1(t) where the top of the axis meets the cavity wall when we make the substitution y=y1≡(0, 0, r1) yieldingEmbedded Image(3.11)Equating this with (3.5) givesEmbedded Image(3.12)and we will see later that conservation of angular momentum implies that b3=0.

The momentum associated with an arc length dθ of the ring isEmbedded Image(3.13)Integrating this over θ we see that, as expected, the total momentum of the ring at time t is just the momentum of the centre of mass:Embedded Image(3.14)Adding this to the momentum of the shellEmbedded Image(3.15)gives the total momentum of the body:Embedded Image(3.16)

The total angular momentum of the ring isEmbedded Image(3.17)where we have used the fact that dQ(ω0t)/dt=ω0Q(ω0t+π/2). For any matrix A and angle α, we have the identityEmbedded Image(3.18)where eijk is the alternating (Levi-Civita) symbol andEmbedded Image(3.19)This identity allows us to simplify the expression for each element of qr(t) and we haveEmbedded Image(3.20)where H is the antisymmetric complex matrixEmbedded Image(3.21)Substituting this back in equation (3.20) givesEmbedded Image(3.22)whereEmbedded Image(3.23)in which we have set b3=0 because the torque dqr/dt has no third component assuming there is no friction at the contacts between the top and the cavity. The angular momentum of the shell isEmbedded Image(3.24)and so therefore the total angular momentum of the body isEmbedded Image(3.25)whereEmbedded Image(3.26)Thus, the effective moment of inertia of this seemingly rigid body depends on the frequency ω. In addition, the constitutive laws in equations (3.16) and (3.25) take the formEmbedded Image(3.27)to order ϵ withEmbedded Image(3.28)

4. Effective relations for a random medium

This section summarizes an approach that was initiated by Willis (1981a,b). The equation of motion of any classical continuum can be expressed in the formEmbedded Image(4.1)where σ denotes stress, p is momentum density and f is body-force density. This has to be supplemented by constitutive equations that relate stress to strain and momentum density to velocity: in the context of small deformations of a (visco-) elastic medium,Embedded Image(4.2)Here, the symbol * denotes convolution with respect to time. The constitutive operator C and the mass density ρ may depend on position. For a random medium, this dependence is obtained from the underlying stochastic process.

The analysis will also extend to ensembles of locally periodic media having a very small spatial period ϵ: for each realization r=(r1, r2, r3), where ri∈(0, 1), we could takeEmbedded Image(4.3)where the convolution operator C(x, y) and the function ρ(x, y) where y=(y1, y2, y3) are periodic in the each variable yi with period 1. Ensemble averaging then corresponds to integrating r over the unit cube.

For such media, ensemble averaging the equation of motion (4.1) gives1Embedded Image(4.4)This equation could be solved directly if it could be supplemented by ‘effective constitutive relations’ relating mean stress and mean momentum density to mean strain and mean velocity. It was shown explicitly by Willis (1997), though it was implicit from earlier work starting with Willis (1981a,b), that they took the formEmbedded Image(4.5)Here, * still denotes convolution with respect to time. The ‘effective’ operators are also non-local (integral) operators in the space variables. Willis (1985) demonstrated rigorously, for a model composite in which only the density varied, that the effective density operator was indeed a second-order tensor. The form of equation (4.5) is explained briefly below.

Introduce a comparison material with constitutive properties C0 and ρ0 in place of the actual C and ρ, and express the constitutive relations (4.2) in the formEmbedded Image(4.6)so thatEmbedded Image(4.7)The equations of motion (4.1) now implyEmbedded Image(4.8)Taken together with boundary and initial conditions, equation (4.8) can be solved in terms of the relevant Green's function for the comparison medium, to give (in the kind of symbolic notation already introduced)Embedded Image(4.9)where u0 is the solution of the problem for the comparison medium. It follows thatEmbedded Image(4.10)where the operators Sx, etc. are obtained from appropriate derivatives of G. It can be shown that St and Mx are formal adjoints (Willis 1981b): Embedded Image in the sense thatEmbedded Image(4.11)Equations (4.7) now require thatEmbedded Image(4.12)In addition, by using the ensemble averages of equations (4.10),Embedded Image(4.13)Formally, equations (4.13) have solution2Embedded Image(4.14)Finally, by ensemble averaging and employing the ensemble average of equations (4.6), the effective constitutive relations (4.5) are obtained, withEmbedded Image(4.15)The operator T cannot be found exactly but approximations can be generated—see the references.

We note that to evaluateEmbedded Imagefor arbitrary-independent fields a and b, we of course just replace 〈e〉 with a and Embedded Image with b in equation (4.13) and solve for τ and π. By making this replacement and ensemble averaging equation (4.13), we see that T satisfiesEmbedded Image(4.16)

5. The new elastodynamic equations in the absence of spinning subelements

To motivate this section, we first give a simple example which shows that the average velocity Embedded Image is not always a suitable macroscopic variable. Imagine the models of figures 1 or 2 with the rigid material replaced by an elastic phase and with the springs absent, so that the cores are disconnected from the matrix. The propagation of elastic waves in such structures will clearly not be influenced by the mass of the core. In addition, it is unclear what value one should assign to the displacement field in the void phase. If one fills the void phase by a highly compressible light isotropic material then the displacement in the void phase, and hence Embedded Image and the effective operators which enter the Willis equations, will still be dependent on the Poisson's ratio of this highly compressible material even in the limit as its moduli vanish. For such materials, it seems more natural that the ‘effective density’ corresponds approximately to the average density excluding the mass of the cores and the macroscopic velocity corresponds approximately to that of the matrix.

To allow for such a description, we look for equations which apply to the ensemble average of a weighted displacement field w(x)u(x, t), where the weight function w(x) may be zero in some hidden regions where there are voids or which is not accessible to measurement: here w(x)=wr(x) is dependent on the particular realization r in the ensemble of materials. To derive our equations, we first assume that there are no voids and w is fairly close to 1. Then we argue that the same effective equations generally apply when voids are present and when w(x) is zero in some hidden regions.

To allow for the case that only a weighted mean 〈uw〉 of 〈u〉 may be observable, where uw(x,t)=w(x)u(x,t), note first thatEmbedded Image(5.1)Hence, employing also (4.14),Embedded Image(5.2)where H denotes the Heaviside step function. (We assume u=0 for t≤0.)

Define ew as the strain associated with uw:Embedded Image(5.3)It follows thatEmbedded Image(5.4)and soEmbedded Image(5.5)whereEmbedded Image(5.6)From equation (4.16), we see that Rw is the identity operator when w(x)=1. Thus, we expect that it should have an inverse at least when w(x) is close to 1. Using this inverse allows us to express the Willis constitutive relations (4.5) in terms of 〈ew〉 and Embedded ImageEmbedded Image(5.7)whereEmbedded Image(5.8)This constitutive relation (5.7) together with ensemble average of equation (5.3),Embedded Image(5.9)and equation (4.4) govern the behaviour of 〈uw〉. The behaviour will not be dependent on the choice of C0 or ρ0.

Although this derivation was based on the Willis equations which assume that there were no voids present in the structure, we may now consider the situation in which voids are present. More specifically, we assume that f=0 in any regions where voids are present for any realization in the ensemble (otherwise (4.1) has no solution) and we assume that we can fill the voids with highly compliant light material and obtain a material satisfying the assumptions of continuum elasticity. (If f is non-zero in a region where the material is highly compliant, then the homogenized equations can take a different form: see eqn (4) in Ávila et al. (2005).) We make this replacement and set w(x) close to zero in the highly compliant phase. Although the derivation of equations (4.5) and (5.7) assumed that the fluctuations in C(x) and ρ(x) were not large, and w(x) was close to 1, we expect no dramatic change in the form of these equations to occur when C(x), ρ(x) and w(x) approach zero in some reasonable way in some reasonably shaped regions. (At some later time this needs to be made mathematically precise.) Hence, equations (4.5), (4.14) and (5.7) should still hold. Even though 〈e〉 and Embedded Image depend on the way that C(x) and ρ(x) approach zero in these regions, we expect that the fields 〈σ〉, 〈p〉, 〈ew〉 and Embedded Image should not. Therefore, it seems reasonable to expect the constitutive relation (5.7) to hold generally even when there are voids (provided w(x) is taken to be zero in the void regions).

By the same argument as made in §2, the relations (5.7) should be the correct ones to apply to any continuum body (without internal masses) where there is microstructure below the limit of observation or where only the ensemble average of a weighted displacement field w(x)u(x, t) is accessible to measurement, where the weight function w(x) may be zero in some hidden regions.

6. The stress energy tensor and obtaining the ensemble averaged energy density

The preceding equations allow one to determine the ensemble-averaged stress tensor and the ensemble-averaged momentum. But from the viewpoint of special relativity, these together with the ensemble-averaged energy density are just components of the ensemble-averaged stress energy tensor, whose components mix under Lorentz transformations. Therefore, for completeness it makes sense to find the equation which determines the ensemble-averaged energy density, especially if one is interested in how the equations transform under space-time coordinate transformations (see section 3 of Milton et al. (2006)). Let us assume the force, momentum, stress, displacement and velocity are all proportional to some small parameter ϵ. Then, to first order in ϵ, the stress energy tensor Embedded Image and the four-force density F as they appear in special relativity with x4=ict take the formEmbedded Image(6.1)where E is the energy density. The equation Embedded Imageij,j=Fi is then equivalent to the equations of conservation of momentum (with i=1, 2, 3) and conservation of energy (with i=4). Conservation of momentum implies the equation of motion (4.1). Conservation of energy impliesEmbedded Image(6.2)where the right-hand side is zero because the work Embedded Image done by the forces is of order ϵ2. At first sight, this equation seems a little strange since E includes not only the rest energy of the body, but also mechanical and thermal energies in part created as the result of viscous damping. To calculate the thermal energy density and hence E, one clearly needs to have an equation governing heat transport which has not been included in our analysis. The answer to this apparent paradox is that phonons (which are responsible for heat transport) carry momentum. To avoid these complexities, we assume that there is no heat transport, i.e. zero thermal conductivities.

If the ‘conventional’ equations hold at the microscopic scale we can write,Embedded Image(6.3)where the latter is the continuum version of E=mc2, and the first equation neglects the momentum associated with the flux of elastic energy, because it is of order ϵ2. (By including this, one would see that in special relativity, the momentum depends on the stress and hence on the strain.) Then substituting equation (6.3) into equation (6.2) givesEmbedded Image(6.4)which is the familiar law of conservation of mass, and mass is equivalent to energy, which is why equation (6.2) is the law of conservation of energy.

Ensemble averaging (6.2) givesEmbedded Image(6.5)which with equation (4.5)2 yields (to first order in ϵ)Embedded Image(6.6)Having solved the new elastodynamic equations for 〈uw〉, one could use (6.6) to solve for the energy density 〈E〉 assuming its value is known at t=0, and neglecting heat transport.

Another point is that the symmetry of the stress energy tensor Embedded Image can be called into question as we will see in §7. The standard argument (Misner et al. 1973) for establishing its symmetry is that if it were not symmetric then the torque created on an arbitrarily small cube would set that cube into arbitrarily great angular acceleration, which rightly seems absurd. However, the loophole in this argument is that the continuum equations break down at the level at which the underlying microstructure becomes apparent, so one cannot consider arbitrarily small cubes. For a discussion of modifications to Einstein's theory of general relativity which allow for non-symmetric stress energy tensors see Hehl et al. (1976). The antisymmetric part of the stress–energy tensor is known as the spin tensor. There is also experimental support for theories where the stress is not symmetric: see http://silver.neep.wisc.edu/∼lakes/Coss.html.

7. The new elastodynamic equations with microinertia

The Willis equations and the extensions (5.7) apply to the response of material with ‘conventional’ tensor of moduli C and mass density ρ. By contrast, §2 developed the proposition that the apparent densities of mass and inertia would both be tensor-valued convolution operators with respect to time; in addition, the macroscopic velocity and spin were taken to be those of the rigid matrix. Again for a composite containing inclusions of this type, the average displacement field 〈u〉 might not be observable, but perhaps the ensemble average of a weighted displacement field w(x)u(x, t), where the weight function w(x) may be zero in some hidden regions, would be accessible to measurement. These two aspects will be incorporated sequentially in the equations now to be developed.

(a) Allowing for inclusions with microstructure

Assume that the composite has conventional elastic (or viscoelastic) moduli, but (at the chosen scale of modelling) it has a density matrix ρ and also displays microinertia so that its total kinetic energy density isEmbedded Imagewhere Embedded Image and ω=(1/2)curl(v). As will be discussed in a subsequent paper, this formula for the total kinetic energy is not applicable to all materials which below the chosen scale of modelling contain spinning masses. For example, the introduction of the local moment of inertia tensor Embedded Image assumes that the spinning mass elements are rigid. Nevertheless, we believe the above formulation should be a good approximation when the medium below the chosen scale of modelling is a composite with seemingly rigid inclusions of the type discussed in §§2 and 3 embedded in an elastic matrix. We also remark that there is a large body of literature regarding continua which below the level of modelling contain rotating elements: see the review article of Lakes (1995) and references therein.

The coupling between v and ω is admitted because, in line with the description of the kinematics of the microstructure in §§2 and 3, v might not be the velocity of the centre of mass of the element. Equivalently, its linear and angular momentum densities are, respectively,Embedded Image(7.1)It seems natural to take the stress in such a medium to be non-symmetric. Therefore, let sji represent the i-component of traction on an element whose normal is in the j-direction. In addition, let s=σ+a, where σ is symmetric and a is antisymmetric. The equations of motion for the medium may then be expressedEmbedded Image(7.2)assuming conventional body-force density f but no body-couples. The symmetric part of the stress is given by σ=C*e; by contrast, the antisymmetric part a is not given constitutively, but has to be such that (7.2)2 is satisfied. The contribution from a to the equation of linear momentum (7.2)1 can be taken into account by making a constitutive relation for the total stress s. The governing equations then comprise the equation of motion (7.2)1, together with the linear momentum relation (7.1)1 and the new ‘constitutive relation’Embedded Image(7.3)where dkl=ul,k (d=u) denotes the displacement gradient andEmbedded Image(7.4)It is convenient also to express the linear momentum relation in the formEmbedded Image(7.5)where Embedded Image.

A composite made of such ‘unconventional’ material can be treated in the same way as in §6. A comparison material with constitutive properties C0 and ρ0 is chosen, relative to which the following equations are obtained:Embedded Image(7.6)(Note that if microinertia is absent and C′ satisfies the usual symmetries of elasticity tensors then to avoid singularities it could be necessary to choose a C0, which does not satisfy these additional symmetries.) Expressing their solution (cf. (4.14)) asEmbedded Image(7.7)the corresponding effective constitutive relations take the formEmbedded Image(7.8)withEmbedded Image(7.9)

(b) Employing a weighted average

Finally, to allow for the case that only a weighted mean 〈uw〉 of 〈u〉 may be observable, where uw(x, t)=w(x)u(x, t), note first thatEmbedded ImageHence, employing also equation (7.7),Embedded Image(7.10)

Define dw as the displacement gradient associated with uw:Embedded Image(7.11)It follows thatEmbedded Image(7.12)whereEmbedded Image(7.13)in which δ denotes the Dirac delta-function so that these convolutions reduce to simple products.

It is now possible to express the effective constitutive relations (7.8) in terms of 〈dw〉 and Embedded Image asEmbedded Image(7.14)from which an equation governing 〈uw〉 follows. We remark that although we took w(x) to be a scalar function, the analysis would still go through if w(x) was replaced by a matrix valued function w(x); this would allow giving more weight to certain displacement components in some areas of the composite.

8. Discussion

For a time-harmonic disturbance at a radian frequency ω≠0, the effective constitutive relations (7.14) have the general formEmbedded Image(8.1)and the associated equation of motion for u isEmbedded Image(8.2)in which C, S, D, ρ are non-local operators in space, depending on ω as a parameter. (Even if they were local and isotropic, S could conceivably be proportional to the alternating tensor.) Their explicit calculation is difficult if not impossible, but they can nevertheless be approximated. Equation (4.12) can be conditionally averaged to generate a statistical hierarchy of equations, and then approximations can be developed by making some closure assumption such as the quasi-crystalline approximation (Lax 1952). An attractive alternative to this is to note that the hierarchy is equivalent to a single ‘stochastic variational principle’ (Willis 1981a,b, 1985) and to substitute suitable classes of trial fields into this. Possibilities of this type exist also for the generalized formulation leading to equations (7.8) and (7.14); these will be developed elsewhere.

The Willis equations and their extension (5.7) apply to the response of material with ‘conventional’ tensor of moduli C and mass density ρ. One can develop the same formalism for treating composites where the equations (5.7) hold at the microlevel. The predicted ‘macroscopic equations’ will again take the same form as equation (5.7). This implies that the general form of the relations (5.7) is closed under the operation of building a ‘hierarchical composite’ in which the constituent materials already conform to relations of this type, i.e. these equations are closed under homogenization (with suitable restrictions on the microstructure). This is expected since one can take an ensemble average by subaveraging within subensembles and then averaging the subaverages, with appropriate weighting. It is possible also to demonstrate that the form (7.14) is closed under coordinate transformations (Milton et al. 2006) whereas the form of the conventional continuum elastodynamic equations is not invariant under these transformations. Such closure properties for conductivity and electromagnetism (Post 1962; Kohn & Vogelius 1984; Ward & Pendry 1996) are a consequence of the geometric structure of electromagnetism (e.g. section 3.2 of Zolla et al. 2005) and are important for the development of electrical and electromagnetic cloaking materials. Corresponding possibilities for the generalized elastodynamics considered here remain to be explored.

Again let us assume the force, momentum, stress, displacement and velocity are all proportional to some small parameter ϵ. Then the introduction of a non-symmetric stress, necessarily entails that the stress energy tensor Embedded Image is non-symmetric and (with x4=ict) takes the formEmbedded Image(8.3)to order ϵ, where E(x, t) is the energy density and c2ϕ(x, t) is the energy flux density in the absence of forces. The equations of conservation of momentum and energy are contained in the equation Tij,j=Fi where F is the four-force density given by equation (6.1)2. When the stress energy tensor was symmetric, we could determine all elements of its ensemble average to first order in ϵ (assuming the force, momentum, stress and displacement were all proportional to ϵ and that there was no thermal conduction) as follows from §§5 and 6. When microinertia is present then §7 shows how to obtain 〈s〉 and 〈p〉, but it remains unclear how to obtain 〈ϕ〉 and 〈E〉 in this case.

It is remarked finally that a formulation for inhomogeneous electrodynamics, exactly parallel to the elastodynamic formulation of Willis (1981a,b), was developed by Willis (1984). It leads to a general set of effective constitutive relations in which the electric displacement Embedded Image and magnetic induction field Embedded Image are coupled to both the electric field Embedded Image and magnetic field intensity Embedded Image through the equationsEmbedded Image(8.4)When Embedded Imageeff, λeff and μeff act locally in space (though not necessarily locally in time) these are the constitutive equations for bi-anisotropic media (see, for instance, Serdikukov et al. 2001 and references therein). The presence of the adjoint operator λeff† in (8.4)2 is a consequence of the formal self-adjointness of the underlying equations of electrodynamics. This symmetry would be broken if the more general procedure of §7, of generating equations for weighted averages, were adopted. We also remark that such couplings appear for electromagnetic fields in the presence of a gravitational field: see eqn (4) in the problem of section 90 of Landau & Lifshitz (1975). In isotropic media (called bi-isotropic media), the coupling λeff is non-zero only in chiral materials (i.e. materials with mirror asymmetric structures) and leads to novel effects such as optical activity. We expect that many novel physical effects will similarly be found to be associated with the equations (4.5), (5.7), (7.8) and (7.14).


Special thanks go to Marc Briane whose joint collaboration with us on the companion paper (Milton et al. 2006) stimulated this research. Also they thank Ping Sheng and Michael Vogelius for their useful conversations, and Andrey Lagarkov, Alexei Efros, Ross McPhedran, Evgenii Narimanov and Nikolay Zheludev for pointing out relevant references. Luc Tartar is thanked for his inspiring remark that homogenization theory is key to understanding physical laws. The referees are thanked for their comments. Graeme Milton is grateful for support from the National Science Foundation through grant DMS-0411035, and from the Institut de Recherche Mathématique de Rennes (IRMAR) and the Université de Rennes and thanks the Institut National des Sciences Appliquées de Rennes for its hospitality.


  • It is assumed here that the body force f is sure. For a discussion of the more general case, see Luciano & Willis (2000).

  • This general structure can be confirmed by obtaining a series solution by iteration, for a weakly-heterogeneous medium. See, for instance, Willis (1997).

    • Received August 7, 2006.
    • Accepted November 16, 2006.


View Abstract