Velocity-Stress Equations for Wave Propagation in Anisotropic Elastic Media

where ρ is the density of the medium, w the displacement, and c[4] the fourth-order stiffness tensor [2]. Equation (1.1) has been derived based on the equation of motion in conjunction with the elastic constitutive equation. Equation (1.1) has been solved by the finite-difference methods, e.g., [21], and the time-domain finite-element methods, e.g., [33], for propagating waves, and the frequency-domain finite-element methods, e.g., [6], for normal mode analysis of standing waves.


Introduction
Conventionally, the second-order elastodynamics equation has been employed to model waves in elastic solids: [4] ∇w , (1.1) where ρ is the density of the medium, w the displacement, and c [4] the fourth-order stiffness tensor [2]. Equation (1.1) has been derived based on the equation of motion in conjunction with the elastic constitutive equation. Equation (1.1) has been solved by the finite-difference methods, e.g., [21], and the time-domain finite-element methods, e.g., [33], for propagating waves, and the frequency-domain finite-element methods, e.g., [6], for normal mode analysis of standing waves.
Alternatively, one can use a modern numerical method [11,12] to solve the first-order velocity-stress equations to obtain the transient solution of wave propagation in elastic solids. For example, Virieux [25] modeled waves in earth crust by using a finite-difference method to solve the velocity-stress equations. LeVeque [12,13] simulated stress waves in isotropic, elastic solids by using a upwind scheme. Shorr [22] developed a specialized formulation of the finite-element method to simulate propagating waves. Käser and Dumbser [10] used the discontinuous Galerkin method [20] to model waves in anisotropic earth crust. Yu et al. [32] simulated non-linear stress-wave propagation in hypo-elastic media by using the space-time Conservation Element and Solution Element (CESE) method [4]. In the setting of the second-order wave equation, Eq. (1.1), ample theoretical analyses about anisotropic elasticity can be found in the literature, e.g., [2,24]. In general, mathematical properties of the first-order velocity-stress equations must be similar to that of the second-order wave equation Eq. (1.1). However, only limited discussions, e.g., Puente et al. [18] and Yang et al. [30], about the eigen structure of the velocity-stress equations can be found in the literature.
The objective of the present chapter is to analyze the mathematical properties of the first-order velocity-stress equations. In particular, their connection to the conventional second-order wave equations, Eq. (1.1) is of interest. To this end, this paper will be focused on analyzing the eigen structure of the velocity-stress equations. The approach employed is identical to that demonstrated in Warming et al. [29]. To be general, the derivation will be based on the properties of a triclinic solid, which have 21 independent stiffness constants and no symmetry. As such, the result will be applicable to all anisotropic, elastic solids. Moreover, in an one-dimensional space, the Jacobian matrix of the velocity-stress equations is shown to be equivalent to the classical Christoffel matrix of the second-order wave equations. In contrast to the conventional approach, however, a harmonic solution of plane waves is not assumed. To demonstrate the capabilities of the new formulation, the CESE method is applied to solve the velocity-stress equations for modeling wave propagation in a block of beryl, a solid with hexagonal symmetry. Wave expansion from a point source is demonstrated. The calculated group velocity profile compared well with the analytical solution. Additional result of waves interacting with interfaces separating regions with different lattice orientations are also reported. All results show salient features of wave propagation.
In the remainder of this section of Introduction, the equations of elastodynamics are reviewed in order to stage the further derivation in the following sections. First, the equation of motion without the body force is considered: where v is the velocity vector and T is the Cauchy stress tensor. T can be related to the strain tensor S by the elastic relation: where c [4] is the fourth-order stiffness tensor. To proceed, Eq. (1.3) is differentiated by time t to yield: where the superscript t denotes transpose. The model equations include Eqs. (1.2) and (1.4). In an three-dimensional space, there are nine independent equations with the velocity and stress components as the unknowns. Next, the Voigt notation [2] is applied and the second-order tensors in Eq. (1.4) is changed to be 6-component vectors: ⎛ (1.5) In the following, boldfaced T = (T 1 , T 2 , T 3 , T 4 , T 5 , T 6 ) t and S = (S 1 , S 2 , S 3 , S 4 , S 5 , S 6 ) t denote the 6-component stress and strain vectors. The non-boldfaced T and S in Eqs. (1.2) (1.3), and (1.4) denote the second-order stress and strain tensors. Aided by the Voigt notation, the elastic constitutive relation is rewritten as T = CS, where C is the second-order 6 × 6 stiffness matrix.
The rest of this chapter is organized as in the following. Section 2 shows the first-order velocity-stress equations in a vector-matrix form. Section 3 proves that the velocity-stress equations are hyperbolic. Section 4 shows that the Jacobian matrix of the one-dimensional velocity-stress equations is equivalent to the classic Christoffel matrix. As a concrete example, Section 5 illustrates a special form of the velocity-stress equations for elastodynamics in beryl, a solid of hexagonal symmetry. Section 6 reports the numerical solution of wave propagation in a block of beryl. The solution is obtained by applying the CESE method to solve the velocity-stress equations. In Section 7 we provide the concluding remarks. At the end of the chapter, a list of cited references and several appendices are attached.

