## Abstract

Slowness surface for bulk wave propagation in anisotropic media can be divided into concave, saddle and convex regions by the parabolic lines. When a parabolic line crosses a symmetry plane, it leaves either an inflection point or a parabolic point. Surface normal at these points is associated with cuspidal point and swallowtail point, respectively, on the wave surface and in phonon focusing patterns. By examining the degeneracies in the Stroh eigenvalue equation, we have calculated the cuspidal points in cubic crystals analytically. In this work, the parabolic point and its surface normal are discussed. The main idea is to establish a connection between the parabolic point and the extraordinary transonic state that is related to a degeneracy with a multiplicity of four in the Stroh eigenvalue equation. Such a connection yields a series of simple expressions, which determine the locations of parabolic points and the corresponding swallowtail points. The result is demonstrated using phonon focusing patterns of cubic crystals, and the method also provides a tool for general discussion of the slowness surface geometry.

## 1. Introduction

Slowness surface, defined by the Christoffel equation for the bulk wave propagation in anisotropic elastic media, exhibits many interesting features (Musgrave 1970; Auld 1973). It is well known that each sheet of the slowness surface can be divided into concave, saddle and convex regions separated by the so-called parabolic line along which the Gaussian curvature is zero (Every 1981; Wolfe 1998). There are also various kinds of acoustic axes where two slowness sheets meet. A proper way to examine the geometry of slowness surfaces is to study the signatures left by the parabolic lines at symmetry planes. In a symmetry plane, the slowness surface may have points with either in-plane zero curvature or ex-plane zero curvature. The former refers to an inflection point marking a fold on the slowness surface and the latter refers to a parabolic point as shown in figure 1.

Wave surface, defined by the envelope of the slowness surface, describes group velocity for bulk waves. The geometry of the wave surface can be studied by examining either the Gauss map theoretically or phonon focusing patterns experimentally (Northrop & Wolfe 1980; Every 1981, 1986; Hurley & Wolfe 1985; Wolfe 1998). The phonon focusing pattern is just a polar projection of the wave surface, which consists mainly of two types of caustic lines: (i) a parabola-formed line centred at a *cuspidal point* originated from the inflection point and (ii) a pair of semicubical formed lines centred at a *swallowtail point* resulting from the parabolic point.

Geometrically, it was suggested that transformation from the slowness surface to the wave surface can be addressed by the Legendre transformation in catastrophe theory (Sewell 1982). However, an explicit scheme remains lacking for dealing with the slowness surface in anisotropic media. By examining the Christoffel equation, Every (1981, 1986) derived a series of critical existence conditions for the cuspidal/swallowtail points in cubic crystals (see also lines A, B, D, E and G in the figures 5 and 9), and Shuvalov & Every (1996, 2000) succeeded in identifying some transverse zero-curvature points in symmetry planes using the perturbation method. But exact expressions for positions of cuspidal/swallowtail points remain unsolved explicitly. Other studies on the wave surface and its geometry have been limited in numerical simulations in relation to phonon imaging (Northrop & Wolfe 1980; Hurley & Wolfe 1985; Wolfe 1998).

In a recent study (Wang 2008), the present author established a scheme to calculate the positions of the cuspidal points in cubic crystals based on the Stroh (1962) formalism for elastodynamics. The main idea was to recognize the connection between the inflection points on the slowness surface and the triple degeneracies in the sextic Stroh eigenvalue equation. In the case of the parabolic points, the scheme has to be extended. In §2, we will explore the relationship between the parabolic point and the so-called zero-curvature transonic state, which is associated with a degeneracy of multiplicity of four in the Stroh eigenvalue equation (Barnett *et al*. 1990). Such a relationship will make the determination of the parabolic points easier compared to the perturbation method (Shuvalov & Every 1996, 2000). The swallowtail point on the wave surface, originated from the surface normal at the parabolic point, will then be calculated analytically. The method developed here can be applied to crystals with other symmetries.

## 2. Stroh formalism and degeneracies

