Geometrical characteristics of phase and group velocity surfaces in anisotropic media

The phase and group velocity surfaces are essential for wave propagation in anisotropic media. These surfaces have certain features that, especially, for shear waves result in complications for modelling and inversion of recorded wavefields. To analyse wave propagation in an anisotropic model, it is important to identify these features in both the phase and group domains. We propose few characteristics for this analysis: the energy flux angle, decomposed in the polar and azimuth angle correction angles and enhancement factor, which is able to characterize both singularity points and triplication zones. The very simple equation that controls the triplications is derived in the phase domain. The proposed characteristics are illustrated for elastic and acoustic anisotropic models of different symmetry classes.


I N T RO D U C T I O N
The wave propagation in anisotropic elastic media is very complex subject of research. In particular, the very common features of S waves in low symmetry anisotropic models are singularity points and triplications of group velocity surfaces. These features make the S wave and converted wave data processing practically impossible in the vicinity of triplication points and triplication cusps. Ray theory byČervený (2001) among others has problems in these points, and the wave amplitudes are not defined. At finite frequencies, these points result in complex waveforms that can be an issue for elastic imaging and Full Wave Inversion (FWI), specifically with linearized gradient (Oh and Alkhalifah, 2016). The details of the connection between the complexity of the group velocity surface and practical applications can be found in Grechka (2017).
In this paper, we introduce a few geometrical characteristics of transformations from the phase velocity surface to the group one and vice versa. One of these characteristics, the enhancement factor, results in a very simple equation defining the triplications for anisotropic media of arbitrary symmetry. The proposed characteristics are evaluated for anisotropy models of different symmetry: transversely isotropic with a vertical symmetry axis, elastic, acoustic and elliptic orthorhombic model, monoclinic model with a horizontal symmetry plane and triclinic model.

