International
Tables for Crystallography Volume D Physical properties of crystals Edited by A. Authier © International Union of Crystallography 2013 
International Tables for Crystallography (2013). Vol. D, ch. 1.3, pp. 72100
doi: 10.1107/97809553602060000902 Chapter 1.3. Elastic properties^{a}Institut de Minéralogie et de Physique des Milieux Condensés, 4 Place Jussieu, 75005 Paris, France, and ^{b}Laboratoire de Physique des Milieux Condensés, Université P. et M. Curie, 75252 Paris CEDEX 05, France In this chapter, the strain and stress tensors are defined and their main properties are derived. The elastic tensors, elastic stiffnesses and elastic compliances are then introduced. Their variation with orientation, depending on the crystal class, is given in the case of Young's modulus. The next part is devoted to the propagation of waves in continuous media (linear dynamic elasticity and the Christoffel matrix); the relation between the velocity and the elastic constants is given for the cubic, hexagonal and tetragonal classes. The experimental determination of elastic constants and their pressure and temperature dependence are discussed in separate sections. The last two sections of the chapter concerns nonlinear elasticity (second and higherorder elastic constants) and nonlinear dynamical elasticity. 
Let us consider a medium that undergoes a deformation. This means that the various points of the medium are displaced with respect to one another. Geometrical transformations of the medium that reduce to a translation of the medium as a whole will therefore not be considered. We may then suppose that there is an invariant point, O, whose position one can always return to by a suitable translation. A point P, with position vector , is displaced to the neighbouring point P′ by the deformation defined by The displacement vector constitutes a vector field. It is not a uniform field, unless the deformation reduces to a translation of the whole body, which is incompatible with the hypothesis that the medium undergoes a deformation. Let Q be a point that is near P before the deformation (Fig. 1.3.1.1). Then one can write
After the deformation, Q is displaced to Q′ defined by
In a deformation, it is more interesting in general to analyse the local, or relative, deformation than the absolute displacement. The relative displacement is given by comparing the vectors and PQ. Thus, one has Let us set Replacing by its expansion up to the first term gives
If we assume the Einstein convention (see Section 1.1.2.1 ), there is summation over j in (1.3.1.2) and (1.3.1.3). We shall further assume orthonormal coordinates throughout Chapter 1.3; variance is therefore not apparent and the positions of the indices have no meaning; the Einstein convention then only assumes repetition of a dummy index. The elements and are the components of dr and , respectively. Let us put where represents the Kronecker symbol; the 's are the components of matrix unity, I. The expressions (1.3.1.2) can also be written using matrices M and B: The components of the tensor are nonzero, unless, as mentioned earlier, the deformation reduces to a simple translation. Two cases in particular are of interest and will be discussed in turn:
If the components are constants, equations (1.3.1.3) can be integrated directly. They become, to a translation,
The fundamental property of the homogeneous deformation results from the fact that equations (1.3.1.4) are linear: a plane before the deformation remains a plane afterwards, a crystal lattice remains a lattice. Thermal expansion is a homogeneous deformation (see Chapter 1.4 ).
Some crystals present a twin microstructure that is seen to change when the crystals are gently squeezed. At rest, the domains can have one of two different possible orientations and the influence of an applied stress is to switch them from one orientation to the other. If one measures the shape of the crystal lattice (the strain of the lattice) as a function of the applied stress, one obtains an elastic hysteresis loop analogous to the magnetic or electric hysteresis loops observed in ferromagnetic or ferroelectric crystals. For this reason, these materials are called ferroelastic (see Chapters 3.1 to 3.3 and Salje, 1990). The strain associated with one of the two possible shapes of the crystal when no stress is applied is called the macroscopic spontaneous strain.
Let be the basis vectors before deformation. On account of the deformation, they are transformed into the three vectors The parallelepiped formed by these three vectors has a volume V′ given by where is the determinant associated with matrix B, V is the volume before deformation and represents a triple scalar product.
The relative variation of the volume is It is what one calls the cubic dilatation. gives directly the volume of the parallelepiped that is formed from the three vectors obtained in the deformation when starting from vectors forming an orthonormal base.
1.3.1.2.4. Expression of any homogeneous deformation as the product of a pure rotation and a pure deformation
Let us project the displacement vector on the position vector OP (Fig. 1.3.1.2), and let be this projection. The elongation is the quantity defined by where , , are the components of r. The elongation is the relative variation of the length of the vector r in the deformation. Let A and S be the antisymmetric and symmetric parts of M, respectively:
Only the symmetric part of M occurs in the expression of the elongation:
The geometrical study of the elongation as a function of the direction of r is facilitated by introducing the quadric associated with M: where is a constant. This quadric is called the quadric of elongations, Q. S is a symmetric matrix with three real orthogonal eigenvectors and three real eigenvalues, , , . If it is referred to these axes, equation (1.3.1.7) is reduced to
One can discuss the form of the quadric according to the sign of the eigenvalues :
In order to follow the variations of the elongation with the orientation of the position vector, one associates with r a vector y, which is parallel to it and is defined by where k is a constant. It can be seen that, in accordance with (1.3.1.6) and (1.3.1.7), the expression of the elongation in terms of y is
Thus, the elongation is inversely proportional to the square of the radius vector of the quadric of elongations parallel to OP. In practice, it is necessary to look for the intersection p of the parallel to OP drawn from the centre O of the quadric of elongations (Fig. 1.3.1.3a):

Equally, one can connect the displacement vector directly with the quadric Q. Using the bilinear form the gradient of , , has as components
One recognizes the components of the displacement vector u, which is therefore parallel to the normal to the quadric Q at the extremity of the radius vector Op parallel to r.
The directions of the principal axes of Q correspond to the extremal values of y, i.e. to the stationary values (maximal or minimal) of the elongation. These values are the principal elongations.
If the deformation is a pure rotation Hence we have
The quadric Q is a cylinder of revolution having the axis of rotation as axis.
If the deformation is small but arbitrary, i.e. if the products of two or more components of can be neglected with respect to unity, one can describe the deformation locally as a homogeneous asymptotic deformation. As was shown in Section 1.3.1.2.4, it can be put in the form of the product of a pure deformation corresponding to the symmetric part of , , and a pure rotation corresponding to the asymmetric part, : Matrix B can be written where I is the matrix identity. As the coefficients of are small, one can neglect the product and one has is a symmetric matrix that represents a pure deformation. is an antisymmetric unitary matrix and, since A is small, Thus, represents a rotation. The axis of rotation is parallel to the vector with coordinates which is an eigenvector of . The magnitude of the rotation is equal to the modulus of this vector.
In general, one is only interested in the pure deformation, i.e. in the form of the deformed object. Thus, one only wishes to know the quantities and the symmetric part of M. It is this symmetric part that is called the deformation tensor or the strain tensor. It is very convenient for applications to use the simplified notation due to Voigt:One may note that The Voigt strain matrix S is of the form
Let us consider an orthonormal system of axes with centre P. We remove nothing from the generality of the following by limiting ourselves to a planar problem and assuming that point P′ to which P goes in the deformation lies in the plane (Fig. 1.3.1.4). Let us consider two neighbouring points, Q and R, lying on axes and , respectively (, ). In the deformation, they go to points Q′ and R′ defined by

