# High-order perturbations of the phase velocity and polarization of qP and qS waves in anisotropic media

Véronique Farra

► To cite this version:

Véronique Farra. High-order perturbations of the phase velocity and polarization of qP and qS waves in anisotropic media. Geophysical Journal International, 2001, 147, pp.93-104. 10.1046/j.0956-540X.2001.00510.x. insu-03597750

HAL Id: insu-03597750

https://insu.hal.science/insu-03597750

Submitted on 4 Mar 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# High-order perturbations of the phase velocity and polarization of $qP$ and $qS$ waves in anisotropic media

Véronique Farra

Département de Sismologie, Institut de Physique du Globe de Paris, 4 Place Jussieu, 75252 Paris Cedex 05, France. E-mail: farra@ipgp.jussieu.fr

Accepted 2001 May 15. Received 2001 May 15; in original form 2000 October 16

## SUMMARY

An approximate expression of the eikonal equation and the polarization vector may be obtained in weakly anisotropic media from first-order perturbation theory. The advantage of this approximation for $qP$ wave is that the squared phase velocity is linear in the elastic parameters. For $qS$ waves, the first-order approximation is more complicated and can be expressed in terms of the square root of a quadratic function in the elastic parameters. Higher order perturbations can be obtained by an iterative procedure which improves the accuracy of the approximations. Explicit analytic formulae of the approximate squared phase velocities and polarizations are given for orthorhombic and transversely isotropic symmetries. Numerical comparisons between the exact and the approximate phase velocities and polarization vectors obtained at different orders illustrate the accuracy of the approximate formulae presented. For realistic anisotropy, the second-order expressions of the squared phase velocities are accurate approximations which do not cost much more with respect to the first-order computations. Third order expressions of the squared phase velocities are very accurate and need only computation of the first-order approximations. Second order expressions should be used to have good approximations of polarization vectors outside the vicinity of singularities. Higher order approximations of the $qS$-waves eigenvectors should be applied in neighbouring directions of singularities.

**Key words:** phase velocity, seismic anisotropy, perturbation method, polarization vector.

## 1 INTRODUCTION

Three elastic waves (the $qP$ and the two $qS$ waves) can propagate along any direction in an unbounded anisotropic medium. The slowness surface represents the directional dependence of the inverse phase velocity of the three waves. It consists of three sheets; the inner sheet associated with the $qP$ wave is convex, but the other two sheets, associated with the $qS$-waves, can display concave and saddle-shaped regions in addition to convex regions. There are special directions (singularities) for which the two quasi-shear waves sheets come into contact. The shape of the slowness surface plays a central role in the interpretation of a wide range of wave phenomena in anisotropic solids.

In an anisotropic material, not only the phase velocities depend on the propagation direction but, in general, the polarization directions are neither parallel nor perpendicular to the propagation direction. The phase velocities and the polarization directions correspond to the eigenvalues and eigenvectors of the so-called Christoffel matrix. Analytic expressions of the eigenvalues and the eigenvectors can be found only for simple symmetries (isotropy and hexagonal symmetry); in most symmetries,

the eigenvalues and the eigenvectors have to be calculated by numerical algorithms.

Based on the observation that most anisotropic media are weakly anisotropic, several researchers have proposed approximations for wave speeds. Approximate equations for phase velocities in weakly anisotropic media were first derived by Backus (1965). Thomsen (1986) derived expressions of the phase velocities in the case of transverse isotropy. Sayers (1994), Mensch & Rasolofosaon (1997) and Pšenčík & Gajewski (1998) obtained approximate relations for the phase velocity of the $qP$ wave in arbitrary symmetry. As shown by Backus (1965) the approximate formulae for the phase velocity in weakly anisotropic media follow simply from the first-order perturbation theory for anisotropic media (see also Červený & Jech (1982) and Jech & Pšenčík (1989)). The advantage of this approximation for $qP$ waves is that it is linear in the elastic parameters. For $qS$ waves, the first-order approximation is more complicated and can be expressed in terms of the square root of a quadratic function in the elastic parameters. Second order expressions have been obtained by Farra (1999) in the case of transverse isotropy in order to improve the accuracy of traveltime computation.
The approximate formulae obtained for the qP wave have found important applications in the approximate evaluation of kinematic and dynamic quantities such as rays, traveltimes, polarization vectors and amplitudes which are used in the computation of the qP-wave Green's function by the ray method in inhomogeneous weakly anisotropic media (Červený & Jech 1982; Farra 1989; Nowack & Pšenčík 1991; Pšenčík & Gajewski 1998; Mensch & Farra 1999). Higher order expressions may be useful to improve the accuracy of such computations (Druzhinin 1996; Farra 1999).

Shear-wave modelling is complicated by the presence of the shear-wave singularities. Shear-wave singularities cause anomalies in the polarization of the wavefield and in the geometry of wave surfaces (Crampin & Yedlin 1981; Helbig 1994; Rúmpker & Thomson 1994; Vavryčuk 1999). They can cause breakdown of modelling algorithms connected to numerical instabilities that arise whenever the velocity sheets of two waves are close to each other. Analytic approximative formulae of the phase velocities may be useful to provide explicit solutions to various problems of wave propagation in weakly anisotropic media or in the neighbourhood of singularities.

In this paper, an iterative procedure is used to obtain explicit analytic formulae of the Christoffel matrix eigenvalues and eigenvectors corresponding to higher and higher orders of perturbation. The eigenvalues and the eigenvectors of the Christoffel matrix give the squared phase velocities and the polarization vectors of the waves propagating in the corresponding homogeneous anisotropic medium. For a realistic anisotropy, numerical comparisons between the exact and the approximate quantities (phase velocities and polarization vectors) obtained at different orders illustrate the accuracy of the approximate formulae.

## 2 CHRISTOFFEL EQUATION

Let us introduce the Christoffel matrix $\Gamma(\mathbf{p})$, whose elements are dependent of a vector $\mathbf{p}$ and given by:

$$ \Gamma_{jk} = p_i p_l a_{ijkl}, \quad (1) $$

with implicit summation on repeated indices. The parameters $a_{ijkl} = c_{ijkl}/\rho$ are the density normalized elastic parameters and $p_i$ are the components of the vector $\mathbf{p}$. The matrix $\Gamma(\mathbf{p})$ is a symmetric matrix with three positive eigenvalues $G_m(\mathbf{p})$. The corresponding eigenvectors $\mathbf{g}_m(\mathbf{p})$ (defined as unit vectors) are mutually orthogonal and satisfy the Christoffel equations:

$$ (\Gamma - G_m \mathbf{I}) \mathbf{g}_m = 0, \quad (2) $$

where $\mathbf{I}$ is the $3 \times 3$ identity matrix.

Three waves (the quasi-P wave and two quasi-S waves denoted by $qP$, $qS_1$ and $qS_2$, respectively) can propagate in the anisotropic solid defined by the elastic parameters $a_{ijkl}$. Each wave is associated with one of the eigenvalues $G_m$ denoted by $G_P$, $G_{S_1}$ and $G_{S_2}$, respectively, and defined by ascending order $G_{S_2}(\mathbf{p}) \le G_{S_1}(\mathbf{p}) < G_P(\mathbf{p})$ for every $\mathbf{p}$. The corresponding eikonal equation may be written in the following form (Červený 1972):

$$ G_m(\mathbf{p}_m) = 1, \quad m = P, S_1, S_2, \quad (3) $$

where the vector $\mathbf{p}_m = \mathbf{n}/V_m(\mathbf{n})$ is the slowness vector of the considered wave for the phase propagation direction defined by

the unit vector $\mathbf{n}$. The phase velocity squared $V_m^2$ is given by

$$ V_m^2(\mathbf{n}) = G_m(\mathbf{n}). \quad (4) $$

The polarization vector of the wave is parallel to the corresponding eigenvector $\mathbf{g}_m$.

In an isotropic medium, the eigenvalues of the two shear-waves are equal and the polarization vector of the shear wave may take an arbitrary orientation in the plane orthogonal to the P-wave polarization vector $\mathbf{g}_P$.

## 3 FIRST-ORDER PERTURBATIONS FOR THE qP AND qS WAVES

Analytic expressions of the Christoffel matrix eigenvalues can be found only for simple symmetries (isotropy or hexagonal symmetry). Fortunately, geological media are often weakly anisotropic which makes perturbation techniques relevant. In this section, first-order perturbation theory is used to solve approximately the Christoffel equation. Eigenvalues and eigenvectors are expanded into perturbation series due to the perturbation of the elastic parameters with respect to a reference isotropic medium. The procedure is well-known in many domains of physics, such as quantum mechanics (see for example Morse & Feshbach 1953; Landau & Lifshitz 1966). For the seismological applications of first-order perturbation approaches, Backus (1965) Červený (1982) and Hanyga (1982) suggested linear formulae for the determination of phase velocities of the body waves propagating in anisotropic media. For a reference isotropic medium, these formulae can be used for qP waves only. To obtain formulae applicable to qS waves as well, degenerate perturbation theory must be used (see Landau & Lifshitz 1966; Jech & Pšenčík 1989).

Let us assume that we have a reference medium characterized by elastic parameters $a_{ijkl}^{(0)}$. We denote $\Gamma^{(0)}(\mathbf{p})$ the reference Christoffel matrix, $G_m^{(0)}(\mathbf{p})$ and $\mathbf{g}_m^{(0)}(\mathbf{p})$ the reference eigenvalues and corresponding eigenvectors. Let us consider a perturbation of the model, such that the elastic parameters are changed from $a_{ijkl}^{(0)}$ to $a_{ijkl}^{(0)} + \Delta a_{ijkl}$. In the perturbed model, the Christoffel matrix $\Gamma(\mathbf{p})$ is given by $\Gamma = \Gamma^{(0)} + \Delta\Gamma$, where the elements of $\Delta\Gamma(\mathbf{p})$ are given by $\Delta\Gamma_{jk} = p_i p_l \Delta a_{ijkl}$. The expansion of the Christoffel equation (2) to first-order gives:

$$ (\Delta\Gamma - \Delta G_m \mathbf{I}) \mathbf{g}_m^{(0)} + (\Gamma^{(0)} - G_m^{(0)} \mathbf{I}) \Delta\mathbf{g}_m = 0, \quad m = P, S_1, S_2, \quad (5) $$

where $\Delta G_m(\mathbf{p})$ is the perturbation of the considered eigenvalue and $\Delta \mathbf{g}_m(\mathbf{p})$ is the perturbation of the associated eigenvector. Moreover, from the requirement that $\mathbf{g}_m$ is a unit vector, $\mathbf{g}_m \cdot \mathbf{g}_m = 1$, one gets to first order:

$$ \mathbf{g}_m^{(0)} \cdot \Delta \mathbf{g}_m = 0. \quad (6) $$

Let us assume that the reference medium is isotropic. Let us introduce the reference vector system ($\mathbf{g}_P^{(0)}, \mathbf{g}_1^{(0)}, \mathbf{g}_2^{(0)}$), where $\mathbf{g}_P^{(0)}(\mathbf{p})$ is the eigenvector corresponding to the P-wave eigenvalue $G_P^{(0)}(\mathbf{p})$ in the isotropic reference medium, $\mathbf{g}_1^{(0)}(\mathbf{p})$ and $\mathbf{g}_2^{(0)}(\mathbf{p})$ are two mutually perpendicular unit vectors situated in the plane orthogonal to $\mathbf{g}_P^{(0)}$. The freedom in the choice of vectors $\mathbf{g}_1^{(0)}$ and $\mathbf{g}_2^{(0)}$ is a consequence of the coincidence of the eigenvalues of the two quasi-shear waves in the reference isotropic medium.
In the isotropic reference medium, we can write the eigenvalues of the Christoffel matrix as:

$$G_P^{(0)}(\mathbf{p}) = \mathbf{g}_P^{(0)\prime} \Gamma^{(0)}(\mathbf{g}_P^{(0)}) = A_{33}^{(0)} \mathbf{p}^2 \quad (7)$$

$$G_S^{(0)}(\mathbf{p}) = G_{S_1}^{(0)}(\mathbf{p}) = G_{S_2}^{(0)}(\mathbf{p}) = \mathbf{g}_1^{(0)\prime} \Gamma^{(0)}(\mathbf{g}_1^{(0)}) = \mathbf{g}_2^{(0)\prime} \Gamma^{(0)}(\mathbf{g}_2^{(0)}) = A_{44}^{(0)} \mathbf{p}^2 \quad (8)$$

where we use the classical Voigt notation of contracted indices for the components of the fourth-rank tensor $a_{ijkl}$ (Auld 1973) and the superscript t to denote the transposed vector. $A_{33}^{(0)}$ and $A_{44}^{(0)}$ correspond to the P and S-waves velocities squared of the isotropic medium.

The eigenvector $\mathbf{g}_P^{(0)}(\mathbf{p})$ corresponding to the eigenvalue $G_P^{(0)}(\mathbf{p})$ is:

$$\mathbf{g}_P^{(0)} = \frac{1}{p} \begin{bmatrix} p_x \\ p_y \\ p_z \end{bmatrix}, \quad (9)$$

where $p_x, p_y$ and $p_z$ are the components of the vector $\mathbf{p}$ and $p = \sqrt{p_x^2 + p_y^2 + p_z^2}$.

The unit vectors $\mathbf{g}_1^{(0)}$ and $\mathbf{g}_2^{(0)}$ may be chosen as

$$\mathbf{g}_1^{(0)} = \frac{1}{pp_r} \begin{bmatrix} -p_x p_z \\ -p_y p_z \\ p_x^2 + p_y^2 \end{bmatrix}, \quad \mathbf{g}_2^{(0)} = \frac{1}{p_r} \begin{bmatrix} -p_y \\ p_x \\ 0 \end{bmatrix}, \quad (10)$$

with $p_r = \sqrt{p_x^2 + p_y^2}$.

Moreover, let us note the following properties:

$$\mathbf{g}_1^{(0)\prime} \Gamma^{(0)}(\mathbf{g}_P^{(0)}) = \mathbf{g}_2^{(0)\prime} \Gamma^{(0)}(\mathbf{g}_P^{(0)}) = \mathbf{g}_1^{(0)\prime} \Gamma^{(0)}(\mathbf{g}_2^{(0)}) = 0. \quad (11)$$

Let us specify the expression (5) for the quasi-P wave. The scalar product of eq. (5) written for $m=P$, with $\mathbf{g}_P^{(0)}$ gives:

$$\Delta G_P = \mathbf{g}_P^{(0)\prime} \Delta \Gamma \mathbf{g}_P^{(0)}. \quad (12)$$

The perturbation vector $\Delta \mathbf{g}_P$ can be obtained from the scalar products of eq. (5), written for $m=P$, with $\mathbf{g}_n^{(0)}$, $n=1, 2$, taking into account eq. (6):

$$\Delta \mathbf{g}_P = \sum_{n=1,2} \frac{\mathbf{g}_n^{(0)\prime} \Delta \Gamma \mathbf{g}_P^{(0)}}{G_P^{(0)} - G_S^{(0)}} \mathbf{g}_n^{(0)}. \quad (13)$$

The first-order expression of the eigenvalue $G_P = G_P^{(0)} + \Delta G_P$, denoted by $G_P^{(1)}$, can be expressed as, see eqs (7) and (12):

$$G_P^{(1)} = \mathbf{g}_P^{(0)\prime} \Gamma \mathbf{g}_P^{(0)} \quad (14)$$

and the first-order expression of the eigenvector $\mathbf{g}_P = \mathbf{g}_P^{(0)} + \Delta \mathbf{g}_P$ is given by:

$$\mathbf{g}_P^{(1)} = \mathbf{g}_P^{(0)} + \sum_{n=1,2} \frac{\mathbf{g}_n^{(0)\prime} \Gamma \mathbf{g}_P^{(0)}}{G_P^{(0)} - G_S^{(0)}}, \quad (15)$$

where we use the property (11), $\mathbf{g}_P^{(0)\prime} \Gamma^{(0)}(\mathbf{g}_n^{(0)}) = 0$, $n=1, 2$.

Let us specify the expression (5) for the quasi-shear waves. The $qS_1$ and $qS_2$ zero-order eigenvectors $\mathbf{g}_{S_1}^{(0)}$, $k=1, 2$, are situated in the plane orthogonal to the zero-order P-wave polarization vector $\mathbf{g}_P^{(0)}$. Though the choice of the vectors $\mathbf{g}_1^{(0)}$ and $\mathbf{g}_2^{(0)}$ is not unique, the eigenvectors $\mathbf{g}_{S_1}^{(0)}$ and $\mathbf{g}_{S_2}^{(0)}$ can not be

chosen arbitrarily because they are subject to the requirement that the change in them, denoted by $\Delta \mathbf{g}_{S_1}$ and $\Delta \mathbf{g}_{S_2}$, should be small under the action of the perturbation $\Delta a_{ijkl}$. One can write the zero-order eigenvectors $\mathbf{g}_{S_1}^{(0)}$ as linear combinations of the vectors $\mathbf{g}_1^{(0)}$ and $\mathbf{g}_2^{(0)}$:

$$\mathbf{g}_{S_k}^{(0)} = \sum_{n=1,2} \alpha_{S_k}^n \mathbf{g}_n^{(0)}. \quad (16)$$

The scalar products of (5), written for $m=S_k$, with $\mathbf{g}_n^{(0)}$, $n=1, 2$, give:

$$\begin{align}
(B_{11} - \Delta G_{S_k}) \alpha_{S_k}^1 + B_{12} \alpha_{S_k}^2 &= 0 \nonumber \\
B_{12} \alpha_{S_k}^1 + (B_{22} - \Delta G_{S_k}) \alpha_{S_k}^2 &= 0 \tag{17}
\end{align}$$

where

$$B_{ij} = \mathbf{g}_i^{(0)\prime} \Delta \Gamma \mathbf{g}_j^{(0)}. \quad (18)$$

The condition of solvability of system (17) gives (see Jech & Pšenčík 1989):

$$\Delta G_{S_k} = \frac{1}{2} \left[ B_{11} + B_{22} \pm \sqrt{(B_{11} - B_{22})^2 + 4B_{12}^2} \right]. \quad (19)$$

The first-order expression of the eigenvalue $G_{S_k} = G_{S_k}^{(0)} + \Delta G_{S_k}$, denoted by $G_{S_k}^{(1)}$, can be expressed as, see (8), (11) and (19):

$$G_{S_k}^{(1)} = \frac{1}{2} \left[ M_{11} + M_{22} \pm \sqrt{(M_{11} - M_{22})^2 + 4M_{12}^2} \right], \quad (20)$$

where

$$M_{ij} = \mathbf{g}_i^{(0)\prime} \Gamma \mathbf{g}_j^{(0)}. \quad (21)$$

The eigenvalue of the $qS_1$-wave is given by eq. (20) with a positive sign in front of the square root. The components of the corresponding vector $\mathbf{g}_{S_1}^{(0)}$ can be deduced from eq. (17) by taking into account that $\mathbf{g}_{S_1}^{(0)}$ is a unit vector (Pšenčík 1998):

