Soft metamaterials with dynamic viscoelastic functionality tuned by pre-deformation

The small amplitude dynamic response of materials can be tuned by employing inhomogeneous materials capable of large deformation. However, soft materials generally exhibit viscoelastic behaviour, i.e. loss and frequency-dependent effective properties. This is the case for inhomogeneous materials even in the homogenization limit when propagating wavelengths are much longer than phase lengthscales, since soft materials can possess long relaxation times. These media, possessing rich frequency-dependent behaviour over a wide range of low frequencies, can be termed metamaterials in modern terminology. The sub-class that are periodic are frequently termed soft phononic crystals although their strong dynamic behaviour usually depends on wavelengths being of the same order as the microstructure. In this paper we describe how the effective loss and storage moduli associated with longitudinal waves in thin inhomogeneous rods are tuned by pre-stress. Phases are assumed to be quasi-linearly viscoelastic, thus exhibiting time-deformation separability in their constitutive response. We illustrate however that the effective incremental response of the inhomogeneous medium does not exhibit time-deformation separability. For a range of nonlinear materials it is shown that there is strong coupling between the frequency of the small amplitude longitudinal wave and initial large deformation. This article is part of the theme issue ‘Rivlin's legacy in continuum mechanics and applied mathematics’.

essentially equivalent to the theory of Simo, based on internal variables and evolution equations in the incompressible limit [46]. As referred to above, small-on-large (SOL) theory is a classical theory developed predominantly over the second half of the twentieth century in the context of nonlinear elasticity, allowing the study of stability under large deformation and also wave propagation in prestressed states [7]. Similarly, theories have been presented for viscoelasticity but the more complex constitutive response makes such a study rather complex and the conclusions depend strongly on the viscoelastic constitutive response chosen. Rivlin was a pioneer in this field, developing a general theory of SOL viscoelasticity with Pipkin for materials with fading memory [47], which allows one to write down a general expression for the incremental stress in a medium subject to large deformation in either a non-steady or steady configuration. Rivlin's subsequent work with Hayes employed the Rivlin-Ericksen constitutive form to lay down the foundations for the study of the propagation of small-amplitude viscoelastic waves in a homogeneously deformed medium [48][49][50].
Zapas and Wineman considered small oscillatory torsional deformations superposed on a uniform extension with the BKZ constitutive model [51], motivated by a classical result of Rivlin [52], who showed that in the elastic case the superposed torsional stiffness is independent of the strain energy density function and can be expressed entirely in a known relationship between the axial force and the axial stretch ratio. Zapas and Wineman showed that viscoelastic oscillations associated with the BKZ model breaks this independence.
For reasons of ease of implementation, as described above, there is great interest in using the Simo/QLV models of nonlinear viscoelasticity. Recall that these theories exhibit strain independent relaxation and therefore the incremental equations give rise to incremental stress-strain relations that are time-deformation separable [53,54]. In a related theory, Morman and Nagtegaal proposed a simplified theory of SOL viscoelasticity which was time-deformation separable [55]. It is known however that although this assumption is often reasonable for homogeneous materials, heterogeneous materials respond rather differently and the so-called effective relaxation function of the composite can be strongly dependent on stretch. Posing a functional form for this relaxation function however is non-trivial. In order to deal with this issue, Kim et al. introduced a so-called static deformation correction factor into the effective relaxation function [53,54], although this manner of introduction and the choice of its functional dependence are rather arbitrary. Here, assuming time-deformation separability for homogeneous phases of a composite rod with simple geometry (piecewise constant properties), we illustrate how the effective incremental response depends on both pre-stretch and time (frequency). Using homogenization theory in this pre-stressed state one can derive the manner in which the effective incremental linear properties depend on the pre-deformation.
It should be clear from the above discussion that Ronald S. Rivlin provided the foundational theoretical work in the area of viscoelastic waves in pre-stressed materials. We hope that the present work illustrates just one aspect of the legacy of Rivlin's work, i.e. the importance of theoretical work associated with wave propagation in pre-stressed nonlinear materials. It shows that his work has had a lasting influence today and impacts on fields such as nonlinear viscoelastic materials and metamaterials that have the potential for far-reaching impact in a number of fields of modern science and technology.
The paper proceeds as follows: in §2, we describe the formulation of linear viscoelasticity in the time and subsequently in the frequency domain. Following this in §3, we illustrate how one can employ homogenization theory in order to derive the effective viscoelastic Young's modulus and in particular that this (complex) modulus can have strong frequency dependence even in the homogenization regime due to its relaxation spectrum. This is particularly the case for soft materials, which tend to have long relaxation times. In §4, we describe the finite deformation of nonlinear viscoelastic homogeneous materials and the relation of the quasi-static long time limit of a uniaxial deformation to the associated hyperelastic deformation problem. In §5, we then formally derive the incremental equation governing longitudinal waves in a homogeneous, thin incompressible quasi-linearly viscoelastic rod that has been subjected to uniaxial tension, inducing large time-dependent deformations. We show that in this homogeneous case the incremental equation has coefficients that are frequency-deformation separable. In §6, we then study inhomogeneous rods to determine the effective viscoelastic Young's modulus in the pre-stressed state, assuming that each homogeneous phase has time-deformation separability. Importantly, it is shown that the effective (homogenized) incremental response does not behave in this separable manner with strong coupling between frequency and stretch.