P H A S E A N D G RO U P V E L O C I T I E S
The phase velocity (velocity of a plane wave) in anisotropic medium can be obtained from solving the Christoffel equation, where I is the unit matrix, is the Christoffel matrix with the elements (Helbig, 1994) ik = c i jkl n j n l , where c i jkl is the reduced (density normalized) stiffness tensor, and n = (n 1 , n 2 , n 3 ) T is the phase directional vector. The characteristic equation (1) is the cubic equation for phase velocity squared, where coefficients in (3) are invariants of symmetric positive definite matrix , that is, P = tr( ), Q = tr(cof( )) and R = det( ) (Helbig, 1994;Červený, 2001), where tr(), cof(), det() stand for trace, cofactor and determinant of Christoffel matrix. Equation (3) has three real positive roots which are shown in Appendix A. The velocity of the wave packet formed by the superposition of plane waves is called the group (ray) velocity V. In a homogeneous non-dispersive medium (Musgrave, 1970;Cervený, 2001), with the group velocity projections given by where p i is the slowness projection, g i are the corresponding normalized eigenvector projections, and ∇ = ∂ . The equation for the group velocity can be found iň Cervený (2001). The group velocity vector V (from equation (4) is defined by The projections n 1 , n 2 and n 3 are given by n 1 = sin θ cos φ, n 2 = sin θ sin φ, n 3 = cos θ , with θ and φ being tilt and azimuth phase angles, respectively. We can arrive to equations for group velocity projections (Fedorov, 1968;Ben-Menahem and Sena, 1990), From equations (7), we can compute the group velocity magnitude From equations (7), we can also compute the group velocity direction vector, N with projections N i = V i /V , N 1 = sin cos , N 2 = sin sin , N 3 = cos , where and are the tilt and azimuth group angles.
The expression for phase velocity is defined through the slowness notations (the duality between the phase and group domains is exploited), The phase velocity projections can be explicitly given by where From comparison of equations (8) and (11), we have the relation between velocity derivatives in the phase and group domains (12)

G E O M E T R I CA L C H A R AC T E R I S T I C S
The geometrical characteristics of the phase velocity surface v = v(θ, φ) (solution of equation (3)) are defined by the derivatives of velocity upon the polar and azimuth phase angles. By using first-order derivatives, we can compute different geometrical characteristics.
The simplest one defines the angle between the phase and group velocity vectors and called the energy flux (or power flow) angle ψ, where cos ψ = n T N.
From equation (8), it follows that This angle can be interpreted as a measure of local anisotropy in a sense that in case of weak-anisotropy approximation, the energy flux angle is assumed to be zero. In the case of a singularity point, the phase velocity derivatives go to infinity (Appendix B), and the energy flux angle ψ = π/2. Following from equation (12), the energy flux angle can also be defined in a group domain, Other geometrical characteristics defined by the firstorder derivatives are the polar and azimuth angle corrections, θ and φ. They define the difference between the phase and group angles. For the longitudinal direction according to Helbig (1993), both velocity derivatives are zero, and, thus, the energy flux angle being computed for P wave is zero.

Correction angles in phase domain
The polar group angle and the azimuth group angle can be computed from equations (7) as follows: where correction angles θ and φ are defined by The correction angles can be interpreted as the magnitude of anisotropy being decomposed into the polar and azimuth angles.
In vertical symmetry planes of an anisotropic medium with orthorhombic symmetry (ORT), ∂v/∂φ = 0, therefore, φ = 0, equation (17) reduces to the a transversely isotropic medium with a vertical symmetry axis (VTI), and θ = ψ. In the case of singularity point (see Appendix B), the polar correction angle tends to −π/2, since tan θ → −∞. The azimuth correction angle is given by where the phase velocity v is defined by the root of characteristic equation (6), see Appendix A, and indices θ and φ stand for corresponding derivatives (Appendix B).

Correction angles in group domain
The correction angles in the group domain can be introduced in a similar way, From equations (17), (18), (21) and (22), we can compute the difference between group and phase azimuth angles defined in group and phase domains, and the difference between group and phase polar angles defined in group and phase domains,

Enhancement factor
A more complex geometrical characteristic is the enhancement factor (Maris, 1971;Jaeken and Cottenier, 2016) A, which can be introduced as follows: where dω and d are the solid angles spanned by a set of phase velocity vectors n and group velocity vectors N, respectively. Characteristic A from equation (26) is proportional to the ratio of the Gaussian curvatures of the phase and group velocity surfaces at corresponding points and allows to control the abnormal behaviour of these surfaces (Dubrovin et al., 1984). The enhancement factor can also be associated with the ray intensity parameter proposed byČervený (1972), geometrical spreading (Kendall and Thomson, 1989) and wave propagation metric tensor (Klimeš, 2002).
The angle dω has a very simple expression ( Fig. 1): while the angle d is given by where × denotes the vector product. Substituting equations (27) and (28) into equation (26) The enhancement factor defined in equation (29) can be computed in the phase domain by taking the derivatives of   (D4). P wave is shown to the left, S1 wave in the middle and S2 wave to the right. The position of the singularity point is indicated by a black point. the group velocity projections over phase angles (Appendix C). After tedious algebraic manipulations, the final expression for the enhancement factor is given by where the terms R 1 , R 1 and R 1 are defined by the second-order derivatives of the phase velocity: The equation for D controls the abnormal behaviour of group velocity vectors. If D → 0 (A → ∞, d → 0), this point results in a cusp shaping the triplication zone; if D → ∞ (A → 0, d → ∞), this point it is a singularity one. Equation is controlling the triplications in an arbitrary anisotropic medium. For the azimuthally independent medium, ∂v/∂φ = 0, equation (32) reduces to a well-known condition for triplication in a VTI medium (Dellinger, 1991;Thomsen and Dellinger, 2003;Roganov and Stovas, 2010),  In the vicinity of a singularity point, being proportional to the determinant of the matrix of secondorder derivatives of the phase velocity upon phase angles. The enhancement factor in vicinity of the singularity point is given by

Transversely isotropic medium with a vertical symmetry axis
In case of the transversely isotropic medium with a vertical symmetry axis (VTI), ∂v/∂φ = 0, we have = φ = 0, and the energy flux angle (16) is reduced to The enhancement factor (30) is simplified to where the group velocity (8) is defined by In this case, sin θ + 1 v ∂v ∂θ cos θ ≥ 0, and the term 1 v ∂ 2 v ∂θ 2 + 1 controls the SV wave triplications. There are no singularity points since the SH wave is decoupled from equation (3).

Symmetry planes in orthorhombic medium
For vertical symmetry planes in the ORT medium (the firstorder azimuth angle derivative is zero, ∂v/∂φ = 0), the angle equations are the same as for the VTI medium (equation 36), but the second-order derivative over the azimuth angle is not zero. Equation (B6) takes the form Therefore, the term R 3 = 0, and the enhancement factor is different VTI equation (37), The expression in the denominator of (40) is equivalent to the one obtained in Stovas et al. (2019) for SV-SH wave triplications in vertical symmetry planes of the ORT medium.
In this case, the term R 1 = 1 v ∂ 2 v ∂θ 2 + 1 controls in-plane triplication, and the term R 2 = sin 2 θ + 1 v ∂ 2 v ∂φ 2 + 1 v ∂v ∂θ sin θ cos θ controls out-of-plane triplication. The singularity points are achieved when the derivatives of the phase velocity are not defined.
For a horizontal symmetry plane in the ORT medium (the first-order polar angle derivative is zero), the polar angles are the same. θ = 0 and the azimuth angle equation are the same as for VTI case: The second-order derivative over the polar angle, ∂φ 2 and R 3 = 0. The enhancement factor is given by In this case, the term 1 + 1 v ∂ 2 v ∂φ 2 controls in-plane triplication, and the term 1 + 1 v ∂ 2 v ∂θ 2 controls out-of-plane triplication.
When velocity derivatives are not defined, we have the singularity points.

Ellipsoidal case
Explicit equations can be obtained for the ellipsoidal case. In this case, the phase and group velocity equations are given, respectively, by v 2 = v 2 1 sin 2 θ cos 2 φ + v 2 2 sin 2 θ sin 2 φ + v 2 3 cos 2 θ   The correction angles φ, θ and the energy flux ψ are, respectively, defined as The enhancement factor is given by In this case, we do not have any triplications. In Figure 2, we show the geometrical characteristics computed from equations (46), (48) and (49) for (P wave-based) ellipsoidal velocity model with parameters given in equation (D5). One can see that the behaviour of all characteristics is smooth and simple.

Acoustic orthorhombic medium
This model was introduced in Alkhalifah (2003), Stovas (2015) and Masmoudi and Alkhalifah (2016). This special case is an acoustic ORT model (on-axes S wave phase Figure 12 Azimuth correction angle computed for the triclinic model with parameters defined in equation (D2). P wave is shown to the left, S1 wave in the middle and S2 wave to the right. The positions of singularity points are indicated by black points. Note singularity points for P-S1 and S1-S2 waves.

Figure 13
Energy flux angle computed for the triclinic model with parameters defined in equation (D3). P wave is shown to the left, S1 wave in the middle and S2 wave to the right. The positions of singularity points are indicated by black points. Note singularity points for P-S1 and S1-S2 waves. velocities are zero). In this case, the stiffness coefficient matrix takes the form (Abedi et al., 2019) where v j are the on-axis P wave velocities, and η 1 , η 2 and η 3 are the anellipticity parameters defined in the (x,y), (x,z) and (y,z) symmetry planes, respectively. In this case, the P wave is similar to elastic ORT case, while two S wave artefacts are completely different from their counterparts in the elastic case (Xu and Stovas, 2019b).
In any symmetry plane, R = 0, and characteristic equation reduces to which indicates that one of S wave artefact does not exist in the symmetry plane, and another one is coupled with P wave. . P wave is shown to the left, S1 wave in the middle and S2 wave to the right. The positions of singularity points are indicated by black points. Note singularity points for P-S1 and S1-S2 waves. For example, the non-trivial solution of equation (52) in the horizontal symmetry plane is given by where the signs + and -correspond to P wave and S wave artefacts, respectively. For the S wave artefact, in the horizontal symmetry plane, we have R 3 = 0, R 1 ≥ 0 and R 2 < 0. Therefore, the S wave artefact in the symmetry plane is given by the out-of-plane triplication branch. The caustic point is defined by equation R 1 = 0 and is located at an azimuth angle of The caustic condition is given by inequalities The singularity point is located at θ = 0, and there are no other singularity points in symmetry planes since only one S wave artefact exists there. Similar analysis can be performed for vertical symmetry planes.
The azimuth correction angle (equation (19)) computed for the acoustic ORT model (standard) with parameters defined in equation (D4) is shown in Figure 3, and the energy flux angle (equation (15)) is shown in Figure 4. The enhancement factor for acoustic ORT can be seen in Figure 5. The singularity point is clearly seen at phase angles θ = 62.1 • and φ = 16.8 • (the position can be computed analytically, see Xu and Stovas (2019b)) for both S1 and S2 waves. The cusps are identified at φ = 0 and θ = 90 • and φ = 90 • and θ = 90 • for the S1 wave. The S1 and S2 waves in the acoustic orthorhombic model are considered as artefacts. These artefacts are given by triplication sheets.
For the elastic ORT model, the azimuth correction angle, energy flux angle and enhancement factor are shown in Figures 6, 7 and 8, respectively. From the azimuth correction angle and energy flux angle plots (Fig. 6), we can see that the accuracy of the weak-anisotropy approximation is very low for this model. The maximum energy flux angle for P waves reaches 16 degrees, and the one for S1 waves reaches 25 degrees. The singularity points (Table 1) are clearly seen on all characteristics computed for S waves by abrupt increase in density of the contour lines. One can see that results are very different for elastic and acoustic models. The largest difference is observed for S1 and S2 waves (since the S waves in acoustic ORT models are considered as artefacts). For low symmetry models (monoclinic and triclinic), we still plot (for simplicity) only the first quadrant keeping in mind that monoclinic model is not symmetric in azimuth angle, and triclinic model has no symmetry in the polar and azimuth angles. Azimuth correction angle, energy flux angle and enhancement factor computed for the monoclinic model are shown in Figures 9, 10 and 11, respectively. The accuracy of weak-anisotropy approximation is even worse for this model. The maximum energy flux angle is 30 degrees (P wave) and 45 degrees (S waves). Two singularity points (Table 2) are clearly seen in all the S wave plots. The azimuth angle correction factor (Fig. 9) also assists in identifying the triplication zone clearly seen at around θ = 30 • (there is a jump in the azimuth correction factor from −90 to 90 degrees). Azimuth correction angle, energy flux angle and enhancement factor computed for the triclinic model are shown in Figures 12, 13 and 14, respectively. The week-anisotropy approximation is extremely inaccurate (maximum energy flux angles are 45, 60 and 70 degrees for P, S1 and S2 waves, respectively. For this model, we notice the singularity points for P and S1 waves and the singularity points for S1 and S2 waves ( Table 3). The triplication zones are very complex in shape (especially for S2 waves).

C O N C L U S I O N S
We derived geometrical characteristics of the phase and group velocity surfaces in anisotropic media. The energy flux angle is decomposed into the polar and azimuth correction angles. The new characteristic, enhancement factor is derived and used to control different features of the surfaces like singularity points and triplications. This factor provides an easy way to identify the complications in the group velocity surface. From the enhancement factor, a simple equation controlling triplications is derived. This equation is split into three factors responsible for in-plane and out-of-plane triplications in the group velocity surface.

AC K N OW L E D G E M E N T S
Alexey Stovas acknowledges the GAMES project at NTNU. Yuriy Roganov and Vyacheslav Roganov acknowledge Tesseral Technologies Inc. for financial support.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request.

R E F E R E N C E S
In a similar way, the velocity derivative on the azimuth angle can be computed as The second-order derivatives can be obtained from the differentiation of equation (B2): In a similar way, The cross-derivative is given by In a singularity point, all the velocity derivatives given in equations (B3)-(B7) go to infinity. This condition is reduced to where d is defined in equation (A3). Velocity v 2 satisfies the condition v S2 ≤ v 2 ≤ v S1 and can approximately be considered as average S wave velocity. Eliminating phase velocity from equations (B1) and (B8) results in the following expression: which is the well-known condition for multiple roots of the cubic equation.

A P P E N D I X C : G RO U P V E L O C I T Y D E R I VAT I V E S
The group velocity projections are given in equation (7), The derivatives of these projections over the phase angles can be directly computed from equations (C1): After algebraic manipulations, we obtain equation (30)

A P P E N D I X D : A N I S OT RO P I C M O D E L S
In order to illustrate the derived geometrical characteristics, we select the models with different symmetries.
The ORT model (model standard) is defined by the density normalized stiffness matrix (Schoenberg and Helbig, 1997): The triclinic model is adopted from Vavryčuk (2005a): and monoclinic (with a horizontal symmetry plane) is obtained from triclinic one by setting corresponding stiffness coefficients to zero, The acoustic orthorhombic model is obtained from the standard orthorhombic (defined in equation (D1), based on equations (43), The ellipsoidal model (for P wave) can be obtained from the acoustic orthorhombic model given by equation (D4) by setting all anelliptic parameters to zero,