Generalization of the Zabolotskaya equation to all incompressible isotropic elastic solids

TL;DR

This paper generalizes the Zabolotskaya equation to all incompressible isotropic elastic solids, revealing the impossibility of linear polarization in nonlinear shear waves and analyzing harmonic generation in Gaussian beams.

Contribution

It derives a new coupled nonlinear equation system for shear waves in all incompressible isotropic elastic solids, extending the Zabolotskaya equation and exploring polarization and harmonic phenomena.

Findings

Coupling between anti-plane and in-plane motions cannot be removed in general.

Pure shear beams generate only odd harmonics.

Slight in-plane noise induces second harmonic generation.

Abstract

We study elastic shear waves of small but finite amplitude, composed of an anti-plane shear motion and a general in-plane motion. We use a multiple scales expansion to derive an asymptotic system of coupled nonlinear equations describing their propagation in all isotropic incompressible non-linear elastic solids, generalizing the scalar Zabolotskaya equation of compressible nonlinear elasticity. We show that for a general isotropic incompressible solid, the coupling between anti-plane and in-plane motions cannot be undone and thus conclude that linear polarization is impossible for general nonlinear two-dimensional shear waves. We then use the equations to study the evolution of a nonlinear Gaussian beam in a soft solid: we show that a pure (linearly polarised) shear beam source generates only odd harmonics, but that introducing a slight in-plane noise in the source signal leads to a second harmonic, of the same magnitude as the fifth harmonic, a phenomenon recently observed experimentally. Finally, we present examples of some special shear motions with linear polarisation.

Tables

Table 1. Table 1: Coefficients of quadratic and cubic nonlinearity for two-dimensional shear wave motion in some soft incompressible solids. The values are calculated from elastic wave measurements of μ 𝜇 \mu (initial shear modulus) and A 𝐴 A , D 𝐷 D (Landau coefficients).

Reference	Material	${\beta }_2=1+\frac{A}{4\mu }$	${\beta }_3=\frac{3}{2}\left(1+\frac{A/2+D}{\mu }\right)$
Catheline et al.	phantom gel 1	$-0.78\pm 0.14$
(2003) [2]	phantom gel 2	$-2.98\pm 0.19$
	phantom gel 3	$-0.76\pm 0.03$
Rénier et al.	5% gelatine gel	$-0.43\pm 0.07$	$4.03\pm 0.66$
(2008) [22]	7% gelatine gel	$0.33\pm 0.01$	$2.50\pm 1.28$
Latorre et al.	soft phantom gel 1	$-3.07\pm 0.10$
(2012) [13]	soft phantom gel 2	$0.44\pm 0.07$
	soft phantom gel 3	$0.50\pm 0.05$
	beef liver 1	$-9.03\pm 2.68$
	beef liver 2	$-14.45\pm 5.47$
	beef liver 3	$-17.02\pm 4.57$
Jiang et al.	pig brain (left)	$-0.43\pm 0.11$	$1.41\pm 0.02$
(2015) [12]	pig brain (right)	$-0.46\pm 2.01$	$0.46\pm 0.07$

Equations

$$\mathbf{T}=-p\mathbf{I}+2\left(\partial W/\partial I_1\right)\mathbf{B}-2\left(\partial W/\partial I_2\right){\mathbf{B}}^{-1},$$

$$\text{div }\mathbf{T}=\rho \mathbf{a},$$

$$\text{Div }\mathbf{P}=\rho {\partial }^2\mathbf{x}/{\partial }^{2}t.$$

$$x=X+\tilde{u}(X,Y,t),y=Y+\tilde{v}(X,Y,t),z=Z+\tilde{w}(X,Y,t),$$

$$\mathbf{F}=\left[\begin{array}{ccc}1+\tilde{u}{,}_X& \tilde{u}{,}_Y& 0\\ \tilde{v}{,}_X& 1+\tilde{v}{,}_Y& 0\\ \tilde{w}{,}_X& \tilde{w}{,}_Y& 1\end{array}\right],$$

