Character of Wave Propagation in Anisotropic Media

With Emphasis on Transverse Anisotropy

David H. von Seggern

The classic relation between stress and strain in an elastic medium involves a rank-4 tensor of 81 coefficients.  Invoking various symmetry relations can lead to as few as two coefficients, as in the limiting case of isotropy.  The equation for wave propagation in anisotropic media was sufficiently complex that pre-computer investigators had to make simplifying assumptions on the tensor in order to evaluate, by hand calculation, even a few points for direction-dependent velocity, slowness, and wavefront.   Mathematica enables one to fully plot these quantities even in the most general case.  This article demonstrates the power of Mathematica to handle previously daunting hand calculations and proceeds to visualize the transverse anisotropic case (5 coefficients), equivalent to hexagonal symmetry in crystals, in some detail.  It also contains a novel approach to investigating arbitrary transverse anisotropy through parameters controlled in the Manipulate function of Mathematica.

Introduction

Within physics, basic analysis of materials assumes that their properties are uniform with position; this is the homogeneous assumption, and many results can be formulated easily and deterministically for homogeneous materials.  Introducing inhomogeneity entails more complexity; and results, unless statistical, are usually dependent on numerical solution.  Again within physics, when a certain property of a material has no dependence on direction, then the material is described as  isotropic with respect to that property.  When the property is directionally dependent, then we describe the material as anisotropic with respect to that property.  Note that having one directionally dependent property does not necessarily imply that other properties are also anisotropic.  Anisotropy is present in most real media and is therefore a concern in the analysis of many properties; for instance, velocity of acoustic or elastic waves, velocity of light waves, rate of heat diffusion, and permeability to fluid transport.  

This article deals with anisotropic, but homogeneous, elastic media with application to mineral materials and to Earth structure.   It is largely inspired by the seminal papers of M. J. P. Musgrave ([1] and [2]; hereafter, Musgrave I and Musgrace II) and by the associated paper of Miller and Musgrave ([3]; hereafter, Musgrace III) on anisotropic elastic media.  Those papers graphically showed the behavior of wavefronts in anisotropic media through what was, at the time, laborious hand calculations.  At nearly the same time, [4] and [5] were also treating the velocity and wavefronts and performing similarly arduous calculations to produce figures representing them.  Musgrave I laid out the basic equations for wave propagation in a general anisotropic media; dealt with the concepts of velocity surface, slowness surface (inverse of velocity surface), and wave surface; and showed the complex relations among rays and wavefronts.  Musgrave II and III followed up the general theory with applications to hexagonal and cubic crystals, respectively.  These homogeneous, yet anisotropic, materials allow for some simplifications in the equations.   In these latter two papers, figures based on very arduous calculations at the time were shown for the velocity, slowness, and wave surfaces.  Most striking are cusps which appear in the wave surface, and these were related mathematically to certain characteristics of the slowness surface.  

We assume a Cartesian coordinate system where (x, y, z) = (1, 2, 3) for subscripts denoting direction when defining the algebraic functions representing elastic stress and strain and the elastic waves arising from changes in these quantities.  We will work in the polar and spherical coordinate systems for 2-D and 3-D respectively when employing Mathematica functions to plot figures.  

Letting σ denote stress and τ denote strain, one starts with the general stress-strain relation for an elastic medium [6]:

	(Subscript[σ, ij])= (Subscript[c, ijkl])(Subscript[τ, kl])

The elastic coefficients c are a rank-4 tensor of 81 elements, often called the stiffness tensor.  Using symmetry conditions on (Subscript[c, ijkl]) stated in [7],  the 81-coefficient general tensor can be reduced to a 6 x 6 stiffness coefficient matrix expressed as (Subscript[c, pq]) and the equation simplifies thus:

	(Subscript[σ, p]) = (Subscript[c, pq])(Subscript[τ, q]) 

where σ and τ are now simple vectors of length six.  This reduction is often called the Voigt convention or Voigt notation.  The strain components (Subscript[τ, q]) are the six independent components of the strain matrix τ where (Subscript[τ, kl]) = (Subscript[τ, lk]):

	(Subscript[τ, 1]) = (Subscript[τ, 11])
	(Subscript[τ, 2]) = (Subscript[τ, 22])
	(Subscript[τ, 3]) = (Subscript[τ, 33])
	(Subscript[τ, 4]) = (Subscript[τ, 23])
	(Subscript[τ, 5]) = (Subscript[τ, 31])
	(Subscript[τ, 6]) = (Subscript[τ, 12])

and similarly for (Subscript[σ, p]).  The 6-by-6 coefficient matrix (Subscript[c, pq]) is symmetric; and, due to symmetry conditions again stated in [7], there are only 21 independent elastic coefficients for the general elastic medium.  The most general anisotropic medium would be described by these 21 coefficients, but traditionally authors have studied only media with additional symmetry so that the mathematics is tractable.  

We follow Musgrave I in defining the most general elastic case, using a Cartesian coordinate system (x, y, z).  A general elastic medium supports three types of waves: compressional (P), longitudinal shear (SV), and transverse shear (SH).  We first seek the velocity surface, which is that surface defined by the velocity in the (l, m, n) direction for a given plane wave, where (l, m, n) are the direction cosines comprising the unit vector.   We define a “ray” as along the (l, m, n) direction, and clearly the velocity will vary with the ray direction.  We define a “velocity surface” to be formed by the radii equal to velocity along the (l, m, n) directions.   The velocity surface is defined by solutions to equation 30 in Musgrave I, which has three solutions for the three waves.  P waves have their displacement along a ray in a vertical plane.  SV waves have their displacements in a vertical plane coincident with that containing the P-wave displacements while SH waves have their displacement in the plane normal to that containing P and SV displacements and containing the ray.   We begin by algebraically solving the determinant equation 30 of Musgrave I for ϵ:

	(Det[({{A-e, H, G}, {H, B-e, F}, {G, F, C-e}})]) = 0

using the three parameters A, B, C and substituting parameters F, G, H for his βγ, γα, αβ, respectively.  These parameters will be defined below.  One obtains a cubic equation in e = ρ(v^2) where ρ is density and v is the phase velocity.   The quantity e has units of energy density, or energy per unit volume; later we must remove ρ and take the square root to obtain the phase velocity v.

The cubic equation resulting from the determinant operation is

	A B C-A (F^2)-B (G^2)+2 F G H-C (H^2)+(-A B-A C-B C+(F^2)+(G^2)+(H^2))e+(A+B+C) (e^2)-(e^3)
	
We solve the cubic equation in terms of the parameters (A, B, C, F, G, H).  This is done in the initialization cell, and the solutions for ρ(v^2) for P, SV, and SH waves can be examined by removing the comment brackets.

Solutions for cubic equations can take one of three forms:

	1) three real roots
	2) a single “triple point” real root
	3) a real root with two complex roots, being complex conjugates of one another
	
In the case of an isotropic medium, all three roots are real as for form 1 and quite simple as shown below.  Form 2 will not occur for real elastic media.  The roots for anisotropic cases will generally take form 3, with the latter two roots being complex conjugates.  The first root for the P wave is always real because only even powers of the constants are contained in the radical.  The latter two solutions are the SV and SH waves.  In this cell, Mathematica returns the roots with the P first; the second and third are SV and SH, or possibly reversed, depending on the exact coefficients.   This order cannot be assumed for generally anisotropic media when the exact type of anisotropy is specified in the coefficient matrix.  Thus, one must look graphically at the solutions before assigning them to P, SV, or SH.  