Geometrical interpretation of the components of the strain tensor. , , : axes before deformation; , , : axes after deformation. 
As the coefficients are small, the lengths of P′Q′ and P′R′ are hardly different from PQ and PR, respectively, and the elongations in the directions and are
The components , , of the principal diagonal of the Voigt matrix can then be interpreted as the elongations in the three directions , and . The angles α and β between PQ and P′Q′, and PR and P′R′, respectively, are given in the same way by One sees that the coefficient of Voigt's matrix is therefore The angle is equal to the difference between angles before deformation and after deformation. The nondiagonal terms of the Voigt matrix therefore represent the shears in the planes parallel to , and , respectively.
To summarize, if one considers a small cube before deformation, it becomes after deformation an arbitrary parallelepiped; the relative elongations of the three sides are given by the diagonal terms of the strain tensor and the variation of the angles by its nondiagonal terms.
The cubic dilatation (1.3.1.5) is (taking into account the fact that the coefficients are small).
Matrix M has only one coefficient, , and reduces to (Fig. 1.3.1.5a) The quadric of elongations is reduced to two parallel planes, perpendicular to , with the equation .
Let us consider a solid C, in movement or not, with a mass distribution defined by a specific mass ρ at each point. There are two types of force that are manifested in the interior of this solid.

Now consider a volume V within the solid C and the surface S which surrounds it (Fig. 1.3.2.2). Among the influences that are exterior to V, we distinguish those that are external to the solid C and those that are internal. The first are translated by the body forces, eventually by volume couples. The second are translated by the local contact forces of the part external to V on the internal part; they are represented by a surface density of forces, i.e. by the stresses that depend only on the point Q of the surface S where they are applied and on the orientation of the normal n of this surface at this point. If two surfaces S and S′ are tangents at the same point Q, the same stress acts at the point of contact between them. The equilibrium of the volume V requires:
The equilibrium of the solid C requires that:
Using the condition on the resultant of forces, it is possible to show that the components of the stress can be determined from the knowledge of the orientation of the normal n and of the components of a ranktwo tensor. Let P be a point situated inside volume V, , and three orthonormal axes, and consider a plane of arbitrary orientation that cuts the three axes at Q, R and S, respectively (Fig. 1.3.2.3). The small volume element PQRS is limited by four surfaces to which stresses are applied. The normals to the surfaces PRS, PSQ and PQR will be assumed to be directed towards the interior of the small volume. By contrast, for reasons that will become apparent later, the normal n applied to the surface QRS will be oriented towards the exterior. The corresponding applied forces are thus given in Table 1.3.2.1. The volume PQRS is subjected to five forces: the forces applied to each surface and the resultant of the volume forces and the inertial forces. The equilibrium of the small volume requires that the resultant of these forces be equal to zero and one can write (including the inertial forces in the volume forces).

As long as the surface element dσ is finite, however small, it is possible to divide both terms of the equation by it. If one introduces the direction cosines, , the equation becomes When dσ tends to zero, the ratio tends towards zero at the same time and may be neglected. The relation then becomes This relation is called the Cauchy relation, which allows the stress to be expressed as a function of the stresses , and that are applied to the three faces perpendicular to the axes, , and . Let us project this relation onto these three axes: The nine components are, by definition, the components of the stress tensor. In order to check that they are indeed the components of a tensor, it suffices to make the contracted product of each side of (1.3.2.4) by any vector : the lefthand side is a scalar product and the righthand side a bilinear form. The 's are therefore the components of a tensor. The index to the far left indicates the face to which the stress is applied (normal to the , or axis), while the second one indicates on which axis the stress is projected.
Let us return to equation (1.3.2.1) expressing the equilibrium condition for the resultant of the forces. By replacing by the expression (1.3.2.4), we get, after projection on the three axes, where and the inertial forces are included in the volume forces. Applying Green's theorem to the first integral, we have
The equilibrium condition now becomes In order that this relation applies to any volume V, the expression under the integral must be equal to zero, or, if one includes explicitly the inertial forces, This is the condition of continuity or of conservation. It expresses how constraints propagate throughout the solid. This is how the cohesion of the solid is ensured. The resolution of any elastic problem requires solving this equation in terms of the particular boundary conditions of that problem.
Let us now consider the equilibrium condition (1.3.2.2) relative to the resultant moment. After projection on the three axes, and using the Cartesian expression (1.1.3.4) of the vectorial products, we obtain (including the inertial forces in the volume forces). is the permutation tensor. Applying Green's theorem to the first integral and putting the two terms together gives
In order that this relation applies to any volume V within the solid C, we must have or
Taking into account the continuity condition (1.3.2.5), this equation reduces to
A volume couple can occur for instance in the case of a magnetic or an electric field acting on a body that locally possesses magnetic or electric moments. In general, apart from very rare cases, one can ignore these volume couples. One can then deduce that the stress tensor is symmetrical:
This result can be recovered by applying the relation (1.3.2.2) to a small volume in the form of an elementary parallelepiped, thus illustrating the demonstration using Green's theorem but giving insight into the action of the constraints. Consider a rectangular parallelepiped, of sides , and , with centre P at the origin of an orthonormal system whose axes , and are normal to the sides of the parallelepiped (Fig. 1.3.2.4). In order that the resultant moment with respect to a point be zero, it is necessary that the resultant moments with respect to three axes concurrent in this point are zero. Let us write for instance that the resultant moment with respect to the axis is zero. We note that the constraints applied to the faces perpendicular to do not give rise to a moment and neither do the components , , and of the constraints applied to the faces normal to and (Fig. 1.3.2.4). The components and alone have a nonzero moment.

