# First-order ray tracing for qS waves in inhomogeneous weakly anisotropic media

Véronique Farra

► To cite this version:

Véronique Farra. First-order ray tracing for qS waves in inhomogeneous weakly anisotropic media. Geophysical Journal International, 2005, 161, pp.309-324. 10.1111/j.1365-246X.2005.02570.x. insu-03601119

HAL Id: insu-03601119

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

Submitted on 8 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
# First-order ray tracing for $qS$ waves in inhomogeneous weakly 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 2004 December 23. Received 2004 October 14; in original form 2004 April 14

## SUMMARY

First-order ray tracing equations and paraxial ray tracing system are derived for quasi-shear waves propagating in a smooth inhomogeneous weakly anisotropic medium. The first-order equations are based on explicit formulae in terms of slowness vector and elastic parameters. They only depend on 15 independent combinations of elastic parameters. For isotropic media, they reduce to exact equations. The behaviour of the approximate equations is studied in the vicinity of singularities. Ray tracing equations behave properly if the wave mode is changed when the ray crosses the conical or the intersection singularity, but the paraxial ray tracing system may be singular. A second-order correction of the traveltime is obtained by integration along the approximate rays. Comparison of the exact and approximate traveltimes in transversely isotropic and orthorhombic homogeneous media shows small relative errors except in the close vicinity of cusps.

**Key words:** anisotropy, perturbation method, $qS$ waves, ray tracing, second-order traveltime.

## 1 INTRODUCTION

Perturbation theory is a useful tool for approximate study of wave properties in weakly anisotropic media. It gives approximate but very simple and transparent expressions for physical quantities such as the phase velocity and the polarization vector of individual waves (see, among others, Mensch & Rasolofosaon 1997; Pšenčík & Gajewski 1998; Farra 2001; Farra & Pšenčík 2003). Such approximations are important as they offer physical insight into the dependence of wave attributes on the elastic parameters of the medium. In anisotropic media of arbitrary symmetry, exact expressions depend on all 21 elastic parameters. In weakly anisotropic media, the individual seismic waves only depend on a limited set of their combinations, the so-called weak anisotropy parameters (see Farra & Pšenčík 2003). The combinations of elastic parameters corresponding to the $qP$ wave are different from combinations for the $qS$ waves. Standard ray tracing equations for anisotropic media do not reflect this fact since they appear identical for all three types of wave. Moreover, the standard ray equations are quite complicated. The right-hand sides of the ray tracing equations contain up to 81 terms consisting of individual elastic parameters. Many of these terms are mutually dependent and the calculation of the right-hand sides of the ray tracing equations requires many unnecessary operations.

In this paper, approximate ray tracing equations are presented for the $qS$ waves. The approximation is of the first order with respect to the size of deviation of the medium from isotropy. A similar approach is used for $qP$ wave in Pšenčík & Farra (2004). The application of the standard ray method (Červený 1972) to the propagation of the $qS$ waves in inhomogeneous anisotropic media is more complicated than for the $qP$ wave, because of difficulties related to singularities. Singularities can cause breakdowns of the ray tracing algorithms (Shearer & Chapman 1989; Gajewski & Pšenčík 1990; Vavryčuk 2001). Singularities often appear with triplications of the wave front. Triplications complicate the geometry of the wave front but do not pose complications for ray tracing equations. Moreover, in the vicinity of shear wave singularities, the two $qS$ waves do not propagate independently but are mutually coupled (see, for example, Chapman & Shearer 1989).

In Sections 2 and 3, basic equations and standard (exact) ray equations are reviewed. In Section 4, the matrix $B$, whose elements control various attributes of elastic waves, is introduced and first-order approximations for phase velocities are obtained. In Sections 5 and 6, first-order ray tracing (FORT) equations and paraxial ray tracing systems are derived for $qS$ waves. These equations are based on explicit formulae which make the dependence on parameters of the medium transparent. Their behaviour in the vicinity of singularities is studied analytically. In Section 7, a second-order traveltime correction to be evaluated along the first-order rays is derived. In Section 8, explicit expressions are given for transverse isotropic (TI) media. The accuracy of approximate traveltime formulae is illustrated with numerical examples in Section 9.

In the following, all the lowercase indices range over the values 1, 2 and 3 and the uppercase indices range over the values 1 and 2. The subscript with brackets $[m]$ characterizes the wave mode, $m = 1, 2$ for the $qS$ waves and $m = 3$ for the $qP$ wave. The superscript in parenthesis
indicates the approximation order of the quantity. Voigt notation $A_{\alpha\beta}$ for density-normalized elastic parameters, with $\alpha$ and $\beta$ running from 1 to 6, is used in parallel with the tensor notation $a_{ijkl}$. The Einstein summation convention is used for the repeated subscripts.

## 2 GENERAL EXPRESSIONS

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

$$ \Gamma_{jk} = p_i p_l a_{ijkl}. \qquad (1) $$

The matrix $\Gamma(\mathbf{p})$ is called the generalized Christoffel matrix in contrast to the standard Christoffel matrix in the definition of which the vector $\mathbf{p}$ is defined as a unit vector. 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 $3 \times 3$ matrix with three positive eigenvalues $G_{[m]}(\mathbf{p})$ and corresponding eigenvectors $g_{[m]}(\mathbf{p})$. Since the elements of matrix $\Gamma$ are homogeneous functions of the second degree in $\mathbf{p}$, $G_{[m]}(\mathbf{p})$ and $g_{[m]}(\mathbf{p})$ are homogeneous functions of degree two and zero, respectively (see Červený 2001).

Three wave modes (the quasi-$P$ wave and two quasi-$S$ waves) can propagate in the anisotropic solid defined by the density-normalized elastic parameters $a_{ijkl}$. Each wave mode is associated with one of the eigenvalues $G_{[m]}$. For each of the wave modes, the corresponding slowness vector, denoted $\mathbf{p}_{[m]}$, satisfies the eikonal equation which can be written in the following form (Červený 1972):

$$ G_{[m]}(\mathbf{p}_{[m]}) = 1. \qquad (2) $$

The slowness vector $\mathbf{p}_{[m]} = \mathbf{n}/V_{[m]}(\mathbf{n})$ is related to the phase velocity $V_{[m]}(\mathbf{n})$ in the wave normal direction defined by the unit vector $\mathbf{n}$. Since $G_{[m]}(\mathbf{p})$ is a homogeneous function of degree 2 with respect to $\mathbf{p}$, it can be deduced from eq. (2) that the phase velocity squared is given by:

$$ V_{[m]}^2(\mathbf{n}) = G_{[m]}(\mathbf{n}). \qquad (3) $$

The polarization vector of the wave is identical with the corresponding eigenvector $\mathbf{g}_{[m]}$.

In the following sections, the eigenvalues $G_{[m]}$ are ordered as follows: $G_{[2]} \le G_{[1]} \le G_{[3]}$. The eigenvalues $G_{[M]}$, $M = 1, 2$, correspond to the $qS_M$ waves, the $qS_1$ wave being defined as the faster quasi-shear wave in terms of phase velocity; the remaining eigenvalue $G_{[3]}$ belongs to the $qP$ wave. The three phase-velocity sheets can be separated, but they can also come into contact along so-called singularity directions.

Singularities are very common in all kinds of anisotropy. They are defined as directions where two waves have coincident phase velocities (or eigenvalues $G_{[m]}$). For most geological materials, the waves with a coincident phase velocity are the quasi-shear waves. We will assume that only shear wave singularities are present in the medium. Three types of singularities can be distinguished: the kiss, the intersection and the conical singularity (see, for example, Crampin & Yedlin 1981). All these singularities can appear in weakly as well as in strongly anisotropic media. In singular directions, it is not possible to specify uniquely the polarization vectors $\mathbf{g}_{[M]}$ of the $qS_M$ waves. It is only possible to find the plane in which the polarization vectors $\mathbf{g}_{[M]}$ are situated, i.e. the plane perpendicular to the third polarization vector $\mathbf{g}_{[3]}$.

## 3 RAY EQUATIONS

The ray equations discussed in the following sections can be used for any 3-D continuous distributions of the elastic parameters. The Christoffel matrix $\Gamma$, the eigenvalues $G_{[m]}$ and the eigenvectors $\mathbf{g}_{[m]}$ are functions of the position vector $\mathbf{x}$.

The ray equations can be written in a compact way by using the Hamiltonian formulation. Let us introduce the Hamiltonian:

$$ H(\mathbf{x}, \mathbf{p}) = \frac{1}{2} (G(\mathbf{x}, \mathbf{p}) - 1), \qquad (4) $$

where the symbol $G$ denotes one of three eigenvalues of the Christoffel matrix. In the following sections the subscript $[m]$, which characterizes the wave mode, will be omitted for simplification in the ray equations.

The ray tracing equations in inhomogeneous anisotropic media can be written as:

$$ \begin{aligned} \frac{dx_i}{d\tau} &= \frac{\partial H}{\partial p_i} = \frac{1}{2} \frac{\partial G}{\partial p_i}, \\ \frac{dp_i}{d\tau} &= -\frac{\partial H}{\partial x_i} = -\frac{1}{2} \frac{\partial G}{\partial x_i}, \end{aligned} \qquad (5) $$

(see, for example, Červený 1989; Farra 1989). Here $x_i$ are coordinates of the trajectory of the ray and $p_i$ are components of the slowness vector at corresponding points of the ray. The variable $\tau$ is a parameter along the ray, which has the meaning of propagation time. Let us mention that the group velocity is given by the first set of eqs (5).

At each point of a ray, the slowness vector satisfies the eikonal equation:

$$ H(\mathbf{x}, \mathbf{p}) = 0, \qquad (6) $$

which yields (2). As $H$ is constant along any solution of (5), it is sufficient to satisfy the eikonal equation (6) at the source in order to satisfy it along the whole ray.
Ray tracing equations (5) have been expressed compactly in terms of partial derivatives of the Christoffel matrix with respect to $x_i$ and $p_i$ and the polarization vector $g$ of the considered wave mode (Červený 1972; Gajewski & Pšenčík 1987, 1990):

$$ \frac{dx_i}{d\tau} = \frac{1}{2} g^T \frac{\partial \Gamma}{\partial p_i} g, \qquad \frac{dp_i}{d\tau} = -\frac{1}{2} g^T \frac{\partial \Gamma}{\partial x_i} g \qquad (7) $$

where the superscript T is used to denote the transposed vector. This requires the computation of many terms at each step of integration of the ray equations, including the 21 elastic parameters and their spatial derivatives. The right-hand side of the ray tracing equations (7) contains up to 81 terms consisting of individual elastic parameters. Many of these terms are mutually dependent, so that the calculation of the right-hand side of ray tracing equations often requires many unnecessary operations, which may even increase numerical errors.

Let us describe a ray by the position vector $\mathbf{x}(\tau)$ and the slowness vector $\mathbf{p}(\tau)$. The traveltime between $\mathbf{x}(\tau_0)$ and $\mathbf{x}(\tau)$ is obtained by integration along the ray