Consider an infinite homogeneous elastic medium. The equation of motion for steady states is given by(2.1)where *C*_{ijkl} is the elastic stiffness tensor; *ρ* is the density; and *u* is the displacement. We define a reference plane *R*=(** m**,

**), where**

*n***and**

*m***are two unit vectors and**

*n***⊥**

*m***, and consider solutions described by**

*n***(**

*u***)=**

*x***exp{i**

*a**k*[(

**+**

*m**p*

**).**

*n***−**

*x**vt*]} with the wavevector

**=**

*k**k*(

**+**

*m**p*

**n**). The polarization vector

**will obey(2.2)where the term (**

*a**mn*), etc. are given by (

*mn*)

_{ik}

*=m*

_{j}

*c*

_{ijkl}

*n*

_{l}, and the repeated subscripts denote summation. Traction vector

**, which describes the traction exerted by the wave on the planes parallel to**

*b***.**

*n***=0, is defined by**

*x***=−[(**

*b**nm*)+

*p*(

*nn*)]

**. By defining a six-dimensional vector**

*a**ξ*=[

**,**

*a***]**

*b*^{T}, Stroh (1962) reformulated (2.2) into a six-dimensional eigenvalue equation,(2.3)whereand

*I*is a 3×3 identity matrix. Note that the characteristic equations for (2.2) and (2.3) are identical. Generally, the Stroh eigenvalues

*p*'s appear as three pairs of complex conjugate numbers when

*v*is sufficiently low. As

*v*increases, certain pair(s) of eigenvalues become real, giving bulk wave solutions with real wavevector

**=**

*k**k*(

**+**

*m**p*

**). Transition of the eigenvalues from complex to real, or vice versa, is usually referred as a transonic state, and it can be interpreted geometrically as follows (figure 2**

*n**a*): construct a vertical line

*L*parallel to

**and place it at large**

*n**v*(or small

*v*

^{−1}), so that it intersects a slowness curve leaving two intersection points marking two bulk wave solutions (

*k*_{1},

*k*_{2}) associated with two real Stroh eigenvalues

*p*

_{1}=tan

*ϕ*

_{1}and

*p*

_{2}=tan

*ϕ*

_{2}, where

*ϕ*

_{i}are the inclination angles of wavevectors

*k*_{i}from

**. By reducing**

*m**v*, the line

*L*moves towards the right until it contacts tangentially the outermost point of the slowness curve, and it leads to a degeneracy

*p*

_{1}=

*p*

_{2}=

*p*,(2.4)or

*k*_{1},

*k*_{2}→

**. Decreasing**

*k**v*further will yield six complex eigenvalues. Such a velocity

*v*is denoted as the transonic velocity. Geometrically, since the line L forms a tangent at the contact point, the transonic velocity

*v*represents then the group velocity

*v*

_{g}which is parallel to

**. In other words, the degeneracy in the Stroh eigenvalue equation defines the curve normal and the group velocity simultaneously.**

*m*The transonic states have been studied extensively and they have been categorized into two groups: there are six ordinary transonic states (types 1–6) (Chadwick & Smith 1977; Ting 1996) and four extraordinary transonic states (types E1–E4) (Barnett *et al*. 1990). The ordinary transonic state corresponds to simple duplex degeneracy as shown in figure 2*a*,*b*, while the extraordinary transonic state is related to a degeneracy with multiplicity of four or six. Consider a symmetric type 4 transonic state shown in figure 2*b*. There are four Stroh eigenvalues coalescing into two duplex-degenerated ones when the line L contacts tangentially two points at the slowness curve with the decreasing *v*. If such a transonic state evolves into a symmetric type 1 transonic state, the transition can take place only via an intermediate state with the straight slowness curve, i.e. an E1 transonic state as shown in figure 2*c*. In other words, an E1 transonic state in fact represents a parabolic point. At such an E1 transonic state, *four* Stroh eigenvalues coalesce into one real number and the slowness curve then exhibits zero curvature. When ** m** lies in a symmetry plane, the existence of E1 transonic state or parabolic point requires the inclination angles

