M Deschamps

G Huet

Complex surface rays associated with inhomogeneous skimming and Rayleigh waves

Complex surface rays associated with inhomogeneous skimming and Rayleigh waves

A c c e p t e d m a n u s c r i p t 1. Introduction

On the extension beyond the cuspidal edge of the acoustic wave surface of an elastically anisotropic solid is a sharp but non-singular feature in the dynamic response functions.

The arrival times of such features coincide with maxima of the energy flux. These quasi-singular features are shown in many measurements of the point source/point receiver mentioned in numerous publications. This behaviour in the geometric shadow has been predicted by Maris in 1983 [START_REF] Maris | Effect of time phonon wavelenght on phonon focusing[END_REF], by applying to phonons analysis the pioneerly work performed by Pearcey in electromagnetism [START_REF] Pearcey | The structure of an electromagnetic field in the neighbourhood of a caustic[END_REF]. It has been experimentally observed at ultrasonic frequencies by making use of quasi-monochromatic sources [START_REF] Hauser | Internal diffraction of ultrasound in crystal: phonon focusing at long wavelength[END_REF][START_REF] Weaver | Acoustic flux imaging in anisotropic media[END_REF][START_REF] Kim | Ultrasonic imaging of the group velocity surface about the cubic axis in silicon[END_REF]. This phenomenom of diffraction at an edge, which is observed although no actual screen spatially limits the radiation, is called internal diffraction or diffraction at the cusp edge. It cannot be explained by the theory of real rays.