$$ T = \int_{\tau_0}^{\tau} \left( \mathbf{p} \cdot \frac{d\mathbf{x}}{d\tau} - H(\mathbf{x}, \mathbf{p}) \right) d\tau. \qquad (8) $$

The ray satisfies the eikonal equation (6); thus the second term in (8) is zero. Moreover, since the eigenvalue $G$ is a homogeneous function of degree 2 in $\mathbf{p}$, it can be shown that:

$$ \mathbf{p} \cdot \frac{d\mathbf{x}}{d\tau} = \frac{1}{2} p_i \frac{\partial G}{\partial p_i} = 1, \qquad (9) $$

so that (8) can be simply written as:

$$ T = \tau - \tau_0. \qquad (10) $$

Along the ray $(\mathbf{x}(\tau), \mathbf{p}(\tau))$, we can compute the paraxial propagator matrix $\mathbf{P}(\tau, \tau_0)$. The paraxial propagator matrix is the propagator matrix of the paraxial ray tracing system:

$$ \frac{d\mathbf{P}}{d\tau} = \mathbf{D}(\tau)\mathbf{P}, \qquad (11) $$

where

$$ \mathbf{D}(\tau) = \begin{pmatrix} \mathbf{S} & \mathbf{T} \\ -\mathbf{R} & -\mathbf{S}^T \end{pmatrix} \qquad (12) $$

is a 6 × 6 matrix computed on the reference ray at $\tau$. The elements of the 3 × 3 matrices $\mathbf{R}$, $\mathbf{S}$ and $\mathbf{T}$ are defined by:

$$ R_{ij} = \frac{\partial^2 H}{\partial x_i \partial x_j} = \frac{1}{2} \frac{\partial^2 G}{\partial x_i \partial x_j}, \qquad S_{ij} = \frac{\partial^2 H}{\partial p_i \partial x_j} = \frac{1}{2} \frac{\partial^2 G}{\partial p_i \partial x_j}, \qquad T_{ij} = \frac{\partial^2 H}{\partial p_i \partial p_j} = \frac{1}{2} \frac{\partial^2 G}{\partial p_i \partial p_j}. \qquad (13) $$

The propagator matrix $\mathbf{P}(\tau, \tau_0)$ is the 6 × 6 matrix with initial condition $\mathbf{P}(\tau_0, \tau_0) = \mathbf{I}_6$ ($\mathbf{I}_6$ is the 6 × 6 identity matrix). The propagator matrix of the paraxial system (11) has many applications, for example in the computation of the paraxial rays around the reference ray, the second derivatives of the traveltime and the ray theoretical amplitude, in the two-point ray tracing as well as in ray tracing in perturbed media (see, among others, Červený et al. 1988; Červený 1989; Farra & Madariaga 1987; Farra 1993).

The second partial derivatives of $G$ in (13) can be written as (see eq. 20 in Gajewski & Pšenčík 1990):

$$ \frac{\partial^2 G_{[m]}}{\partial z_i \partial z_j} = \mathbf{g}_{[m]}^T \frac{\partial^2 \Gamma}{\partial z_i \partial z_j} \mathbf{g}_{[m]} + \sum_{m'=1}^{m'=3, m' \neq m} \frac{2}{G_{[m]} - G_{[m']}} \left( \mathbf{g}_{[m]}^T \frac{\partial \Gamma}{\partial z_i} \mathbf{g}_{[m']}ight) \left( \mathbf{g}_{[m]}^T \frac{\partial \Gamma}{\partial z_j} \mathbf{g}_{[m']}ight), \qquad (14) $$

where $z_i$ may be either $x_i$ or $p_i$.

The ray tracing equations (7) and the paraxial ray tracing system (11–13) are identical in general form for all three wave modes propagating in inhomogeneous anisotropic media. The type of wave has to be specified by the initial conditions and does not change along a ray in a smooth anisotropic medium (however, it may change at interfaces); the only exception is related to the rays of $qS_1$ and $qS_2$ waves passing through shear-wave singularities. It may happen when integrating eqs (7) that the ray crosses a singularity, i.e. the direction of the local slowness vector $\mathbf{p}(\tau)$ is a singularity direction for the medium at $\mathbf{x}(\tau)$. For the $qS_M$ waves, the polarization vector cannot be determined uniquely at that point and the right-hand side of the ray equations (7) cannot be evaluated. Moreover, the paraxial ray tracing system (11–13) is singular since (14) is infinite when $G_{[1]} = G_{[2]}$. Singularities cause difficulties in tracing rays in inhomogeneous anisotropic media (Shearer & Chapman 1989; Gajewski & Pšenčík 1990; Vavryčuk 2001).

## 4 APPROXIMATE EXPRESSIONS FOR THE PHASE VELOCITY IN WEAKLY ANISOTROPIC MEDIA

In the following sections, we assume that the medium is weakly anisotropic, so that the parameters $a_{ijkl}$ can be expressed as follows:

$$ a_{ijkl} = a_{ijkl}^{(0)} + \Delta a_{ijkl} = a_{ijkl}^{(0)} + \epsilon b_{ijkl}. \qquad (15) $$
The $a_{ijkl}^{(0)}$ are the density-normalized elastic parameters in a reference isotropic medium:

$$a_{ijkl}^{(0)} = (\alpha^2 - 2\beta^2)\delta_{ij}\delta_{kl} + \beta^2(\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk}) \quad (16)$$

and $\Delta a_{ijkl} = \epsilon b_{ijkl}$ are their perturbations. $\epsilon$ is a quantity measuring the size of the deviation of the anisotropic medium from the isotropic reference medium. In (16), $\alpha$ and $\beta$ denote P- and S-wave velocities of the reference isotropic medium and $\delta_{ij}$ is the Kronecker delta.

Assuming $\epsilon$ to be a small quantity, Farra (2001) and Farra & Pšenčík (2003) used a perturbation approach to derive approximations of the phase velocities and polarization vectors for weakly anisotropic media. We briefly review some of their results which are needed in the following sections.

Let us write the vector **p** as **p** = **p**n, where **n** is a unit vector and **p** is the length of vector **p**. We introduce three mutually perpendicular unit vectors **e**_{[k]}(**n**), $k = 1, 2, 3$, where **e**_{[3]} = **n**. The vectors **e**_{[1]} and **e**_{[2]}, situated in the plane perpendicular to **n**, can be chosen arbitrarily but they should vary smoothly with variation of **n**. Moreover, we state that **e**_{[k]}(**p**) = **e**_{[k]}(**n**), so that the components of vectors **e**_{[k]} are homogeneous functions of degree 0 with respect to **p**.

Farra & Pšenčík (2003) use the vectors **e**_{[k]} to define the symmetric 3 × 3 matrix **B**(**p**) whose elements are

$$B_{kl} = \mathbf{e}_{[k]}^{\mathrm{T}} \mathbf{\Gamma} \mathbf{e}_{[l]}, \quad k = 1, 2, 3, \quad l = 1, 2, 3. \quad (17)$$

The elements of $\mathbf{\Gamma}$ and the components of the vectors $\mathbf{e}_{[k]}$ being homogeneous functions of degree 2 and 0, respectively, with respect to **p**, it is easy to show that the elements $B_{ij}$ are homogeneous functions of degree 2.

In isotropic media, the matrix **B**(**p**) is diagonal with diagonal terms $B_{11}(\mathbf{p}) = B_{22}(\mathbf{p}) = \beta^2 p^2$ and $B_{33}(\mathbf{p}) = \alpha^2 p^2$, where $\alpha$ and $\beta$ denote P- and S-wave velocities of the isotropic medium. In weakly anisotropic media, the first-order approximation of the phase velocity depends on elements of the matrix **B** (Farra & Pšenčík 2003). The matrix **B** is independent of the choice of the reference isotropic medium. It can be shown that the quantity $B_{11} - B_{22}$ and the off-diagonal elements of the matrix **B** are of first order in the parameter $\epsilon$ (Farra & Pšenčík 2003).

The first-order approximations $G_{[m]}^{(1)}(\mathbf{p})$ of the eigenvalues of the Christoffel matrix can be written in terms of the elements $B_{ij}(\mathbf{p})$, see Farra (2001):

$$G_{[1]}^{(1)} = \frac{1}{2} \left( B_{11} + B_{22} + \sqrt{(B_{11} - B_{22})^2 + 4B_{12}^2} \right), \quad G_{[2]}^{(1)} = \frac{1}{2} \left( B_{11} + B_{22} - \sqrt{(B_{11} - B_{22})^2 + 4B_{12}^2} \right), \quad G_{[3]}^{(1)} = B_{33}. \quad (18)$$

The approximate formulae (18) are valid to first order in $\epsilon$. One can notice that $G_{[1]}^{(1)}(\mathbf{p}) \ge G_{[2]}^{(1)}(\mathbf{p})$.

Since the elements of the matrix **B** are independent of the velocities $\alpha$ and $\beta$ of the reference isotropic media, $G_{[m]}^{(1)}$ in (18) are independent of them, too. As the elements $B_{ij}$, $G_{[m]}^{(1)}$ are homogeneous functions of degree 2 with respect to **p**.

In the wave normal direction defined by the unit vector **n**, the first-order approximation of the phase velocity squared, denoted $V_{[m]}^{(1)2}$, is given by

$$V_{[m]}^{(1)2}(\mathbf{n}) = G_{[m]}^{(1)}(\mathbf{n}), \quad m = 1, 2, 3. \quad (19)$$

In the first-order approximation, the singularity directions are characterized by coincident first-order phase velocities, so that $G_{[1]}^{(1)} = G_{[2]}^{(1)}$, and therefore,

$$B_{12} = 0, \quad B_{11} = B_{22} \quad (20)$$

(see eq. 18). Let us remark that the first-order singularity directions may deviate from the exact singularity directions (see Farra & Pšenčík 2003).

In the following, we frequently use two specifications of the vectors $\mathbf{e}_{[k]}$ and, thus, of the matrix **B**. In one, we use the vectors $\mathbf{e}_{[k]}$, denoted $\hat{\mathbf{e}}_{[k]}$, with $\hat{\mathbf{e}}_{[3]} = \mathbf{n}$ and $\hat{\mathbf{e}}_{[K]}$ chosen so that the corresponding **B** matrix, denoted $\hat{\mathbf{B}}$, satisfies the conditions

$$\hat{B}_{12} = 0, \quad \hat{B}_{11} > \hat{B}_{22}. \quad (21)$$

The second condition excludes singularity directions from our consideration (see eq. 20). For singularity directions, the vectors $\hat{\mathbf{e}}_{[1]}$ and $\hat{\mathbf{e}}_{[2]}$ cannot be specified uniquely.

The diagonal elements of the matrix $\hat{\mathbf{B}}(\mathbf{p})$ specify the first-order approximations $G_{[m]}^{(1)}(\mathbf{p})$ of the eigenvalues of the Christoffel matrix (Farra 2001) (see eqs 18 and 21):