$$I_1=(1+\tilde{u}{,}_X{)}^2+\tilde{u}{,}_Y^2+\tilde{v}{,}_X^2+(1+\tilde{v}{,}_Y{)}^2+1+\tilde{w}{,}_X^2+\tilde{w}{,}_Y^2,$$

$$\displaystyle I_{2}=I_{1}\,+(\tilde{u},_{X}\tilde{w},_{Y}-\tilde{u},_{Y}\tilde{w},_{X}{\color[rgb]{0,0,0}+}\tilde{w},_{Y})^{2}$$

$$\displaystyle\qquad\quad\,+(\tilde{v},_{X}\tilde{w},_{Y}-\tilde{v},_{Y}\tilde{w},_{X}-\tilde{w},_{X})^{2}-\tilde{w},_{X}^{2}-\tilde{w},_{Y}^{2}{\color[rgb]{0,0,0},}$$

$$\tilde{u}{,}_X+\tilde{v}{,}_Y+\tilde{u}{,}_X\tilde{v}{,}_Y-\tilde{u}{,}_Y\tilde{v}{,}_X=0.$$

$$\tilde{u}=ϵu,\tilde{v}=ϵv,\tilde{w}=ϵw,$$

$$\chi ={ϵ}^{2}cX,\eta =ϵY,\tau =t-X/c,\tilde{p}={ϵ}^{2}p,$$

$$u=c\psi {,}_Y,v=-c\psi {,}_X.$$

$$\tilde{u}={ϵ}^{2}c\psi {,}_{\eta },\tilde{v}=ϵ\psi {,}_{\tau }-{ϵ}^3{c}^2\psi {,}_{\chi },$$

$$I_1=I_2=3+{ϵ}^{2}J,\text{where}J=[(\psi {,}_{\tau \tau }{)}^2+(w{,}_{\tau }{)}^2]/c^2,$$

$$\partial W/\partial I_1={\gamma }_0+{ϵ}^2{\gamma }_{1}J,\partial W/\partial I_2={\lambda }_0+{ϵ}^2{\lambda }_{1}J,$$

$${\gamma }_0={\frac{\partial W}{\partial I_1}|}_{I_1=I_2=3},$$

$${\lambda }_0={\frac{\partial W}{\partial I_2}|}_{I_1=I_2=3},$$

$${\gamma }_0+{\lambda }_0=\mu /2,$$

$${\beta }_2=-\frac{{\lambda }_0}{{\gamma }_0+{\lambda }_0}=-\frac{2}{\mu }{\frac{\partial W}{\partial I_2}|}_{I_1=I_2=3},$$

$${\beta }_3=\frac{3}{2}\frac{{\gamma }_1+{\lambda }_1}{{\gamma }_0+{\lambda }_0}=\frac{3}{\mu }{\left(\frac{{\partial }^{2}W}{\partial I_1^2}+2\frac{{\partial }^{2}W}{\partial I_1\partial I_2}+\frac{{\partial }^{2}W}{\partial I_2^2}\right)|}_{I_1=I_2=3}.$$

$$\psi {,}_{\tau \tau \eta \eta }+2\psi {,}_{\tau \tau \tau \chi }+\frac{2{\beta }_3}{3c^2}(J\psi {,}_{\tau \tau }){,}_{\tau \tau }-\frac{{\beta }_2}{c^2}(w{,}_{\tau \tau \tau }w{,}_{\eta }-w{,}_{\tau }w{,}_{\eta \tau \tau })=0,$$

$$w{,}_{\eta \eta }+2w{,}_{\tau \chi }+\frac{2{\beta }_3}{3c^2}(Jw{,}_{\tau }){,}_{\tau }$$

$$-\frac{{\beta }_2}{c^2}[\psi {,}_{\tau \tau \tau }w{,}_{\eta }-\psi {,}_{\eta \tau \tau }w{,}_{\tau }+2(\psi {,}_{\tau \tau }w{,}_{\eta \tau }-\psi {,}_{\eta \tau }w{,}_{\tau \tau })]=0.$$

$$W=\sum _{i+j=1}^{\infty }{C}_{ij}(I_1-3{)}^i(I_2-3{)}^j,$$