*ϕ*

_{i}for the four bulk waves

*k*_{i},

*i*=1, 2, 3, 4 being zero:(2.5)Such a relationship can then be regarded as an existence condition for the parabolic points in the symmetry planes. (Note that although the curve normal at an E1 transonic state lies in the symmetry plane, the surface normal might point sideways with respect to the basis vector

**as shown in figure 3.)**

*m*For cubic elastic media, the matrix *Γ* in (2.2) often has a simple form when the basis vector ** n** is normal to a symmetry plane. This enables us to reformulate the condition (2.5) into a set of simple equations determining locations of the parabolic points in the principal symmetry planes: (001) and . With parabolic points known, their surface normal will be determined using the Stroh formalism. For the sake of brevity, we will use the relative elastic constants

*a*=

*c*

_{11}/

*c*

_{44}and

*b*=

*c*

_{12}/

*c*

_{44}and the anisotropy factor

*Δ*=

*a*−

*b*−2 in our discussion.

## 3. Location of parabolic points

In order to locate the parabolic points, let us define the reference plane *R* with (** m**,

**)=(**

*n***,**

*k*

*n*_{1}) as shown in figure 1. By changing the propagation direction

**or**

*k**θ*, we scan slowness sections so that E1 zero-curvature transonic states can be identified using the condition (2.5).

As long as ** n** is normal to a symmetry plane, the characteristic equation for (2.2), ‖

*Γ*‖=0, can always be expressed in terms of a bicubic equation of

*p*,(3.1)A symmetric E1 extraordinary transonic state (with respect to

**), if it exists, requires a fourfold degeneracy with degenerated eigenvalue**

*m**p*=0. This condition for the zero-curvature transonic state, or the parabolic point, can now be reformulated as(3.2)The first condition,

*a*

_{0}=‖(

*mm*)−

*ρv*

^{2}

*I*‖=0, is in fact the characteristic equation of the Christoffel equation, and its solutions (in terms of

*ρv*

^{2}) define then the three

*transonic state branches*belonging to the elliptic sheet

*S*

_{b}and two non-elliptic sheets

*S*

_{a}and

*S*

_{c}. The second condition,

*a*

_{2}=0, ensures the degeneracy in

*p*with multiplicity of four, or a parabolic point. Algebraically, the condition (3.2) can be replaced by the resultant between

*a*

_{0}and

*a*

_{2}(see, e.g. Bôcher 1936),(3.3)because the vanishing resultant guarantees a common root between the two polynomials,

*a*

_{0}(

*ρv*

^{2},

*θ*) and

*a*

_{2}(

*ρv*

^{2},

*θ*), and it leaves Res(

*a*

_{0},

*a*

_{2}) dependent only on

*θ*.

### (a) Parabolic point in (001) plane

We will now calculate directions of the parabolic points in the (001) plane. Let ** m**=cos

*θ*

*e*_{x}+sin

*θ*

*e*_{y}and

**=**

*n*

*e*_{z}. The matrix

*Γ*will have the following elements:(3.4)The term

*a*

_{0}in the characteristic equation ‖

*Γ*‖=0 from (3.4) can be factorized into two parts: one for the non-elliptic transonic branches (

*S*

_{a,}

*S*

_{c}) and the other for the elliptic transonic branch (

*S*

_{b}). The condition (3.3) can thus be rewritten as(3.5)where(3.6)after a leading factor

*Δ*sin

^{2}2

*θ*is omitted. The direction of the parabolic point

*θ*along off-symmetry direction will then be determined from vanishing

*R*

_{a}and

*R*

_{b}. The condition

*R*

_{a}=0 suggests that there is at most one parabolic point (

*θ*∈(0,

*π*/4]) on the transonic state branch

*S*

_{a}and it exists only in the crystals with

*Δ*<0. On the transonic state branch

*S*

_{b}, however,

*R*

_{a}=0 implies that the entire branch can have zero transverse curvature if

*a*

^{2}+

*ab*−2(

*b*+1)

^{2}=0.

### (b) Parabolic points in plane

