{'fraction_non_alphanumeric': 0.050710177250851037, 'fraction_numerical': 0.022968267010995032, 'mean_word_length': 4.128009630818619, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Supershear surface waves reveal prestress and anisotropy of soft materials

Surface waves play important roles in many fundamental and applied areas from seismic detection to material characterizations. Supershear surface waves with propagation speeds greater than bulk shear waves have recently been reported, but their properties are not well understood. In this Letter, we describe theoretical and experimental results on supershear surface waves in rubbery materials. We find that supershear surface waves can be supported in viscoelastic materials with no restriction on the shear quality factor. Interestingly, the effect of prestress on the speed of the supershear surface wave is opposite to that of the Rayleigh surface wave. Furthermore, anisotropy of material affects the supershear wave much more strongly than the Rayleigh surface wave. We offer heuristic interpretation as well as theoretical verification of our experimental observations. Our work points to the potential applications of supershear waves for characterizing the bulk mechanical properties of soft solid from the free surface.

Surface waves play important roles in many fundamental and applied areas from seismic detection to material characterizations. Supershear surface waves with propagation speeds greater than bulk shear waves have recently been reported, but their properties are not well understood. In this Letter, we describe theoretical and experimental results on supershear surface waves in rubbery materials. We find that supershear surface waves can be supported in viscoelastic materials with no restriction on the shear quality factor. Interestingly, the effect of prestress on the speed of the supershear surface wave is opposite to that of the Rayleigh surface wave. Furthermore, anisotropy of material affects the supershear wave much more strongly than the Rayleigh surface wave. We offer heuristic interpretation as well as theoretical verification of our experimental observations. Our work points to the potential applications of supershear waves for characterizing the bulk mechanical properties of soft solid from the free surface.
Surface wave motion in solids is a classical problem in mechanics, acoustics and seismology, and has found broad applications to nondestructive testing (NDT) of materials [1][2][3], surface acoustic wave devices [4][5][6] and seismic activity monitoring [7][8][9]. The Rayleigh surface wave is the most well studied wave propagating on the free surface of solid. Besides the Rayleigh wave, numerical simulations [10,11] and laboratory experiments [12,13] confirmed the existence of another type of surface wave that propagates faster than bulk shear waves [14]. This so-called supershear Rayleigh wave allows the surface elastic energy to propagate away with a greater speed and is responsible for supershear crack propagation and supershear earthquake [15][16][17]. Recently, growing needs for mechanical characterization of soft matters, such as elastomers, hydrogels and biological tissues, have spurred renewed interests in surface wave motion [18][19][20][21][22][23][24][25]. While the Rayleigh surface wave in soft matters is well established and explored for applications, much less is known about supershear waves in soft materials [21]. Supershear surface waves are inherently leaky in order to manifest a higher speed than bulk shear waves. Interestingly, supershear surface waves were thought to be supported only in highly viscoelastic materials with high wave attenuation, for example, when the shear wave quality factor (i.e., storage modulus/loss modulus) is less than 6.29 [14,26].
In this letter, we reveal several underappreciated properties of supershear waves in soft materials via experimental and theoretical investigations. Firstly, we show the supershear surface waves can be supported over a broad frequency range, with no restriction on dissipation of the material. This invalidates the previous notion that supershear surface waves are only supported in highly dissipative materials. We visualize supershear surface waves in three dimensions in transparent rubber using optical coherence elastography. The measured wave profiles are in excellent agreement with our theoretical model. Secondly, we show that the stress on the soft material alters the speed of the supershear surface wave opposite to what is commonly known for Rayleigh surface waves. We account for this apparently surprising result using the acoustoelastic theory, finite-element (FE) simulation, and simple heuristic explanation. Extending this finding, we show how material anisotropy affects supershear surface waves distinctly different from Rayleigh surface waves.
In Fig. 1(a) we show the schematic of our experiment. The sample used in this study is a silicone rubber (EcoFlexTM 00-50, Smooth-On Inc.). The approximate size is 8×4 cm in the lateral extent and 4 cm along depth, which is large enough to avoid wave reflections at the edges. Mechanical waves were excited at the surface by a vertically vibrating, flat tip with a circular contact area with a radius a of ∼ 0.75 mm, which is driven by a PZT actuator (Thorlabs, PA4CEW). The driving frequency of the PZT, f s , was step-tuned from 2 to 20 kHz. The vertical component of the displacement (on the order of 10 nm) was measured by using a home-built sweep-source optical coherence tomography (OCT) system [23]. The OCT system measures the optical phase delay modulation of a reflected laser beam at different lateral locations (∼ 8 ms per location) with a sampling rate of ∼ 43 kHz. By conducting Fourier transformation of the acquired phase trace, we obtain the amplitude and phase of the displacement (Supplemental Material). To understand the observation, we perform a theoretical study on the near-field displacement of the surface waves. Consider a semi-finite solid which occupies the space z ≥ 0, as shown in Fig. 1(a). Since the problem is axisymmetric, we introduce a cylindrical coordi- Amplitude (nm) Phase (rad) x 1 x 3 x 1 nate system (r, θ, z) and suppose the tip imposes a uniform time-harmonic pressure p(r, t) = p 0 e i2πfst (r ≤ a) to the contact surface, where p 0 denotes the amplitude of the pressure, a the radius of the cylindrical tip, and t the time. The equilibrium equation for isotropic materials is (λ + 2µ)∇∇ · u − µ∇ × ∇ × u = ρ∂ 2 u/∂t 2 , where ρ is the density. λ and µ are Lamé constants, and λ µ for soft materials. Assuming a time-harmonic wave u = u 0 e i2πfst , we get an explicit equation for z component of the surface displacement u 0z (r > a) (Supplementary Material): where R and SS in the sub-and super-scripts indicate the Rayleigh and supershear waves, respectively, J 1 is the Bessel function of the first kind, and H 0 (−kr) → e −ikr r −1/2 ; Equation (1) describes the two propagating waves with wavenumbers, k R and k SS , which are the solutions of a secular equation: where k S = 2πf s ρ/µ is the wavenumber of the bulk shear wave, and F = dF/dk. The two physical solutions are k R = 1.047k S and Here the negative sign in the imaginary coefficient results in an exponential decay of the supershear wave amplitude over r.
The phase velocity c of a wave is related to its wavenumber k, via c = 2πf s / Re(k). For pure elastic materials with a real value of k S , the supershear wave is leaky with a complex value of k SS and has c SS = 2.13c T 0 , where c T 0 = [Re( ρ/µ)] −1 is the phase velocity of the bulk shear wave. For general viscoelastic materials, we find that the ratio c SS /c T 0 increases slightly as the shearwave quality factor Q (= Re(µ)/ Im(µ)) decreases; for example, it is 2.18 at Q = 6 and 2.99 at Q = 0.
The two terms in Eq. (1) results in the interference patterns observed in our experiments. To quantitatively compare the theoretical model and experiments, in Fig.  1(d) we show the experimental normalized displacement data (markers) along the radial direction at f s = 12 kHz and also plot the theoretical curves ordained by fitting the experimental data with Eq. (1) with k S as the only fitting parameter. The best fitting gives k S ≈ 6.93−1.03i mm −1 . Using ρ ≈ 1, 070 kg/m 3 , we obtain the storage modulus µ ≈ 110.8 kPa and the loss modulus µ ≈ 33.7 kPa (µ = µ + iµ ). By performing the Fourier transformation we move the experimental data from the spatial domain to the wavenumber domain (see Fig. 1(e)) and observe two peaks that correspond to the Rayleigh and supershear surface waves, respectively. The theoretical model shows a good agreement with the experiment in the wavenumber domain. Strictly speaking, an inverse Hankel transformation is required to extract the wavenumber from the displacement profiles. However, We find that the Fourier transform is a convenient approximation with errors less than 4% for all our experi-mental cases (see Supplementary Material, Fig. S1). In Figs. 1(f) and (g) we plot the experimentally measured surface displacements u R z (r) and u SS z (r). As shown in Fig. 1(h)), the wave amplitudes decrease exponentially over r, much steeper than the r −1/2 dependence expected for radially propagating surface waves in pure elastic materials [27]. The fast exponential decay is due to the relatively low Q of the rubber sample with µ /µ = 3.3 at 12 kHz.
We further investigated the effect of viscoelasticity by varying the stimulus frequency. Figure 2(a) and (b) show the near-field surface displacements obtained at f s from 2 to 20 kHz in the space and wavenumber domains, respectively (see Supplementary Material, Fig. S2, for the comparison with the theoretical solution). The relative amplitudes of the two modes are slightly frequency dependent. This is because the excitation efficiency of the wave is sensitive to the size, a, of pressure loading. From From the measured wavenumbers, the phase velocities of the two waves are computed. The results are shown in Fig. 2(c) and (d). The velocities were measured to increase with frequency. This indicates the viscoelasticity of the material. The experimental frequency dependence is fit well into the Kelvin-Voigt model using µ = µ e + i2πf s η, where µ e = 101.3 kPa and η = 0.65 Pa·s (see inset in Fig. 2(c)). It has been previously proposed that the supershear surface wave was supported only when Q < 6.29 [26]. However, at f s = 2 kHz we measure k S ≈ 1.25 − 0.052i mm −1 from the experimental data using Eq. (1) and find Q ≈ 12.1. Certainly, Eq.
(1) predicts the existence of the supershear surface wave without restriction on the quality factor in soft materials [14,28].
Typical soft materials including rubber has nonlinear elastic properties or hyperelasticity. Both Young's modulus along the direction of strain and shear modulus facing the direction of the strain increase with the magnitude of strain. In these hyperelastic materials, an extensional prestress increases the velocity of waves propagating along the direction of stretch, and compression decreases the wave speed. This general characteristic has been observed and well documented for the Rayleigth surface wave both experimentally and theoretically. However, little is known about the effect of prestress on the supershear surface wave.
To describe the supershear wave, we use the acoustoelastic theory in continuum mechanics [29,30] and consider a nominal uniaxial stress N 11 in a lateral coordinate, x 1 , and its resulting nominal strain ε along the direction. For isotropic, incompressible, hyperelastic materials, the deformation gradient tensor is given by F = diag(1 + ε, (1 + ε) −1/2 , (1 + ε) −1/2 ). Among several different hyperelastic models, we consider the Mooney-Rivlin model that is widely used for rubbery materials under modest deformation (0 < ε < 100%) [31,32]. It describes the strain energy function as Here, µ is shear modulus of the material at the stress-free condition, ζ is a material-dependent parameter (0 ≤ ζ ≤ 1), I 1 = tr (C), and I 2 = I 2 1 − tr (C 2 ) /2, where C = tr (F T F ) is the right Cauchy-Green deformation tensor. The nominal stress is related to the energy function by N = ∂W/∂F −pF −1 , wherep is the Lagrange multiplier for the incompressibility constraint (or pressure to maintain constant volume). For a small uniaxial stress, The constitutive equation for acoustoelastic waves in the pre-stressed, deformed material can be written with the Eulerian elasticity tensor A 0 defined as where the subscripts i, j, p, q, α, β represent different coordinates. The infinitesimal elastic wave u in the deformed material coordinates is governed by the incremental dynamic equation [33,34] where the Einstein summation convention has been adopted.p denotes the increment ofp. Utilizing the plane wave solution of the dynamic equation and the stress free boundary conditions at the surface (x 3 = 0), we obtain the following secular equation: where c T 0 = µ/ρ is the shear wave speed in the stressfree configuration. Note that the coefficient of the linear ε term for the supershear wave is always negative for ζ = [0, 1]. Using a custom-built mechanical setup and the OCT system ( Fig. 3(a) and Supplementary Material), we measured the velocities of the Rayleigh and supershear surface waves in a rubber sample at different magnitudes of prestress. Figure 3 (b)-(d) shows the experimental results for a range of ε from −0.2 to 0.2. The Rayleigh wave velocity increases with lateral strain as expected. Interestingly, we find that the supershear wave velocity actually decreases with the strain, the opposite behavior. The best fit to Eq. (5) was obtained with ζ ≈ 0 (Supplemental Material), which is shown in Fig. 3. No specific physical interpretation of ζ = 0 in the phenomenological Mooney-Rivlin model is known. The ratio of c SS /c R ≈ 2.23 − 3.30ε + 5.46ε 2 is uniquely related to strain.
To understand the negative sensitivity of c SS to uniaxial extension, we performed FE simulation (Supplemental Material). Figure 3(e) illustrates wave motion in the x 1x 3 cross-sectional plane in the stress-free configuration. Three distinct waves from the excitation point are seen: the Rayleigh wave with phase planes normal to the surface, the supershear wave with phase planes tilt at a angle θ SS , and a spherical shear wave. The supershear wave is a leaky surface wave, whose energy is radiated into the medium in the form of a planar shear wave. This shear wave has a speed of c T 0 . The leaky angle θ SS satisfies Snell's law: c SS cos(θ SS ) = c T 0 , from which θ SS = 62 • . The supershear surface wave can be viewed as a bulk shear wave created at the surface and propagating with the steep angle into the medium. Now, we have a qualitative explanation for the negative dependence of c SS on ε.
As the medium is stretched in x 1 , it is compressed in x 3 by the Poisson's ratio. This compression decreases wave speeds along x 3 , just like compression in x 1 decreases wave speeds along x 1 . Since the propagation direction of the supershear wave is more vertical than horizontal, the effect of prestress on c SS is opposite to that of c R . c T = (α cos 4 θ + 2β sin 2 θ cos 2 θ + γ sin 4 θ)/ρ, (6) and c SS = c T (θ SS ) / cos(θ SS ). Over −0.5 ≤ ε ≤ 0.5, we find the variation of θ SS is < 1.4 • (Supplementary Material). Therefore, c SS ≈ 2.13 c T (θ = 62 • , ε). Whereas the speed of the Rayleigh surface wave c R ≈ 0.955 c T (θ = 0 • , ε).
From the angular interpretation, we expect c SS to be sensitive to material anisotropy. In an orthotropic material with the symmetry orientation aligned with the coordinate system, the shear wave speed is given by Eq. (6) with α = γ = C 55 and 2β = C 11 + C 33 − 2C 13 − 2C 55 , where C ij are components of the stiffness matrix [35] (Supplemental Material). We introduce an anisotropy indexĀ = E * 1 /µ 13 − 4, where E * 1 = (C 11 + C 33 − 2C 13 ) is the plane strain Young's modulus, and µ 13 = C 55 , a shear modulus. Then, c T (θ) = c T 0 [ (cos 4 θ + (Ā + 2) sin 2 θ cos 2 θ + sin 4 θ) ] 1/2 , where c T 0 = µ 13 /ρ (see Fig.  4(a)). For isotropic materials,Ā = 0 and c T (θ) = c T 0 . As shown in Fig. 4, the ratio of the supershear surface wave velocity to the bulk shear speed shows a strong dependence onĀ while the ratio of the Rayleigh surface wave speed to the shear wave speed remains nearly unchanged, especially whenĀ > 0. The leaky angle θ SS has a weak dependence onĀ (Supplementary Material). Our FE simulations agree well with the theoretical curves ( Fig. 4(b)). In conclusion, we have reported on the properties of the supershear surface wave in soft materials. We provided theoretical and experimental descriptions of the effects of viscoelasticity, prestress and anisotropy of the material on the velocity of the supershear wave. The specific properties of the supershear surface wave are quite distinct from those of the Rayleigh surface waves. The difference gives us a new opportunity to use the two wave modes to characterize the material anisotropy or mechanical stress [36] of bulk materials from measurements at their free surface. One possible application is to characterize soft biological tissues using the two surface waves with OCT elastography [37]. The destructive interference of the supershear and Rayleigh surface waves in the near field may be applied to trap and manipulate particle across length scales [38]. Finally, our findings may be useful in the investigations of supershear dynamics of solids, such as ultra-fast dynamic ruptures, where the stress concentration gives rise to the hyperelastic stiffening.