$${\beta }_2=-\frac{2}{\mu }{C}_{01},{\beta }_3=\frac{6}{\mu }\left(C_{20}+C_{11}+C_{02}\right).$$

$$W=\mu \text{tr}\left({\mathbf{E}}^2\right)+\frac{A}{3}\text{tr}\left({\mathbf{E}}^3\right)+D\text{tr}{\left({\mathbf{E}}^2\right)}^2,$$

$${\beta }_2=1+\frac{A}{4\mu },{\beta }_3=\frac{3}{2}\left(1+\frac{A/2+D}{\mu }\right).$$

$$w=\epsilon w_1+{\epsilon }^2{w}_2+{\epsilon }^3{w}_3+{\epsilon }^4{w}_4+{\epsilon }^5{w}_5+\dots ,$$

$$\psi =\epsilon {\psi }_1+{\epsilon }^2{\psi }_2+{\epsilon }^3{\psi }_3+{\epsilon }^4{\psi }_4+{\epsilon }^5{\psi }_5+\dots ,$$

$$w(0,\eta ,\tau )=\epsilon A_1\exp (-{\omega }^2{\eta }^2/c^2)\sin (\omega \tau ),$$

$$\psi (0,\eta ,\tau )\equiv 0,$$

$$w_{1,\eta \eta }+2w_{1,\tau \chi }=0,{\psi }_{1,\eta \eta \tau \tau }+2{\psi }_{1,\tau \tau \tau \chi }=0,$$

$$w_1(\chi ,\eta ,\tau )=W_1(\chi ,\eta )\exp (\text{i}\omega \tau )+\text{c.c.},$$

$$W_{1,\eta \eta }+2\text{i}\omega W_{1,\chi }=0.$$

Full text

keywords: non-linear elasticity, non-linear waves, harmonics, Zabolotskaya equation, multiple scales.

1 Introduction

Anti-plane shear motions are some of the simplest motions around to investigate the nonlinear equations of elastodynamics [10]. This framework is quite general to start with, but its main limitation becomes apparent as soon as we write down the equations of motion, because it turns out that there are very few materials that may sustain a pure state of anti-plane shear in the absence of body forces. In general the equations of elastodynamics reduce to an overdetermined system of partial differential equations. As summarized by Pucci and Saccomandi [19, 20] and Saccomandi [24], this compatibility problem can only be resolved under some special circumstances, for certain classes of materials.

In nonlinear acoustics, the evolution of an initial disturbance with a well-defined direction of propagation has been studied (among others) for longitudinal waves by Zabolotskaya [31] in compressible isotropic nonlinear elasticity and by Norris and Kostek [16] in compressible anisotropic nonlinear elasticity (see the review by Norris [17]).

For transverse waves in compressible isotropic nonlinear elasticity, the situation is more complex because transverse beams are undistorted in the second-order approximation. As shown by Zabolotskaya [31], the strain energy has to be expanded up to the fourth order to obtain a nonlinear motion. The resulting scalar Zabolotskaya equation (Z-equation) describing the propagation of small-but-finite amplitude shear waves in nonlinear elasticity [1, 32] is now a standard model equation [9].

In nonlinear acoustics, there is a widespread belief that the Z-equation is a model equation valid for any constitutive equation, although it is not, simply because the Z-equation is based on anti-plane shear motion which, as we recalled above, is sustainable only by some restricted classes of theoretical models for solids. Indeed, Destrade et al. [5] showed that for nonlinear isotropic incompressible solids, the scalar Z-equation is valid only in the so-called ‘generalized neo-Hookean solids’ and for some other special subclasses of constitutive equations. The strain energy of generalized neo-Hookean solids is restrictive and does not reflect real solids, because it depends on one strain invariant only, instead of the two required for general nonlinear isotropic incompressible solids (see also Horgan and Saccomandi [11]).

Hence, nonlinear shear waves with linear polarization do not exist in general isotropic incompressible solids, only in solids with a very special form of strain energy density which might not be representative of any real-world material. In general, nonlinear shear waves necessarily couple in-plane to anti-plane motion.