$$\begin{align}
\alpha_{S_1}^1 &= \sqrt{\frac{1}{2} \left( 1 + \frac{M_{11} - M_{22}}{\sqrt{\Delta}} \right)} \nonumber \\
\alpha_{S_1}^2 &= \operatorname{sign}(M_{12}) \sqrt{\frac{1}{2} \left( 1 - \frac{M_{11} - M_{22}}{\sqrt{\Delta}} \right)}, \tag{22}
\end{align}$$

where we introduce the notation

$$\Delta = (M_{11} - M_{22})^2 + 4M_{12}^2. \quad (23)$$

The eigenvalue of the $qS_2$-wave is given by eq. (20) with a negative sign in front of the square root. The corresponding vector $\mathbf{g}_{S_2}^{(0)}$ has the following components:

$$\begin{align}
\alpha_{S_2}^1 &= -\operatorname{sign}(M_{12}) \sqrt{\frac{1}{2} \left( 1 - \frac{M_{11} - M_{22}}{\sqrt{\Delta}} \right)} \nonumber \\
\alpha_{S_2}^2 &= \sqrt{\frac{1}{2} \left( 1 + \frac{M_{11} - M_{22}}{\sqrt{\Delta}} \right)}. \tag{24}
\end{align}$$

Some instabilities may appear in the computation of the vectors $\mathbf{g}_{S_1}^{(0)}$ and $\mathbf{g}_{S_2}^{(0)}$ in regions in which the two quasi-shear waves propagate with nearly the same velocity.

Let us remark that the first-order expressions (14) and (20) of the eigenvalues are independent of the choice of the reference medium. Moreover, the zero order eigenvectors $\mathbf{g}_{S_1}^{(0)}$ and $\mathbf{g}_{S_2}^{(0)}$ are specified uniquely by the Christoffel matrix $\Gamma$ and satisfy the
relations:

$$g_{S_1}^{(0)t} \Gamma g_{S_2}^{(0)} = g_{S_1}^{(0)t} \Delta \Gamma g_{S_2}^{(0)} = 0, \qquad (25)$$

and

$$G_{S_k}^{(1)} = g_{S_k}^{(0)t} \Gamma g_{S_k}^{(0)}, \quad k = 1, 2. \qquad (26)$$

For the clarity of the following parts, let us remark that the quantities $g_{S_n}^{(0)t} \Gamma g_{P}^{(0)}$, $g_{S_n}^{(0)t} \Gamma g_{P}^{(0)}$ ($n=1, 2$) and $g_{1}^{(0)t} \Gamma g_{2}^{(0)}$ are first-order terms.

## 4 HIGH ORDER PERTURBATIONS FOR THE qP WAVE

In order to obtain a general expression of the eigenvalue $G_P$, it is convenient to write eq. (2) in the base ($g_P^{(0)}, g_{S_1}^{(0)}, g_{S_2}^{(0)}$). The scalar products of (2), written for $m=P$, with $g_{S_n}^{(0)}, n=1, 2$ give:

$$g_{S_1}^{(0)} \cdot g_P = \frac{g_{S_1}^{(0)t} \Gamma g_P^{(0)}}{G_P - G_{S_1}^{(1)}} \cdot g_P^{(0)} \cdot g_P$$ $$g_{S_2}^{(0)} \cdot g_P = \frac{g_{S_2}^{(0)t} \Gamma g_P^{(0)}}{G_P - G_{S_2}^{(1)}} \cdot g_P^{(0)} \cdot g_P \qquad (27)$$

where $G_{S_1}^{(1)}$ and $G_{S_2}^{(1)}$ are defined by eq. (26) and we use the property (25).

By definition, $g_P$ is a unit vector, $g_P \cdot g_P = 1$; therefore, using (27), one obtains:

$$g_P^{(0)} \cdot g_P = \frac{1}{\sqrt{1 + \left( \frac{g_{S_1}^{(0)t} \Gamma g_P^{(0)}}{G_P - G_{S_1}^{(1)}} \right)^2 + \left( \frac{g_{S_2}^{(0)t} \Gamma g_P^{(0)}}{G_P - G_{S_2}^{(1)}} \right)^2}} \qquad (28)$$

The scalar product of (2), written for $m=P$, with $g_P^{(0)}$ gives:

$$(G_P - G_P^{(1)}) g_P^{(0)} \cdot g_P = g_P^{(0)t} \Gamma g_{S_1}^{(0)} g_{S_1}^{(0)} \cdot g_P + g_P^{(0)t} \Gamma g_{S_2}^{(0)} g_{S_2}^{(0)} \cdot g_P, \qquad (29)$$

where $G_P^{(1)}$ is given by (14).

Inserting relations (27) into (29), one finds:

$$G_P = G_P^{(1)} + \frac{(g_{S_1}^{(0)t} \Gamma g_P^{(0)})^2}{G_P - G_{S_1}^{(1)}} + \frac{(g_{S_2}^{(0)t} \Gamma g_P^{(0)})^2}{G_P - G_{S_2}^{(1)}} \qquad (30)$$

Let us remark that the eigenvalue $G_P$ appears on both sides of the eq. (30). From eq. (30), one can obtain $G_P$ by an iterative procedure or a perturbation approach at any order. In the perturbation approach, $G_P$ and $\Gamma$ are expanded into perturbation series on both sides of equation (30). Taking into account that $g_{S_n}^{(0)t} \Gamma g_P^{(0)}$, $n=1, 2$, are first-order terms, one gets $G_P^{(1)}$ as the first-order expression of $G_P$. From (30), we notice that the first-order expression $G_P^{(1)}$ is always smaller than the exact eigenvalue $G_P$.

The second-order expression of $G_P$, denoted by $G_P^{(2)}$, is

$$G_P^{(2)} = G_P^{(1)} + \frac{(g_{S_1}^{(0)t} \Gamma g_P^{(0)})^2}{G_P^{(0)} - G_S^{(0)}} + \frac{(g_{S_2}^{(0)t} \Gamma g_P^{(0)})^2}{G_P^{(0)} - G_S^{(0)}} \qquad (31)$$

which can also be written using the vectors $g_1^{(0)}$ and $g_2^{(0)}$:

$$G_P^{(2)} = G_P^{(1)} + \frac{(g_1^{(0)t} \Gamma g_P^{(0)})^2}{G_P^{(0)} - G_S^{(0)}} + \frac{(g_2^{(0)t} \Gamma g_P^{(0)})^2}{G_P^{(0)} - G_S^{(0)}} \qquad (32)$$

The third-order expression of $G_P$ can be written as

$$G_P^{(3)} = G_P^{(1)} + \frac{(g_{S_1}^{(0)t} \Gamma g_P^{(0)})^2}{G_P^{(1)} - G_{S_1}^{(1)}} + \frac{(g_{S_2}^{(0)t} \Gamma g_P^{(0)})^2}{G_P^{(1)} - G_{S_2}^{(1)}} \qquad (33)$$

It should be emphasized that eq. (33) is not exactly an expansion into perturbation series because the denominators contain perturbed terms $G_P^{(1)} - G_{S_1}^{(1)}$ and $G_P^{(1)} - G_{S_2}^{(1)}$. The perturbation series may be obtained by making approximate expansion of the denominators. It is clear that this should give a more complicated formula and less convenient calculation.

The expression of the qP-wave eigenvalue at order $l, l \ge 3$, denoted by $G_P^{(l)}$, can be obtained via the following recursion relation between the $(l-2)$ approximation and the $l$th:

$$G_P^{(l)} = G_P^{(1)} + \frac{(g_{S_1}^{(0)t} \Gamma g_P^{(0)})^2}{G_P^{(l-2)} - G_{S_1}^{(1)}} + \frac{(g_{S_2}^{(0)t} \Gamma g_P^{(0)})^2}{G_P^{(l-2)} - G_{S_2}^{(1)}} \qquad (34)$$

In terms of iterative procedure, the expressions obtained at uneven order (even order, respectively) correspond to the sequence obtained by inserting the approximate trial solution $G_P = G_P^{(1)}$ ($G_P = G_P^{(0)}$, respectively) in the right-hand side of relation (30) and using the recursion relations (34). The expressions obtained at uneven order are independent of the choice of the reference medium. This is not the case for the expressions obtained at even order.

Using the expression of $G_P$ obtained at order $l$ in (27) and (28), one obtains the vector $g_P$ at order $l+1$. For example, the first-order expression of $g_P$ is:

$$g_P^{(1)} = \frac{1}{\sqrt{1 + \left( \frac{g_{S_1}^{(0)t} \Gamma g_P^{(0)}}{G_P^{(0)} - G_S^{(0)}} \right)^2 + \left( \frac{g_{S_2}^{(0)t} \Gamma g_P^{(0)}}{G_P^{(0)} - G_S^{(0)}} \right)^2}} \times \left[ g_P^{(0)} + \frac{g_{S_1}^{(0)t} \Gamma g_P^{(0)}}{G_P^{(0)} - G_S^{(0)}} g_P^{(0)} + \frac{g_{S_2}^{(0)t} \Gamma g_P^{(0)}}{G_P^{(0)} - G_S^{(0)}} g_P^{(0)} \right], \qquad (35)$$

which can also be written in the base ($g_P^{(0)}, g_1^{(0)}, g_2^{(0)}$):

$$g_P^{(1)} = \frac{1}{\sqrt{1 + \left( \frac{g_1^{(0)t} \Gamma g_P^{(0)}}{G_P^{(0)} - G_S^{(0)}} \right)^2 + \left( \frac{g_2^{(0)t} \Gamma g_P^{(0)}}{G_P^{(0)} - G_S^{(0)}} \right)^2}} \times \left[ g_P^{(0)} + \frac{g_1^{(0)t} \Gamma g_P^{(0)}}{G_P^{(0)} - G_S^{(0)}} g_1^{(0)} + \frac{g_2^{(0)t} \Gamma g_P^{(0)}}{G_P^{(0)} - G_S^{(0)}} g_2^{(0)} \right]. \qquad (36)$$

