One emerging research area within the fields of acoustic and elastic metamaterials involves designing subwavelength structures that display elastic instabilities in order to generate an effective medium response that is strongly nonlinear. To capture the overall frequency-dependent and dispersive macroscopic response of such heterogeneous media with subwavelength heterogeneities, a theoretical framework is developed that accounts for higher-order stiffnesses of a resonant, nonlinear inclusion that varies with a macroscopic pre-strain, and the inherent inertia associated with an inclusion embedded in a nearly incompressible elastic matrix material. Such a model can be used to study varying macroscopic material properties as a function of both frequency and pre-strain and the activation of such microscale instabilities due to an external, macroscopic loading, as demonstrated with a buckling metamaterial inclusion that is of interest due to its tunable and tailorable nature. The dynamic results obtained are consistent with similar static behavior reported in the literature for structures with elastic instabilities.

## I. INTRODUCTION

Acoustic and elastic metamaterials have attracted significant attention in recent years due to their ability to overcome limitations of conventional materials by exploiting engineered sub-wavelength structures.^{1} Traditionally, the fields of acoustic and elastic metamaterials focus on linear, periodic media with negative effective mass density or stiffness through design of subwavelength resonant structures, and the homogenization schemes developed to describe the effective properties. However, there is now an increasing interest in an effective material response that considers nonlinearity. Most studies of nonlinear acoustic metamaterials rely on spatial periodicity. Examples of periodic, one-dimensional nonlinear metamaterials include mass-spring chains where springs possess cubic nonlinearity,^{2} chains of bistable springs^{3} or beams,^{4} structured waveguides with periodic plates and side holes,^{5} and thin rods with finite internal mass-spring oscillators.^{6} Example nonlinear acoustic and elastic metamaterials in two dimensions include arrays of granular chains^{7} and deformable lattices with voided cavities,^{8} active components,^{9} or resonators containing coated masses, thin beams, and voids.^{10} Proposed applications include non-reciprocal devices,^{9,11,12} phononic switches and rectifiers,^{8,13–15} and tunable nonlinear lenses for filtering, focusing, or improved resolution of acoustic imaging devices.^{7,16,17}

The acoustic metamaterials of interest generate the strong nonlinearity as a result of an engineered elastic compliance rather than by exploiting resonance phenomena. Some initial research on the nonlinear behavior of this type of metamaterial includes buckling beams in the form of negative stiffness honeycombs,^{18,19} woodpile structures,^{20} and compliant, three-dimensional (3D) elements known as buckliballs that are fabricated from soft polymers.^{21} Previous nonlinear acoustic metamaterial research has focused on large elastic deformation and dispersion analysis using Bloch–Floquet methods, which are limited to linear wave phenomena propagating through a structure subjected to large static deformation. However, heterogeneous media with randomly dispersed hyperelastic inclusions are potentially useful for exploiting traditional nonlinear wave phenomena such as phase conjugation, parametric amplification, Raman scattering, and enhanced energy dissipation.^{22} A model to study the homogenized frequency-dependent response of a metamaterial consisting of a dilute concentration of highly nonlinear subwavelength inclusions embedded in a viscoelastic material is developed here.

Example small-scale inclusions that gives rise to large microstructure-induced macroscopic nonlinearity are hollow encapsulated shells that exhibit buckling and are embedded in a host medium. These types of heterogeneous media have been studied for cases of fluid^{23,24} and soft elastic solid^{25–27} host media. Experiments have verified that large values of quadratic-order nonlinearity can be achieved as the microspheres buckle due to an external loading pressure in castor oil^{23} and nearly incompressible soft plastic.^{26} These results for buckling microspheres inspire work in buckling metamaterials. The advantage of buckling metamaterials is their potential to open up a larger design space for effective material properties, one that is more controllable and tailorable than the microspheres, and which may in turn push the limits of achievable nonlinearity.

The present work is influenced by the study of bubbly liquids and porous elastic media. Quasi-static homogenization methods for fluids and fluid-like media can be described using iso-pressure assumptions for both the effective medium and its constituents, such as the immiscible mixture law.^{28,29} Limits of the mixture law have been derived where the constituents possess small and approximately identical nonlinearity.^{30,31} For the case of a soft inclusion within a stiffer matrix, such as porous elastomers^{30} or bubbly liquids,^{28,32} it has been proven that only a small volume fraction is necessary to induce large nonlinearity of the effective medium. Therefore, it is anticipated that the strong nonlinearity of the metamaterial inclusions will exhibit strong macroscopic nonlinearity for dilute concentrations even when frequency dependence is considered.

Quasi-static assumptions are insufficient for describing an effective medium that contains viscoelastic constituents. Early work that took into account relaxation and nonlinearity associated with the viscoelastic elements led to wave equations that predict attenuation, dispersion, and harmonic generation in the resulting effective medium.^{33,34} In the present work, it is of interest to account for the effects of microscale inertia due to the dynamic oscillations of subwavelength inclusions, which can also influence the macroscopic linear and nonlinear properties. Similar dynamic homogenization methods have been applied to cases of air bubbles in water,^{35,36} as well as a soft elastomer matrix containing gas^{37} or voids.^{25} For air bubbles in water, viscous losses have also been taken into account.^{38}

These previously studied models, however, do not capture a more general nonlinear inclusion undergoing purely volumetric expansion, which is of interest to study small-scale metamaterial elements, nor do they account for the influence of microscale inertia of the heterogeneity, or a varying microscopic response due to a macroscopic loading condition. The present work develops a theoretical framework for a dynamic homogenization model that quantifies a dispersive effective medium comprised of compliant, nonlinear elements embedded in a nearly incompressible viscoelastic matrix material. The model is influenced by an effective medium model developed by Zabolotskaya and Soluyan for bubbly liquids.^{35,38} The present work not only extends the previously developed theory to account for an arbitrary spherical inclusion with strong material nonlinearity and non-negligible mass, but is also applied to an incremental response to account for finite deformations.

Expressions for incremental deformation theory are applied to the constitutive relationships of the inclusion and viscoelastic matrix in Sec. II. Two coupled equations are then derived in Sec. III to model the nonlinear wave propagation through a heterogeneous medium by following the analysis for gas bubbles in a fluid by Zabolotskaya and Soluyan.^{38} The frequency-dependent effective medium behavior is then obtained in Sec. IV by solving the coupled model equations at first- and second-order. Only second-harmonic generation is considered in this paper, but analyses for cubic nonlinearity are straightforward extensions of the present work.^{39} The low-frequency limit considered in Sec. V yields an evolution equation in the form of the Korteweg–deVries–Burgers equation. It is worth noting that the theory developed here is not only applicable to metamaterials. Instead, it is relevant for any microscale inclusions described by expanding the strain energy density to at most fourth-order in radial strain, which accounts for up to cubic nonlinearity in the dependence of stress on strain. For example, the present model could be relevant to buckling microspheres if an analytic constitutive relationship were given.