Here we derive a generalization of the Z-equation which works for any isotropic incompressible elastic solid (Section 2). As expected, the resulting system of equations couples an in-plane motion to the anti-plane shear motion [21]. We use a multiple scale expansion to derive the system of equations, placing ourselves in the general framework of exact nonlinear hyperelasticity. We then make the link with weakly nonlinear elasticity, and recover the system of equations established by Wochner et al. [29], albeit in a different form, with two equations instead of three, involving a ‘stream’ function, from which the in-plane motion and the transverse component of the motion can be deduced. We also peruse the literature to demonstrate that so far, no real-world incompressible solid has been found such that in-plane motion is decoupled from out-of-plane motion in general. Note that such a decoupling is possible for some special motions, as we show at the end of the paper, see also [11, 19] for a detailed discussion on this point.

In Section 3, we analyse the data recently provided by Espíndola et al. [6] in their investigation of nonlinear waves propagating in gelatine. They showed experimentally that an initially linearly polarized transverse Gaussian beam generated odd harmonics, but they did not remark on the small peak observed at the second harmonic. Here we show that it can be explained by considering that the initial Gaussian beam possesses a small amount of in-plane noise, in a sense clarified in Section 3, of one order of magnitude smaller than the out-of-plane component.

Finally in Section 4 we present examples of anti-plane linear polarisation that can be achieved for some special nonlinear motions.

2 The generalized Zabolotskaya (GZ) system

2.1 Derivation in exact nonlinear elasticity

We denote by $\boldsymbol{F}=\partial\boldsymbol{x}/\partial\boldsymbol{X}$ the deformation gradient associated with the motion $\boldsymbol{x}=\boldsymbol{x}(\boldsymbol{X},t)$, where $\boldsymbol{X}$ and $\boldsymbol{x}$ are the positions in the reference and current configurations, respectively. The left Cauchy-Green strain tensor is $\boldsymbol{B}=\boldsymbol{FF}^{T}$ and the first two principal strain invariants are $I_{1}=\text{tr}\boldsymbol{B}$ and $I_{2}=[(\text{tr}\boldsymbol{B})^{2}-\text{tr}(\boldsymbol{B}^{2})]/2$.

In all generality, the strain energy density $W$ for isotropic incompressible hyperelastic materials is a function of those two invariants only: $W=W(I_{1},I_{2})$, say. Then the Cauchy stress tensor $\boldsymbol{T}$ is [23]

[TABLE]

where $p$ is a Lagrange multiplier associated with the incompressibility constraint, which reads locally as $\text{det}(\boldsymbol{F})=1$.

The equations of elastodynamics in the absence of body forces read

[TABLE]

where $\rho$ is the (constant) mass density and $\boldsymbol{a}$ the acceleration vector. Equivalently, in terms of the nominal stress tensor $\boldsymbol{P}\equiv\boldsymbol{TF}^{-T}$, we can write them as

[TABLE]

In this paper we take $X$ as the direction of propagation and consider the following class of two-dimensional motions

[TABLE]

describing an anti-plane shear motion $\tilde{w}(X,Y,t)$ superimposed onto an in-plane motion with components $\tilde{u}(X,Y,t)$ and $\tilde{v}(X,Y,t)$.

For these motions we compute the components of the deformation gradient as

[TABLE]

in the $\boldsymbol{e}_{i}\otimes\boldsymbol{E}_{j}$ basis, where ($\boldsymbol{E}_{1},\boldsymbol{E}_{2},\boldsymbol{E}_{3}$) and ($\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}$) are the unit vectors along the Cartesians axes ($X,Y,Z$) and ($x,y,z$), respectively, and the comma denotes partial differentiation. From this expression we deduce

[TABLE]

and the incompressibility condition $\det\boldsymbol{F}=1$ reads

[TABLE]

Note that we used this condition in computing $I_{2}$ above.

Our goal is to derive an asymptotic system able to describe two-dimensional shear waves. To this end we introduce a small parameter $\epsilon$ such that the amplitudes can be written as

[TABLE]

where $u$, $v$, $w$ are functions of $X,Y,t$ only, and are of order zero. Then we introduce the following scalings

[TABLE]

where $c$ is the speed of linear (infinitesimal) transverse elastic waves.