$$G_{[m]}^{(1)} = \hat{B}_{mm} = \hat{\mathbf{e}}_{[m]}^{\mathrm{T}} \mathbf{\Gamma} \hat{\mathbf{e}}_{[m]}, \quad m = 1, 2, 3. \quad (22)$$

The vectors $\hat{\mathbf{e}}_{[1]}$ and $\hat{\mathbf{e}}_{[2]}$ specify the zero-order polarization vectors $\mathbf{g}_{[1]}^{(0)}$ and $\mathbf{g}_{[2]}^{(0)}$ of the $qS$ waves.

In the other specification we use the vectors $\mathbf{e}_{[k]}$, denoted $\bar{\mathbf{e}}_{[k]}$ and defined in Appendix A. We denote the corresponding matrix by a bar: $\bar{\mathbf{B}}$. For orthorhombic media, explicit expressions for the elements of the matrix $\bar{\mathbf{B}}(\mathbf{n})$ can be found in Appendix A. Though the equations written in the following sections are obtained for arbitrary vectors $\mathbf{e}_{[K]}$, $K = 1, 2$, we use the vectors $\bar{\mathbf{e}}_{[K]}$ and the matrix $\bar{\mathbf{B}}$ for practical applications.

Elements of a matrix **B** specified for arbitrarily chosen vectors $\mathbf{e}_{[k]}$ are related to elements of $\hat{\mathbf{B}}$ by simple relations. If we denote by $\Phi^{(0)}(\mathbf{p})$, $0 \le \Phi^{(0)} < \pi$, the angle between the vectors $\mathbf{e}_{[1]}(\mathbf{p})$ and $\hat{\mathbf{e}}_{[1]}(\mathbf{p})$, the vectors $\hat{\mathbf{e}}_{[1]}$ and $\hat{\mathbf{e}}_{[2]}$ can be expressed in terms of unit vectors $\mathbf{e}_{[K]}$ in the following way:

$$\hat{\mathbf{e}}_{[1]} = \mathbf{e}_{[1]} \cos \Phi^{(0)} + \mathbf{e}_{[2]} \sin \Phi^{(0)}, \quad \hat{\mathbf{e}}_{[2]} = -\mathbf{e}_{[1]} \sin \Phi^{(0)} + \mathbf{e}_{[2]} \cos \Phi^{(0)}. \quad (23)$$
The elements $B_{IJ}$ transform into $\hat{B}_{IJ}$ in the following way (see eqs 17 and 23):

$$ 
\begin{aligned}
\hat{B}_{11} &= B_{11} \cos^2 \Phi^{(0)} + 2B_{12} \cos \Phi^{(0)} \sin \Phi^{(0)} + B_{22} \sin^2 \Phi^{(0)}, \\
\hat{B}_{22} &= B_{11} \sin^2 \Phi^{(0)} - 2B_{12} \cos \Phi^{(0)} \sin \Phi^{(0)} + B_{22} \cos^2 \Phi^{(0)}, \\
\hat{B}_{12} &= (B_{22} - B_{11}) \cos \Phi^{(0)} \sin \Phi^{(0)} + B_{12}(\cos^2 \Phi^{(0)} - \sin^2 \Phi^{(0)}).
\end{aligned}
\quad (24) $$

The angle $\Phi^{(0)}$ can be determined from the equation:

$$ \tan 2\Phi^{(0)} = \frac{2B_{12}}{B_{11} - B_{22}}. \quad (25) $$

Eq. (25) follows from the first condition in (21), taking into account the relation for $\hat{B}_{12}$ in (24). Two solutions $\Phi^{(0)}, 0 \le \Phi^{(0)} < \pi$, differing by $\pi/2$, satisfy (25). The second condition in (21) guarantees unique determination of the angle $\Phi^{(0)}$.

The components of the zero-order polarization vectors $g_{[M]}^{(0)} = \hat{e}_{[M]}$ in the basis $e_{[k]}$ are denoted as follows (see eq. 23):

$$ g_{[1]}^{(0)} = \begin{bmatrix} \cos \Phi^{(0)} \\ \sin \Phi^{(0)} \\ 0 \end{bmatrix}, \quad g_{[2]}^{(0)} = \begin{bmatrix} -\sin \Phi^{(0)} \\ \cos \Phi^{(0)} \\ 0 \end{bmatrix}. \quad (26) $$

Thus, from (22) and (24), one can write for $M = 1, 2$:

$$ G_{[M]}^{(1)} = g_{[M]}^{(0)T} \mathbf{B} g_{[M]}^{(0)}. \quad (27) $$

Moreover, the following relation can be deduced from the first condition in (21), taking into account the relation for $\hat{B}_{12}$ in (24):

$$ g_{[1]}^{(0)T} \mathbf{B} g_{[2]}^{(0)} = \hat{B}_{12} = 0. \quad (28) $$

## 5 FIRST-ORDER RAY TRACING EQUATIONS

In the following, we derive an approximate, first-order ray tracing (FORT) for the $qS$ waves propagating in a smooth inhomogeneous weakly anisotropic medium. The matrix $\mathbf{B}$, the approximate eigenvalues $G_{[M]}^{(1)}$, the angle $\Phi^{(0)}$ and the vectors $g_{[M]}^{(0)}$ are functions of the position vector $\mathbf{x}$.

Instead of the exact eigenvalue $G$, we use the approximate expression (18) and insert it into the Hamiltonian (4):

$$ H^{(1)}(\mathbf{x}, \mathbf{p}) = \frac{1}{2} (G^{(1)}(\mathbf{x}, \mathbf{p}) - 1). \quad (29) $$

The approximate rays $(\mathbf{x}^{(1)}, \tau)$, $(\mathbf{p}^{(1)}, \tau))$ are solutions of the approximate ray tracing equations:

$$ \frac{dx_i^{(1)}}{d\tau} = \frac{\partial H^{(1)}}{\partial p_i} = \frac{1}{2} \frac{\partial G^{(1)}}{\partial p_i}, \quad \frac{dp_i^{(1)}}{d\tau} = -\frac{\partial H^{(1)}}{\partial x_i} = -\frac{1}{2} \frac{\partial G^{(1)}}{\partial x_i}. \quad (30) $$

Let us mention that the first set of eqs (30) corresponds to the first-order group velocity (Farra 2004).

At each point of the ray, the approximate slowness vector satisfies the approximate eikonal equation:

$$ H^{(1)}(\mathbf{x}^{(1)}, \mathbf{p}^{(1)}) = 0. \quad (31) $$

As $H^{(1)}$ is constant along any solution of (30), i.e.

$$ \frac{dH^{(1)}}{d\tau} = \frac{\partial H^{(1)}}{\partial x_i} \frac{dx_i^{(1)}}{d\tau} + \frac{\partial H^{(1)}}{\partial p_i} \frac{dp_i^{(1)}}{d\tau} = 0, \quad (32) $$

it is sufficient to satisfy the eikonal equation (31) at some initial value $\tau_0$ in order to satisfy it along the whole ray.

The application of expression (30) for the $qS$ waves is more complicated than for the $qP$ wave because of singularities. Using the expression (27) for $G_{[M]}^{(1)}$, $M = 1, 2$, one can write (30) as:

$$ \frac{dx_i^{(1)}}{d\tau} = \frac{1}{2} \mathbf{g}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_i} \mathbf{g}^{(0)}, \quad (33) $$

$$ \frac{dp_i^{(1)}}{d\tau} = -\frac{1}{2} \mathbf{g}^{(0)T} \frac{\partial \mathbf{B}}{\partial x_i} \mathbf{g}^{(0)}. \quad (34) $$

In order to obtain (33) and (34), we use the property (28) and the following relations deduced from (26):

$$ \frac{\partial \mathbf{g}_{[1]}^{(0)}}{\partial z_i} = \mathbf{g}_{[2]}^{(0)} \frac{\partial \Phi^{(0)}}{\partial z_i}, \quad \frac{\partial \mathbf{g}_{[2]}^{(0)}}{\partial z_i} = -\mathbf{g}_{[1]}^{(0)} \frac{\partial \Phi^{(0)}}{\partial z_i}, \quad (35) $$

where $z_i$ may be either $x_i$ or $p_i$, so that

$$ \mathbf{g}_{[M]}^{(0)T} \mathbf{B} \frac{\partial \mathbf{g}_{[M]}^{(0)}}{\partial z_i} = 0. \quad (36) $$
In (33) and (34), the partial derivatives of the elements of matrix **B** are calculated at $(\mathbf{x}^{(1)}, \mathbf{p}^{(1)})$. The partial derivatives $\partial B_{jk}/\partial p_i$ can be expressed in terms of the spherical coordinates $(\rho, \theta, \phi)$ of vector **p** (see Appendix B). The partial derivatives $\partial B_{jk}/\partial x_i$ can be determined from the expression of $B_{jk}(\mathbf{x}, \mathbf{p})$, keeping in mind that the only spatially dependent parameters are the elastic parameters.

Eqs (33) and (34) can be compared with the exact ray tracing equations (7). The Christoffel matrix $\Gamma$ and the exact polarization vector $\mathbf{g}_{[M]}$ of the $qS_M$ wave are substituted by the matrix **B** and the zero-order polarization vector $\mathbf{g}_{[M]}^{(0)}$ whose components are expressed in the basis $\mathbf{e}_{[k]}$. Noticing that $g_{[M]}^{(0)} = 0$ (see eq. 26), one can see that the first-order ray tracing equations (33) and (34) for $qS$ waves are related to the three elements $B_{11}$, $B_{12}$ and $B_{22}$ of matrix **B** and their partial derivatives with respect to $x_i$ and $p_i$. Therefore, they are only controlled by 15 combinations of density-normalized elastic parameters (see Farra & Pšenčík 2003). These combinations are the following: $A_{44}$, $A_{55}$, $A_{66}$, $2A_{12} - A_{11} - A_{22}$, $2A_{13} - A_{11} - A_{33}$, $2A_{23} - A_{22} - A_{33}$, $A_{14} - A_{24}$, $A_{14} - A_{34}$, $A_{15} - A_{25}$, $A_{15} - A_{35}$, $A_{16} - A_{26}$, $A_{16} - A_{36}$, $A_{45}$, $A_{46}$, $A_{56}$. Note that they are different from combinations for the $qP$ wave. The integration of the FORT equations should be faster, and perhaps with fewer numerical errors, than the integration of the exact equations which depend on the 21 elastic parameters.