S1. MECHANICAL MODEL FOR SURFACE WAVE EXCITATION WITH A TIME HARMONIC STIMULUS
As shown in Fig. 1(a), let the plane z = 0 be the boundary of the semi-infinite soft solid which occupies the space z ≥ 0. A harmonic pressure p(r, t) = p 0 e iωt is applied on a circular region (r ≤ a) of the free surface where p 0 is the amplitude of the pressure and t denotes the time. ω = 2πf s , where f s is the stimulus frequency. Since the problem is axisymmetric we have u = u 0 e iωt , where u 0 = u 0r e r + u 0z e z . e r and e z denote the unit vectors along r and θ directions, respectively. Inserting u into the equilibrium equation (λ + 2µ)∇∇ · u − µ∇ × ∇ × u = ρ∂ 2 u/∂t 2 , we can get where ρ is the density. λ and µ are Lamé constants and λ µ for soft solids studied here. To solve this problem, we introduce φ = ∇ · u 0 and ψe θ = ∇ × u 0 , where e θ denotes the unit vector along θ direction.

S3
whereσ zz,0 andσ zr,1 are the 0 and 1-order Hankel transformations of σ zz and σ zr , respectively. Solving Eq.(S-7) with the boundary conditions given by Eqs. (S-9) and (S-10), for Rayleigh surface wave we get In Eq. (S-11), we have taken e −z √ k 2 −k 2 S forψ 1 to make sureψ 1 → 0 when z → +∞. For supershear surface wave, we should take e z √ k 2 −k 2 S forψ 1 [2], which results in a plus sign for the second term of F(k). To get a universal form for the secular equation, we replace the equation (S-12) with Substitution of Eq. (S-11) in Hankel transformation of Eq. (S-2), we can get the Hankel transformations of u 0z and u 0rū The displacement u 0z and u 0r can be obtained by performing inverse Hankel transformations to Eq. (S-14).
Here we focus on u z , which is measured in our experiments.
Accordong to Royston et al. [3], an equivalent form of Eq.(S-15) is where K 0 is the modified Bessel function of the second kind. The superiority of Eq. (S-16) is that the Cauchy principal value theorem is applicable. There are some comprehensive discussions on how to deal with the contributions from the poles and the branch cuts [4][5][6] when performing the integration in Eq. (S-16). Since we are only interested in the surface waves, we can neglect the contributions of the branch integrals and only consider the contributions of the residues [4]. Therefore we can get where F = ∂F/∂k, and k denotes the root of F(k) = 0. For soft materials (λ µ) we considered in this study, k, k S k L . In this case Eq. (S-13) reduce to Taking the two roots of Eq.(S-18) that correspond to the Rayleigh and supershear surface waves, i.e., k R = 1.047k 2 and k SS = (0.4696 − 0.1355i)k 2 , we get the Eq. (1) in the main text, where we have used the Hankel function of the first kind H 0 (x) to take the place of K 0 (x) (note that K 0 (x) = iπ/2H (1) 0 (ix)). The other root of of Eq.(S-18) has a positive imagery part, which denotes an exponential increase in wave amplitude along the radial direction. This root thus is excluded.