The model equation
To proceed, a Cartesian coordinate system is employed to formulate the velocity-stress equations Eqs. (1.2) and (1.4): where K (μ) with μ = 1, 2, 3 are 3 × 6 matrices: In Eq. (2.2), δ ij is the Kronecker delta. The rank, i.e., the dimension of the column and row spaces [23], of K (μ) with μ = 1, 2, 3 is three. In the above equations, the value of each entry in C (or c [4] ) must be measured based on the same Cartesian coordinates employed in formulating Eq. (2.1). This coordinate system could be arbitrary. However, the entries of C reported in the literature, e.g., Auld [2], are usually tabulated based on a chosen coordinate system, which is aligned with the lattice orientation of the solid of interest.
To proceed, the velocity-stress equations, Eq. (2.1), is recast into a vector-matrix form: where the unknown vector u (2) , and A (3) are the three Jacobian matrices: where 0 m×n denotes a m × n null matrix. Since the rank of matrices K (1) , K (2) , and K (3) is 3, matrices A (1) , A (2) , and A (3) can be easily reduced to their Echelon forms [23] and the result shows that the rank of each of the three 9 × 9 matrices A (1) , A (2) , and A (3) is 6.

Hyperbolicity of the velocity-stress equations
In this section, the hyperbolicity of the velocity-stress equations Eq. (2.3) is proved based on analyzing a linear combination of the three Jacobian matrices: Aided by Eq. In the following, we will show that all eigenvalues of B are real. Then we will show that B is diagonalizable by its eigenvector matrix. To proceed, we introduce the following definition: Recall that |h| = 1. Since h 1 , h 2 , and h 3 are not all-zero, the rank of matrix G is 3. Aided by this definition of G, matrix B can be written as Because the rank of G is 3, the rank of the 6 × 3 matrix CG t must be less than or equal to 3, such that rank(B) ≤ 6. The eigenvalues of B can be obtained by solving its characteristic equation: where I n is a n × n identity matrix and λ the eigenvalue. We assume that eigenvalues of B are either zero or non-zero. If all non-zero eigenvalues of B are real, all eigenvalues of B are real. The existence of zero eigenvalues of B is self-evident because rank(B) ≤ 6 [23]. For non-zero eigenvalues, aided by the Schur complements, Eq. (3.6) can be transformed to be is a 3 × 3 matrix. For completeness, the Schur complements are provided in Appendix 8. To proceed, we recall that the stiffness matrix C is symmetric and positive-definite [17]. Therefore, matrix Γ must be symmetric and positive-definite. As a result, the eigenvalues of Γ, i.e., ρλ 2 , are real and positive. Since the density ρ must be positive, λ 2 is positive. Therefore, the non-zero eigenvalues of B must be real and in positive-negative pairs, i.e., for a positive eigenvalue λ + , there must be a negative counterpart λ − = −λ + .
To recap, the spectrum of B includes the null eigenvalue and three non-zero, positive-negative pairs of eigenvalues. Eigenvalues in pairs could be repeated such that there could be at least one and at most three pairs of distinct non-zero eigenvalues. For isotropic solids, B has two pairs of distinct non-zero eigenvalues. The above proof of the real spectrum of B is underpinned by the property of matrix C, which is symmetry and positive-definiteness. As will be shown in the following, this property of C is equally important in constructing the eigenvector matrix of B.
To proceed, we show that matrix B can be diagonalized by its eigenvector matrix. Essentially, we need to find nine linearly independent eigenvectors for B. First, we show that there exist six linearly independent eigenvectors associated with the non-zero eigenvalues. We then show the existence of three linearly independent eigenvectors associated with the null eigenvalue. To proceed, we consider the eigen problem of matrix B: where λ + is a positive eigenvalue of B and the associated eigenvector containing a 3-component vector q v and a 6-component vector q T . Aided by the definition of matrix B, Eq. (3.5), the eigen problem, Eq. (3.9), can be recast to which gives: By arbitrarily specifying a value to an entry of q v , q v can be obtained by solving Eq. (3.11). q T can then be obtained by substituting q v into Eq. (3.12). As such, the eigenvector q + can be fully constructed.
Essentially, the eigen problem of the 9 × 9 matrix B in Eq. (3.9), is reduced to a eigen problem for the 3 × 3 matrix Γ in Eq. (3.11). Since Γ is a 3 × 3, symmetric, positive-definite matrix, the spectral theorem [23] of linear algebra asserts that Γ always has three linearly independent eigenvectors, regardless whether Γ has degenerate eigenvalues or not. These three eigenvectors of Γ are denoted by q v,1 , q v,2 , and q v, 3 . The three associated eigenvectors of B with positive eigenvalues are denoted by q +,1 , q +,2 , and q +,3 . Aided by Eq. (3.12), q +,1 , q +,2 , and q +,3 are readily determined by q v,1 , q v,2 , and q v, 3 . Because q v,1 , q v,2 , and q v,3 are linearly independent, q +,1 , q +,2 , and q +,3 must also be linearly independent. By following exactly the same procedure, one can derive the three independent eigenvectors q −,1 , q −,2 , and q −,3 associated with the three negative eigenvalues λ − = −λ + . To proceed, since are the six linearly independent eigenvectors of B. The corresponding non-zero eigenvalues of B can be written as 3 . We note that it is possible for B to have degenerate non-zero eigenvalues. Nevertheless, all six eigenvectors associated with the non-zero eigenvalues are linearly independent. This concludes the discussions about the eigenvectors of B associated with the non-zero eigenvalues.
The remaining three eigenvectors of B are associated with the null eigenvalue. Recall that B is a 9 × 9 matrix with rank(B) ≤ 6. Thus, the nullity of B must be greater than or equal to 3, i.e., there are at least three linearly independent vectors spanning the null space of B [23]. Aided by the definition of null space and eigenvector, the basis of the null space of B consists of the eigenvectors associated with the null eigenvalue of B. As such, there must be exactly three eigenvectors associated with the null eigenvalue of B, because (i) there are at most nine linearly independent eigenvectors of B, (ii) there are at least three linearly independent eigenvectors associated with the null eigenvalues of B, and (iii) there are exactly six linearly independent eigenvectors associated with the non-zero eigenvalues of B. These three eigenvectors are denoted by q 0,1 , q 0,2 , and q 0,3 . This concludes the discussion of three linearly independent eigenvectors associated with the null eigenvalue of B.
Together, there are nine linearly independent eigenvectors of B: q ±,i , q 0,i , i = 1, 2, 3. The diagonalizing eigenvector matrix of B can be constructed as: where Q is a 9 × 9 matrix composed of nine column vectors. We then perform similarity , λ +,2 , and λ +,3 with λ 0 = 0. This concludes the discussions of diagonalizability of B by the above nine linearly independent eigenvectors.
As an aside, each of Jacobian matrices A (μ) , μ = 1, 2, 3, is a special case of B with h ν = δ μν , ν = 1, 2, 3. For these special cases, we let q If we let μ = 1, it can be straightforwardly shown that the fifth, sixth, and seventh columns of A (1) are all-zero. For μ = 2, the fourth, sixth, and eighth columns of A (2) are all-zero. For μ = 3, the fourth, fifth, and ninth columns of A (3) are all-zero. For completeness, the eigenvectors associated with the trivial eigenvalue for A (μ) , μ = 1, 2, 3, are tabulated in the following: Note that a 9 × 9 identity matrix I 9 = QQ −1 is inserted in the second term of Eq. (3.16). Since Q −1 is a constant matrix, it can be moved into the partial differentiation operator. Aided by the eigenvalue matrix Λ and the characteristic variableŝ Eq. (3.16) can be rewritten in the characteristic form as where Λ is a diagonal matrix consisting of the eigenvalues of B. Equation (3.17) contains nine decoupled equations:

Coordinate transformation and the Christoffel matrix
In this section, we show how matrix B is connected to the classic Christoffel matrix. Essentially, we will show that matrix Γ in Eq. To connect the present formulation to the classical Christoffel matrices tabulated in the literature, we proceed to transform matrix Γ formulated in global coordinates to beΓ in local coordinates. We first define the following matrices associated with the Jacobian matrices A (μ) , μ = 1, 2, 3: which are special cases of Γ with h ν = δ μν , ν = 1, 2, 3. Γ (1) , Γ (2) , and Γ (3) are tensors in E 3 and they can be transformed to local coordinates by using the rotation matrix R shown in Eq. (8.1) in Appendix 8:Γ v,i are the same vector defined in the global and local coordinate systems, respectively. The algebraic relation between the Γ (μ) andΓ (μ) in Eq. (4.2) needs to be further elaborated. To proceed, Bond's matrix [2] is used to transform the stiffness matrix from the local coordinates to the global coordinates: Details of the transformation are provided in Appendix 8. For a wide range of media, entries ofC in the local coordinate system are widely available in the literature, e.g., [2]. Substituting Eqs. (4.1) and (4.3) into Eq. (4.2) gives: where is a 3 × 6 matrix with its entries as the components of the direction cosine corresponding to x μ -axis. To proceed, Appendix 8 shows that: The algebra involved is lengthy but straightforward. Aided by Eq. (4.6) and the orthogonality of matrix R, Eq. (4.5) becomes To proceed, the relation between Γ andΓ is derived in the following. Let L be the linear combination of L (1) , L (2) , and L (3) : Aided by Eq. (4.5), Eq. (4.8) can be expanded and becomes Aided by Eqs. (3.8), (3.4), and (4.9), Γ can be readily transformed toΓ: Next, the unit vector h is also transformed into the local coordinates: By comparing Eqs. (4.10) and (4.12) to the Christoffel matrix tabulated in the literature, e.g., Auld [2], it can be seen thatΓ in Eq. (4.10) is indeed the Christoffel matrix. Special cases of Γ in the directions along the axes of the local coordinate system, i.e.,x 1 -,x 2 -, andx 3 -axis, arē Γ (μ) , μ = 1, 2, 3, respectively.