The existence conditions for the parabolic points in the plane can also be reached in a similar fashion as in the (001) plane. In order to avoid complex calculation, we rotate the crystal by *π*/4 about *e*_{x}-axis and keep the earlier definition for the reference plane (** m**,

**) unaltered. The elastic stiffness constants after rotation are given by**

*n**c*

_{22}=

*c*

_{33}=(

*a*+

*b*+2)/2,

*c*

_{23}=(

*a*+

*b*−2)/2 and

*c*

_{44}=(

*a*−

*b*)/2, while the rest remains unchanged. The matrix

*Γ*for the ‘new’ crystal is then defined by(3.7)Note that

*a*

_{0}in the characteristic equation ‖

*Γ*‖=0 can still be factorized into two parts. The condition (3.3) can be expressed as(3.8)with(3.9)where

*f*

_{E}=2

*a*

^{2}+

*ab*−3

*a*−3

*b*

^{2}−9

*b*−4, after a leading factor

*Δ*sin

^{2}

*θ*is omitted. The vanishing resultants,

*T*

_{a}=0 and

*T*

_{b}=0, will define the directions of the parabolic points

*θ*on

*S*

_{a}and

*S*

_{b}, respectively.

To sum up, we have located all the three off-symmetry directions in the cubic crystals admitting the parabolic points given explicitly by (3.6)_{1}, (3.9)_{1} and (3.9)_{2}. The parabolic points defined by (3.6)_{1} and (3.9)_{2} agree with the results obtained from the Gaussian curvature calculation using the perturbation method (see eqn (25) and (21) in Shuvalov & Every 1996). The advantage of the present method is to locate all the parabolic points direct from the characteristic equation for (2.2).

## 4. Swallowtail point: surface normal at the parabolic point

With the parabolic point *θ* clarified, we can now determine its surface normal ** g** as shown in figure 3. Let us construct the reference plane (

**,**

*m***)=(**

*n***,**

*g*

*n*_{2}) (figure 1), or explicitly,

**=cos**

*m**φ*

*e*_{x}+sin

*φ*

*e*_{y}and

**=−sin**

*n**φ*

*e*_{x}+cos

*φ*

*e*_{y}as shown in figure 4.

By adjusting *φ*, ** m** is aligned with the surface normal

**, and by changing the velocity**

*g**v*, we move the line L towards the parabolic point at

*θ*forming a tangent at it. This leads to an ordinary transonic state (similar to the case shown in figure 2

*a*) with

**‖**

*m***and**

*g***‖L. The angle**

*n**φ*will then describe the surface normal

**at the parabolic point, namely, the swallowtail point on the wave surface and in the phonon focusing pattern. Geometrically, from figure 4, one can observe a simple relationship among the directions of the parabolic point**

*g**θ*and its surface normal

*φ*,(4.1)where

*ϕ*is the inclination angle of the wavevector

**from the basis vector**

*k***. Since**

*m**ϕ*defines actually a degenerated Stroh eigenvalue (

*p*=tan

*ϕ*), the relationship between

*θ*,

*φ*and

*p*will be given by

*p*=tan(

*θ*−

*φ*), which in turn can be reformulated as(4.2)The purpose of this reformulation is to make later deductions easier although such a relationship will produce a spurious solution. Fortunately, the spurious solution can be easily identified and hence removed afterward.

Now, we will analyse the degenerated Stroh eigenvalue associated with the transonic state above. With ** m** and

**defined within a symmetry plane, the matrix**

*n**Γ*in (2.2) can always be diagonalized into two blocks and its characteristic equation can be decomposed into two parts ‖

*Γ*‖=

*Q*

_{4}

*Q*

_{2}=0,(4.3)The quartic part

*Q*

_{4}=0 yields eigenvalues

*p*

_{i}(

*i*=1, 2, 3, 4) associated with the non-elliptic slowness branch

*S*

_{a}, and the quadratic part

*Q*

_{2}=0 defines the eigenvalues

*p*

_{i}(

*i*=5, 6) associated with the elliptic slowness branch