The $(l+1)$th order expression of $g_P$, $l \ge 1$, is given by:

$$g_P^{(l+1)} = \frac{1}{\sqrt{1 + \left( \frac{g_{S_1}^{(0)t} \Gamma g_P^{(0)}}{G_P^{(l)} - G_{S_1}^{(1)}} \right)^2 + \left( \frac{g_{S_2}^{(0)t} \Gamma g_P^{(0)}}{G_P^{(l)} - G_{S_2}^{(1)}} \right)^2}} \times \left[ g_P^{(0)} + \frac{g_{S_1}^{(0)t} \Gamma g_P^{(0)}}{G_P^{(l)} - G_{S_1}^{(1)}} g_{S_1}^{(0)} + \frac{g_{S_2}^{(0)t} \Gamma g_P^{(0)}}{G_P^{(l)} - G_{S_2}^{(1)}} g_{S_2}^{(0)} \right]. \qquad (37)$$

The even order approximations of $g_P$ depend on the first-order approximation of the eigenvalues via the recursion relation (34) and are independent of the choice of the reference medium.
## 5 HIGH ORDER PERTURBATIONS FOR THE qS-WAVES

In order to obtain a general expression of the qS-waves eigenvalues $G_{S_k}$, $k=1, 2$, one can write eqs (2) in the base $(g_P^{(0)}, g_1^{(0)}, g_2^{(0)})$. The scalar product of eq. (2), written for $m=S_k$, with $g_P^{(0)}$ gives:

$$ (G_{S_k} - G_P^{(1)}) g_P^{(0)} \cdot g_{S_k} = g_P^{(0)\prime} \Gamma g_1^{(0)} g_1^{(0)} \cdot g_{S_k} + g_P^{(0)\prime} \Gamma g_2^{(0)} g_2^{(0)} \cdot g_{S_k}, \quad (38) $$

where $G_P^{(1)}$ is given by eq. (14).

The scalar products of eq. (2), written for $m=S_k$, with $g_n^{(0)}$, $n=1, 2$ give:

$$ \begin{aligned} & (g_1^{(0)\prime} \Gamma g_1^{(0)} - G_{S_k}) g_1^{(0)} \cdot g_{S_k} + g_1^{(0)\prime} \Gamma g_2^{(0)} g_2^{(0)} \cdot g_{S_k} \\ & + g_1^{(0)\prime} \Gamma g_P^{(0)} g_P^{(0)} \cdot g_{S_k} = 0 \\ & g_2^{(0)\prime} \Gamma g_1^{(0)} g_1^{(0)} \cdot g_{S_k} + (g_2^{(0)\prime} \Gamma g_2^{(0)} - G_{S_k}) g_2^{(0)} \cdot g_{S_k} \\ & + g_2^{(0)\prime} \Gamma g_P^{(0)} g_P^{(0)} \cdot g_{S_k} = 0. \end{aligned} \quad (39) $$

Taking relation (38) into account, one obtains from (39) the following system:

$$ \begin{aligned} & (M_{11} - G_{S_k}) g_1^{(0)} \cdot g_{S_k} + M_{12} g_2^{(0)} \cdot g_{S_k} = 0 \\ & M_{12} g_1^{(0)} \cdot g_{S_k} + (M_{22} - G_{S_k}) g_2^{(0)} \cdot g_{S_k} = 0 \end{aligned} \quad (40) $$

where

$$ M_{ij} = g_i^{(0)\prime} \Gamma g_j^{(0)} + \frac{g_P^{(0)\prime} \Gamma g_i^{(0)} g_P^{(0)\prime} \Gamma g_j^{(0)}}{G_{S_k} - G_P^{(1)}}. \quad (41) $$

The condition of solvability of system (40) yields:

$$ G_{S_k} = \frac{1}{2} \left[ M_{11} + M_{22} \pm \sqrt{(M_{11} - M_{22})^2 + 4M_{12}^2} \right]. \quad (42) $$

The eigenvalue of the $qS_1$-wave (respectively, of the $qS_2$-wave) is given by eq. (42) with a positive sign (respectively negative sign) in front of the square root. Let us remark that the eigenvalue $G_{S_k}$ appears on both sides of eq. (42), because of the definition of $M_{ij}$.

From expression (42), one can obtain $G_{S_k}$ by using an iterative procedure or a perturbation approach at any order. In the perturbation approach, $G_{S_k}$ and $\Gamma$ are expanded into perturbation series in the expression of $M_{ij}$ given by (41). The expression of $G_{S_k}$ at order $l$, denoted by $G_{S_k}^{(l)}$, is given by eq. (42) with the $M_{ij}$ replaced by their $l$th order expressions. Therefore, the first-order expression $G_{S_k}^{(1)}$ is given by eq. (42) with

$$ M_{ij}^{(1)} = g_i^{(0)\prime} \Gamma g_j^{(0)}, \quad (43) $$

which is identical to expressions (20–21).

The second-order expression $G_{S_k}^{(2)}$ is given by eq. (42) with

$$ M_{ij}^{(2)} = g_i^{(0)\prime} \Gamma g_j^{(0)} + \frac{g_P^{(0)\prime} \Gamma g_i^{(0)} g_P^{(0)\prime} \Gamma g_j^{(0)}}{G_S^{(0)} - G_P^{(0)}}. \quad (44) $$

The third-order expression $G_{S_k}^{(3)}$ is given by eq. (42) with

$$ M_{ij}^{(3)} = g_i^{(0)\prime} \Gamma g_j^{(0)} + \frac{g_P^{(0)\prime} \Gamma g_i^{(0)} g_P^{(0)\prime} \Gamma g_j^{(0)}}{G_{S_k}^{(1)} - G_P^{(1)}}. \quad (45) $$

The expression at order $l$, $l \ge 3$, of the $qS_k$-wave eigenvalue, denoted by $G_{S_k}^{(l)}$, is given by eq. (42) with

$$ M_{ij}^{(l)} = g_i^{(0)\prime} \Gamma g_j^{(0)} + \frac{g_P^{(0)\prime} \Gamma g_i^{(0)} g_P^{(0)\prime} \Gamma g_j^{(0)}}{G_{S_k}^{(l-2)} - G_P^{(1)}}. \quad (46) $$

It should be emphasized that eq. (46) is not exactly an expansion into perturbation series because the denominator contains the perturbed term $G_{S_k}^{(l-2)} - G_P^{(1)}$. However, expression (46) is equivalent to the expansion into perturbation series at order $l$ and is much more convenient for calculation.

This is a recursion relation between the $(l-2)$ approximation and the $l$th. In terms of iterative procedure, the expressions obtained at uneven order (even order, respectively) correspond to the sequence obtained by inserting the approximate trial solution $G_{S_k}^{(1)}$ ($G_{S_k}^{(0)}$, respectively) in the right-hand side of relation (42) and using the iteration process via the recursion relations (42) and (46). The expressions obtained at uneven order are independent of the choice of the reference medium. This is not the case for the expressions obtained at even order.

Let us now consider the expression of the eigenvector $g_{S_k}$. By definition, $g_{S_k}$ is a unit vector, therefore

$$ (g_1^{(0)} \cdot g_{S_k})^2 + (g_2^{(0)} \cdot g_{S_k})^2 + (g_P^{(0)} \cdot g_{S_k})^2 = 1. \quad (47) $$

The approximation of the eigenvector $g_{S_k}$ at order $l$, denoted by $g_{S_k}^{(l)}$, can be obtained by using the approximations of order $(l+1)$, $G_{S_k}^{(l+1)}$ and $M_{ij}^{(l+1)}$, in relations (38) and (40), taking into account relation (47). The expressions obtained at even order are independent of the choice of the reference medium. Note that the evaluation of the term $g_P^{(0)} \cdot g_{S_k}^{(0)}$ requires only the approximation $G_{S_k}^{(l-1)}$. However, in practical applications, it is more convenient for calculation to use approximations computed at order $(l+1)$, $G_{S_k}^{(l+1)}$ and $M_{ij}^{(l+1)}$.

In order to get more insight into the leading perturbation orders of the polarization approximations, let us assume that the vectors $g_1^{(0)}$ and $g_2^{(0)}$ correspond to the vectors $g_{S_k}^{(0)}$ and $g_{S_2}^{(0)}$, respectively. The components of $g_{S_k}$ in the base $(g_P^{(0)}, g_S^{(0)}, g_{S_2}^{(0)})$ satisfy similar equations to (38–40), where the coefficients $M_{ij}$, denoted in this case by $M_{S_i S_j}$, are given by:

$$ M_{S_i S_j} = g_{S_i}^{(0)\prime} \Gamma g_{S_j}^{(0)} + \frac{g_P^{(0)\prime} \Gamma g_{S_i}^{(0)} g_P^{(0)\prime} \Gamma g_{S_j}^{(0)}}{G_{S_k} - G_P^{(1)}}. \quad (48) $$

Let us remark that in terms of leading perturbation order, the $M_{S_i S_j}$, $i=1, 2$, have a first-order leading term and $M_{S_1 S_2}$ has a second-order leading term because of relation (25). Therefore, if $M_{S_1 S_1}$ and $M_{S_2 S_2}$ are not equal to first-order, $G_{S_k}$ is equal to one of the $M_{S_i S_j}$ at second order (This can be obtained from relation (42) by expanding the square root into perturbation series):

$$ G_{S_1}^{(2)} = \max_{i=1,2} M_{S_i S_1}^{(2)} \quad (49) $$

