Notes on ANISOTROPY & CUBIC CRYSTALS

• The stress-strain relation is the three dimensional generalization of Hooke's law for an elastic solid relating the force (stress) to the displacement (strain). Usually the displacements are sufficient small that the linear relation s_ij = c_ijkle_kl is valid. However, the displacements considered here are large enough that the nonlinear terms contribute significantly and, in fact, give rise to shock formation.
  • For example, a recent experiments with waves generated by laser-generated, thermoelastic excitation in silicon generated particle velocities of over 40 m/s [A. Lomonosov and P. Hess, Nonlinear Acoustics in Perspective, edited by R. J. Wei (Nanjing University Press, Nanjing, 1996), pp. 106-111]. Because the surface acoustic wave speed for silicon in this direction is 4730 m/s, the Mach number equals 0.0085. To compare, the equivalent Mach number in air at room temperature is equivalent to a SPL of 155 dB (re 20 mPa).
  • The same experiment generated peak-to-peak particle displacements of over 700 nm in silicon. Because the lattice constant for silicon is only 0.543 nm, this means that the displacement is nearly 1300 times the size of the atoms in the crystal.
• Elastic constants specify the material properties of a particular crystal. The materials considered here are assumed to be fully elastic (i.e., no energy is lost during the deformation).
  • The energy of deformation has the form

    F = F0+k1 bij eij + k2 cijkl eijekl + k3 dijklmn eijeklemn +¼ .

    The elastic constants are named from this expression, i.e., the c_ijkl are the coefficient for the second order terms. Assuming the initial energy of the crystal is zero with zero deformation and applying the definition of the stress

    sij = ¶F ¶(¶ui/¶xj)  ,

    then gives rise to the stress-strain relation shown above.
  • In general, there are 81 possible SOE constants and 729 possible TOE constants. But because the c_ijkl and d_ijklmn are tensors and because the strain tensor is symmetrical, there is a maximum of 21 independent SOE constants and a maximum of 56 independent TOE constants (e.g., triclinic crystal).
• This presentation will only discuss cubic crystals, not because this is a limitation of the approach used here, but because this is the crystal type that has been most actively studied previous researchers. As a result, cubic crystals have the most empirical elastic constant data available for both the SOE and TOE constants.
• Because of the symmetries of cubic crystals, the number of independent SOE and TOE constants is reduced to 3 and 6, respectively. In contrast, isotropic materials have only 2 SOE and 3 TOE constants, respectively.
• One way to characterize cubic crystals is to construct the anisotropy ratio h. It is defined in such a way that isotropic materials have an anisotropy ratio of one. As will be shown later, the ratio allows the surface acoustic wave speed variation of the crystals to be conveniently grouped.
• Several specific examples of cubic crystals will be given in this presentation including potassium chloride (KCl), nickel (Ni), and silicon (Si). These materials were chosen because they exhibit a range of anisotropy ratios and crystal types.
• A feeling for the structure and symmetries of these crystals can be obtained from the simple structural models of the face-centered cubic (fcc) and diamond cubic crystals shown here.
  • Face centered cubic: The spheres each represent one atom in the crystal. For nickel, all the atoms are Ni atoms while for potassium chloride the atoms in the diagram alternate between K and Cl. Hence the latter structure can also be considered to be either a simple cubic lattice in which two different atoms alternate or two fcc lattices that inter-penetrate. The lattice spacing is 0.630 nm for KCl and 0.352 nm for Ni.
  • Diamond: Si has a diamond lattice which can also be considered to be two fcc lattices, one displaced relative to the other by (1/4, 1/4, 1/4). Note also that every atom has four nearest neighbors. The lattice spacing for Si is 0.543 nm.
  • While these crystals look different, they are considered cubic because of their symmetry properties. In a cubic crystal,

    The cube is invariant under rotations of 90^° about any of three perpendicular axes, and also invariant about the three planes perpendicular to these axes. It has three-fold rotation-inversion axes about lines between opposite corners, and two-fold rotation and mirror symmetries about the axis [110]. [from lecture notes for Condensed Matter Physics by Michael Marder]

    These are then the symmetries that are taken into account to determine how many independent elastic constants are needed to fully describe these systems.