*S*

_{b}.

The transonic state defining the surface normal ** g** is in fact an ordinary transonic state of type 1 (figure 2

*a*). Since such a transonic state refers to a duplex-degenerated Stroh eigenvalue, it will require vanishing

*discriminant*of the determinant ‖

*Γ*‖: (i) for a type 1 transonic state on the

*S*

_{a}branch (figure 4), the degenerated Stroh eigenvalues, say,

*p*

_{1}=

*p*

_{2}, requires

*D*[

*Q*

_{4}]=0, where

*D*[

*Q*

_{4}] denotes the discriminant for the polynomial

*Q*

_{4}, or equivalently Res(

*Q*

_{4}, d

*Q*

_{4}/d

*p*) (see, e.g. Bôcher 1936). (ii) For a transonic state on the

*S*

_{b}branch, the degenerated Stroh eigenvalues (

*p*

_{5}=

*p*

_{6}) demands

*D*[

*Q*

_{2}]=0, which simply means .

Up to now, we have a set of conditions defining the relationship between the location of the parabolic point *θ* (3.6) and (3.9), its surface normal *φ* (4.2) and its corresponding transonic state with the degenerated *p*'s (given by *D*[*Q*_{4}]=0 and/or *D*[*Q*_{2}]=0). The next step is to substantiate *D*[*Q*_{4}] and *D*[*Q*_{2}] in the symmetry planes in cubic crystals so that the direction of the surface normal ** g**, defined by

*φ*, can be calculated analytically.

### (a) Swallowtail point in (001) plane

Since ** m** and

**are now defined in the (001) plane, the matrix**

*n**Γ*can be expressed with the following elements:(4.4)Note that the matrix

*Γ*here is identical to the one as eqn (14) in Wang (2008). Here, we will examine the condition for a duplex-degenerated Stroh eigenvalue rather than a triple degenerated one.

In the (001) plane, the parabolic point *θ* is located only on the non-elliptic branch *S*_{a} and is defined by (3.6)_{1}; we can then focus on the quartic part of the characteristic equation. Since it is a type 1 transonic state associated with the parabolic point, it will be subjected to the following three constraints:(4.5)The first constraint defines the parabolic point; the second one ensures that the Stroh eigenvalue *p* is a degenerated one as a transonic state requires (i.e. the line *L* approaches the parabolic point in figure 4*a* by changing *v*); and the last one guarantees that the transonic state does refer to the parabolic point. Here, we choose to calculate the discriminant *D*[*Q*_{4}] with respect to *ρv*^{2} in order to eliminate the variable *ρv*^{2} so that *D*[*Q*_{4}] becomes a function of *p* and *φ*.

By eliminating *θ* and *p* from (4.5), we can get a condition for the direction of the surface normal ** g** described by

*φ*. This is done by calculating the resultants between

*D*[

*Q*

_{4}],

*P*and

*R*

_{a}in (4.5) successively. The final resultant, Res(Res(

*D*,

*P*),

*R*

_{a}), is only dependent on

*φ*. As a consequence of the reformulation of (4.2), the final resultant can always be factorized into two parts, where one of them is spurious. By omitting the spurious part, we can reach a condition for the surface normal direction

*φ*at the parabolic point as a function of the relative elastic constants

*a*and

*b*. The final result, Res(Res(

*D*,

*P*),

*R*

_{a})=0, turns out to be a simple trigonometric equation,(4.6)where

*f*

_{D}=(

*a*−1)(

*a*+

*b*)+

*Δ*(

*b*+1)

^{2}, after some trivial pre-factors are stripped off. Its solution, , describes hence the swallowtail point in the (001) plane. Figure 5 illustrates the variation of as function of

*a*and

*b*. The line

*D*defined by

*f*

_{D}=0 (related to ) marks the critical existence condition for this swallowtail point, which is consistent with the condition

*D*in Every (1981).