Elastic solids with hexagonal symmetry
In this section, we illustrate how the above formulation is applied to model waves in elastic solids of hexagonal symmetry. A solid of hexagonal symmetry is shown in Figure 1. The Cartesian coordinates are chosen such that the x 3 -axis is aligned with the hexagonal cylinder of the medium, and the x 2 -axis is aligned with the crystal lattice. In this coordinate system,

3) become
Similarly, for A (2) and A (3) , the related matrices are The above Jacobian matrices in the model equations are valid only if the local coordinates (x 1 , x 2 ,x 3 ) are employed, in which the entries of the stiffness matrix such as that shown in Eq. (5.1) can be directly taken from the literature. To proceed, we transform the model equations from the local coordinates (x 1 ,x 2 ,x 3 ) to the global coordinates (x 1 , x 2 , x 3 ). In general, the rotation can be defined by the three Euler angles [1]: α, β, and γ. The coordinate transformation matrix R, Eq. (8.1), can be readily specified: cos γ cos β cos α − sin γ sin α cos γ cos β sin α + sin γ cos α − cos γ sin β − sin γ cos β cos α − cos γ sin α − sin γ cos β sin α + cos γ cos α sin γ sin β sin β cos α sin β sin α cos β ⎞ ⎠ .

Modeling wave propagation in a block of beryl
In this section, we apply the space-time Conservation Element and Solution Element (CESE) method [4] to solve the velocity-stress equations for modeling wave propagation in a block of beryl, a elastic solid of hexagonal symmetry. The CESE method is a modern numerical method for time accurate solutions of linear and non-linear hyperbolic partial differential equations. The CESE method treats space and time in a unified way in calculating the space-time flux. The method is simple, robust, and its operational count is comparable to that of a second-order central difference scheme. The backbone of the CESE method is the a scheme [4], which is neutrally stable without artificial dissipation for solving a scalar convection equation with a constant coefficient. For more complex wave equations, the a scheme is extended with added artificial damping, including the a-α-scheme [4] for shock capturing and the c-τ scheme [5] for Courant-Friedrichs-Lewy (CFL) number insensitive calculations. To model multi-dimensional problems, the CESE method uses unstructured meshes [27] with triangles and tetrahedra as the basic elements in two-and three-dimensional spaces. Extensions to use quadrilaterals and hexahedra are also available [36]. In the past, the CESE method has been successfully applied to solve a wide range of various fluid dynamics problems, including aero-acoustics [14], cavitations [19], complex shock waves [8], detonations [26], magnetohydrodynamics (MHD) [34,35], etc. Moreover, the method has been applied to model electromagnetic waves [28] and nonlinear stress waves in isotropic solids [3,31,32]. In the present chapter, the following results are obtained by using the a-α-scheme.
The two-dimensional velocity-stress equations are employed to model wave propagation in a block of beryl, which is a solid of hexagonal symmetry. The material properties are taken from [16]: ρ = 2.7 g/cm 3 , c 11 = 26.94 × 10 11 dynes/cm 2 , c 12 = 9.61 × 10 11 dynes/cm 2 , c 13 = 6.61 × 10 11 dynes/cm 2 , c 33 = 23.63 × 10 11 dynes/cm 2 , and c 44 = 6.53 × 10 11 dynes/cm 2 . The computational domain is −1 m ≤ x (1) ≤ 1 m and −1 m ≤ x (2) ≤ 1 m. A unstructured mesh composed of triangular elements is used. The square domain is divided into 2.2 millions of triangular elements. The maximum length of the element edge is 2.68 × 10 −3 m. The initial condition is a two-dimensional step function v = aH(σ − |x|), where a represents a constant initial velocity to initiate the expanding wave from the origin. The constant σ in the initial condition is chosen to be 5 × 10 −3 m.
The two-dimensional code is parallelized by domain decomposition. First, a graph of element connectivity is built based on the unstructured mesh. The connectivity graph is processed by METIS [7] to partition the domain. As shown in Fig. 3, the overall spatial domain is decomposed into 16 sub-domains for parallel computing. Numerical calculations for elements in each sub-domain is distributed to a workstation as a part of a networked cluster. In all cases, Δt = 65 × 10 −9 s and the maximum CFL number ν is about 0.95. Three cases of wave propagation are calculated: (1) calculation of the group velocity profile to assess  In Case 1, two-dimensional numerical simulations are performed to calculate the density of the total energy, which is the summation of the kinetic energy and the strain energy normalized by area. For elastic solids, the group velocity is identical to the energy velocity [2]. Here, we consider wave propagation in 100 (x 2 -x 3 ) or 010 (x 1 -x 3 ) plane of the hexagonal solid. The orientation of the solid is set to be θ = 0 • and φ = 90 • , to have x (1) aligned to x 3 , and x (2) aligned to x 2 . The group velocities have the analytical solution. The group velocities of pure-shear (SH), quasi-longitudinal (qL), and quasi-shear (qS) waves in a block of beryl are plotted in Fig. 4(a) [16]. The simulation result, shown in Fig. 4(b), compares well with the analytical solution.
In Case 2, the orientation of the computational plane is in the direction of θ = 20 • and φ = 70 • . Similar to that in Case 1, the initial condition is a velocity source at the origin with σ = 0.00 m. Figure 5 shows two snapshots of the calculated energy profiles at t = 26 and 130 μs. With the non-reflecting boundary condition implemented on the four boundaries of the computational domain, waves propagate out of the domain without spurious reflection.
For Case 3, we calculate wave propagation in a domain composed of three blocks of beryl with three different lattice orientations. The partition of the domain is shown in Fig. 6. Shown in Fig. 7(a), an initial, cylindrical wave is generated at the origin. This is because of the lattice orientation of the central block such that an isotropic wave expansion occurs. This wave expands and interacts with the interfaces separating the neighbouring blocks of beryl with different lattice orientations. Figure 7 shows four snapshots of the calculated energy profiles at different times. Due to complex wave/interface interactions, all wave modes, including SH, qL and qS are excited at the interface. As a results, refracted waves at different speeds can be seen in Figs. 7(b) and 7(c). The original circular wave front is fractured after passing through the interfaces, as shown in Fig. 7(d). Complex wave features have been successfully captured by the numerical results.