Symmetry of the stress tensor: the moments of the couples applied to a parallelepiped compensate each other. 
For face 1, the constraint is if is the magnitude of the constraint at P. The force applied at face 1 is and its moment is
Similarly, the moments of the force on the other faces are
Noting further that the moments applied to the faces 1 and 1′ are of the same sense, and that those applied to faces 2 and 2′ are of the opposite sense, we can state that the resultant moment is where is the volume of the small parallelepiped. The resultant moment per unit volume, taking into account the couples in volume, is therefore It must equal zero and the relation given above is thus recovered.
We shall use frequently the notation due to Voigt (1910) in order to express the components of the stress tensor: It should be noted that the conventions are different for the Voigt matrices associated with the stress tensor and with the strain tensor (Section 1.3.1.3.1).
The Voigt matrix associated with the stress tensor is therefore of the form
1.3.2.5.2. Interpretation of the components of the stress tensor – special forms of the stress tensor

If the surface of the solid C is free from all exterior action and is in equilibrium, the stress field inside C is zero at the surface. If C is subjected from the outside to a distribution of stresses (apart from the volume forces mentioned earlier), the stress field inside the solid is such that at each point of the surface where the 's are the direction cosines of the normal to the surface at the point under consideration.

Consider a medium that is subjected to a stress field . It has sustained a deformation indicated by the deformation tensor S. During this deformation, the forces of contact have performed work and the medium has accumulated a certain elastic energy W. The knowledge of the energy density thus acquired is useful for studying the properties of the elastic constants. Let the medium deform from the deformation to the deformation under the influence of the stress field and let us evaluate the work of each component of the effort. Consider a small elementary rectangular parallelepiped of sides , , (Fig. 1.3.2.8). We shall limit our calculation to the components and , which are applied to the faces 1 and 1′, respectively.

Determination of the energy density in a deformed medium. PP′ represents the displacement of the small parallelepiped during the deformation. The thick arrows represent the forces applied to the faces 1 and 1′. 
In the deformation , the point P goes to the point P′, defined by A neighbouring point Q goes to Q′ such that (Fig. 1.3.1.1) The coordinates of are given by
Of sole importance is the relative displacement of Q with respect to P and the displacement that must be taken into account in calculating the forces applied at Q. The coordinates of the relative displacement are
We shall take as the position of Q the point of application of the forces at face 1, i.e. its centre with coordinates (Fig. 1.3.2.8). The area of face 1 is and the forces arising from the stresses and are equal to and , respectively. The relative displacement of Q parallel to the line of action of is and the corresponding displacement along the line of action of is . The work of the corresponding forces is therefore
The work of the forces applied to the face 1′ is the same (, and change sign simultaneously). The works corresponding to the faces 1 and 1′ are thus and for the two stresses, respectively. One finds an analogous result for each of the other components of the stress tensor and the total work per unit volume is
Let us consider a metallic bar of length loaded in pure tension (Fig. 1.3.3.1). Under the action of the uniaxial stress (F applied force, area of the section of the bar), the bar elongates and its length becomes . Fig. 1.3.3.2 relates the variations of Δl and of the applied stress T. The curve representing the traction is very schematic and does not correspond to any real case. The following result, however, is common to all concrete situations:

Young's modulus is not sufficient to describe the deformation of the bar: its diameter is reduced, in effect, during the elongation. One other coefficient, at least, is therefore necessary. In a general way, let us consider the deformation of a continuous anisotropic medium under the action of a field of applied stresses. We will generalize Hooke's law by writing that at each point there is a linear relation between the components of the stress tensor and the components of the strain tensor: The quantities and are characteristic of the elastic properties of the medium if it is homogeneous and are independent of the point under consideration. Their tensorial nature can be shown using the demonstration illustrated in Section 1.1.3.4 . Let us take the contracted product of the two sides of each of the two equations of (1.3.3.2) by the components and of any two vectors, x and y: The lefthand sides are bilinear forms since and are secondrank tensors and the righthand sides are quadrilinear forms, which shows that and are the components of fourthrank tensors, the tensor of elastic compliances (or moduli) and the tensor of elastic stiffnesses (or coefficients), respectively. The number of their components is equal to 81.
Equations (1.3.3.2) are Taylor expansions limited to the first term. The higher terms involve sixthrank tensors, and , with coefficients, called thirdorder elastic compliances and stiffnesses and eighthrank tensors with coefficients, called fourthorder elastic compliances and stiffnesses. They will be defined in Section 1.3.6.4. Tables for thirdorder elastic constants are given in Fumi (1951, 1952, 1987). The accompanying software to this volume enables these tables to be derived for any point group.
It is convenient to write the relations (1.3.3.2) in matrix form by associating with the stress and strain tensors column matrices and with the tensors of the elastic stiffnesses, c, and of the elastic compliances, s, square matrices (Section 1.1.4.10.4 ); these two matrices are inverse to one another. The number of independent components of the fourthrank elastic tensors can be reduced by three types of consideration:

Expression (1.3.2.7) of the strain energy stored per unit volume in a medium for a small deformation can be integrated when the medium is strained under a stress according to linear elasticity. Applying relation (1.3.3.2), one gets for the density of strain energy
Let us apply a hydrostatic pressure (Section 1.3.2.5.2). The medium undergoes a relative variation of volume (the cubic dilatation, Section 1.3.1.3.2). If one replaces in (1.3.3.8) the stress distribution by a hydrostatic pressure, one obtains for the components of the strain tensorFrom this, we deduce the volume compressibility, , which is the inverse of the bulk modulus, κ: This expression reduces for a cubic or isotropic medium to
Under the action of a hydrostatic pressure, each vector assumes a different elongation. This elongation is given by equation (1.3.1.6): where the 's are the direction cosines of r. The coefficient of linear compressibility is, by definition, . Replacing by its value , we obtain for the coefficient of linear compressibility In the case of a cubic or isotropic medium, this expression reduces to
The coefficient of linear compressibility is then equal to one third of the coefficient of volume compressibility. We note that the quadric of elongations is a sphere.
If the applied stress reduces to a uniaxial stress, , the strain tensor is of the form In particular, We deduce from this that Young's modulus (equation 1.3.3.1) is
The elongation of a bar under the action of a uniaxial stress is characterized by and the diminution of the cross section is characterized by and . For a cubic material, the relative diminution of the diameter is One deduces from this that is necessarily of opposite sign to and one calls the ratio Poisson's ratio.
Putting this value into expression (1.3.3.12) for the coefficient of compressibility in cubic or isotropic materials gives
As the coefficient of compressibility, by definition, is always positive, we have
In practice, Poisson's ratio is always close to 0.3. It is a dimensionless number. The quantity represents the departure from isotropy of the material and is the anisotropy factor. It is to be noted that cubic materials are not isotropic for elastic properties. Table 1.3.3.2 gives the values of , , , ν and for a few cubic materials.
It is interesting to calculate Young's modulus in any direction. For this it is sufficient to change the axes of the tensor . If A is the matrix associated with the change of axes, leading to the direction changing to the direction , then Young's modulus in this new direction is with The matrix coefficients are the direction cosines of with respect to the axes , and . In spherical coordinates, they are given by (Fig. 1.3.3.3) where θ is the angle between and , and ϕ is the angle between and . Using the reduction of for the various crystal classes (Section 1.1.4.9.9 ), we find, in terms of the reduced twoindex components, the following.
The representation surface of , the inverse of Young's modulus, is illustrated in Figure 1.3.3.4 for crystals of different symmetries. As predicted by the Neumann principle, the representation surface is invariant with respect to the symmetry elements of the point group of the crystal but, as stated by the Curie laws, its symmetry can be larger. In the examples of Fig. 1.3.3.4, the symmetry of the surface is the same as that of the point group for aluminium (Fig. 1.3.3.4a), tungsten (Fig. 1.3.3.4b) and sodium chloride (Fig. 1.3.3.4c), which have as point group, for tin (Fig. 1.3.3.4e, ) and for calcite (Fig. 1.3.3.4f, ). But in the case of zinc (Fig. 1.3.3.4d, ), the surface is of revolution and has a larger symmetry. It is interesting to compare the differences in shapes of the representation surfaces for the three cubic crystals, depending on the value of the anisotropy factor, which is larger than 1 for aluminium, smaller than 1 for sodium chloride and close to 1 for tungsten (see Table 1.3.3.2). In this latter case, the crystal is pseudoisotropic and the surface is practically a sphere.
The isotropy relation between elastic compliances and elastic stiffnesses is given in Section 1.3.3.2.3. For reasons of symmetry, the directions of the eigenvectors of the stress and strain tensors are necessarily the same in an isotropic medium. If we take these directions as axes, the two tensors are automatically diagonalized and the second relation (1.3.3.7) becomes
These relations can equally well be written in the symmetrical form
If one introduces the Lamé constants, the equations may be written in the form often used in mechanics:
Two coefficients suffice to define the elastic properties of an isotropic material, and , and , μ and λ, μ and ν, etc. Table 1.3.3.3 gives the relations between the more common elastic coefficients.

We saw in Section 1.3.2.3 that the condition of equilibrium is
If we use the relations of elasticity, equation (1.3.3.2), this condition can be rewritten as a condition on the components of the strain tensor: Recalling that the condition becomes a condition on the displacement vector, : In an isotropic orthonormal medium, this equation, projected on the axis , can be written with the aid of relations (1.3.3.5) and (1.3.3.9): This equation can finally be rearranged in one of the three following forms with the aid of Table 1.3.3.3.
The elastic properties of materials have been considered in the preceding section in the static state and the elastic constants have been defined in terms of the response of the material to particular static forces. It is effectively the way the elastic constants have been measured in the past, although the measurements could not be very precise. A way of proceeding frequently used now is to excite a mechanical wave in the crystal and measure its propagation velocity or the wavelength associated with a particular frequency. One method consists in sending a train of ultrasonic waves through the crystal; one uses a pulse generator and a piezoelectric transducer glued to the crystal. The elapsed time between the emission of the train of waves and its reception after reflection from the rear face of the sample is then measured. Another method involves producing a system of standing waves after reflection at the inner surface of the crystal and determining the set of resonance frequencies. The experimental techniques will be described in Section 1.3.4.6.
The purpose of the next sections is to establish relations between the wavelength – or the velocity of propagation – and the elastic constants.
Consider the propagation of a wave in a continuous medium. The elongation of each point will be of the form where ν is the frequency and q is the wavevector. The velocity of propagation of the wave is
We saw in Section 1.3.3.6 that the equilibrium condition is Here the only volume forces that we must consider are the inertial forces:
The position vector of the point under consideration is of the form where only depends on the time and defines the mean position. Equation (1.3.4.3) is written therefore
Replacing u by its value in (1.3.4.1), dividing by and using orthonormal coordinates, we get
It can be seen that, for a given wavevector, appears as an eigenvalue of the matrix of which the vibration vector u is an eigenvector. This matrix is called the dynamical matrix, or Christoffel matrix. In order that the system (1.3.4.5) has a solution other than a trivial one, it is necessary that the associated determinant be equal to zero. It is called the Christoffel determinant and it plays a fundamental role in the study of the propagation of elastic waves in crystals.
Let be the direction cosines of the wavevector q. The components of the wavevector are With this relation and (1.3.4.2), the system (1.3.4.5) becomes Putting in (1.3.4.6), the condition that the Christoffel determinant is zero can be written
On account of the intrinsic symmetry of the tensor of elastic stiffnesses, the matrix is symmetrical.
If we introduce into expression (1.3.4.7) the elastic stiffnesses with two indices [equation (1.3.3.6)], we find, for instance, for and
The expression for the effective value, , of the `stiffened' elastic stiffness in the case of piezoelectric crystals is given in Section 2.4.2.2 .
Equation (1.3.4.7) may be written
This shows that in a dynamic process only the sums can be measured and not and separately. On the contrary, can be measured directly. In the cubic system therefore, for instance, is determined from the measurement of on the one hand and from that of on the other hand.
The Christoffel determinant has three roots and the Christoffel matrix, being Hermitian with real coefficients, has three real eigenvalues and three orthogonal eigenvectors. The wavevector q, therefore, encompasses three waves with vibration vectors , , which are perpendicular to one another. In the general case, there is no particular angular relationship between the vibration vectors (or polarization vectors). However, if the latter are parallel to certain symmetry directions in the crystal, one of the vibration vectors is along this direction. The corresponding wave is called longitudinal. The two other waves have their polarization direction perpendicular to the wavevector and are thus transverse. If one of the polarization vectors is almost parallel to the wavevector, which often happens, then one speaks of the vibration as being quasilongitudinal.
We shall limit ourselves to cubic, hexagonal and tetragonal crystals and consider particular cases.

In hexagonal crystals, there are five independent elastic stiffnesses, , , , , and (Section 1.1.4.10.4 ).

In tetragonal crystals, there are six independent elastic stiffnesses, , , , , and (Section 1.1.4.10.4 ).

As mentioned in Section 1.3.4.1, the elastic constants of a material can be obtained by the elastic response of the material to particular static forces; however, such measurements are not precise and the most often used approach nowadays consists of determining the velocity of ultrasonic waves propagating along different directions of the crystal and calculating the elastic constants from the Christoffel determinants (1.3.4.8). The experimental values are often accurate enough to justify the distinction between static and dynamic values of the elastic constants and between phase and group velocities, and the careful consideration of the frequency range of the experiments.

The use of the resonance technique is a well established approach for determining the velocity of sound in a gas by observing nodes and antinodes of a system of standing waves produced in the socalled Kund tube. In the case of transparent solids, optical means allow us to visualize the standing waves and to measure the wavelength directly (Zarembowitch, 1965). An easier procedure can be used: let us consider a transparent crystal in the shape of a parallelepiped (Fig. 1.3.4.1). A piezoelectric transducer is glued to the crystal and excited at varying frequencies. If the bonding between the transducer and the crystal is loose enough, the crystal can be considered as free from stress and the sequence of its resonance frequencies is given bywhere n is an integer, V the phase velocity of the wave in the direction orthogonal to the parallel faces and l the distance between these faces.
The looseness of the bonding can be checked by the regularity of the arithmetic ratio, . On account of the elastooptic coupling, a phase grating is associated with the elastic standingwave system and a light beam can be diffracted by this grating. The intensity of the diffraction pattern is maximum when resonance occurs. A large number of resonance frequencies can be detected, usually more than 100, sometimes 1000 for nonattenuating materials. Consequently, in favourable cases the absolute value of the ultrasonic velocity can be determined with an uncertainty less than .
Pulseecho techniques are valid for transparent and opaque materials. They are currently used for measuring ultrasonic velocities in solids and can be used in very simple as well as in sophisticated versions according to the required precision (McSkimmin, 1964). In the simplest version (Fig. 1.3.4.2), an electronic pulse generator excites the mechanical vibrations of a piezoelectric transducer glued to one of two planeparallel faces of a specimen. An ultrasonic pulse whose duration is of the order of a microsecond is generated and transmitted through the specimen. After reflection at the opposite face, it returns and, when it arrives back at the transducer, it gives rise to an electronic signal, or echo. The whole sequence of such echos is displayed on the screen of an oscilloscope and it is possible to measure from them the time interval for transit. Usually, Xcut quartz crystals or ferroelectric ceramics are used to excite longitudinal waves and Ycut quartz is used to excite transverse waves. In many cases, a circulator, or gate, is used to protect the receiver from saturation following the main `bang'. This method is rough because the beginning and the end of a pulse are not well characterized. Several improvements have therefore been made, mainly based on interferometric techniques (pulsesuperposition method, `sing around' method etc.). Nevertheless, if the absolute value of the ultrasonic velocity is not determined with a high accuracy by using pulseecho techniques, this approach has proved valuable when relative values of ultrasonic velocities are needed, e.g. temperature and pressure dependences of ultrasonic velocities.