• Because these systems are anisotropic, the wave propagation is different depending on how the crystal is cut and the direction that the wave is traveling.
  • The surfaces of cut crystals have traditionally been described using a crystallographic convention called Miller indices. Miller indices are defined by finding three non-collinear atoms on the surface that intersect the crystal axes and then applying the following method:

    1. Find the intercepts of the three basis axes in terms of the lattice constants.
    2. Take the reciprocals of these numbers and reduce to the smallest three integers having the same ratio. The result is enclosed in parentheses (hkl). [from Kittel, Introduction to Solid State Physics, 2nd ed., John Wiley & Sons, New York, p. 34 (1965).]

    Note that if the Miller indices are interpreted as a vector components, the resulting vector is normal to the surface of the cut.
  • Directions are specified in a different way:

    The indices of a direction in a crystal are expressed as the set of the smallest integers which have the same ratios as the components of a vector in the desired direction referred to the axis vectors. The integers are written in square brackets, [uvw]. The x axis is the [100] direction; the -y axis is the [0,-1,0] direction. A full set of equivalent directions is denoted this way: áuvwñ. [from Kittel, Introduction to Solid State Physics, 2nd ed., John Wiley & Sons, New York, p. 34 (1965).]

    This presentation will use both of these notations frequently.
• The diagrams show some typical crystal cuts and how they are specified using the Miller index notation. For example, the (001) intersects the crystal axes at, x = ¥ (1/¥ = 0), y = ¥, and z = 1. Similarly, the (111) plane intersects the crystal axes at x = 1, y = 1, and z = 1. Due to length constraints, this presentation will focus only on waves propagating in the (001) and (111) planes.

  Notes on THEORY

• Briefly, the approach used here involves the calculation of the Hamiltonian energy function through cubic order in the wave variables, choosing appropriate generalized coordinates, applying the equations of motion in canonical form, and deriving evolution equations for the slowly varying amplitudes in a suitable retarded time frame.
  • Note that computing the Hamiltonian the quadratic order would only give rise to linear terms in the model equations. Thus, the potential energy terms to at least cubic order in the strain must be included to model nonlinear effects.
  • Note also that this method is very general. It is applicable to any elastic material for which the SOE and TOE constants are known and to any cut and direction in such a material.
• Assumptions:
  1. It is assumed that the nonlinear solution is close to the linear solution; in particular the depth dependence of each frequency is the same as in the linear solution.
  2. It is assumed that the wave travels only in one direction, i.e., no compound waves.
• The velocity waveforms in the solid are written as Fourier series where the u_ni are the eigenfunctions of the linearized wave equation.
  • The coordinate system for the solution is chosen such that the the z-axis is perpendicular to the surface of the solid and the x-axis is in the direction of the propagation of the wave.
  • Note that on the surface the waveforms simplify to

    vi(x,z,t) = ¥ å n = -¥  vn(x) 3 å s = 1  b(s)i eint        [ v*n = v-n]

    where t = k₀x-w₀t is the retarded time and the b_i^(s) are determined from the linear problem.
  • Note that surface acoustic waves are non-dispersive, i.e., their wave speed is not frequency dependent.
• The coupled, nonlinear spectral evolution equations that result from this approach are shown above. Here v_n is the complex amplitude of the nth harmonic, a_n is the attenuation coefficient for the nth harmonic, w₀ is the characteristic angular frequency of the wave, r is the density of the material, c is the surface acoustic wave speed derived, and R_mn is the nonlinearity matrix.
  • The ad hoc attenuation term a_n = n²a₁ is added to the left-hand side for purposes of numerical stability when solving the equations. For all the cases shown here the dimensionless value of a₁ is 0.025. This attenuation is sufficiently weak that its main effect is to stabilize the portion of the waveform in the neighborhood of the shock without significantly the remainder of the waveform. Note that a₁ here is the analog of the the Goldberg number G for nonlinear acoustic waves in fluids.
  • Physically, the nonlinearity coefficients R_mn represent the strength of the coupling between different harmonics in the wave. They are given by a complicated analytical expression which can be determined completely by knowing the SOE and TOE constants of the material. (See the Supplement below for the equation for R_mn and more information.)
  • For the case of isotropic materials, these equations can be shown to reduce to the evolution equations previously derived by Zabolotskaya [E. A. Zabolotskaya, ``Nonlinear propagation of plane and circular waves in isotropic solids,'' J. Acoust. Soc. Am. 91, 2569-2575 (1992)].
