The Directional Dependence of Elastic Stiffness and Compliance Shear Coefficients and Shear Moduli in Cubic Materials

The mathematics associated with the representation surfaces for elastic stiffness and compliance shear coefficients and shear moduli for planes in the standard 001–011–111 stereographic triangle for cubic materials is examined as a function of the degree of anisotropy, A=2c44/(c11−c12), of the material under consideration. Similarities and differences between the surfaces for the elastic stiffness and compliance shear coefficients and the shear moduli are highlighted.

direction under consideration and the length of which is proportional to the shear modulus. These graphical plots either had the direction of observation in the shear plane (Figs. 6 and 7 of their paper) or parallel to the normal to the shear plane ( Fig. 8 of their paper) and clearly show the directional dependence of G r for planes other than {001} or {111}.
Later, in a series of papers, Turley and Sines [6][7][8] reported the results of numerical computations of the directional dependence of elastic stiffness and compliance shear coefficients and shear moduli in silicon, copper and molybdenum using the same graphical representations as Workman and Evans, with, in addition, a very powerful graphical representation using a standard stereographical triangle on which to represent c rr and G r for silicon for various (hkl). The mathematical expressions given by these authors for the shapes of c 44 and G 4 as a function of the shear direction and the normal to the shear plane are highly complex, although readily amenable to numerical computation.
Turley and Sines made a number of interesting observations about their numerical results without showing formal mathematical proofs. For example, it is apparent from Fig. 2 of [6] that for planes of the form (0kl), (hhl) and (hkk) defining the borders of the standard 001-011-111 stereographic triangle for cubic materials, the polar plots of c rr as a function of shear vector in the plane are symmetrical about [100] and [0lk], [110] and [llh], and [011] and [khh] respectively, and have 2mm point symmetry at the origin. However, for a general plane (hkl) within this standard stereographic triangle, such as (156), the two perpendicular directions about which the polar plots are symmetrical are not specified in terms of directions [uvw] lying in (hkl). For the specific plane (156) Turley and Sines show graphically in [6] that these directions will be rotated by what they term a phase angle φ from the two perpendicular directions [510] and [3 15 13] lying in (156).
In addition, they comment in [8] that the shapes (or loci) of the shear stiffness coefficient, c 44 , for silicon when plotted on the standard stereographic triangle closely resemble those of the shear modulus, G 4 , but are not identical, without quantifying this degree of resemblance.
Subsequent to the work of Turley and Sines, Walpole [9] independently examined the variation of shear modulus as a function of the shear direction and the normal to the shear plane and identified mathematically the extremum values of the shear modulus as (i) when either the shear direction is parallel to 001 or the shear plane is {001}, and (ii) when the shear direction is 110 and the shear plane is {110}. These extrema are apparent in the numerical simulations of Turley and Sines [7], e.g., their Fig. 3(a).
In the work reported here we have revisited the mathematics associated with the representation surfaces for c rr , s rr and G r for (hkl) in the standard 001-011-111 stereographic triangle for cubic materials, paying particular attention to the specification of principal axes of these polar plots for different choices of (hkl) and the forms of the loci of these plots with respect to these principal axes. We have also extended our analysis to the examination of polar plots for c qr and s qr for q, r = 4, 5 and 6 with q = r, the coefficients of which are zero when the rotated axis system within which c qr and s qr are specified aligns with the cubic axes defining the lattice of the material under consideration. As others have noted [6,8], these considerations are relevant to the study of imperfections in anisotropic crystals, such as the strain fields and stress fields associated with point defects and dislocations (e.g., [10][11][12][13][14]).