$$ G_{S_2}^{(2)} = \min_{i=1,2} M_{S_i S_2}^{(2)} \quad (50) $$
where the second-order expressions $M_{S_i S_j}^{(2)}$ are given by:

$$M_{S_i S_j}^{(2)} = G_{S_i}^{(1)} + \frac{(g_{\mathbf{P}}^{(0)})^T \Gamma g_{\mathbf{S}_i}^{(0)})^2}{G_S^{(0)} - G_P^{(0)}}. \qquad (51)$$

Let us introduce the indices $k_1, k_2$, such that $k_1 \neq k_2$ and $G_{S_k}^{(2)} = M_{S_k S_k}^{(2)}$ at second order. From system (40), one can see that $g_{S_k}^{(0)} \cdot g_{S_k}^{(0)}$ should have a zero order leading term and $g_{S_k}^{(0)} \cdot g_{S_k}^{(0)}$ should have a first-order leading term. From system (40), one gets:

$$g_{S_k}^{(0)} \cdot g_{S_k} = \frac{M_{S_1 S_2}}{G_{S_k} - M_{S_k S_k}} g_{S_k}^{(0)} \cdot g_{S_k} \qquad (52)$$

Introducing eq. (52) into (38), one obtains:

$$g_{\mathbf{P}}^{(0)} \cdot g_{\mathbf{S}_k} = \frac{g_{\mathbf{P}}^{(0)} \Gamma g_{\mathbf{S}_k}^{(0)}}{G_{S_k} - G_{\mathbf{P}}^{(1)}} \frac{G_{S_k} - G_{S_k}^{(1)}}{G_{S_k} - M_{S_k S_k}} g_{\mathbf{S}_k}^{(0)} \cdot g_{\mathbf{S}_k}. \qquad (53)$$

From the requirement that $g_{\mathbf{S}_k}$ is a unit vector, one gets the component $g_{S_k}^{(0)} \cdot g_{S_k}$ by inserting eqs (52) and (53) into the relation:

$$(g_{S_k}^{(0)} \cdot g_{S_k})^2 + (g_{S_k}^{(0)} \cdot g_{S_k})^2 + (g_{\mathbf{P}}^{(0)} \cdot g_{\mathbf{S}_k})^2 = 1. \qquad (54)$$

The component of $g_{\mathbf{S}_k}$ along the vector $g_{S_k}^{(0)}$ is inversely proportional to the factor $G_{S_k} - M_{S_k S_k}$ which is close in absolute value to $G_{S_1} - G_{S_2}$. High order effects should be observed in regions in which the two quasi-shear waves propagate with nearly the same phase velocity. In these regions, the accuracy of the eigenvector approximation will be strongly dependent on the perturbation order. It is important to remark that the index $k_1$ may not be equal to $k$ in some region, so that the eigenvector, say $g_{S_k}$, can have a larger component along $g_{S_k}^{(0)}$ than along $g_{S_k}^{(0)}$. In these regions, $M_{S_k S_k}^{(2)}$ is larger than $M_{S_k S_k}^{(2)}$, though $G_{S_1}^{(1)}$ is larger than $G_{S_2}^{(1)}$. Such an effect can only be seen at an order $l \ge 2$.

# 6 ORTHORHOMBIC AND TRANSVERSELY ISOTROPIC SYMMETRIES

## 6.1 Orthorhombic symmetry

We shall focus our analysis on the orthorhombic symmetry. Orthorhombic media (Schoenberg & Helbig 1997) arise, for example, from combinations of vertically aligned cracks and TIV anisotropy (transverse isotropy with vertical symmetry axis). An orthorhombic medium is defined by nine independent density normalized elastic parameters $A_{I,I}$ and three mutually perpendicular planes of symmetry. We shall use the coordinate system $(x, y, z)$ with coordinate planes coinciding with the symmetry planes of anisotropy. In such a coordinate system, called the ‘crystal’ coordinate system, the elastic matrix is given by:

$$A = \begin{pmatrix} A_{11} & A_{12} & A_{13} & 0 & 0 & 0 \\ A_{12} & A_{22} & A_{23} & 0 & 0 & 0 \\ A_{13} & A_{23} & A_{33} & 0 & 0 & 0 \\ 0 & 0 & 0 & A_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & A_{55} & 0 \\ 0 & 0 & 0 & 0 & 0 & A_{66} \end{pmatrix} \qquad (55)$$

and the elements of the Christoffel matrix $\Gamma$ are given by:

$$
\begin{align}
\Gamma_{11} &= A_{11} p_x^2 + A_{66} p_y^2 + A_{55} p_z^2, \notag \\
\Gamma_{22} &= A_{66} p_x^2 + A_{22} p_y^2 + A_{44} p_z^2, \notag \\
\Gamma_{33} &= A_{55} p_x^2 + A_{44} p_y^2 + A_{33} p_z^2, \notag \\
\Gamma_{12} &= (A_{12} + A_{66}) p_x p_y, \notag \\
\Gamma_{13} &= (A_{13} + A_{55}) p_x p_z, \notag \\
\Gamma_{23} &= (A_{23} + A_{44}) p_y p_z. \tag{56}
\end{align}
$$

Let us introduce the following parameters (Mensch & Farra 1999):

$$\begin{align}
\hat{A}_{12} &= A_{12} - \frac{A_{11} + A_{22}}{2} + 2A_{66}, \quad \hat{A}_{13} = A_{13} - \frac{A_{11} + A_{33}}{2} + 2A_{55}, \notag \\
\hat{A}_{23} &= A_{23} - \frac{A_{22} + A_{33}}{2} + 2A_{44}. \tag{57}
\end{align}$$

Using the expressions (9) and (10) of the vectors $g_{\mathbf{P}}^{(0)}$, $g_{\mathbf{I}}^{(0)}$ and $g_{\mathbf{S}}^{(0)}$, one can write:

$$\begin{align}
g_{\mathbf{P}}^{(0)\mathrm{T}} \Gamma g_{\mathbf{P}}^{(0)} &= A_{11} p_x^2 + A_{22} p_y^2 + A_{33} p_z^2 + 2 \hat{A}_{12} \frac{p_x^2 p_y^2}{\mathbf{p}^2} \notag \\
&\quad + 2 \hat{A}_{13} \frac{p_x^2 p_z^2}{\mathbf{p}^2} + 2 \hat{A}_{23} \frac{p_y^2 p_z^2}{\mathbf{p}^2}, \tag{58}
\end{align}$$

$$\begin{align}
g_{\mathbf{P}}^{(0)\mathrm{T}} \Gamma g_{\mathbf{I}}^{(0)} &= \frac{-p_z}{\mathbf{p}^2 p_r} \left[ 2 \hat{A}_{12} p_x^2 p_y^2 + (\hat{A}_{13} p_x^2 + \hat{A}_{23} p_y^2)(p_z^2 - p_r^2) \right. \notag \\
&\qquad \left. + \frac{\mathbf{p}^2}{2} ((A_{11} - A_{33})p_x^2 + (A_{22} - A_{33})p_y^2) \right], \tag{59}
\end{align}$$

$$\begin{align}
g_{\mathbf{P}}^{(0)\mathrm{T}} \Gamma g_{\mathbf{I}}^{(0)} &= \frac{p_x p_y}{2 p p_r} [2 \hat{A}_{12} (p_x^2 - p_y^2) + 2 (\hat{A}_{23} - \hat{A}_{13}) p_z^2 \notag \\
&\quad + (A_{22} - A_{11}) \mathbf{p}^2], \tag{60}
\end{align}$$

$$\begin{align}
g_{\mathbf{I}}^{(0)\mathrm{T}} \Gamma g_{\mathbf{I}}^{(0)} &= \left[ 2 \hat{A}_{12} \frac{p_x^2 p_y^2}{\mathbf{p}^2 \mathbf{p}_r^2} - 2 (\hat{A}_{13} p_x^2 + \hat{A}_{23} p_y^2) \frac{p_z^2}{\mathbf{p}^2} \right. \notag \\
&\qquad \left. + (A_{55} p_x^2 + A_{44} p_y^2) \frac{\mathbf{p}^2}{\mathbf{p}_r^2} \right], \tag{61}
\end{align}$$

$$g_{\mathbf{I}}^{(0)\mathrm{T}} \Gamma g_{\mathbf{S}}^{(0)} = \frac{p_x p_y p_z}{\mathbf{p}_r^2 p} [(\hat{A}_{23} - \hat{A}_{13})\mathbf{p}_r^2 - \hat{A}_{12}(p_x^2 - p_y^2) + (A_{55} - A_{44})\mathbf{p}^2], \qquad (62)$$

$$g_{\mathbf{S}}^{(0)\mathrm{T}} \Gamma g_{\mathbf{I}}^{(0)} = \left[ A_{66}(p_x^2 + p_y^2) + (A_{44}p_x^2 + A_{55}p_y^2)\frac{p_z^2}{\mathbf{p}_r^2} - 2\hat{A}_{12}\frac{p_x^2 p_y^2}{\mathbf{p}_r^2} \right]. \qquad (63)$$

Some of the perturbation expressions depend on the reference eigenvalues $G_{S_k}^{(0)}$ and $G_{S_k}^{(0)}$. Among the isotropic reference media, the so-called isotropic replacement medium (IRM) of the orthorhombic medium (Fedorov 1968; Sayers 1994; Mensch & Rasolofosaon 1997) can be used. Let us introduce the ‘intrinsic’ deviation from a reference medium $\frac{\|A - A_0\|}{\|A\|}$ (see definition in Arts et al. (1991) and Mensch & Rasolofosaon (1997)), where $A$ and $A_0$ are the elastic tensors of the orthorhombic medium and the reference medium, respectively. $\| \cdot \|$ denotes the norm
of the considered fourth-rank tensor (for example, $||\mathbf{A}||^2 = A_{ijkl}A_{ijkl}$, with implicit summation on repeated indices). The so-called isotropic replacement medium (IRM) of the orthorhombic medium is a solution of the minimization of the norm $||\mathbf{A} - \mathbf{A}_{\text{IRM}}||$. In the IRM medium, the eigenvalues $G_P^{(0)}$ and $G_S^{(0)}$ are given by eqs (7–8) with:

$$ 
\begin{aligned}
A_{33}^{(0)} &= (3A_{11} + 3A_{22} + 3A_{33} + 2A_{12} + 2A_{13} + 2A_{23} + 4A_{44} \\
&\quad + 4A_{55} + 4A_{66})/15, \\
A_{44}^{(0)} &= (A_{11} + A_{22} + A_{33} - A_{12} - A_{13} - A_{23} + 3A_{44} \\
&\quad + 3A_{55} + 3A_{66})/15,
\end{aligned} 
\quad (64) $$

which can be also written:

$$ 
\begin{aligned}
A_{33}^{(0)} &= (A_{11} + A_{22} + A_{33})/3 + 2(\hat{A}_{12} + \hat{A}_{13} + \hat{A}_{23})/15, \\
A_{44}^{(0)} &= (A_{44} + A_{55} + A_{66})/3 - (\hat{A}_{12} + \hat{A}_{13} + \hat{A}_{23})/15.
\end{aligned} 
\quad (65) $$

Using expressions (58–63) in relations (30) and (42), one can calculate at any order the phase velocity squared $V_m^2(\mathbf{n})$ of the three waves in the direction defined by the unit vector $\mathbf{n}$ (see eq. 4). In a weakly orthorhombic medium, the qP-wave phase slowness surface is controlled by the 6 parameters $A_{11}$, $A_{22}$, $A_{33}$, $\hat{A}_{12}$, $\hat{A}_{13}$, $\hat{A}_{23}$ in its first-order expression (14) and by seven parameters in its second-order expression (32). The qS-waves phase slowness surfaces are controlled by the six parameters $A_{44}$, $A_{55}$, $A_{66}$, $\hat{A}_{12}$, $\hat{A}_{13}$, $\hat{A}_{23}$ in their first-order expressions (given by eqs 20 and 21) and by nine parameters in their second-order expressions (given by eqs 42 and 44).

## 6.2 Transverse isotropy

In a transversely isotropic medium with symmetry axis along the z-axis, one has $\hat{A}_{12}=0$, $\hat{A}_{13}=\hat{A}_{23}$, $A_{11}=A_{22}$ and $A_{44}=A_{55}$. Therefore, $g_P^{(0)}g_2^{(0)} = g_1^{(0)}g_2^{(0)} = 0$. The eigenvalue corresponding to one of the quasi-S-waves (the so-called qSH-wave) is given by:

$$ G_{SH}(\mathbf{p}) = A_{66}\mathbf{p}_r^2 + A_{44}\mathbf{p}_z^2, \quad (66) $$

the corresponding eigenvector is $g_2^{(0)}$.

The first-order expressions of the qP wave and the other qS-wave (the so-called qSV-wave) eigenvalues are given by:

$$ 
\begin{aligned}
G_P^{(1)}(\mathbf{p}) &= A_{11}\mathbf{p}_r^2 + A_{33}\mathbf{p}_z^2 + 2\hat{A}_{13} \frac{\mathbf{p}_r^2 \mathbf{p}_z^2}{\mathbf{p}^2}, \\
G_{SV}^{(1)}(\mathbf{p}) &= A_{44}\mathbf{p}^2 - 2\hat{A}_{13} \frac{\mathbf{p}_r^2 \mathbf{p}_z^2}{\mathbf{p}^2}.
\end{aligned} 
\quad (67) $$

The third-order expressions are given by:

$$ 
\begin{aligned}
G_P^{(3)}(\mathbf{p}) &= G_P^{(1)}(\mathbf{p}) + \frac{\mathbf{p}_r^2 \mathbf{p}_z^2}{G_P^{(1)}(\mathbf{p}) - G_{SV}^{(1)}(\mathbf{p})} \\
&\quad \times \left[ \hat{A}_{13} \frac{\mathbf{p}_z^2 - \mathbf{p}_r^2}{\mathbf{p}^2} + \frac{1}{2} (A_{11} - A_{33}) \right]^2, \\
G_{SV}^{(3)}(\mathbf{p}) &= G_{SV}^{(1)}(\mathbf{p}) - \frac{\mathbf{p}_z^2 \mathbf{p}_r^2}{G_P^{(1)}(\mathbf{p}) - G_{SV}^{(1)}(\mathbf{p})} \\
&\quad \times \left[ \hat{A}_{13} \frac{\mathbf{p}_z^2 - \mathbf{p}_r^2}{\mathbf{p}^2} + \frac{1}{2} (A_{11} - A_{33}) \right]^2.
\end{aligned} 
\quad (68) $$

The second-order expressions of the eigenvectors $g_P$ and $g_{SV}$ are parallel to the vectors (see expressions 37 and 38)

$$ g_P^{(0)} - \frac{p_z p_r}{G_P^{(1)}(\mathbf{p}) - G_{SV}^{(1)}(\mathbf{p})} \left[ \hat{A}_{13} \frac{p_z^2 - \mathbf{p}_r^2}{\mathbf{p}^2} + \frac{1}{2} (A_{11} - A_{33}) \right] g_I^{(0)} \quad (69) $$

and

$$ g_I^{(0)} + \frac{p_z p_r}{G_P^{(1)}(\mathbf{p}) - G_{SV}^{(1)}(\mathbf{p})} \left[ \hat{A}_{13} \frac{p_z^2 - \mathbf{p}_r^2}{\mathbf{p}^2} + \frac{1}{2} (A_{11} - A_{33}) \right] g_P^{(0)}, \quad (70) $$

respectively.

The qP and qSV-waves velocities squared in the phase normal direction defined by the unit vector $\mathbf{n}$ are given by eq. (4). Denoting by $\theta$ the angle between the phase propagation direction and the z-axis, one gets from eq. (67) the first-order expressions of the velocities squared:

$$ 
\begin{aligned}
V_P^{2(1)}(\theta) &= A_{11} \sin^2 \theta + A_{33} \cos^2 \theta + 2 \hat{A}_{13} \sin^2 \theta \cos^2 \theta, \\
V_{SV}^{2(1)}(\theta) &= A_{44} - 2 \hat{A}_{13} \sin^2 \theta \cos^2 \theta,
\end{aligned} 
\quad (71) $$

and from eq. (68) the third order expressions of the velocities squared:

$$ 
\begin{aligned}
V_P^{2(3)}(\theta) &= V_P^{2(1)}(\theta) + \frac{\sin^2 \theta \cos^2 \theta}{V_P^{2(1)}(\theta) - V_{SV}^{2(1)}(\theta)} \\
&\quad \times \left[ \hat{A}_{13}(\cos^2 \theta - \sin^2 \theta) + \frac{1}{2} (A_{11} - A_{33}) \right]^2, \\
V_{SV}^{2(3)}(\theta) &= V_{SV}^{2(1)}(\theta) - \frac{\sin^2 \theta \cos^2 \theta}{V_P^{2(1)}(\theta) - V_{SV}^{2(1)}(\theta)} \\
&\quad \times \left[ \hat{A}_{13}(\cos^2 \theta - \sin^2 \theta) + \frac{1}{2} (A_{11} - A_{33}) \right]^2.
\end{aligned} 
\quad (72) $$

## 7 NUMERICAL EXAMPLE

We consider the orthorhombic medium with matrix of density normalized elastic parameters (in $\text{km}^2 \text{s}^{-2}$):

$$ A = \begin{pmatrix}
10.8 & 2.2 & 1.9 & 0.0 & 0.0 & 0.0 \\
2.2 & 11.3 & 1.7 & 0.0 & 0.0 & 0.0 \\
1.9 & 1.7 & 8.5 & 0.0 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 3.6 & 0.0 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 3.9 & 0.0 \\
0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 4.3
\end{pmatrix}. \quad (73) $$

The corresponding IRM medium has elastic parameters $A_{33}^{(0)}=10.04 \text{ km}^2\text{/s}^2$ and $A_{44}^{(0)}=4.01 \text{ km}^2\text{/s}^2$. We recall that the IRM medium is the closest isotropic approximation of the anisotropic model. The ‘intrinsic’ deviation of the orthorhombic medium (73) with respect to the IRM medium is 10.6 per cent.

Let us illustrate the accuracy of the approximate formulae (14) and (31–34) of the qP-wave phase velocity squared. The phase velocity squared in the phase normal direction defined by the unit vector $\mathbf{n}$ is given by eq. (4). The vector $\mathbf{n}$ is defined by its polar angle $\theta$ ($\theta=0^{\circ}$ for propagation along the z-axis)
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Map of the exact *qP*-wave phase velocity (in km s⁻¹) as a function of the polar angle θ and the azimuth angle φ of the phase propagation direction. The model is the orthorhombic medium with elastic parameters given by (73).</figcaption></figure>

and azimuth φ (φ=0° for propagation in the *Oxz* plane). Because the orthorhombic symmetry has three mutually perpendicular planes of symmetry, all the calculations were made for 0°≤θ≤90° and 0°≤φ≤90°. Let us denote by *V*<sub>*P*</sub> the true *qP*-wave phase velocity and *V*<sub>*P*</sub><sup>(l)</sup> the approximated phase velocity at order *l*. Fig. 1 represents the true *qP*-wave phase velocity as a function of the phase angles θ and φ. The maximum variation with respect to the IRM model is 8 per cent. Table 1 gives the maximum relative error