Also of interest in the present work is the frequency-dependent response and resulting harmonic generation, for example, metamaterial inclusion with designed elastic instabilities embedded in a nearly incompressible viscoelastic material proposed in Sec. VI. Such an element is not only tunable with an external activation, but tailorable with the geometric design and choice of material properties and may offer more controllable macroscopic material properties than its more conventional counterparts.

## II. NONLINEAR CONSTITUTIVE BEHAVIOR

Consider a spherical inclusion that has a nonlinear constitutive relationship and is embedded in a viscoelastic matrix material. Of interest is the motion and deformation when accounting for both geometric and material nonlinearity of the inclusion and the matrix material. Finite deformations are therefore necessary to properly account for nonlinear effects. However, the constitutive behavior in nonlinear elasticity is more complicated than the linear elastic stress-strain relationship. To simplify the equations, incremental deformations from an applied pre-strain are introduced to describe both the inclusion and the matrix material. The small-on-large behavior at each increment represents a small perturbation about an imposed pre-strain.

The inclusion of interest is spherical, undergoing purely radial deformation due to uniform pressure imposed on the surface. Assume that the inclusion is at its global equilibrium configuration defined by radius *R*_{0} and pressure *P*_{I0}. The inclusion is then deformed to a new reference state, defined as the pre-strain configuration through application of a static pressure *P*_{I1}. The new configuration is defined by

which is a normalized change in radius from *R*_{0} to the pre-strained radius *R*_{1}. In the present work, subscripts “1” will refer to the pre-strain states. Now, let the inclusion be further deformed a small amount to the total, instantaneous configuration, defined by pressure *P*_{I}, radius *R*, and normalized change in radius

The difference between the total configuration and pre-strain is defined as the incremental deformation field. The incremental dimensionless radius is identical to the infinitesimal strain tensor for an object undergoing purely radial deformation. Therefore, let $\epsilon =\xi \u2212\xi 1$ be defined from Eqs. (1) and (2) and represent the strain increment from the reference state *R*_{1}.

Since the deformation is radially symmetric, the stress tensor is only non-zero on the diagonal. For this particular case, the Cauchy stress tensor is related to the hydrostatic pressure with the relation $\sigma ij=\u2212P\delta ij$.^{40} The same definition holds when considering incremental deformations. Therefore, let the total pressure be expressed as a Taylor series expansion about the pre-strain. The incremental pressure *p*_{I} may then be defined as

where the Taylor series coefficients, *K*_{1}, *K*_{2}, and *K*_{3}, represent the linear, quadratically nonlinear, and cubically nonlinear stiffness moduli local to the pre-strain configuration and are defined as

The stiffness moduli *K*_{1}, *K*_{2}, and *K*_{3} all vary as a function of the applied pre-strain. Although Eq. (3) is obtained by assuming *ε* is a small quantity, the pressure-strain relationship has been expanded to cubic order because the proposed microscale inclusion possesses magnitudes of *K*_{2} and *K*_{3} that may be large enough to render the higher-order terms non-negligible.

It is assumed that the matrix is a nearly incompressible viscoelastic material. In this case, the shear stress in the matrix that interacts with the surface of the inclusion is accounted for through an effective pressure *P*_{M}.^{41} The effective matrix pressure associated with the shear stress is formulated in terms of incremental deformations and evaluated at the surface of the inclusion as

where *μ*_{M} and *A*_{M} refer to the local shear modulus and third-order elastic coefficient, which is a quadratic order term, respectively. These local moduli are defined explicitly based on Ref. 41 as

and *P*_{M1} denotes a pre-strain imposed on the matrix material. The following analysis is limited to $\mu M0/KM0\u226a1$, where $\mu M0$ and $KM0$ are static values of the shear and bulk moduli, respectively, and $AM0$ is a static nonlinear elastic modulus at quadratic order.

## III. EFFECTIVE MEDIUM MODEL

A theoretical model is now derived to study nonlinear propagation of a pressure wave. It is assumed that the only source of acoustic nonlinearity is induced via the frequency-dependent oscillations of the inclusions, which accounts for material and geometric effects both of the inclusion and matrix material local to the inclusion surface. The change in volume of the inclusions due to an incident pressure wave represents a nonlinear source term within the otherwise linearized inhomogeneous wave equation. The derivation most closely follows that developed by Zabolotskaya and Soluyan as presented in Refs. 32 and 38 for bubbly liquids insofar as the main assumption is that the nonlinearity is dominated by the presence of the inclusions. As such, the inclusions appear as a distribution of volume (monopole) sources in an otherwise linear wave equation that is coupled to a nonlinear dynamical model of Rayleigh–Plesset form describing the response of the inclusions. While nonlinear propagation effects within the matrix itself may be worth including in subsequent work, they lie beyond the scope of the present work.

The volume of the entire effective medium $V*$ is defined by the sum of the volumes of its constituents as $V*=NVI+\u2009VM$, where *N* is the number of inclusions, *V*_{I} is the volume of the inclusion, and *V*_{M} is the volume of the matrix material. Assume that the total volumes may be defined as the summation of a pre-strain, denoted with a subscript 1, and an increment. Divide the volume conservation equation by $V*1$ to obtain

The volume fraction at each pre-strain is defined as $\varphi =NVI1/V*1$, such that $1\u2212\varphi =VM1/V*1$. The present analysis is limited to dilute volume fractions whereby $\varphi \u226a1$.

To obtain an approximation for the density of the effective medium, substitute the total and pre-strain densities for the effective medium volumes in Eq. (10). The total density may also be expressed as the summation of a pre-strain state and a small increment that represents an acoustic perturbation, where $\rho =\rho 1+\rho \u2032$. The same substitution for density is made for volume of the matrix. Then, Eq. (10) may be expressed in terms of volume fraction, incremental effective densities of the effective medium and matrix, and incremental volume of the inclusion as

Only the linear terms in the binomial expansion of the density are needed to derive a linear wave equation with a source distribution characterized by the volume oscillations *v* of the inclusions. Binomially expand $(1+\rho M\u2032/\rho M1)\u22121$ and $(1+\rho *\u2032/\rho *1)\u22121$ up to linear order. The quadratic order terms in the expansions, $(\rho *\u2032/\rho *1)2$ and $(\rho M\u2032/\rho M1)2$ are neglected because their contributions are small relative to the strong material nonlinearity assumed for the inclusions. Such approximations have been reasonably made for bubbly liquids and should be valid for the current inclusions.^{32,38}

