Antiplane Stoneley waves propagating at the interface between two couple stress elastic materials

We investigate antiplane Stoneley waves, localized at the discontinuity surface between two perfectly bonded half-spaces. Both half-spaces are elastic linear isotropic and possess a microstructure that is described within the theory of couple stress materials with micro-inertia. We show that the microstructure deeply affects wave propagation, which is permitted under broad conditions. This outcome stands in marked contrast to classical elasticity, where antiplane Stoneley waves are not supported and in-plane Stoneley waves exist only under very severe conditions on the material properties of the bonded half-spaces. Besides, Stoneley waves may propagate only beyond a threshold frequency (cuton), for which an explicit expression is provided. For a given frequency above cuton, this expression lends the admissible range of material parameters that allows propagation (passband). In particular, significant contrast between the adjoining materials is possible, provided that Stoneley waves propagate at high enough frequency. Therefore, micro-inertia plays an important role in determining the features of propagation. Considerations concerning existence and uniqueness of antiplane Stoneley waves are given: it is found that evanescent and decaying/exploding modes are also admitted. Results may be especially useful when accounting for the microstructure in non-destructive testing (NDT) and seismic propagation.


Introduction
The quest for proving the existence of new types of localized waves, similar in nature to Rayleigh waves occurring at a free surface, begun shortly after the discovery of these by Lord Rayleigh [32]. Love [18, p.165] investigated the possibility of waves propagating at the free surface of a layer perfectly bonded to a half-space, in an attempt to explain the problematic (from the theoretical standpoint) appearance of shear horizontal Rayleigh waves in seismograms [22]. Later, Stoneley [30,31] investigated the existence of waves localized at the surface of a discontinuity between two materials. As he points out "Whereas, however, Prof. Love's problem is concerned with a disturbance confined chiefly to the free surface, the present paper deals with a wave motion that is greatest at the surface of separation of the two media" [31]. Indeed, these waves go under the name of generalized Rayleigh waves or, quite fittingly, Stoneley waves. Stoneley concludes that "we can definitely assert that when the wave-velocities are not too widely different for the two media, a wave of the Rayleigh type can exist at the interface". This finding is somewhat surprising, in that it allows localized waves to exist only at "weak" boundaries, whereas the exact opposite could be expected. At the time of Stoneley, the main motivation behind the investigation was geophysics: namely to determine whether seismic energy could escape at the Gutenberg-Wiechert boundary between the Earth's mantle and the core. Indeed, existence conditions were supposed (by Love and Stoneley) to correspond to Wiechert conditions, expressing equality of the shear wave speed across the surface of discontinuity. Precise quantification of the range of existence of Stoneley waves came much later, by Scholte [28]. Scholte conditions are very severe, so much that "the existence of such modes is the exception rather than the rule" [2]. Indeed, as pointed out in Owen [25] and Hsieh et al. [13], they are satisfied by merely 31 combinations among 900 isotropic materials.
Existence and uniqueness of Stoneley waves have been investigated in the monograph by Cagniard [4,Chap.4] and by Chadwick and Borejko [5]. (Notably, Chadwick has been Robert Stoneley's research student.) Extension to anisotropic materials, which admit generalized Stoneley waves, was given by Lim and Musgrave [17] and by Barnett et al. [2].
Recently, the role of material microstructure has attracted considerable attention in connection with wave propagation and related phenomena [9-11, 21, 22, 24, 29]. One way of encompassing material microstructure into the models is by means of polar theories, among which couple stress (CS) theory is perhaps the simplest [19,26]. Couple stress theory builds on top of classical elasticity (CE) by adding an extra kinematical field, named the micro-rotation. However, in contrast to micro-polar theory, micro-rotation is taken to be entrained by the displacement field [15]. As a result, a length scale is introduced in the system, and the theory is no longer self-similar. Consequently, some unphysical results appearing in CE are repaired, such as the non-dispersive nature of bulk and Rayleigh waves and the lack of support for antiplane (or shear horizontally polarized) Rayleigh and interfacial (Stoneley) waves. It is precisely this last matter that is addressed in this paper.
In recent times, a number of contributions have appeared in the literature concerning the propagation of Stoneley waves, mostly within CE. The speed of Stoneley waves guided by a perfect interface between two elastic half-spaces is determined analytically by Vinh et al. [34], for the case of equal bulk moduli. Consideration of Stoneley waves propagating at a loose interface between two elastic half-spaces is given by Vinh and Giang [33]. Similarly, guided propagation occurring between two half-spaces in elastically constrained contact is considered by Anh et al. [1]. Since CE does not support antiplane Stoneley waves, every listed contribution deals with waves polarized in the sagittal plane, to which the term Stoneley waves is traditionally attached. As notable exceptions, Eremeyev et al. [7] show that antiplane Stoneley waves occur at the interface between two half-spaces in the presence of surface elasticity, while propagation of in-plane Stoneley waves is considered by Kumar et al. [16] within the context of the modified couple stress theory.
To the best of our knowledge, no contribution appears in the literature investigating antiplane Stoneley waves within the couple stress theory. We show that incorporation of microstructural details into the model dramatically alters the propagation features. Indeed, as illustrated in Nobili et al. [22], microstructure provides new pathways for energy transport, which take the form of novel wave propagation patterns being supported. In this paper, we show that, in marked contrast to CE, antiplane Stoneley waves are supported in CS materials under very general conditions concerning the elastic contrast between the media in contact. Furthermore, propagation is only permitted beyond a threshold frequency, named cuton frequency, that is an increasing function of this contrast. Consequently, the role of rotational inertia is especially important for determining the range of admissible parameters for propagation to occur. In this context, a proof of existence and uniqueness of antiplane Stoneley waves is also given.
Traditionally, Stoneley waves have been exploited in borehole seismics to determine the shear-wave velocities at different depths. Also, thin-film applications are possible, as it is demonstrated by Rokhlin et al. [27]. A pioneering application to non-destructive testing (NDT) is experimentally investigated in Hsieh et al. [13]. The recent monograph by Dal Moro [6] collects 14 real-life case study applications of surface and near-surface waves, ranging from geophysics to civil engineering, from geotechnics to moonquakes. Today, the investigation of Stoneley waves is sustained by modern promising developments in the field of acoustic NDT, as recently suggested by Ilyashenko [14]. This paper aims to add a further tool to the investigator, in the form of antiplane interfacial waves reflecting the microstructure underneath the continuous media.

Two couple stress elastic perfectly bonded half-spaces
We consider two half-spaces, named A and B, perfectly bonded along a plane surface. We introduce a righthanded Cartesian coordinate system (O, x 1 , x 2 , x 3 ), whose axes are directed along the relevant unit vectors (e 1 , e 2 , e 3 ). The coordinate system is located in such a way that the plane x 2 0 corresponds to the contact surface between A and B; see Fig. 1. Both half-spaces possess a microstructure, which is described within the theory of linear couple stress (CS) elasticity. Their relevant properties are henceforth denoted by the superscript k ∈ {A, B}.
The stress state in each half-space depends not only on the classical Cauchy force stress tensor s k , but also on the couple stress tensor μ k . The latter characterizes the polar behaviour of the material such that, for any directed surface of unit normal n k , it determines the internal reduced couple vector q k , acting across that surface, where the superscript T denotes the transposed tensor. s k is conveniently decomposed into its symmetric and skew-symmetric part, whereas μ k is decomposed into deviatoric and spherical parts, wherein a dot denotes the scalar product, and 1 is the rank-2 identity tensor.
In each half-space, occupying the volume B k , the internal virtual power may be expressed as follows (see, e.g. Koiter [15] and Ottosen [24]): where grad denotes the gradient operator and a superposed dot denotes time differentiation. Here, u k and ϕ k are the displacement and the micro-rotation vector fields, which are related as follows: where it is understood that a subscript comma denotes partial differentiation, i.e. u i, j (grad u) i j ∂u i /∂ x j , and E is the rank-3 Levi-Civita permutation tensor. Summation over twice repeated subscripts is assumed throughout. We introduce the classical linear strain tensor alongside the torsion-flexure (or wryness) tensor By combining Eqs. (4), (5) and (7), we infer that the torsion-flexure tensor is purely deviatoric, and the couple stress tensor, to all intents and purposes, may be replaced by the sole deviatoric part μ D . To light notation, hereinafter we write μ k with the understanding that μ k D is meant.

Constitutive equations
We assume hyperelastic isotropic material behaviour for both the half-spaces A and B. Accordingly, in each of these we define a free-energy density U k (ε k , χ k ), such that the following constitutive relations hold: where four material parameters are introduced for each space, namely the classical Lamé moduli, Λ k and G k > 0, alongside l k > 0 and −1 < η k < 1, characterizing the microstructure. Positive definiteness of the free energy demands as illustrated by Gourgiotis and Bigoni [12], where necessary conditions for wave propagation are given, with emphasis on antiplane deformations.

Equations of motion
In each half-space, the equations of motion, in the absence of body forces, read div s k ρ kük , where (div s k ) i s k pi, p . Here, ρ k and J k ≥ 0 are the mass density and the rotational inertia per unit volume, respectively. Besides, (axial τ k ) i 1 2 E i j p τ k j p denotes the axial vector attached to the skew-symmetric tensor τ k , such that τ k i j E i j p (axial τ k ) p .

Antiplane shear deformations
We assume antiplane shear deformations, such that, in each half-space, the displacement field u k reduces to the out-of-plane component only , and there is no dependence on the x 3 coordinate. Within this framework and for homogeneous media, following Nobili et al. [22], we get the governing equation in terms of displacement, in which the symbolˆ indicates the two-dimensional Laplace operator in x 1 and x 2 . Equation (10) correctly reduces to the corresponding expressions developed in Fan and Xu [8], Zhang et al. [35], and Zisis [36] in the context of statics and in the absence of rotational inertia. Precisely in this form, Eq. (10) was first given by Clebsch to describe the motion of a plate, including prestress and rotational inertia effects; see Gourgiotis and Bigoni [12] and references therein.

Reduced force traction vector and tangential part of couple stress traction
With Koiter [15], we introduce the reduced force traction vector p k acting across a surface with unit normal n, where μ k nn n k · μ k n k q k · n k and × denotes the cross-product between vectors. Similarly, for the couple stress traction vectorq k , only the tangential part can be really imposed on a boundary, Considering now the boundary surface x 2 0 separating the two half-spaces, we have so that the out-of-plane component of the reduced force traction and the in-plane components of the couple stress traction read, respectively, for medium A, and for medium B.

Nondimensional form of governing equations
We are now in a position to bring in the nondimensional form the governing equations of Sect.
whereupon υβ c s B /c s A . We observe that the limiting case in which the half-space B is constituted by a classical elastic isotropic material, i.e. in the absence of a microstructure for B, can still be retrieved by taking Substituting the thus introduced nondimensional variables in Eq. (10) lends the nondimensional governing equations holding in A and B, respectively. Here, the symbol indicates the two-dimensional Laplace operator with respect to ξ 1 and ξ 2 , whereas the dimensionless parameters k 0 are defined as (see also Mishuris et al. [20] Nobili et al. [21]):

Time-harmonic solution
We consider time-harmonic and straight-crested antiplane waves moving in the sagittal plane (ξ 1 , ξ 2 ): where ı is the imaginary unit and Ω ωT A > 0 indicates the nondimensional time frequency. Substituting the solution form (17) into Eqs. (16), we obtain the pair of meta-biharmonic partial differential equations (PDEs): provided that we make the proper choice for Θ, namely It is easily seen that Θ is a bounded function of Ω, while δ is the wave number of shear horizontal (SH) travelling bulk waves. The latter may be rewritten as Combining Eqs. (20) and (21), we obtain the connection This connection becomes simply A 0 0 A cr in the special case δ δ cr . With such definitions, Eq. (18.2) may be factored as where we have let the dimensionless wave numbers for bulk travelling and bulk evanescent waves in medium B In particular, ψ may be interpreted as a generalization to material B of the parameter δ, Indeed, in the special case of a single homogeneous full-space, that is for β υ 1, we have ψ δ. The following asymptotics hold: We introduce the shorthands whereby we obtain the bulk travelling and bulk evanescent (dimensional) wavespeeds for material B, The corresponding wave numbers δ 1,2 may be used in (25) instead of κ 1,2 when expressing the corresponding bulk wavespeeds in terms of c A s . In the absence of a microstructure for B, that is in the limit (15), we have κ 1 → √ δδ cr , κ 2 /υ → √ δ/δ cr , and we retrieve the classical bulk SH wavespeed For the boundary conditions, Eq. (13.1) takes on the form In the absence of a microstructure for B, Eq. (27) reproduces the classical result

Antiplane Stoneley waves
For guided waves propagating along the interface ξ 2 0, we take wherein we have introduced the shorthand κ Θ K , with K kl A denoting the dimensionless (spatial) wave number in the propagation direction ξ 1 . Moreover, we define the dimensional phase speed in the propagation direction, The general decaying solution of Eq. (19) is and for Eq. (23) where the coefficients a 1,2 and b 1,2 are four amplitudes to be determined. Here, we have let the decay indices for material A, and, similarly, for material B. Here, the square root is made definite by taking the branch that corresponds to the positive square root of any positive real argument. (Equivalently, we may demand that A 1,2 and B 1,2 → √ p, as κ p → +∞ real and positive.) Branch cuts for the square root are taken parallel to the imaginary axis in anti-symmetric fashion; see the discussion in Noble [23, §1.1]. It is emphasized that this is not the choice for the cuts taken in Cagniard [4], where a finite cut is considered instead. However, according to this choice, we get an odd real-valued function along the real axis, and this jeopardizes decay. At any rate, Rayleigh and Stoneley waves need to be slower than the slowest bulk mode.
Clearly, ı A 1,2 and ı B 1,2 are the (dimensionless) wave numbers in the thickness direction ξ 2 , in the relevant material. Consideration of the dimensional wave numbers squared, matches favourably with the corresponding results given by Fan and Xu [8]. As discussed in Nobili et al. [22], for κ > δ, the solution (30) corresponds to a localized travelling wave moving slower than the corresponding bulk wave in material A. Similarly, for κ > δ 1 , the solution (31) corresponds to a localized travelling wave moving slower than the corresponding bulk wave in material B.

Rayleigh function
We define the general form of the Rayleigh function that is valid for either half-space, provided that we substitute λ 1,2 with the relevant decay index along ξ 2 . For instance, for the half-space A, we have λ 1,2 A 1,2 , η η A , and we get R A 0 (κ) R 0 (κ, A 1 , A 2 , η A ). Multiplication by (A 1 − A 2 ) lends the Rayleigh function in the form already exposed in Nobili et al. [21,22], In similar fashion, for the half-space B, we have R B 0 (κ) R 0 (κ, B 1 , B 2 , η B ), and multiplication by ( which is a generalization of Eq. (35), wherein the role of the material parameters is not concealed behind the scaling.
We now prove the existence of at least a real zero for the Rayleigh function. On the real axis, for the travelling bulk wave number κ δ, while we have the asymptotics whereupon the Rayleigh function eventually becomes real negative. We conclude, by continuity, that at least a real root is admitted. This already proves the existence of antiplane Rayleigh waves. In fact, we can prove, by the argument principle, that three pairs of central-symmetric roots are present: one real pair corresponding to Rayleigh waves, one purely imaginary pair, corresponding to Rayleigh-like waves, and a third pair of complex conjugated roots, as discussed in Nobili et al. [22]. This proof is given in "Appendix A".

Frequency equation for antiplane Stoneley waves
For perfect adhesion between the half-spaces at the joining surface ξ 2 0, we enforce the boundary conditions Plugging the solutions (30), (31) into the boundary conditions (39) lends a homogeneous system of linear algebraic equations in the unknown amplitudes a 1,2 , b 1,2 . This system admits non-trivial solutions inasmuch as the following secular (or frequency) equation is satisfied: in which Δ is the determinant of the linear system. Introducing the ratio Γ G B /G A , the determinant in Eq. (40) may be written as with where R A 0 (κ) and R B 0 (κ) are the Rayleigh functions for the relevant half-space and D 1 (κ) is the coupling term. For the latter, we have Here, dependence on Ω by δ, δ 1 , and δ 2 is implicitly assumed. Equation (42) is the CS counterpart of the Rayleigh function given by Cagniard [4,], valid for isotropic CE media. The Rayleigh function (42) exhibits symmetry with respect to A ↔ B inversion, recalling that we also have Γ ↔ Γ −1 and β ↔ β −1 . Moreover, when A B, that is whence propagation is possible only for A 1,2 0, that amounts to finding bulk waves. Similarly, in the absence of either half-space, that is for Γ 0 or Γ → ∞, we find the Rayleigh function of the remaining half-space, that is R A 0 (κ) or R B 0 (κ), respectively.

Antiplane Stoneley waves
We now prove the condition for the existence of a real root for the Stoneley frequency function. We have the asymptotics whence, on the real axis, the frequency equation eventually becomes negative. For a given triple δ, δ 1 and δ 2 , D 0 (κ) is monotonically decreasing, and the possibility of a real root for the frequency equation hinges on the fact that Condition (45) is necessary for the existence of travelling antiplane Stoneley waves. Equality provides a cuton frequency Ω cuton beyond which propagation is possible. In "Appendix B", this simple analysis is put into the wider perspective of determining existence and uniqueness of antiplane Stoneley waves. From it, we see that the following scenarios are possible: -for Ω < Ω cuton , condition (45) is not satisfied: the number of roots is 4, two complex-conjugated, located in the second/fourth quadrant, and two opposite purely imaginary. As discussed in Nobili et al. [22], complex roots represent waves decaying/exploding in every direction and have little significance in unbounded media. Conversely, purely imaginary roots represent Stoneley-like waves travelling in the interior of the medium and decaying/exploding along the interface. Such roots are important in semi-infinite situations. -For Ω > Ω cuton , condition (45) is satisfied: the number of roots is 6, and, alongside the previous four zeros, a pair of real opposite roots, corresponding to travelling Stoneley waves, appears. Figure 2 shows that the cuton frequency is a monotonically increasing function of Γ exhibiting a vertical asymptote. Consequently, a critical value Γ c exists for the ratio Γ , beyond which propagation is blocked. Indeed, for large values of Γ , the root landscape, as it appears from the argument principle, becomes more involved, and, for instance, real (and purely imaginary) roots are eventually lost. A full analysis of all possible scenarios rests outside the scope of this paper. We merely observe that, for a given Ω, the condition (45) demands positivity of a quadratic function of Γ , generally concave upwards, which intersects the Γ -axis to the right of the origin, in the light of (37), at Γ Ω < Γ c ; see Fig. 3. As a consequence, a real interval of admissible shear modulus ratios is highlighted, 0 < Γ < Γ Ω , which accommodates propagation at and beyond the specified frequency Ω. Figure 3 supports the observation that this admissible interval increases with Ω up to the asymptotic value Γ c . In fact, for large Ω, the quadratic positive real root stabilizes very close to Γ c . A different situation develops when rotational inertia in medium A is smaller than in medium B, for example, in our parameter set, we consider A 0 0.3. Then, δ 1 grows with Ω faster than δ and eventually overtakes it. This behaviour reflects itself in that the cuton function exhibits a horizontal asymptote at Ω cuton,asym ≈ 3.20663, corresponding to this overtaking; see Fig. 4. Upon reaching Ω cuton,asym from below, the quadratic function of Γ eventually reverses convexity and moves above the Γ -axis: thus, propagation is admitted for any Γ , as in Fig. 5a. Beyond Ω cuton,asym , the assumption δ > δ 1 is violated: a new admissibility interval may be determined considering the roots Γ 1Ω and Γ 2Ω of D 0 (δ 1 ). However, given that convexity has reversed, propagation now occurs outside the interval [Γ 1Ω , Γ 2Ω ]. We thus see, in contrast to CE, that propagation of antiplane Stoneley waves is largely possible, even when the material properties are significantly different.
Finally, we observe that the case η A 0 is special, for then intersection with the ordinate axis occurs at the origin, that is a double root for Γ , i.e. Γ Ω 0. However, we have already proved that, in this situation, Stoneley waves collapse into bulk waves.

Dispersion curves
Travelling wave solutions are possible inasmuch as a set of real solution pairs (κ S , Ω S ) may be found for the frequency Eq. (41). This is possible in the open interval κ > δ M , whence we retrieved the well-known fact, already pointed out by Cagniard [4], that Stoneley waves are slower than the slowest bulk wave. In our 0.5 , we find that Stoneley waves propagate (a) within the finite admissible range (0, Γ Ω ] for Γ at Ω 1 (solid, black); this interval grows and eventually becomes unbounded upon reaching the asymptotic frequency Ω 3.21 (dashed, red). Beyond this frequency (b), we have δ 1 > δ and admissibility demands D 0 (δ 1 ) ≥ 0, which sets the unbounded admissibility interval (0, Γ 1Ω ] ∪ [Γ 2Ω , ∞) (color figure online) example, δ > δ 1 , whence, in the light of (25), δ 1 is the wave number of the fastest bulk wave while, in its neighbourhood, sits the fastest Rayleigh wave κ 1R . Also, we can show that Stoneley waves are faster than the slowest Rayleigh wave, whose wave number is κ R δ. Indeed, looking at (42), we see that, for κ κ R , the first term drops out (by definition of κ R ) and only negative terms remain, in the light of the Rayleigh function R B (κ) being monotonically decreasing (and zero at κ 1R δ 1 ). A solution κ S for (42), thought of as a function of Γ , is possible only inasmuch as R A (κ S ) > 0, and this occurs only for κ S < κ R . We thus prove the result already observed in Hsieh et al. [13] and in Lim and Musgrave [17]. Indeed, following the latter, "for all geometries examined [in the context of anisotropic CE], the interface wave velocity is found to lie between the higher free surface (generalized Rayleigh) wave velocity and the slowest bulk wave velocity".
In the following, for the sake of definiteness, when plotting dispersion curves we assume the parameters Γ 0.1, υ β 1.1, A 0 ≤ B 0 0.5, and η B 0.5, whereby bulk (and Rayleigh) waves are faster in B, i.e. δ 1 < δ. Figure 6, presents dispersion curves for Stoneley waves expressing the wavespeed c S as a multiple of the shear bulk wave speed of CE, c A s , as in Eq. (29). Similarly to what occurs for Rayleigh waves [21,22], Stoneley waves are initially dispersive, yet they soon develop a horizontal asymptote that is located above (below) the shear wave speed of CE, according to A 0 ≷ A 0 cr . Therefore, the dispersive nature of the wave is restricted to small wave numbers. However, in the absence of rotational inertia in the half-space A (recall Γ is small), that is for A 0 0 (see Fig. 6a)), the curve is monotonically increasing and dispersion is always warranted. As already discussed, bulk non-dispersive waves, moving with the constant speed c A s , are found for either η A 0 or A 0 A 0 cr . So far, the behaviour of Stoneley waves is very similar to that of Rayleigh waves (see Nobili et al. [21, §5]), of which they are perturbations through Γ . However, in contrast to Rayleigh waves, a cuton frequency is met below which propagation is prevented. Therefore, Stoneley waves present a zero-frequency block-band, and their propagation follows an optical branch. Figure 7 illustrates the role of Γ and η A in determining the cuton frequency. Also, propagation is possible below a critical value Γ Ω (case δ > δ 1 ) and outside the finite interval Γ 1Ω < Γ < Γ 2Ω (case δ 1 > δ), as discussed in §6. Particularly, in significant contrast to CE, Stoneley wave propagation occurs in CS elasticity under pretty general conditions, well beyond Wiechert's conditions, and they appear to be the rule rather than the exception.
Still, similarly to Rayleigh waves, Stoneley waves are perturbations of the travelling bulk modes. Indeed, following Nobili et al. [22], a convenient approach to the determination of Stoneley wave numbers is obtained expanding Eq. (40) around δ (or δ 1 , depending which is the largest) setting Solutions of the resulting linear equation in S , are plotted in Fig. 8 against the numerically evaluated frequency spectra. It is seen that approximated and "exact" spectra perfectly overlap in the whole range considered. The coefficients d 2 and a 1 are given in "Appendix C".

Conclusions
Propagation of antiplane Stoneley waves is investigated within the context of couple stress theory, also accounting for rotational inertia. This investigation, which, to the best of our knowledge, has no counterpart in the available literature, attempts to discuss the role of material microstructure in developing new pathways for energy transport at the interface between two media in contact. Indeed, it is shown that antiplane Stoneley The question of existence and uniqueness of such waves is also addressed by the argument principle. It is found that, besides travelling waves, evanescent and decaying/exploding modes are also admitted, in a complex wave pattern. Interestingly, propagation is possible only beyond a threshold frequency, for which an explicit expression is given. In general, this threshold frequency is an increasing function of the elastic contrast between the media in contact. Thus, it appears that lack of propagation for dissimilar materials, typical of classical elasticity, is relaxed into high-frequency propagation by the presence of the microstructure. As a result, rotational inertia plays an important role as it affects the admissibility range for propagation.
We show that Stoneley waves, just like Rayleigh waves, are perturbations of the relevant bulk modes. Thus, antiplane Stoneley waves are perturbations of SH bulk waves. Consequently, an approximated linear (in the wave number) expression for locating Stoneley roots is given, which proves extremely accurate, when compared to the numerical solution, in a wide frequency range. The possibility that antiplane Stoneley waves may propagate under very broad conditions possesses great importance in seismology and in non-destructive testing of materials. Indeed, it provides a simple tool for probing the microstructural properties of the materials in contact.
We determine the number of zeros of (34) in the cut complex plane by the argument principle; see Beardon [3]. To fix ideas, we give the proof for R A (κ), but the argument easily extends to R B (κ). We construct the mapping of the simple curve γ by the Rayleigh function, R A (γ ), and count its index (or winding number). Looking at Fig. 9, we see that γ consists of the circle γ R of arbitrarily large radius R, together with the loops γ ±δ around the centrally symmetric pair of cuts [±δ, ±δ ∓ ı∞) and the loops γ ±ı around [±ı, ±ı∞). By the asymptotics (38), we infer that, when moving along γ R , the image point makes four complete turns around the origin. We now turn to the loops around the cuts, and, in the light of the central symmetry property, only loops sitting in either half-plane are considered, and the resulting winding number is then doubled. On the loop γ −ı , we have A 1 and A 2 purely imaginary, whence Eq. (34) remains in the same form but now in terms of real numbers. In the limit as this loop shrinks down to the cut, γ −ı is mapped onto an open curve approaching the real line from above, i.e. from positive imaginary numbers. Conversely, the loop γ δ is mapped onto an S-shaped open curve  10 Simple curve γ (green, solid) whose mapping by the Stoneley frequency equation D 0 (γ ) is used to determine existence and uniqueness of antiplane Stoneley roots. Here, to fix ideas, we have assumed δ 1 < δ and δ 2 < 1 (color figure online) as in Fig. 11, which intersects the real axis three times, named d 1 < d 2 < d 3 . In particular, d 1 < 0 is located to the left of the origin, while d 2 R A (δ) δ 4 η 2 ≥ 0 is always to the right. Together, D(γ −ı ) ∪ D(γ −ı δ) form a non-simple curve winding once around the origin, that is closed when including the points at infinity. We conclude that six order 1 roots are expected.
Basically, this is a singularly perturbed polynomial equation inasmuch as η is assumed to be small. In this context, we observe that for the case η 0, corresponding to the strain gradient theory, Rayleigh waves collapse into bulk waves as Eq. (34) reduces to λ 1 λ 2 λ 2 2 − λ 2 1 , whose real roots correspond to bulk waves λ 1,2 0. In fact, Rayleigh roots are generally perturbations around either bulk wave speed; see Nobili et al. [22]. It is observed that in Nobili et al. [22] a different choice is made for the cuts, according to which the complex-conjugated pair of zeros may fall outside the Riemann sheet. In fact, with our choice for the cuts, existence of all roots is always warranted.

Appendix B: Proof for the number of roots of the Stoneley frequency function
To this aim, we enlarge our viewpoint and think of D 0 as a function of the complex variable s. Then, D 0 (s) appears centrally symmetric, i.e. D 0 (s) D 0 (−s). We determine the number of zeros of D 0 (s) in the cut complex plane by the argument principle. Accordingly, we determine the index (winding number) of the curve D 0 (γ ), where γ γ R ∪ γ ±δ ∪ γ ±δ 1 ∪ γ ±ı is the simple curve shown in Fig. 10. Here, to fix ideas, we assume δ 1 < δ and δ 2 < 1.
When Γ is small enough, the following analysis resembles that given for the Rayleigh function. By the asymptotics (44), as the point κ moves on the curve γ R , its mapping D 0 (κ) makes four complete turns about the origin, whence the index is 4.
As in Fig. 11, γ δ is mapped into an open loop having three intersections with the real axis, d 1 < 0, d 2 and d 3 > 0, with d 2 D 0 (δ). The explicit expression for d 2 is given in the Appendix. In contrast, d 1 and d 3 may  be found numerically imposing the condition I[D 0 (δ ∓ ε − ı y)] 0, respectively, with ε → 0 + and y > 0.
When Γ is small enough, this loop looks just like the S-shaped curves encountered in the Rayleigh case, but, unlike there, its intersection d 2 with the real axis is not necessarily positive. Indeed, this loop has index − 1 2 inasmuch as d 2 < 0, that occurs for small values of Ω. In this situation, D 0 (s) possesses two pairs of roots: a complex-conjugated pair and a purely imaginary pair. Upon increasing Ω, the cuton frequency Ω cuton is reached such that d 2 0 and the real root κ S is located precisely at the bulk wave number δ. In consideration of the fact that δ is a monotonically increasing function of Ω and so is D 0 (Ω), for Ω > Ω cuton we have that D 0 (γ δ ) winds around the origin as in Fig. 11b. Thus, we find three pairs of roots: a complex-conjugated pair, a purely imaginary pair, and a real pair. Similarly, γ δ 1 is mapped onto a loop closed at infinity which never encircles the origin and contributes nothing to the index. Finally, the loop γ −ı is mapped onto the real axis from above (i.e. from the side of positive imaginary part) moving from left to right; see Fig. 12. This curve brings an index 1 2 regardless of Ω.