General Stiffness and Compliance Tensor Transformation Relations
For an arbitrary rotation of axes from one axis system to another, tensors of the fourth rank T ij kl (in full tensor notation) transform as for i, j , k, l, m, n, p and q all taking values from 1 to 3 [3,15], where the a im are direction cosines specifying the angle between the ith axis of the 'new' axis system and the mth axis of the 'old' axis system. The stiffness tensor C ij kl and the compliance tensor S ij kl relate a stress σ ij to a strain ε kl through the equations σ ij = C ij kl ε kl and ε ij = S ij kl σ kl .
For our purposes here the 'old' axis system is the orthonormal axis system aligned with respect to the 100 directions of the cubic crystal. Sixty coefficients of tensors of the fourth rank are zero in this axis system for cubic crystals. The twenty one non-zero C ij kl and S ij kl for cubic crystals are (i = j ) [3,15], in which the contracted two suffix Voigt notation is denoted by c and s for stiffness and compliance to distinguish these from the full tensor C and S components. The non-zero c ij and s ij are related to one another through the equations with the shear moduli G 4 = G 5 = G 6 = 1/s 44 for cubic crystals in this axis system. Thomas [16] has shown that the C ij kl and S ij kl transform from axes 1, 2 and 3 aligned with respect to the 100 directions of the cubic crystal to a new orthonormal set of axes 1 , 2 and 3 , so that in the new axis system the C ij kl and S ij kl can be written in the succinct forms C ij kl = c 12 δ ij δ kl + c 44 (δ ik δ jl + δ il δ jk ) + (c 11 − c 12 − 2c 44 )a iu a ju a ku a lu , 4S ij kl = 4s 12 δ ij δ kl + s 44 (δ ik δ jl + δ il δ jk ) + (4s 11 − 4s 12 − 2s 44 )a iu a ju a ku a lu , (5) in which δ is the Kronecker delta and the dummy suffix u takes the values 1, 2 and 3.

Directional Dependence of the Elastic Stiffness and Compliance Shear Coefficients
To establish how the shear modulus varies as a function of shear plane and shear direction, the co-ordinate transformation from the 'old' set of axes to the 'new' set of axes can be chosen without loss of generality so that axis 3 is parallel to the normal of the plane of shear and 2 is parallel to the direction of shear. Hence, It is apparent from Eqs. (6) and (7) that interchanging directions 3 and 2 , i.e., interchanging the direction of shear and the normal to the plane of shear, has no effect on either c 44 or s 44 , so that these are invariant with an exchange of plane of shear and direction of shear. It is also evident that a further simplification can be made to these two coefficients.
It also follows that both can be expressed in terms of c 44 and the anisotropy ratio [3,15,17]: by the Weiss zone law [15].
If A = 1, as will be the case for isotropic cubic materials such as W ( [12], p. 837), Eqs. (6) and (7) reduce to the equations c 44 = c 44 and s 44 = s 44 = 1/c 44 respectively. It is also evident from Eqs. (11) and (12) that c 44 , s 44 and G 4 = 1/s 44 can all be normalised with respect to c 44 so that they are only functions of A.

Plane of Shear {001} or Direction of Shear 001
If 3 is the plane (001) in the standard stereographic triangle, a 31 = a 32 = 0 and a 33 = 1. Since 2 lies in (001), a 23 = 0, so that c 44 = c 44 and s 44 = s 44 = 1/c 44 , and so G 4 = c 44 , i.e., if the plane of shear is (001) the elastic stiffness coefficient c 44 and the shear modulus are both identical and independent of the direction of shear in (001) and A. By symmetry this result is true for {001} or on any plane when the direction of shear is 001 , as we have already noted. c 44 and G 4 therefore plot as identical circles on the representation of the shear modulus within a plane used by Wortman and Evans [1] and Turley and Sines [6].

Plane of Shear {111} or Direction of Shear 111
If 3 is the plane (111) in the standard stereographic triangle, a 31 = a 32 = a 33 = 1/ √ 3. Hence, Once again, it is evident that the elastic stiffness coefficient c 44 and the shear modulus are both independent of the direction of shear in (111), and that they will be identical if A = 1. By symmetry this result is true for {111} and when the direction of shear is 111 . However, the circles representing the loci of c 44 and G 4 will have different radii unless A = 1 because For a material such as Si for which A = 1.564, this ratio is 1.045, and for Cu with A = 3.21, this is 1.338. For highly anisotropic cubic materials such as δ-Pu at room temperature with A = 7.03 and Li at 78 K with A = 9.39 [18], c 44 /G 4 rises to 2.149 and 2.666 respectively.