Since the matrix material is nearly incompressible, $\rho M1\u2248\rho M0\u2261\rho M$ and the local stiffness moduli at each increment are identical to those at the global equilibrium. The density of the matrix is related to the pressure at linear order as $\rho M\u2032/\rho M1\u2261p/\rho McMl2$, where $cMl2=[ KM0+(4/3)\mu M0 ]/\rho M$ refers to the quasi-static, longitudinal sound speed of the surrounding matrix material. The nonlinearity in the equation of state of the matrix material is neglected because it is assumed to be much smaller than the nonlinearity of the inclusions. Common values of the parameter of nonlinearity at quadratic order, expressed as *B*/*A*, are on the order of 5–12 for fluids or fluid-like media.^{42} For example, *B*/*A* = 5 for water, and it is even smaller for air, for which *B*/*A* = 0.4. It is assumed that the strong nonlinearity of the inclusions of interest yields a value of *B*/*A* that is many orders of magnitude larger than that of the matrix material.^{39} In addition, it is also assumed that the matrix is much stiffer than the inclusions. For these assumptions, the contributions to the overall nonlinearity of the effective medium from the parameter of nonlinearity for the matrix material are negligible in comparison to that of the inclusion. Such an approximation is consistent with approximations made for bubbly liquids^{32,35} and is realized most simply for the immiscible mixture law for fluids in the quasi-static limit.^{29} It is therefore reasonable to assume a linearized equation of state for the matrix material.

Furthermore, it is assumed that the effective density is approximately that of the matrix material, such that $\rho *1\u2248\rho *\u2248\rho M$. Although $\rho I$ is not necessarily small compared to $\rho M$, reasonable values of $\rho I$ have little influence on the macroscopic density because $\varphi \u226a1$. With these simplifying approximations, Eq. (11) becomes

Since nonlinear convection terms are neglected in the derivation of Eq. (12), the linearized equations of mass and momentum conservation are combined to obtain a linear wave equation in the form of

Substitution of Eq. (12) into Eq. (13) to eliminate the density yields the following inhomogeneous wave equation:

The change in volume of the inclusion is determined by an augmented form of the Rayleigh–Plesset equation including up to quadratic order nonlinearity. The current form describes the motion of an inclusion as defined in Sec. II and is derived following the formalism used in Ref. 41. This derivation appears in Appendix A, where it is assumed that the surrounding matrix is a nearly incompressible viscoelastic material. The derived dynamics ignore interactions between inclusions, requiring the concentration to be dilute enough that each inclusion behaves independently. Given that $\varphi \u226a1$ is required to derive Eq. (14), it is reasonable that these interactions are ignored.

In terms of change in volume, the dynamic model equation up to quadratic order may be expressed as

Although Eq. (15) is of a similar form to the dynamical equation utilized for bubbly liquids, e.g., in Refs. 32 and 38, the coefficients defined in Eqs. (16)–(20) also account for the inertia and damping of the inclusion, contributions of a nearly incompressible elastic matrix, and the effects of the pre-strain.

The left-hand side of Eq. (15) describes a linear oscillator with undamped natural frequency (squared) given by

and the viscous loss factor

Both the undamped natural frequency *ω*_{1} and viscous loss factor *δ* are functions of pre-strain states and will vary as a function of the deformation through the variables *ξ*_{1}, *R*_{1}, *ρ*_{I1}, *K*_{1}, and *μ*_{M}.

The right-hand side of Eq. (15) contains both the acoustic pressure acting on the inclusion and the quadratic nonlinearity associated with the oscillation of the inclusion. The coefficient of the acoustic pressure *p* in Eq. (15) is

The other two terms on the right-hand side of Eq. (15) account for quadratic nonlinearity. The coefficient *b* represents the nonlinearity in the constitutive relation acting on the surface of the inclusion, and includes both nonlinearity of the inclusion and the matrix material, where

In the linear, fluid limit for which *μ*_{M} = *A*_{M} = 0, *b* is identical to the coefficient of nonlinearity *β* of the inclusion^{39,43} and it is similar to that obtained for bubbly liquids.^{38} Therefore, *b* can be considered to be an effective coefficient of nonlinearity that includes an increased stiffness and quadratic-order nonlinearity due to the effects of the shear stress of the matrix on the surface of the inclusion.

The remaining quadratic term is associated with inertia, the coefficient for which is defined as

The appearance of both *ρ*_{I1} and *ρ*_{M} indicates that the inertial nonlinearity is associated with the moving mass of both the inclusion and the matrix material surrounding the inclusion, respectively. Although only quadratic nonlinearity is retained in the present analysis, accounting also for cubic order nonlinearity in Eq. (15) is relatively straightforward.^{39}

The coupled Eqs. (14) and (15) for *p* and *v* comprise the effective medium model. These equations can be solved together numerically or they can be solved iteratively to describe harmonic generation as in Sec. IV, where a specified form of the inclusion is considered. For acoustic waves with frequencies well below the natural frequency of the inclusion, a single nonlinear evolution equation in terms of the acoustic pressure can be obtained as in Sec. V.

## IV. SECOND-HARMONIC GENERATION

Equations (14) and (15) can be solved iteratively to describe a wave at frequency *ω* that generates a secondary wave at frequency 2*ω*, following a standard expansion for the pressure and volume as done in Ref. 38 and subsequently reviewed in Ref. 32. The linear solution at frequency *ω* provides insight into the frequency-dependent attenuation and dispersion of the effective medium. The second-order solution at frequency 2*ω* describes the strength of the acoustic nonlinearity.

The pressure and volume are expressed in terms of the fundamental wave and its second harmonic as

It is assumed that the harmonic amplitudes *p _{n}* and

*v*are both on the order of

_{n}*ϵ*, where

^{n}*ϵ*is typically defined as the acoustic Mach number and is a small parameter, such that $\u03f5\u226a1$. Substitute Eqs. (21) and (22) into Eq. (14) and equate terms of order

*ϵ*to obtain

Equation (15) then becomes

Now solve Eq. (24) for *v*_{1} and substitute that expression into Eq. (23) to obtain the following equation for the first-order pressure:

where

and

The quantity $c\u03031$ in Eq. (25) is obtained by setting *n* = 1 in Eq. (26), but more generally $c\u0303n$ is the complex phase speed at frequency *nω*. The real phase speed *c _{n}* and attenuation coefficient

*α*at frequency

_{n}*nω*are given by

Now consider second-harmonic generation. The wave equation at order *ϵ*^{2} is obtained by substituting Eqs. (21) and (22) into Eq. (14) to obtain

Likewise, Eq. (15) becomes

where $c\u03032$ is defined by Eq. (26) with *n* = 2 and

The magnitude of the complex function $\beta \u0303(\omega )$ provides a measure of the increase in acoustic nonlinearity introduced by the dynamic oscillations of the inclusions in connection with second-harmonic generation. There are two components of $\beta \u0303$, one involving the coefficient of constitutive nonlinearity *b* and the other involving the coefficient of inertial nonlinearity *d* multiplied by $\omega 2$. As noted by Zabolotskaya,^{36} these components may cancel if constitutive nonlinearity offsets inertial nonlinearity.