It may happen when integrating eqs (33) and (34) that the approximate ray crosses a singularity, i.e. the direction of the local slowness vector $\mathbf{p}^{(1)}(\tau)$ is a first-order singularity direction for the medium at $\mathbf{x}^{(1)}(\tau)$. The eigenvalue should be changed ($G_{[1]}^{(1)}$ by $G_{[2]}^{(1)}$ or $G_{[2]}^{(1)}$ by $G_{[1]}^{(1)}$) when the ray crosses the conical or the intersection singularity so that the zero-order polarization vector $\mathbf{g}^{(0)}$ and therefore the first-order group velocity change smoothly along the ray (see Appendix C). This change is required otherwise the ray tracing can produce an unphysical abrupt change of the ray direction when crossing the singularity (see Vavryčuk 2001). Therefore, the wave mode can be $qS_1$ on some parts of the ray and $qS_2$ on some other parts. The presence of the first-order singularity can be detected by looking at the behaviour of $B_{11} - B_{22}$ and $B_{12}$ along the ray (see eq. 20), the type of the singularity being determined from the partial derivatives of the $B_{JK}$ (see Appendix C).

# 6 FIRST-ORDER PARAXIAL RAY TRACING SYSTEM

In the first-order approximation, we can compute the propagator matrix $\mathbf{P}^{(1)}$ of the approximate paraxial ray tracing system calculated along the approximate ray $(\mathbf{x}^{(1)}(\tau), \mathbf{p}^{(1)}(\tau))$:

$$ \frac{d\mathbf{P}^{(1)}}{d\tau} = \mathbf{D}^{(1)}(\tau)\mathbf{P}^{(1)}, \qquad (37) $$

where $\mathbf{D}^{(1)}$ is given by the expressions (12) and (13) with the Hamiltonian $H$ substituted by $H^{(1)}$ (see eq. 29).

The elements of the matrix $\mathbf{D}^{(1)}$ are simply related through expressions similar to (13) to the second partial derivatives of $G^{(1)}$. From expressions (27), (28) and (35), one can write the second partial derivatives of $G_{[M]}^{(1)}$ as:

$$ \frac{\partial^2 G_{[M]}^{(1)}}{\partial z_i \partial z_j} = \mathbf{g}_{[M]}^{(0)T} \frac{\partial^2 \mathbf{B}}{\partial z_i \partial z_j} \mathbf{g}_{[M]}^{(0)} + \sum_{M'=1}^{M'=2, M' \neq M} \frac{2}{G_{[M]}^{(1)} - G_{[M']}} \left( \mathbf{g}_{[M]}^{(0)T} \frac{\partial \mathbf{B}}{\partial z_i} \mathbf{g}_{[M']}(0) \right) \left( \mathbf{g}_{[M]}^{(0)T} \frac{\partial \mathbf{B}}{\partial z_j} \mathbf{g}_{[M']}(0) \right), \qquad (38) $$

where $z_i$ may be either $x_i$ or $p_i$. In order to obtain (38), we use the following relation:

$$ \left( G_{[1]}^{(1)} - G_{[2]}^{(1)} \right) \frac{\partial \Phi^{(0)}}{\partial z_i} = \mathbf{g}_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial z_i} \mathbf{g}_{[2]}^{(0)}, \quad i = 1, 2, 3, \qquad (39) $$

derived from (28). In (38), all the quantities are calculated at $(\mathbf{x}^{(1)}(\tau), \mathbf{p}^{(1)}(\tau))$.

Expression (38) can be compared with the second partial derivatives of the exact eigenvalue (see eq. 14). The Christoffel matrix $\Gamma$ and the exact polarization vectors $\mathbf{g}_{[m]}$ of the three wave modes are substituted by the matrix **B** and the zero-order polarization vectors of the two $qS$ waves whose components are expressed in the basis $\mathbf{e}_{[k]}$. This should reduce the number of operations necessary for the evaluation of the right-hand side of the paraxial ray tracing system.

The right-hand side of (38) is singular wherever $G_{[1]}^{(1)} = G_{[2]}^{(1)}$, i.e. when the approximate ray $(\mathbf{x}^{(1)}, \mathbf{p}^{(1)})$ crosses a first-order singularity. The behaviour of the second partial derivatives of $G_{[M]}^{(1)}$ with respect to $\mathbf{p}$ components can be studied in the vicinity of singularities (see Appendix C). For the kiss and intersection singularities, the last term in (38) has a finite limit when the wave normal direction approaches the singularity direction. For the conical singularity, this term tends to infinity. The other second partial derivatives of $G_{[M]}^{(1)}$ can tend to infinity at singularities depending on the evolution of the singularity direction with position. The first-order equations have the same behaviour as the exact equations in the vicinity of singularities. The first-order formulae being explicit, their behaviour can be studied analytically.

Noticing that $g_{[13]}^{(0)} = g_{[23]}^{(0)} = 0$ (see eq. 26), one can see that, as for the first-order ray tracing equations (33) and (34), the first-order paraxial ray tracing system (37) is related to the three elements $B_{11}$, $B_{12}$ and $B_{22}$ of matrix **B** and their partial derivatives with respect to $x_i$ and $p_i$. Therefore, it only depends on 15 combinations of density-normalized elastic parameters.

# 7 SECOND-ORDER TRAVELTIME

Let us introduce $\Delta H = H - H^{(1)}$ the difference between the exact Hamiltonian $H$ and its approximation (29). Farra & Le Bégat (1995) show that an approximate traveltime can be obtained by integration of the term $\Delta H$ along an approximate trajectory. Using the FORT ray $(\mathbf{x}^{(1)}, \mathbf{p}^{(1)})$ as the approximate trajectory, one can write the approximate traveltime, denoted $T^{(2)}$, between $\mathbf{x}^{(1)}(\tau_0)$ and $\mathbf{x}^{(1)}(\tau)$ as:

$$ T^{(2)} = T^{(1)} - \int_{\tau_0}^{\tau} \Delta H(\mathbf{x}^{(1)}, \mathbf{p}^{(1)}) d\tau, \qquad (40) $$
where $T^{(1)}$ is the traveltime along the FORT ray:

$$T^{(1)} = \int_{\tau_0}^{\tau} \left( \mathbf{p}^{(1)} \cdot \frac{d\mathbf{x}^{(1)}}{d\tau} - H^{(1)}(\mathbf{x}^{(1)}, \mathbf{p}^{(1)}) \right) d\tau. \quad (41)$$

The FORT ray being the first-order approximation of a ray, the traveltime (41) is of the first order and the approximation (40) is of the second order.

The FORT ray satisfies the approximate eikonal (31); thus the second term in (41) is zero. Moreover, the eigenvalue $G^{(1)}$ being a homogeneous function of degree 2 in $\mathbf{p}$, it can be shown that:

$$\mathbf{p}^{(1)} \cdot \frac{d\mathbf{x}^{(1)}}{d\tau} = \frac{1}{2} p_i^{(1)} \frac{\partial G^{(1)}}{\partial p_i} = 1, \quad (42)$$

so that expression (41) can be simply written as:

$$T^{(1)} = \tau - \tau_0. \quad (43)$$

Therefore, $T^{(1)}$ is obtained by solving the FORT equations.

An estimation of $\Delta H = \frac{1}{2}(G - G^{(1)})$ can be obtained from the second-order approximation $G^{(2)}$ of the eigenvalues of the qS-waves, i.e. $\Delta H = \frac{1}{2}(G^{(2)} - G^{(1)})$ in (40). The second-order approximation of the eigenvalues of the qS waves can be written as (Farra 2001; Farra & Pšenčík 2003):

$$G^{(2)} = \mathbf{g}^{(1)\mathrm{T}} \mathbf{M} \mathbf{g}^{(1)}. \quad (44)$$

The elements of the matrix $\mathbf{M}$ specified for the basis $\mathbf{e}_{[k]}$ are given by:

$$M_{I,J}(\mathbf{x}, \mathbf{p}) = B_{I,J}(\mathbf{x}, \mathbf{p}) + \frac{B_{I3}(\mathbf{x}, \mathbf{p})B_{J3}(\mathbf{x}, \mathbf{p})}{(\beta^2(\mathbf{x}) - \alpha^2(\mathbf{x}))p^2}, \quad (45)$$

where $\alpha$ and $\beta$ are the $P$ and $S$ velocities of the reference isotropic medium. The $P$ and $S$ velocities of the so-called isotropic replacement medium (IRM) (see e.g. Mensch & Rasolofosaon 1997) can be used for $\alpha$ and $\beta$.

In (44), the vector $\mathbf{g}^{(1)}$ is specified by its components in the basis $\mathbf{e}_{[k]}$:

$$\mathbf{g}^{(1)} = \begin{bmatrix} \cos \Phi^{(1)} \\ \sin \Phi^{(1)} \\ 0 \end{bmatrix}, \quad (46)$$

where the angle $\Phi^{(1)}$, $0 \le \Phi^{(1)} < \pi$, is determined from the equation:

$$\tan 2\Phi^{(1)} = \frac{2M_{12}}{M_{11} - M_{22}}. \quad (47)$$

Let us remark that the second-order approximation (44) of the eigenvalues of the qS waves is similar to the first-order approximation (27) with matrix $\mathbf{B}$ substituted by matrix $\mathbf{M}$ and vector $\mathbf{g}^{(0)}$ by vector $\mathbf{g}^{(1)}$. Each of the two solutions $\Phi^{(1)}$, $0 \le \Phi^{(1)} < \pi$, of eq. (47) is associated with one of the qS waves, the vector $\mathbf{g}^{(1)}$ being the projection of the first-order polarization vector in the plane perpendicular to $\mathbf{n}$ (see Farra & Pšenčík 2003). Among the two solutions of eq. (47), the angle $\Phi^{(1)}$ should be chosen so that the vector $\mathbf{g}^{(1)}$ used in (44) is the closest to $\mathbf{g}^{(0)}$. This guarantees that we are dealing with first- and second-order wave-surface elements which are close to each other. This criterion is important in the neighbourhood of singularities.

## 8 FORT FOR VTI MEDIA

We now specialize to a transversely isotropic medium with symmetry axis along the $x_3$-axis (VTI medium). The phase velocity surface of the qS waves can be separated in two regular sheets in terms of polarizations: the $qSV$ wave is associated with polarization vector $\mathbf{g}_{[SV]}$ (polarization vector within the incident vertical plane) and the $qSH$ wave with polarization vector $\mathbf{g}_{[SH]} = \bar{\mathbf{e}}_{[2]}$ (polarization vector normal to the incident vertical plane). The separation in terms of polarizations is impossible for general anisotropy.

In a VTI medium, the elements $\bar{\mathbf{B}}_{ij}$ are given by:

$$\begin{align}
\bar{\mathbf{B}}_{11}(\mathbf{x}, \mathbf{p}) &= A_{44} p^2 - 2 \hat{A}_{13} \frac{(p_1^2 + p_2^2) p_3^2}{p^2}, \nonumber \\
\bar{\mathbf{B}}_{22}(\mathbf{x}, \mathbf{p}) &= A_{66} (p_1^2 + p_2^2) + A_{44} p_3^2, \nonumber \\
\bar{\mathbf{B}}_{13}(\mathbf{x}, \mathbf{p}) &= p_3 \sqrt{p_1^2 + p_2^2} \left( \hat{A}_{13} \frac{p_3^2 - (p_1^2 + p_2^2)}{p^2} + \frac{1}{2}(A_{11} - A_{33}) \right), \nonumber \\
\bar{\mathbf{B}}_{12}(\mathbf{x}, \mathbf{p}) &= \bar{\mathbf{B}}_{23}(\mathbf{x}, \mathbf{p}) = 0 \tag{48}
\end{align}$$