• The solution procedure is outlined above. By first solving the linear problem for the eigenvalues, eigenvectors, and small-signal wave speed, the nonlinearity matrix can be constructed. Once the nonlinearity matrix is determined, the model equations can be integrated. A fourth order Runge-Kutta routine was used for the results shown here. Typically 200 harmonics were employed in each calculation.

  Notes on WAVE SPEED VARIATION

• One difference between the isotropic and anisotropic cases that exhibits itself even in the linear approximation is the variation of the wave speed as a function of the direction of the propagation of the wave.
• Here the surface acoustic wave speed relative to the bulk shear speed is plotted versus direction for several different materials. The top graph shows the wave speed for the (001) cut as a function of the angle from the á100ñ direction (periodic every 90^°) while the bottom graph shows the wave speed for the (111) cut as a function of the angle from the á1,1,-2ñ (periodic every 60^°). For comparison, the wave speed of fused quartz is also plotted. (fused quartz is isotropic and hence has a constant wave speed in all directions).
• Note that the wave speed curves group by anisotropy ratio. Anisotropy ratios greater than one give slower wave speeds while anisotropy ratios less than one give faster wave speeds. See the review paper by Farnell for similar graphs with more materials [G. W. Farnell, ``Properties of elastic surface waves,'' Physical Acoustics, Vol. 6, ed. W. P. Mason, Academic Press, New York, 1970, pp. 109-166].

• This diagram shows the particle motion of a typical surface acoustic wave in an anisotropic medium. Consider the case of a surface acoustic wave with an initially sinusoidal velocity waveform in an isotropic and anisotropic material. Assume that the x-axis is in the propagation direction and that the z-axis is normal to the surface cut.
• In the isotropic case, the particle motion will be elliptical, and the plane of the ellipse will be the xz-plane (indicated by the dashed rectangle in the side view above). There is no transverse component to the particle motion.
• In the anisotropic case, the particle motion is still elliptical, but most generally the plane of the ellipse is rotated by some angle (indicated here by the letter f) out of the xz-plane. Different views of this transverse motion are shown in the top and front views.
• To summarize, surface acoustic wave motion has generally has longitudinal (x-axis), vertical (z-axis), and transverse (y-axis) components.

  Notes on WAVEFORM DISTORTION: CRYSTALS

• These diagrams show the nonlinear waveform distortion and shock formation for KCl, Ni, and Si in the (111) plane in the á1,1,-2ñ direction.
• The left column contains the longitudinal velocity waveform while the right column contains the vertical velocity waveform. Due to the symmetry of this particular case, there is no transverse velocity component.
• In each waveform, a non-dimensional velocity component is plotted versus a non-dimensional retarded time. The initial waveform is a single frequency, continuous wave. The velocity v is scaled such that the total initial velocity amplitude is unity, and the time t is scaled by the period of the initial sinusoidal signal. The propagation distance X is scaled such that X = 1 corresponds to approximately one shock formation distance.
• KCl: In its longitudinal waveform, cusps form near the shock front while in its vertical waveform. a peak forms. Note that the wave distorts asymmetrically in both cases. Interestingly, the longitudinal waveform looks similar to the waveform distortion of an acoustic wave in a fluid except that here the coefficient of nonlinearity here is negative. Hence peaks move slower and troughs, faster in contrast to the behavior in fluids where peaks move faster, and troughs, slower.
• Ni: In contrast to KCl, here the vertical waveform forms a shock and longitudinal waveform forms a peak. However, in this case the coefficient of nonlinearity is positive so that the peaks move faster and troughs, slower.
• Si: Here the waveforms look qualitatively similar to Ni but there also seem to be some mixed features. For example, the vertical waveform forms a shock but the cusping is so asymmetrical that one looks more like a peak. The vertical waveform forms a distinct peak but the peak is so asymmetrical that one side looks like a shock.
• As can be seen, the distortion of the waveforms is very complicated. The cause of the specific features of the distortion in the various crystals has not yet been explained and will be a source of future work.

  WAVEFORM DISTORTION: DIRECTIONS