For propagation of a progressive plane wave in the +*x* direction with boundary conditions $p1(x=0)=p0$ and $p2(x=0)=0$ and source amplitude *p*_{0}, the solution of Eq. (25) is

where $k\u0303n=n\omega /c\u0303n$ is the complex wave number. Substitution of Eq. (34) into Eq. (32) yields a forced, inhomogeneous ordinary differential equation for *p*_{2}, the solution of which is

The amplitude of the second harmonic varies with propagation distance from the source, frequency, and pre-strain. In the limit of weak dispersion, $|k\u03032|\u223c2|k\u03031|$, the amplitude of the second harmonic grows linearly with *x*. The regime of linear growth is the same as that identified for propagation through bubbly liquids.^{32,38}

## V. LOW-FREQUENCY PROPAGATION

It may also be of interest to consider the low-frequency behavior that corresponds to the response of the effective medium at frequencies below the resonance frequency of the embedded inclusions. In the low-frequency limit, the higher-order inertial effects are much smaller than both the linear terms and the constitutive nonlinearity and may reasonably be discarded. Furthermore, the smallness of the quadratic terms and terms with time derivatives in Eq. (15) permit it to be rewritten recursively in the form of $v=v(p)$ as follows:^{32}

Let the low-frequency sound speed be defined as

Recall that the parameters defining *c*_{*} in Eq. (38) vary as a function of the pre-strain deformation and are therefore not constant. The term within the square brackets in Eq. (38) is the harmonic average of the compressibilities, defined as the inverse of stiffness, with $\varphi \u226a1$ and an effective stiffness of the inclusion $(1+\xi 1)[ K1+(4/3)\mu M ]$. Therefore, *c*_{*} represents the sound speed obtained by a simple mixture law. This approximation is applicable only to media subjected to iso-stress or iso-pressure, such as fluids with zero shear modulus or fluid-like media with shear a modulus much smaller than the bulk modulus.

A one-dimensional evolution equation for progressive plane waves is now derived. A slow scale with $\tau =t\u2212x/c*$ and $x1=\u03f5x$ is introduced to reduce, at leading order, the wave equation in Eq. (37) to an evolution equation in the form of the Korteweg–deVries–Burgers equation:^{32}

where

denotes attenuation,

corresponds to dispersion, and

represents the quadratic acoustic nonlinearity in the low-frequency limit. As previously discussed, *b* represents an effective coefficient of nonlinearity that accounts for both the nonlinearity of the matrix associated with shear stress on the surface of the inclusion, and of the inclusion itself. This parameter is identical to *β*_{I} in the fluid limit. Therefore, Eq. (42) is in the form that would be obtained for an immiscible mixture, with effective stiffness of the inclusion $(1+\xi 1)[ K1+\u2009(4/3)\mu M ]$, effective nonlinearity of the inclusion *b*, and negligible quasi-static nonlinearity of the matrix.^{29}

## VI. EXAMPLE: A NONLINEAR METAMATERIAL INCLUSION

As an application of the theoretical framework developed above, the frequency-dependent effective medium behavior is obtained for hyperelastic inclusions embedded in a nearly incompressible viscoelastic medium that induces a nonlinear elastic response due to buckling of the inclusion structure. Sections VI A–VI C offer a qualitative understanding of the model capabilities and the behavior achievable for one example nonlinear inclusion with elastic instabilities. The inclusion geometry was chosen because the elastic deformation is known from previous studies to demonstrate strongly nonlinear behavior.^{44,45} However, the acoustic response due to propagating waves through a composite medium has not yet been explored, making the chosen design a prime candidate for the present work. The proposed metamaterial inclusion is appealing because the instabilities stem from the geometric design. Therefore, these instabilities are more tailorable and controllable than other buckling inclusions, such as buckling microspheres in a fluid^{23,24} or an elastic^{25–27} matrix material.

It should be emphasized that the theoretical model is valid not only for the chosen example inclusion, but for any class of inclusions, metamaterial or otherwise, that exhibit purely radial deformations and is reasonably described using Eqs. (1)–(6). It may also be used to solve the inverse problem to design, tailor, and optimize the inclusions for specific applications, which should be the focus of future work.

### A. Microscale behavior

The chosen example metamaterial inclusion is shown in Fig. 1(a), where a vertical slice through the center of the sphere reveals a two-dimensional cut in Fig. 1(b). The inclusion is symmetric about both the cut plane and the highlighted symmetry lines. The inclusion features six sets of two orthogonal beams, where the thickness of each beam is a symmetric function of the arc length distance from its center. The beam is thinner in the center and on the ends that attach to the body of the inclusion in order to displace the beam more at those positions than at the thicker segments.^{44} The thickness modulated beam therefore behaves in the same manner as double beam elements that constrain the second buckling mode as described in Refs. 44 and 45. Elastic instabilities are induced by concentrating the imposed pressure at the inclusion-matrix interface at the center of the thickness modulated beam via the pressure transformer, as shown in Fig. 1. The inclusion is assumed to have an initial radius of *R*_{1} = 8.3 mm and is constructed out of nylon. The volume-averaged strain energy density is calculated using the commercial finite element software Comsol, where a displacement-controlled loading condition is imposed on the beam elements shown in Fig. 1 to induce the radial displacement of the inclusion.^{45} As is common practice in effective medium theory in general and in the field of metamaterials in particular, the complicated inclusion geometry depicted in Fig. 1 may be modeled as a homogeneous spherical inclusion with effective material properties subject to the restriction that the radius is much smaller than the wavelengths of interest. Therefore, the present work is also restricted to the long-wavelength limit. The strains, pressures, and incremental stiffness coefficients are then obtained numerically based on Eqs. (2)–(6) using techniques first described in Klatt and Haberman,^{45} which in turn describe the effective behavior of the homogenized spherical inclusion.

The unstable nature of an element with buckled beams is briefly discussed here, but a more detailed description of the modeling of a similar geometry was provided by Klatt and Haberman.^{45} Although Klatt and Haberman considered double beam elements, the behavior is analogous to that of a thickness modulated beam.^{44} The pressure-versus-strain response is shown in Fig. 2(a), where the internal pressure of the inclusion *P*_{I} appears in the solid line and the total pressure on the surface of the inclusion *P*_{I} + *P*_{M} is denoted with a dotted line. The stiffness versus strain appears in Fig. 2(b), where the local linear stiffness of the inclusion *K*_{1} is shown in the solid line and the effective stiffness when also accounting for contributions from the viscoelastic matrix $K1+(4/3)\mu M$ is shown with a dotted line.