In a solid, the elastic constants are temperature and pressure dependent. As examples, the temperature dependence of the elastic stiffnesses of an aluminium single crystal within its stability domain (the melting point is 933 K) and the pressure dependence of the elastic stiffnesses of the ternary compound KZnF_{3} within its stability domain (the crystal becomes unstable for a hydrostatic pressure of about 20 GPa) are shown in Figs. 1.3.5.1 and 1.3.5.2, respectively.

Temperature dependence of the elastic stiffnesses of an aluminium single crystal (after Every & McCurdy, 1992). 

Pressure dependence of the elastic stiffness of a KZnF_{3} crystal. Reproduced with permission from Ultrasonics Symposium Proc. IEEE (Fischer et al., 1980). Copyright (1980) IEEE. 
We can observe the following trends, which are general for stable crystals:

These observations can be quantitatively justified on the basis of an equation of state of a solid: where represents the stress tensor, the strain tensor, X the position of the elementary elements of the solid and Θ the temperature.
Different equations of state of solids have been proposed. They correspond to different degrees of approximation that can only be discussed and understood in a microscopic theory of lattice dynamics. The different steps in the development of lattice dynamics, the Einstein model, the Debye model and the Grüneisen model, will be presented in Section 2.1.2.7 . Concerning the temperature and the pressure dependences of the elastic constants, we may notice that rather sophisticated models are needed to describe correctly the general trends mentioned above:

Table 1.3.5.1 gives typical values of for some cubic crystals considered within their stability domain. In column 6, the `elastic Debye temperature' of the crystal, , has been calculated according to the formula where h is the Planck constant, is the Boltzmann constant, v is an average velocity (see for instance De Launay, 1956) and n is the number of atoms per unit volume.

It is interesting to compare , the `elastic Debye temperature', with , the `calorimetric Debye temperature'. The definition of will be given in Section 2.1.2.7 . It results from the attempt at founding a universal description for the thermal properties of solids when the temperature is expressed as a reduced temperature, ; is obtained from calorimetric measurements at low temperature. It is worth noting that accurate values of lowtemperature elastic constants and lowtemperature calorimetric measurements lead to an excellent agreement between and [better than 2 or 3 K (De Launay, 1956)]. This agreement demonstrates the validity of the Debye model in the vicinity of 0 K. From Table 1.3.5.1, we can observe that for ionic crystals is, in general, greater than . This remark is not valid for covalent and metallic crystals. Typical orders of magnitude are given in Table 1.3.5.2. These statements concern only general trends valid for stable crystals.

In the case of temperatureinduced phase transitions, some elastic constants are softened in the vicinity and sometimes far from the critical temperature. As an example, Fig. 1.3.5.3 shows the temperature dependence of in RbCdF_{3}, CsCdF_{3} and TlCdF_{3} single crystals. RbCdF_{3} and TlCdF_{3} undergo structural phase transitions at 124 and 191 K, respectively, while CsCdF_{3} remains stable in this temperature range. The softening of when the temperature decreases starts more than 100 K before the critical temperature, . In contrast, Fig. 1.3.5.4 shows the temperature dependence of in KNiF_{3}, a crystal that undergoes a para–antiferromagnetic phase transition at 246 K; the coupling between the elastic and the magnetic energy is weak, consequently decreases abruptly only a few degrees before the critical temperature. We can generalize this observation and state that the softening of an elastic constant occurs over a large domain of temperature when this constant is the order parameter or is strongly coupled to the order parameter of the transformation; for instance, in the cooperative Jahn–Teller phase transition in DyVO_{4}, is the soft acoustic phonon mode leading to the phase transition and this parameter anticipates the phase transition 300 K before it occurs (Fig. 1.3.5.5).

