Phase shift at caustics along rays in anisotropic media

Abstract

In isotropic ray tracing, the ray approximation to the wavefield undergoes a phase shift when the ray crosses a caustic. The cumulative number of such phase shifts along a ray is usually called the KMAH index. The sign of these phase shifts is prescribed by the sign of the angular frequency in combination with the sign convention used for the Fourier transformation. In isotropic media the KMAH index always increases by one or by two, depending on the type of caustic crossed. For (quasi-)shear waves in anisotropic media the KMAH index may decrease. This is the case if the associated slowness sheet is locally concave in one or two of its principal directions of curvature.

1 Introduction

In isotropic ray tracing, it is well known that a phase shift occurs in the ray approximation to the wavefield when a ray propagates through a caustic. At a caustic the cross-sectional area of an infinitesimally thin ray tube around a central ray vanishes. Such a cross-sectional area is modelled by the determinant det [Q(s) ] of a 2 × 2 matrix which occurs in the dynamic ray-tracing formalism of Červený (1985), where s is the curve length along the ray. In isotropic media, assuming a time-harmonic phase factor of exp (−iωt), the phase is increased by −sgn (ω)π/2 for each zero encountered in det [Q(s) ], where double zeros are counted twice. For positive frequencies all phase shifts are in the same direction.

For (quasi-)shear waves in anisotropic media ‘anomalous’ phase shifts of opposite sign may occur. This phenomenon has been discussed previously by several authors, either implicitly [e.g. a very general treatment by Lewis (1965) ] or explicitly [e.g. Garmany (1988a,b), dealing with WKBJ approximations in stratified media, or Kravtsov & Orlov (1993), in a context of stationary-phase approximations]. Only recently, Klimeš (1997, 1998) described the phase shifts at caustics in anisotropic media in terms of traces and determinants of the fundamental 2 × 2 matrix solutions, computed in Červený’s formalism of dynamic ray tracing. The analysis was given for a bundle of rays which starts from a point-source. Klimeš’s analysis rests on several specific symplectic properties of the solution of the dynamic ray-tracing system. This paper gives an alternative, more geometrical description of the phase shifts at caustics in anisotropic media, fully in terms of the curvature of the wave front in the vicinity of the caustic. The sign of the phase shift turns out to be related to the local convexity or concavity of the associated slowness sheet. The criteria for computing the sign of the phase shift, as presented here, are shown to be consistent with Klimeš’s results, and these are valid for rather general ray bundles, rather than only for ray bundles which start at a point-source.

2 Fundamentals Of Anisotropic Dynamic Ray Tracing

Consider a bundle of rays, parametrized by two ray parameters γ₁ and γ₂ (e.g. shooting angles). We adopt the notation and concepts of isotropic dynamic ray tracing as formulated by Červený (1985). Let x(s) and p(s) denote the position of a ray in this bundle and the associated slowness vector (i.e. the gradient of traveltime τ), respectively, as a function of curve length s along the ray. In isotropic dynamic ray tracing, it is convenient to use a ray-centred system of orthonormal basis vectors {e₁ (s),e₂ (s),t(s) } Popov & Pšenčík 1978), where e₁ (s) and e₂ (s) span the tangent plane of the wave front at x(s), and t(s) is a unit vector tangent to the ray. Here we prefer to interpret t(s) as the normalized slowness vector, which allows a generalization of such a ray-centred formulation to anisotropic dynamic ray tracing (Červený 1997; Bakker 1996). Generally, t(s) is not tangent to the ray in the anisotropic case, but it is perpendicular to the wave front. Along a ray we denote the group (or ray) velocity by v_gr , and the phase velocity (i.e. the reciprocal of the length of p) by v_ph . In isotropic media one has v = v_gr = v_ph .

The dynamic ray-tracing system can be formulated in terms of components along these ray-centred basis vectors. Let Q(s) and P(s) be the 2 × 2 matrices with entries Q_IJ (s) = 〈e_I (s),∂x/∂γ_J (s) 〉 and P_IJ (s) = 〈e_I (s),∂p/∂γ_J (s) 〉, respectively. Then the dynamic ray-tracing system reads

The specific form of the matrices S_IJ [which can be computed on the basis of results in Gajewski & Pšenčík (1990) ] is not of importance here, except for the symmetry of S₁₂ and S₂₁ , and the equality S₁₁ = −S^T₂₂ . These properties imply the usual symplecticity of the dynamic ray-tracing system. Since the trace of the 4 × 4 matrix in eq.(1) vanishes, it follows from Liouville’s theorem that Q(s) and P(s) cannot have a common null vector for any s if they do not have a common null vector at the start of the ray. It can be shown that this remains true in cases where interfaces are crossed (Červený 1997). For purposes of interpretation, it is worthwhile to mention that S₁₁ , S₂₁ and S₂₂ vanish in a homogeneous medium. Furthermore, in an isotropic medium one has S₁₂ = v² I.