Figure 6 shows the focusing patterns for four hypothetical media marked by filled circles in figure 5. One can observe that in the region just below the *D* line, and the swallowtail *points* inward (figure 6*b*). With decreasing *a*, the swallowtail point approaches to the edge of the outer square formed pattern, a cuspidal point, at certain *a* (marked by the curve M in figure 5*a*). When *a* decreases further, the swallowtail point withdraws from the edge, while the swallowtail *points* outward forming inner square-like pattern (figure 6*c*,*d*). Such a transition takes place when the inflection point and the parabolic point coincide and form a monkey saddle point.

One can define a boundary line marking the transition. Since the parabolic point is defined by (3.6)_{1}, and the inflection point can be determined by calculating the zero in-plane curvature *κ* for the slowness branch *S*_{a} ( with *s*=1/*v* defined by *ρv*^{2}=1/2{*a*+1−[(*a*−1)^{2}+(2*a*−*a*^{2}+2*b*+*b*^{2})sin 2*θ*]^{1/2}}) (e.g. Musgrave 1970), condition for the coalesce between these two points (or monkey saddle point) can be deduced and it is given by(4.7)which is represented by the curve M in figure 5.

The curve M divides the figure 5*a* into two regions. Above the curve M, the parabolic point lies near the symmetry axis. Below the curve M, the parabolic point switches with the inflection point and becomes the centre of the outer parabolic line (figure 7), which results in four swallowtails forming an inner square-like pattern. The inner parabolic line consisting of inflection points produces outer circular pattern (or an octagon pattern if *b* is large enough).

Experimentally, a lot of attention was paid to the central square pattern without differentiating its origins. The size of this pattern should be determined by the cuspidal point *φ*_{1} defined in Wang (2008) above the curve M and by the swallowtail point under the curve M. Figure 5*b* combines *φ*_{1} (fig. 3 in Wang 2008) and (figure 5*a*) into one plot where the two sets of results meet exactly along the curve. The earlier work based on simulations (Hurley & Wolfe 1985) agrees with this plot, although it did not explicitly distinguish the difference between the cuspidal points and the swallowtail points.

### (b) Swallowtail points in plane

For the plane, referring to the new crystal defined in §3*b*, the matrix *Γ* can be written with the following elements:(4.8)Its characteristic equation can be written as *Q*_{4}*Q*_{2}=0 with and *Q*_{4}=*Γ*_{33}. Since there are parabolic points in both *S*_{a} and *S*_{b} branches defined by *T*_{a}=0 and *T*_{b}=0, respectively, we will consider them separately (figure 8).

#### (i) The swallowtail point pertaining to *S*_{a}

The derivation of the swallowtail point on the transonic state branch *S*_{a} (figure 8*a*) is similar to the case in the (001) plane because it involves the quartic equation from (4.8). Similar to (4.5), the swallowtail point is determined by(4.9)The final result, Res(Res(*D*[*Q*_{4}], *P*),*T*_{a})=0, determines the swallowtail point on *S*_{a}, , and it can be expressed in terms of a simple trigonometric equation,(4.10)where *c*_{n}(*a*, *b*)'s are four simple but lengthy polynomials of *a* and *b*. The variation of the swallowtail point is illustrated in figure 9*a*.

Concerning the critical existence condition for this caustic point, leads to two limiting curves in the *a*−*b* plot. By setting , we get(4.11)which is represented by the curve A in figure 9*a*. Note that (4.11) is identical with the condition A in Every (1981). When *Δ*=0, we have(4.12)giving (or [111] direction) consistent with the isotropic case. Between the curve A and *Δ*=0, we found that does not yield any physical solution.

#### (ii) The swallowtail point pertaining to *S*_{b}

Along the transonic state branch *S*_{b}, deduction of the swallowtail point *φ* (figure 8*b*) is less complicated. Starting *Q*_{2}=*Γ*_{33} from (4.8), the swallowtail point is determined by(4.13)The final result, Res(Res(*D*[*Q*_{2}], *P*),*T*_{b})=0, can be derived easily and it leaves a very simple condition for the swallowtail point :(4.14)where *f*_{E}=2*a*^{2}+*ab*−3*a*−3*b*^{2}−9*b*−4.

The variation of is illustrated in figure 9*b*. For this swallowtail point, we can reach four critical conditions from (4.14): leads to *a*−*b*=0 and *f*_{E}=0 (see the curve E); leaves *Δ*=0 (consistent with the isotropic case); and gives rise to *f*_{G}=2*a*^{2}−*ab*−5*a*−*b*^{2}−3*b*=0 (see the curve G). We also observe that no physical solution is found in the region between the curve G and *Δ*=0. The curves E and G are also consistent with the conditions derived by Every (1981).

#### (iii) Examples

In order to demonstrate the results for the plane, we calculated the focusing patterns for four hypothetical cubic media with relative elastic constants satisfying *a*+*b*=4.0 with *Δ* varying from −1.5 to 1.5 with increment 1.0 (figure 9). Their phonon focusing patterns are displayed in figure 10 in two columns.

The left column shows the variation of and the right column . The right column indicates that the swallowtail point moves from ST to FT mode as *Δ* increases. The swallowtails (figure 10*b*,*d*) reduce to small diffuse segments in the perimeter of the anticaustic circle in FT mode (figure 10*h*). The left column indicates that the swallowtail point exists only on ST mode when *Δ*>0 (figure 10*e*,*g*). In ST mode, figure 10*b*,*d*,*e*,*g* indicate that the swallowtail point is determined by and in the *Δ*<0 and *Δ*>0 regions, respectively. Compared to the earlier works based on simulations (Hurley & Wolfe 1985), our result gives a rigorous description of the swallowtail points.

Finally, we will conclude our discussion by revisiting two examples used in the earlier investigation on cuspidal points (Wang 2008). Figure 11 shows locations of the swallowtail points in the phonon focusing pattern of CaF_{2} and GaAS. In CaF_{2}, two swallowtail points in the plane are calculated analytically and they agree with the simulations. In GaAs (which lies above the curve M in figure 5), it can be shown that in the (001) plane the inflection and the parabolic points are located along *θ*=12.19° and *θ*=11.78°, respectively. Their corresponding caustic points, a cuspidal and a swallowtail point, are along *φ*_{1}=18.07° and , respectively, which are very close to each other. With the results presented in figure 11 and Wang (2008, figs. 6 and 7), we have determined all the characteristic points on the wave surface and phonon focusing patterns in CaF_{2} and GaAS explicitly, and therefore in cubic crystals in general.

## 5. Concluding remarks

Although many aspects of the slowness/wave surface geometry have been well recognized and qualitatively understood, an effective method to determine the characteristic points, such as inflection/parabolic points on the slowness surface and cuspidal/swallowtail points on the wave surface, remain lacking. The closed-form solution based on the Christoffel equation (Every 1981) provided an analytical scheme to determine the wave surface in general, and it yielded a series of critical conditions for the cuspidal/swallowtail points. The present scheme, based on the Stroh formalism, is shown to be an effective tool to calculate the locations of the cuspidal/swallowtail points.

The common foundation in treating the cuspidal point (Wang 2008) and the swallowtail point is that, they are associated, highly degenerated Stroh eigenvalues, and can then be studied in terms of the extraordinary transonic states. The difference is that the cuspidal point can be determined by employing the Stroh formalism once, while the swallowtail point requires applying the Stroh formalism twice: first to locate the parabolic point and then to calculate its surface normal. By confining our discussion to the symmetry planes, we can use the strong geometric and algebraic constraints to reach simple and explicit/implicit expressions for the locations of the cuspidal and swallowtail points.

The present scheme, established in this work and the earlier one (Wang 2008), makes it possible to calculate all the singular points on the slowness/wave surface of cubic crystals. Moreover, the explicit results can also be employed in recovering the relative elastic constants from a given phonon focusing pattern efficiently.

## Acknowledgments

The author would like to thank the referees for their critical reading and valuable comments.

## Footnotes

- Received October 28, 2008.
- Accepted December 3, 2008.

- © 2009 The Royal Society