S2. ACOUSTOELASTIC MODEL
Here the incremental dynamic theory is briefly revisited and then adopted to derive the secular equation for the surface waves that incorporates the effect of prestress. Consider, a nominal stress N homogeneously deforms the sample from the stress-free configuration to the current configuration in which the surface waves propagate. The infinitesimal elastic wave u in the current configuration is governed by the incremental dynamic equation [7] A 0 where the Einstein summation convention has been adopted. A 0 is the Eulerian elasticity tensor and defined as A 0 piqj = F pα F qβ ∂ 2 W/∂F iα ∂F jβ (i, j, p, q, α, β = 1, 2, 3), where W (F ) is the strain energy function, F iα = ∂x i /∂X α is the deformation gradient, and x i and X α denote the Cartesian coordinates of the points in the current and stress-free configurations. The nominal stress is related to the strain energy function by N = ∂W/∂F −pF −1 , wherep is the Lagrange multiplier for the incompressibility constraint.p denotes the increment ofp. Denote the principle stretch ratios by λ i (i = 1, 2, 3), then we have F = diag(λ 1 , λ 2 , λ 3 ). For the incompressible materials we are interested in this study, λ 1 λ 2 λ 3 = 1 and ∇ · u = 0. (S-20) We consider the plane wave in x 1 − x 3 (i.e., u 2 = 0 and ∂()/∂x 2 = 0). From Eq. (S-20) we can introduce a stream function χ(x 1 , x 3 , t) such that: u 1 = ∂χ/∂x 3 and u 3 = −∂χ/∂x 1 . Inserting these into Eq. (S-19) and eliminatingp we can get αχ ,1111 + 2βχ ,1133 + γχ ,3333 = ρ(χ ,11tt + χ ,33tt ), (S-21) where () ,i and () ,t denote the partial deviates with respect to coordinate x i (i = 1, 3) and t, and On the surface (x 3 = 0) the stress-free boundary conditions apply. So the incremental nominal stress, denoted bŷ N , which can be expressed with χ by [7]  For the supershear surface wave, we should take, without loss of generality, s 1 and −s 2 . In this way we get a sign change in Eq. (S-29). Similar as Eq. (S-18) we get Eq. (S-30) is the secular equation that incorporates the effect of the prestress, from which the phase velocity of the surface wave c can be obtained by c = Re(C −1 ) −1 .
While different notations has been adopted, Eq. (S-30) will reduce to Eq. (S-18) in the absence of prestress. In the stress-free configuration α = β = γ = µ, where µ is the linear shear modulus. Recall the definitions C = ω/k and k S = ω/ µ/ρ, we can rewrite Eq. (S-29) as which is equivalent to F(k) = 0 after some simple derivations.