(see Appendix A).
The angles $\Phi^{(0)}$ and $\Phi^{(1)}$ are equal to 0° or 90° (see eqs 25 and 47); it yields $\tilde{e}_{[1]}$ for the $qSV$-wave zero-order polarization vector and $\tilde{e}_{[2]}$ for the $qSH$ wave. The first-order phase velocities squared of the two $qS$ waves are $V_{[SV]}^{(1)2}(\mathbf{x}, \mathbf{n}) = \tilde{B}_{11}(\mathbf{x}, \mathbf{n})$ and $V_{[SH]}^{(1)2}(\mathbf{x}, \mathbf{n}) = \tilde{B}_{22}(\mathbf{x}, \mathbf{n})$, respectively.

The first-order Hamiltonians (29) can be written as:

$$ 
\begin{aligned}
H_{[SV]}^{(1)}(\mathbf{x}, \mathbf{p}) &= \frac{1}{2} \left( A_{44} p^2 - 2 \hat{A}_{13} \frac{(p_1^2 + p_2^2) p_3^2}{p^2} - 1 \right), \\
H_{[SH]}^{(1)}(\mathbf{x}, \mathbf{p}) &= \frac{1}{2} \left[ A_{66} (p_1^2 + p_2^2) + A_{44} p_3^2 - 1 \right].
\end{aligned} 
\quad (49) $$

The FORT equations can be obtained from (30). For the $qSH$ wave, the polarization vector in the zero-order approximation coincides with the exact one. Thus, the FORT equations are exact and can be written as

$$ 
\begin{aligned}
\frac{dx_1}{d\tau} &= A_{66}p_1, \\
\frac{dx_2}{d\tau} &= A_{66}p_2, \\
\frac{dx_3}{d\tau} &= A_{44}p_3, \\
\frac{dp_i}{d\tau} &= -\frac{1}{2} \left( \frac{\partial A_{66}}{\partial x_i} (p_1^2 + p_2^2) + \frac{\partial A_{44}}{\partial x_i} p_3^2 \right).
\end{aligned} 
\quad (50) $$

Moreover, as $\Delta H_{[SH]} = 0$, the traveltime along the ray can be determined directly from (43).

For the $qSV$ wave, the FORT equations can be written as

$$ 
\begin{aligned}
\frac{dx_1}{d\tau} &= \left( A_{44} - 2\hat{A}_{13} \frac{p_3^4}{p^4} \right) p_1, \\
\frac{dx_2}{d\tau} &= \left( A_{44} - 2\hat{A}_{13} \frac{p_3^4}{p^4} \right) p_2, \\
\frac{dx_3}{d\tau} &= \left( A_{44} - 2\hat{A}_{13} \frac{(p_1^2 + p_2^2)^2}{p^4} \right) p_3, \\
\frac{dp_i}{d\tau} &= -\frac{1}{2} \left( \frac{\partial A_{44}}{\partial x_i} p^2 - 2 \frac{\partial \hat{A}_{13}}{\partial x_i} \frac{(p_1^2 + p_2^2) p_3^2}{p^2} \right).
\end{aligned} 
\quad (51) $$

The traveltime can be obtained to second order from (40) with the Hamiltonian perturbation given by

$$ \Delta H_{[SV]}(\mathbf{x}, \mathbf{p}) = \frac{1}{2} \frac{p_3^2 (p_1^2 + p_2^2)}{(\beta^2 - \alpha^2) p^2} \left( \hat{A}_{13} \frac{p_3^2 - (p_1^2 + p_2^2)}{p^2} + \frac{1}{2} (A_{11} - A_{33}) \right)^2. \quad (52) $$

In the isotropic case, the phase velocities of the two shear-waves are equal and do not depend on the wave normal direction:

$$ V_{[1]}^2 = V_{[2]}^2 = A_{44}. \quad (53) $$

The first-order Hamiltonian is given by $H(\mathbf{x}, \mathbf{p}) = \frac{1}{2}(A_{44}p^2 - 1)$ and the FORT equations reduce to:

$$ 
\begin{aligned}
\frac{dx_i}{d\tau} &= A_{44}p_i, \\
\frac{dp_i}{d\tau} &= -\frac{p^2}{2} \frac{\partial A_{44}}{\partial x_i}.
\end{aligned} 
\quad (54) $$

The above equations are standard ray tracing equations for isotropic media with S-wave velocity $\sqrt{A_{44}}$.

## 9 NUMERICAL EXAMPLES

In order to illustrate the performance of the second-order traveltime formula we consider the simple case of a homogeneous transversely isotropic medium with symmetry axis along the $x_3$-axis (VTI medium). The model is characterized by the parameters $A_{ij}$, in $(\text{km/s})^2$, with values $A_{11} = A_{22} = 20.16$, $A_{33} = 19.63$, $A_{12} = 7.40$, $A_{13} = A_{23} = 7.26$, $A_{44} = A_{55} = 3.48$, $A_{66} = 6.38$. This model corresponds to model B in Pšenčík & Vavryčuk (2002) and Farra (2004). Anisotropy of this model is about 30 per cent for $qS$ waves. Due to the axial symmetry of the medium, it is sufficient to investigate just a quadrant of a vertical plane containing the axis of symmetry. Fig. 1 shows the phase velocities of the two $qS$ waves as a function of angle $\theta$. The angle $\theta$ specifies the direction of the wave normal, $\theta = 0^{\circ}$ corresponds to the direction along the axis of symmetry and $\theta = 90^{\circ}$ corresponds to the direction perpendicular to it. The phase velocity is shown by circles for the $qS_1$ wave and by crosses for the $qS_2$ wave. The $qS_1$ wave is $qSV$ for $\theta < 59^{\circ}$ and $qSH$ for $\theta > 59^{\circ}$.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Phase velocity section of *qS* waves for the TI medium specified in the text. *θ* is the angle of the wave normal with the axis of symmetry. The *qS*₁ wave is shown by circles, the *qS*₂ wave by crosses.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Normalized traveltimes of *qS* waves as functions of the ray angle with the symmetry axis. The normalized traveltime residual Δ*T*/(*T* − *T*ᵢ)/*T*ᵢ is computed with respect to the traveltime *T*ᵢ obtained in the isotropic medium with *S* velocity *V*<sub>S</sub> = 1.87 km s<sup>−1</sup>. *T*<sub>e</sub> (black lines) is the exact traveltime of the *qSV* and *qSH* waves in the homogeneous TI medium specified in the text. The residuals corresponding to first- and second-order traveltimes *T*₁ and *T*₂ are shown by blue and red crosses, respectively.</figcaption></figure>

In homogeneous anisotropic media, dpᵢ⁽¹⁾/dτ = 0 in (34) and thus the solution of the first-order ray tracing equations can be obtained analytically. Let us introduce the normalized traveltime residual Δ*T* = (*T* − *T*ᵢ)/*T*ᵢ with respect to the traveltime *T*ᵢ calculated in the isotropic medium with *S* velocity *V*<sub>S</sub> = 1.87 km s<sup>−1</sup>. We denote by *T*<sub>e</sub> the exact traveltime and *T*₁ and *T*₂ the first- and second-order approximations of the traveltime (see eqs 43 and 40). The traveltimes are calculated for the two *qS* waves, the exact and approximate rays being computed for regularly specified wave-normal angles *θ*. In expression (40) for the second-order traveltime, the *P* and *S* velocities of the IRM medium were used. In Fig. 2, the three quantities (*T*<sub>e</sub> − *T*ᵢ)/*T*ᵢ, (*T*₁ − *T*ᵢ)/*T*ᵢ and (*T*₂ − *T*ᵢ)/*T*ᵢ are plotted as functions of the ray angle with the symmetry axis. In homogeneous media, these quantities are independent of the distance between the source and the station. The variability of (*T*<sub>e</sub> − *T*ᵢ)/*T*ᵢ with the ray direction is only due to anisotropy. One can see the triplication of the *qSV* wave which is due to the anisotropy. The second-order traveltime approximates very well the exact traveltime except in the close vicinity of the cusps. There is a mispositioning of the cusps due to the first-order approximation of the ray direction. One can notice that the first-order and second-order traveltimes corresponding to the *qSH* wave are identical to the exact traveltimes.

Let us now consider a more complicated anisotropy, the orthorhombic medium used by Crampin (1991). The model is specified by the density-normalized elastic parameters *A*<sub>ij</sub>, in (km/s)<sup>2</sup>, with values *A*<sub>11</sub> = 16.26, *A*<sub>22</sub> = 16.36, *A*<sub>33</sub> = 13.61, *A*<sub>12</sub> = 4.41, *A*<sub>13</sub> = 3.61, *A*<sub>23</sub> = 3.63, *A*<sub>44</sub> = 5.15, *A*<sub>55</sub> = 4.10, *A*<sub>66</sub> = 4.73. Fig. 3 shows the map of relative differences (in per cent) of *qS*₁ and *qS*₂ phase velocities as function of polar angle *θ* and azimuth *φ* specifying the wave normal. Because the orthorhombic symmetry has three mutually perpendicular planes of symmetry, all the calculations are made in the octant defined by 0° ≤ *θ* ≤ 90° and 0° ≤ *φ* ≤ 90°. This model has four singularities in the map shown: two singularities for *φ* = 0°, *θ* = 26° and *θ* = 66°, one at *θ* = 90° and *φ* = 73°, and another one for *θ* = 52° and *φ* = 52°. The singularities often appear with triplications of the wave front. The surface of the *qS*₁ wave has a triplication in the vicinity of the *Ox*₁*x*₃ propagation plane (*φ* = 0°) for the ray direction polar angle around 45°. The *qS*₂-wave surface has two triplications: the one crossing most
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Map of relative difference (in per cent) of $qS_1$ and $qS_2$ wave phase velocities as function of azimuth $\phi$ and polar angle $\theta$ of the wave normal for the orthorhombic medium specified in the text.</figcaption></figure>

of the octant is related to the pinch connecting three of the singularities in Fig. 3, the other triplication is related to the singularity at $\theta = 66^\circ$ and $\phi = 0^\circ$ (see Fig. 4 in Crampin 1991).

In Fig. 4, the exact and second-order quantities $(T_e - T_i)/T_i$ and $(T_2 - T_i)/T_i$ are plotted as functions of the ray angle $\theta_g$ with the $x_3$-axis for three wave normal azimuths (0°, 45° and 90°). The normalized traveltime residuals are computed with respect to the traveltime $T_i$ calculated in the isotropic medium with $S$ velocity $V_S = 2.27 \text{ km s}^{-1}$. The traveltimes were calculated for the two $qS$ waves, the exact and approximate rays being computed for regularly specified wave normal angles $\theta$. In expression (40) for the second-order traveltime, the $P$ and $S$ velocities of the IRM medium were used. Let us remark that for the wave-normal azimuth $\phi = 45^\circ$, the rays are not contained in the same vertical plane. The loop seen on the corresponding panel is related to one of the $qS_2$ wave triplications. The second-order traveltime approximates very well the exact traveltime except in the close vicinity of the $qS_1$ cusp in the $0^\circ$ azimuth.

## 10 CONCLUSION

The first-order ray tracing equations and paraxial ray tracing system are obtained for $qS$ waves propagating in inhomogeneous weakly anisotropic media. The equations are written using the Hamiltonian formulation with the first-order Hamiltonian substituting the exact one. A second-order approximation of the traveltime is obtained by integration of a correcting term along the approximate rays. In contrast to other methods used to compute approximate rays (Farra 1989; Nowack & Pšenčík 1991; Mensch & Farra 1999), the proposed approach does not require the calculation of reference rays in a reference medium. The FORT rays are obtained directly by solving the FORT equations. A reference medium is only required for the calculation of the second-order traveltime.

In the perturbation formulae, an important role is played by the matrix **B** whose elements control various attributes of elastic waves (see Farra & Pšenčík 2003). The first-order ray equations for $qS$ waves depend on three elements $B_{11}$, $B_{12}$ and $B_{22}$ of the matrix **B** and their partial derivatives. They are based on explicit formulae which make the dependence on parameters transparent. The equations are only controlled by 15 independent combinations of the elastic parameters. This reduces number of operations necessary for their evaluation. For isotropic media, the first-order equations reduce to standard exact ray equations.

The first-order formulae being explicit, their behaviour can be studied analytically in the vicinity of shear wave singularities. The presence of a singularity along the ray can be detected from the quantities $B_{11} - B_{22}$ and $B_{12}$. Ray tracing equations behave properly if the wave mode is changed ($qS_1$ into $qS_2$ or vice versa) when the ray crosses a conical or intersection singularity. However, the paraxial ray tracing system may be singular at singularities. Ray theory and its extensions, such as the Gaussian beam method (Červený et al. 1982) and Maslov asymptotic theory (Chapman & Drummond 1982), should meet difficulties for the computation of synthetic seismograms, as they need to compute the paraxial quantities.

As shown by tests made in homogeneous TI and orthorhombic media, the use of the FORT rays together with the second-order traveltime correction gives negligible errors except in the vicinity of cusps. For TI anisotropy of about 30 per cent, the relative errors of the second-order traveltimes are under 0.3 per cent.

The next step in this study is the generalization of the approach to layered media. This will allow calculation of synthetic seismograms in laterally varying layered weakly anisotropic media. For $qS$ waves, it will be then necessary to address the problem of the $qS$-wave coupling.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Normalized traveltimes for *qS* waves as functions of the ray angle θg with the *x*3-axis. Each panel corresponds to a constant wave normal azimuth. The normalized traveltime residual Δ*T*/Tᵢ = (*T* − Tᵢ)/Tᵢ is computed with respect to the traveltime Tᵢ obtained in the isotropic medium with *S* velocity *V*<sub>S</sub> = 2.27 km s<sup>-1</sup>. T<sub>e</sub> (black crosses) is the exact traveltime of the *qS* waves in the homogeneous orthorhombic medium specified in the text. The residuals corresponding to second-order traveltimes T<sub>2</sub> are shown by red circles.</figcaption></figure>