By extending the slowness vector of plane waves to complex values, the propagation of the acoustic energy along a prescribed direction does not correspond to a finite number of arrivals (i.e. the geometrical arrivals) but has a continuum of solutions given by the full set of inhomogeneous plane waves whose energy flux is oriented in the direction of observation. This generalized plane wave is represented by a bivector [START_REF] Hayes | Inhomogeneous plane waves[END_REF][START_REF] Poirée | Les ondes planes évanescentes dans les fluides parfaits et les solides élastique[END_REF][START_REF] Ph | Energy flux for damped inhomogeneous plane waves in viscoelastic fluids[END_REF]. Their propagation has been studied theoretically by Ph. Boulanger and M. Hayes [START_REF] Ph | Bivectors and waves in mechanics and optics[END_REF][START_REF] Ph | Electromagnetic Plane Waves in Anisotropic Media: An Approach Using Bivectors[END_REF] and a general study of their reflection and transmission through the plane interfaces can be found in [START_REF] Deschamps | Reflection and refraction of the inhomogeneous plane wave, Acoustic Interaction with Submerged Elastic Structures, Part I[END_REF][START_REF] Van Den Abeele | Inhomogeneous plane wave scattering and mode stimulation on periodic rough surfaces[END_REF]. This inhomogeneous plane waves obviously do not necessarily refer to wave packet arrivals (in the strict sense of wave propagation they, for instance, do not satisfy Fermat's principle). However specific inhomogeneous plane waves have interesting properties with regards to the acoustic rays since they are close to satisfy the criterion of stationary phase. For these particular complex plane waves, the arrival time corresponding to a complex ray can be associated. This study has been done for a radiating source in an infinite anisotropic medium [START_REF] Deschamps | Complex rays and internal diffraction at the cusp edge. in Surface waves in anisogtropic and laminated bodies and defect detection[END_REF].

These results can immediately be applied to the skimming waves that propagate on the free surface of a semi-infinite anisotropic solid. In addition, a generalisation to the Rayleigh waves is possible which is the purpose of this paper.

In Section 3 we introduce the method of skimming waves and extend the results to the Rayleigh equation. Finally in Section 4, the waves arrival times of all the rays, complex or not, are compared with the Green function giving the dynamic response to a point source located on the free surface.

A c c e p t e d m a n u s c r i p t

Page 3

Problem statement

Let us consider a semi-infinite anisotropic elastic medium, which occupies a region x 3 > 0 , with its surface x 3 = 0 in cartesian coordinates (x 1 , x 2 , x 3 ) . Let us assume that this surface is subjected to the action of a point like source and a point like receiver. This Lamb's problem can be solved by different methods presented in various papers, see for example [START_REF] Bescond | Dynamical surface response of a semi-infinite anisotropic elastic media to an impulsive force[END_REF][START_REF] Every | The elastodynamic response of a semi-infinite anisotropic solid to sudden surface loading[END_REF]. The components of the displacement field in the x 3 -direction is plotted in Fig. 1.

Figure 1 Waveform for the copper crystal at the distance of 5 mm between the source and the observation point when the observation direction is oriented at 60° from the principal direction x 1 .

The normal of the surface coincide with the crystallographic axis [0 0 1].

More details concerning the calculations will be given later on. From the typical waveform chosen, it can be clearly shown the complexity of the waveform due to numerous wavefront arrivals. The associated discontinuities in the temporal shape are in agreement with the arrival times of the wavefront that propagate at the group velocity (see [START_REF] Bescond | Dynamical surface response of a semi-infinite anisotropic elastic media to an impulsive force[END_REF] for this calculation). All types of wave arrivals are observed. These different points are marked as follows: (R) for the Rayleigh wave, (Ps) for the pseudo Rayleigh wave, (i) for the surface skimming waves which have wave vectors on the surface, (o) for the surface skimming waves which have wave vectors out of the surface, (Rc) and (Sc) for the particular contributions associated with the waves that come from the diffraction by the cusps of the Rayleigh and the surface skimming waves, respectively. The folding wavefronts for both the Rayleigh wave and the surface skimming wave are made clear.

In fact, in what follows, we focus our interest only on the wave arrival times which can be observed by the waveform, and in particular on those associated with (Rc) and (Sc) waves. To this end, let us study the inhomogeneous plane wave solutions of both the Christoffel equation and the Rayleigh equation. In the next two chapters, we show how to use such solutions to interpret the arrival time associated with (Rc) and (Sc) times.

The skimming waves

The Christoffel equation

The skimming waves are the plane bulk waves satisfying the wave equation, whose the energy flux is located in the plane (x 1 , x 2 ) . Thus, they are solution of the Christoffel equation without any other conditions. Let us write down this problem with our notations. The wave equation is given by

A c c e p t e d m a n u s c r i p t

Page 4

σ ij, j = ρ ∂ t v i , ( 1 )
where ρ is the mass density and the stress tensor σ ij , for a lossless anisotropic medium is given by the Hooke law

ij ijkl kl C σ = ε , ( 2 )
with the elastic stiffness constants ijkl C and the strain tensor ε kl . Although some numerical results will be presented in Fig. 5 for skimming waves for which the slowness vector is not in the free surface, this case is not treated in this paper for the simplicity.

Let us then restrict the analysis to inhomogeneous plane waves that propagate along the surface. At any point X and time t the acoustic velocity field acquires the form [START_REF] Ph | Bivectors and waves in mechanics and optics[END_REF]:

V = a o P exp iω(t -S / / ⋅ X ) , ( 3 )
where a o is the complex amplitude, P represents the polarisation vector and the variable ω is the angular frequency. The complex vector S / / denotes the slowness vector (the so-called slowness bi-vector). The subsript // indicates that the associated vector belong to the surface. The quantities P and S / / are given by the Christoffel eigenproblem following from the equation of motion

Γ ij v j = ρ δ ij v j , ( 4 )
where Γ represent the Christoffel matrix with components

Γ ij = C iklj s k s l (i, j, k, l = 1...3)
and δ ij is the Kronecker delta. The variables s i denote the slowness components in the x i -axis (i = 1...3) . The condition of vanishing the determinant

| | 0 ij ij Γ -ρδ = ( 5 )
yields the implicit Christoffel characteristic equation.

Our study is restricted to a monoclinic anisotropic semi-infinite medium with a plane of symmetry parallel to the boundary. By virtue of these restrictions, the left hand side of Eq. ( 5) is reduced to a product of two polynomial functions of variables s 1 2 and s 2 2 . One function, corresponding to the plane waves polarized in the (x 1 , x 2 ) plane (noted P and SV), is of the second order. The other, associated with the plane waves polarized orthogonaly to the (x 1 , x 2 ) plane (noted SH), is of the first order.

A c c e p t e d m a n u s c r i p t

The energy equation

The energy equations provide important relations between the complex slowness bivector and the energy flux. In [START_REF] Hayes | Energy flux for trains of inhomogeneous plane waves[END_REF], the following relations between the energy velocity c n , the real part ′ S // and the imaginary part ′′ S // of the slowness bi-vector are obtained

′ S // ⋅ c n = 1 , (6)
′′ S // ⋅ c n = 0 , [START_REF] Bescond | Dynamical surface response of a semi-infinite anisotropic elastic media to an impulsive force[END_REF] where the energy velocity, the time-averaged Poynting vector < ′ J // > and the timeaveraged kinetic and potential energy densities < E c > and < E p > are respectively:

c n = < ′ J // > < E c + E p > , ( 8
) < ′ J // >= - 1 2 Re σ ⋅V ( ) , (9)
< E c >= 1 4 ρV ⋅V , (10)
< E p >= 1 4 ε ij C ijkl ε kl , (11)
here the bar superscript means the complex conjugation and < > implies the averaging by time period. Note that all these new quantities are real. Two other important remarks may be drawn from Eqs. ( 6) and ( 7), which will be used extensively in the next section. First, the damping vector given by ′ S // is always normal to the energy flux direction.