C. Results for Mooney-Rivlin model
To capture the material response, we introduce the incompressible Mooney-Rivlin model where 0 ≤ ζ ≤ 1 is a dimensionless parameter. I 1 = tr (C) and I 2 = I 2 1 − tr (C 2 ) /2, where C = tr (F T F ) is the right Cauchy-Green deformation tensor. The Mooney-Rivlin model reduce to the neo-Hookean model [8] when ζ = 1.
The first-order approximation of c R is c R = (0.955 + 0.615ε) µ/ρ for ζ = 0, c R = (0.955 + 1.09ε) µ/ρ for ζ = 1. (S-37) As shown in Fig. S3, the experimental data approximately follows the theoretical curve given by ζ = 0. We thus used ζ = 0 in Fig. 4 of the main text. The supershear shear surface wave is leaky because of the greater wave speed than the bulk shear wave. According to the Snell's law, the leaky angle θ SS (see Fig. 4 Figure S4 shows the effect of the prestress on the leaky angle. θ SS varies from 62.88 • to 60.64 • over a broad range of nominal strain (−0.5 ≤ ε ≤ 0.5).

S4. FINITE ELEMENT ANALYSIS
The finite element analysis was performed using Abaqus/standard 6.13 (Dassault Systèmes Simulia Corp.). A square domain 40×40 mm 2 was built in the analysis, which was large enough to avoid wave reflection. The left side of the domain was a symmetric boundary. A compression/tension along the horizontal direction was applied in a static analysis step. In the subsequent analysis we used an implicit dynamic analysis step and applied a time-harmonic, local pressure to the surface to excite elastic waves. We adopted a gradient mesh (∼ 0.025 mm to ∼ 0.5 mm from top to bottom of the model) and the CPE8RH element type (8-node biquadratic, reduced integration, hybrid with linear pressure).

S9
S5. EXPERIMENTS Figure S6 shows the experimental setup. The sample was clamped and then compressed/stretched manually using a linear translation stage. We drew a mesh grid on the surface of the sample and measured the stretch ratio by the deformation of the mesh grid. The sample was prepared by mixing the two parts with at a 50:50 mass ratio (Ecoflex 0030, Smooth-On Inc.), pouring the mixture into a mold, and then curing the sample at room temperature overnight. Figure S7 shows the schematic of the M-B scan. The laser beam scanned synchronously with the stimulus signal sent to the PZT. At each lateral location, we acquired ∼ 350 A lines (M scan) at a sampling rate of ∼ 43 kHz. Then the laser beam moved to the next localization (B scan). From each M scan, the vibration on the surface was Fourier transformed to obtain the amplitude A and phase ϕ, i.e., u z (t) = Ae i(ωt+ϕ) . Then we reported the real and imagery parts of the surface displacement.