The non-monotonic inclusion pressure versus strain is calculated numerically from the strain energy density and imposed displacement obtained from the finite element analysis via a displacement-controlled loading. However, a more realistic scenario for the inclusions is that of pressure loading. State A (solid circle) represents the global equilibrium state. The deformation of the inclusion itself, as opposed to one embedded in a matrix material, is determined by an external loading pressure that changes the internal pressure of the inclusion. As the externally imposed pressure increases, the inclusion compresses until it reaches the local maximum, denoted by the solid triangle at state B. This point denotes a strain state where the local linear stiffness is approximately zero, as evident by state B in Fig. 2(b).

The gray box in both Figs. 2(a) and 2(b) represents the microscale negative stiffness regime and corresponds to strain states where $K1\u22640$. When the inclusion is unconstrained, an additional increase in pressure would cause the inclusion to undergo a large and rapid reduction in radius, where it snaps to a higher pressure state rather than follow the solid black line. The behavior for unloading is similar. This so-called snap-through behavior is due to the elastic instabilities associated with the inclusion geometry. The inclusion configuration denoted by position C is unstable. Therefore, the inclusion cannot exist in that deformed state without a reaction force to counter the decreased pressure associated with an increased strain of the negative stiffness regime. For a fluid matrix where *μ*_{M0} = 0 Pa, the total stiffness on the surface of the inclusion is only that of the inclusion.

In the context of wave propagation, these instabilities are often undesirable on the macroscale. Previous work has indicated that it is possible to obtain macroscopic stability for an unstable microscopic element if it is constrained in the microscale negative stiffness regime, which occurs when the shear modulus of the surrounding matrix material is sufficiently large relative to the stiffness of the unstable element.^{45–51} Macroscopic stability means that the curvature of the strain energy density function is positive for all macroscopic strains, which requires that the linear stiffness of the effective medium always be positive. This is achieved when $K1+(4/3)\mu M>0$. Therefore, the present example includes a matrix with a static shear modulus that is large enough in magnitude to obtain macroscopic stability of the heterogeneous medium at all deformation states, which requires that $\mu M0>(3/4)K1/(1+\xi 1)$ for all deformation states. If the shear modulus were not large enough, then the microscale negative stiffness would induce macroscale negative stiffness and the effective medium would exhibit instabilities manifested as snap-through phenomena similar to that of the inclusion. The present work not only confirms the static stability observed in other works^{45–51} but also offers a model to study the dynamic response at some stabilized static deformation. Macroscopic stability is therefore desired for the present example so that wave propagation can be considered for an inclusion constrained to any deformation state.

This model is also valid for a fluid matrix and will yield similar results for an inclusion with varying stiffness. However, the current case of constrained negative stiffness is considered more interesting because it offers the ability to tailor the macroscopic stiffness in a way that cannot be achieved with conventional, continuous materials and emphasizes the appeal of the proposed inclusion design.

Let the initial volume fraction of inclusion be only *ϕ*_{0} = 5 × 10^{−4}. Although this volume fraction is very small, it is sufficient to induce large nonlinearity on the macroscale as discussed for soft inclusion within either a fluid matrix, such as for bubbly liquids,^{36} or within an elastic matrix with comparable quadratic nonlinearity to the inclusion.^{30,31} The matrix material is assumed to be a soft, rubber-like elastomer that is nearly incompressible,^{27,52,53} with the following material properties chosen to illustrate the behavior of interest: the bulk modulus is *K*_{M} = 2.2 GPa, the density is *ρ*_{M} = 1000 kg/m^{3}, and the shear modulus is *μ*_{M0} = 1.5 MPa. To focus explicitly on the influence of the inclusion on the macroscopic behavior, $AM0=0$ Pa, such that the matrix has only linear material properties. Recall, however, that the effects of geometric nonlinearity have been included in the model and are still accounted for in the analysis presented here.

The effective pressure due to the shear stress of the matrix material is calculated numerically using Eqs. (7)–(9), where the strain increment *ε* in Eq. (7) is that of the inclusion. The total pressure acting on the surface of the inclusion *P*_{I} + *P*_{M} in Fig. 2(a) is now monotonic, which is indicative of macroscopic stability, and has an increased magnitude relative to *P*_{I}. Thus, the value of *μ*_{M0} is large enough to constrain the inclusions for all strain states, including those where $K1\u22640$. An imposed pressure results in loading and unloading that follows the dotted line, such that the inclusion is moving through the negative stiffness regime along the solid line. At pre-strain C, the shear stress of the matrix provides a sufficient reaction force to balance the decreased internal pressure of the inclusion for an increased external pressure, such that the beam elements of the inclusion remain in a deformed configuration that would be unstable in the absence of the restoring force provided by the deformed matrix. If the shear modulus were not sufficiently large, the reaction force would not constrain the inclusion, and the entire heterogeneous medium would exhibit snap-through phenomena associated with macroscopic instability. The effective local linear stiffness when accounting for the shear effects of the matrix becomes $K1+(4/3)\mu M$, as shown in the dotted line in Fig. 2(b), is always positive.

The frequency-dependent behavior is now calculated numerically based on the developed equations in Secs. III and IV for three pre-strain values depicted in Figs. 2(a) and 2(b), where an imposed pressure *P*_{0} constrains the inclusions to the states defined by A (circle), B (triangle), and C (square). Pre-strain A represents the global equilibrium. Although the pressure and strain at pre-strain A are identical, the stiffness differs in Fig. 2(b), where just the inclusion response is shown with the filled circle and the stiffness when accounting for the matrix appears with an open circle. Pre-strain B represents a point approximately at the zero linear stiffness state, but is a small, positive value. Pre-strain C denotes a state where the microscale inclusions are constrained to remain in a configuration that would be unstable without the presence of the matrix. This range is called the microscale negative stiffness regime.

### B. Frequency-dependent effective properties

In the present analysis, let the shear viscosity of the matrix be $\eta M=5\u2009Pa\u2009s$ and the bulk viscosity of the inclusion be $\zeta I=0\u2009Pa\u2009s$. Again, the matrix material property is within the range of values available for soft rubber-like materials in the literature for nearly incompressible media,^{52,53} and the viscosity of the inclusion is neglected to focus on other effects. From Eq. (17), it is apparent that the viscous damping due to the matrix material varies as a function of pre-strain even for a matrix displaying constant shear viscosity. Additionally, the undamped natural frequency from Eq. (16) also varies as a function of pre-strain. For ease of comparing the frequency-dependent properties of the effective medium at each pre-strain, it is advantageous to define a dimensionless frequency relative to one reference frequency. Let the undamped natural frequency at the zero pre-strain state, defined as $\omega 0=\omega 1|\xi 1=0$, serve as the reference frequency for the discussion that follows. The frequency-dependent behavior for all pre-strain cases considered will be analyzed as a function of the dimensionless frequency *ω*/*ω*_{0}.

In Fig. 3, the phase speed as defined in Eq. (26) is normalized by the quasi-static longitudinal sound speed of the matrix *c*_{M}_{l} and shown on a log scale. Well below the resonance frequency for all three cases, the phase speed *c*_{1} approaches the low-frequency sound speed *c*_{*} defined by Eq. (38), which varies as a function of pre-strain. In the low-frequency limit, the effective medium stiffness is lower than that of the matrix material due to the presence of compliant inclusions oscillating in phase with the pressure wave. There is a decrease in stiffness at the pre-strain B (dotted-dashed line) relative to A (solid line), and thus the low-frequency phase speed for pre-strain B is less than A. Similarly, pre-strain C (dotted line) indicates a further decrease in overall stiffness and thus a lower phase speed.

As the frequency increases, the phase speed decreases until the frequency equals the undamped natural frequency at each pre-strain, which corresponds to an antiresonance with minimum phase speed. As anticipated from Fig. 2(b), pre-strain C has the smallest natural frequency, then the natural frequency at pre-strain B, and pre-strain A case represents the largest natural frequency as exhibited in Fig. 3. For frequencies beyond the antiresonance, the phase speed variation becomes large, culminating in a maximum value at a frequency greater than *ω*/*ω*_{0} = 1. The large variation corresponds to frequencies in the range