We will conduct asymptotic expansions up to $\mathcal{O}(\epsilon^{3})$ and neglect terms of higher orders. We point out that this expansion procedure is the usual expansion of nonlinear acoustics, see Norris [16], instead of the slow time scaling $t^{\prime}=\epsilon^{2}t$, which is often found in Continuum Mechanics [15]. These choices are equivalent, as they correspond to mapping the initial conditions into a source term and vice-versa.

It is now convenient to introduce a “stream function” $\psi=\psi(X,Y,t)$ such that [21]

[TABLE]

Then, using the chain rule we obtain

[TABLE]

and the incompressibility equation (7) is automatically satisfied at order $\mathcal{O}(\epsilon^{3})$. Notice how the scalings (8)-(9) ensure that the shear motion is dominant and the in-plane motion is small in comparison, with amplitudes at least one order of magnitude smaller, see (11).

Further, we find the following expansions for the principal invariants,

[TABLE]

is a non-dimensional quantity. We can then expand the derivatives of the strain energy density as

[TABLE]

where the constants $\gamma_{0}$, $\gamma_{1}$, and $\lambda_{0}$, $\lambda_{1}$ are defined as follows,

[TABLE]

Enforcing continuity with linear isotropic incompressible elasticity, we find that [23]

[TABLE]

where $\mu=\rho c^{2}$ is the infinitesimal shear modulus (the second Lamé coefficient).

Now we introduce the following non-dimensional coefficients $\beta_{2}$ and $\beta_{3}$,

[TABLE]

We then find (details not reproduced here) that to order $\epsilon^{3}$, the equations of motion read

[TABLE]

This is the Generalised Zabolotzkaya system (GZ system), describing transverse waves travelling in any incompressible isotropic solid. Once a solid is specified by a given strain energy density $W$, the constants $\beta_{2}$ and $\beta_{3}$ are computed from the formulas above, and the GZ equations form a system of two coupled nonlinear partial differential equations for $\psi$ and $w$. Once solved, it yields the displacement components $u$, $v$ (from (10)) and $w$ and the motion is described in its entirety.

It was first established by Wochner et al. [29] in the context of weakly nonlinear elasticity, with which we now connect.

2.2 Connection with weakly nonlinear elasticity

In all generality, we may expand the strain energy density $W$ in a Rivlin series, as [23]

[TABLE]

where the $C_{ij}$ are constants. We then find that

[TABLE]

Now, at the same level of approximation, the $i+j=2$ Rivlin expansion of $W$ is equivalent [3] to the following fourth-order Landau expansion of weakly-nonlinear elasticity [33, 4],

[TABLE]

where $\boldsymbol{E}=(\boldsymbol{F}^{T}\boldsymbol{F}-\boldsymbol{I})/2$ is the Green-Lagrange strain tensor, and $A$, $D$ are the third- and fourth-order nonlinear Landau constants, respectively. These constants are linearly connected [3] to the $C_{ij}$, and eventually we find that

[TABLE]

Hence $\beta_{2}$ invokes third-order elasticity only, and $\beta_{3}$, fourth-order elasticity. With respect to stress-strain relationships, $\beta_{2}$ invokes quadratic nonlinearities, and $\beta_{3}$, cubic nonlinearities. These constants were introduced in papers by Zabolotskaya and collaborators [29, 26, 27]. With this connection it is a simple matter to identify our formulation (17) of the equations of motion with that of Wochner et al. [29].

If we were to study one-dimensional plane shear waves (depending on only one space variable), then the derivatives with respect to $\eta$ would vanish from the GZ system (17) and the coefficient $\beta_{2}$ would play no role in the motion: there would be no quadratic nonlinearity for the wave and we would then recover the result of Zabolotskaya et al. [33], with the same coefficient of cubic nonlinearity $\beta_{3}$.