Temperature dependence of the elastic constant in RbCdF_{3}, CsCdF_{3} and TlCdF_{3} crystals; the crystals of RbCdF_{3} and TlCdF_{3} undergo structural phase transitions (after Rousseau et al., 1975). 

Temperature dependence of the elastic constant in KNiF_{3}, which undergoes a para–antiferromagnetic phase transition. Reprinted with permission from Appl. Phys. Lett. (Nouet et al., 1972). Copyright (1972) American Institute of Physics. 
As mentioned above, anharmonic potentials are needed to explain the stress dependence of the elastic constants of a crystal. Thus, if the strainenergy density is developed in a polynomial in terms of the strain, only the first and the second elastic constants are used in linear elasticity (harmonic potentials), whereas higherorder elastic constants are also needed for nonlinear elasticity (anharmonic potentials).
Concerning the pressure dependence of the elastic constants (nonlinear elastic effect), considerable attention has been paid to their experimental determination since they are a unique source of significant information in many fields:

In a solid body, the relation between the stress tensor T and the strain tensor S is usually described by Hooke's law, which postulates linear relations between the components of T and S (Section 1.3.3.1). Such relations can be summarized by (see equation 1.3.3.2)where the 's are the elastic stiffnesses.
For a solid under finite strain conditions, Hooke's law, valid for infinitesimal deformations, does not hold, and the fundamental definitions for stress and strain must be revisited.
Finite elastic strains may be treated from two different viewpoints using either the Lagrangian (material) or the Eulerian (spatial) descriptions.
Let us consider a fixed rectangular Cartesian coordinate system with axes (). Any particular position vector r of components (, , ) denotes a point in space. A point that always moves with the material is called a particle or material point. Let every particle be identified by its coordinates at some reference time . These reference coordinates, referred to the same Cartesian system, will be denoted by (, , ) and the corresponding position vector a. A particular vector a can serve as a name for the particle located at that position at the reference time .
The vectors r and a both specify a position in a fixed Cartesian frame of reference. At any time, we associate each r with an a by the rule that r is the present position vector of the particle initially at a. This connection between r and a is written symbolically as where
The coordinates that identify the particles are called material coordinates. A description that, like (1.3.6.1), uses (t, , , ) as independent variables is called a material or Lagrangian description.
The converse of (1.3.6.1) and (1.3.6.2) may be written where
A spatial description or Eulerian description uses the independent variables (t, , , ), the being called spatial coordinates.
Now, for the sake of simplicity, we shall work with the Lagrangian formulation exclusively. For more details see, for instance, Thurston (1964) and Wallace (1970, 1972).
The displacement vector from the reference position of a particle to its new position has as components
The term strain refers to a change in the relative positions of the material points in a body. Let a final configuration be described in terms of the reference configuration by setting t equal to a constant in (1.3.6.1). Then t no longer appears as a variable and (1.3.6.1) can be written where the are the independent variables. It follows that
Let now the particle initially at (, , ) move to (, , ). The square of the initial distance to a neighbouring particle whose initial coordinates were is The square of the final distance to the same neighbouring particle is
In a material description, the strain components are defined by the following equations: Substituting (1.3.6.6) into (1.3.6.7), it follows that Hence
If the products and squares of the displacement derivatives are neglected, the strain components reduce to the usual form of `infinitesimal elasticity' [see equation (1.3.1.8)]:
It is often useful to introduce the Jacobian matrix associated with the transformation (a, x). The components of this matrix are where
From the definition of matrix J, one has and where , and are the transpose matrices of da, dx and J, respectively, and δ is the Kronecker matrix.
The Lagrangian strain matrix S may then be written symbolically:
When finite strains are concerned, we have to distinguish three states of the medium: the natural state, the initial state and the final or present state: The natural state is a state free of stress. The initial state is deduced from the natural state by a homogeneous strain. The final state is deduced from the initial state by an arbitrary strain.
Concerning the stress tensor, as pointed out by Thurston (1964), the stressdeformation relation is complicated in nonlinear elasticity because `the strain is often referred to a natural unstressed state, whereas the stress is defined per unit area of the deformed body'. For this reason, the differential of work done by the stress is not equal to the stress components times the differentials of the corresponding strain components. So, following Truesdell & Toupin (1960), we shall introduce a thermodynamic tension tensor defined as the first derivative of the energy with respect to strain. If the internal energy U per unit mass is considered, the thermodynamic tension refers to an isentropic process. Then where σ is the entropy and the volumic mass in the initial state.
If the Helmholtz free energy F is considered, the thermodynamic tension refers to an isothermal process. Thenwhere Θ is the temperature. It will be shown in Section 1.3.7.2 that
Following Brugger (1964), the strainenergy density, or strain energy per unit volume Φ, is assumed to be a polynomial in the strain: where , , X denotes the configuration of the initial state and the 's are the Lagrangian finite straintensor components.
If the initial energy and the deformation of the body are both zero, the first two terms in (1.3.6.9) are zero. Note that is a stress and not an intrinsic characteristic of the material. In this expression, the elastic stiffnesses and are the second and thirdorder stiffnesses, respectively. Since the strain tensor is symmetric, pairs of subscripts can be interchanged [see equation (1.3.3.4)]:
More accurately, the isentropic and the isothermal elastic stiffnesses are defined as the nth partial derivatives of the internal energy and the Helmholtz free energy, respectively. For example, the thirdorder isentropic and isothermal stiffnesses are, respectively, where the internal energy, U, is a function of X, and σ, and the Helmholtz free energy, F, is a function of X, and Θ.
From these definitions, it follows that the Brugger stiffness coefficients depend on the initial state. When no additional information is given, the initial state is the natural state.
The thirdorder stiffnesses form a sixthrank tensor containing 3^{6} = 729 components, of which 56 are independent for a triclinic crystal and 3 for isotropic materials. The independent components of a sixthrank tensor can be obtained for any point group using the accompanying software to this volume. They are listed in Table 1.3.6.1 for each Laue class. They can be reduced to threeindex coefficients. Table 1.3.6.2 gives the values of the thirdorder stiffnesses of some materials.