Concluding remarks
In this chapter, we have presented the first-order velocity-stress equations as an alternative to the conventional second-order wave equations for modeling wave propagation in anisotropic elastic solids. The eigen structure of the equations has been thoroughly studied by analyzing the composed Jacobian matrix of the first-order equations, i.e., B = ∑ 3 μ=1 h μ A (μ) . The hyperbolicity of the equations has been rigorously proved by showing the real spectrum of matrix B and diagonalizability of B. B has a degenerate null eigenvalue. For non-zero eigenvalues, B has at least one and at most three positive-negative pairs of real eigenvalues. Moreover, B has nine linearly independent eigenvectors, which can be used to diagonalize B. The eigenvalues represent the wave speeds, and the eigenvectors are related to the polarization of the waves. Since the first-order velocity-stress equations are hyperbolic, the Cauchy problem is well-posed [9]. For numerical simulation, the proven hyperbolicity justifies the use of modern finite-volume methods, originally developed for solving conservation laws, to solve the velocity-stress equations. For upwind methods [12] and discontinuous Galerkin methods [18], information about eigenvalues and eigenvectors of matrix B are critically important. The derived characteristic relations of the governing equations can be directly used to derive the Riemann solver and the flux function as the building blocks of the methods. Moreover, by clearly defining the local and global coordinate systems, we also show that matrix B is directly connected to the classic Christoffel matrix. Therefore, nine coupled first-order velocity-stress equations are shown to be directly related to the classic second-order elastodynamic equations. To demonstrate the application of the velocity-stress equations, we apply the space-time CESE method to solve the equations for modeling wave propagation in a block of beryl. The calculated energy profiles compare well with the analytical solution. We also simulate wave propagation in a heterogeneous solid composed of three blocks of beryl with different lattice orientations. Numerical results show complex wave/interface interactions.