For both isotropic and anisotropic media, one defines the matrix M = PQ⁻¹ of the second-order derivatives of the traveltime in the directions tangent to the wavefront, and v_ph (s)M(s) is the curvature matrix of the wavefront at x(s).

Assuming a phase factor exp (−iωt) for Fourier transformation from the frequency domain to the time domain, the ray approximation to the wavefield along a ray in a smoothly varying medium is expressed by

where ψ is a constant vector, calibrated at the start, and ρ is the density. In an isotropic medium u should be interpreted as the coordinates of the particle displacement along the (varying) ray-centred basis vectors (e.g. Červený 1985), and in an anisotropic medium u should be interpreted as the coordinates along the (varying) polarization vectors associated with each of the three wave modes (e.g. Pšenčík & Teles 1996; Červený 1997). In isotropic media it is well known that 1/{det [Q(s) ] }^1/2 is to be interpreted as (1/ |det [Q(s) ] |^1/2 )exp [−isgn (ω) (π/2)k(s,s₀ ) ], where the KMAH index k(s,s₀ ) counts the number of zeros in det [Q(s) ] in the interval (s₀ ,s) (double zeros counted twice), assuming that the ray starts at s₀ .

Hence, for ω > 0, in isotropic media, the phase is decreased by π/2 for each zero encountered in det [Q(s) ]. In Section 3 it is shown, however, that this is not always the case for (quasi-)shear waves in anisotropic media.

3 Phase Shift At Caustics

Assuming that P(s₀ ) and Q(s₀ ) do not have a common null vector at the initial point s₀ , we can have the following three situations at a point s_c which is on a caustic:

(1)case 1: Q(s_c ) = 0 (double zero, i.e. point caustic) and P(s_c ) is of full rank;

(2)case 2: Q(s_c ) has rank 1 (single zero, i.e. line caustic) and P(s_c ) is of full rank;

(3)case 3: Q(s_c ) has rank 1 (single zero, i.e. line caustic) and P(s_c ) has rank 1.

The appropriate treatment of the phase shift at a caustic can be based on the paraxial approximation of the wavefield at x(s,q) = x(s) + q₁ e₁ (s) + q₂ e₂ (s) in a paraxial vicinity of a central ray. With q^T = (q₁ ,q₂ ), the wavefield is approximated by u(ω,s,q) = U(ω,s,q)exp [iωτ(s) ], where

In the vicinity of a caustic this paraxial approximation has to be interpreted as a limiting case for analytic continuation of the traveltime τ(s), where Im [M(s) ] > 0 is used for its second-order derivative, and ω > 0 is assumed. This requirement is consistent with the Gaussian-beam approach as presented by Červený & Pšenčík (1983) and Červený (1985).

By rewriting det (Q) as det (P)det (M⁻¹ ) in eq.(3), one finds that the singularities in the vicinity of the point s_c at a caustic are contained in the factor {det [M(s) ] }^1/2exp [iω(1/2)q^T M(s)q], provided that det [P(s_c ) ] ≠ 0. Expressions of this kind typically occur as stationary-phase approximations in Maslov theory. For a 2-D situation, it was stated by Kravtsov & Orlov (1993) that the phase factor of the square root depends on the sign of the quadratic form in the exponent. A similar conclusion can be drawn from Appendix F of Lewis (1965). Note that in anisotropic media the polarization vector depends smoothly on s, unless the ray accidentally hits the caustic with a slowness vector pointing into a singular direction of the associated slowness sheet. Such exceptions are not considered in this paper.

In the following, we shall analyse the singularity in the phase factor in more detail. Let μ₁ (s) and μ₂ (s) be the eigenvalues of M(s), with m₁ (s) and m₂ (s) being the associated normalized eigenvectors. We consider cases 1–3 listed above.

Case 1

As Q(s_c ) = 0, one has Q(s) ≈ [ (s − s_c )/v_gr (s_c ) ]S₁₂(s_c )P(s_c ) in the vicinity of the caustic, which implies M⁻¹ (s) ≈  [ (s − s_c )/v_gr (s_c ) ]S₁₂ (s_c ). This shows that the eigensystem of M⁻¹ (s) is smoothly varying if S₁₂ (s_c ) is not rank-deficient, with

where the σ_I are the eigenvalues of S₁₂ (s_c ). Consequently, one obtains