The three independent constants for isotropic materials are often taken as , and and denoted respectively by , , , the `thirdorder Lamé constants'.
The `thirdorder Murnaghan constants' (Murnaghan, 1951), denoted by , are given in terms of the Brugger constants by the relations
Similarly, the fourthorder stiffnesses form an eighthrank tensor containing components, 126 of which are independent for a triclinic crystal and 11 for isotropic materials (the independent components of a sixthrank tensor can be obtained for any point group using the accompanying software to this volume).
For a solid under finite strain conditions, the definition of the elastic compliance tensor has to be reconsidered. In linear elasticity, the secondorder elastic compliances were defined through the relations (1.3.3.2): while, in nonlinear elasticity, one has where
In most experiments, the initial stress is small compared with the secondorder elastic constants (for example, 1 GPa hydrostatic pressure compared with the usual value c_{ijkl} = 100 GPa). Consequently, the deformation between the initial (stressed) state and the natural (unstressed) state is small compared with 1. For this reason, it is convenient to expand the elastic constants in the initial state as a power series in the strain about the natural state. To avoid confusion, we introduce new notations: now represents the coordinates in the natural or unstressed state; X represents the coordinates in the initial or homogeneously strained state; are the components of displacement. All letters with superscript bar refer to the natural state; for example, denotes the Lagrangian strain in the natural state; .
Now, in order to relate the properties at X to those at , we need to specify the strain from to X. Let Consequently, The secondorder elastic constants at X can be expressed in terms of the second and thirdorder elastic constants at : or This expression holds for both isentropic and isothermal elastic constants.
The elastic strainenergy density has appeared in the literature in various forms. Most of the authors use the Murnaghan constants as long as isotropic solids are concerned. However, most of the literature uses Brugger's thermodynamic definition when anisotropic media are under consideration (Brugger, 1964).
The elastic strainenergy density for an isotropic medium, including thirdorder terms but omitting terms independent of strain, may be expressed in terms of three strain invariants, since an isotropic material is invariant with respect to rotation: where λ and μ are the secondorder Lamé constants, are the thirdorder Murnaghan constants, and , , are the three invariants of the Lagrangian strain matrix. These invariants may be written in terms of the strain components as The elastic strainenergy density for an anisotropic medium (for example a medium belonging to the most symmetrical groups of cubic crystals) is (Green, 1973)
In recent years, the measurements of ultrasonic wave velocities as functions of stresses applied to the sample and the measurements of the amplitude of harmonics generated by the passage of an ultrasonic wave throughout the sample are in current use. These experiments and others, such as the interaction of two ultrasonic waves, are interpreted from the same theoretical basis, namely nonlinear dynamical elasticity.
A first step in the development of nonlinear dynamical elasticity is the derivation of the general equations of motion for elastic waves propagating in a solid under nonlinear elastic conditions. Then, these equations are restricted to elastic waves propagating either in an isotropic or in a cubic medium. The next step is the examination of two important cases:
Finally, the concept of natural velocity is introduced and the experiments that can be used to determine the third and higherorder elastic constants are described.
For generality, these equations will be derived in the X configuration (initial state). It is convenient to obtain the equations of motion with the aid of Lagrange's equations. In the absence of body forces, these equations are or where L is the Lagrangian per unit initial volume and are the elements of the Jacobian matrix.
For adiabatic motionwhere U is the internal energy per unit mass.
Combining (1.3.7.2) and (1.3.7.3), it follows that which can be written since
Using now the equation of continuity or conservation of mass: and the identity of Euler, Piola and Jacobi: we get an expression of Newton's law of motion: with becomes since , the thermodynamic tensor conjugate to the variable , is generally denoted as the `second Piola–Kirchoff stress tensor'.
Using Φ, the strain energy per unit volume, Newton's law (1.3.7.4) takes the form and
As an example, let us consider the case of a plane finite amplitude wave propagating along the axis. The displacement components in this case become Thus, the Jacobian matrix reduces to
The Lagrangian strain matrix is [equation (1.3.6.8)] The only nonvanishing strain components are, therefore, and the strain invariants reduce to
In this case, the strainenergy density becomes Differentiating (1.3.7.6) with respect to the strains, we get All the other .
From (1.3.7.5), we derive the stress components: Note that this tensor is not symmetric.
For the particular problem discussed here, the three components of the equation of motion are
If we retain only terms up to the quadratic order in the displacement gradients, we obtain the following equations of motion:
In this case, the strainenergy density becomes Differentiating (1.3.7.8) with respect to the strain, one obtains All other . From (1.3.7.5), we derive the stress components: In this particular case, the three components of the equation of motion are
If we retain only terms up to the quadratic order in the displacement gradients, we obtain the following equations of motion: which are identical to (1.3.7.7) if we put
The coordinates in the medium free of stress are denoted either a or . The notation is used when we have to discriminate the natural configuration, , from the initial configuration X. Here, the process that we describe refers to the propagation of an elastic wave in a medium free of stress (natural state) and the coordinates will be denoted .
Let us first examine the case of a pure longitudinal mode, i.e.
The equations of motion, (1.3.7.7) and (1.3.7.9), reduce to for an isotropic medium or for a cubic crystal (most symmetrical groups) when a pure longitudinal mode is propagated along [100].
For both cases, we have a onedimensional problem; (1.3.7.7) and (1.3.7.9) can therefore be written
The same equation is also valid when a pure longitudinal mode is propagated along [110] and [111], with the following correspondence: Let us assume that ; a perturbation solution to (1.3.7.10) is where with
If we substitute the trial solutions into (1.3.7.10), we find after one iteration the following approximate solution: which involves secondharmonic generation.
If additional iterations are performed, higher harmonic terms will be obtained. A well known property of the firstorder nonlinear equation (1.3.7.10) is that its solutions exhibit discontinuous behaviour at some point in space and time. It can be seen that such a discontinuity would appear at a distance from the origin given by (Breazeale, 1984) where is the initial value for the particle velocity.
We now consider the propagation of smallamplitude elastic waves in a homogeneously strained medium. As defined previously, or a are the coordinates in the natural or unstressed state. X are the coordinates in the initial or homogeneously strained state. are the components of displacement from the initial state due to the wave.
Starting from (1.3.7.4), we get Its partial derivative is If we expand the state function about the initial configuration, it follows that
The linearized stress derivatives become If we let , the equation of motion in the initial state is The coefficients do not present the symmetry of the coefficients except in the natural state where and are equal.
The simplest solutions of the equation of motion are plane waves. We now assume plane sinusoidal waves of the form where k is the wavevector.
Substitution of (1.3.7.14) into (1.3.7.13) results in or with .
The quantities and A are, respectively, the eigenvalues and eigenvectors of the matrix . Since is a real symmetric matrix, the eigenvalues are real and the eigenvectors are orthogonal.
The main experimental procedures for determining the third and higherorder elastic constants are based on the measurement of stress derivatives of ultrasonic velocities and on harmonic generation experiments. Hydrostatic pressure, which can be accurately measured, has been widely used; however, the measurement of ultrasonic velocities in a solid under hydrostatic pressure cannot lead to the whole set of thirdorder elastic constants, so uniaxial stress measurements or harmonic generation experiments are then necessary.
In order to interpret wavepropagation measurements in stressed crystals, Thurston (1964) and Brugger (1964) introduced the concept of natural velocity with the following comments:
`According to equation of motion, the wave front is a material plane which has unit normal k in the natural state; a wave front moves from the plane to the plane in the time . Thus W, the natural velocity, is the wave speed referred to natural dimensions for propagation normal to a plane of natural normal k.
In a typical ultrasonic experiment, plane waves are reflected between opposite parallel faces of a specimen, the wave fronts being parallel to these faces. One ordinarily measures a repetition frequency F, which is the inverse of the time required for a round trip between the opposite faces.'
In most experiments, the thirdorder elastic constants and higherorder elastic constants are deduced from the stress derivatives of . For instance, Table 1.3.7.1 gives the expressions for and for a cubic crystal. These quantities refer to the natural state free of stress. In this table, p denotes the hydrostatic pressure and the 's are the following linear combinations of thirdorder elastic constants:


References
Breazeale, M. A. (1984). Determination of thirdorder elastic constants from ultrasonic harmonic generation. Physical Acoustics, Vol. 17, edited by R. N. Thurston, pp. 2–75. New York: Academic Press.Brillouin, L. (1932). Propagation des ondes électromagnétiques dans les milieux matériels. Congrès International d'Électricité, Vol. 2, Section 1, pp. 739–788. Paris: GauthierVillars.
Brugger, K. (1964). Thermodynamic definition of higherorder elastic coefficients. Phys. Rev. 133, 1611–1612.
De Launay, J. (1956). The theory of specific heats and lattice vibrations. Solid State Physics, Vol. 2, edited by F. Seitz & D. Turnbull, pp. 219–303. New York: Academic Press.
Every, A. G. & McCurdy, A. K. (1992). Second and higher order elastic constants. In LandoltBörnstein, Group III, Condensed Matter, Volume 29, Subvolume A. Berlin/Heidelberg/New York: Springer.
Fischer, M. (1982). Third and fourthorder elastic constants of fluoperovskites CsCdF_{3}, TlCdF_{3}, RbCdF_{3}, RbCaF_{3}. J. Phys. Chem. Solids, 43, 673–682.
Fischer, M., Zarembowitch, A. & Breazeale, M. A. (1980). Nonlinear elastic behavior and instabilities in crystals. Ultrasonics Symposium Proc. IEEE, pp. 999–1002.
Fumi, F. G. (1951). Thirdorder elastic coefficients of crystals. Phys. Rev. 83, 1274–1275.
Fumi, F. G. (1952). Thirdorder elastic coefficients in trigonal and hexagonal crystals. Phys. Rev. 86, 561.
Fumi, F. G. (1987). Tables for the thirdorder elastic tensors in crystals. Acta Cryst. A43, 587–588.
Green, R. E. (1973). Treatise on Material Science and Technology, Vol. 3. New York: Academic Press.
McSkimmin, H. J. (1964). Ultrasonic methods for measuring the mechanical properties of solids and fluids. Physical Acoustics, Vol. IA, edited by W. P. Mason, pp. 271–334. New York: Academic Press.
Melcher, R. L. & Scott, B. A. (1972). Soft acoustic modes at the cooperative Jahn–Teller transition in DyVO_{4}. Phys. Rev. Lett. 28, 607–610.
Michard, F., Zarembowitch, A., Vacher, R. & Boyer, L. (1971). Premier son et son zéro dans les nitrates de strontium, barium et plomb. Phonons, edited by M. A. Nusimovici, pp. 321–325. Paris: Flammarion.
Murnaghan, F. D. (1951). Finite Deformation in an Elastic Solid. New York: John Wiley and Sons.
Nouet, J., Zarembowitch, A., Pisarev, R. V., Ferré, J. & Lecomte, M. (1972). Determination of T_{N} for KNiF_{3} through elastic, magnetooptical and heat capacity measurements. Appl. Phys. Lett. 21, 161–162.
Rousseau, M., Gesland, J. Y., Julliard, J., Nouet, J., Zarembowitch, J. & Zarembowitch, A. (1975). Crystallographic, elastic and Raman scattering investigations of structural phase transitions in RbCdF_{3} and TlCdF_{3}. Phys. Rev. 12, 1579–1590.
Salje, E. K. H. (1990). Phase Transitions in Ferroelastic and Coelastic Crystals. Cambridge University Press.
Thurston, R. N. (1964). Wave propagation in fluids and normal solids. Physical Acoustics, Vol. 1A, edited by W. P. Mason, pp. 1–109. New York: Academic Press.
Truesdell, C. & Toupin, R. (1960). The classical field theories. Handbuch der Physik, Vol. III/1, edited by S. Flügge. Berlin, Göttingen, Heidelberg: SpringerVerlag.
Voigt, W. (1910). Lehrbuch der Kristallphysik. 2nd ed. (1929). Leipzig: Teubner. Photoreproduction (1966). New York: Johnson Reprint Corp.
Wallace, D. C. (1970). Thermoelastic theory of stressed crystals and higherorder elastic constants. Solid State Physics, Vol. 25. New York: Academic Press.
Wallace, D. C. (1972). Thermodynamics of Crystals. New York: John Wiley and Sons.
Zarembowitch, A. (1965). Etude théorique et détermination optique des constantes élastiques de monocristaux. Bull. Soc. Fr. Minéral. Cristallogr. 28, 17–49.