Linear viscoelasticity
Consider the most general linear viscoelastic constitutive relation, in the form [39,56] whereσ andê are the (second order) linearized stress and strain tensors, respectively, whileĜ(t) is the (fourth-order) stress relaxation tensor, noting that ':' indicates double contraction. One can also write down this constitutive relation in its inverse form, i.e.
whereĴ(t) is the creep compliance tensor. In this paper however attention will be focused on (2.1). The constitutive equations (2.1) and (2.2) are written in the form of a time-convolution, which takes into account any jump discontinuities in the arguments of strain or stress. When the motion starts at the instant t = 0, (2.1) can be written as [56] Integrating by parts, this becomeŝ Let us consider now the constitutive law governing isotropic linear viscoelastic materials. The relaxation tensor then takes the form G(t) = 3κ(t)I 1 + 2μ(t)(I 2 − I 1 ), (2.5) whereκ(t) andμ(t) are two independent relaxation functions and the fourth-order basis tensors I 1 and I 2 have Cartesian components defined as where a superscript + denotes the half-range Fourier transform Experimentally, it is often found that to a good approximation the Poisson ratio ν of a wide class of viscoelastic materials is very weakly dependent on frequency [57]. This motivates the following constitutive relation in the frequency domain as an alternative to (2.8): In (2.11), only one of the material moduli depends on ω: the transform domain Young's modulus E(ω) = −iωE + (ω), with E + (ω) being the Fourier half-range transform of the extensional relaxation functionÊ(t). Furthermore, we see that (2.11) is also in a hydrostatic/deviatoric split format.
In the incompressible limit, κ → ∞, ν → 1/2 and together with e H → 0, (2.8) and (2.11) become, respectively, σ (ω) = −P(ω)I + 2M(ω)e D (ω) (2.12) and where P(ω) is the so-called Lagrange multiplier that accommodates the additional constraint of incompressibility. It is clear that for the consistency of (2.12) and (2.13) in the incompressible limit E(ω) = 3M(ω), which is consistent given that E → 3μ in this limit. Traditionally, Prony series have been used extremely successfully to model the relaxation behaviour of a wide array of viscoelastic materials [58]. These usually take the form, e.g. for the extensional relaxation functionÊ (2.14) where τ j are the associated relaxation times of the medium in question. The simplest models accommodate a single relaxation time τ and a convenient choice of amplitudes giveŝ where E I and E ∞ are known as the instantaneous and long-term Young's moduli and τ is the sole relaxation time of the material. In physical terms, E I is Young's modulus that would be measured from the initial load curve when a material is subjected to rapid extension (on timescales t τ ) up to a strain e max and then held at this strain (termed a 'stress relaxation test'). On the other hand, E ∞ is Young's modulus that would be measured from the load curve when the medium is subjected to a very slow (quasi-static) extension on timescales t τ . The medium is perfectly elastic if E I = E ∞ . The relaxation time τ is obtained by fitting the model to the relaxation test data [58]. We note that for QLV the relaxation functions can be obtained in the linear elastic regime since they are independent of deformation [40]. Employing the form (2.15), we find that where  (1 -f 0 )a q a Figure 1. Illustrating the geometry of the periodic inhomogeneous rod. The periodic cell is of length a q, where q is the cross-sectional lengthscale of the rod. Phase 1 (black) has volume fraction φ 0 and phase 2 (white) has volume fraction (1 − φ 0 ).
The latter inequality arises given that materials are 'glassy' (stiff) at high frequency and 'rubbery' (soft) at low frequency. The shear modulus can also be written in this form in the frequency domain, i.e. M(ω) = M S − iM L .