The upper and lower bounds of the inequality in Eq. (43) are derived by considering the lossless case, where *δ* = 0. The lower frequency limit in Eq. (43) is obtained as $c\u03031\u21920$ in Eq. (26), and the upper frequency bound corresponds to where $c\u03031\u2192\u221e$ in Eq. (26). When damping is included, as in Fig. 3, the magnitude of $c\u03031$ is finite, but may still become large within the bounds defined in Eq. (43). In the high-frequency limit beyond the upper bound of Eq. (43), $c1(\omega )\u2192cMl$.

The frequency range defined in Eq. (43) represents a regime characterized by high attenuation. As the phase speed approaches infinity, which is the limit for an incompressible medium, the propagating waves attenuate more rapidly. Attenuation of the fundamental wave times the wavelength from Eq. (29) is shown in Fig. 4 for pre-strains A (solid line), B (dotted-dashed line), and C (dotted line). The wavelength $\lambda 1=c1/f$ is a function of the dispersive phase speed *c*_{1} from Eq. (28). In the frequency range where the phase speed variance is the largest, $\alpha 1(\omega )\lambda 1$ is non-zero and becomes large, but is very small outside of this frequency range. In Fig. 4, the peak magnitudes occur at the undamped natural frequency of each respective pre-strain. Pre-strain C possesses an attenuation coefficient with the largest width, indicating that the negative stiffness regime attenuates waves over the largest frequency range. Pre-strain B has the second largest bandwidth for high attenuation due to favorable damping properties of a system with zero linear stiffness state. The smallest bandwidth for high attenuation occurs for pre-strain A.

The complex wavenumber $k\u03031$ can be expressed as the summation of real and complex components, such that $k\u03031=\omega /c1\u2212i\alpha 1$. In the absence of loss in the matrix material or inclusion, the phase speed goes to zero at the natural frequency of each pre-strain, and the frequency range defined in Eq. (43) corresponds to no propagation. When $\eta M$ is small but non-zero, the phase speed is bounded but may still become large in magnitude, such that $\omega /c1\u21920$ and $k\u03031\u2248\u2212i\alpha 1$. In this case, waves are purely evanescent. However, if *α*_{1} is too large, the expansion for pressure and volume defined in Eqs. (21) and (22) are no longer valid since the assumptions $|p2|\u2009\u226a\u2009|p1|$ and $|v2|\u2009\u226a\u2009|v1|$ are no longer true. These regimes of high attenuation will be studied in future work.

Second-harmonic generation is used to illustrate the frequency and pre-strain dependence of the acoustic nonlinearity, which may become very large in magnitude. Let the following be defined as the normalized coefficient of quadratic nonlinearity:

where $c\u03032$ and $\beta \u0303$ are the frequency and deformation dependent coefficients of nonlinearity and phase speed of the second harmonic defined in Eq. (26) (with *n* = 2) and Eq. (33), respectively. The low-frequency limits of *c*_{*} and *β*_{*}, as defined in Eqs. (38) and (42), respectively, are evaluated for the zero pre-strain case to provide a constant reference value when illustrating dependence on frequency and deformation.

The normalized coefficient of quadratic nonlinearity $B(\omega )$ appears as a function of frequency in Fig. 5. At all three pre-strains there exists two resonance peaks. The resonances are defined as the frequencies for which Δ_{1} = 0 and Δ_{2} = 0, where *ω _{n}* appears in Eq. (27). For pre-strain A (solid line), one resonance occurs when the frequency equals the undamped natural frequency of the fundamental wave, such that

*ω*/

*ω*

_{0}= 1. Another resonance peak occurs at a lower frequency, where

*ω*/

*ω*

_{0}= 0.5, corresponding to the second harmonic frequency 2

*ω*equaling the undamped natural frequency. The same qualitative behavior occurs for pre-strain B (dotted-dashed line) and pre-strain C (dotted line), but the frequencies of the resonances are shifted downward due to the decrease in undamped natural frequency at those pre-strains.

For the case of gas bubbles in a liquid, analytic expressions for the second-harmonic amplitude were developed by Zabolotskaya and Soluyan^{38} for the cases of *ω* = *ω*_{0} and *ω* = 2*ω*_{0}. Those expressions illustrate that there is growth of the second-harmonic amplitude with propagation distance up to a maximum value, and then a decrease in amplitude beyond that value as the gas bubbles absorb energy. Similar expressions are not derived in the present work, but the behavior for the present nonlinear inclusions will be qualitatively similar to that described for bubbly liquids and readers are referred to Ref. 38 for more information regarding such phenomena.

Although the same number of resonance peaks occurs at all three pre-strains, antiresonance exists only for pre-strain C. An antiresonance occurs at frequencies where $\beta \u0303(\omega )=0$. Thus, the antiresonance occurs when the constitutive nonlinearity at quadratic order offsets the coefficient of inertial nonlinearity, i.e., for frequencies at which $b=3(\omega /\omega 1)2d$ in Eq. (33). When the coefficient of nonlinearity defined in Eq. (33) is negative, *b* is negative and only an imaginary frequency will result in $b=3(\omega /\omega 1)2d$. Since an imaginary frequency is non-physical, no antiresonance exists when $\beta \u0303(\omega =0)$ is negative. The value of *d* is always positive, but no such restriction exists for *b*, which therefore determines the sign of $\beta \u0303(0)$.