Second, the energy velocity c n equals the phase velocity in the direction of the Poynting vector < ′ J // > .

Energy velocity calculation

In this section, the energy velocity of the skimming waves that propagate along the surface in the fixed direction of the unit vector n is obtained. As it follows from Eq. ( 7) that the component of the slowness vector of the inhomogeneous plane waves is real in the flux direction. Hence, the following relation can be written

S // ⋅ n = ′ s n (12)
where ′

s n = ′ c n -1
is the slowness component in n -direction, which is real in virtue of Eq. [START_REF] Bescond | Dynamical surface response of a semi-infinite anisotropic elastic media to an impulsive force[END_REF]. Noting by θ n the angle between the vectors n and x 1 , the relation ( 12) can be

A c c e p t e d m a n u s c r i p t

Page 6

rewritten as follows

′ s n = s 1 cos Ź θ n + s 2 sinŹ θ n , (13)
where

s 1 = ′ s 1 + i ′′ s 1 and s 2 = ′ s 2 + i ′′
s 2 are the complex slowness component in x 1 -direction and x 2 -direction, respectively. Let us study in detail the restrictive case of wave polarised in the surface. Expressing the variable s 2 from Eq. ( 13) and substituting it in the characteristic polynomial we obtain a new polynomial equation in the form

4 1 0 ( ) 0 k k n k b s s = ′ = ∑ . ( 14
)
The symmetry is broken since this polynomial is of the fourth order of s 1 . The coefficients b k depend on ′ s n . Solving Eq. ( 14) for a fixed value of ′ s n and for a given direction of the observation, i. e. a given value of θ n , is equivalent to find the intersection of the slowness curves and the line, normal to n (see Fig. 2). This parameter has been introduced to calculate the various Green functions using the Cagniard method [START_REF] Van Der Hijden | Propagation of transient elastic waves in stratified anisotropic media[END_REF].

For the illustration, the surface with the normal [0 0 1] of a Copper crystal is selected. The non zero stiffness constants of this cubic crystal are: C 11 = 170 , C 11 = 12.3 and C 44 = 75.5 (GPa). The mass density is 8.9 (g / cm 3 ) . The real part of the slowness vector ′ S // is plotted in Fig. 2. The black and grey lines correspond to the homogeneous and inhomogeneous plane waves, respectively. The dashed line represents the normal to the direction of the observation. As the parameter ′ s n increases from zero (line A) to a large value (line D), different observation can be made graphically. From the line A to the line B two solutions are real. The line B corresponds to the wave arrival since this line is exactly tangent to the slowness curve and the direction n coincide with the direction of the energy flux. At this point two real solutions become identical and beyond this point they become complex conjugated. Other solutions are real until the last arrival time (line C). Finally both solutions are complex. At the arrival time, the additional relation

d ′ s 1 d ′ s n → ∞ (15)

A c c e p t e d m a n u s c r i p t

Page 7 is satisfied. These results are classical and they provide the group velocities of homogeneous plane wave by solving the following system of equations

b k ( ′ s n ) s 1 k k=0 4 ∑ = 0 d ′ s n d ′ s 1 = 0 . ⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪ (16)
Solving this system with respect to ′ s n we obtain the energy velocity for the homogeneous plane waves as follows

′ c n = 1/ ′ s n . (17)
Singularities of the displacements at the wave arrival times can be observed in the Green function [START_REF] Van Der Hijden | Propagation of transient elastic waves in stratified anisotropic media[END_REF]. These arrival times classically defines the wave surfaces by the ray theory (or real rays) [START_REF] Musgrave | Crystal acoustics[END_REF][START_REF] Every | Ballistic phonons and the shape of the ray surface in cubic crystals[END_REF].