$$ \frac{V_P^{(l)} - V_P}{V_P} $$

(in per cent) for *l*=1, 2, 3, 4. The even order approximations *V*<sub>*P*</sub><sup>(2)</sup> and *V*<sub>*P*</sub><sup>(4)</sup> were computed with the reference eigenvalues *G*<sub>*P*</sub><sup>(0)</sup> and *G*<sub>*S*</sub><sup>(0)</sup> of the IRM medium. One can note the increasing accuracy with the order *l* of the approximation.

Fig. 2 shows the angular deviation of the exact *qP*-wave polarization vector from the phase propagation direction. It reaches maximum values of about 8°. Table 2 gives the maximum angular deviation γ<sub>*P*</sub><sup>(l)</sup> in degrees between the exact polarization vector and the approximated polarization vector at order *l*=0, 1, 2. The first and second-order expressions of the polarization vector are given by eqs (36) and (37),

<table><caption>**Table 1.** Maximum error (in per cent) of the *qP*-wave phase velocity as a function of the perturbation order. The steps used in the calculations were 1° in θ and φ, for 0°≤θ≤90° and 0°≤φ≤90°.</caption><thead><tr><th><sup>V<sub>P</sub><sup>(1)</sup></sup>/<sub>V<sub>P</sub></sub> × 100</th><th><sup>V<sub>P</sub><sup>(2)</sup></sup>/<sub>V<sub>P</sub></sub> × 100</th><th><sup>V<sub>P</sub><sup>(3)</sup></sup>/<sub>V<sub>P</sub></sub> × 100</th><th><sup>V<sub>P</sub><sup>(4)</sup></sup>/<sub>V<sub>P</sub></sub> × 100</th></tr></thead><tbody><tr><td>0.6</td><td>0.06</td><td>0.011</td><td>0.0009</td></tr></tbody></table>

<table><caption>**Table 2.** Maximum deviation (in degrees) between the exact and the approximated *qP*-wave eigenvector as a function of the perturbation order. The steps used in the calculations were 1° in θ and φ, for 0°≤θ≤90° and 0°≤φ≤90°.</caption><thead><tr><th>γ<sub>P</sub><sup>(0)</sup></th><th>γ<sub>P</sub><sup>(1)</sup></th><th>γ<sub>P</sub><sup>(2)</sup></th></tr></thead><tbody><tr><td>8°</td><td>1°</td><td>0.15°</td></tr></tbody></table>

respectively. The first-order approximation was computed with the reference eigenvalues of the IRM medium. The second-order approximation gives very accurate results.

Let us illustrate the accuracy of the computation of the two *qS*-waves phase velocities obtained from expression (42) where the coefficients *M*<sub>*ij*</sub> are computed with eqs (43–46), respectively. Let us denote by *V*<sub>*S*<sub>*k*</sub></sub> (*k*=1, 2) the true *qS*-waves phase velocities and *V*<sub>*S*<sub>*k*</sub></sub><sup>(l)</sup> the approximated phase velocities at order *l*. Fig. 3 represents the true *qS*<sub>*l*</sub> and *qS*<sub>*k*</sub> phase velocities as functions of the phase angle θ. Each curve corresponds to a different azimuth φ. There is a point singularity, also called single point degeneracy (Schoenberg & Helbig 1997) in the symmetry plane *Oxz* at θ=34.4° and φ=0°, where the two quasi-shear waves velocity sheets come into contact at *V*<sub>*S*<sub>*l*</sub></sub>=*V*<sub>*S*<sub>*k*</sub></sub>=1.955 km/s. The cross-section of the quasi-shear waves velocity sheets by any plane passing through the singular direction (θ=34.4°, φ=0°) gives two curves with discontinuous tangents at the point singularity, except for the plane tangent to the cone θ=34.4°, in which the curves are smooth at the point singularity. Fig. 4 represents the map of the normalized difference (in per cent) of the two exact *qS*-waves phase velocities. One can see that the two quasi-shear waves propagate with nearly the same phase velocity for

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Map of the angular deviation (in degrees) of the phase normal and the exact polarization vector *g*<sub>*P*</sub> as a function of the polar angle θ and the azimuth angle φ of the phase propagation direction. The model is the orthorhombic medium with elastic parameters given by eq. (73).</figcaption></figure>
<table><caption>**Table 3.** Maximum error (in per cent) of the $qS_1$ and $qS_2$-waves phase velocities as a function of the perturbation order. The steps used in the calculations were 1° in $\theta$ and $\phi$, for $0° \le \theta \le 90°$ and $0° \le \phi \le 90°$.</caption><thead><tr><th>S<sub>n</sub></th><th><sup>V<sub>Sn</sub><sup>(1)</sup></sup>&frasl;<sub>V<sub>Sn</sub></sub> &times; 100</th><th><sup>V<sub>Sn</sub><sup>(2)</sup></sup>&frasl;<sub>V<sub>Sn</sub></sub> &times; 100</th><th><sup>V<sub>Sn</sub><sup>(3)</sup></sup>&frasl;<sub>V<sub>Sn</sub></sub> &times; 100</th><th><sup>V<sub>Sn</sub><sup>(4)</sup></sup>&frasl;<sub>V<sub>Sn</sub></sub> &times; 100</th></tr></thead><tbody><tr><td>S<sub>1</sub></td><td>0.7</td><td>0.08</td><td>0.007</td><td>0.0008</td></tr><tr><td>S<sub>2</sub></td><td>1.4</td><td>0.15</td><td>0.027</td><td>0.002</td></tr></tbody></table>

$32° \le \theta \le 38°$ and $0° \le \phi \le 25°$. In this region, the two phase velocity sheets approach each other in a pinch. Along the pinch, the velocity sheets come into contact smoothly at the point singularity. Table 3 gives the maximum relative errors

$$ \frac{V_{S_k}^{(l)} - V_{S_k}}{V_{S_k}} $$

(in per cent) for $l=1, 2, 3, 4$. The even order approximations $V_{S_k}^{(2)}$ and $V_{S_k}^{(4)}$ were computed with the reference eigenvalues $G_P^{(0)}$ and $G_S^{(0)}$ of the IRM medium. The accuracy of the approximations increases with the order $l$.

Fig. 5 shows the map of the angular deviation between the exact polarization vector $g_{S_1}$ and the vector $g_1^{(0)}$ given by expression (10). One can see a large rotation of the direction of the polarization vector $g_{S_1}$ as the phase propagation direction

<figure><img src="image_4.png" /><figcaption>**Figure 3.** $qS_1$ and $qS_2$ true phase velocities as functions of the polar angle $\theta$ of the phase propagation direction. Each curve corresponds to a constant azimuth $\phi$, $0° \le \phi \le 90°$. The step between two azimuths is 10°. There is a point singularity in the symmetry plane $Oxz$ at $\theta=34.4°$ and $\phi=0°$, where the two phase velocity sheets come into contact at $V_S=1.955$ km/s. The model is the orthorhombic medium with elastic parameters given by (73).</figcaption></figure>

$\phi$ (deg)

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Map of the normalized difference (in per cent) of the exact $qS_1$ and $qS_2$ phase velocities as functions of the polar angle $\theta$ and the azimuth angle $\phi$ of the phase propagation direction. The model is the orthorhombic medium with elastic parameters given by (73).</figcaption></figure>

passes through the region where the two quasi-shear waves propagate with nearly the same phase velocity (see Fig. 4). As the direction of phase propagation passes through this region, the polarizations of the two quasi-shear waves sheets are exchanged. Such a phenomena is described by Crampin (1981) for orthorhombic orthopyroxene. Fig. 6 shows the

$\phi$ (deg)

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Map of the angular deviation (in degrees) of the $qS_1$-wave exact polarization vector $g_{S_1}$ and the vector $g_1^{(0)}$ given by eq. (10) as a function of the polar angle $\theta$ and the azimuth angle $\phi$ of the phase propagation direction. The steps used in the calculations were 0.75° in $\theta$ and $\phi$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Map of the angular deviation (in degrees) of the $qS_1$-wave exact polarization vector $g_{S_1}$ and the zero-order eigenvector $g_{S_1}^{(0)}$ as a function of the polar angle $\theta$ and the azimuth angle $\phi$ of the phase propagation direction. The steps used in the calculations were 0.75° in $\theta$ and $\phi$.</figcaption></figure>

map of the angular deviation between the exact eigenvector $g_{S_1}$ and the zero-order eigenvector $g_{S_1}^{(0)}$. The deviation is large for $25^\circ \le \theta \le 50^\circ$ and is above 70 degrees in the region defined by $34.4^\circ \le \theta \le 38.9^\circ$ and $0^\circ \le \phi \le 25^\circ$. In this region, the exact polarization vector $g_{S_1}$ has a larger component along the vector $g_{S_2}^{(0)}$ than along the vector $g_{S_1}^{(0)}$. Fig. 7 shows the map of the

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Map of the angular deviation (in degrees) of the zero-order eigenvector $g_{S_1}^{(0)}$ and the vector $g_{S_1}^{(0)}$ given by eq. (10) as a function of the polar angle $\theta$ and the azimuth angle $\phi$ of the phase propagation direction. The steps used in the calculations were 0.75° in $\theta$ and $\phi$.</figcaption></figure>