A: The Schur complements
Consider matrix M, which can be divided into 4 sub matrices: where A, B, C, and D are p × p, p × q, q × p, and q × q matrices, respectively. If A −1 and D −1 exist, then matrices D − CA −1 B and A − BD −1 C are the Schur complements of A and D [15]. In this case, matrix M can be factored out as This factorization is useful to calculate det(M).

B: Rotation of Cartesian coordinate system
This appendix provides the transformation relations for the global and local coordinate systems. Let e 1 , e 2 , and e 3 be the natural basis of the global coordinate system, and let e 1 ,ē 2 , andē 3 be the natural basis of the local coordinate system. An arbitrary vector b can be represented by using either the global coordinates or the local coordinates as: is the rotation matrix and it is orthonormal, i.e., R −1 = R t or ∑ 3 k=1 r ik r jk = ∑ 3 k=1 r ki r kj = δ ij , i, j = 1, 2, 3. The natural bases of the global and local coordinate systems are also related through R defined in Eq. (8.1): r ji e j , i = 1, 2, 3. (8.2) To proceed, consider the transformation of a tensor E: where E is defined in the global coordinates, andĒ in the local coordinates.

C: Rotation of governing equations
This appendix shows that the velocity-stress equations, Eq. (1.4), remain in the same form when formulated in different Cartesian coordinate systems. To proceed, we consider an arbitrary tensor E. In the following, we show that ∇ · E =∇ · (ER), where∇· is the divergence operator defined in the local coordinate system, and R is the rotation matrix R defined in Appendix 8. Consider the ith component of ∇ · E with i = 1, 2, 3: The above equation shows that ∇ · E =∇ · (ER) because R is a constant spatial tensor.
Next, apply coordinate transformation to the equation of motion, Eq. (1.2): Aided by Eq. (8.5), Eq. (8.6) can be rewritten as: Then perform the coordinate transformation for constitutive equations, Eq. (1.3): R t TR = R t c [4] SR. Recast c [4] in Eq. (8.7) into the index form: T mn = r im T ij r jn = r im c ijkl S kl r jn = r im c ijkl r koSop r l p r jn = r im r jn r ko r l p c ijklSop =c mno pSop , (8.8) where c ijkl andc mno p are the components of the fourth-order stiffness tensor c [4] andc [4] , respectively. For conciseness, Eq. (8.8) uses Einstein summation convention. By using the direct notation, Eq. (8.8) can be recast asT =c [4]S .