COMPLEX RAY IN ANISOTROPIC SOLIDS: EXTENDED FERMAT’S PRINCIPLE

. In contrast to homogeneous plane waves, solutions of the Chris-toﬀel equation for anisotropic media, for which a determined number of rays can be observed in a ﬁxed direction of observation, inhomogeneous plane waves provide a continuum of “rays” that propagate in this direction. From this con- tinuum, some complex plane waves can be extracted for verifying a deﬁnition of quasi-arrivals, based on the condition that the time of ﬂight would vary the less in extension to the Fermat’s principle that stipulates a stationary time of ﬂight for wave arrivals. The dynamic response in some angular zones contain prominent, although not singular, features whose arrivals cannot be described by the classical ray theory. These wave packet’s arrivals can be described by quasi-fronts associated to speciﬁc inhomogeneous plane waves. The extent of the phenomena depends on the degree of anisotropy. For weak anisotropy, such quasi-fronts can be observed. For stronger anisotropy, the use of inhomogeneous plane waves, due to their complex slowness vector, permits a simple description of quasi-arrivals that refer to the internal diﬀraction phenomenon. Some examples are given for diﬀerent wave surfaces, showing how the wave fronts can be extended beyond the cuspidal edges for forming closed wave surfaces.


1.
Introduction. Wave surfaces can be defined as the envelope formed by the normals of the slowness surfaces and vice versa, for example, see [15]. This duality, that takes place in the context of geometrical acoustics, provides the concept of rays referring to the acoustical energy velocity (or group velocity). The extreme point of their energy velocity vector lies on the wave surfaces and this latter one is therefore, by definition, normal to the slowness sheets. Wave surfaces for anisotropic solids may exhibit cuspidal behaviours in certain directions (while for isotropic solids these surfaces reduce to spheres). For such media, if we restrict the consideration to the two-dimensional problem of waves polarised and propagating in a symmetry plane, the number of wavefronts can be either two or four. They are called real rays in this paper. Four solutions appear when direction of observation lies in the folded area of the wave surfaces. If not, only two solutions are present. The junction of the two domains takes place at the cuspidal edges where two solutions of the quasitransverse type merge. Away from this point, the two real solutions disappear. In practice, extending beyond the cuspidal edges there are sharp but non-singular features in the dynamic response functions. The arrival times for these features coincide with the maxima in the acoustic energy flux. These quasi-singular features show up prominently in point source/point receiver measurements, and have been commented on in a number of publications, cf. for exemple [4,7,8,11,12]. This has been experimentally observed at ultrasonic frequencies by making use of quasimonochromatic sources [8,11,17]. 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 is not confined 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. They obviously do not all necessarily refer to wave packet arrivals (in the strict sense of wave propagation they, for instance, do not satisfy the Fermat's principle). However specific inhomogeneous plane waves have very interesting properties with regard to acoustic rays since they are close to satisfying the criterion of minimum time of flight. For these particular complex plane waves, among the continuum of solutions, the time of flight is quasi-stationary. Based on this remark, the equation can be obtained to calculate the "proper" associated energy velocity beyond the cuspidal edges. This equation only refers to complex Christoffel's equation. Limiting the study to planes of symmetry and plotting the associated complex wave surfaces, it can be shown that four rays always exist in any direction for both quasi-isotropic and anisotropic media. Depending on the cusp presence, these rays are all real or two are real or two complex. In other words, it is always possible to define four closed wave surfaces. Depending on the angles of energy propagation, either four surfaces are associated to real rays (cusp area), or the two surfaces are associated with two real rays and two others with two complex rays. These calculations can easily explain the physical phenomena usualy observed in these particular areas.
2. Inhomogeneous plane waves: Background. Let us consider inhomogeneous plane waves. At any space point X and any time t, the acoustic velocity field U is then of the form [9,10,13]: where a 0 is the complex amplitude factor, P represents the complex polarisation vector. The real variable ω is the angular frequency. The complex vector S denotes the slowness vector (the so-called slowness bi-vector). For a lossless anisotropic medium, the stress tensor σ ij is given by: where C ijkl stand for the elastic stiffness constants and ε kl is the strain tensor. The parameters s i and u i (for i = 1 · · · 3 ) denote, respectively, the slowness and the particle velocity components on the X i -axes. The quantities P and S are given by the Christoffel's eigenproblem following from the equation of motion where G represents the Christoffel's matrix. Its components are G ij = C ijkl s k s l -ρδ ij (for i, j, k, l = 1. . . 3 ), where ρ denotes the mass density and δ ij is the Kronecker symbol. The condition of vanishing determinant of G yields the Christoffel's implicit characteristic equation From Eqs. 2 and 3, the following relations between the energy velocity and the slowness bi-vector can be obtained where the new quantity c e is such that: where < J >= 1 2 Re σŪ , < E c >= 1 4 ρUŪ and < E p >= 1 4 ε ij C ijklεkl . The superposed bar means the complex conjugate, Re {x} expresses the real part of the quantity x and <> implies time averaging. The vector < J > represents the time-averaged Poynting vector. The quantities < E c > and < E p > stand for the time-averaged kinetic and potential energy densities, respectively. Two important remarks may be drawn from Eqs. 5 and 6, which will be used extensively in the next section. First, the damping vector given by S is always normal to the energy flow direction given by n e , (see Eq. 6). Second, the energy velocity c e is equal to the phase velocity of the wave in the direction of the Poynting vector. If this direction is given by the unit vector n e , Eq. 5 takes the form: where s e = |c e | -1 is the real part of the slowness component in the n e -direction. A more detailed discussion on energy equations for inhomogeneous plane waves may be found in reference [3,10].
It is important to note that the inhomogeneous planes waves refer, in most cases, to surface waves. However, in this context, they are considered as bulk inhomogeneous plane waves, in a sence that they do not popagate along any surface [10]. Such bulk complex waves can be observed experimentally, in agreement with theory, in certain problems of reflection-refraction [6]. 3. Energy velocity from the Fermat's principle. Let us assume now that the observation point is located at a distance r from the origin and that the vector X is oriented in the n e -direction. The question under consideration in this section is: what are the bulk inhomogenous (or not) plane waves for which the energy velocity is oriented in the n e -direction? For an inhomogeneous wave, the arrival time of the phase front at the point X is given by the relation t = rS · n e . Setting s e = t/r, which is of course real, it is then obtained which is in agreement with the conditions of energy direction Eqs. 5 and 6. In contrast to homogeneous plane waves, i.e. S = 0, the energy velocity direction n e is parallel to the fixed observation direction only for slowness vectors for which this unique direction is normal to the slowness curves, as is well known. This can be obtained, for example, by using the Fermat's principle, which can be written in the general form as follows: where τ (η) represents the wave arrival time and is some vector parameter that affects the propagation. In our problem, the function τ (η) can be replaced by s e and the parameter η can be chosen as the real slowness component s 2 . Then, Fermat's principle gives This condition is satisfied only for the homogenous plane waves for which the normal to the slowness curve is parallel to the n e -direction, as already mentioned. In other words, this condition provides the waves for which the energy direction is given by n e . Singularities, i.e. infinite displacement values at the wave arrival time, are then observed in the Green's function [12]. This classically defines the wave surfaces of the ray theory (or real rays). For inhomogeneous plane waves, Eq. 11 is not satisfied. This is not problematic, since under the condition 9, the energy velocity vector is always oriented along the n e -direction, due to the damping vector has been chosen normal to this direction. In [5], it has been found that, for certain inhomogeneous plane waves, i.e. for a specific value of s e , the derivative function introduced in Eq. 11 can exhibit maxima instead of infinity. Based on this observation, associated wave arrival times can be defined by the following equation: which is solved simultaneously with the Christoffel's equation 3. This equation defines quasi arrival times, in the sense that Fermat's principle is not satisfied. Keeping in mind that the energy direction of an inhomogeneous plane wave is necessarily along the n e -direction, see Eq. 8, this is not inconvenient. The wavefronts associated to Eq. 12 will be called complex wavefronts (or complex rays), in contrast to those, given by solutions of Eq. 11, which will be called real wavefronts (or real rays).
For numerical computations, let us now consider the propagation of two in-plane polarised waves in the principal plane (X 1 , X 2 ) of a model cubic crystal. The stiffness constants are (in GPa): C 11 = 110.6, C 66 = 27.66 and C 12 is obtained from the anisotropy parameter ε = C 11 (2 C 66 + C 12 ) -1 . The mass density is ρ = 2.33 g/cm 3 . In Fig. 1, the absolute value of the derivative function of the real slowness component on the X 2 -axis with respect to the slowness component on the n e -direction is plotted versus the phase velocity along the n e -direction, for three typical cases.
Firstly, both the longitudinal and the transversal modes satisfy Eq. 12 for specific inhomogeneous plane waves, see Figs. 1a and 1b. The anisotropy parameter has been chosen to be ε = 0.935 and the angle of observation is θ = 1°. Clearly, in addition to the velocities of the real waves, which are about 7.0 and 3.6 mm/µs, respectively for the L and T modes at an energy velocity of about 3.3 mm/µs, a sharp peak appears on the derivative function for both modes. Secondly, keeping the same anisotropy parameter but taking the angle of observation of θ = 15°(see Figs. 1c and 1d), only the L mode satisfies this equation and the cusp is not present due to the low anisotropy, see next Figs. 2a and 2b. In such a situation, the single peak is smoothed. Finally, for the anisotropy parameter of ε = 0.761, again only the mode satisfies Eq. 12 but for the chosen angle of observation the cusp exists, see Figs. 1e and 1f, and Fig. 2d in the next section.   They are plotted on two separate graphs (see Figs. 2a and 2b), where the real solutions (L and T ) are, of course, the same on both figures. Note that, even for weak anisotropy, there exist two complex solutions, which almost satisfy Fermat's principle, as previously observed in Figs. 1a and 1b. In addition, it should be remarked that, for the T e mode and for specific angles of observation, the complex solution disappears since the corresponding velocity cannot be greater than the velocity of the homogeneous plane wave. When this mode disappears, i.e. for instance after a small angle of observation, the velocity of the L e mode becomes larger than the T mode and ensures mode continuity. For the two other anisotropies, only the longitudinal mode has a complex solution. When the anisotropy increases (see Fig. 2c), the future cusp influence can be discerned through the L e mode, although the cusp does not really exist. Finally, when the cusp appears (see Fig.  2d), due to higher anisotropy, the T e mode completely vanishes and the L e mode takes place progressively in the continuity of the cusp.

5.
Green's function. Regarding the above results, the question naturally arises: is the new complex solution always observable? To answer this question, let us consider the calculation of the time response of an infinite anisotropic continuum submitted to a line force normal to the n e -direction. The time dependence is described by a delta function. The solution of such a Green problem can be expressed by one Fourier integral and one Laplace integral. Skipping details, which may be found in [12,16], it can be shown, by using the Cagniard-de Hoop method with the time t = rs e , that these two intergrals can be evaluated analytically, and the displacement field d(r, θ, t) then takes the form: where r is the line source-observation point distance, Im {x} indicates the imaginary part of the quantity x, θ stands for the angle of observation (or the angle between the normal to the line source orientation n e and the X 1 -axis) and the variable s 2 represents the slowness on the X 2 -direction given at any time t (or any s e ) by the Christofell's equation 3 such that g(s 1 , s 2 , 0) = 0, with s 1 = (s e -s 2 sin θ) / cos θ. The index n refers to the two solutions of the Christoffel's equation. Solving this equation reveals that the partial wave contribution to the waveform is effective when its associated slowness becomes and stays complex. The terms A n (s e ) are obtained by solving the inhomogeneous problem of the initial differential equation system. They are associated with each partial wave n = L, T and they depend on the type of source. This result is limited for propagation in a principal plane but a more general case can be also decribed by such an approach [16]. Clearly, when Eq. 11 is satisfied, the waveform has a discontinuity which can be associated with a wavefront. It is then not surprising that when ds 2 /ds e exibits a maximum its contribution can also be directly visible on the waveform, as will be further illustrated in Fig. 5. However, this contribution is not necessarily well defined, in terms of time of flight, due to the influence of the source through the functions A n (s e ).
Let us now compare the waveform deduced from Eq. 13 with the extra solutions defined by Eq. 12. To this end, consider a line source along the X 3 -axis that generates a delta force oriented on the X 1 -direction. Fig. 3 shows polar 2D image  that represents the radial particle velocity field calculated from Eq. 13 for various angles of observation θ. The time associated with the waveforms is transformed into a velocity by the relation: c = r/t. The amplitude of the wave, given by Eq. 13, is converted to gray scale, such that zero amplitude corresponds to black pixels. The associated polar energy velocity curves, shown in Fig. 2, are also plotted for comparison. In this figure, it is made clear that, on the one hand, the extension of the cusp is well observed for high anisotropy, as already described in reference [5], and on other hand, even though the cusp does not exist for low anisotropy, a wave packet can be associated to the complex ray noted L e . This can be emphasised by plotting the velocity corresponding to the maximum amplitude in the velocity field extracted from Fig. 3, as illustrated in Fig. 4. Once again very good agreement    between the simple model of the inhomogeneous plane waves and the calculated wave arrival can be observed. For more explanation, let us examine in detail two typical waveforms corresponding to the two anisotropy parameters, as shown on Fig. 5. The graphs represent the total displacement field along the direction, as well as its decomposition into both  i.e., for a energy velocity approximately greater than 3.4 mm/µs. Consequently, the extra longitudinal wave L e cannot be cancelled since this is the only contribution at that time. A broad peak is then visible. Note that the broadening of the peak is directly connected to the imaginary part of the complex solutions. To emphasise this point, let us consider a larger anisotropy, i.e. ε = 0.761 (cf. Fig. 2d) and two angles of observation: θ = 50 • and θ = 25 • , as reported in Figs. 5c and 5d. The first angle is such that the observation point is far from the cusp; the imaginary part is large and the maximum in the waveform is considerably broadened, while for the second angle, the observation point is close to the cusp; the imaginary part is small and a narrow peak is observed. (b) Polar waveforms (radius: r/t; angle: θ). Black: zero amplitude; White: highest amplitude. Figure 6. Wavefronts and 3D polar waveforms for a long fibers composite material.
As a final example, in Fig. 6, let us inspect the case of a long fiber composite material, for which the stiffness constants are (in GPa): C 11 = 188, C 22 = 9.7, C 66 = 6.7 and C 12 = 3.5. The mass density is ρ = 1.7 g/cm 3 . The fibers are oriented along the X 1 -direction. Both ends of the cusp are extended distinctly by the complex solutions. In Fig. 6a, the new solution stops at the real longitudinal wavefront L. In 6b, it continues the circle described by the real shear wavefront T . Note that for this material, paper [14] contains a good agremment between the complex wavefronts and experiments.
In this reference, a comparison between the results presented in this paper and those obtained with the stationary phase method is also discussed. The essence of the stationary phase technique, for instance when approximating the acoustic field in the vicinity of a cusp, is to consider the constructive interferences of the infinite summation of homogeneous plane waves that propagate in different directions. The energy velocity of each wave is, of course, not in the correct energy direction, i.e. along the n e -direction, as well as the wave resulting from the summation. In our case, the maximum of energy in the waveform can be interpreted as constructive interferences between inhomogeneous plane waves, that all propagate in the correct direction, i.e. the damping vector is normal to the n e -direction, but with different amplitude decays. In fact, since the wave vector is complex, more parameters are necessary to describe the wave. As a result, the solutions of Eq. 12 are richer in the domain of wave diffraction.
6. Conclusion. This paper proposes a method to describe extension of the cuspidal edges of wave surfaces in anisotropic media. The extension is deduced, without acoustic field calculation, by searching for the solution to the Christoffel's equation, in terms of inhomogeneous plane waves, that almost satisfies Fermat's principle. For homogeneous plane waves, it is necessary to satisfy Fermat's principle exactly in order to have the Poynting vector oriented in the direction of observation. But for inhomogeneous plane waves, this is not necessary since this condition of energy orientation is taken into account by the fact that the damping vector is chosen normal to observation point.