Here we are dealing with two-dimensional shear waves, and the system shows a strong coupling between the three components of the wave in general. Mathematically speaking, there are several ways to simplify the Generalised Z-system (GZ-system) of Equation (17) into decoupled equations, by playing on, and taking special values of the constants $\beta_{2}$ and $\beta_{3}$. For instance [29, 5] by taking $\beta_{2}=0$: in that special case, the GZ-system decouples into an equation for $\psi$ and an equation for $w$. But solids with the special property $\beta_{2}=0$ do not exist in the real world. We show this in Table 1, where we computed the constants $\beta_{2}$ and $\beta_{3}$ from several experimental sources. The conclusion is that in general isotropic incompressible solids, quadratic nonlinearities cannot be separated from cubic nonlinearities when it comes to two-dimensional shear wave motion, contrary to what is postulated in the original paper by Wochner et al. [29], and pursued by several works that followed, see for example [30, 18, 7, 8, 28].

However, for some special motions, linear or plane polarisation in the transverse plane can be decoupled from the longitudinal motion: we present such examples in Section 4. But first we study the effects of the coupling in the GZ-system on the evolution of a shear Gaussian beam and the generation of higher-order harmonics.

3 Gaussian beams in soft incompressible solids

Espíndola et al. [6] recently generated and measured nonlinear shear waves in gelatine. Figure 1 reproduces their data in the case of a high excitation amplitude generating a typical cubically nonlinear shock profile. The transmitted fundamental frequency at 100 Hz has the largest amplitude in the normalised spectrum, and some energy is generated at the (smaller) third harmonic (300 Hz) and at the (minor) fifth harmonic (500 Hz). Something that went un-noticed in the paper is that there is also a peak at the second harmonic (200 Hz), of comparable amplitude to the one of the fifth harmonic at $500$ Hz.

In this section we show that a shear wave beam source condition which is linearly polarised in the $Z-$direction (i.e. $\psi(0,\eta,\tau)\equiv 0$) does not produce a second harmonic in general, even when $\beta_{2}\neq 0$ as is the case for real incompressible solids (in their modelling, Espíndola et al. [6] take $\beta_{2}=0$). On the other hand, we show that if the polarisation of the shear wave beam is slightly misaligned, and therefore in our language ‘noisy’, and allows for some small amplitude variations in the $XY-$plane (i.e. $\psi(0,\eta,\tau)\neq 0$), then a second harmonic is generated, with an amplitude comparable to that of the fifth harmonic.

We consider a regular perturbative solution of the GZ system (17) via a new small parameter $\varepsilon$,

[TABLE]

where $w_{i}$, $\psi_{i}$ are functions to be determined at each order. For a general reference on this perturbation method, we refer to the books by Blackstock and Hamilton [9] or by Naugolnykh and Ostrovski [14].

3.1 Pure anti-plane shear beam

First we take the source data (at $\chi=0$) to be a pure anti-plane Gaussian beam, i.e.

[TABLE]

where $A_{1}$ is a constant and $\omega$ is the frequency of the beam (or equivalently, $c/\omega$ is the effective source radius).

By substitution in the GZ system (17) we obtain at the first order in $\varepsilon$,

[TABLE]

two equations that are independent of the material parameters $\beta_{2}$ and $\beta_{3}$. Clearly we may take $\psi_{1}\equiv 0$ to satisfy the second equation, keeping a pure anti-plane shear beam at first order. We look for a solution to the first equation in the form

[TABLE]

where $\textrm{i}=\sqrt{-}1$, $W_{1}$ is a complex function, and ‘c.c.’ stands for ‘complex conjugate’. Then (23)1 reduces to

[TABLE]

The solution to this parabolic equation, subject to the condition (3.1), is [14]

[TABLE]

where $H$ and $Q$ are the following real functions

[TABLE]

Now we move on to order $\mathcal{O}\left(\varepsilon^{2}\right)$. We notice that our first order solution $w_{1}(\xi,\eta,\tau)$ is such that $w_{1,\tau\tau}=-\omega^{2}w_{1}$; it follows that

[TABLE]

Then the first equation of the GZ system (17) at that order reduces to

[TABLE]

for which we may take $\psi_{2}\equiv 0$ as a solution, maintaining the pure anti-plane shear beam at the second order. Then the second equation of the GZ system (17) at second order reduces to $w_{2,\eta\eta}+2w_{2,\tau\chi}=0$, for which we may also take the trivial solution $w_{2}\equiv 0$.