For these rays, it is necessarily to satisfy the system [START_REF] Deschamps | Reflection and refraction of the inhomogeneous plane wave, Acoustic Interaction with Submerged Elastic Structures, Part I[END_REF] to have an energy flux along the direction n , but for all complex solutions presented in Fig. 2 the energy flux is along direction in virtue of Eq. ( 7). As shown in [START_REF] Deschamps | Complex rays and internal diffraction at the cusp edge. in Surface waves in anisogtropic and laminated bodies and defect detection[END_REF], for a specific value of ′ s n , among this continuum of solutions, the derivative in Eq. ( 15) exhibits a finite maximum instead of infinite value. In Fig. 3 the function d

′ s 1 / d ′ s n is plotted versus ′ s n for the solutions
shown by the dashed black line and the dashed grey line in Fig. 2. The discontinuity and the maximum are marked by black and grey squares, respectively. The wave arrival times can be found from the following system of equations:

b k ( ′ s n ) s 1 k k=0 4 ∑ = 0 d 2 ′ s 1 d ′ s n 2 = 0 . ⎧ ⎨ ⎪ ⎪ ⎩ ⎪ ⎪ (18)
If the solution exists, then, a concentration of energy is observed for this wave arrival time. In our case, this will define the complex rays since they are associated with specific inhomogeneous plane waves. These inhomogeneous plane waves quasi satisfy the criterion of stationary phase. In other words, they quasi satisfy the Fermat's principle.

A c c e p t e d m a n u s c r i p t

Page 8

Let us now use this result to inspect the skimming waves of the problem described in section 2. To highlight the physical phenomena arising close to the cusps existing on the polar curves of the energy velocity, consider the copper crystal again. Both solutions of Eqs. ( 6) and ( 18) are calculated for various angles of observation and presented by polar plot of ′ c n in Fig. 4. The homogeneous and inhomogeneous plane waves are identified by the black and grey solid lines, respectively. Clearly, the ends of the cusp are extended and connected two by two by the complex solutions. These complex rays generalize the real rays and give us a new explanation to the diffraction phenomena near the cusp denoted by (Sc) in Fig. 1, whose description can be found in the literature [START_REF] Every | Time domain dynamic response functions of elastically anisotropic solids[END_REF]. To make this point clear, the arrival time and the energy velocity associated the black and grey points in Fig. 3 are shown in Figs. 1 and4, respectively. Finally, the important point to note is that the skimming wavefronts start at the ends of the cusp and connect two by two these extreme points. The polar curves of the energy velocity of the surface skimming wave whose wave vector is out of the surface plane are plotted in Fig. 5. For such rays, the theoretical calculations are not presented because the generalisation to complex rays is not performed at the moment and remain to be inspected, only the real wavefront are plotted for the comparison with the Green function, which will be presented in section 5.

The successful use of such complex rays to study the phenomena close and beyond the cusp edge is made clear. Naturally, this study can be generalized for the case of Rayleigh waves.

The dispersion equation

The displacement field associated with the Rayleigh wave is given by