• Next consider what happens when the crystal, cut, and propagation distance are held constant and only the direction of propagation is changed.
• Here the longitudinal, vertical, and transverse velocity waveforms are plotted for KCl in the (111) plane for a variety of directions relative to the á1,1,-2ñ direction. Each individual waveform shown is the waveform that results after an initially sinusoidal signal is propagated to a distance of X = 2 for that case. The scaling is the same as in the previous set of graphs.
• Here the mixing of features between the various components is seen even more than in the previous plots. The longitudinal waveform at 0^° starts in a shock, and turns into more of a peaked waveform as the angle is increased to 5^°. Over the same change in direction, the peak in the vertical waveform turns into an asymmetric peak with one side steepened. Finally, the transverse component (which is zero at 0^°) appears as a peak and grows in amplitude as the angle increases.
• The cause of these changes as a function of direction will also be subject of future work.

  Notes on ATYPICAL ENERGY TRAPPING

• A particularly interesting case of atypically energy trapping in the lowest order harmonics arises for wave propagation in KCl in the (001) plane in the á100ñ direction. While this type of distortion has been seen before in dispersive media (e.g., nonlinear optical media for light, bubbly liquids for sound), we believe that this is the first time that it has been predicted for a non-dispersive medium.
• First look at the top pair of plots. The left plot is the longitudinal velocity waveform and the right plot is the vertical velocity waveform. Due to the symmetry of this case, no transverse velocity component exists. The scaling is the same as in the previous waveform plots. Note that no shock formation occurs- only low frequency ``wiggles'' appear.
• This is not a numerical effect but a physical effect caused by the fact that the nonlinearity matrix element that describes the strength of the coupling between the first harmonic and higher order harmonics is approximately an order of magnitude smaller than neighboring matrix elements. As a result, it is more difficult for energy to transfer from lower to higher harmonics as required for shock formation.
• This atypical energy trapping in the lowest harmonics can be seen by comparing the harmonic propagation curves for a more typical case with the case described above. The plots in the ``Comparison'' section show the amplitudes of the first five harmonics as a function of propagation distance. The harmonic components are scaled relative to the initial amplitude of the first harmonic and the propagation distance is scaled as described previously.
• The top plot shows the harmonic propagation curves for KCl in the (111) plane in the á1,1,-2ñ direction (the waveforms shown on the previous pages). The first harmonic amplitude falls while the other harmonic amplitudes rise as nonlinear interactions cause energy to flow into the higher frequency terms. At each propagation distance, the amplitudes decrease monotonically from the first to the fifth harmonic.
• In contrast, the harmonic propagation curves for KCl in the (001) plane in the á100ñ direction are very different. While the first harmonic falls, the the second harmonic rises much more quickly and even exceeds the amplitude of the first harmonic around X = 2. The third and fourth harmonics rise much more slowly initially than in the previous case but eventually grow to amplitudes comparable to the first harmonic before attenuation causes all the harmonic components to decrease.
• As can be seen from this harmonic propagation curve, most of the energy is contained within the first two harmonics at least close to the source. If the system of model equations is truncated after just the first two terms and no dissipation is assumed, then by applying the boundary conditions shown above allows an analytical solution to be derived.
• A graphical comparison of the analytical and numerical solution shown above. As expected, the agreement is good close to the source X < 1 for both harmonics and reasonably good for the first harmonic out to about X = 4. Part of the discrepancy between the two solutions comes from the fact that the numerical solution contains dissipation that the analytical solution does not.
  • If dissipation terms are added to the truncated system of evolution equations and then that system is integrated numerically, it can be shown that agreement is virtually identical near the source. The dissipation also accounts for about half of the difference between the curves at the further propagation distances.

  Notes on CONCLUSION & FUTURE WORK

• To conclude, we investigated a set of real cubic crystal waveforms and found that they exhibit asymmetrical distortion, shocks with cusped spikes, and a phase shift of the waveforms as a function of direction. We also predicted atypical energy trapping in the lowest order harmonics in a non-dispersive medium for the case of KCl for the (001) plane in the á100ñ and derived an approximate analytical solution valid close to the source.
• Future work will probably include studying the relationship between the nonlinearity matrix and waveform distortion, deriving an approximate analytical expression for the shock formation distance, comparing pulsed waveform modeling to experimental data, and investigating non-cubic crystals.

  Supplement: NONLINEARITY MATRIX & LINEAR THEORY

The nonlinearity matrix is given by

Start with linearized wave equation