Linear viscoelastic wave propagation through inhomogeneous rods
Consider now viscoelastic wave propagation in a thin rod whose cross section has characteristic length scale q as depicted in figure 1. Pose the propagation of longitudinal elastic waves of the form subject to lateral stress-free conditions σ 22 = σ 33 = 0 and with time-harmonic dependence ensuring that we work in the Fourier domain. Assume that u 1 (x 1 ) = e ikx 1 where k = 2π/λ is the wavenumber and λ is the wavelength and assume that kq 1 which is the 'thin rod' regime.
In this asymptotic regime, u 2 ≈ −νx 2 u 1 (x 1 ), u 3 ≈ −νx 3 u 1 (x 1 ) [59] and shear strains arising due to u 2,1 and u 3,1 (where we have defined the notation u i,j = ∂u i /∂x j ) are then an order of magnitude smaller than u 1 (x 1 ) given that x 1 scales naturally on 1/k, whereas x 2 and x 3 scale on q. In the incompressible regime of interest here ν = 1/2. Under these conditions, u 1 (x 1 ) satisfies the following governing equation in the frequency domain [60]: where ρ is the mass density. The displacement u 1 (x 1 ) is the fundamental mode that would propagate in the rod at low frequencies, i.e. kq 1 where k = ρω 2 /E(ω). Suppose now that, although the rod has uniform cross-section, it is non-uniform in its mechanical properties, possessing a periodic microstructure with characteristic length scale of O(a). We assume that q a so that the rod can certainly be considered to be thin both with respect to the microstructure and the wavelength. Dispersive effects due to the thickness of the rod can therefore be considered to be negligible and longitudinal waves satisfy (3.2) now where E(ω) and ρ depend on x 1 [60]. With reference to figure 1, we will consider a two phase, periodic composite rod consisting of a repeated periodic cell of length a, within which there are two isotropic, viscoelastic materials. Note that for ease of discussion we shall call phase 1 the inclusion material and denote its volume fraction in the material by φ 0 . In the following section, we will describe the influence of large, quasi-static pre-stress on viscoelastic wave propagation, while in this section we shall describe the situation when there is no pre-stress; this illustrates numerous aspects of viscoelastic wave propagation in heterogeneous media and sets the scene before we extend the study to incorporate the influence of pre-stress. While there has been significant focus on the propagation of viscoelastic waves in layered media [61,62] and the role of damping in viscoelastic phononic crystals [63][64][65][66][67][68], the regime of a thin inhomogeneous rod does not appear to have been discussed before. It therefore provides a convenient asymptotic regime in which to study the influence of viscoelasticity in the large deformation regime, as we will do in later sections once the un-stressed scenario has been considered here.
Referring to (3.2), the equation governing dimensional displacements u 1 (x 1 ) in the inhomogeneous rod is where ρ 0 is the inhomogeneous (and real) mass density in the unstressed state (hence the superposed 0), i.e.
where n ∈ Z. Furthermore E 0 is the (inhomogeneous and complex) frequency domain Young's modulus in the unstressed state, i.e. (3.5) We assume that both materials can be described adequately by a single relaxation time (τ 1 and τ 2 in phases 1 and 2, respectively) Prony series of the form (2.15) and with associated instantaneous and long-term Young's moduli E I 1,2 , E ∞ 1,2 , respectively. Define the non-dimensional variablex = x/a, and the scaled displacementũ(x) = u(ax)/U for some typical displacement magnitude U. Then the governing equation we have defined the non-dimensional parameter = aω/c c = k c a where c c = (E c /ρ c ) 1/2 is a lowfrequency characteristic wave speed given that we have introduced the characteristic elastic Young's modulus, relaxation time and density E c , τ c and ρ c respectively. In particular, we note that the transform domain Young's modulus is piecewise constant, i.e.
and where we have definedẼ ∞ Note that we require k c q 1 in order for (3.6) to be valid; we are yet to choose = k c a but this is limited by the assumption k c q 1, i.e. also cannot be 'too large', although provided q/a is small enough we can take as large as we wish. The scaled Young's modulusẼ 0 r , r = 1, 2 depends on numerous parameters but we stress the dependence onD since this is the dependence on the non-dimensional frequency parameter that generates a low-frequency dynamic response.
Let us now restrict attention to the case of main interest in this article: low-frequency propagation in the regime where frequency dependence can still arise from the viscoelastic behaviour of the phases that comprise the medium. When wavelengths are much longer than the microscale a, we are in the so-called separation-of-scales regime, so that = k c a 1. Using asymptotic homogenization on (3.6), with the assumption that the parameters in the storage and  1.  1), it is straightforward to show that [11,28] the effective Young's modulus in the unstressed state (scaled on E c ) takes on the following harmonic mean form: . (3.10) Recall once again thatẼ 0 1 andẼ 0 2 depend on the non-dimensional frequencyD even for 1. This dependence on ω is solely due to the presence of the (re-scaled) Deborah numberD = ωτ c in the viscoelastic storage and loss moduli. The effective modulus will, of course, also depend on the relaxation time parametersτ r but as noted above we explicitly note its dependence onD only due to our interest on low-frequency dispersive behaviour. What this means is that even for 1, where in the purely elastic case we would be in the homogenization regime with a constant effective Young's modulus, we now have a frequency-dependent effective Young's modulusẼ 0 * =Ẽ 0 * (D), even in this homogenization regime. With phase behaviour as described above, the effective (complex) Young's modulus (3.10) is therefore written asẼ Given that this is the frequency (Fourier) domain effective modulus, the time domain relaxation function associated with extensional deformation can be obtained by Fourier inverting (3.11). This is of practical interest, because it is clear that even in this simple geometry the relaxation function of an inhomogeneous medium is not a simple linear superposition of the relaxation functions of its constituents. Since here we are primarily interested in time-harmonic wave propagation problems, we are mainly interested in the frequency domain response.
In figure 2, we illustrate the effective response by plotting the effective storage and loss moduli for the specific example wheñ where we have taken E c = 10 5 Pa, τ c = 1s, given that these properties are typical of polyurethanes when varying composition and cross-linked densities [69,70], noting in particular that relaxation times for polyurethane can vary between 10 and 1000s. In particular, given thatD = ωτ c = ω, it is clear that there is a strong frequency dependence even in the separation of scales (homogenization) regime. Note the blending of loss modulus amplitudes, and therefore a broader band attenuative effect at low frequency with the two phase medium at φ 0 = 0.5 in figure 2b. The frequency at which the storage modulus increases (usually significantly) is usually termed the glass transition region, referring to the associated increase in stiffness.