V = P(x 3 )expŹ iω(t -S // . X ) , (19)

A c c e p t e d m a n u s c r i p t

Page 9

where the slowness vector of this surface wave in the plane (x 1 , x 2 ) is denoted by S // = (s 1 , s 2 ) . The function

P(x 3 ) = P k exp(-i ω s 3k x 3 ) k=1 3 ∑ , ( 20
)
consists of three partial mode contributions with the component of the slowness vector s 3k in x 3 -axis. The variables are given by solving the Christoffel equation. Skipping detailed calculations of the Christoffel equation, it can be reduced to the following form [START_REF] Taziev | Dispersion relation for acoustic waves in an anisotropic elastic half-space[END_REF]:

Y 2 (s 3k ) + B 1 Y (s 3k ) + B 0 = 0 , (21)
where

Y (s 3k ) = s 31 s 32 + s 33 s 32 + s 31 s 33 , (22)
is a combination of three unknown slownesses. Two others, s 1 and s 2 , are assumed to be known. In the above expression, the coefficients B 1 and B 0 are polynomial functions of s 1 and s 2 . To complete the description of the surface waves, the boundary condition must be taken into account. In this case the free traction conditions acquires the following implicit form:

Y 2 (s 3k ) + A 1 Y (s 3k ) + A 0 = 0 , (23)
where the coefficients A 1 and A 0 are also polynomial functions of s 1 and s 2 . For more details see [START_REF] Taziev | Dispersion relation for acoustic waves in an anisotropic elastic half-space[END_REF]. The polynomial dispersion equation of the Rayleigh wave is obtained by eliminating the variable Y (s 3k ) from the polynomial functions (in Eqs. ( 21) and ( 23)). This is done by calculating the resultant. After simplification, a polynomial equation in the implicit form is obtained

12 2 2 2 1 0 ( ) 0 k k k s s α = = ∑ , (24)
where the coefficients α 2k (s 2 ) depend on s 2 . The degree of the Rayleigh polynomial depends on the anisotropy and on the surface orientation. In the case under consideration, i. e., when the surface coincides with the plane of symmetry, the degree with respect to s 1 2 is 12, for the more general case of triclinic materials the degree of Rayleigh polynomial is 27 [START_REF] Taziev | Two-partial acoustic surface waves[END_REF][START_REF] Taylor | Surface Waves in Anisotropic Media: the Secular Equation and its Numerical Solutions[END_REF][START_REF] Ting | Explicit secular equations for surface waves in monoclinicwith the simmetry plane at x1=0, x2=0 or x3=0[END_REF]. This degree is reduced to 4, when considering the propagation in a symmetry axis direction [START_REF] Destrade | The explicit secular equation for surface waves in monoclinic elastic crystal[END_REF]. Finally, for the isotropic case this degree is 3.

A c c e p t e d m a n u s c r i p t

Page 10

For a fixed value of the direction of propagation, the real solutions of Eq. ( 24), i.e., for a real couple (s 1 ,s 2 ) , which also satisfy the condition of free surface, correspond to the Rayleigh wave. This real solution is unique, however an infinite number of complex solutions exists for a complex slowness vector on the surface S // . This generalization introduces the Rayleigh inhomogeneous waves in the sense that waves are not only attenuated normal to this interface due to the attenuation of three partial modes but also in the surface in the direction given by the imaginary part ′′ S // of the slowness vector.

Energy relation

As for the skimming waves, the energy equation provides relations between the slowness vector and the Pointing vector of Rayleigh waves. To get these relations, let us multiply Eq. ( 4) by the velocity conjugation and take the time-averaging introduced in section 3.2. This yields the relation

S // ⋅ < J // >= 2 < E c > + i 2ω v j ∂ 3 σ j3 , with j = 1...3 (25)
between the complex Poynting vector < J // >= -

1 2
σ ⋅V on the surface, the time-averaged kinetic density < E c > and the slowness vector S // . The additional term on the right hand side results from the interaction of three partial waves in x 3 -axis. In the similar way, multiplying Eq. ( 2) by the complex conjugation of the strain, the vector J // can be married with the time-averaged potential energy density < E p > ,

S // ⋅ < J // >= 2 < E p > - i 2ω σ j3 ∂ 3 v j , with j = 1...3 (26)
The densities < E c > and < E p > are defined by Eqs. ( 10) and [START_REF] Poirée | Les ondes planes évanescentes dans les fluides parfaits et les solides élastique[END_REF], respectively, while the field is given by Eq. [START_REF] Every | Time domain dynamic response functions of elastically anisotropic solids[END_REF]. Finally, the sum of Eq. ( 25) and the complex conjugation of Eq. ( 26) yield

S // ⋅ < ′ J // >Ź =Ź< E c + E p > -< E ω > -i < E l > , (27)
with

< E ω >= 1 2ω Im(v j ∂ 3 σ j3 ) , (28)
< E l >= -1 2ω ∂ 3 J 3 , (29)

A c c e p t e d m a n u s c r i p t

Page 11

where J 3 = -σ 3 j v j . The relation given by Eq. ( 27) depends on x 3 . To obtain an invariant energy relation, let us integrate this relation along the depth. So, noting

<< (.) >>= < (.) > dx 3 0 +∞ ∫ , ( 30
)
we obtain the following equation:

( ′ S // + i ′′ S // ) ⋅ ′ c n = 1- << E ω >> << E c + E p >> -i << E l >> << E c + E p >> , ( 31
)
with << E l >>= (2ω) -1 (J 3 (x 3 = +∞) -J 3 (x 3 = 0)) which represents the energy leakage out of and in the material and the quantity

′ c n = << ′ J // >> << E c + E p >> (32)
defines the energy velocity of the surface wave. The condition of free surface yields zero normal component of the Poynting vector at the surface, i.e. J 3 (x 3 = 0) = 0 . Moreover, the amplitudes of the physical solutions are exponentially decreasing in the depth and then J 3 (x 3 = +∞) = 0 . In addition, the term << E ω >> can be neglected for small values of inhomogeneities, i.e., under small value of the slowness damping vector ′′ S // , which is our case under consideration. Taking into account all these comments from Eq. ( 27), we obtain

′ S // ⋅ ′ c n → 1, ( 33
) ′′ S // ⋅ ′ c n → 0 . ( 34
)
Figure 6. Representation of the phase, damping and energy velocity vectors of a partial mode of the inhomogeneous Rayleigh wave.