At order $\mathcal{O}\left(\varepsilon^{3}\right)$ we obtain an equation for $w_{3}$ from the second equation of the GZ system (17), where the coefficient of cubic nonlinearity $\beta_{3}$ now plays a role:

[TABLE]

Substituting our solution at order one, we have

[TABLE]

with $W_{1}$ given in (26). Because of the linearity of this equation, we may look for solutions in the form

[TABLE]

and obtain independent equations for the unknown first and third harmonic amplitude functions $W^{(1)}$ and $W^{(3)}$. For example, the determining equation for the first harmonic of the solution is the forced heat equation,

[TABLE]

A similar forced heat equation determines $W^{(3)}_{3}$ and the solution is complete at order 3 for $w$. Here we see that, as is usual for transverse waves, the harmonic introduced by the initial condition is triplicated. At this same order we have the following equation for the in-plane components, coming from the first equation of the GZ system (17),

[TABLE]

for which we may take $\psi_{3}\equiv 0$, maintaining a pure anti-plane shear beam at third order.

Moving on now to order $\mathcal{O}\left(\varepsilon^{4}\right)$, we find for the first equation of the GZ system (17),

[TABLE]

and now the in-plane motion is excited (we checked by a trivial direct computation that the coupling term in the brackets is not null.) We point out that at this order ($\varepsilon^{4}$), the lowest order where the in-plane motion manifests itself, it is composed of even harmonics only. Then the second equation of the GZ system (17) at fourth order reduces to $w_{4,\eta\eta}+2w_{4,\tau\chi}=0$, which we may solve by taking $w_{4}\equiv 0$.

Here we see that the presence of the $\beta_{2}$ coefficient does eventually lead to a coupling between in-plane and anti-plane wave components. The coupling is weak, in the sense that the anti-plane shear component $w$ of amplitude $\varepsilon$ is coupled to an in-plane function $\psi$ of magnitude $\varepsilon^{4}$.

The interaction between the $\psi_{4}$ term of the in-plane motion and the anti-plane motion $w$ manifests itself at the next order, again thanks to the presence of the $\beta_{2}$ coupling term. Hence at order $\mathcal{O}\left(\varepsilon^{5}\right)$ we find that the first equation of the GZ system (17) gives

[TABLE]

(and we checked that the bracketed term is not zero.) The same structure and sequences are repeated at higher orders. Clearly, this means that only odd harmonics are generated: a pure (hypothetical) anti-plane Gaussian beam cannot generate even harmonics.

3.2 Noisy shear beam

Now let us consider the case of a source data which is not a pure anti-plane shear Gaussian beam; indeed, it is impossible to achieve perfect beam focusing experimentally, and we expect a slightly noisy initial data, which we now model as

[TABLE]

where $B_{2}$ is a constant.

In this case, in the first order of $\varepsilon$ we find the same solution as for the pure anti-plane shear beam initial data, see previous section. At second order, the identity (28) holds and the first equation of the GZ system (17) gives the following equation for $\psi_{2}$,

[TABLE]

as in the previous section, but now with source value (3.2)2. The solution is

[TABLE]

where $H$ and $Q$ were defined in (3.1). It follows by comparison with (24) and (26) that $w_{1}$ and $\psi_{2}$ are proportional to each other:

[TABLE]

As before, $w_{2}=0$.

At order $\mathcal{(}\varepsilon^{3})$, the determining equation for $w_{3}$ reads

[TABLE]

Here it is an easy matter to check that the bracketed term is identically zero. Therefore $w_{3}$ is given by the solution of (30) obtained in the previous section. Then the equation for $\psi_{3}$ is (34), which we solve with $\psi_{3}\equiv 0$.

At order $\mathcal{O}(\varepsilon^{4})$, the first equation of the GZ system (17) determines $\psi_{4}$ as the solution to

[TABLE]

The difference between this equation and the equation (35) obtained at this order when the initial data is a pure anti-plane shear is the additional term proportional to $\beta_{3}$. However, we see from Equation (40) that

[TABLE]

Now by a direct computation (not reproduced here) we find that the solution $\psi_{4}$ of (42) contains the first, third and fourth harmonics.