Quasi-static deformation of a quasi-linear viscoelastic medium
Note that with the configuration considered above, the effective properties of the periodic rod are fixed once the material properties of the constituent phases and their volume fractions are specified. One can envisage many situations where one would like to modify the effective properties of a material in situ. As discussed in the Introduction, over the last decade a significant amount of research has been conducted into the use of nonlinear elastic materials as a means to modify the incremental dynamic response of a medium. Here we study the influence of viscoelastic effects on longitudinal wave propagation in thin rods when the rod is subject to uniaxial deformation via a tensile stress T. As we shall describe in §4c, we shall assume that the initial large deformation of the thin rod is quasi-static and piecewise homogeneous. Strictly, given that the rod is inhomogeneous, complex local deformations would result close to interfaces. However, these will be confined to small regions within the vicinity of the interfaces and would certainly not have a dominant response in the low-frequency regime of interest here. Before analysing the structure of the inhomogeneous rod let us consider what it means to be in the quasi-static regime of a viscoelastic material that has been subjected to finite deformation. For convenience, we do not append a hat to time varying variables in this section as we have done above, because we are not envisaging taking Fourier transforms of these variables -they are associated with the large deformation problem.

(a) Quasi-static deformations
We shall study nonlinear viscoelastic media that behave in a manner described by the QLV theory [40,45,71]. The general constitutive expression for anisotropic media takes the form where Π is the second Piola-Kirchhoff stress with superscripts 've' and 'e' referring to the viscoelastic and elastic stresses, respectively, and G is a fourth-order reduced relaxation tensor, where reduced refers to the fact that it is non-dimensional, unlike the relaxation tensor in linear viscoelasticity, which has dimensions of stress, see e.g. (2.3). The stress Π e is the (instantaneous) elastic stress, derived from a strain energy W as is standard in nonlinear elasticity. The form prescribed in (4.1) preserves objectivity and provides a balance between realistic modelling and ease of implementation in computational simulations. Of specific note and of importance in what is to follow in this article, for isotropic, incompressible QLV materials, relaxation is independent of strain. Recall that the Cauchy stress T is related to Π by where F is the deformation gradient while the superscript 'T' denotes its transpose and J = detF.
Restricting attention now to the scenario of interest in the present paper, isotropic incompressible media, equation (4.1) takes the form [40] T where P(t) is the Lagrange multiplier associated with satisfying the constraint of incompressibility. D(t) is the reduced deviatoric scalar relaxation function, which without loss of generality is defined subject to the condition that D(0) = 1 and furthermore It should be noted that is the deviatoric part of the Cauchy stress. C = F T F is the right Cauchy-Green strain tensor and W i , (i = 1, 2) is the derivative of the strain energy function W with respect to the invariants I i , (i = 1, 2) of C, and furthermore I 3 = J 2 = det C = 1 due to the constraint of incompressibility. In a homogeneous, isotropic and incompressible material subjected to a simple uniaxial homogeneous extension in the uniaxial direction x 1 , a point in the undeformed configuration, prescribed by Cartesian coordinates (X 1 , X 2 , X 3 ) moves to in the deformed configuration, where λ is the uniaxial stretch along the x 1 -direction. Here we suppose that such a deformation has arisen given the uniaxial stress distribution of the form where T(t) is therefore the uniaxial tensile load. Under the quasi-linear viscoelastic stress-strain law (4.3), it is straightforward now to determine the relationship between T(t) and λ(t), starting by writing down the expressions for the diagonal stresses, i.e.
T ve (4.10) Thus, by substracting (4.10) from (4.9), we obtain where from (4.4) and Π e D22 = 2 At this point, let us be more specific about the type of material under study. We consider Mooney-Rivlin materials, with strain energy function of the form [72] W = μ I 2 where μ ∞ is the long-term shear modulus and τ denotes the relaxation time. The scaling used in the reduced relaxation function with μ I factored out in order to appear in the strain energy function (4.14), ensures that D(0) = 1. However, if one wished, μ ∞ could be employed in the strain energy function in place of μ I . This would yield a different scaling of the reduced relaxation function, leading to D(0) = μ I /μ ∞ , which would then be a factor appearing in (4.3). Finally, equation (4.11) becomes (b) Relaxation and creep tests: the equivalence of the quasi-static limit with hyperelastic theory We now describe how the large time limit of the QLV theory, which we shall term as the quasistatic limit, is equivalent to a purely hyperelastic deformation having the same strain energy function as that used in the QLV analysis but with the μ I in (4.14) interchanged with μ ∞ given that this is associated with long-term deformation. Note that one can impose the stretch λ(t) and determine the resulting time-dependent stress T(t) or vice-versa. The latter is more challenging mathematically given that (4.16) is a nonlinear integral equation for the unknown λ(t), although an efficient method to determine the solution of this equation was introduced in [40]. As an example of the resulting deformations, suppose that the relaxation time of a homogeneous material is τ = 2 s and therefore define τ c = τ andt = t/τ . We can consider imposing either the deformation in the 'ramp' form as depicted in figure 3a or we can pose the non-dimensional tensionT = T/μ I in the ramp form as depicted in figure 3c. Imposing stretch as in (4.18) gives rise to the classical relaxation test, with the resulting tension predicted in figure 3b for three different ratest = 10, 20, 50, noting that this is a time scaled on τ . In this instance, the tension follows directly by computing the integrals on the right-hand-side of (4.16) at successive times. Of particular note is that the slower the deformation the smaller the relaxation effect, although regardless of this,T tends to a fixed value ast → ∞. Analogously in figure 3d, the respective stretch predictions are plotted when the ramp stress (4.19) is imposed witht = 10, 20, 50, by solving the integral equation that results from (4.16) as described in [40]. In all plots here a Mooney-Rivlin strain energy function is taken with γ = 0 and we also took μ ∞ /μ I = 0.3. We note that the horizontal 'dot-dashed' lines in figure 3b,d are the long-term viscoelastic deformations and can be determined from standard hyperelastic theory via the relations (see e.g.  where B = FF T and β 1 = 2(W 1 + I 1 W 2 ) and β 2 = −2W 2 , (4.21) noting that for the calculation leading to the horizontal dot-dashed lines in figures 3b,d, the Mooney-Rivlin strain energy function (4.14) has been employed, but with μ I interchanged with μ ∞ , i.e.
where here, γ = 1/2 recovers the neo-Hookean model W NH . This interchange of moduli in the case of hyperelasticity is consistent given that we wish to determine the long-term viscoelastic deformation, and its accuracy is confirmed by the agreement of the long-time QLV and hyperelastic results depicted in figure 3b,d. In particular, for the uniaxial elastic deformation described by relations (4.20) give T = −P + β 1 λ 2 + β 2 λ 4 and 0 = −P + β 1 λ + β 2 λ 2 (4.24) so that For the scenario depicted in figure 3b when γ = 0, the long-term elastic limit is

(c) Quasi-static pre-stress of an inhomogeneous incompressible rod
We now consider the nonlinear deformation of an inhomogeneous incompressible rod which has periodic microstructure as in §3. We assume that each phase of the rod is quasi-linear viscoelastic but is subjected to a slow ramp deformation and we are interested only in the large time limit, which as we have shown above is equivalent to the hyperelastic theory with an appropriate strain energy function. In the deformed state, the length and width of the periodic cell of the rod will change to a and q , respectively (deformed from a and q in the unstressed state) and hence if we non-dimensionalize these on a we find that in the deformed configuration, the nth periodic cell is defined by the domain where η = q /a 1 and noting that we have used a tilde on spatial variables given that we have non-dimensionalized.
The inhomogeneous rod will be subjected to uniaxial tension as in (4.8) with T 11 = T and lateral stress-free conditions. Since the cell is long and thin (η 1), we neglect the effects of necking close to the phase interfaces and assume a simple planar deformation where each crosssectionalx 2x3 plane within each phase deforms identically. This appears to be a reasonable approximation for η 1 and will certainly capture the leading-order behaviour. This assumption is similar to the approximation made in [37,74,75], for example (although the nature of the approximation is not discussed there). As a result of these assumptions together with the constraint of incompressibility, the resulting deformation in the rth phase (r = 1, 2) is described bỹ where Λ r and 1/ √ Λ r are the principal stretches in the longitudinal and lateral directions in the rth phase, r = 1, 2. The constants γ (n) r are associated with a rigid body motion in the rth phase and nth cell. Such rigid body displacements are required in order to satisfy the continuity of displacement boundary conditions on the interfaces between the phases, as will be described below.
Given continuity of tension across phase interfaces and assuming T is imposed, the stress-stretch relations for each phase, derived from (4.25), yield an equation for the principle stretch Λ r in phase r = 1, 2. The resulting Lagrangian forms of displacements are where as noted above, the constants γ (n) r are deduced by insisting on continuity of displacements atX 1 = n, n + φ 0 , i.e. (4.33) and These are coupled with the assumption that (without loss of generality) γ (0) 2 = 0, which gives a fixed reference. We note again that the specific values of γ (n) r are merely rigid body displacements required in order to satisfy the continuity of the body. The nonlinear deformation above now serves as an equilibrium state from which to perturb via the consideration of superposed longitudinal waves.

Longitudinal waves in pre-stressed elastic and quasi-linear viscoelastic thin rods
Before considering wave propagation in the pre-stressed inhomogeneous medium, let us consider how one determines the incremental equations in a homogeneous quasi-linear viscoelastic rod. We first summarize the existing theory for hyperelastic rods before deriving the theory in the QLV scenario.

(a) Incremental elastic waves
The classical theory of 'small-on-large' allows one to write down the equations governing smallamplitude elastic vibrations superposed on nonlinear elastic deformation and here we briefly summarize the theory for the case of hyperelastic materials [7]. In the context of a homogeneous incompressible thin rod, consider an initial uniaxial deformation as in (4.23) and (4.24), before assuming the presence of a superposed longitudinal wave of the form (3.1), working in the same asymptotic regime as that following (3.2) so that we neglect shear strains and lateral displacements are such that u 2,2 = u 3,3 = −νu 1,1 = − 1 2 u 1,1 (5.1) for the incompressible case studied here. Assuming that the incremental longitudinal stress is Σ we therefore have the incremental relations (see e.g. (3.16) in ch. 1 of [7]) where P is the Lagrange multiplier derived from (4.24) and p is the incremental Lagrange multiplier. The incremental moduli A ij are defined as  (5.6) and the incremental stress-strain relationship is then given by where ρ = ρ 0 given that the medium is incompressible. For the specific Mooney-Rivlin strain energy function, we find that and we note that this agrees with the form deduced in the same asymptotic (thin rod) limit in (3.10) of [76]. When λ = 1, A = 3μ ∞ = E ∞ .

(b) Incremental viscoelastic waves
Although the small-on-large theory associated with hyperelastic materials is well established, the analogous theory associated with QLV materials is not. We shall therefore derive such a theory, tailored to the specific configuration of interest here. A more general theory shall be developed elsewhere in the future. In addition to the large time-dependent uniaxial deformation (4.7), let us consider a separate deformation that is 'close' to the deformation (4.7), denoted byx = χ (X). We define the difference between position vectors in the deformed configurations as The total deformation gradient can be written as where γ = grad u and Grad and grad denote the gradient operator with respect to X and x, respectively. Furthermore,J = detF = J + tr(γ )J to first order in the perturbation and hence incompressibility requires tr(γ ) = u k,k = 0 as expected. Let us write the total Cauchy stressT and Piola-Kirchhoff stressΠ in the form T = T + τ andΠ = Π + π (5.13) so that τ and π denote the difference in Cauchy and second Piola-Kirchhoff stresses between the two deformation states. From (4.2), 14) The total viscoelastic Cauchy stress can therefore be written as where and where p is the incremental lagrange multiplier such thatP = P + p, and is determined by imposing τ ve 22 = 0. Note that Next, denoting the divergence operator with respect to x (x) by div (div) and neglecting inertia associated with the initial finite deformation, the conservation of momentum equation takes the form divT ve = ρü, (5.19) where ρ = ρ 0 , recalling that ρ 0 is the mass per unit of volume in the undeformed configuration. At leading order, (5.19) is div T ve = 0, (5.20) whilst upon defining Σ ve = τ ve − γ T ve , the equation governing the perturbation is div Σ ve = ρü. (5.21) Assume that the incremental displacement field takes the form (consistent with the longitudinal waves described above) where H denotes the Heaviside function, indicating that this perturbation is 'switched on' at time t = t 1 , where t 1 0 and Now consider the QLV theory for the specific case of a Mooney-Rivlin material defined by (4.14) and for the uniaxial deformation defined in (4.7). It can be shown that the resulting incremental viscoelastic stresses take the form (2 + λ 3 (s))λ 3/2 (t) e −iωs − λ 3/2 (s)(λ 3 (s) − 1) e −iωt u 2,1 ds (5.25) and the shear stresses are therefore set to zero upon neglecting shear strains. The contribution to (5.21) is then only due to Σ ve 11 so that the equation of motion reduces to in the domain of integration. Further, consider t t 1 and t τ , so that any transients on the lefthand side of (5.26) decay to zero. In this quasi-static regime, the equation governing perturbations is then where A is as defined for a Mooney-Rivlin material in (5.10) and we have also defined where we recall thatẼ = E I /E ∞ = μ I /μ ∞ in the incompressible regime and D = τ ω is the Deborah number. Finally, with reference to §3, as it may be convenient to scale time or frequency on a characteristic time τ c one can also define (analogously to (3.9)) R S (D) = 1 +Ẽ(Dτ Compare (5.31) with the incremental case without viscoelasticity, (5.9). Notably, in this quasistatic regime one can therefore conclude that homogeneous QLV materials lead to small amplitude perturbations that exhibit frequency-deformation separability [53,54]. The deformation dependence is purely associated with the incremental modulus A(λ) which can be determined from the hyperelastic theory and the frequency dependence can be determined from the unstressed scenario. This is a consequence of the assumptions on t t 1 t and the QLV model. In figure 4, we illustrate this by showing that pre-stress merely scales the storage and loss moduli associated with longitudinal vibrations of a QLV material that is under uniaxial tension.

Effective low-frequency viscoelastic waves in a pre-stressed inhomogeneous rod
The equations governing the incremental response of a thin homogeneous rod have been derived in both the hyperelastic and QLV scenarios. We now illustrate the effective low-frequency response of an inhomogeneous incompressible thin rod of the same geometry as that described in §3 but now under the quasi-static finite uniaxial deformation as described in §4c. We first consider the response in the absence of loss before moving on to the problem of main interest in this paper. Given that the material is incompressible ρ 0 = ρ throughout and furthermore the volume fraction of each phase remains fixed during deformation, i.e. φ 0 = φ.

(a) Elastic waves
The equation governing longitudinal elastic waves in each phase is (5.9) with phase modulus A r and density ρ r . Non-dimensionalizing as per §3, the governing equation for incremental longitudinal waves is then given by Since the form of (6.1) is entirely equivalent to that in (3.6), the homogenization process can proceed in precisely the same manner so that we can determine the effective Young's modulus, stressing that this holds in the elastic regime but now in a pre-stressed state. This allows us to determine the influence of pre-stress on the effective wavespeed of the superposed wave. We can employ formula (3.10) for the effective Young's modulus with the only modifications being that we now use the incremental modulus in place of E, noting that (for now at least) the incremental Young's modulus is purely elastic. Therefore, the effective incremental Young's modulus (scaled on E c ) for a two-phase periodic composite rod is straightforwardly defined as The effective mass density in the stressed state is simply the arithmetic mean: which we note is scaled on ρ c .
We consider examples where phases are Mooney-Rivlin as already defined in (4.22) as well as the Yeoh model, i.e.
where α > 0. In this Yeoh strain energy case, the incremental modulus A as defined in (5.8) is Upon denotingμ r = μ r /E c , r = 1, 2, we set the dimensionless elastic constants and we also note again that in the incompressible scenario E ∞ = 3μ ∞ . In figure 5, we plotÃ * as a function of φ, the volume fraction of phase 1, for a range of strain energy function combinations. We note in particular the sensitivity of the behaviour to the choice of strain energy. For the parameter regime illustrated, it is noteworthy that for all cases except one, pre-stress tends to stiffen the effective material response as additional phase 1 material is added to the inhomogeneous medium. The anomaly is the combination of Mooney-Rivlin (phase 1) and Yeoh (phase 2) for which pre-stress appears to have a softening effect for increasing φ. This indicates that the Yeoh model becomes relatively stiffer (compared to Mooney-Rivlin) for a givenT and therefore additional Mooney-Rivlin material only serves to soften the response. This macroscopic effect arises due to the combination of incremental moduli (5.10) and (6.5) in the harmonic mean form (6.2).
In figure 6, for the same range of strain energy combinations as in figure 5, we illustrate the specific response to tension for the case when φ = 0, 0.5 and 1 illustrating the ability to tune the incremental Young's modulus extremely effectively by combining nonlinear materials with predeformation. We limit our discussion here given that this effect is well known and move on to the case of dominant interest-the influence of viscoelasticity.

(b) Viscoelastic waves
Proceeding now to the viscoelastic scenario and for the regime of interest here, we employ (5.31) as the equation governing longitudinal wave propagation in each phase of the inhomogeneous rod.
Employing the notation used above, we can therefore define the pre-stressed incremental (frequency dependent) viscoelastic modulus in the rth phase as in (5.34), i.e. E r = A rRr (noting now the absence of the superscript '0' as was employed for the unstressed modulus) and thus write (5.31) in an appropriate form for the inhomogeneous rod  noting that now E depends on x 1 and also on pre-stress. This equation is then nondimensionalized by choosing an appropriate E c and ρ c , to yield where = k c a is defined as above. We can now determine the effective properties given that we have the equation in our canonical form. This allows us to write down the effective Young's modulus in the transform domain in the form (1 − φ)Ẽ 1 + φẼ 2 (6.9) and we note that deformation and relaxation effects are clearly not separable in (6.9). This simple case serves to illustrate that even for the simplest of microstructural geometries one finds that  for some deformation dependent function C(λ). This form has been employed as a phenomenological expression for the form of effective strain-dependent relaxation of inhomogeneous materials by e.g. [53,54]. Expression (6.9) is exact in the asymptotic regime of interest here for low-frequency propagation in thin rods. In figure 7a, we illustrate the stretch dependence in the case of an inhomogeneous rod whose phases behave quasi-statically as neo-Hookean materials with long-term moduli defined in (3.12), noting that μ ∞ = E ∞ /3 in the definition of the strain energy functions. Viscoelastic properties of the phases are as defined in (3.12). We plot the effective incremental loss and storage moduli associated with the unstressed and stressed (T = ±10) curves. It should be noted that the pre-stress stiffens the effect across a broad range of frequencies (as should be expected) but more importantly, this effect is nonuniform across the frequency spectrum, leading to the conclusion that frequency-deformation separability does not hold for the inhomogeneous material considered here. In figure 7b, we plot the value of the Deborah numberD at which the loss modulus achieves its local maximum close tõ D = 0.1 with varyingT and for increasing phase 1 shear modulus:μ ∞ 1 = 5 (green),μ ∞ 1 = 10 (red) andμ ∞ 1 = 50 (blue) whenμ ∞ 2 = 3. These results are in stark contrast to the homogeneous bar case as depicted in figure 4 where the spatial homogeneity ensures that the frequency-deformation separation remains.
To further illustrate the lack of frequency-deformation separability, the effective incremental loss and storage moduli are plotted in figure 8 for the cases when the quasi-static nonlinear elastic behaviour of phase 2 is of Yeoh type and phase 1 is (a) neo-Hookean and (b) Mooney-Rivlin with γ = 0, together with the parameter set (3.12). The striking effect of the lack of frequencydeformation separability is particularly evident: this effect is emphasized here when the phases possess different strain energy functions.

Conclusion
A methodology and formulation has been presented in order to determine the effective incremental dynamic (loss and storage) moduli associated with pre-stressed inhomogeneous nonlinear materials. In particular, it is noted that frequently relaxation times of such soft materials are long and therefore these media are very often frequency dependent even at 'low' frequencies where classical homogenization methods would normally give properties that are independent of frequency. It is shown that the effective incremental moduli are very sensitive to the choice of strain energy functions of the constituent phases. Furthermore, and perhaps most importantly, it is shown that even when constituent phases are assumed to be time-deformation separable in their individual constitutive response, the resulting effective behaviour of the thin inhomogeneous rod has strong frequency-deformation (and hence time-deformation) coupling. We note that at higher frequencies, such periodic materials would act as soft-phononic crystals. One could extend the above analysis to that context but the fully dispersive rod theory would then have to be developed, taking into account lateral inertia and other effects that become important at higher frequencies.
Data accessibility. All results can be obtained from information provided in the manuscript. Authors' contributions. Both authors contributed equally. They have both read and approved the manuscript. Competing interests. The authors declare that they have no competing interests. Funding. W.J.P. acknowledges the EPSRC for funding his research fellowship (EP/L018039/1).