The three solutions are for the P, SV, and SH waves.  In evaluating certain specific cases below, the three solutions were arranged to correspond to the three waves using graphic information developed after the solutions were obtained in Mathematica.  We have not yet introduced the full complexity contained in the parameters (A, B, C, F, G, H).   It is useful to continue the treatment by defining these six parameters in terms of the 21 independent elastic coefficients as done in Musgrave I (eq. 24-29).

	A = (Subscript[c, 11])(l^2)+2 (Subscript[c, 16]) l m+(Subscript[c, 66]) (m^2)+2 (Subscript[c, 15])l n+2(Subscript[c, 56])m n+(Subscript[c, 55]) (n^2)
	B = (Subscript[c, 66]) (l^2)+2 (Subscript[c, 26]) l m+(Subscript[c, 22]) (m^2)+2 (Subscript[c, 46]) l n+2 (Subscript[c, 24]) m n+(Subscript[c, 44]) (n^2)
	C = (Subscript[c, 55]) (l^2)+2 (Subscript[c, 45])l m+(Subscript[c, 44]) (m^2)+2 (Subscript[c, 35]) l n+2 (Subscript[c, 34]) m n+(Subscript[c, 33]) (n^2)
	F = (Subscript[c, 56]) (l^2)+((Subscript[c, 25])+(Subscript[c, 46])) l m+(Subscript[c, 24]) (m^2)+((Subscript[c, 36])+(Subscript[c, 45])) l n+((Subscript[c, 23])+(Subscript[c, 44])) m n+(Subscript[c, 34]) (n^2)
	G = (Subscript[c, 15]) (l^2)+((Subscript[c, 14])+(Subscript[c, 56])) l m+(Subscript[c, 46]) (m^2)+((Subscript[c, 13])+(Subscript[c, 55]) l n+((Subscript[c, 36])+(Subscript[c, 45])) m n+(Subscript[c, 35]) (n^2)
	H = (Subscript[c, 16]) (l^2)+((Subscript[c, 12])+(Subscript[c, 66])) l m+(Subscript[c, 26]) (m^2)+((Subscript[c, 14])+(Subscript[c, 56])) l n+((Subscript[c, 25])+(Subscript[c, 46])) m n+(Subscript[c, 45])(n^2)

One gets a feeling for the complexity of the solutions by substituting these expansions in the cubic equation above.  Different types of anisotropy can then be modeled by setting values for different sets of coefficients, taking into account symmetries and relations in the different types.  Substitution of these expanded parameters into the solutions for ρ(v^2) creates very long expressions, containing hundreds of terms.  The following cell will produce these lengthy solutions for P, SV, and SH waves.   It must be evaluated to be able to do further calculations below.

Unprotect[C];
Clear[A,B,C,F,G,H];
equation=Collect[Expand[Det[{{A-e,H,G},{H,B-e,F},{G,F,C-e}}]],e];
solution=Solve[equation==0,e];
rhov2P=e/. solution[[2]];
rhov2SV=e/. solution[[3]];
rhov2SH=e/. solution[[1]];
awave={"P","SV","SH"};
Clear[c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66];
Clear[l,m,n];
A=c11 (l^2)+2 c15 l n+2 c16 l m+c55 (n^2)+2 c56 m n+c66 (m^2);
B=c22 (m^2)+2 c24 m n+2 c26 l m+c44 (n^2)+2 c46 l n+c66 (l^2);
C=c33 (n^2)+2 c34 m n+2 c35 l n+c44 (m^2)+2 c45 l m+c55 (l^2);
F=m n (c23+c44)+c24 (m^2)+l m (c25+c46)+c34 (n^2)+l n (c36+c45)+c56 (l^2);
G=l n (c13+c55)+l m (c14+c56)+c15 (l^2)+c35 (n^2)+m n (c36+c45)+c46 (m^2);
H=l m (c12+c66)+l n (c14+c56)+c16 (l^2)+m n (c25+c46)+c26 (m^2)+c45 (n^2);

One can show the lengthy solutions, which take up several pages, by evaluating the following cell.

Print["P-wave general solution:\n",rhov2P];
Print["SV-wave general solution:\n",rhov2SV];
Print["SH-wave general solution:\n",rhov2SH];

Visualization of the Completely General Case

With a completely general expression for the velocity (actually ρ(v^2)), a Manipulate function can be set up to allow manipulation of all 21 elastic constants, plus density, while plotting the resultant velocity in planar cross-sections.  The evaluation takes several tens of seconds for any change in elastic constants and is not provided here.   It is included as an appendix with this article.  When performed, this can show the capability of Mathematica to visualize extremely complex equations though.   

Hand computation for any setting of the parameters, except for specific cases of anisotropy as treated below, would be nearly impossible, thus requiring earlier investigators to treat only simplified cases of anisotropy, and then only in planes intersecting the 3-D solutions, as described below.  Musgrave II and III treated just two tractable cases of hexagonal and cubic symmetry in crystals.  The expressions involving the parameters (A, B, C, F, G, H) are significantly simplified in those papers, but still remain daunting to evaluate without computer aid for given elastic coefficients of the crystals.  We will proceed to look at such simplified cases.  One of our goals is to check the calculations in Musgrave II and III.

Types of Symmetry

In this section, we provide a means to impose further symmetries and simplifications in the general 21-coefficient equation above.    Let us describe the simplest of cases before turning to the cases treated in Musgrave II and III.  This would be for an isotropic medium where the velocity is independent of direction.  The 19 rules shown in Table 1  are from [7], (pp. 160-161); they leave only two coefficients: (Subscript[c, 11]) and (Subscript[c, 12]).

({{, rule}, {1, (Subscript[c, 13])->(Subscript[c, 12])}, {2, (Subscript[c, 14])->0}, {3, (Subscript[c, 15])->0}, {4, (Subscript[c, 16])->0}, {5, (Subscript[c, 22])->(Subscript[c, 11])}, {6, (Subscript[c, 23])->(Subscript[c, 12])}, {7, (Subscript[c, 24])->0}, {8, (Subscript[c, 25])->0}, {9, (Subscript[c, 26])->0}, {10, (Subscript[c, 33])->(Subscript[c, 11])}, {11, (Subscript[c, 34])->0}, {12, (Subscript[c, 35])->0}, {13, (Subscript[c, 36])->0}, {14, (Subscript[c, 44])->(((Subscript[c, 11])-(Subscript[c, 12]))/2)}, {15, (Subscript[c, 45])->0}, {16, (Subscript[c, 46])->0}, {17, (Subscript[c, 55])->(((Subscript[c, 11])-(Subscript[c, 12]))/2)}, {18, (Subscript[c, 56])->0}, {19, (Subscript[c, 66])->(((Subscript[c, 11])-(Subscript[c, 12]))/2)}})

This solution is given in terms of the two remaining elastic coefficients, and there is no directional dependence of velocity.  This corresponds to the “three real roots” form of solution to a cubic equation.  However, two of the roots, for SV and SH, are equal and there are really only two propagation velocities.  The SV and SH motions are orthogonal and each orthogonal to that of the P motion.  Thus, along any ray, one has local displacements along three orthogonal directions, the P motion along the ray and the two S motions perpendicular to it.  It is common in seismology to use (Subscript[c, 12]) = λ and (Subscript[c, 11]) = λ + 2μ so that the P velocity becomes (Sqrt[(λ+2μ)/ρ]) and the S velocity (SV and SH both) becomes (Sqrt[μ/ρ]) where ρ is the density of the medium.  λ and μ are called Lame’s constant and rigidity, respectively, in the seismological literature.  An isotropic medium can have a range of behavior depending on the values of the two constants, and this behavior is encapsulated in the parameter called Poisson’s ratio ν where

	ν = (λ/(2(λ + μ)))
	
The value of ν is in the range [0, 1/2] where the lower bound is obtained when λ = 0 and the upper bound is obtained when μ = 0 (liquid).  The special case of ν = 1/4 (λ = μ) is given the name Poisson solid, and much of the solid Earth has a value of ν near that of a Poisson solid.  In this case, the ratio of P to S velocities is  (Sqrt[3]).  The graphical illustrations of velocity, slowness, and wave surfaces will be trivial for the isotropic case, being spherical surfaces for all three waves.  The only difference in surfaces is caused by varying (Subscript[c, 11]) and (Subscript[c, 12]) (equivalently, λ and μ) to change the velocities of P and S waves.

We don’t treat the isotropic case any further due to its simplicity.  We proceed now to the case of hexagonal symmetry in Musgrave II.  Hexagonal symmetry is common for many crystals and demonstrates anisotropy.  The text cell below shows the rules that reduce the number of independent elastic coefficients to five ((Subscript[c, 11]), (Subscript[c, 12]), (Subscript[c, 13]), (Subscript[c, 33]), (Subscript[c, 44])) in the case of hexagonal symmetry.  Evaluating the input cell that follows will set up expressions for ρ(v^2), to be used later in graphics displays.  Other symmetry cases can be studied; for instance, cubic symmetry, treated in Musgrave II, needs only three coefficients ((Subscript[c, 11]), (Subscript[c, 12]), (Subscript[c, 44])).  We treat only hexagonal symmetry in detail here because it leads naturally into transverse anisotropy common in much of the real Earth.  

Consider a transversely anisotropic medium in which the velocity in the horizontal direction (x-y plane) differs from that in the vertical (z) direction.    This is fulfilled by crystals of hexagonal symmetry (Musgrave II) and is affirmed in several other studies [8, 9, 10].  For this type of anisotropy, the z axis is an axis of symmetry with the five coefficients determining the degree of transverse anisotropy.  The 16 rules to convert the general anisotropic case to the hexagonal symmetry case are given in Table 2.   Transverse anisotropy is a common model applied to the Earth.  This type of anisotropy naturally represents much of the Earth’s interior due to the mode of deposition for sediments, non-hydrostatic stresses, layering, fracturing systems with preferred orientations, or some combination of all these.  Under some circumstances, transverse isotropy may become azimuthally dependent in the x-y plane, as for spreading ocean floors [11], leading to a slightly more complex stiffness matrix.

({{, rule}, {1, (Subscript[c, 14])->0}, {2, (Subscript[c, 15])->0}, {3, (Subscript[c, 16])->0}, {4, (Subscript[c, 22])->(Subscript[c, 11])}, {5, (Subscript[c, 23])->(Subscript[c, 13])}, {6, (Subscript[c, 24])->0}, {7, (Subscript[c, 25])->0}, {8, (Subscript[c, 26])->0}, {9, (Subscript[c, 34])->0}, {10, (Subscript[c, 35])->0}, {11, (Subscript[c, 36])->0}, {12, (Subscript[c, 45])->0}, {13, (Subscript[c, 46])->0}, {14, (Subscript[c, 55])->(Subscript[c, 44])}, {15, (Subscript[c, 56])->0}, {16, (Subscript[c, 66])->(((Subscript[c, 11])-(Subscript[c, 12]))/2)}})

Next consider a medium in which the velocity in the three axial directions is identical but varies off-axis in a way such that the velocity surface is point symmetric about (0, 0, 0).  This is the case for a cubic crystal (Musgrave III).  The 18 simplifying rules are given in Table 3.  The velocity surface is not a rotational surface although the three principle planes are all planes of reflection for the surface.  The curve formed by cutting the surface with the x-y plane reflects at the x or y axis but, contrary to hexagonal symmetry, is not a circle.  We will not examine the cubic symmetry case in this article; but by making the substitutions shown in the table above, one can produce solutions such as illustrated in Musgrave III when the constants given in that paper are used.  Cubic symmetry is displayed by certain crystals but is not a form of anisotropy contained in the real earth on a macro scale.

({{, rule}, {1, (Subscript[c, 13])->(Subscript[c, 12])}, {2, (Subscript[c, 14])->0}, {3, (Subscript[c, 15])->0}, {4, (Subscript[c, 16])->0}, {5, (Subscript[c, 22])->(Subscript[c, 11])}, {6, (Subscript[c, 23])->(Subscript[c, 12])}, {7, (Subscript[c, 24])->0}, {8, (Subscript[c, 25])->0}, {9, (Subscript[c, 26])->0}, {10, (Subscript[c, 33])->(Subscript[c, 11])}, {11, (Subscript[c, 34])->0}, {12, (Subscript[c, 35])->0}, {13, (Subscript[c, 36])->0}, {14, (Subscript[c, 45])->0}, {15, (Subscript[c, 46])->0}, {16, (Subscript[c, 55])->(Subscript[c, 44])}, {17, (Subscript[c, 56])->0}, {18, (Subscript[c, 66])->(Subscript[c, 44])}})

Rules for additional types of anisotropy exemplified by crystals are given in [7] (pp. 159-160).  The cubic and other types are also not considered in this notebook because the hexagonal type demonstrates adequately the nature of anisotropy in the Earth.  It also demonstrates the ease with which Mathematica can compute and display meaningful results.  To be precise, the cases we investigate should be referred to as transverse anisotropy with a vertical axis of symmetry, but we assume the symmetry about the vertical axis henceforth.  

Set the replacement rules for elastic coefficients in the hexagonal-symmetry case and simplify the solutions in this cell.  Solutions for ρ(v^2) are returned in terms of the remaining five coefficients.   Compare these relatively short expressions to those obtained for the general case above.  Note there is no dependence on l or m for hexagonal symmetry, and so the solutions are azimuthally independent.

symmetry="hexagonal";
rules={c14->0,c15->0,c16->0,c22->c11,c23->c13,c24->0,c25->0,c26->0,c34->0,c35->0,c36->0,c45->0,c46->0,c55->c44,c56->0,c66->((c11-c12)/2)};
AA=A/. rules;BB=B/. rules;CC=C/. rules;FF=F/. rules;GG=G/. rules;HH=H/. rules;
equation=Det[{{AA-e,HH,GG},{HH,BB-e,FF},{GG,FF,CC-e}}];
solution=Solve[equation==0,e];
eSH=Simplify[e/. {solution[[1]]},(l^2)+(m^2)+(n^2)==1][[1]];
eSV=Simplify[e/. {solution[[2]]},(l^2)+(m^2)+(n^2)==1][[1]];
eP=Simplify[e/. {solution[[3]]},(l^2)+(m^2)+(n^2)==1][[1]];
Print["((ρv)^2) SOLUTIONS for symmetry = ",symmetry];
Print["P:",eP];
Print["SV:",eSV];
Print["SH:",eSH];

Now we define functions to express the three velocities for hexagonal crystals using the three solutions just obtained.  Similar functions could be set up for any symmetry case.  These functions allow specification of actual values for remaining coefficients after symmetry conditions such as Table 2 or 3 are applied.

vPhex[c11in_,c12in_,c13in_,c33in_,c44in_,ρin_]:=Sqrt[eP/ρ/.{c11->c11in,c12->c12in,c13->c13in,c33->c33in,c44->c44in,ρ->ρin}];
vSVhex[c11in_,c12in_,c13in_,c33in_,c44in_,ρin_]:=Sqrt[eSV/ρ/.{c11->c11in,c12->c12in,c13->c13in,c33->c33in,c44->c44in,ρ->ρin}];
vSHhex[c11in_,c12in_,c13in_,c33in_,c44in_,ρin_]:=Sqrt[eSH/ρ/.{c11->c11in,c12->c12in,c13->c13in,c33->c33in,c44->c44in,ρ->ρin}];

Visualization of Velocity for Media With Specified Coefficients

It can be seen by examining the hexagonal-symmetry solutions that ρ(v^2) for the three waves take the third form given above for solutions of a cubic equation -- one real root and two, possibly, complex-conjugate roots.  The SH wave is the real root and does not include either (Subscript[c, 13]) or (Subscript[c, 33]).  This is a consequence of the particle motion always being normal to the plane of the ray.  Furthermore, velocity is independent of l and m, thus circular about the z axis.    Note that both P and SV are independent of (Subscript[c, 12]).  This is a consequence of the particle motion being in the vertical plane holding the ray.  The solutions for P and SV waves are paired and would be complex conjugates if the expression under the radical sign is negative.  However, the fact that the expression ρ(v^2) is intrinsically real appears to dictate that the expression under the radical is always positive.  That it is always positive has been verified by running a Monte Carlo algorithm wherein the 5 elastic constants are drawn from [0, 1] ranges while the unit vector (l, m, n) is varied randomly over a sphere.   No cases were found where the expression was negative.   This does not guarantee that the full expressions for ρ(v^2) for P and SV are positive and real though.  This was checked by repeating the Monte Carlo algorithm on the full expressions.  In the case of P, ρ(v^2)was shown to always be positive real, but SV had some cases where ρ(v^2) was negative real.  We dismiss such combinations of constants as physically improbable in actual materials.  Similarly, the structure of the expression for ρ(v^2) for the SH wave can theoretically admit cases where that expression would be negative; such cases are again physically improbable.

The transversely anisotropic media to be considered require five elastic coefficients to be specified.   Specific cases are simply evaluated by entering the desired elastic coefficients in the cell below.  We are working in MKS units while the reference frame of Musgrave II was in cgs units.  Musgrave’s elastic coefficients need to be multiplied by (10^4) and the density by (10^3).  Thus, our working units are MPa for the elastic coefficients and kg/(m^3) for density. 

Use the next cell to evaluate algebraically the expressions for a specific transversely anisotropic medium and then graphically plot the velocity surface for one of the chosen waves (P, SV, or SH).  The constants are entered in MKS units, and they are equivalent to stress in units of MPa.  Density ρ is entered as kg/(m^3).  The velocity then is obtained in units of km/s.  Here are some typical value ranges:

100,000 ≤ c11 ≤ 300,000 MPa
10,000 ≤ c12 ≤ 150,000 MPa
10,000 ≤ c13 ≤ 100,000 MPa
50,000 ≤ c33 ≤ 300,000 MPa
25,000 ≤ c44 ≤ 150,000 MPa
2,000 ≤ ρ ≤ 4,000 kg/(m^3)

Two important sets of parameters are those used in Musgrave II for the hexagonal crystals beryl and zinc.  Note that there has been a change in values to convert to MKS units of MPa and kg/(m^3).

c11		c12		c13		c33		c44		ρ
	beryl		269,400	96,100		66,100		236,300	65,300		2,700
	zinc		143,000	17,000		33,000	  	50,000		40,000		7,100

Coefficients for cubic crystals are given in Musgrave III.  Both hexagonal and cubic configurations are studied rigorously in [12].  All that remains to do now is to specify the actual values for the remaining coefficients, whose number depends on the type of symmetry chosen.  Figure 1 shows the result of specifying one set of exact coefficients.  The following cell should be evaluated in order to enter values for the elastic constants, according to the symmetry specified.  Dialog boxes will request numerical values for the elastic constants.  After the constants are entered, a Manipulate function will allow viewing of the P, SV, and SH velocity surfaces in 3-D.  Note that an isotropic medium can be specified while setting the hexagonal symmetry coefficients by first setting (Subscript[c, 11]) and (Subscript[c, 12]) and further applying

		(Subscript[c, 13]) = (Subscript[c, 12])
		(Subscript[c, 33]) = (Subscript[c, 11])
		(Subscript[c, 44])= (((Subscript[c, 11])- (Subscript[c, 12]) )/2)
	
Enter values for the coefficients into the functions defined above for hexagonal symmetry (vPhex, etc.) and compute the 3-D velocity surfaces (vP, etc.) in the following cell.  Set the hexagonal-symmetry coefficients arbitrarily, but according to guidelines in Table 4.  Note that a truer surface can be plotted by setting PlotPoints to 50 or more, and MaxRecursion to 2 or more but that the plotting time increases as the square of PlotPoints and linearly as MaxRecursion.

c11in=270000;
c12in=96000;
c13in=66000; 
c33in=236000;
c44in=65000;
ρin=2700.;
vP=vPhex[c11in,c12in,c13in,c33in,c44in,ρin];
vSV=vSVhex[c11in,c12in,c13in,c33in,c44in,ρin];
vSH=vSHhex[c11in,c12in,c13in,c33in,c44in,ρin];
coefLabel= "\nc11 = "<>ToString[c11in]<>", c12 = "<>ToString[c12in]<>", c13 = "<>ToString[c13in]<>"\nc33 = "<>ToString[c33in]<>", c44 = "<>ToString[c44in]<>", ρ = "<>ToString[ρin];
(*Using polar angle θ and azimuthal angle ϕ in SphericalPlot3D.*)
Clear[θ,ϕ];
Manipulate[
temp =ToExpression["v"<>wave]/. {l->Cos[ϕ]Sin[θ],m->Sin[ϕ]Sin[θ],n -> Cos[θ]};
SphericalPlot3D[temp,{θ,0,π},{ϕ,0,2π},AxesLabel->{X,Y,Z},PlotRange->Automatic,
PlotPoints->Automatic, MaxRecursion->Automatic,
PlotLabel->"velocity for "<>symmetry<>" symmetry -- " <>wave<>coefLabel,BaseStyle->12,LabelStyle->12,ImageSize->400],
{wave,awave}  (*Choose a wave -- P, SV, or SH.*),
TrackedSymbols->{wave},
SaveDefinitions->True
]

Slowness Surface

We now introduce the slowness surface, which is just the geometric inverse of the velocity surface with respect to the unit sphere such that the slowness s in terms of the unit vector is s(l, m, n) = 1/v(l, m, n).  The shape of the slowness surface will be critical in predicting the shape of the wave surface, which tracks the wavefront of constant phase.  The mathematical importance of the slowness surface is discussed in [13] in a clear and concise manner for anisotropic media, with emphasis on transverse anisotropy.   

The velocity and slowness in cross-sections containing the z axis and in the x-y plane can be plotted in the subsequent cells.  All three waves (P, SV, SH) are treated together, provided the velocity-squared expressions (v2P, v2SV, v2SH) have had values for coefficients assigned in the previous input cell.  Figure 2 shows typical results of this evaluation for vertical cross-sections.  An additional cell can be evaluated to graph the sections in the equatorial plane.  

Note that the plots shown in Figure 2 agree with those in Musgrave II when Table 5 coefficient values are used.  In Musgrave II, due to the arduous nature of hand calculations, only a few points were calculated by that author in order to define the curves.  The scales on the x and y axes for the slowness surfaces are in error in Musgrave II, and are shown here correctly.  The next two cells allow one to plot the velocity and slowness in vertical cross-section or on the equatorial plane.

Evaluate the following cell to plot the velocity and slowness in vertical cross-section using 3-D velocity expressions developed above (vP, etc).  The angle ψ is the azimuthal angle of the vertical plane, counterclockwise from the positive x axis.  For hexagonal symmetry (commonly call transverse anisotropy), the profile will not change with ψ while it may for other crystal symmetries.   Here we use the Cartesian angle ϑ, counterclockwise from the horizontal plane.

aψ =0;
ψx = Cos[aψ];
ψy = Sin[aψ];
Clear[ϑ];

vPz = Simplify[vP/.{ l ->ψx*Cos[ϑ],m->ψy*Cos[ϑ],n ->Sin[ϑ]}];
sPz= 1/vPz;
vSVz=Simplify[vSV/.{ l ->ψx*Cos[ϑ],m->ψy*Cos[ϑ],n ->Sin[ϑ]}];
sSVz=1/vSVz;
vSHz=Simplify[vSH/.{ l ->ψx*Cos[ϑ],m->ψy*Cos[ϑ],n ->Sin[ϑ]}];
sSHz=1/vSHz;

p1=PolarPlot[{vPz,vSVz,vSHz},{ϑ, 0, 2π},PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->12,LabelStyle->12,PlotLabel->"velocity in vertical plane\n"<>symmetry<>" symmetry"<>coefLabel,AspectRatio->Automatic,ImageSize->300];
p2=PolarPlot[{sPz,sSVz,sSHz},{ϑ, 0, 2π},PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->12,LabelStyle->12,PlotLabel->"slowness in vertical plane\n"<>symmetry<>" symmetry"<>coefLabel,AspectRatio->Automatic,ImageSize->300];
GraphicsRow[{p1,p2},ImageSize->{750,500}]

By setting θ = 0, we can see the horizontal cross-section of the velocity and slowness surfaces.  They will, of course, be circles for hexagonal symmetry.  Evaluate this cell to plot those projections.  We use the variable ψ to represent the CCW angle within the x-y plane.

Clear[ψ];

vPh=Simplify[vP/.{ l ->Cos[ψ],m->Sin[ψ],n ->0}];
sPh= 1/vPh;
vSVh=Simplify[vSV/.{ l ->Cos[ψ],m->Sin[ψ],n ->0}];
sSVh=1/vSVh;
vSHh=Simplify[vSH/.{ l ->Cos[ψ],m->Sin[ψ],n ->0}];
sSHh=1/vSHh;

p1=PolarPlot[{vPh,vSVh,vSHh},{ψ,0,2π},PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->12,LabelStyle->12,
AspectRatio->Automatic,PlotLabel->"velocity in equatorial plane\n"<>symmetry<>" symmetry"<>coefLabel,ImageSize->300];
p2=PolarPlot[{sPh,sSVh,sSHh},{ψ,0,2π},PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->12,LabelStyle->12,
AspectRatio->Automatic,PlotLabel->"slowness in equatorial plane\n"<>symmetry<>" symmetry"<>coefLabel,ImageSize->300];
GraphicsRow[{p1,p2},ImageSize->{800,700}]

Wave Surface

We now introduce the wave surface, also called the group-velocity surface or wavefront.  It is a snapshot in time of the leading edge of the propagating wavefront.  Again, the reader is referred to the clear discussion in [13] on the relation of the three surfaces treated here.  For waves in transversely anisotropic media, the slowness and wave surfaces are polar reciprocals of one another, with the z axis being an axis of symmetry.

The wave surface is defined as equal phase of the plane waves at some finite, constant time T when each plane wave at (l, m, n) emanates through (0, 0, 0) at t = 0.  Physically, it shows how the energy propagates in the medium.  For convenience, let T = 1 second.  Equivalently, one can assume a point source at (0, 0, 0) and draw the wavefront at T = 1.  It is pointed out in Musgrave II and III that the wave surface is dependent on the curvature of the slowness surface.  When this curvature changes sign (convex to concave or vice-versa), a cusp is generated in the wave surface.  These are illustrated in Musgrave II and III for the media studied there.  The cusps occur in pairs, and the midpoint between them is where the slowness surface reaches a relative minimum, with small curvature.  The shape of the wave surface can be predicted by plotting the curvature of the slowness surface over one quadrant.  Curvature r is defined thus

	r = ((((Subscript[r, 0]^2) + 2(Subscript[r, 1]^2) - (Subscript[r, 0])(Subscript[r, 2])))/((((Subscript[r, 0]^2) + (Subscript[r, 1]^2)))^(3/2)))
	
where
	(Subscript[r, 0])= curve itself,
	(Subscript[r, 1])= 1st derivative of the curve wrt angle θ,
	(Subscript[r, 2])= 2nd derivative of the curve wrt angle θ.
	
This is an extremely lengthy expression but is computed in Mathematica quickly and can be plotted as in the next cell for all three waves.  It plots the curvature r of the slowness curve as a function of θ, the polar angle.  Curvature r = 1/ρ where ρ is the radius of curvature, and curvature = 0 when radius of curvature = ∞.  Where the curvature changes sign indicates that the wave surface will have a cusp (they occur in pairs) with triplication of the wave arrival in time between the two cusps.  The angle where the curvature reaches a minimum will indicate the mid-point between cusps of the wave surface.  For an isotropic medium, the fact that velocity is constant means the curvature is constant.  

Evaluate the following cell to graph the curvature of the slowness vertical profile for the three waves.  We plot the curvature for only one quadrant because the graph is the same in all 4 quadrants.

Clear[η];
rwave={};
Do[
r0 =ToExpression["s"<>awave[[i]]<>"z"]/.ϑ->π/2 - θ;  (*Change back to polar angle.*)
r1=D[r0,θ];
r2=D[r0,{θ,2}];
r= (r0^2 + 2*r1^2 - r0*r2)/(r0^2 + r1^2)^(3/2);  (*This is the curvature, not radius of curvature.*)
AppendTo[rwave, r] 
,{i,1,3}
];
rwavedegree=rwave/.θ->η Degree;

(*PlotLegends does not work here unless Evaluate is wrapped around the expression.*)
Plot[Evaluate[rwavedegree],{η, 0, 90},PlotRange->{{0,90},All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,
MaxRecursion->6,Frame->True,FrameStyle->Thick,FrameTicks->{{Automatic,None},{Table[i*10,{i,0,9}],None}},FrameLabel->{"polar angle (degrees)","curvature ϱ of slowness"},GridLines->{Table[i*10,{i,1,8}],Automatic},GridLinesStyle->{Thin,Thin},BaseStyle->14,LabelStyle->14,PlotRangePadding->{0,0},AspectRatio->1/2,PlotLabel->symmetry<>" symmetry"<>coefLabel,ImageSize->400]

The derivation of the wave surface in Musgrave II or III for hexagonal and cubic symmetry, respectively, is done with some elegant geometry and algebra specific to those media, allowing tractable calculation of the points on the wave surface.  The general case is easily represented algebraically though as in [5] (eq. 1 and 2), where expressions were derived directly from the original work in [4] (eq. 17 and 18).  The result avoids any manipulation or computation using the stiffness coefficients.  Here we treat only the transverse anisotropic case where the wave surface is invariant about the z axis.

Let (OverVector[n]) be the normal to the slowness surface at angle θ.  Then the point on the wave surface that corresponds to (OverVector[n])(θ) is parallel to it and at a time T = 1 second.  The surface can be scaled by a constant to represent it at different times.  Let (OverVector[u]) be the vector to the wave surface while (OverVector[v]) is the vector to the velocity surface, and let (Subscript[u, h]) and (Subscript[u, v]) be the horizontal and vertical components of the wave surface.  Then it can be shown, taking the formulation in [5], that by changing from polar angle θ to the angle ϑ counterclockwise from the x-y plane, we have these vector components of (OverVector[u]):

		(Subscript[u, h]) = |v|cos(ϑ) - ((d|v|)/(d ϑ))sin(ϑ)
		(Subscript[u, z]) = |v|sin(ϑ) + ((d|v|)/(d ϑ))cos(ϑ)	
	
The form of these expressions indicates that the wave surface is the velocity surface with added quantities which relate to the rate of change of the velocity surface in the angular direction.  These quantities, (Subscript[u, h]) and  (Subscript[u, z]), can easily be evaluated for small increments of ϑ to provide a list of points to plot in joined fashion in order to graph the projection of the wave surface on a vertical plane holding the z axis as in the next cell.  Some error messages may be generated, but are inconsequential regarding the plotted result.  

Evaluate this cell to plot the wave surfaces in vertical cross-section through the z axis at a time of 1 second.  Depending on the elastic coefficients, there may be cusps in the SV wave profile.

(*Use the Cartesion polar angle ϑ (CCW from the horizontal plane) for the angle in the evaluation.*)
(*Note that the exact azimuthal angle aψ was set above in plotting the velocity and slowness in vertical plane.*)
inc=π/1800;  (*A small inc is required to define cusps well.*)
angle=Table[aϑ//N,{aϑ,0,2π,inc}];
na=Length[angle];
del={};
waveSurface={};
(*Loop over three waves.*)
Do[
temp = ToExpression["v"<>awave[[i]]<>"z"];
dvdϑ=0;
dvdϑ=D[temp,ϑ];
uh=temp*Cos[ϑ] - dvdϑ*Sin[ϑ];
uz=temp*Sin[ϑ] + dvdϑ*Cos[ϑ];
hvector=Table[Chop[uh],{ϑ,0.,2π,inc}];
zvector=Table[Chop[uz],{ϑ,0.,2π,inc}];
AppendTo[del,angle - ArcTan[hvector,zvector]];
AppendTo[waveSurface,Transpose[{hvector,zvector}]]
,{i,1,3}
];
ListPlot[waveSurface,Joined->True,PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->12,LabelStyle->12,AspectRatio->Automatic,PlotLabel->symmetry<>" symmetry\nwave surface at 1 s in vertical plane"<>coefLabel,ImageSize->300]

An interesting graph is that shown by Musgrave II for beryl or zinc (Figures 5 and 3, respectively, in that paper).  This is the deviation of the ray from the direction to the slowness surface from (0, 0, 0).  The ray for a given point on the slowness surface is associated with the wave surface at a deviated angle which is derived when computing the wave surface.  Sweeping down from the positive z axis as in Figure 4, a positive deviation followed by a negative deviation is associated with cusps, as in the case of the SV wave surface for the beryl crystal.  

Evaluate the following cell to graph the ray deviations in the wave surface.  For a given angle to the slowness surface, one adds this deviation to perturb the ray direction to the wave surface.

angles=Reverse[Take[angle,1+Floor[na/4]]];
delp={};
Do[
AppendTo[delp,Take[del[[i]],1+Floor[na/4]]];
Reverse[delp[[i]]]
,{i,1,3}
];
delp=Chop[delp];
ListPlot[{Transpose[{angles/Degree,delp[[1]]/Degree}],Transpose[{angles/Degree,delp[[2]]/Degree}],Transpose[{angles/Degree,delp[[3]]/Degree}]},PlotRange->{{0,90},All},Joined->True,PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,Frame->True,FrameStyle->Thick,FrameTicks->{{Automatic,None},{Table[i*10,{i,0,9}],None}},FrameLabel->{"polar angle (degrees)","deviation of ray (degrees)"},GridLines->{Table[i*10,{i,1,8}],Automatic},GridLinesStyle->{Thin,Thin},PlotRangePadding->0,BaseStyle->14,LabelStyle->14,AspectRatio->1/2,PlotLabel->symmetry<>" symmetry"<>coefLabel,ImageSize->500]

When doing anisotropic, hexagonal-symmetry calculations, the results bear some comparison with those done by Musgrave II.  Using the parameter values given in Table 5 above for beryl, the Mathematica evaluation here shows the cusps of the SV wave surface somewhat more accurately.  Also, the ray deviation plot is slightly different than the one in Musgrave II.  These differences are not surprising, given the complexity of the expressions.  For the zinc case, a small difference between the SV-wave surface computed here at the cusps on the horizontal axis, being very tiny, and the larger cusps shown in Musgrave II suggests perhaps a small error crept into the hand calculations.  

A similar approach can be used to obtain the equations for the wave surface in the x-y plane (θ = π/2 or z = 0):

	(Subscript[u, x]) = |v|cos(ψ) - ((d|v|)/(d ψ))sin(ψ)
	(Subscript[u, y]) = |v|sin(ψ) + ((d|v|)/(d ψ))cos(ψ)	
	
where ψ is the counterclockwise angle from the positive x axis in the x-y plane.  The wave surface for the hexagonal symmetry case will, of course, be circular in the x-y plane.  Again, the evaluation may produce inconsequential error messages.  

Evaluate the following cell to graph the wave surfaces in the equatorial plane at a time of 1 second.  For the case of hexagonal symmetry, this horizontal profile of the wave surface is a circle for all three waves.

xvector={};
yvector={};
delz0={};
waveSurfacez0={};
Clear[ψ];
(*Loop over three waves.*)
Do[
temp = ToExpression["v"<>awave[[i]]<>"h"];
dvdψ=0;
dvdψ=D[temp,ψ];
ux=temp*Cos[ψ] - dvdψ*Sin[ψ];
uy=temp*Sin[ψ] + dvdψ*Cos[ψ];
xvector=Table[Chop[ux],{ψ,0.,2π,inc}];
yvector=Table[Chop[uy],{ψ,0.,2π,inc}];
AppendTo[delz0,angle - ArcTan[xvector,yvector]]; (*can use ArcTan on lists*)
AppendTo[waveSurfacez0,Transpose[{xvector,yvector}]]
,{i,1,3}
];
ListPlot[waveSurfacez0,Joined->True,PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->12,LabelStyle->12,AspectRatio->Automatic,PlotLabel->symmetry<>" symmetry\nwave surface at 1 s in equatorial plane"<>coefLabel,ImageSize->300]

Arbitrary Transverse Anisotropy

We will treat arbitrary transverse anisotropy (ATA) in terms of deviations from the isotropic case.  Recall that, for transverse anisotropy (hexagonal symmetry), only five elastic coefficients ((Subscript[c, 11]), (Subscript[c, 12]), (Subscript[c, 13]), (Subscript[c, 33]), (Subscript[c, 44])) remain.  In this case, there is no x-y dependence, and the solution at any polar angle is circular about the z axis.  Three measures of transverse anisotropy can be formed by varying the coefficients (Subscript[c, 13]), (Subscript[c, 33]), and (Subscript[c, 44]) relative to the isotropic ones thus:
   	
	(Subscript[c, 44]) = f·((Subscript[c, 11]) - (Subscript[c, 12]))/2,
	(Subscript[c, 13]) = g·(Subscript[c, 12]),
	(Subscript[c, 33]) = h·(Subscript[c, 11]).

We call (f, g, h) the ATA deviation parameters.  Note that, when f = g = h = 1, we have the case of isotropy.  In that case the constant (Subscript[c, 44]) reduces to μ, the constant (Subscript[c, 13]) (= (Subscript[c, 12])) to λ, and the constant (Subscript[c, 11]) (= (Subscript[c, 33])) to λ + 2μ; the solution is then spherically symmetric.  If λ = μ, the elastic solid becomes a Poisson solid with the compressional velocity (Sqrt[3])times the shear velocity.  

Evaluate the next cell to set up expressions for manipulating the deviation parameters (f, g, h).  It defines the expressions for velocity and slowness for arbitrary transverse anisotropy.  Note that the initial setup of the expressions A, B, C, F, G, H must have been made in the very first cell of this article.

Clear[f,g,h];
Clear[vPATA,vSVATA,vSHATA];
(*Enter 2 coefficients to define isotropic medium; the remaining 3 elastic constants will be defined with f,g,h as shown above.*)
coef11=Input["Enter c11 in MPa: "];
coef12=Input["Enter c12 in MPa: "];
density=Input["Enter ρ in kg/m^3: "];

(*Set up the equations for hexagonal symmetry to represent transverse anisotropy.*)
rules={c14->0,c15->0,c16->0,c22->c11,c23->c13,c24->0,c25->0,c26->0,c34->0,c35->0,c36->0,c45->0,c46->0,c55->c44,c56->0,c66->(c11 - c12)/2};
AA=A/.rules;BB=B/.rules;CC=C/.rules;FF=F/.rules;GG=G/.rules;HH=H/.rules;
equation =Det[{{AA - e ,HH,GG},{HH,BB - e,FF},{GG,FF,CC - e}}];
(*Now solve the cubic equation algebraically for e = ((ρv)^2) and apply the unit vector relation.*)
solution = Solve[equation == 0, e];
eP=Simplify[e/.{solution[[3]]},l^2 + m^2 + n^2 ==1];
eP=eP[[1]];
eSV=Simplify[e/.{solution[[2]]},l^2 + m^2 + n^2 ==1];
eSV=eSV[[1]];
eSH=Simplify[e/.{solution[[1]]},l^2 + m^2 + n^2 ==1];
eSH=eSH[[1]];

vPA=(Sqrt[(eP/ρ)/.{c11->coef11,c12->coef12,c13->g*coef12,c33->h*coef11,c44->f*(coef11 - coef12)/2,ρ->density}])/.{l->Cos[ϑ],m->0,n -> Sin[ϑ]};
vSVA=(Sqrt[(eSV/ρ)/.{c11->coef11,c12->coef12,c13->g*coef12,c33->h*coef11,c44->f*(coef11 - coef12)/2,ρ->density}])/.{l->Cos[ϑ],m->0,n -> Sin[ϑ]};
vSHA=(Sqrt[(eSH/ρ)/.{c11->coef11,c12->coef12,c13->g*coef12,c33->h*coef11,c44->f*(coef11 - coef12)/2,ρ->density}])/.{l->Cos[ϑ],m->0,n -> Sin[ϑ]};

Now we define functions which can be used in a Manipulate in which the ATA parameters (f, g, h) can be varied.

(*tempP =vPA/.{l->Cos[ϑ],m->0,n -> Sin[ϑ]};*)
vPATA[f_,g_,h_]:= Evaluate[vPA];
sPATA[f_,g_,h_]:= Evaluate[1/vPA];

(*tempSV = vSVA/.{l->Cos[ϑ],m->0,n -> Sin[ϑ]};*)
vSVATA[f_,g_,h_]:=Evaluate[vSVA];
sSVATA[f_,g_,h_]:=Evaluate[1/vSVA];

(*tempSH = vSHA/.{l->Cos[ϑ],m->0,n -> Sin[ϑ]};*)
vSHATA[f_,g_,h_]:=Evaluate[vSHA];
sSHATA[f_,g_,h_]:=Evaluate[1/vSHA];

Note that the function Evaluate was wrapped around the three defined functions for P, SV, and SH.  Using these expressions, we can plot the velocity surface’s intersection in a vertical plane that includes the z axis as a function of ϑ, the Cartesian angle, which is counterclockwise from the x-y (n = 0) plane.  Similarly, we can plot the slowness surface’s intersection in the same vertical plane.  Again, we use the early [4, 5] approach to computing the wave surface.  The following cell allows one to display the vertical slice of velocity, slowness, and wave surfaces through the z axis as the ATA parameters (f, g, h) are allowed to vary.  Again, the full surface is invariant in azimuthal angle.  In the following Manipulate cell, the velocity and slowness expressions are local functions which can be manipulated using their parameters (f, g, h) as Manipulate variables.   ContinuousAction is set to False because the expressions to be evaluated are complicated and not amenable to rapid changes.  This same cell shows the intersection of the wave surface in the vertical plane.  

Evaluate this cell to allow one to plot profiles of velocity, slowness, and wavefront while manipulating the ATA deviation parameters (f, g, h).  By manipulating these parameters, one can relate the shape of the wave surface to the amount of deviation of elastic coefficients from an isotropic case.  Cusps, and their associated triplication in arrivals, appear and disappear for the SV wave as these parameters are manipulated.  Try f < 1 and h > 1 to produce cusps in the four quadrants or try f > 1 and h < 1 to produce cusps on the x and z axes.  The ratio of the two isotropic coefficients ((Subscript[c, 11]) and (Subscript[c, 12])), as input, will affect how this wave surface appears for specific (f, g, h).  For values of f, g, or h much different than unity, anomalous wave surfaces can appear, unlikely to be of application in the real Earth.

(*Use the Cartesion angle ϑ (CCW from the x-y plane) rather than the polar angle θ.*)
(*Start with the isotropic assumption f = g = h = 1.*)
inc=1.0*π/180;  (*A small increment is required to define cusps well.*)

(*Include all three waves in the Manipulate.  Start with the isotropic assumption e = f = g = 1.*)
Manipulate[
p1=PolarPlot[{vPATA[f,g,h],vSVATA[f,g,h],vSHATA[f,g,h]},
{ϑ, 0, 2π},PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->12,LabelStyle->12,PlotLabel->"intersection of velocity surface\n with vertical plane containing z axis.",AspectRatio->Automatic,ImageSize->250];
p2=PolarPlot[{sPATA[f,g,h],sSVATA[f,g,h],sSHATA[f,g,h]},
{ϑ, 0, 2π},PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->12,LabelStyle->12,PlotLabel->"intersection of slowness surface\n with vertical plane containing z axis.",AspectRatio->Automatic,ImageSize->250];
(*Compute the wave surface.*)
waveSurfaceATA={};
Do[
temp = ToExpression["v"<>awave[[i]]<>"ATA"];
vgh=temp[f,g,h]*Cos[ϑ] - D[temp[f,g,h],ϑ]*Sin[ϑ];
vgz=temp[f,g,h]*Sin[ϑ] + D[temp[f,g,h],ϑ]*Cos[ϑ];
(*Compute the horizontal and vertical components of the group velocity vg.*)
hvector=Table[Chop[vgh],{ϑ,0.,2π,inc}];
zvector=Table[Chop[vgz],{ϑ,0.,2π,inc}];
AppendTo[waveSurfaceATA,Transpose[{hvector,zvector}]]
,{i,1,3}
];
(*Start with the isotropic assumption f = g = h = 1.*)
p3=ListPlot[waveSurfaceATA,Joined->True,PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->12,LabelStyle->12,AspectRatio->Automatic,PlotLabel->"intersection of wave surface\n with vertical plane containing z axis.",AspectRatio->Automatic,ImageSize->250];
coeflabel= "coefficients:"<>"\nc11 = "<>ToString[coef11]<>"\nc12 = "<>ToString[coef12]<>"\nc13 = "<>ToString[g*coef12]<>"\nc33 = "<>ToString[h*coef11]<>"\nc44 = "<>ToString[f*(coef11 - coef12)/2]<>"\nρ = "<>ToString[density];
p4=Graphics[Text[coeflabel,{0,0},BaseStyle->16]];
GraphicsGrid[{{p1,p2},{p3,p4}},Spacings->{0,50},ImageSize->700],
{{f,1.0},0.5,1.5,Appearance->"Labeled",LabelStyle->16},
{{g,1.0},0.5,1.5,Appearance->"Labeled",LabelStyle->16},
{{h,1.0},0.5,1.5,Appearance->"Labeled",LabelStyle->16},
LabelStyle->16,
ContinuousAction->False,
SaveDefinitions->True,
TrackedSymbols->{f,g,h}
]

We can visualize the rays in a transverse anisotropic medium with variable coefficients using much of the code given above for the wave surface.  The wave surface, or wavefront, is actually the envelope of endpoints of rays at a given time T.  In this case, we again use T = 1.   The following cell allows one to see the ray diagram while manipulating the ATA deviation parameters as above and also while manipulating the number of rays to track.  Each ray set (P, SV, or SH) can be visualized separately.  Note the concentration or focusing of SV rays at the double cusps in each quadrant (try f < 1 and h > 1) or on the axes (try f > 1 and h < 1).

(*Use the Cartesion angle ϑ (CCW from the x-y plane) rather than the polar angle θ.*)
(*Start with the isotropic assumption e = f = g = 1.*)

Manipulate[
inc=2.0*π/nrays;  (*A small increment is required to define cusps well.*)
temp = ToExpression["v"<>wave<>"ATA"];
vgh=temp[f,g,h]*Cos[ϑ] - D[temp[f,g,h],ϑ]*Sin[ϑ];
vgz=temp[f,g,h]*Sin[ϑ] + D[temp[f,g,h],ϑ]*Cos[ϑ];
hvector=Table[Chop[vgh],{ϑ,inc,2π,inc}];
zvector=Table[Chop[vgz],{ϑ,inc,2π,inc}];
endpoints= Transpose[{hvector,zvector}];
origin={0,0};
raysATA=Riffle[endpoints,{origin},{1,-1,2}];
(*Display the rays for specified wave.*)
p1=ListPlot[raysATA,Joined->True,PlotRange->{All,All},AxesLabel->{"x","z"},PlotStyle->Directive[Black,Thin],BaseStyle->12,LabelStyle->12,AspectRatio->Automatic,PlotLabel->wave<>"-wave rays at T = 1 in\n the vertical plane containing z axis.",AspectRatio->Automatic,ImageSize->250];
coeflabel= "coefficients:"<>"\nc11 = "<>ToString[coef11]<>"\nc12 = "<>ToString[coef12]<>"\nc13 = "<>ToString[g*coef12]<>"\nc33 = "<>ToString[h*coef11]<>"\nc44 = "<>ToString[f*(coef11 - coef12)/2]<>"\nρ = "<>ToString[density];
p2=Graphics[Text[coeflabel,{0,0},BaseStyle->16]];
GraphicsRow[{p1,p2},ImageSize->700],
{wave,awave}  (*Choose a wave -- P, SV, or SH.*),
{{nrays,180,"# rays"},2,720,2,Appearance->"Labeled",LabelStyle->16},
{{f,1.0},0.5,1.5,Appearance->"Labeled",LabelStyle->16},
{{g,1.0},0.5,1.5,Appearance->"Labeled",LabelStyle->16},
{{h,1.0},0.5,1.5,Appearance->"Labeled",LabelStyle->16},
LabelStyle->16,
ContinuousAction->False,
SaveDefinitions->True,
TrackedSymbols->{f,g,h,wave,nrays}
]

Summary

Anisotropy means that velocity varies with direction of propagation in an elastic medium.  We started with the most general case of anisotropy which involves 21 independent elastic coefficients in a 6×6 matrix relating stress to strain.  The solution to the wave equation in anisotropic material leads to three waves (P, SV, SH), one compressional wave and two shear waves, each dependent on the 21 coefficients.   Although it is possible to visualize the wave character from a point source in the most general case, further assumptions on symmetry make the problem more tractable and more rapidly visualized.   These assumptions follow from the nature of perfect crystals, with the limit being an isotropic medium with velocity constant in all directions.  Anisotropy for particular crystals exhibiting hexagonal or cubic symmetry were studied in early papers, with illustrations of hand-calculated 2-D profiles of velocity, slowness, and wave surfaces.  We have enabled checking of these early calculations and have enabled visualization of the completely general case of anisotropy with 21 coefficients.  Hexagonal crystals are related to a common anisotropic configuration called transverse anisotropy which is a widely applied assumption for the elastic properties of the Earth, rather than simply isotropic.  Transverse anisotropy with a vertical axis of symmetry entails a different velocity in the vertical direction than in the horizontal direction.  We have formulated a way to visualize wave character in media of arbitrary transverse anisotropy, starting from the limiting case of isotropy.  The wave surface, in particular, responds to variations in transverse anisotropy with a non-spherical shape, often with cusps and triplication of arrivals.

Conclusion

What was only once visualized with painstaking hand calculations can now be done in Mathematica in a few seconds or fraction of a second.  We found that the hand calculations made in seminal papers for the character of elastic wave propagation in anisotropic media were remarkably accurate even though they would have required hours or days to complete.  With Mathematica, anisotropy can be studied in even more detail and is a good example of how mathematical simplifications are often no longer needed.

Appendix

The cell in this appendix provides a way to visualize velocity for the general anisotropic case (21 elastic coefficients) in the vertical plane through the z axis at an arbitrary rotation angle ψ.  The initial display is for an isotropic medium.  The evaluations take several tens of seconds for any change in one or more of the 21 elastic constants, and so this cell is provided for when one is interested in some very specific case of anisotropy.   A similar plot for slowness can easily be made by inverting the velocity expression.

Unprotect[C];  (*C is used intrinsically in Mathematica, but for expressions not relevant to this work.*)
Clear[A,B,C,F,G,H];
Clear[c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66];
Clear[l,m,n];
awave={"P","SV","SH"};

equation= Collect[Expand[Det[{{A - e,H,G},{H,B - e,F},{G,F,C - e}}]],e];
solution = Solve[equation == 0, e];
(*We rearrange the solutions into the order P,SV,SH.*)
rhov2P=e/.solution[[2]];
rhov2SV=e/.solution[[3]];
rhov2SH=e/.solution[[1]];

A=l^2*c11+m^2*c66+n^2*c55+2*m*n*c56 + 2*n*l*c15 + 2*l*m*c16;
B=l^2*c66+m^2*c22+n^2*c44+2*m*n*c24 + 2*n*l*c46 + 2*l*m*c26;
C=l^2*c55+m^2*c44+n^2*c33+2*m*n*c34 + 2*n*l*c35 + 2*l*m*c45;
F=l^2*c56+m^2*c24+n^2*c34+m*n*(c23 + c44) + n*l*(c36 + c45) + l*m*(c25 + c46);
G=l^2*c15+m^2*c46+n^2*c35+m*n*(c36 + c45) + n*l*(c13 + c55) + l*m*(c14 + c56);
H=l^2*c16+m^2*c26+n^2*c45+m*n*(c25 + c46) + n*l*(c14 + c56) + l*m*(c12 + c66);
(*Make starting elastic coefficients be for an isotropic medium.*)
(*Note we are using the angle ϑ in the vertical Cartesian plane, measured CCW from the horizontal axis.*)
Clear[ϑ];
Clear[aψ];
tempP =Sqrt[ rhov2P/.{l->Cos[aψ]Cos[ϑ],m -> Sin[aψ]Cos[ϑ],n->Sin[ϑ]}]/ρ; (*vertical plane at angle aψ*)
vPplot[c11_,c12_,c13_,c14_,c15_,c16_,c22_,c23_,c24_,c25_,c26_,c33_,c34_,c35_,c36_,c44_,c45_,c46_,c55_,c56_,c66_,ρ_,aψ_]:=Evaluate[tempP];

tempSV = Sqrt[rhov2SV/.{l->Cos[aψ]Cos[ϑ],m -> Sin[aψ]Cos[ϑ],n->Sin[ϑ]}]/ρ; (*vertical plane at angle aψ*)
vSVplot[c11_,c12_,c13_,c14_,c15_,c16_,c22_,c23_,c24_,c25_,c26_,c33_,c34_,c35_,c36_,c44_,c45_,c46_,c55_,c56_,c66_,ρ_,aψ_]:=Evaluate[tempSV];

tempSH = Sqrt[rhov2SH/.{l->Cos[aψ]Cos[ϑ],m -> Sin[aψ]Cos[ϑ],n->Sin[ϑ]}]/ρ; (*vertical plane at angle aψ*)
vSHplot[c11_,c12_,c13_,c14_,c15_,c16_,c22_,c23_,c24_,c25_,c26_,c33_,c34_,c35_,c36_,c44_,c45_,c46_,c55_,c56_,c66_,ρ_,aψ_]:=Evaluate[tempSH];

(*Include all three waves in the Manipulate.  This takes several tens of seconds to evaluate the first time and for any change in slider values after that.*)
Manipulate[
PolarPlot[
{
vPplot[c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66,ρ,aψ],
vSVplot[c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66,ρ,aψ],
vSHplot[c11,c12,c13,c14,c15,c16,c22,c23,c24,c25,c26,c33,c34,c35,c36,c44,c45,c46,c55,c56,c66,ρ,aψ]
},
{ϑ, 0, 2π},PlotRange->{All,All},PlotStyle->{{Blue,Thick},{Red,Thick},{Green,Thick}},PlotLegends->awave,BaseStyle->20,PlotLabel->"Intersection of velocity surface \nwith plane at ψ = "<>ToString[aψ,FormatType->TraditionalForm]<>"\ncontaining z axis.",AspectRatio->Automatic],
{{c11,300000.},100000.,500000.,Appearance->"Labeled"},
{{c12,100000.},50000.,150000.,Appearance->"Labeled"},
{{c13,100000.},50000.,150000.,Appearance->"Labeled"},
{{c14,0.0},0.0,100000.,Appearance->"Labeled"},
{{c15,0.0},0.0,100000.,Appearance->"Labeled"},
{{c16,0.0},0.0,100000.,Appearance->"Labeled"},
{{c22,300000.},100000.,500000.,Appearance->"Labeled"},
{{c23,100000.},50000.,150000.,Appearance->"Labeled"},
{{c24,0.0},0.0,100000.,Appearance->"Labeled"},
{{c25,0.0},0.0,100000.,Appearance->"Labeled"},
{{c26,0.0},0.0,100000.,Appearance->"Labeled"},
{{c33,300000.},100000.,500000.,Appearance->"Labeled"},
{{c34,0.0},0.0,100000.,Appearance->"Labeled"},
{{c35,0.0},0.0,100000.,Appearance->"Labeled"},
{{c36,0.0},0.0,100000.,Appearance->"Labeled"},
{{c44,100000.},50000.,150000.,Appearance->"Labeled"},
{{c45,0.0},0.0,100000.,Appearance->"Labeled"},
{{c46,0.0},0.0,100000.,Appearance->"Labeled"},
{{c55,100000.},50000.,150000.,Appearance->"Labeled"},
{{c56,0.0},0.0,100000.,Appearance->"Labeled"},
{{c66,100000.},50000.,150000.,Appearance->"Labeled"},
{{ρ,3000.0},1000.,5000.,Appearance->"Labeled"},
{{aψ,0.0,"ψ"},0,π,π/12,Appearance->"Labeled"},
ContinuousAction->False,
SynchronousUpdating->False,
TrackedSymbols->True
]

Acknowledgments

The author wishes to acknowledge the early investigators of anisotropy, many of them in the reference list below.  With skillful mathematical detail, they unlocked the character of anisotropy in elastic media.  Employing painstaking calculations, they illustrated the complex behavior of wavefronts in anisotropic media before the computer era.

References

M. J. P. Musgrave, ”On the propagation of elastic waves in aeolotropic media: I. General principles”, Proc. R. Soc. Lond. A,  226, 1954 pp. 339-355. doi: 10.1098/rspa.1954.0258

M. J. P. Musgrave, ”On the propagation of elastic waves in aeolotropic media: II. Media of hexagonal symmetry”, Proc. R. Soc. Lond. A,  226, 1954 pp. 356-366. doi: 10.1098/rspa.1954.0259

G. F. Miller and M. J. P. Musgrave, ”On the propagation of elastic waves in aeolotropic media: III. Media of cubic symmetry”, Proc. R. Soc. Lond. A,  236, 1956 pp. 352-383. doi: 10.1098/rspa.1956.0142

G. W. Postma, “Wave propagation in a stratified medium”, Geophysics, 20, 1955 pp. 780-806.  doi: 10.1190/1.1438187

T. Krey and K. Helbig, “A theorem concerning anisotropy of stratified media and its significance for reflection seismics”, Geophys. Prosp., 4, 1956 pp. 294-302.  doi: 10.1111/j.1365-2478.1956.tb01411.x

D. Frederick and T. S. Chang, Continuum Mechanics, Boston: Allyn and Bacon Inc., 1965.

A. E. H. Love, Treatise on the Mathematical Theory of Elasticity, 4th ed.,  New York: Dover Publications, 1944.

F. Gassmann, “Introduction to seismic travel time methods in anisotropic media”, Pure Appl. Geophys., 58,1964 pp. 63-112.  doi: 10.1007/BF00879140

P. F. Daly and F. Hron, “Reflection and transmission coefficients for seismic waves in ellipsoidally anisotropic media”, Geophys., 44, 1979 pp. 27-38.  doi: 10.1190/1.1440920

J. G. Berryman, “Long-wave elastic anisotropy in transversely isotropic media”, Geophys., 44, 1979 pp. 896-917. doi: 10.1190/1.1440984

D. Forsyth, “The early structural evolution and anisotropy of the oceanic upper mantle”, Geophys. J. Int., 43, 1975 pp. 103-162.  doi: 10.1111/j.1365-246X.1975.tb00630.x

V. Vavrycuk, “Parabolic lines and caustics in homogeneous weakly anisotropic solids”, Geophys. J. Int., 152, 2003 pp. 318-334.  doi: 10.1046/j.1365-246X.2003.01845.x

K. Helbig, “Discussion on ‘The reflection, refraction, and diffraction of waves in media with an elliptical velocity dependence’”, Geophysics, 44, 1979 pp. 987-990.  doi: 10.1190/1.1440990