As seen the results (33) and (34) are similar to those for complex skimming inhomogeneous plane waves, cf. Eqs. ( 6) and [START_REF] Bescond | Dynamical surface response of a semi-infinite anisotropic elastic media to an impulsive force[END_REF]. For one partial mode these vector relations are illustrated in Representation of the phase, damping and energy velocity vectors of a partial mode of the in where the surface is Π // and planes ′ Π // and ′′ Π // (orthogonal to n ) contain the phase and the damping vectors, respectively. Of course, three partial modes differ by the bi-vector component in the x 3 -axis only. So, these planes remain identical.

To conclude this part, it is of interest to underline that the notion of Rayleigh wave is generalized due to the presence of attenuation on the surface orthogonal to the direction of the energy flux. Indeed, this last point remains in force for usual Rayleigh wave, for

A c c e p t e d m a n u s c r i p t

Page 12

which ′′ S // = 0 and for which the direction of attenuation coincides with the x 3 -axis. In addition, in this case the term < E ω > is rigorously zero.

Energy velocity calculation

In this section let us use the relations (33) and (34). Since their forms are identical to Eqs. ( 6) and ( 7) for skimming wave, the energy velocity for real and complex Rayleigh waves can be obtained in a similar way. Introduce new slowness variable defined by Eq. ( 13) and express new relation form Eq. ( 24) as follows