From Eqs. (4) and (5), the local quadratic stiffness moduli may be defined as a function of the local linear stiffness moduli through the following relation: $K2\u221d\u2212\u2202K1/\u2202\xi $. Thus, *K*_{2} can be approximated from the slope of the *K*_{1} versus *E*_{I} curve shown in the solid line in Fig. 2(b). The slope of *K*_{1} is negative in the strain region between zero pre-stain and the local minimum of *K*_{1} and positive for negative strains greater in magnitude than the local minimum. Note that *K*_{2} is always negative for positive strains. For the specific example inclusion and matrix properties chosen, the sign of *K*_{2} determines the sign of both *b* from Eq. (19) and *β*_{*} in Eq. (42). However, this will not always be the case; for example *A*_{M} may offset *K*_{2} such that *b* is positive when *K*_{2} is negative. At pre-strain A and B, *K*_{2} is negative and thus $\beta \u0303(0)$ is also negative. However, for the pre-strain C, both *K*_{2} and $\beta \u0303(0)$ are positive and an antiresonance at $\omega =\omega 1b/3d$ is present.

In the low-frequency limit, the quadratic nonlinearity is dominated purely by the quasi-static nonlinearity of the inclusion and matrix defined by *b* in Eq. (19). The normalization in Eq. (44) is defined to ensure that the zero frequency limit of *B*(*ω*) at the global equilibrium is equal to 1, which is the case for pre-strain A. At pre-strain B, the normalized coefficient of quadratic nonlinearity at zero frequency is more than an order of magnitude larger than at pre-strain A, while pre-strain C is approximately 3 orders of magnitude larger than at pre-strain A in the low-frequency limit. Additionally, the peak amplitudes are largest for pre-strain C, second largest for pre-strain B, and smallest for pre-strain C. This illustrates that very large macroscopic acoustic nonlinearity can be induced due to the elastic instabilities on the microscale. Furthermore, the large range of variability possibly highlights the tunable nature of such metamaterial composites. In the high-frequency limit, where $\omega \u2192\u221e$, the amplitude of the coefficient of nonlinearity decays as $\omega \u22124$ and $\beta \u0303(\omega )\u21920$. Therefore, beyond the undamped resonance at each pre-strain, the significance of the acoustic nonlinearity decreases as a function of increasing frequency.

### C. Second-harmonic generation

Last, it is of interest to consider the amplitude of the fundamental and second harmonic as a function of distance from the source. The amplitude of the fundamental wave $| p1 |$ (solid line), as calculated from Eq. (34), and second harmonic $| p2 |$ (dotted line), obtained from Eq. (35), are shown in Pa over two propagation distances in Fig. 6. It is important that the behavior occurs over multiple wavelengths, and the lower *x* axis shows the dimensionless distance relative to the wavelength *λ*_{1}. The upper *x* axis shows the distance from the source in meters. Although the dimensionless propagation distance $x/\lambda 1$ is identical for Figs. 6(a)–6(c), the actual distance from the source differs for each case and the upper *x* axes are not identical for Figs. 6(a)–6(c). Propagation distance for which there is strong harmonic generation is relevant to the design and optimization of the inclusion.

Since the phase speed, attenuation, and coefficient of nonlinearity vary with both frequency and strain, specific frequencies are chosen at each pre-strain of interest to explore the behavior of $| p1 |$ and $| p2 |$. In Fig. 5, it was ascertained that the most prominent nonlinearity occurs at the undamped natural frequency of each pre-strain. However, high attenuation also occurs at the undamped natural frequency. Therefore, frequencies were chosen at *ω*/*ω*_{1} = 0.85 for each respective pre-strain to be far enough below resonance that attenuation does not completely dominate the response.

For pre-strain A, shown in Fig. 6(a), this frequency is *ω*/*ω*_{0} = 0.85 and the corresponding nonlinearity is $|\beta \u0303|=1.7\xd7102$. The amplitude of the fundamental appears relatively constant as it propagates from the source. The magnitude of the second harmonic exhibits an oscillatory behavior with a peak amplitude more than 4 orders of magnitude smaller than $| p1 |$. The periodicity in $| p2 |$ corresponds to the dispersion length $2\pi /| k2\u22122k1 |$, where *k _{n}* =

*ω*/

*c*is the real part of the complex wave number. At distances less than one dispersion length, dispersion effects are negligible.

_{n}For pre-strain B, shown in Fig. 6(b), the chosen frequency of *ω*/*ω*_{1} = 0.85 now corresponds to *ω*/*ω*_{0} = 0.34 with $|\beta \u0303|=7.3\xd7103$. The second harmonic still exhibits an oscillatory response over the period of one dispersion length, which is damped out due to the strong attenuation. The amplitude of the fundamental wave is also attenuated over the propagation distance shown. The peak amplitude of $| p2 |$ is slightly more than 3 orders of magnitude smaller than $| p1 |$. The amplitude of the second harmonic at pre-strain B is larger than at pre-strain A, despite the higher attenuation, because $|\beta \u0303|$ is an order of magnitude larger at pre-strain B than at A.

The lowest frequency considered is for pre-strain C, shown in Fig. 6(c), where *ω*/*ω*_{1} = 0.85 corresponds to *ω*/*ω*_{0} = 0.23 with nonlinearity $|\beta \u0303|=2.2\xd7104$. This pre-strain also possesses the largest attenuation, and both the fundamental and second harmonic are attenuated as they propagate from the source more than at pre-strains A and B. The most attenuation is present at pre-strain C, as evident by the decaying amplitudes of both the fundamental and second harmonic. For $| p2 |$, the oscillatory response is damped out over the propagation distance shown more quickly than for pre-strain B. Pre-strain C also corresponds to the largest coefficient of nonlinearity, which is 1 order of magnitude larger than at pre-strain B and 2 orders of magnitude larger than at pre-strain A, and in turn, an increased amplitude of $| p2 |$ in Fig. 6(c) relative to Figs. 6(a) and 6(b). At $x/\lambda 1=8,\u2009| p2 |$ is approximately 1 order of magnitude smaller than $| p1 |$. Thus, $| p2 |$ from pre-strain C is closer in magnitude to $| p1 |$ than for pre-strains A and B, indicating that pre-strain C generates the strongest second harmonic.

Even though $|\beta \u0303|$ at pre-strain C is very large, $| p2 |\u226a| p1 |$ is still satisfied. For the case of strong dispersion, the phase speed of the second harmonic diverges from that of the fundamental, such that the quantity $1\u2212(c\u03032/c\u03031)2$ in Eq. (35) may become large in magnitude. Thus, weak second harmonic generation relative to the fundamental wave may still be satisfied even for very large $|\beta \u0303|$. For weak dispersion, i.e., where $|1\u2212(c\u03032/c\u03031)2|\u223c1,\u2009\u2009|\beta \u0303|$ is more restricted. With the current inclusion design, it may be possible to induce macroscopic nonlinearity that violates $| p2 |\u226a| p1 |$, and future work should consider higher-order nonlinear effects.

## VII. CONCLUSION