ACKNOWLEDGMENTS

The author is grateful to Ivan Pšenčík for his constructive comments on the original manuscript. This is IPGP contribution no 2027.

REFERENCES

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

Červený, V., 1989. Ray tracing in factorized anisotropic inhomogeneous media, *Geophys. J. Int.*, **99**, 91–100.

Červený, V., 2001. *Seismic Ray Theory*, Cambridge University Press, Cambridge.

Červený, V., Popov, M.M. & Pšenčík, I., 1982. Computation of wave fields in inhomogeneous media—Gaussian beam approach, *Geophys. J. R. astr. Soc.*, **70**, 109–128.

Červený, V., Klimeš, L. & Pšenčík, I., 1988. Applications of dynamic ray tracing, *Phys. Earth planet. Inter.*, **51**, 25–35.
Chapman, C.H. & Drummond, R., 1982. Body-wave seismograms in inhomogeneous media using Maslov asymptotic theory, *Bull. seism. Soc. Am.*, **72**, S277-S317.

Chapman, C.H. & Shearer, P.M., 1989. Ray tracing in azimuthally anisotropic media—2. Quasi-shear wave coupling, *Geophys. J.*, **96**, 65–83.

Crampin, S., 1991. Effects of point singularities on shear-wave propagation in sedimentary basins, *Geophys. J. Int.*, **107**, 531–543.

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

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

Farra, V., 1993. Ray tracing in complex media, *J. appl. Geophys.*, **30**, 55–73.

Farra, V., 2001. High order perturbations of the phase velocity and polarization of $qP$ and $qS$ waves in anisotropic media, *Geophys. J. Int.*, **147**, 93–104.

Farra, V., 2004. Improved first-order approximation of group velocities in weakly anisotropic media, *Stud. Geoph. Geod.*, **48**, 199–213.

Farra, V. & Le Bégat, S., 1995. Sensitivity of $qP$-wave traveltimes and polarization vectors to heterogeneity, anisotropy and interfaces, *Geophys. J. Int.*, **121**, 371–384.

Farra, V. & Madariaga, R., 1987. Seismic waveform modeling in heterogeneous media by ray perturbation theory, *J. geophys. Res.*, **92**, 2697–2712.

Farra, V. & Pšenčík, I., 2003. Properties of the zero-, first- and higher-order approximations of attributes of elastic waves in weakly anisotropic media, *J. acoust. Soc. Am.*, **114**, 1366–1378.

Gajewski, D. & Pšenčík, I., 1987. Computation of high-frequency seismic wavefields in 3-D laterally inhomogeneous anisotropic media, *Geophys. J. R. astr. Soc.*, **91**, 383–411.

Gajewski, D. & Pšenčík, I., 1990. Vertical seismic profile synthetics by dynamic ray tracing in laterally varying layered anisotropic structures, *J. geophys. Res.*, **95**, 11 301–11 315.

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 $\epsilon$, $\delta$ and $\gamma$, *Geophys. J. Int.*, **128**, 43–64.

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. & Farra, V., 2004. First-order ray tracing for $qP$ waves in inhomogeneous weakly anisotropic media, *Geophysics*, submitted.

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.

Pšenčík, I. & Vavryčuk, V., 2002. Approximate relation between the ray vector and the wave normal in weakly anisotropic media, *Stud. Geophys. Geod.*, **46**, 793–807.

Shearer, P.M. & Chapman, C.H., 1989. Ray tracing in azimuthally anisotropic media—1. results for models of aligned cracks in the upper crust, *Geophys. J.*, **96**, 51–64.

Shuvalov, A.L., 1998. Topological features of the polarization fields of plane acoustic waves in anisotropic media, *Proc. R. Soc. Lond. A*, **454**, 2911–2947.

Vavryčuk, V., 2001. Ray tracing in anisotropic media with singularities, *Geophys. J. Int.*, **145**, 265–276.

Vavryčuk, V., 2003. Parabolic lines and caustics in homogeneous weakly anisotropic solids, *Geophys. J. Int.*, **152**, 318–334.

## APPENDIX A: EXPRESSIONS FOR THE ELEMENTS OF THE MATRIX $\tilde{B}$ IN ORTHORHOMBIC, TI AND ISOTROPIC MEDIA

In this appendix we restrict ourselves to orthorhombic media. An orthorhombic medium is defined by nine independent density-normalized elastic parameters $A_{L}$ and three mutually perpendicular planes of symmetry.

In the 'crystal' coordinate system, the elastic matrix is given by:

$$ \mathbf{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 (A1) $$

We introduce the following parameters (Mensch & Farra 1999):

$$ \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}, \quad \hat{A}_{23} = A_{23} - \frac{A_{22} + A_{33}}{2} + 2A_{44}. \qquad (A2) $$

Let $\mathbf{n}$ be the unit wave vector. We define three mutually perpendicular unit vectors $\tilde{e}_{[1]}$, $\tilde{e}_{[2]}$ and $\tilde{e}_{[3]} = \mathbf{n}$ in the following way:

$$ \tilde{e}_{[1]} = \frac{1}{n_r} \begin{bmatrix} n_1 n_3 \\ n_2 n_3 \\ -n_r^2 \end{bmatrix}, \quad \tilde{e}_{[2]} = \frac{1}{n_r} \begin{bmatrix} -n_2 \\ n_1 \\ 0 \end{bmatrix}, \qquad (A3) $$

with $n_r = \sqrt{n_1^2 + n_2^2}$.

Denoting $\theta$ and $\phi$, the polar angle and the azimuth specifying the vector $\mathbf{n}$, the components of the wave normal are $n_1 = \sin \theta \cos \phi$, $n_2 = \sin \theta \sin \phi$ and $n_3 = \cos \theta$. Moreover, we introduce the notation $\tilde{B}_{kl}(\mathbf{x}, \theta, \phi) = \tilde{B}_{kl}(\mathbf{x}, \mathbf{n})$.
In orthorhombic media, the elements of the symmetric matrix $\tilde{B}(\mathbf{x}, \theta, \phi)$ can be written as:

$$
\begin{align*}
\tilde{B}_{11} &= 2 (\hat{A}_{12} C_{\phi}^2 S_{\phi}^2 - \hat{A}_{13} C_{\phi}^2 - \hat{A}_{23} S_{\phi}^2) S_{\theta}^2 C_{\theta}^2 + A_{55} C_{\phi}^2 + A_{44} S_{\phi}^2, \\
\tilde{B}_{12} &= -C_{\phi} S_{\phi} C_{\theta} \{[\hat{A}_{23} - \hat{A}_{13} - \hat{A}_{12} (C_{\phi}^2 - S_{\phi}^2)] S_{\theta}^2 + A_{55} - A_{44}\}, \\
\tilde{B}_{22} &= (A_{44} C_{\phi}^2 + A_{55} S_{\phi}^2) C_{\theta}^2 + (A_{66} - 2 \hat{A}_{12} C_{\phi}^2 S_{\phi}^2) S_{\theta}^2, \\
\tilde{B}_{13} &= C_{\theta} S_{\theta} \left[ 2 \hat{A}_{12} C_{\phi}^2 S_{\phi}^2 S_{\theta}^2 + (\hat{A}_{13} C_{\phi}^2 + \hat{A}_{23} S_{\phi}^2) (C_{\theta}^2 - S_{\theta}^2) + \frac{1}{2} (A_{11} - A_{33}) C_{\phi}^2 + \frac{1}{2} (A_{22} - A_{33}) S_{\phi}^2 \right], \\
\tilde{B}_{23} &= S_{\theta} C_{\phi} S_{\phi} \left[ \hat{A}_{12} (C_{\phi}^2 - S_{\phi}^2) S_{\theta}^2 + (\hat{A}_{23} - \hat{A}_{13}) C_{\theta}^2 + \frac{1}{2} (A_{22} - A_{11}) \right], \\
\tilde{B}_{33} &= (A_{11} C_{\phi}^2 + A_{22} S_{\phi}^2) S_{\theta}^2 + A_{33} C_{\theta}^2 + 2 \hat{A}_{12} C_{\phi}^2 S_{\phi}^2 S_{\theta}^4 + 2 (\hat{A}_{13} C_{\phi}^2 + \hat{A}_{23} S_{\phi}^2) S_{\theta}^2 C_{\theta}^2,
\end{align*}
\tag{A4}
$$