Plane of Shear {hkl}
Without loss of generality, we can consider the specific plane (hkl) with h ≤ k ≤ l within the standard stereographic triangle. We will consider plotting the loci of c 44 and s 44 first. Suppose two orthonormal directions β 1 and β 2 are identified within (hkl) that together with the direction cosines corresponding to the normal to (hkl), 3 , form a right-handed set of axes. The table of direction cosines between the axis set 1, 2 and 3 and the axis set β 1 , β 2 and 3 will be of the form Hence, with respect to β 1 and β 2 , a shear direction of unit length making an angle θ with β 1 and (90 • − θ ) with β 2 will have direction cosines  6) and (7) readily shows that there will be no sin θ cos θ terms in these equations if a 2 31 α 11 α 21 + a 2 32 α 12 α 22 + a 2 33 α 13 α 23 = 0.
Under these circumstances, Eqs. (6) and (7) are both of the general form where p = P + (Q+R) Consideration of this polar equation shows that, for non-zero values of q, it represents the locus of the solutions of a sextic equation; defining tan θ = y/x, Eq. (20) becomes This is one form of the sextic equation The locus defined by Eq. (20) is symmetric about the two perpendicular directions defined by θ = 0 • and θ = 90 • . Hence, for two randomly chosen perpendicular directions, β 1 and β 2 , Eq. (19) enables the phase angle φ defined by Turley and Sines [6] to be determined. Equation (19) can be further simplified. Since β 1 , β 2 and 3 form an orthonormal righthanded set of vectors, so that the condition expressed by Eq. (19) becomes This together with the condition that β 2 lies in (hkl), enables the two perpendicular directions about which the polar plots are symmetrical for a general plane (hkl) within this standard stereographic triangle to be found. If the required directions are taken to be of the form [1vw], then v is a solution of the quadratic equation  [6], with axes 1 , 2 and 3 in our notation permuting cyclically to become 2 , 3 and 1 , in their notation. That these are the perpendicular directions defining the principal axes of the loci of c 44 , s 44 and G 4 for (hkl) within the standard stereographic triangle can be confirmed by a second consideration of Eqs. (6) and (7) with the substitution in Eq. (18). Defining the function under these circumstances, an examination of f (θ) shows that it has turning points whenever For planes on the border of the standard 001-011-111 stereographic triangle, identification of the two perpendicular directions about which the polar plots are symmetrical is more straightforward. Planes on the great circle connecting 001 and 011 are of the form (0kl) and so in terms of their directions cosines, a 31 = 0, giving solutions from Eqs. (24) and (25)  and ±[l l 2h]. Finally, for planes of the form (hkk) a 32 = a 33 , and so these directions are ±[011] and ±[2khh]. Having determined the two perpendicular directions about which the polar plots are symmetrical, the forms of the polar plots can be examined for planes of the various special and general forms.

Plane of Shear (0kl)
For (0kl) the table of direction cosines between the axis set 1, 2 and 3 and the axis set β 1 , β 2 and 3 can be chosen to be of the form with a 2 32 + a 2 33 = 1. Substituting from Eq. (18), Eqs. (11) and (12) become and respectively. The polar plots of c 44 and s 44 depart from circles most noticeably for planes on the great circle connecting 001 and 011 at (011), taking the forms For A > 1, c 44 will take its minimum value and s 44 its maximum value for (0kl) when θ = 0 • and a 33 = 1/ √ 2, i.e., for (011) when the shear direction is [011], when these equations simplify to and so for this particular shear system G 4 = c 44 . Clearly, for A < 1, c 44 will take its maximum value and s 44 its minimum value for shear on (011)

Plane of Shear (hhl)
For (hhl) the table of direction cosines between the axis set 1, 2 and 3 and the axis set β 1 , β 2 and 3 can be chosen to be of the form with 2a 2 31 + a 2 33 = 1. Substituting from Eq. (18), Eqs. (11) and (12) become and respectively. The polar plots of c 44 and s 44 are circles for (001) and (111) is a maximum. This leads to the following condition on a 31 : For Si with A = 1.564, this function takes a maximum value of 1.159 when a 31 = 0.422 between (001) and (111), so that (hhl) is close to (112); for Li, this function has a maximum of 1.802 when a 31 = 0.451, so that (hhl) is close to (335).

Plane of Shear (hll)
For (hll) the table of direction cosines between the axis set 1, 2 and 3 and the axis set β 1 , β 2 and 3 can be chosen to be of the form respectively. These polar plots vary from a circle at (111) to the polar plot at (011) which accentuates anisotropy the most.