Now suppose that M(s), and hence its eigenvalues, are perturbed into the complex z domain, with positive imaginary parts. Then Fig. 1 shows the path travelled by [μ_I (s) ]⁻¹ when s passes the point s_c at the caustic in the forward direction, depending on the sign of the σ_I . One reads off the phase changes from the plotted paths of the principal branches of the square roots [μ_I (s) ]^−1/2, which should not cross the branch line formed by the negative real axis. In the limit for vanishing imaginary parts of the eigenvalues, the square roots [μ_I (s) ]^−1/2 should be multiplied by a phase factor exp [i(π/2)sgn (σ_I ) ] if s passes the point s_c in the forward direction.

Hence, the amplitude should be multiplied by a factor

Rays for which the matrix S₁₂ is rank-deficient at the caustic (i.e. the signs of its eigenvalues are questionable) are not considered in this paper.

Case 2

Now M⁻¹ (s_c ) has two distinct eigenvalues, and hence the eigensystem of M⁻¹ (s) depends smoothly on s in the vicinity of s_c . Again eq.(5) is valid, but now only one of the eigenvalues μ_I (s), say μ₁ , has a singularity. Starting with the equalities dM/ds = (dP/ds)Q⁻¹+P(dQ⁻¹/ds) and dQ⁻¹/ds = −Q⁻¹ (dQ/ds)Q⁻¹ , we investigate the singular behaviour of μ₁ . By substitution of the dynamic ray-tracing system (eq.1) into these relations, one obtains

and hence

where S₁₁ = −S^T₂₂ , and where S₁₂ and S₂₁ are symmetric. As the eigenvalues of M⁻¹ (s_c ) are well separated, one obtains v_gr [d(1/μ₁ )/ds]=v_gr m^T₁ (dM⁻¹ /ds)m₁ → (m^T₁ S₁₂ m₁ ) (s_c ) for s → s_c . Consequently, as an analogue to eq.(4), one obtains

which implies, similarly to case (1), that the amplitude has to be multiplied by

when s passes the point s_c in the forward direction. The exceptional case where (m^T₁ S₁₂ m₁ ) (s_c ) = 0 is not treated in this paper.

Case 3

The situation in which both Q(s_c ) and P(s_c ) have rank 1 may occur as an exception in 3-D, but if 2-D modelling is embedded in 3-D modelling (when P has always rank 1, i.e. no geometrical spreading out of the plane of the cross-section) this is just a normal case of a caustic.

In order to analyse this situation, let a be a null vector of P(s_c ), and let b be a null vector of Q(s_c ). Since M(s)Q(s)a = P(s)a→0 for s→s_c , and Q(s_c )a ≠ 0, one of the eigenvalues of M(s), say μ₂ , takes the value 0 at s_c . Similarly, one obtains M(s)Q(s)b→P(s_c )b ≠ 0, while Q(s)b→0 for s→s_c , which implies 1/μ₁ (s)→0. Furthermore, the eigenvectors of M(s) vary smoothly with s, and at s_c we can make the identification m₂ (s_c ) = Q(s_c )a and m₁ (s_c ) = P(s_c )b, where appropriate scaling of a and b is assumed such that the m_I are orthonormal.

In this case, rewriting det (Q) as det (P)det (M⁻¹ ) in eq.(3) is not convenient. Rather, we investigate Q in terms of the orthonormal eigenvectors m_I . Since μ₁ m^T₁ Q = m^T₁ MQ = m^T₁ P, we obtain

The determinant of the latter matrix is non-zero at s_c . Otherwise, this matrix would have a null vector m for which m^T₁ Pm = 0 and m^T₂ Qm = 0, which in turn would imply that m is a common null vector of P(s_c ) and Q(s_c ), since m₁ (s_c ) spans the range of P(s_c ), and m₂ (s_c ) spans the range of Q(s_c ) at s_c .

In this case the factor exp [iω(1/2)μ₂ (q^T m₂ )² ] behaves smoothly in the vicinity of s_c . Therefore, the phase jump at the caustic is determined by the phase jump of μ^1/2₁exp [iω(1/2)μ₁ (q^T m₁ )² ], and eqs (9) and (10) are again valid as a consequence of eq.(7), provided (m^T₁ S₁₂ m₁ ) (s_c ) ≠ 0. Thus, the present case 3 leads to the same phase jump as case 2, i.e. the vanishing of det (P) at the caustic does not influence the phase jump.

4 KMAH Index For Anisotropic Media

The results of the analysis in Section 3 can be summarized in terms of the KMAH index for anisotropic media. For ω > 0, and with the phase factor exp (−iωt) in the Fourier transformation, the factor 1/[det (Q) ]^1/2 is rewritten as {1/ |det [Q(s) ] |^1/2 }exp [−i(π/2)k(s,s₀ ) ], where the KMAH index k(s,s₀ ) accumulates the phase shifts in the following way.