where $C_{\theta}, S_{\theta}, C_{\phi}$ and $S_{\phi}$ mean cos $\theta$, sin $\theta$, cos $\phi$ and sin $\phi$, respectively.

Since the elements $\tilde{B}_{kl}$ are homogeneous functions of degree 2 with respect to $\mathbf{p}$, one can use the relation $\tilde{B}_{kl}(\mathbf{x}, \mathbf{p}) = p^2 \tilde{B}_{kl}(\mathbf{x}, \theta, \phi)$ to obtain their expressions from (A4).

In a VTI medium, one has $A_{11} = A_{22}$, $\hat{A}_{12} = 0$, $\hat{A}_{13} = \hat{A}_{23}$ and $A_{44} = A_{55}$. The elements $\tilde{B}_{ij}$ only depend on the angle $\theta$ and are given by:

$$
\begin{align*}
\tilde{B}_{11} &= A_{44} - 2 \hat{A}_{13} S_{\theta}^2 C_{\theta}^2, \\
\tilde{B}_{22} &= A_{66} S_{\theta}^2 + A_{44} C_{\theta}^2, \\
\tilde{B}_{13} &= C_{\theta} S_{\theta} \left[ \hat{A}_{13} (C_{\theta}^2 - S_{\theta}^2) + \frac{1}{2} (A_{11} - A_{33}) \right], \\
\tilde{B}_{12} &= \tilde{B}_{23} = 0, \\
\tilde{B}_{33} &= A_{11} S_{\theta}^2 + A_{33} C_{\theta}^2 + 2 \hat{A}_{13} S_{\theta}^2 C_{\theta}^2.
\end{align*}
\tag{A5}
$$

In isotropic media, $\hat{A}_{12} = \hat{A}_{13} = \hat{A}_{23} = 0$, $A_{44} = A_{55} = A_{66}$ and $A_{11} = A_{22} = A_{33}$. The corresponding elements $\tilde{B}_{ij}$ are given by:

$$
\tilde{B}_{11} = \tilde{B}_{22} = A_{44}, \quad \tilde{B}_{33} = A_{33}, \quad \tilde{B}_{12} = \tilde{B}_{13} = \tilde{B}_{23} = 0. \tag{A6}
$$

## APPENDIX B: EXPRESSIONS FOR THE PARTIAL DERIVATIVES OF ELEMENTS OF MATRIX B

Because the elements of matrix **B** are homogeneous functions of degree 2 with respect to **p**, their partial derivatives with respect to $p_i$ can be expressed in a simple way by using the spherical coordinates $(p, \theta, \phi)$ of vector **p**.

Let us introduce the notation $\mathcal{B}_{kl}(\mathbf{x}, \theta, \phi) = \mathcal{B}_{kl}(\mathbf{x}, \mathbf{n})$. The elements $\mathcal{B}_{kl}$ satisfy the relation $\mathcal{B}_{kl}(\mathbf{x}, \mathbf{p}) = p^2 \mathcal{B}_{kl}(\mathbf{x}, \theta, \phi)$. The first partial derivatives of element $\mathcal{B}_{kl}$ are given by:

$$
\begin{align}
\frac{\partial \mathcal{B}_{kl}}{\partial p_i}(\mathbf{x}, \mathbf{p}) &= p \left( 2 \mathcal{B}_{kl} n_i + \frac{\partial \mathcal{B}_{kl}}{\partial \theta} \bar{e}_{[1]i} + \frac{1}{\sin \theta} \frac{\partial \mathcal{B}_{kl}}{\partial \phi} \bar{e}_{[2]i} \right), \notag \\
\frac{\partial \mathcal{B}_{kl}}{\partial x_i}(\mathbf{x}, \mathbf{p}) &= p^2 \frac{\partial \mathcal{B}_{kl}}{\partial x_j}, \quad i=1, 2, 3, \tag{B1}
\end{align}
$$

where the vectors $\bar{e}_{[K]}$ are defined in Appendix A.

The second partial derivatives of element $\mathcal{B}_{kl}$ are given by:

$$
\begin{align*}
\frac{\partial^2 \mathcal{B}_{kl}}{\partial x_i \partial x_j}(\mathbf{x}, \mathbf{p}) &= p^2 \frac{\partial^2 \mathcal{B}_{kl}}{\partial x_i \partial x_j}, \\
\frac{\partial^2 \mathcal{B}_{kl}}{\partial p_i \partial x_j}(\mathbf{x}, \mathbf{p}) &= p \left( 2 \frac{\partial \mathcal{B}_{kl}}{\partial x_j} n_i + \frac{\partial^2 \mathcal{B}_{kl}}{\partial x_j \partial \theta} \bar{e}_{[1]i} + \frac{1}{\sin \theta} \frac{\partial^2 \mathcal{B}_{kl}}{\partial x_j \partial \phi} \bar{e}_{[2]i} \right), \\
\frac{\partial^2 \mathcal{B}_{kl}}{\partial p_i \partial p_j}(\mathbf{x}, \mathbf{p}) &= 2 \mathcal{B}_{kl} \delta_{ij} + \frac{\partial^2 \mathcal{B}_{kl}}{\partial \theta^2} \bar{e}_{[1]i} \bar{e}_{[1]j} + \frac{1}{\sin \theta} \left( \frac{\partial^2 \mathcal{B}_{kl}}{\partial \theta \partial \phi} - \frac{\cos \theta}{\sin \theta} \frac{\partial \mathcal{B}_{kl}}{\partial \phi} \right) (\bar{e}_{[1]i} \bar{e}_{[2]j} + \bar{e}_{[2]i} \bar{e}_{[1]j}) \\
&\quad + \frac{\partial \mathcal{B}_{kl}}{\partial \theta} (\bar{e}_{[1]i} n_j + n_i \bar{e}_{[1]j}) + \frac{1}{\sin \theta} \left( \frac{1}{\sin \theta} \frac{\partial^2 \mathcal{B}_{kl}}{\partial \phi^2} + \cos \theta \frac{\partial \mathcal{B}_{kl}}{\partial \theta} \right) \bar{e}_{[2]i} \bar{e}_{[2]j} \\
&\quad + \frac{1}{\sin \theta} \frac{\partial \mathcal{B}_{kl}}{\partial \phi} (\bar{e}_{[2]i} n_j + n_i \bar{e}_{[2]j}). \tag{B2}
\end{align*}
$$
# APPENDIX C: BEHAVIOUR OF THE FIRST-ORDER GROUP VELOCITIES AND ANGLE $\Phi^{(0)}$ IN THE NEIGHBOURHOOD OF A SINGULARITY

In order to study the behaviour of the first-order group velocities in the vicinity of singularities, we follow a similar approach to that used by Shuvalov (1998) (see also Vavryčuk 2003 for the exact expressions).

We shall use the following vectors $\mathbf{q}(\mathbf{p})$, $\mathbf{r}(\mathbf{p})$ and $\mathbf{s}(\mathbf{p})$ defined by their components:

$$q_i = \frac{1}{2} \left( \frac{\partial B_{11}}{\partial p_i} + \frac{\partial B_{22}}{\partial p_i} \right), \quad r_i = \frac{1}{2} \left( \frac{\partial B_{11}}{\partial p_i} - \frac{\partial B_{22}}{\partial p_i} \right), \quad s_i = \frac{\partial B_{12}}{\partial p_i}, \qquad (C1)$$

and the 3 × 3 matrices $\mathbf{F}(\mathbf{p})$ and $\mathbf{G}(\mathbf{p})$ whose elements are:

$$F_{ij} = \frac{1}{4} \left( \frac{\partial^2 B_{11}}{\partial p_i \partial p_j} - \frac{\partial^2 B_{22}}{\partial p_i \partial p_j} \right), \quad G_{ij} = \frac{1}{2} \frac{\partial^2 B_{12}}{\partial p_i \partial p_j}. \qquad (C2)$$

The components $q_i$, $r_i$, $s_i$ and the matrix elements $F_{ij}$ and $G_{ij}$ are continuous functions of $\mathbf{p}$. Let us remark that the vector $\mathbf{q}$ is independent of the choice of the vectors $\mathbf{e}_{[K]}$ used to define the elements $B_{KL}$, but $\mathbf{r}$ and $\mathbf{s}$ depend on it.

The first-order group velocities $v_{[M]}^{(1)}$ defined by (33) and the quantity $g_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_i} g_{[2]}^{(0)}$ needed in the evaluation of the second partial derivatives of $G_{[M]}^{(1)}$ (see eq. 38), can be written in terms of the angle $\Phi^{(0)}$ and the vectors $\mathbf{q}$, $\mathbf{r}$ and $\mathbf{s}$:

$$\mathbf{v}_{[M]}^{(1)}(\mathbf{n}) = \frac{1}{2V_{[M]}^{(1)}(\mathbf{n})} \left[ \mathbf{q}(\mathbf{n}) - (-1)^M (\cos 2\Phi^{(0)}(\mathbf{n})\mathbf{r}(\mathbf{n}) + \sin 2\Phi^{(0)}(\mathbf{n})\mathbf{s}(\mathbf{n})) \right], \quad M = 1, 2, \qquad (C3)$$

and

$$\mathbf{g}_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_i} \mathbf{g}_{[2]}^{(0)} = -\sin 2\Phi^{(0)} r_i + \cos 2\Phi^{(0)} s_i, \quad i = 1, 2, 3. \qquad (C4)$$

Let us denote by $\mathbf{n}_0$ the unit vector specifying the first-order singularity direction. In the direction $\mathbf{n}_0$, the eqs (20) are satisfied. Denoting $\mathbf{q}_0 = \mathbf{q}(\mathbf{n}_0)$, $\mathbf{r}_0 = \mathbf{r}(\mathbf{n}_0)$ and $\mathbf{s}_0 = \mathbf{s}(\mathbf{n}_0)$ the vectors calculated at $\mathbf{n}_0$, one can show from (20) and (B1) that:

$$\mathbf{r}_0 \cdot \mathbf{n}_0 = \mathbf{s}_0 \cdot \mathbf{n}_0 = 0, \qquad (C5)$$

so that the vectors $\mathbf{r}_0$ and $\mathbf{s}_0$ lie in the same plane orthogonal to $\mathbf{n}_0$.