Plane of Shear (hkl)
For (hkl) planes within the standard stereographic triangle the full forms of Eqs. (11) and (12) have to be used with Eq. (18). If the methodology used by Turley and Sines [6] is used for the polar plots of polar plots of c 44 s 44 and G 4 , so that axes are referred to the meridian tangent in the plane and the direction in the plane perpendicular to the meridian direction, as for example in their Fig. 3, the principal axes of the polar plots will be rotated with respect to the axes chosen by a predictable angle θ defined by a suitable solution of Eq. (28).

Polar Plot Representations for c 4s 4and G 4
As has been noted in Sect. 3, both C 3232 and S 3232 are expressions which can be rearranged to be of the general form r = p + q cos 2θ (Eq. (20)) with respect to the two perpendicular axes about which the polar plots are symmetrical for the (hkl) plane under consideration. The most extreme departure from a circle occurs for shear on the (011) plane in the 001-011-111 stereographic triangle, but it is evident from Eq. (32) that even for this plane p > |q| for both C 3232 and S 3232 , and so there is a defined limit with respect to which the cos 2θ term dominates the polar plot of a circle represented by r = p.
Since C 3232 and S 3232 are quantities defined by both a plane and a direction, it is not possible to plot them on a single surface, as is routinely done for Young's modulus and Poisson's ratio in cubic materials (see, for example, [1]). It is therefore necessary to explore alternative methods of representing these quantities graphically. One method of presenting a tensor component T 3232 and its reciprocal 1/T 3232 introduced by both Wortman and Evans [1] and developed by Turley and Sines [6] is to view the polar plot in the 1 -2 plane, i.e., viewed down the antiparallel direction to the normal to the plane of shear, 3 . The development of Turley and Sines was to group together a number of such plots on a standard 001-011-111 stereographic triangle. Examples of such representations are shown in Figs. 1, 2, 3 for c 44 s 44 and G 4 for (a) Cu, (b) δ-Pu and (c) Nb, for which A = 0.51 [19], using the representation of Turley and Sines on the 001-011-111 stereographic triangle. It can also be instructive to examine the polar plots for (0kl), (hhl) and (hll) planes as groups, as in the example shown for c 44 for Cu in Fig. 4.
For a general (hkl) in Figs. 1, 2, 3, the polar plots are oriented with respect to the vector [uv0] common to (hkl) and (001). Hence, for (0kl) planes, the plots are oriented with respect to [100]; they are therefore symmetric about the vertical and horizontal plotting axes. For (hhl) planes the plots are oriented with respect to the common [110] direction, and so the principal axes are at 45°to the horizontal plotting axes.
For the (hll) planes, the common [011] direction, and therefore the common principal axis, follows the great circle between 111 and 011. Finally, for a plane such as (123), the vector [210] common to (hkl) and (001) defines the direction and sense with respect to which the principal axes of the polar plot must be oriented.
It is evident from a comparison of Figs. 1 and 3 for c 44 and G 4 that there are clear differences in the shapes of the polar plots for c 44 and G 4 if the plane of interest is not (001).  Fig. 1(a). The locus at 001 has unit radius in units of c 44 . (c) Stereographic projection of polar plots of c 44 for Nb (A = 0.51) viewed as in Fig. 1(a). The locus at 001 has unit radius in units of c 44 The differences are most apparent at and around the 011 pole, i.e., one of the likely slip planes for dislocations in b.c.c. metals, as is shown in Fig. 5 for the 011 pole for δ-Pu. These differences between c 44 and G 4 for highly anisotropic cubic materials need to be borne in mind when considering their mechanical behaviour. For situations involving the imposition of external shear stresses to a single crystal of a cubic material, s 44 is the relevant shear coefficient to consider, and therefore G 4 is the effective shear modulus. For local strain fields around dislocations in single crystals, c 44 is the relevant shear coefficient usually considered [12], although others use compliance coefficients for their analysis when deriving expressions such as the energy of a dislocation or its displacement field [20,21].