Then the second equation of the GZ system (17) at the fourth order may be solved by taking $w_{4}\equiv 0$, as in the previous section.

Now we write down the second equation of the GZ system (17) at order $\mathcal{O}\left(\varepsilon^{5}\right)$ for $w_{5}$. It reads

[TABLE]

Our goal here is to show that the solution $w_{5}$ of Equation (3.2) contains the second harmonic. Clearly, it can only arise from the first bracketed term on the right-hand side, which is

[TABLE]

Using (40), we see that this forcing term is proportional to

[TABLE]

Further, taking the solution {\color[rgb]{0,0,0}W}_{1} in Equation (24) and the {\color[rgb]{0,0,0}W}^{(1)}_{3} component of (32), we find the forcing term above to be

[TABLE]

where the ellipsis refers to higher harmonics terms. Our claim is that

[TABLE]

Indeed if the left-hand side above were zero, then integration would yield $W^{(1)}_{3}=\Gamma(\chi)W_{1}$ for some arbitrary function $\Gamma$. Substitution into the forced heat equation (33) and taking (25) into account, we would then obtain

[TABLE]

which is absurd in view of the dependence of $H$ on $\eta$.

The conclusion is that $w_{5}$ contains a second harmonic according to (46), in addition to the expected fifth harmonic coming out of the full solution. This modeling result is perfectly aligned with the experimental data presented by Espíndola et al. [6], see also Figure 1.

We could of course go further in our asymptotic expansion. At higher orders, we would then expect the next even harmonics to be expressed, but with rapidly diminishing amplitudes, as can roughly be inferred from the experimental results in Figure 1.

4 Linear polarisation for special motions

In the final section we explore some avenues to attain linear polarisation, simply by imposing that $\psi\equiv 0$ (and then $u=v=0$ by (10), and only $w$ remains).

With that assumption, the GZ-system (17) reduces to the overdetermined system,

[TABLE]

In general two equations cannot be solved simultaneously for a single unknown, but it might be possible that these two equations are compatible for special classes of solutions.

For example, consider the following ansatz,

[TABLE]

where $f$ is an arbitrary function. Then the first equation in (50) is the trivial identity and the system is no longer overdetermined. Introducing the ansatz in Equation (50)2, we obtain

[TABLE]

Now we specialise the analysis to the cases where $f,_{\eta\eta}$ and f,_{\eta}^{2}+2f{\color[rgb]{0,0,0},}_{\chi} are functions of $\chi$ only, say

[TABLE]

where $f_{1}$, $f_{2}$ are arbitrary functions. In that case, it follows from (53) that

[TABLE]

where $f_{0}$ is an arbitrary function and $c_{0}$, $c_{1}$ are arbitrary constants. The choice $c_{0}=c_{1}=0$ and $f_{0}\equiv 0$ is remarkable, because it reduces Equation (52) to

[TABLE]

where $\hat{F}=F,_{\zeta}$. With the Sionoid-Cates transformation [25],

[TABLE]

where $\chi_{0}$ is a characteristic length, we rewrite the equation as an inviscid cubic Burger’s equation,

[TABLE]

The ansatz in (51) determines a remarkable class of shear waves for which the overdetermined system (50) admits solutions but this is not the only possibility. Indeed, another possible class is given by ‘line solitary wave’ solutions, where we take $w=W(k^{\prime}\chi+k^{\prime\prime}\eta-\omega\tau)$. It is also interesting to point out that the first equation in (50) is identically satisfied in the case

[TABLE]

Then the approximate solution corresponding to the source data (3.1) ‘uncouples’ the anti-plane and in-plane motions at any order.

Acknowledgements

EP and GS have been partially supported for this work by the Gruppo Nazionale per la Fisica Matematica (GNFM) of the Italian non-profit research institution Istituto Nazionale di Alta Matematica Francesco Severi (INdAM) and the PRIN2017 project “Mathematics of active materials: From mechanobiology to smart devices” funded by the Italian Ministry of Education, Universities and Research (MIUR). We are most grateful to Gianmarco Pinton and David Espíndola for sharing the experimental data used to generate Figure 1.