C k ( ′ s n )s 1 k k=0 24 ∑ = 0 (35)
Here the coefficients C k depend on ′ s n . For a fixed direction of θ n = 80°, the real parts of the 24 solutions for copper are plotted in the plane ( ′ s 1 , ′ s 2 ) . The pure real solution, corresponding to the usual Rayleigh wave, is shown by black, the complex solutions, close to the so-called Pseudo-Rayleigh waves, are shown by dashed lines. This is a well known special surface wave (see, for example [START_REF] Auld | Acoustic fields and waves in solids[END_REF]). All other solutions are plotted by grey lines. The geometric connection between the line of normal n and the solutions, commented in Fig. 2, is still available for the slowness curves associated with the Rayleigh equation. As an example, let us note that when the line is tangent to the Rayleigh slowness curve, the corresponding flux is oriented in the direction n . Consequently, all what has been developed in section 3. The skimming waves can be applied to the Rayleigh equation. In particular, for large enough values of ′ s n , the complex solutions of the Rayleigh wave equation exist and correspond to the inhomogeneous surface waves. We will come back to their physical meaning further. For small values of ′ s n , the solutions have no physical meaning and they correspond to the fictitious solutions introduced to get the polynomial form of the Rayleigh equation. As illustrated in Fig. 8, the same behaviour of the derivative d ′ s 1 / d ′ s n is observed. This function is plotted for the Rayleigh solution and for positive value of ′ s n versus the phase angle θ (angle between ′ S // and the x 1 -axis). At the arrival time the derivative function is infinite and

d ′ s n d ′ s 1 = 0 . (36)
Consequently, for a fixed direction of the observation, the energy velocity of the homogeneous Rayleigh waves can be obtained by solving the system given by the polynomial Eq. ( 35) and those deduced from Eq. (36). The numerical results are presented in Fig. 9, where the black lines correspond to energy velocity associated with the Rayleigh wave in polar coordinate and the angle is the angle of observation θ n . The black points are those associated with the discontinuities presented in Fig. 8. Of course, near the cusp and beyond it the value d ′

s 1 / d ′
s n is large but finite. This is made clear in Fig. 10, where this behaviour can be observed at the angles of observation θ n = 75° and θ n = 64°. As for the skimming waves, this suggests to interpret these phenomena in terms of complex rays. The combination of the dispersion relation, Eq. ( 35), and the polynomial function given by

d 2 ′ s 1 / d ′ s n 2 = 0 , (37)
provides us the solution in terms of ′ s 1 . Referring to Eq. ( 17), we may immediately calculate the energy velocity ′ c n of these complex rays. As pointed out in Fig. 9, the wavefronts (grey line) associated with these inhomogeneous Rayleigh waves shown in Fig. 6, extend the cusp. The grey square and dot indicate two maxima of the curves presented in Fig. 10. The same calculations can be done for the pseudo-Rayleigh wave represented in Fig. 7. The results are shown in Fig. 11. Dashed grey lines correspond to usual Pseudo-Rayleigh wave.

In order to compare our results with the usual solution, the part of curve corresponding to the solution, which is close to the usual pseudo-Rayleigh wave, is plotted by dashed line. This part corresponds exactly to the energy velocity associated with the solution plotted in dashed black in Fig. 7. When the usual pseudo Rayleigh wave exists, i. e., for some specific angles of observation, the energy velocity of the inhomogeneous Rayleigh wave is very close to that of the pseudo Rayleigh wave, but they are different.

As a matter of fact, the attenuation vector associated with this new solution lies on the surface plane (cf. Fig. 6), while for the usual pseudo Rayleigh wave this vector is in the phase plane ′ Π . It is so because the solutions are sought in terms of complex wave number such that the two vectors ′ S // and ′′ S // are colinear for the usual wave. Besides a different description of the usual Pseudo-Rayleigh wave, when the cusp appears due to the anisotropy, the complex rays, first, take place in the continuity of the small cusp, and second, they still exist even when the classical solution stops. This last point is just noted here without any more explanations.

Finally, it is of great interest to note that two cusps of the Rayleigh wavefront are extended as well as the small cusp of the pseudo-Rayleigh wave.

It is important to repeat that all the plotted complex solutions agree with the energy conditions at infinity which imply Eqs. (33) and (34).

Comparison with the Green function calculation

The results obtained in the previous chapters are compared with those obtained by calculation of the Green function. This Green function is the response to a source δ (t)δ (x 1 )δ (x 2 ) located on the free surface Π // . The solution is explained in terms of 3D integral. By using the Cagniard-de Hoop method, two of them can be calculated and one integral remains to be evaluated numerically. For more detail, see [START_REF] Bescond | Dynamical surface response of a semi-infinite anisotropic elastic media to an impulsive force[END_REF].

The Green function gives the evolution of the vertical displacement of a point M with polar coordinates (r,ϕ) on the surface of the medium. For the distance r , the displacement is calculated for the direction of observation ϕ and time t . For an immediate comparison with the wave fronts, the evolution of this displacement is expressed as function of the speed v = r / t rather than the time t . The discontinuities of the vertical displacement which is observed at the wavefront arrival times correspond now to the energy velocities of each waves. Figure 12-A is a 3D representation of these displacements calculated for the parameter ϕ from 0° to 360°. The plane of the figure corresponds to the surface and the amplitude of the displacement is shown by grey scale; black corresponds to a minimum of amplitude and white to a maximum. First of all, let us examine what is well known. On Fig. 12-B and Fig. 12-C, this field calculation is compared with the real rays of skimming, Rayleigh and pseudo-Rayleigh waves, respectively. In all the figures, these contributions are plotted by blue lines. For the pseudo-Rayleigh wave the blue part corresponds to the usual solution, compare Fig. 11 with Fig. 12-C. By analysing these curves, a perfect synchronisation between the arrivals of wavefronts and the discontonuities of the Green function can be seen. Secondly, let us inspect the complex ray arrival times. On Fig. 12-B, the curvature of the complex wavefronts, which connect the cusp two by two, is visible on the Green function. However, the shape of the front is wide and not precisely identified. Since the complex rays are not the true wavefronts and Eq. ( 15) is not satisfied, they do not corresponds to the propagation of discontinuities. The amplitude of the damping vector ′′ S // on the interface indicates the difference between this complex wavefront and the true wavefront. From this point of view at the edge of the cusp its amplitude is zero and it is very small near the cusp. The wavefront is then well defined. In contrast, the damping vector is, for the present example, larger for ϕ = 45 , as well as the front broadening. For the complex Rayleigh wave, the same good agreement between the extension of the small cusp and the field is observed in Fig. 12-C. Finally, Fig. 12-D reports all the complex and real wavefronts and the field calculation. It is made clear that all the feature of the Green function can be explained by the homogeneous skimming and Rayleigh waves, as is well known, completed by the inhomogeneous skimming and Rayleigh waves introduced in this study.

Conclusion

This paper suggests a method to describe extension of the cuspidal edges of wave surfaces. This wave surface is deduced, without acoustic field calculation, by searching the solution of the Christoffel equation in terms of the inhomogeneous plane waves, which almost satisfy the Fermat's principle. For the homogeneous plane waves, it is necessarily to satisfy the Fermat principle exactly in order to have the Poynting vector oriented in the direction of the observation. But for the inhomogeneous plane wave, this is not necessarily since this condition of energy orientation is taken into account by the fact that the damping vector is chosen normally to observation point.

Figure 2 ., ′ s 2 (

 22 Figure 2. Polar curves ′ s 1 , ′ s 2 ( ) of the real part of the bi-vector slowness ′ S // .

Figure 3 .

 3 Figure 3. Absolute value of the derivative d ′ s 1 / d ′ s n versus the phase slowness ′ s n for θ n = 60°.

Figure 4 .

 4 Figure 4. Polar curves ′ c n ,θ n ( ) of the energy velocity in the interface plane. Black and grey lines are the homogeneous and inhomogeneous in-plane skimming waves, respectively. Black and grey squares correspond to θ n = 60°.

Figure 5 .

 5 Figure 5. Polar curves ′ c n ,θ n ( ) of the energy velocity in the interface plane for the homogeneous out of- plane skimming waves.

Figure 7 ., ′ s 2 )

 72 Figure 7. Polar curves ( ′ s 1 , ′ s 2 ) of the real part of the bi-vector slowness S // of the Rayleigh inhomogeneous wave. Solution of Eq. (35) for θ n = 80°. Black line corresponds to the homogeneous Rayleigh wave, dashed line to the Pseudo-Rayleigh wave and grey line to all other solutions.

Figure 8 .A c c e p t e d m a n u s c r i p t Page 13

 813 Figure 8. Function d ′ s 1 / d ′ s n versus the phase angle between ′ S // and x 1 -axis for the homogeneous Rayleigh wave and θ n = 80°.

Figure 9 .

 9 Figure 9. Polar curves ( ′ c n ,θ n ) of the energy velocity of the Rayleigh waves. Black and grey lines correspond to the homogeneous and inhomogeneous Rayleigh waves, respectively, black squares to θ n = 80° and grey squares to θ n = 75° and θ n = 64°.

Figure 10 .A c c e p t e d m a n u s c r i p t Page 14 Figure 11 .

 101411 Figure 10. Function d ′ s 1 / d ′ s n versus the phase slowness ′ s n for two inhomogeneous Rayleigh wave solutions. On the left θ n = 75° and on the right θ n = 64°.

A c c e p t e d m a n u s c r i p t Page 15 Figure 12 .

 1512 Figure 12. A : Polar curve (r / t, θ n ) of the Green function corresponding to the normal displacement at the distance r=5mm from the source. B: Comparison with the energy velocities of the homogeneous (blue) and inhomogeneous (red) skimming waves. C: Comparison with the energy velocities of the homogeneous (blue) and inhomogeneous (red) Rayleigh waves. D: Comparison with all the wavefronts.

A c c e p t e d m a n u s c r i p t Page 16

 16

A c c e p t e d m a n u s c r i p t Page 19 Fig: 1 A c c e p t e d m a n u s c r i p t Page 20

 19120 Fig:1

Fig: 2 A c c e p t e d m a n u s c r i p t Page 21

 221 Fig:2

Fig: 3 A c c e p t e d m a n u s c r i p t Page 22 Fig: 4 A c c e p t e d m a n u s c r i p t Page 23 Fig: 5 A c c e p t e d m a n u s c r i p t Page 24 Fig: 6 A c c e p t e d m a n u s c r i p t Page 25 Fig: 7 A c c e p t e d m a n u s c r i p t Page 26 Fig: 8 A c c e p t e d m a n u s c r i p t Page 27 Fig: 9 A c c e p t e d m a n u s c r i p t Page 28 Fig: 10 A c c e p t e d m a n u s c r i p t Page 29 Fig: 11 A c c e p t e d m a n u s c r i p t

 322423524625726827928102911 Fig:3

Fig: 12

 12 Fig:12

Acknowledgment : The authors wishes to express gratitude to A. Shuvalov and O. Poncelet for helpful discussions.