Further Polar Plot Representations for the Elastic Stiffness and Compliance Shear Coefficients and Shear Moduli
A second method used by Wortman and Evans [1] and Turley and Sines [6] of plotting coefficients is to view down the direction antiparallel to 3 so that both the normal to the plane of shear, 1 , and the shear direction lying in the plane 2 , rotate in the plane of observation.
To consider such representations, we shall use the convention adopted by both these sets of authors; as we have noted in Sect. 3, the elastic stiffness and compliance shear coefficients are invariant with an exchange of the direction of shear and the normal to the plane of shear.
In the general case, the methodology in Sect. 3.3 can be followed, so that the loci for C 1212 ≡ c 66 , 4S 1212 ≡ s 66 and G 6 (with the subscript for G in the contracted notation to denote shear in the 2 direction on the plane whose normal is in the 1 direction) are plotted with respect to principal axes in the plane aligned in the plane. The methodology for C 1212 will be considered here primarily; similar methodologies will pertain for S 1212 and G 6 as noted below.
Using (17) for the relationship between the 1, 2 and 3 and the axis set β 1 , β 2 and 3 , we can then define the perpendicular orthonormal directions 1 and 2 so that  Fig. 1(a). The locus at 001 has unit radius in units of (1/c 44 ). (b) Stereographic projection of polar plots of s 44 for δ-Pu (A = 7.03) viewed as in Fig. 1(a). The locus at 001 has unit radius in units of (1/c 44 ).
A consideration of Eq. (45) shows that there will be no term in sin 2θ cos 2θ , i.e., no term in sin 4θ if This is one form of the polynomial equation of order 10: The locus defined by Eq. (48) is symmetric about the two perpendicular directions defined by θ = 0 • and θ = 90 • . Making the substitutions which follow from the properties of a rotation matrix defined by the table of direction cosines in Eq. (18) [3,15], and recognising that it is apparent that Eq. (47) can be recast into the form 2α 2 21 + a 2 31 α 11 α 21 + 2α 2 22 + a 2 32 α 12 α 22 + 2α 2 23 + a 2 33 α 13 α 23 = 0.
(53) Fig. 3 (a) Stereographic projection of polar plots of G 4 for Cu (A = 3.21) viewed as in Fig. 1(a). The locus at 001 has unit radius in units of c 44 . (b) Stereographic projection of polar plots of G 4 for δ-Pu (A = 7.03) viewed as in Fig. 1(a). The locus at 001 has unit radius in units of c 44 . (c) Stereographic projection of polar plots of G 4 for Nb (A = 0.51) viewed as in Fig. 1(a). The locus at 001 has unit radius in units of c 44 This condition is different from that expressed by Eq. (19). However, for planes on the border of the standard 001-011-111 stereographic triangle, the two perpendicular directions about which these polar plots are symmetrical when viewed down the antiparallel direction to 3 are the same as for C 3232 and S 3232 -the difference only manifests itself for planes within this standard 001-011-111 stereographic triangle. For these planes, defining the function it is evident that g(θ) has turning points whenever confirming that the loci for C 1212 , S 1212 and G 6 have 4mm point symmetry at the origin and that, with respect to the principal axes, the loci for C 1212 , S 1212 are both of the form Thus, for example, for (156), if β is chosen to be parallel to [510], the solution of Eq. (55) confirms that an anti-clockwise rotation of 9.12°is required to align the axes with the principal axes of the loci of C 1212 , S 1212 and G 6 , aligned along maxima of C 1212 and G 6 and minima for S 1212 for A > 1, and vice-versa for A < 1.
Solutions to Eq. (55) on these criteria are shown in Table 2. A comparison of the data in Tables 1 and 2 shows that the directions η 1 and η 2 in the two tables are very similar for (156), the plane considered by Turley and Sines for Figs. 3 and 4 of [6], but that they are noticeably different for other planes such as (125), (123) and (234). A demonstration of this is shown in Fig. 6 for C 3212 and C 1212 for Cu plotted on the (123) plane.
that is, the condition on the two perpendicular directions about which the polar plots are antisymmetrical for a general plane (hkl) is the same for tensors of the form T 3232 and tensors of the form T 3231 . Hence, as is evident in [6], plots of C 3132 and C 3232 (C 1312 and C 1313 in

Average Shear Elastic Constants
As Hirth and Lothe [12] note, the complexity of most problems involving dislocations of differing orientation and Burgers vectors is such that average elastic constants must be used. For such averages, Voigt averages over C ij kl are appropriate for situations in which uniform strain can be assumed, while Reuss averages over S ij kl are appropriate for situations in which uniform stress is a better assumption. Hirth and Lothe state that both averages agree to first order in the anisotropy, but differ to second order (pp. 425, 430 and 431). 1 We note that the representation chosen by Turley and Sines [6] to show C 1312 graphically has 2mm point symmetry and is of the form so that an (arbitrary) offset r 0 > b is added to C 1312 to enable its graphical representation for an angle θ to be always positive definite. An equivalent representation to show the sin 2θ dependence of C 1312 would be simply to plot C 1312 vertically against the polar angle θ horizontally.
For the specific case of the average elastic stiffness and compliance shear coefficients in cubic materials,c 44 ands 44 respectively, the results quoted by Hirth and Lothe can be shown to be:c so thatc Hence, for Cu, for which A = 3.21, this ratio is 1.36, while for δ-Pu (A = 7.03), this ratio is 2.24. While Eq. (62) has the interpretation given by Hirth and Lothe, it is apparent that significant errors can arise in calculations for highly anisotropic cubic materials if the wrong average is unwittingly assumed.

Alternative Graphical Representations of the 'Torsion Modulus' or 'Rigidity Modulus'
An excellent account of the theory of elastic anisotropy and practical considerations of how to measure elastic constants up to the end of 1944 is given by Hearmon [22]. Of particular interest here is the measurement of what Hearmon terms the rigidity modulus, and what others have termed the torsion modulus [23], or 'the shear stress that would develop on the periphery of a cylindrical specimen having a length and diameter = 1 when twisted round an angle of 1 radian' [24] in static methods of measuring the elastic coefficients of single crystals. In general, if a twisting couple is applied to a cylinder of a single crystal of an anisotropic material, bending can also take place in addition to torsion: these are conditions under which the free rigidity modulus, G F , is defined by Hearmon through the equation for a specimen of circular cross-section having its axis along 3 . In full tensor notation the equations for s 44 The formula in Eq. (70) is the formula given in Eq. 10/1 of [24]. This formula is also used by Date and Andrews [25] in their consideration of what they describe as the shear or torsion modulus of a cylindrical bar cut from a single crystal. The advantage of the torsion modulus defined through either Eq. (69) or Eq. (70) is that it depends only on the direction defining the axis of the cylinder, and therefore only on the direction cosines of the cylinder axis with respect to the crystallographic axes of the cubic material under consideration; it is independent of the perpendicular directions of 1 and 2 . This independence of the perpendicular directions of 1 and 2 also holds for the term 1 2 c 44 + c 55 , and so it is evident that 1 2 (c 44 + c 55 ), 1 2 (s 44 + s 55 ) and the torsion modulus for cubic materials can all be represented by a single surface in three dimensions in which the radius vector is in the direction [hkl] under consideration and the length of which is proportional to the term under consideration. This is the principle behind the representation of the torsion modulus in three dimensions in Fig. 4 of [23], reproduced in Fig. 154(b) of [24]. Other representations of 1 2 (s 44 + s 55 ) can be found in §373 of [5] and in Wooster [26]; Wooster uses a different definition of the s ij (see the footnote on p. 237 in Wooster [26]), so that the representations he shows are of what he defines as 2(s 44 + s 55 ) in his notation instead of 1 2 (s 44 + s 55 ) in the Voigt notation.
The clear advantage of these representations is that they can all be represented by a single surface in three dimensions. The clear disadvantage is that they average out differences between moduli in different directions within a plane that are evident in Figs. 1-5.

Discussion and Conclusions
The elastic stiffness coefficients and the shear moduli on high symmetry planes in high symmetry directions of cubic materials are summarised in Table 3. While for cubic materials that are relatively weakly anisotropic, such as Si, it is reasonable to presume that the dependence of c 44 and G 4 on the choice of shear plane and shear direction are very similar, a consideration of the tensor algebra and graphical representations relevant to these two quantities shows that there are clear differences which emerge for cubic materials that are highly anisotropic (A > 3 or A < 1/3) in their elastic coefficients. For materials of high anisotropy, particular differences between c 44 and G 4 are evident on {111} and {110} planes, the former in terms of the magnitudes of c 44 and G 4 , both of which are isotropic on {111} planes, the latter in terms of the differences between c 44 and G 4 that manifest themselves on directions other than the perpendicular 001 and 110 directions within the {110} planes.