(1)If S₁₂ (s_c ) is positive definite, then k(s,s₀ ) is increased by 2 in the case of a point caustic (case 1), and by 1 in the case of a line caustic (cases 2 and 3).

(2)If S₁₂ (s_c ) is negative definite, then k(s,s₀ ) is decreased by 2 in the case of a point caustic (case 1), and by 1 in the case of a line caustic (cases 2 and 3).

(3)If S₁₂ (s_c ) is neither positive definite nor negative definite, then k(s,s₀ ) remains unchanged for a point caustic (case 1), and in the case of a line caustic (cases 2 and 3) the value of sgn [ (m^T₁ S₁₂ m₁ ) (s_c ) ] is added. Some ambiguous situations are not considered here: the situation for a point caustic in which S₁₂ (s_c ) has an eigenvalue equal to 0, and the situation for a line caustic in which (m^T₁ S₁₂ m₁ ) (s_c ) = 0.

Some interpretation of this treatment of the KMAH index is appropriate here. It can be shown (Bakker 1996; Červený 1997) that the slowness sheet of the wave mode of interest associated with the stiffness tensor at x(s_c ) is convex or concave at p(s_c ) if and only if S₁₂ (s_c ) is positive or negative definite, respectively. Thus, the number of positive or negative eigenvalues of S₁₂ (s_c ) equals the number of positive or negative principal curvatures of the relevant slowness sheet at p(s_c ). Moreover, the analysis of Section 3 shows that the values (m^T_I S₁₂ m_I ) (s_c ) determine whether the singular principal curvatures of the wave front change from negative to positive at the caustic, or vice versa.

Note that in an isotropic medium S₁₂ is a scalar multiple of the unit matrix, and therefore positive definite. This implies that the slowness sheets are always convex, and that the singular principal curvatures of a wave front always change from negative to positive values when a caustic is crossed in the forward direction. In anisotropic media the inner slowness sheet (for the P wave) is always convex if it is detached from the other two (for the S waves) (e.g. Duff 1960). Hence anisotropic P waves will have the usual phase-shift behaviour, but anisotropic S waves may have anomalous phase changes, as their slowness sheets can be locally concave in one or two principal directions.

A similar link between the sign of the phase shift and the convexity/concavity of the slowness sheets was made earlier by Garmany (1988a,b) for stratified media. Note also that the results shown here for point caustics are consistent with the so-called index of a point-source, described by Burridge (1967).

Acknowledgements

The author thanks Shell International Exploration and Production BV for permission to publish this paper. Two reviewers are acknowledged for some detailed comments.

References

Appendix

Appendix A:Relation To Klimeš’s Results

Klimeš (1997) has analysed the phase shift at caustics for a bundle of rays generated by a point source.

For a point caustic (case 1) Klimeš’s results are readily seen to be equivalent to the results presented in this paper. For the point caustic Klimeš uses the matrix [ (dQ/ds)P⁻¹ ] (s_c ) instead of S₁₂ (s_c ). This is evidently the same because Q(s_c ) = 0 in such a case, and hence dQ/ds = (1/v_gr ) (S₁₂ P) at s_c .

For a line caustic (cases 2 and 3) the comparison is somewhat more tricky. In these cases Klimeš concludes that sgn [K₁ (s_c )K₂ (s_c ) ] has to be added to the KMAH index, where

with K = Q⁻¹det (Q) (which varies smoothly when s crosses s_c ). Hence, the question is whether sgn (K₁ K₂ ) equals sgn (m^T₁ S₁₂ m₁ ) at s_c .

For case 2, where P(s_c ) has full rank, m₂ (s_c ) spans the range of M⁻¹ (s_c ), which equals (QP⁻¹ ) (s_c ). Since Q(s_c ) has rank 1, m₂ (s_c ) also spans the range of Q(s_c ). This latter property was also shown for case 3, where P(s_c ) has rank 1. Moreover, the range of Q(s_c ) equals the null space of K(s_c ), and vice versa. Consequently, m₂ (s_c ) is in the null space of K(s_c ).

Since tr (Q⁻¹ P) = tr (QQ⁻¹PQ⁻¹ ), we have

Furthermore, K₂ = tr (KS₁₁ Q + KS₁₂ P). Here, the first term vanishes at s_c , because tr (KS₁₁ Q) = tr (KS₁₁ K⁻¹ )det (Q) = tr (S₁₁) det(Q). For the second term we have tr(KS₁₂P) = tr (K⁻¹KS₁₂PK) = tr (S₁₂ PK), which equals m^T₁ S₁₂ PKm₁ at s_c . Hence,

Since m₁ is an eigenvector of PQ⁻¹ , it is also an eigenvector of PK. Now, combination of eqs (A2) and (A3) gives the desired equality