The current work generalized the theoretical dispersion model developed by Zabolotskaya and Soluyan^{38} for bubbly liquids to account for spherical inclusions with a non-negligible moving mass within a viscoelastic matrix material subjected to macroscopic loading conditions. Both the micro- and macroscales are modeled assuming both material and geometric nonlinearity through incremental deformations. In addition to defining the frequency-dependent properties of an effective medium, including phase speed, attenuation, and coefficient of nonlinearity at quadratic order, the fundamental and second harmonic solutions for propagating waves are obtained. Furthermore, a low-frequency evolution equation in the form of the Korteweg–deVries–Burgers equation was derived.

The frequency-dependent wave propagation model also allows identification of three frequency regimes. At frequencies less than the undamped natural frequency of each pre-strain, the system is characterized by high nonlinearity. Harmonic generation is largest at these frequencies. The high-frequency range is where the inclusion motion is considered “frozen” and the effective properties are dominated by those of the matrix. An intermediate range between the low- and high-frequency regimes is where waves are strongly attenuated. The frequencies for each of these regimes depend on the pre-strain imposed on a low volume fraction of an example metamaterial inclusion. This result demonstrates the potential to vary the effective material parameters through mechanical loading of nonlinear acoustic metamaterials. Within the high-frequency limit, the metamaterial requirement that the microstructure be much smaller than the wavelengths of interest will be violated. Thus, the frequency-dependent wave propagation model also provides a means of bounding and verifying the long-wavelength limit, which is paramount to the study of metamaterials.

The present examples chosen reasonably neglected third-harmonic generation by requiring that $| p2 |\u226a| p1 |$. By altering the inclusion design, matrix properties, or source amplitude, for example, cubic nonlinear effects and third-harmonic generation may become more prominent. Extensions to cubic order require accounting for cubic order constitutive and inertial nonlinearity coefficients in the volume formulation of a Rayleigh–Plesset type equation, which is a straightforward augmentation derived in the same way as the coefficient of nonlinearity at quadratic order.

The theoretical model presented in this paper is valid for any nonlinear inclusion that can be modeled following Sec. II and is not limited to the specific inclusion geometry with a thickness modulated beam used as the example here. However, the present model allows the frequency-dependent response for such a metamaterial inclusion to be studied. The frequency dependence of the phase speed, attenuation, and coefficient of nonlinearity all vary as a function of the applied pre-strain. The trends observed for the varying local linear stiffness of the specific inclusion of interest are consistent with behavior previously studied in the literature for a single structure or media that must be periodic. A key difference is that the dynamic response presented here is for a random arrangement of inclusions embedded in a matrix material.

As geometric and material parameters of the inclusion and matrix are altered, the qualitative behavior should remain the same, but the observed magnitudes will vary. For example, the magnitude of the attenuation may decrease, but the frequency range of large attenuation should still be broadest for an inclusion in the negative stiffness regime as compared to a zero linear stiffness state or a state where the constitutive curve is nearly linear. Thus, the theoretical model presented here offers a tool to study the macroscopic frequency- and loading-dependent properties induced by microscale nonlinear inclusions. The elements of interest here offer a means of creating a random effective medium that can yield enhanced attenuation and acoustic nonlinearity that can be tailored for specific applications.

## ACKNOWLEDGMENTS

This work was supported by the Applied Research Laboratories Chester M. McKinney Graduate Fellowship in Acoustics and a Multidisciplinary University Research Initiative from the Office of Naval Research (Grant No. N00014-13-1-0631).

### APPENDIX A: EXTENDED FORM OF THE RAYLEIGH–PLESSET EQUATION

Previous work has focused on the modeling gas bubbles in nearly incompressible elastic^{41} and viscoelastic^{54} media, as well as compliant objects in a nearly incompressible elastic material^{55} using Rayleigh–Plesset type equations. The present derivation combines the methodologies of these three papers, starting with Lagrange's equation,

which is a function of kinetic energy *T*, potential energy *U*, dissipative function Ψ, and the instantaneous radius *R*. The overdot denotes derivatives with respect to time.

The matrix material is assumed to be nearly incompressible; therefore any change in the overall volume of the matrix is neglected. Thus, the kinetic energy is^{41}

where $\rho M=\rho M0$ is the constant, equilibrium density and *R* denotes the instantaneous radius of the sphere. The kinetic energy for a compressible sphere undergoing purely radial deformation is

There are three terms that contribute to the potential energy at the surface of a single inclusion: (i) the internal stored energy of the sphere, (ii) the effective potential energy due to the shear stress of the surrounding matrix, and (iii) the work done on the inclusion by an external force. It is assumed that potential energy is related to pressure such that^{41}

The inclusion pressure is defined in Eq. (3). The effective shear pressure for incremental deformations is defined in Eq. (7) based on Ref. 41. Note that for a fluid, *P*_{M} = 0. The third pressure is the work done by an external forcing function, and may contain time-harmonic and static components.

The dissipative function Ψ is derived from a dissipative function density *ψ*, related through the integral over volume. In the present work, losses are assumed to be linear. The dissipative function density for an isotropic solid with linear losses can be represented by^{56}

where *η* and *ζ* are shear and bulk viscosity coefficients, respectfully, and $\epsilon \u0307ij$ is the infinitesimal strain rate, which neglects nonlinear dissipation terms. For the matrix, the only non-zero velocity is the radial component and only shear viscosity plays a role, where

For the inclusion, only bulk viscosity is present due to the purely radial deformation, where

It is assumed that both $\eta M$ and $\zeta I$ are constants.

An extended Rayleigh–Plesset type equation is obtained by substituting Eqs. (A2)–(A4) and Eqs. (A6) and (A7) into Eq. (A1) and taking the necessary derivatives to obtain

Now it is of interest to derive a volume formulation of Eq. (A8). First, let $VI=VI1+v$, where *v* is a small perturbation. Then, binomially expand the radius in powers of *v* up to quadratic order:

Since the density *ρ*_{I} varies with deformation, it also must be expanded with respect to volume *v*, such that

Note that the volume is expanded only up to linear order because it multiplies terms that are at least linear in order. Thus, for the product to be at most quadratic in order, higher-order terms for density are not necessary. Then, substitute Eqs. (A9) and (A10) into the left-hand side of Eq. (A8) and take the appropriate derivatives to obtain

The loss terms are only retained to linear order because the expressions were derived assuming small strain rates. In terms of volume, the loss is quantified as

For the pressure terms, *ε* must be expanded with respect to *v*, such that

The combined internal pressure of the inclusion and effective shear pressure due to the matrix are defined in Eqs. (3) and (7) to obtain

The example analysis only includes nonlinearity up to quadratic order, and thus *K*_{3} is not necessary for the present work. Extensions considering third-harmonic generation would require maintaining *K*_{3}, as well as expanding all terms to cubic order in volume.^{39} Furthermore, assume that the external pressure can be expressed as the summation of an acoustic pressure wave and an external offset pressure, $Pext=p+P0$. To obtain the desired state $P0=PI1+PM1$.

where *κ*, *ω*_{1}, *δ*, *b*, and *d* are defined in Eqs. (16)–(20).