Let us consider a close neighbourhood of the first-order singularity direction defined by the vector $\mathbf{n}_0$. The deviation $\Delta \mathbf{n} = \mathbf{n} - \mathbf{n}_0$ of the wave normal with respect to the singularity direction $\mathbf{n}_0$ is approximated by the relation:

$$\Delta \mathbf{n} = \psi \mathbf{e}_{\psi} - \frac{\psi^2}{2} \mathbf{n}_0, \qquad (C6)$$

where $\psi$ is the angle between $\mathbf{n}$ and $\mathbf{n}_0$ and $\mathbf{e}_{\psi}$ is the unit vector orthogonal to $\mathbf{n}_0$ specifying the direction from which the singularity is approached.

In the vicinity of the first-order singularity direction, one can write the following approximations:

$$\begin{aligned} B_{11}(\mathbf{n}) - B_{22}(\mathbf{n}) &\approx 2\psi (\mathbf{r}_0 \cdot \mathbf{e}_{\psi} + \psi \mathbf{e}_{\psi}^{\mathrm{T}} \mathbf{F}_0 \mathbf{e}_{\psi}), \\ B_{12}(\mathbf{n}) &\approx \psi (\mathbf{s}_0 \cdot \mathbf{e}_{\psi} + \psi \mathbf{e}_{\psi}^{\mathrm{T}} \mathbf{G}_0 \mathbf{e}_{\psi}), \\ \mathbf{r}(\mathbf{n}) &\approx \mathbf{r}_0 + 2\psi \mathbf{F}_0 \mathbf{e}_{\psi}, \\ \mathbf{s}(\mathbf{n}) &\approx \mathbf{s}_0 + 2\psi \mathbf{G}_0 \mathbf{e}_{\psi}, \end{aligned} \qquad (C7)$$

where the matrices $\mathbf{F}_0$ and $\mathbf{G}_0$ are calculated at $\mathbf{n}_0$. Moreover we denote by $V_S^{(1)}$ the common first-order phase velocity at the singularity.

## Conical singularity

A first-order conical singularity is defined as a direction in which the two first-order slowness sheets of the $qS$ waves touch through the vertices of cone-shaped surfaces. The condition for such a singularity can be written as follows:

$$\mathbf{r}_0 \times \mathbf{s}_0 \neq \mathbf{0}. \qquad (C8)$$

In the vicinity of the singularity, one can write the following approximations from (25), (27) and (C7):

$$\begin{aligned} \tan 2\Phi^{(0)}(\mathbf{n}) &\approx \frac{\mathbf{s}_0 \cdot \mathbf{e}_{\psi}}{\mathbf{r}_0 \cdot \mathbf{e}_{\psi}}, \\ G_{[1]}^{(1)}(\mathbf{n}) - G_{[2]}^{(1)}(\mathbf{n}) &\approx 2\psi \frac{\mathbf{r}_0 \cdot \mathbf{e}_{\psi}}{\cos 2\Phi^{(0)}}, \\ \mathbf{r}(\mathbf{n}) &\approx \mathbf{r}_0, \\ \mathbf{s}(\mathbf{n}) &\approx \mathbf{s}_0. \end{aligned} \qquad (C9)$$
From (C9), one can see that the angle $\Phi^{(0)}$ depends on the direction of approach $\mathbf{e}_{\psi}$ and changes by $\pm\pi/2$ when the wave normal crosses the singularity. In order to demonstrate it, one should remember that $G_{[1]}^{(1)}(\mathbf{n}) \ge G_{[2]}^{(1)}(\mathbf{n})$ (see eq. 18); therefore $\cos 2\Phi^{(0)}$ changes its sign as $\psi$ does when the wave normal crosses the singularity.

For $\mathbf{n}$ approaching $\mathbf{n}_0$, the first-order group-velocity vectors $\mathbf{v}_{[M]}^{(1)}$ tend to the vectors

$$ \frac{1}{2V_S^{(1)}} [\mathbf{q}_0 - (-1)^M (\cos 2\Phi^{(0)}\mathbf{r}_0 + \sin 2\Phi^{(0)}\mathbf{s}_0)], \quad M = 1, 2, \qquad (C10) $$

which envelop a cone (the so-called cone of internal refraction) with the elliptical lid lying in the plane orthogonal to $\mathbf{n}_0$. The $qS_1$ and the $qS_2$ waves exchange their first-order group velocity (because $\Phi^{(0)}$ changes by $\pm\pi/2$) when the wave normal crosses the singularity. Moreover, in the vicinity of the conical singularity, the quantity

$$ \left( \mathbf{g}_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_i} \mathbf{g}_{[2]}^{(0)} \right) \left( \mathbf{g}_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_j} \mathbf{g}_{[2]}^{(0)} \right) (G_{[1]}^{(1)} - G_{[2]}^{(1)})^{-1}, \quad i, j = 1, 2, 3, \qquad (C11) $$

behaves like $1/\psi$ (see C4 and C9), and is therefore singular at the singularity.

## Kiss singularity

A first-order kiss singularity (or tangential singularity) corresponds to a direction in which the first-order slowness sheets of the $qS$ waves touch tangentially at an isolated point. Such a singularity corresponds to the following condition:

$$ \mathbf{r}_0 = \mathbf{s}_0 = \mathbf{0}. \qquad (C12) $$

In the vicinity of the singularity, one can write the following approximations from (25), (27) and (C7):

$$ \tan 2\Phi^{(0)}(\mathbf{n}) \approx \frac{\mathbf{e}_{\psi}^{\mathrm{T}} \mathbf{G}_{0} \mathbf{e}_{\psi}}{\mathbf{e}_{\psi}^{\mathrm{T}} \mathbf{F}_{0} \mathbf{e}_{\psi}}, $$

$$ G_{[1]}^{(1)}(\mathbf{n}) - G_{[2]}^{(1)}(\mathbf{n}) \approx 2\psi^2 \frac{\mathbf{e}_{\psi}^{\mathrm{T}} \mathbf{F}_{0} \mathbf{e}_{\psi}}{\cos 2\Phi^{(0)}}, $$

$$ \mathbf{r}(\mathbf{n}) \approx 2\psi \mathbf{F}_{0} \mathbf{e}_{\psi}, $$

$$ \mathbf{s}(\mathbf{n}) \approx 2\psi \mathbf{G}_{0} \mathbf{e}_{\psi}. \qquad (C13) $$

From (C13), one can see that the angle $\Phi^{(0)}$ depends on the direction of approach $\mathbf{e}_{\psi}$ but is continuous when the wave normal crosses the singularity.

For $\mathbf{n}$ approaching $\mathbf{n}_0$, the first-order group-velocity vectors $\mathbf{v}_{[M]}^{(1)}$ tend to the same vector, see (C3) and (C12):

$$ \frac{1}{2V_S^{(1)}} \mathbf{q}_0, \qquad (C14) $$

whatever the wave mode $M$ and the direction of approach $\mathbf{e}_{\psi}$.

In the vicinity of the kiss singularity, the quantity

$$ \left( \mathbf{g}_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_i} \mathbf{g}_{[2]}^{(0)} \right) \left( \mathbf{g}_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_j} \mathbf{g}_{[2]}^{(0)} \right) (G_{[1]}^{(1)} - G_{[2]}^{(1)})^{-1}, \quad i, j = 1, 2, 3, \qquad (C15) $$

has a finite limit at the singularity, see (C4) and (C13).

## Intersection singularity

An intersection singularity occurs when the first-order slowness sheets of the two $qS$ waves intersect along a line. Such a singularity corresponds to the following condition:

$$ \mathbf{s}_0 = \eta \mathbf{r}_0, \mathbf{r}_0 \neq \mathbf{0}. \qquad (C16) $$

In the vicinity of the singularity, one can write the following approximations from (25), (27) and (C7):

$$ \tan 2\Phi^{(0)}(\mathbf{n}) \approx \eta, $$

$$ G_{[1]}^{(1)}(\mathbf{n}) - G_{[2]}^{(1)}(\mathbf{n}) \approx 2\psi \frac{\mathbf{r}_0 \cdot \mathbf{e}_{\psi}}{\cos 2\Phi^{(0)}}, $$

$$ \mathbf{r}(\mathbf{n}) \approx \mathbf{r}_0 + 2\psi \mathbf{F}_{0} \mathbf{e}_{\psi}, $$

$$ \mathbf{s}(\mathbf{n}) \approx \mathbf{s}_0 + 2\psi \mathbf{G}_{0} \mathbf{e}_{\psi}. \qquad (C17) $$

From eqs (C17), one can see that the angle $\Phi^{(0)}$ changes by $\pm\pi/2$ when the wave normal crosses the singularity.
For $\mathbf{n}$ approaching $\mathbf{n}_0$, the first-order group-velocity vectors $\mathbf{v}_{[M]}^{(1)}$ tend to the vectors

$$ \frac{1}{2V_S^{(1)}} \left( \mathbf{q}_0 - \frac{(-1)^M}{\cos 2\Phi^{(0)}} \mathbf{r}_0 \right). \qquad (C18) $$

The $qS_1$ and the $qS_2$ waves exchange their group velocity (because $\Phi^{(0)}$ changes by $\pm\pi/2$) when the wave normal crosses the singularity. In the vicinity of the intersection singularity, we can write the following approximation of the quantity $g_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_i} g_{[2]}^{(0)}$:

$$ \mathbf{g}_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_i} \mathbf{g}_{[2]}^{(0)} \approx 2\psi (-\sin 2\Phi^{(0)} \mathbf{F}_0 \mathbf{e}_\psi + \cos 2\Phi^{(0)} \mathbf{G}_0 \mathbf{e}_\psi)_i, \qquad (C19) $$

(see C4 and C17), so that

$$ \left( \mathbf{g}_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_i} \mathbf{g}_{[2]}^{(0)} \right) \left( \mathbf{g}_{[1]}^{(0)T} \frac{\partial \mathbf{B}}{\partial p_j} \mathbf{g}_{[2]}^{(0)} \right) (G_{[1]}^{(1)} - G_{[2]}^{(1)})^{-1}, \quad i, j = 1, 2, 3, \qquad (C20) $$

has a finite limit at the singularity.

## Conclusion

Except for the kiss singularity, the angle $\Phi^{(0)}$ is discontinuous when the wave normal crosses the first-order singularity and the zero-order polarization vector $g_{[M]}^{(0)}$ of the $qS_M$ wave changes abruptly. Indeed, the angle $\Phi^{(0)}$ changes by $\pm\pi/2$ when the singularity is crossed, so that the $qS_1$ and $qS_2$ waves exchange the direction of their zero-order polarization vectors and therefore exchange their group velocity. In order to have a continuous variation of the first-order group velocity, when the wave normal crosses a conical or intersection singularity, the eigenvalue should be changed ($G_{[1]}^{(1)}$ to $G_{[2]}^{(1)}$ or vice versa). Moreover, in the case of the conical singularity, the second partial derivatives of the eigenvalues $G_{[M]}^{(1)}$ with respect to $\mathbf{p}$ components are infinite and the paraxial ray tracing system is singular.