angular deviation between the zero-order eigenvector $g_{S_1}^{(0)}$ and the vector $g_{S_1}^{(0)}$. One can see that the region of rapid variation of the direction of the vector $g_{S_1}^{(0)}$ is different from the one in Fig. 5. In this region, the first-order velocity sheets approach each other in a pinch and come into contact in the symmetry plane $Oxz$ at $\theta=38.9^\circ$ and $\phi=0^\circ$. In the area located between these two regions (corresponding to the pinch of the exact velocity sheets and the pinch of the first-order velocity sheets, respectively), the quantity $M_{S_2 S_2}^{(2)}$ is larger than $M_{S_1 S_1}^{(2)}$ and there is a large deviation between the exact eigenvectors and the corresponding zero-order eigenvectors.

Fig. 8 shows the map of the angular deviation (in degrees) between the first-order expression of the polarization vector $g_{S_1}^{(1)}$ and the exact polarization vector $g_{S_1}$ for $15^\circ \le \theta \le 55^\circ$ and $0^\circ \le \phi \le 90^\circ$. The components of the vector $g_{S_1}^{(1)}$ in the base $(g_P^{(0)}, g_I^{(0)}, g_S^{(0)})$ were obtained by using the second-order expressions $G_{S_2}^{(2)}$ and $M_{ij}^{(2)}$, defined by eq. (44), in relations (38) and (40). The first-order approximation $g_{S_1}^{(1)}$ has less than 1 degree error outside two bounded regions. The second-order velocity sheets come into contact at $\theta=34.8^\circ$ and $\phi=0^\circ$. In the area between the pinch of the exact velocity sheets and the pinch of the approximated second-order velocity sheets, $34.4^\circ \le \theta \le 34.8^\circ$ and $0^\circ \le \phi \le 7^\circ$, there is a large deviation between the vectors $g_{S_1}^{(1)}$ and $g_{S_1}$.

The second-order expressions of the polarization vectors were obtained by using the third order expressions $G_{S_2}^{(3)}$ and $M_{ij}^{(3)}$, defined by (45), in relations (38) and (40). The deviation between the approximation $g_{S_1}^{(2)}$ and $g_{S_1}$ is shown on Fig. 9 for $31^\circ \le \theta \le 36^\circ$ and $0^\circ \le \phi \le 24^\circ$. The second-order approximation has less than 1° error outside a region corresponding to the pinch of the exact velocity sheets (see Fig. 4), with approximately 20° length and maximum width of 1°. The third order

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Map of the angular deviation (in degrees) of the $qS_1$-wave exact polarization vector $g_{S_1}$ and the first-order eigenvector $g_{S_1}^{(1)}$ in the propagation directions $15^\circ \le \theta \le 55^\circ$ and $0^\circ \le \phi \le 90^\circ$. The steps used in the calculations were 0.1° in $\theta$ and 1° in $\phi$. The reference medium is the IRM medium.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 9. Map of the angular deviation (in degrees) of the *qS*₁-wave exact polarization vector *g*<sub>S₁</sub> and the second-order eigenvector *g*<sub>S₂</sub>(2) in the propagation directions 31° ≤ *θ* ≤ 36° and 0° ≤ *φ* ≤ 24°. The steps used in the calculations were 0.05° in *θ* and 0.25° in *φ*.</figcaption></figure>

expressions of the polarization vectors has less than 1° error outside a small region, with 6° length and maximum width of 0.1°. High order expressions $G_{S_1}^{(1)}$ and $M_{ij}^{(1)}$ should be used to have good approximations of the qS-waves polarization vectors in neighbouring directions of the singularity.

## 8 CONCLUSIONS

In this paper, an iterative procedure is used to obtain approximate analytic expressions of the Christoffel matrix eigenvalues corresponding to higher order perturbations. Very simple recursion formulae are given between the eigenvalue approximations of order ($l-2$) and $l$. Approximate expressions of the corresponding eigenvectors are obtained at a shifted order ($l-1$) or ($l+1$). The eigenvalues approximations of uneven order and the eigenvectors approximations of even order are independent of the choice of the reference medium. The eigenvalues and the eigenvectors of the Christoffel matrix are associated to the squared phase velocities and the polarization vectors of the three waves (*qP*, *qS*₁, *qS*₂) propagating in an arbitrary direction of a homogeneous medium. Explicit analytic formulae of the approximate squared phase velocities and polarizations are given for orthorhombic and transversely isotropic symmetries. Inspection of the derived formulae makes it possible to study the sensitivity of phase velocities and polarization vectors to elastic parameters.

The example presented in this paper shows the accuracy of the perturbation formulae. For realistic anisotropy (10.6 per cent), the second-order expressions of the squared phase velocities are accurate approximations (errors less than 0.15 per cent if the IRM medium is used as the reference medium) which do not cost much with respect to the first-order computations. Third order expressions of the squared phase velocities are very accurate (errors less than 0.03 per cent) and need only computation of the first-order approximations. The second-order expression of the *qP*-wave polarization has error less than

0.15°. High order expressions of the *qS*-waves polarizations vectors should be used to have good approximations in the vicinity of singularities.

These analytic perturbation formulae could find important applications in the approximate evaluation of kinematic and dynamic quantities of seismic waves propagating in anisotropic media.

## ACKNOWLEDGMENTS

I thank Ivan Pšenčík for his helpful comments. I also thank the two anonymous reviewers for their useful comments which helped me to improve the paper. This is a IPGP contribution No 1757.

## REFERENCES

Arts, R.J., Helbig, K. & Rasolofosaon, P.N.J., 1991. General anisotropic tensor in rocks-approximation, invariants and particular directions, in *61st Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts*, pp. 1534-1535.

Auld, B.A., 1973. *Acoustic Fields and Waves in Solids*, Vol. 1, Wiley, New York.

Backus, G.E., 1965. Possible forms of seismic anisotropy of the uppermost mantle under oceans, *J. Geophys. Res.*, **70**, 3429-3439. (see also comments on this paper by S. Crampin (1982) and reply by the author in *J. Geophys. Res.*, **87**, 4636-4640).

Červený, V., 1972. Seismic rays and rays intensities in inhomogeneous anisotropic media, *Geophys. J. R. astr. Soc.*, **29**, 1-13.

Červený, V., 1982. Direct and inverse kinematic problem for inhomogeneous anisotropic media-linearized approach, *Contr. Geophys. Inst. Slov. Acad. Sci.*, **13**, 127-133.

Červený, V. & Jech, J., 1982. Linearized solutions of kinematic problems of seismic body waves in inhomogeneous slightly anisotropic media, *J. Geophys.*, **51**, 96-104.

Crampin, S., 1981. A review of wave motion in anisotropic and cracked elastic-media, *Wave Motion.*, **3**, 343-391.

Crampin, S. & Yedlin, M., 1981. Shear-wave singularities of wave propagation in anisotropic media, *J. Geophys.*, **49**, 43-46.

Druzhinin, A.B., 1996. Nonlinear ray perturbation theory with its applications to ray tracing and inversion in anisotropic media, *Pure appl. Geophys.*, **148**, 637-683.

Farra, V., 1989. Ray perturbation theory for heterogeneous hexagonal anisotropic medium, *Geophys. J. Int.*, **99**, 723-738.

Farra, V., 1999. Computation of second-order travel-time perturbation by Hamiltonian ray theory, *Geophys. J. Int.*, **136**, 205-217.

Fedorov, F.I., 1968. *Theory of Elastic Waves in Crystals*, Plenum Press, New York.

Hanyga, A., 1982. The kinematic inverse problem for weakly laterally inhomogeneous anisotropic media, *Tectonophysics*, **90**, 253-262.

Helbig, K., 1994. Foundation of anisotropy for exploration seismics, *Handbook of Geophysical Exploration. I. Seismic Exploration*, Vol. 22, pp. xx-xx, eds Helbig, K and Treitel, S., Pergamon Press, Oxford.

Jech, J. & Pšenčík, I., 1989. First-order approximation method for anisotropic media, *Geophys. J. Int.*, **99**, 369-376.

Landau, L.D. & Lifshitz, E.M., 1966. *Mécanique Quantique*, Mir Editions, Moscow.

Mensch, T. & Farra, V., 1999. Computation of P-wave rays, traveltimes and slowness in orthorhombic media, *Geophys. J. Int.*, **138**, 244-256.
Mensch, T. & Rasolofosaon, P., 1997. Elastic-wave velocities in anisotropic media of arbitrary symmetry—generalization of Thomsen parameters *ε*, *δ* and *Γ*, *Geophys. J. Int.*, **128**, 43–64.

Morse, P.M. & Feshbach, H., 1953. *Methods of Theoretical Physics*, McGraw-Hill, New York.

Nowack, R.L. & Pšenčík, I., 1991. Perturbation from isotropic to anisotropic heterogeneous media in ray approximation, *Geophys. J. Int.*, **106**, 1–10.

Pšenčík, I., 1998. Green functions for inhomogeneous weakly anisotropic media, *Geophys. J. Int.*, **135**, 279–288.

Pšenčík, I. & Gajewski, D., 1998. Polarization, phase velocity and NMO velocity of *qP* waves in arbitrary weakly anisotropic media, *Geophysics*, **63**, 1754–1766.

Rümpker, T. & Thomson, C.J., 1994. Seismic-waveform effects of conical points in gradually varying anisotropic media, *Geophys. J. Int.*, **118**, 759–780.

Sayers, C.M., 1994. *P*-wave propagation in weakly anisotropic media, *Geophys. J. Int.*, **116**, 799–805.

Schoenberg, M. & Helbig, K., 1997. Orthorhombic media: Modeling elastic wave behavior in vertically fractured earth, *Geophysics*, **62**, 1954–1974.

Thomsen, L., 1986. Weak anisotropy, *Geophysics*, **51**, 1954–1966.

Vavryčuk, V., 1999. Properties of *S* waves near a kiss singularity: a comparison of exact and ray solutions, *Geophys. J. Int.*, **138**, 581–589.