Resonant ultrasound spectroscopy is used to nondestructively measure the elastic resonances of small solids to elucidate the material's elastic properties or other qualities like size, shape, or composition. Here, we introduce the software *RUScal* for the purpose of determining elastic properties by analyzing the eigenfrequencies of solid specimens with common shapes, such as rectangular parallelepipeds, cylinders (solid and hollow tube), ellipsoids, and octahedrons, as well as irregularly shaped ellipsoids that can be described analytically. All symmetry classes are supported, from isotropic to triclinic, along with the option to add or remove up to three orthogonal mirror planes as well as the ability to reorient the crystal axes with respect the sample edges via Euler angles. Additional features include tools to help find initial sets of elastic constants, including grid exploration and Monte Carlo methods, a tool to analyze frequencies as a function of sample length or crystal orientation, an error analysis tool to assess fit quality, and formatting of the input and output files for batch fitting, e.g., as a function of temperature. This software was validated with published resonant ultrasound spectroscopy data for various materials, shapes, and symmetries with noted improvements in calculation time compared to finite element methods.

## I. INTRODUCTION

Mechanical resonances of solids, like the ringing of a bell, carry information regarding size, shape, and composition, including elastic properties of the resonator. Resonant ultrasound spectroscopy (RUS) is a nondestructive and efficient technique used to measure these mechanical resonances, typically for small, mm-to-cm-sized, solids such as metals, ceramics, and semiconductors. In principle, all adiabatic elastic constants can be obtained from a single RUS measurement. A typical experimental setup is shown in Fig. 1 and consists of a sample whose shape, ideally, can be described analytically (e.g., sphere, cylinder, or rectangular parallelepiped) mounted between piezoelectric transducers: one transducer mechanically excites the sample, and the other picks up the sample's response, which is processed using a lock-in amplifier. When the frequency of the driving transducer matches the eigenfrequency of a full-body mechanical resonance in the sample, the result is a sharp peak in the ultrasound spectrum. Determination of the elastic constants proceeds through modeling in two parts: solving the forward and inverse problems. In the forward problem, the ultrasound spectrum of the material is calculated based on the sample's shape, density, and approximate or estimated elastic constants. The goal of the inverse problem is then to iteratively adjust the elastic constants or other adjustable parameters, such as shape or crystal orientation, until the difference between the calculated and measured resonant frequencies is minimized under certain figures of merit, e.g., the root-mean-squared (rms) error of the residuals. Migliori and Maynard (2005) and Balakirev *et al.* (2019) provide guidance in the assembly and operation of a RUS apparatus.

The RUS technique is relatively new among the nondestructive acoustical methods and traceable to the post-WWII era of geology, seismology, and earth sciences in the 1950s and 1960s. LeCraw (1964) led what amounts to the first experimental RUS investigations for materials science using the resonant sphere technique (RST) for the special case of an isotropic sphere whose free vibrations can be solved analytically (Lamb, 1881). For other cases, the Rayleigh–Ritz method was used to approximate eigenvalues, such as in the work done for isotropic cubes (Holland, 1968). Soga and Anderson (1967) then developed the two-transducer, pump-probe-type RUS technique to exploit the large amplitude surface vibrations of a resonating solid to measure the elastic moduli of anisotropic single crystals up to 1800 K. In the late 1960s and early 1970s, Anderson's then graduate student Harold Demarest developed the theory for cube-shaped samples with cubic symmetry by using orthogonal (Legendre) polynomials rather than trigonometric functions to express the displacement field (Demarest, 1971); this work was known as the cube-resonance method. Then Ohno further developed Demarest's theory for single crystals by extending calculations to crystals with tetragonal, orthorhombic, and trigonal symmetries through the rectangular parallelepiped resonance (RPR) method (Ohno, 1976; Ohno *et al.*, 1986). In the 1990s, Migliori and Sarrao (1997) recognized the critical discovery by Visscher *et al.* (1991) that accommodates arbitrary shapes, albeit at the expense of computing time, by using simple powers of Cartesian coordinates as a choice of basis. Migliori and his colleagues then continued to make strides in automating signal processing and optimization of the measurement system (Migliori, 1991). The increasing power of the modern personal computer has since accelerated the numerical forward and inverse calculations. In the few decades since its inception, RUS has been used widely to study the physical properties of various materials, to probe first- and second-order phase transitions (e.g., Shekhter *et al.*, 2013; Migliori *et al.*, 1990), and to assess specimen quality (Schwarz and Vuorinen, 2000).

Software to accomplish forward and inverse calculations exists freely in the public domain (Balakirev *et al.*, 2019; Fig, 2006; Zadler *et al.*, 2014; Ramshaw, 2015), commercially (e.g., Dynamic Resonance Systems, Magnaflux Quasar Systems, Insight K.K., and Alamo Creek Engineering), or in some cases upon request from research groups. In general, these programs support orthorhombic or higher symmetry and common shapes like cylinders, spheres, and cuboids. Some packages include advanced analysis features, such as genetic algorithms (Remillieux *et al.*, 2015; Ramshaw, 2015), machine learning (Ghosh *et al.*, 2020; Grueso, 2021), finite element analysis for complex shapes or defect detection (Wang *et al.*, 2019; Flynn and Radovic, 2011), or combinations thereof. Perhaps the most widely known programs for routine analysis are *RPR*.exe and *Cyl*.exe by Migliori *et al.* (1993).

In this work, we outline new software for the analysis of resonant ultrasound spectra named *RUScal* written in C language that combines and expands the functionality of existing software with additional shapes, symmetries, and analysis tools for routine RUS analysis. The major additions include hollow cylinder/tube and “potato” shapes (Visscher *et al.*, 1991) as well as all trigonal, tetragonal, monoclinic, and triclinic crystal symmetries, including chiral systems, with the option to turn on/off up to three orthogonal mirror planes. In addition, three analysis tools are included: (1) Monte Carlo search programs with optional genetic algorithm follow-up, to assist in searching for fit input parameters, i.e., elastic constants or Euler angles; (2) a shape evaluation program that calculates ultrasound spectra as a function of sample dimensions or crystalline orientation; and (3) an error analysis program that helps in determining the quality of a fit. The layout of the paper is as follows. First, an introduction to the theory of elasticity is given in Sec. II to familiarize the user with the mathematical methods and physical concepts used in the code. Then an overview of the analysis software is provided, which includes instructions for performing measurements and using the software with a description of its functionalities. In Sec. III, we benchmark many new features of the software against exemplary materials from the literature. The *RUScal* software is available upon request or at https://github.com/RUS-ORNL/RUScal/.

## II. METHODS

### A. The forward problem

The calculation of eigenmodes of a vibrating solid with free boundary conditions is a classical mechanics problem, the formalisms of which can be found in many publications and textbooks [see, e.g., Migliori and Sarrao (1997) and Ashcroft and Mermin (1976)] and are copied here for clarity. In the RUS technique, a sample is placed with point-like contact with two (or more) transducers, one of which mechanically excites the sample with a sinusoidal applied force that is small such that stress and strain are linearly proportional, satisfying Hooke's law,

Here, $\sigma ij$ and $\epsilon kl$ are components of the second-rank stress and strain tensors, respectively, and $Cijkl$ are the 81 components of the fourth-rank elastic stiffness tensor. The Einstein convention for summation notation is used here. The number of elastic constants can be reduced to 21 in the absence of net body torques, regardless of the crystalline symmetry of the material (Leisure and Willis, 1997).

Exact solutions of free oscillators of arbitrary shape and symmetry are not yet realized except for special cases, such as the isotropic sphere (LeCraw, 1964). However, numerical approximations are available. The finite element method (FEM) enables the calculation of eigenmodes of samples with complex shapes and boundary conditions. By comparison, energy minimization (variational) methods are computationally more efficient albeit limited to a smaller number of geometric shapes. Here, we will use the latter, in particular the Rayleigh–Ritz method, to determine the resonance frequencies of a vibrating solid for several conventional shapes, detailed in Sec. II A 2. The following theoretical calculations largely follow those of Migliori and Sarrao (1997). The resonance frequencies of a vibrating three-dimensional (3D) solid are the stationary points of the Lagrangian $L$,

for which the solid has volume *V* that is bound by free surface $S$. The quantities $T$ and $U$ are the kinetic and potential energies per unit volume, respectively, described in terms of a displacement field $u(r)$ (herein shortened to $u$) with components $ui$ (*i*, *j, k, l* = 1, 2, 3 corresponding to *x*, *y*, *z*) (Migliori and Sarrao, 1997),

It is assumed that the time dependence of the displacements is harmonic of the form $ei\omega t$ with angular frequency $\omega $. Using the Rayleigh–Ritz method, each component of the displacement vector is expanded in a complete set of basis functions $\Phi \lambda $,

where $ai,\lambda \u2009(ak,\gamma )$ are the expansion coefficients, and $\lambda (\gamma )=l,m,n$. Each set of basis functions is chosen based on the shape of the sample, which will be discussed further in Sec. II A 2. Following Visscher *et al.* (1991), we expand the function $\Phi \lambda $ in the powers of Cartesian coordinates,

where *l*, *m*, and *n* are non-negative integers. We note that this choice in basis functions, however, fails to describe the discontinuity in strain across an interface between two materials with different densities and elastic constants, such as with coatings or composites. Substituting Eq. (4) into Eqs. (3a) and (3b), the kinetic and potential energy terms become

The Lagrangian can be rewritten in a more compact way,

where the matrices $E$ and $\Gamma $ contain the following elements:

For practical purposes, the series expansion is truncated under the condition

with the number of basis elements $R=3N+1N+2(N+3)/6$. In principle, the limit *N* → ∞ results in the exact solution; however, a polynomial order between 8 and 14 is sufficient to achieve approximate solutions with manageable matrix sizes and computation time to reach convergence (more discussion later in Fig. 3). Generally, the execution time scales as *N*^{3} with some variation when adding mirror planes. In the absence of external forces (Migliori and Sarrao, 1997), stationary points of the Lagrangian are obtained by setting the derivatives of Eq. (7) to zero, resulting in the eigenvalue equation,

which is a generalized symmetric eigenvalue problem with eigenvalues $\rho \omega 2$ and eigenvectors $a$. Equation (10) enables the determination of the resonant frequencies from known elastic constants.

#### 1. Crystal orientation

To achieve the best possible results, crystalline samples should be prepared such that their crystallographic axes align with the edges of the sample body. However, in practice, the crystal can become miscut, or the orientation may be unknown, either of which results in a misaligned elastic tensor. A noteworthy solution may be to polish the sample into a sphere where the crystal orientation is not needed, or rather is ill-defined, with respect to the sample body. In the absence of recutting or polishing, the stiffness tensor needs to be rotated with respect to the sample body using either direct knowledge of the degree of misorientation or fitting of the input parameters. Here, we must use the full tensor notation for the elastic constants $Cijkl$ and apply the rotation in the following way (Arfken *et al.*, 2012; Gladden, 2003):

where $C\u2032\alpha \beta \gamma \delta $ are the elastic constants in the rotated reference frame, and the rotation matrix $R$ is expressed in Kock's convention (Arfken *et al.*, 2012),

The Euler angles psi, theta, and phi refer to rotations (in degrees) about the *z* axis, about the new *y*′ axis, and ($180\xb0\u2212\varphi $) about the new z′′ axis as illustrated in Fig. 2. With these definitions, ($\psi ,0\xb0,\psi )$ is equivalent to $0\xb0,0\xb0,0\xb0$. Gladden (2003) provides an example of utilizing the orientation of a single crystal of alumina (Al_{2}O_{3}) with trigonal symmetry that was intentionally miscut to determine the sign of an elastic constant and to demonstrate the application of Eq. (11); we reproduced Gladden's results in Sec. III D using the *RUScal* code.

#### 2. Shapes

The volume integrals appearing in Eqs. (8a) and (8b) from Sec. II A may be modified to accommodate many shapes by choosing the appropriate set of basis functions ${\Phi \lambda}$. Based on the work by Visscher *et al.* (1991), these matrix elements take the form

where *p*, *q*, and *r* are non-negative integers. Below are the integral forms of Eq. (13) for each shape used in the *RUScal* code under the function *integ*(): the rectangular parallelepiped,

where $(x,y,z)$ correspond to the principal axes $a,b,c$; the sphere and ellipsoid,

the cylinder,

the hollow cylinder,

where $d0$ is the inner diameter, and $A$ is the ratio of the outer-to-inner diameter; the octahedron,

and the potato,

where, as in Visscher *et al.* (1991), the factor in curly brackets is unity if two or three of the integers $p,\u2009q,\u2009r$ are odd and $\pi /2$ otherwise. ξ, υ, and ζ are positive or negative depending on the respective sign of *x*, *y*, and *z* in the considered octant.

### B. The inverse problem

The forward problem enables calculation of resonant frequencies given some physical properties of a sample, such as its shape, density, and elastic constants. However, it is often the case that the interest lies in determining physical properties of the sample from the resonant frequencies, which is the inverse problem. The elastic constants are the usual target of the inverse problem and are determined by iteratively adjusting the initial set of elastic constants until the difference between the set of measured and calculated resonance frequencies is minimized. The figure of merit of the closeness of fit is assessed by chi-squared error—the sum of weighted residuals from $n$ eigenfrequencies (Leisure and Willis, 1997),

where $wn$ is the weight factor between zero and unity for a given level of confidence in a measured resonance, and $fn$ and $gn$ are the measured and calculated resonant frequencies, respectively. Generally, when the correct elastic constants are used, as the polynomial order *N* used in the calculation increases, the difference between calculated and measured eigenfrequencies decreases as illustrated in Fig. 3 for three modes of commercial 4140 steel: lower frequency modes, e.g., first <20 modes, rapidly approach their limiting values with increasing *N*. Because a lower *N* corresponds to faster computation, *N* = 8 or 10 is useful during initial fitting of a few modes. In fitting with more than ∼30 modes, a higher polynomial order is needed. In practice, we find a polynomial order of 12 or 14 is a good compromise between computation time and precision. In general, the number of modes to use should be at least five times the number of free parameters.

#### 1. The measurement

The accuracy of a fit depends on the quality of the measurement. What follows in the remainder of this section is a procedure for RUS measurements adapted from Balakirev *et al.* (2019) to guide new users. The key to a successful RUS experiment is preparation of a good sample. This includes not only accurate measures of the dimensions (±10 μm), mass (±1 mg), and crystallographic alignment with the sample's edges (±1°) but also care in avoiding cracks, voids, or other inhomogeneities that adversely affect the sample's shape, symmetry, or elastic moduli. The acoustical quality of the sample is assessed in the ultrasound spectrum by the quality factor (*Q*), which is a measure of ultrasonic attenuation and calculated as the ratio of the resonant peak's center frequency $f0$ to its width (full width at half-maximum $w$): $Q=f0/w$. High-*Q* materials such as metals typically exhibit *Q* > 10^{3} and appear as sharp peaks in their ultrasound spectra that are well-approximated with Lorentzian line shapes (see Fig. 4). Low-*Q*, broad resonances (*Q* < 100) lead to uncertainty in determining exact frequencies and therefore promote greater uncertainty in the inverse calculation.

Prior to the experiment, it is advised to calculate the ultrasound spectrum, e.g., about 50–100 resonances, using the measured mass, dimensions, and approximate elastic constants. The calculated spectrum serves as a reference for the experiment. Details on how to execute this calculation using *RUScal* are given in Sec. II C 1. A common RUS setup that uses two direct-contact transducers is depicted in Fig. 1. In this case, one carefully places the sample between transducers to avoid (1) damaging the edges or surface and (2) compressing the sample by more than a couple grams-weight of force. The latter is important because the solution to the equation of motion assumes a free-body resonator. Strictly speaking, the exact placement between the transducers is not important with the understanding that some resonances may have surface displacements that are parallel to the transducers' surfaces; the result is a missing mode that will become obvious when comparing against the calculated spectrum and during fitting—see discussion in Secs. 5A and 5B in Balakirev *et al.* (2019). In anticipation of missing a mode, we advise remounting the sample in multiple different orientations and rescanning the frequency range.

While the drive transducer excites the sample, the pickup or receiving transducer detects the sample's response. The receiving signal is a complex voltage response containing both in-phase (real) and quadrature (imaginary) parts—see Fig. 4(b). For a high-*Q* material, these parts can be fit to a complex Lorentzian of the form (Balakirev *et al.*, 2019)

where $A$, $\theta $, $\gamma $, and $Zbg$ are the peak's amplitude, phase, half-width, and background, respectively, and $f$ is frequency. Access to both parts of the signal, along with the calculated spectrum, helps in discarding non-sample artifacts in the ultrasound spectrum, such as from crosstalk between transducers or resonances from buffer rods or other hardware.

### C. The input and output files

The user primarily interacts with the *RUScal* C-program by means of a text (.dat) file with the default name *rusin*, which has a similar format as the text files *rprin* and *cylin* used in the *RPR*.exe and *Cyl*.exe programs, respectively. In viewing these programs to be an industry standard, we formatted our input and output data to be familiar to those users. Figure 5 shows an example of a *rusin* file for a cylindrical sample of 4140 steel whose ultrasound spectrum was measured using the RUS setup described by Balakirev *et al.* (2019). Table I provides descriptions of each parameter used in the *rusin* file along with many options for calculating and fitting data. The pound symbol (#) can be used anywhere in the input file, e.g., for commenting or recordkeeping purposes. The colon symbol (:) can be used in the title line (first line of *rusin* file) to establish a tracking parameter as shown in Fig. 5. For the tracking feature to work properly, place only a float number after the colon, as the software attempts to convert a string (text) to a float (number). If no number is given, the default value of zero is used.

Line . | Column . | Parameter . | Description . |
---|---|---|---|

1 | Header | User comments, e.g., sample name, 128-character maximum. A colon delimiter can be used as a tracking parameter (see text). | |

2 | Comments | Headers for the sample properties in line 3. | |

3 | 1 | Shape | Cylinder = 0, RPR = 1, potato = 4, sphere = 5, octahedron = 6, hollow cylinder = 8 |

3 | 2 | Bravais class | Number of elastic constants (#-cxx) based on Bravais class. |

Isotropic = 2, cubic = 3, hexagonal = 5, trigonal = 306 (32, −3 m, 3 m classes), trigonal = 307 (3, −3 classes), tetragonal = 406 (4 mm, −42 m, 422, 4/mmm classes), tetragonal = 407 (4, −4, 4/m classes), orthorhombic = 9, monoclinic = 13 (all classes, note the diad is parallel to y dimension), triclinic = 21 | |||

3 | 3 | Polynomial | Polynomial order N [see Eq. (9)], (n-poly). Max = 20. |

3 | 4 | Calculation mode #-calc | If 0, then fit free parameters. |

If any integer n > 0, then calculate n modes. | |||

If −1, then execute Monte Carlo search. | |||

If −2, then execute shape evaluation. | |||

If −3, then execute error analysis. | |||

If −4, then execute chi-squared grid. | |||

If −5, then execute genetic algorithm (GenAlg). | |||

3 | 5 | Mirror planes | 0 = no mirror planes, 1 = xy mirror plane, 2 = xy and yz, 3 = xy, yz, and xz |

3 | 6 | Sample mass | Units of grams |

3 | 7 | Chi squared | Chi-squared error from the previous fit. Default is 1.00. |

3 | 8 | Derivatives | Calculate derivatives. 0 = do not calculate, 1 = calculate. |

3 | 9 | Population | GenAlg option: initial population. Default = 500 (integer) |

3 | 10 | Generations | GenAlg option: number of generations. Default = 10 (integer) |

3 | 11 | Evolution | GenAlg option: evolutionary pressure. Default = 3.0 (float) |

3 | 12 | Mixing | GenAlg option: degree of parent mixing. Default = 0.2 (float) |

3 | 13 | Mutation | GenAlg option: probability of mutation. Default = 1.0 (float) |

4 | Comments | Headers for the elastic constants in line 5 | |

5 | Elastic constants | Order of elastic constants for each Bravais class. Units of Mbar or 100 GPa. | |

Isotropic (2) $c11$ and $c44$ | |||

Cubic (3) $c11$, $c12,$ $c44$ | |||

Hexagonal (5) $c33,\u2009c23,\u2009c12,$ $c44$, $c66$ | |||

Trigonal (306) $c11,\u2009c33,\u2009c13,\u2009c12,\u2009c44,\u2009c14$ | |||

Trigonal (307) $c11,\u2009c33,\u2009c13,\u2009c12,\u2009c44,\u2009c14$, $c25$ | |||

Tetragonal (406) $c11,\u2009c33,\u2009c23,\u2009c12,\u2009c44,$ $c66$ | |||

Tetragonal (407) $c11,\u2009c33,\u2009c23,\u2009c12,\u2009c44,$ $c66$, $c16$ | |||

Orthorhombic (9) $c11,\u2009c22,\u2009c33,\u2009c23,\u2009c13,\u2009c12,\u2009c44,\u2009c55,$ $c66$ | |||

Monoclinic (13) c_{11}, c_{22}, c_{33}, c_{23}, c_{13}, c_{12}, c_{44}, c_{55}, c_{66}, c_{15}, c_{25}, c_{35}, c_{46} | |||

Triclinic (21) c_{11}, c_{22}, c_{33}, c_{23}, c_{13}, c_{12}, c_{44}, c_{55}, c_{66}, c_{14}, c_{15}, c_{16}, c_{24}, c_{25}, c_{26}, c_{34}, c_{35}, c_{36}, c_{45}, c_{46}, c_{56} | |||

6 | Constraints | Constraints on the elastic constants | |

If calc mode = 0, then 0 = fit and 1 = fix this constant. | |||

If calc mode = −1, then value equals max deviation of c_{xx} in per thousands or max deviation of Euler angle in degrees. | |||

If calc mode = −3, then values are not used for calculation. | |||

If calc mode = −4, then value equals max deviation of c_{xx} in per thousands. | |||

7 | Comments | Headers for the dimensions in line 8 | |

8 | 1 | Dimension 1 | Units of cm. Cylinder = major axis, RPR = side 1 or a axis, ellipsoid = a axis, hollow cylinder = outer diameter, potato = x^{+}, octahedron = x^{+} |

8 | 2 | Dimension 2 | Units of cm. Cylinder = minor axis, RPR = side 2 or b axis, ellipsoid = b axis, hollow cylinder = inner diameter, potato = y^{+}, octahedron = y^{+} |

8 | 3 | Dimension 3 | Units of cm. Cylinder = height, RPR = side 3 or c axis, ellipsoid = major axis, hollow cylinder = outer height, potato = z^{+}, octahedron = z^{+}. |

8 | 4 | Dimension 4 | Units of cm. Potato = x^{−} |

8 | 5 | Dimension 5 | Units of cm. Potato = y^{−} |

8 | 6 | Dimension 6 | Units of cm. Potato = z^{−}. Note: x^{+} = x^{−} = y^{+} = y^{−} = z^{+}, z^{−} = 0 for a hemisphere. |

9 | Constraints | Constraints on the dimensions | |

If calc mode = 0, then 0 = fit and 1 = fix. Volume remains constant. | |||

If calc mode = −2, then vary one of the following dimensions: 1 = x^{+}, 2 = y^{+}, 3 = z^{+}, 4 = x^{−}, 5 = y^{−}, 6 = z^{−}, 7 = $\psi $, 8 = $\theta $, 9 = $\varphi $. | |||

10 | Comments | Headers for the Euler angles in line 11 | |

11 | 1 | Euler angle 1 | Angle $\psi $ for rotations about the z axis in Kock's convention |

11 | 2 | Euler angle 2 | Angle $\theta $ for rotations about the new y axis in Kock's convention |

11 | 3 | Euler angle 3 | Angle $\varphi $ for rotations of ($\pi $-$\varphi $) about the new z axis in Kock's convention |

12 | Constraints | Constraints on the Euler angles. 0 = fit, 1 = fix | |

If calculation mode = −1, then value equals max deviation of Euler angle in degrees. | |||

If calculation mode = −2, then value corresponds to dimensions or angle to vary: dimensions (1–6) or Euler angles (7 = phi, 8 = theta, psi = 9). | |||

13 | Comments | Headers for the list of frequencies in line 14 | |

14+ | 1 | Frequencies | Measured frequencies in units of MHz |

14+ | 2 | Frequencies | Calculated frequencies in units of MHz. Calculated by the forward problem. |

14+ | 3 | Weights | Weight factors for each measured frequency. See Eq. (20). |

Line . | Column . | Parameter . | Description . |
---|---|---|---|

1 | Header | User comments, e.g., sample name, 128-character maximum. A colon delimiter can be used as a tracking parameter (see text). | |

2 | Comments | Headers for the sample properties in line 3. | |

3 | 1 | Shape | Cylinder = 0, RPR = 1, potato = 4, sphere = 5, octahedron = 6, hollow cylinder = 8 |

3 | 2 | Bravais class | Number of elastic constants (#-cxx) based on Bravais class. |

Isotropic = 2, cubic = 3, hexagonal = 5, trigonal = 306 (32, −3 m, 3 m classes), trigonal = 307 (3, −3 classes), tetragonal = 406 (4 mm, −42 m, 422, 4/mmm classes), tetragonal = 407 (4, −4, 4/m classes), orthorhombic = 9, monoclinic = 13 (all classes, note the diad is parallel to y dimension), triclinic = 21 | |||

3 | 3 | Polynomial | Polynomial order N [see Eq. (9)], (n-poly). Max = 20. |

3 | 4 | Calculation mode #-calc | If 0, then fit free parameters. |

If any integer n > 0, then calculate n modes. | |||

If −1, then execute Monte Carlo search. | |||

If −2, then execute shape evaluation. | |||

If −3, then execute error analysis. | |||

If −4, then execute chi-squared grid. | |||

If −5, then execute genetic algorithm (GenAlg). | |||

3 | 5 | Mirror planes | 0 = no mirror planes, 1 = xy mirror plane, 2 = xy and yz, 3 = xy, yz, and xz |

3 | 6 | Sample mass | Units of grams |

3 | 7 | Chi squared | Chi-squared error from the previous fit. Default is 1.00. |

3 | 8 | Derivatives | Calculate derivatives. 0 = do not calculate, 1 = calculate. |

3 | 9 | Population | GenAlg option: initial population. Default = 500 (integer) |

3 | 10 | Generations | GenAlg option: number of generations. Default = 10 (integer) |

3 | 11 | Evolution | GenAlg option: evolutionary pressure. Default = 3.0 (float) |

3 | 12 | Mixing | GenAlg option: degree of parent mixing. Default = 0.2 (float) |

3 | 13 | Mutation | GenAlg option: probability of mutation. Default = 1.0 (float) |

4 | Comments | Headers for the elastic constants in line 5 | |

5 | Elastic constants | Order of elastic constants for each Bravais class. Units of Mbar or 100 GPa. | |

Isotropic (2) $c11$ and $c44$ | |||

Cubic (3) $c11$, $c12,$ $c44$ | |||

Hexagonal (5) $c33,\u2009c23,\u2009c12,$ $c44$, $c66$ | |||

Trigonal (306) $c11,\u2009c33,\u2009c13,\u2009c12,\u2009c44,\u2009c14$ | |||

Trigonal (307) $c11,\u2009c33,\u2009c13,\u2009c12,\u2009c44,\u2009c14$, $c25$ | |||

Tetragonal (406) $c11,\u2009c33,\u2009c23,\u2009c12,\u2009c44,$ $c66$ | |||

Tetragonal (407) $c11,\u2009c33,\u2009c23,\u2009c12,\u2009c44,$ $c66$, $c16$ | |||

Orthorhombic (9) $c11,\u2009c22,\u2009c33,\u2009c23,\u2009c13,\u2009c12,\u2009c44,\u2009c55,$ $c66$ | |||

Monoclinic (13) c_{11}, c_{22}, c_{33}, c_{23}, c_{13}, c_{12}, c_{44}, c_{55}, c_{66}, c_{15}, c_{25}, c_{35}, c_{46} | |||

Triclinic (21) c_{11}, c_{22}, c_{33}, c_{23}, c_{13}, c_{12}, c_{44}, c_{55}, c_{66}, c_{14}, c_{15}, c_{16}, c_{24}, c_{25}, c_{26}, c_{34}, c_{35}, c_{36}, c_{45}, c_{46}, c_{56} | |||

6 | Constraints | Constraints on the elastic constants | |

If calc mode = 0, then 0 = fit and 1 = fix this constant. | |||

If calc mode = −1, then value equals max deviation of c_{xx} in per thousands or max deviation of Euler angle in degrees. | |||

If calc mode = −3, then values are not used for calculation. | |||

If calc mode = −4, then value equals max deviation of c_{xx} in per thousands. | |||

7 | Comments | Headers for the dimensions in line 8 | |

8 | 1 | Dimension 1 | Units of cm. Cylinder = major axis, RPR = side 1 or a axis, ellipsoid = a axis, hollow cylinder = outer diameter, potato = x^{+}, octahedron = x^{+} |

8 | 2 | Dimension 2 | Units of cm. Cylinder = minor axis, RPR = side 2 or b axis, ellipsoid = b axis, hollow cylinder = inner diameter, potato = y^{+}, octahedron = y^{+} |

8 | 3 | Dimension 3 | Units of cm. Cylinder = height, RPR = side 3 or c axis, ellipsoid = major axis, hollow cylinder = outer height, potato = z^{+}, octahedron = z^{+}. |

8 | 4 | Dimension 4 | Units of cm. Potato = x^{−} |

8 | 5 | Dimension 5 | Units of cm. Potato = y^{−} |

8 | 6 | Dimension 6 | Units of cm. Potato = z^{−}. Note: x^{+} = x^{−} = y^{+} = y^{−} = z^{+}, z^{−} = 0 for a hemisphere. |

9 | Constraints | Constraints on the dimensions | |

If calc mode = 0, then 0 = fit and 1 = fix. Volume remains constant. | |||

If calc mode = −2, then vary one of the following dimensions: 1 = x^{+}, 2 = y^{+}, 3 = z^{+}, 4 = x^{−}, 5 = y^{−}, 6 = z^{−}, 7 = $\psi $, 8 = $\theta $, 9 = $\varphi $. | |||

10 | Comments | Headers for the Euler angles in line 11 | |

11 | 1 | Euler angle 1 | Angle $\psi $ for rotations about the z axis in Kock's convention |

11 | 2 | Euler angle 2 | Angle $\theta $ for rotations about the new y axis in Kock's convention |

11 | 3 | Euler angle 3 | Angle $\varphi $ for rotations of ($\pi $-$\varphi $) about the new z axis in Kock's convention |

12 | Constraints | Constraints on the Euler angles. 0 = fit, 1 = fix | |

If calculation mode = −1, then value equals max deviation of Euler angle in degrees. | |||

If calculation mode = −2, then value corresponds to dimensions or angle to vary: dimensions (1–6) or Euler angles (7 = phi, 8 = theta, psi = 9). | |||

13 | Comments | Headers for the list of frequencies in line 14 | |

14+ | 1 | Frequencies | Measured frequencies in units of MHz |

14+ | 2 | Frequencies | Calculated frequencies in units of MHz. Calculated by the forward problem. |

14+ | 3 | Weights | Weight factors for each measured frequency. See Eq. (20). |

Several parameters are available at the top of the input file. In the third line, these parameters include five shapes, all crystal systems, polynomial order, calculation mode, optional mirror planes, sample mass, rms error from the previous fit (if applicable), and the option to calculate the derivative of each eigenfrequency with respect to each elastic constant, i.e., $df/dcxx$. Input files for each calculation mode are provided along with the *RUScal* code. Two additional parameters can be added: an integer *n* that sets the stopping tolerance for chi-squared variations (10^{−n}; default *n* = 5) and an integer *e* that controls the error bar calculation (0: fast estimation; 1: rigorous—default).

After execution of the of the *RUScal* fitting program (*#-calc* = 0), three text (.dat) files are generated, two of which, *rusio* and *rusout*, are formatted like those generated by the *RPR*.exe and *Cyl*.exe codes (Migliori *et al.*, 1993). The format of *rusio* is identical to that of *rusin* except the free parameters take on values obtained from the last iteration of the fit, and the column of calculated frequencies (F_th) is updated based on the new fitted parameters. The *rusio* file is convenient when one requires multiple fits.

The *rusout* file contains all inputs and fitted data along with several calculations of material properties, as shown in Fig. 6. Note the elastic constants $Cijkl$ are given in the two-index Voigt notation: 11 → 1, 22 → 2, 33 → 3, (12, 21, 23, 32) → 4, (13, 31) → 5, and (12, 21) → 6, where the subscripts *i*, *j, k, l =* 1, 2, 3 correspond to the coordinate axes *x*, *y*, *z*, respectively. Both the Voigt and Reuss values (Cook and Young, 1999) for the bulk $K$ and shear moduli $G$ are provided, which correspond to the upper and lower bounds of these moduli, respectively, i.e., $KV\u2264K\u2264KR$ and $GV\u2264G\u2264GR$,

In Eqs. (22a)–(22d), $sij$ are the components of the compliance matrix, and $cij$ are the elastic constants in the Voigt notation. The true values of the bulk and shear moduli lie between the Voigt and Reuss extremes and often close to the arithmetic or Hill averages (Hill, 1952): $Kavg=KV+KR/2$ and $Gavg=GV+GR/2$. Poisson's ratio $\mu $ as well as the longitudinal $VL$ and transverse $VT$ speeds of sound are calculated from the Hill averages in Eqs. (23a)–(23c):

The average speed of sound $Va$ is given by the following (Anderson, 1963):

The last text file (*summary*.dat) contains a list of fit results that have been ordered in rows, one for each fit that produced a *rusout* file. This summary file is provided to aid analysis and plotting of multiple fit results, such as during a series of temperature-dependent measurements, where the tracking parameter is the independent variable listed in the first column of the file.

#### 1. Forward calculation

Section II A details the mathematical methods used to calculate the eigenfrequencies of normal modes of a free resonator. Using the *RUScal* code, one may generate a number of these frequencies using known properties of a sample, including the shape, density, and approximate elastic constants. To begin, set the *#-calc* parameter in the input file to an integer value greater than zero, which is equal to the number of calculated frequencies. Upon execution of the code, the calculated frequencies will appear both on the active command line and within a newly generated text file named *freqout_calc*.

#### 2. Monte Carlo and genetic algorithm tools

One typically chooses elastic constants that approximate the solution as is required for convergence in the least squares minimization. The elastic constants of many common materials, such as metals, ceramics, and semiconductors, are readily found in published literature or may be obtained separately from experiments, such as pulse-echo ultrasound, neutron, or Brillouin scattering, or from theoretical calculations [see, e.g., the Materials Project (Jain *et al.*, 2013)]. However, it is possible, especially in the case of novel materials, that such data are insufficient or difficult to obtain. We address this issue by providing a Monte Carlo-type routine to probe a wide landscape of possible solutions wherein the values of selected elastic constants are chosen randomly within user-established bounds and the corresponding chi-squared errors are evaluated.

To activate this mode, the user prepares the input (*rusin*) file as usual, providing educated guesses for the elastic constants, and the value of the *#-calc* parameter (see Table I) is set to –1. Then each constraint on the elastic constants is replaced by a value that is the maximum deviation (in per-thousand) from the initial value; this is demonstrated in Fig. 7(a) for a sample of LiYF_{4} (abbreviated as “YLF”), which is a solid-state lasing material with well-known elastic constants [see, e.g., Blanchfield and Saunders (1979)], and five constants are allowed to vary within 20% of their starting values. Note the use of comments in the input file for recordkeeping. Upon execution, two text files named *rms_out* and *rms_sort* are generated and contain a series of rows, each corresponding to an iteration of the chi-squared evaluation of a set of elastic constants, where the latter is sorted from lowest to highest chi squared; the rms is reported in the output files. The sets of elastic constants that correspond to the lines with the lowest rms values should then be used for fitting. We advise fitting with multiple sets of solutions to check that they all converge to the same true minimum—see Sec. II C 4 below for more detail.

Occasionally, two or more elastic constants may become anti-correlated by coincidence of specific dimensions of the sample, which is confirmed by the near identical dependence of each resonant frequency on each elastic constant, i.e., the d*f*/d*c*_{xx} are the same with opposite sign. In this case, it is not straightforward to fit those constants. This occurred for the sample of LiYF_{4} associated with Fig. 7, so the MC tool was used to search for $c33$ and $c23$ values that give the lowest chi-squared error. Figure 7(b) shows the chi-squared values from three consecutive runs of the MC tool, noting that each subsequent run approaches a set of solutions associated with the lowest chi squared. The correlated elastic constants, *c*_{33} and *c*_{23}, were fixed to values measured by Ghani (2013), which agree with those published by other groups (e.g., Blanchfield and Saunders, 1979), and the remaining constants were allowed to vary within 20 ‰ of the starting values. After running the MC tool, we selected about five sets of elastic constants with the lowest chi squared and used each set to fit the data and confirm that each set converged to the same solution, i.e., the solution is robust and lies in the minimum of parameter space.

Supplementary to the Monte Carlo routine above is the option to employ an additional genetic algorithm (“GenAlg”) that mimics the Darwinian theory of natural selection to search for and evaluate generations of solutions iteratively. To use the GenAlg sequence, the value of the *#-calc* parameter (see Table I) is set to −5. The start is the same as the MC program above with the candidate parent population of elastic constants generated via Monte Carlo (500 individuals by default), which are then sorted by the fitness, or chi-squared value, of each parent. The next offspring generation is then created as hybrids between two parents. The parents are chosen at random, with an evolutionary pressure favoring the fittest parents. In addition, there is a small probability of offspring mutation for each specific gene, i.e., individual elastic constant. The results are sorted and serve a population for the next generation. After the number of chosen generations, the algorithm ends with a complete fit using the fittest offspring as starting point. Progress can be monitored by plotting the *rms_sort*.dat file.

The genetic algorithm utilizes default parameters that can be tuned by adding five numbers (two integers and three floats) at the end of the first parameter line in the input file. The population can be set to any integer number above 100 (default: *pop* = 500); the number of generations can be set to any number larger than or equal to 2 (default: *n_gen* = 10); the evolutionary pressure that biases the choice of fittest parents can be set to 1.5 or greater (default: *ev* = 3); the mixing of parents can use strict interpolation or allow for extrapolation (0 ≤ mix ≤ 0.5; default: *mix* = 0.2); and the probability of mutation can be tuned between 0.0/*nc* and 5.0/*nc* (default: *mut_prob* = 1.0/*nc*), were *nc* is the number of free elastic constants. Note that the algorithm does not generally preserve an elite from the parent generation. An elite (*n* individuals) is preserved only if *all* offspring are less fit than these *n* individuals. The equations for choosing and mixing parents *p1* and *p2* to give offspring *o* from the population sorted by rms are

with *r* a random number 0 ≤ *r* ≤ 1, and {*i*} denotes choosing the *i*th element from the population,

with *β* a random number in the range 0-*mix* ≤ *β* ≤ 1+*mix.* If *r < mut_prob*, then

with *r* a random number 0 ≤ *r* ≤ 1, *o _{i}* the

*i*th elastic constant,

*gen*the current generation, and

*range*the deviation range for this constant in the initial Monte Carlo search.

_{i}#### 3. Dimension evaluation

The shape of the sample affects the frequencies of resonances in a way that could cause modes to overlap or become correlated, such as in the case encountered above for YLF. To address this issue, we included the dimension evaluation program, which calculates the eigenfrequencies of resonances systematically as a function of sample length or crystal orientation. Note that only one of the dimensions or Euler angles is varied at a time, and the density is kept constant. To execute this program, set the *#-calc* value to –2, and change the constraint on the dimension to a value between 1 and 9 (see Table I for details). Figure 8 shows both the *rusin* file and results of varying a RPR sample's length. In a practical sense, we find this tool useful during sample preparation to tailor the ultrasound spectrum to have, e.g., wide gaps in frequency between modes, enabling easier mode identification. For example, Fig. 8(b) shows that by cutting the steel sample's height (dimension 3), a gap opens between the first and second modes among others. The generated output text files include *sequence_out* and *sequence_matrix*.

#### 4. Error analysis

While the rms error is important for assessing the precision of a fit, it does not, however, represent the true accuracy of the model, i.e., the error bar of the elastic moduli. Such an analysis must account for systematic errors, such as uncertainty in the measured mass, dimensions, crystal orientation (if applicable), and positions of resonant frequencies. Migliori and Sarrao (1997) outline a set of criteria to achieve an acceptable fit:

The number of modes should be at least 5 times, preferably 8–10 times, the number of free parameters, and the number of missing modes should be less than 10% of the number of measured modes.

The rms error should be less than 0.8%, such as in the 0.2%–0.5% range, but preferably below 0.2%.

Fitted and measured dimensions and Euler angles should not differ by more than a few

*μ*m (e.g., <10*μ*m) and 1°–2°, respectively.The fit should be robust in that changing the starting parameters by a couple percent leads to the same solution.

The model should be able to predict the positions of modes that are not used in fitting, such as a missing mode.

In the absence of a phase transition, the temperature dependence of the elastic moduli should be as smooth as the shift in the measured frequencies.

Once a fit has been obtained, the shape and depth of the minimum in $\chi 2$ versus $cxx$ for each elastic constant can be probed using the included error analysis program. To initiate this program, set the *#-calc* parameter to −3 in the *rusin* file and run *RUScal*. The program will change each elastic constant individually by ±*n*^{2}*0.01% of its starting (fitted) value with *n* between 0 and 5 and calculate the corresponding chi squared. The results are printed in the text file *errorout*. The error on each elastic constant may be estimated by the change in magnitude of each constant that corresponds to a 2% increase in chi squared (Migliori and Sarrao, 1997).

In addition to the error analysis program, a grid search tool is available, by which the chi-squared error is evaluated systematically for each elastic constant (up to three). This mode was included with the intention of probing the parameter landscape around, e.g., adjacent minima of a fit. To use this program, set the *#-calc* parameter to –4, and change the constraint on each elastic constant. The constraint is expressed in per-thousand change relative on each side of the starting value, in a grid of 21 steps. As noted, this program allows one to vary up to three elastic constants by entering nonzero values for the respective constraints.

### D. Graphical user interface

The graphical user interface (GUI) *RUS Analysis* has been provided to help users prepare input files and execute the *RUScal* code. The GUI is written in python using Tkinter libraries compatible with Windows 10, Mac OS, and Linux operating systems. Figure 9 shows the GUI laid out in a 1280 × 1024 frame available with all input parameters and calculation modes. Users must ensure all selections have been made, especially shape and Bravais class, as these affect other options. The check boxes belonging to “calculate derivatives,” Euler angles, dimensions, and elastic constants, when checked, correspond to values of 1 for each parameter in the input file or zero if unchecked. The “bounds” entry boxes for each elastic constant take integer values for the Monte Carlo (–1) and Grid (–4) modes. The list of resonant frequencies must be prepared in a three-column format (measured lines—calculated lines—weight factor) and saved in comma separated (csv), Excel, or text (txt) file format. If the calculated or measured lines are unknown, then substitute their values with zeros with care to avoid missing lines. Once the file has been selected and loaded, the frequencies table will populate, and each line may be selected and modified. Alternatively, the GUI parameters may be automatically populated with values from a prepared input file by clicking either the “Choose input file” or “Read input file” button. The latter will read the *rusio* file located in the same directory as the GUI.

## III. BENCHMARKING RUSCAL ANALYSIS SOFTWARE

The analysis code was tested against published data for materials with various sizes, shapes, and crystal properties that were published with a list of measured resonant frequencies. While we provide capabilities to analyze octahedral shapes and triclinic crystals, to our knowledge there are no published data with which we could have reproduced results and validated the code.

### A. Anisotropic spheres

The first set of examples tests the software for elastically anisotropic single-crystal spheres, both rutile (TiO_{2}) and periclase (MgO), which exhibit tetragonal and cubic crystal symmetries, respectively. Suzuki *et al.* (1992) provided the diameter (8.194 ± 0.001 mm) and density (4.2539 g/cm^{3}) of their rutile TiO_{2} sample as well as a list of the first 43 eigenfrequencies measured via the RST method. Because some modes have the same stiffness, they share the same frequency and become degenerate. Thus, the total corresponding number of calculated modes is 57. Six independent elastic constants are used to describe the elastic properties of TiO_{2}: *c*_{11}, *c*_{33,} *c*_{23}, *c*_{12}, *c*_{44}, and *c*_{66}. Here, we fit the reported data using *N* = 14 polynomial order, which resulted in a rms error of 0.01% and the elastic constants shown in Table II in excellent agreement with results from the authors.

c_{11} (GPa)
. | c_{33} (GPa)
. | c_{23} (GPa)
. | c_{12} (GPa)
. | c_{44} (GPa)
. | c_{66} (GPa)
. | Ref. . |
---|---|---|---|---|---|---|

269.13 | 480.83 | 146.15 | 176.05 | 123.88 | 193.42 | RUScal |

267.59 | 482.84 | 148.09 | 176.42 | 123.43 | 192.85 | Suzuki et al. (1992) |

c_{11} (GPa)
. | c_{33} (GPa)
. | c_{23} (GPa)
. | c_{12} (GPa)
. | c_{44} (GPa)
. | c_{66} (GPa)
. | Ref. . |
---|---|---|---|---|---|---|

269.13 | 480.83 | 146.15 | 176.05 | 123.88 | 193.42 | RUScal |

267.59 | 482.84 | 148.09 | 176.42 | 123.43 | 192.85 | Suzuki et al. (1992) |

The 130.785 mg periclase sample from Oda *et al.* (1994) exhibits cubic symmetry described by three independent elastic constants, *c*_{11}, *c*_{12}, and *c*_{44}, and measures 4.1135 mm in diameter for a density of 3.5886 g/cm^{3}. The authors provided the first 22 eigenmodes, many of which are degenerate, for a total of 51 modes, measured via the RST method. Using *N* = 14 polynomial, our results are in excellent agreement with the authors as shown in Table III with a rms error of ∼0.02%. These results also agree with those from similar studies by Sumino *et al.* (1976) and Sumino *et al.* (1983) and Anderson and Andreatch (1966).

$\rho $ (g/cm^{3})
. | c_{11} (GPa)
. | c (GPa)
. _{12} | c (GPa)
. _{44} | Ref. . |
---|---|---|---|---|

3.589 | 296.58 ± 0.02 | 95.12 ± 0.02 | 156.03 ± 0.01 | RUScal |

3.587 | 296.60 ± 1.30 | 95.90 ± 1.40 | 156.20 ± 0.30 | Sumino et al. (1976) |

3.588 | 298.50 ± 0.80 | 97.50 ± 0.90 | 156.70 ± 0.20 | Sumino et al. (1983) |

3.588 33 | 296.47 | 95.07 | 155.89 | Anderson and Andreatch (1966) |

3.589 | 296.15 ± 0.13 | 95.35 ± 0.12 | 155.89 ± 0.05 | Oda et al. (1994) |

$\rho $ (g/cm^{3})
. | c_{11} (GPa)
. | c (GPa)
. _{12} | c (GPa)
. _{44} | Ref. . |
---|---|---|---|---|

3.589 | 296.58 ± 0.02 | 95.12 ± 0.02 | 156.03 ± 0.01 | RUScal |

3.587 | 296.60 ± 1.30 | 95.90 ± 1.40 | 156.20 ± 0.30 | Sumino et al. (1976) |

3.588 | 298.50 ± 0.80 | 97.50 ± 0.90 | 156.70 ± 0.20 | Sumino et al. (1983) |

3.588 33 | 296.47 | 95.07 | 155.89 | Anderson and Andreatch (1966) |

3.589 | 296.15 ± 0.13 | 95.35 ± 0.12 | 155.89 ± 0.05 | Oda et al. (1994) |

### B. Irregular shape

The sphere represents the shape with the highest symmetry. However, here we seek to test the calculation for the highly irregular, although analytical, “potato” shape by Visscher *et al.* (1991). Six semiaxes (+*x*, +*y*, +*z*, −*x*, −*y*, and −*z*) are used to construct the potato, whose volume is integrated over ellipsoid segments in each octant of the 3D Cartesian coordinate space using the set of basis functions in Eq. (18). The example we use is that of the massive aluminum half-potato in Visscher *et al.* (1991) that weighs 2438.4 g and measures 5.08, 10.16, 7.62, 12.7, 2.54, and 0 cm along each semiaxis in the order given above. The authors provided the first ten modes plus the 20th and 30th modes in Table I of Visscher *et al.* (1991). Using *N* = 9, we achieved a similar rms error of ∼0.15% with the elastic constants shown in Table IV. Upon fitting the elastic constants, we estimated the time to calculate 30 eigenfrequencies to be between less than 1 s for *N* = 6 and ∼6 s for *N* = 9 using a single processor (2.60 GHz) on a local laptop computer.

c_{11} (GPa)
. | c_{44} (GPa)
. | rms error . | Ref. . |
---|---|---|---|

108.53 ± 2.55 | 26.56 ± 0.06 | 0.15% | RUScal |

108 | 27 | 0.16% | Visscher et al. (1991) |

c_{11} (GPa)
. | c_{44} (GPa)
. | rms error . | Ref. . |
---|---|---|---|

108.53 ± 2.55 | 26.56 ± 0.06 | 0.15% | RUScal |

108 | 27 | 0.16% | Visscher et al. (1991) |

### C. Hollow cylinder

The hollow cylinder geometry was tested using data for the silicon carbide (SiC) tube by Donahue *et al.* (2019). Their monolithic SiC tubular specimen was obtained directly from Dow Chemical Company (Midland, MI) from vapor deposited SiC, weighing 0.8844 g and measuring 1.603 cm long with outer and inner diameters of 0.852 and 0.712 cm, respectively. The sample is assumed to be acoustically isotropic, and its elastic properties are described by two independent elastic constants: *c*_{11} and *c*_{44}. Donahue *et al.* solved the forward problem using the FEM, whereas they approached the inverse problem like that of Remillieux *et al.* (2015) by using a genetic algorithm to find the global minimum of the cost function—the weighted difference between measured and calculated eigenfrequencies, similar to Eq. (20) above. Using the provided data and *N* = 13, we find a reasonable fit with a rms error of 0.638% between measured and calculated modes with the elastic moduli shown in Table V. The computation time was ∼1 min using the energy minimization method with a single processor (2.6 GHz) on a laptop computer. Note that at the time of writing this paper, a polynomial order greater than *N* = 13 generates a numerical error. We suspect the root cause is a rounding error during matrix manipulation of the hollow cylinder geometry that creates negative numbers during diagonalization. This may be corrected in the future by increasing the precision using long doubles or libraries like GNU Science library to perform matrix manipulations.

c_{11} (GPa)
. | c (GPa)
. _{44} | E (GPa)
. | K (GPa)
. | rms error . | Ref. . |
---|---|---|---|---|---|

466.36 ± 1.35 | 203.1 ± 0.8 | 452.61 | 195.4 | 0.638% | RUScal |

455.7 ± 3.7 | 206.8 ± 2.9 | 448.5 | 180 | 0.793% | Donahue et al. (2019) |

c_{11} (GPa)
. | c (GPa)
. _{44} | E (GPa)
. | K (GPa)
. | rms error . | Ref. . |
---|---|---|---|---|---|

466.36 ± 1.35 | 203.1 ± 0.8 | 452.61 | 195.4 | 0.638% | RUScal |

455.7 ± 3.7 | 206.8 ± 2.9 | 448.5 | 180 | 0.793% | Donahue et al. (2019) |

### D. Trigonal

The elastic constants of single-crystal corundum (α-Al_{2}O_{3}) are well known and offer the means to test the *RUScal* code for crystals with trigonal symmetry. The data by Ohno *et al.* (1986) were used in this test as this is the first known study using the RPR method to determine the elastic constants of a trigonal crystal. Their sample is a gem-quality crystal measuring 0.3335 ± 0.0004 (parallel to *a* axis), 0.2789 ± 0.0005, and 0.6663 ± 0.0007 cm (parallel to *c* axis) with a mass of ∼0.2468 g determined from the provided density (3.982 g/cm^{3}). The lattice structure of corundum exhibits $R3\xafc$ space group and 3 *m* point group symmetry; thus, its elastic properties are described by six independent elastic constants: *c*_{11}, *c*_{33}, *c*_{23}, *c*_{12}, *c*_{44}, and *c*_{16}. Using *N* = 12 and the 39 eigenmodes provided in Ohno *et al.* (1986), we achieved a rms error of 0.09% with the corresponding elastic constants in Table VI. Error bars were estimated from the scatter in the elastic constants from multiple fits with the same rms error. Results are also compared with those from Wachtman *et al.* (1960), Goto *et al.* (1989), and Gladden (2003), all of whom obtained $cij$ from resonance methods with general agreement.

c_{11} (GPa)
. | c_{33} (GPa)
. | c_{44} (GPa)
. | c_{12} (GPa)
. | c_{13} (GPa)
. | c_{14} (GPa)
. | Ref. . |
---|---|---|---|---|---|---|

499.51 ± 0.95 | 501.78 ± 0.97 | 147.29 ± 0.15 | 164.90 ± 0.97 | 117.29 ± 1.25 | −22.58 ± 0.29 | RUScal |

496.8 ± 1.4 | 500.5 ± 1.6 | 146.8 ± 0.2 | 162.3 ± 1.6 | 115.5 ± 1.6 | −21.9 ± 0.2 | Ohno et al. (1986) |

496.8 ± 1.8 | 498.1 ± 1.4 | 147.4 ± 0.2 | 163.6 ± 1.8 | 110.9 ± 2.2 | −23.5 ± 0.3 | Wachtman et al. (1960) |

497.3 ± 1.5 | 500.9 ± 1.8 | 146.8 ± 0.2 | 162.8 ± 1.7 | 116.0 ± 1.8 | −21.90 ± 0.24 | Goto et al. (1989) |

496.5 | 503.4 | 146.2 | 159.3 | 119.0 | −22.6 | Gladden (2003) |

c_{11} (GPa)
. | c_{33} (GPa)
. | c_{44} (GPa)
. | c_{12} (GPa)
. | c_{13} (GPa)
. | c_{14} (GPa)
. | Ref. . |
---|---|---|---|---|---|---|

499.51 ± 0.95 | 501.78 ± 0.97 | 147.29 ± 0.15 | 164.90 ± 0.97 | 117.29 ± 1.25 | −22.58 ± 0.29 | RUScal |

496.8 ± 1.4 | 500.5 ± 1.6 | 146.8 ± 0.2 | 162.3 ± 1.6 | 115.5 ± 1.6 | −21.9 ± 0.2 | Ohno et al. (1986) |

496.8 ± 1.8 | 498.1 ± 1.4 | 147.4 ± 0.2 | 163.6 ± 1.8 | 110.9 ± 2.2 | −23.5 ± 0.3 | Wachtman et al. (1960) |

497.3 ± 1.5 | 500.9 ± 1.8 | 146.8 ± 0.2 | 162.8 ± 1.7 | 116.0 ± 1.8 | −21.90 ± 0.24 | Goto et al. (1989) |

496.5 | 503.4 | 146.2 | 159.3 | 119.0 | −22.6 | Gladden (2003) |

### E. Miscut crystal

As indicated in Sec. II, an occasion may arise, such as during sample preparation, where the crystallographic axes do not align with the edges of the sample body. In such a case, the elastic tensor needs to be rotated to generate an effective elastic tensor for the specific orientation. In his Ph.D. dissertation, Gladden (2003) sought to determine the sign of the *c*_{14} constant by intentionally miscutting an alumina crystal. His samples consisted of two RPR specimens cut from the same bulk crystal obtained from the Goodfellow Corp. (Boulder City, NV). The measured dimensions of sample 1 and sample 2 are (728.7, 1068.1, and 830.8) *μ*m and (764.7, 1032.2, and 810.2) *μ*m, respectively—sample 2 is used for the calculations below. The +*x* and +*z* directions point along the crystal *a _{h}* (twofold rotational symmetry) and

*c*axis (threefold rotational symmetry). The mass of 2.549 mg was calculated from the given mass and density (3.986 g/cm

^{3}). As the testing of trigonal symmetry is handled primarily in Sec. II, here the focus is on the calculation of eigenfrequencies that depend on the sign of

*c*

_{14}. The starting elastic constants (in GPa) are

*c*

_{11}= 497.3,

*c*

_{33}= 500.9,

*c*

_{23}= 116.0,

*c*

_{12}= 162.8,

*c*

_{44}= 146.8, and

*c*

_{14}= ±21.9. Using the measured frequencies and

*N*= 9, the forward calculation gives the results in Table VII, showing good agreement with Gladden, including some lines that were missed during measurement.

Measured . | Calc. (+c_{14}) (Gladden, 2003) (MHz)
. | Calc. (this work) (MHz) . | Calc. (−c_{14}) (Gladden, 2003) (MHz)
. | Calc. (this work) (MHz) . |
---|---|---|---|---|

2.707 688 | 2.749 469 | 2.749 52 | 2.705 586 | 2.705 63 |

3.672 555 | 3.649 729 | 3.649 79 | 3.680 767 | 3.680 83 |

3.835 520 | 3.867 467 | 3.867 53 | 3.836 131 | 3.836 20 |

3.975 381 | 3.914 977 | 3.915 04 | 3.977 453 | 3.977 52 |

4.169 249 | 4.122 196 | 4.122 27 | 4.155 786 | 4.155 86 |

4.530 27 | 4.669 34 | |||

4.797 780 | 4.627 694 | 4.627 77 | 4.774 484 | 4.774 56 |

4.942 362 | 4.947 461 | 4.947 54 | 4.929 546 | 4.929 63 |

5.086 565 | 5.045 944 | 5.046 03 | 5.097 846 | 5.097 93 |

5.161 420 | 5.139 197 | 5.139 28 | 5.208 598 | 5.208 69 |

5.257 03 | 5.275 40 | |||

5.305 400 | 5.435 047 | 5.435 14 | 5.310 701 | 5.310 79 |

5.323 344 | 5.475 603 | 5.475 70 | 5.347 387 | 5.346 48 |

5.501 625 | 5.479 522 | 5.479 61 | 5.512 956 | 5.513 05 |

5.561 940 | 5.684 228 | 5.684 32 | 5.576 348 | 5.576 44 |

5.739 862 | 5.957 796 | 5.957 90 | 5.734 950 | 5.735 05 |

5.965 31 | 5.806 21 | |||

5.866 725 | 6.153 921 | 6.154 03 | 5.842 906 | 5.843 00 |

6.235 88 | 6.058 46 | |||

6.370 460 | 6.460 665 | 6.460 77 | 6.370 622 | 6.370 73 |

6.488 406 | 6.486 284 | 6.486 39 | 6.462 084 | 6.462 19 |

Measured . | Calc. (+c_{14}) (Gladden, 2003) (MHz)
. | Calc. (this work) (MHz) . | Calc. (−c_{14}) (Gladden, 2003) (MHz)
. | Calc. (this work) (MHz) . |
---|---|---|---|---|

2.707 688 | 2.749 469 | 2.749 52 | 2.705 586 | 2.705 63 |

3.672 555 | 3.649 729 | 3.649 79 | 3.680 767 | 3.680 83 |

3.835 520 | 3.867 467 | 3.867 53 | 3.836 131 | 3.836 20 |

3.975 381 | 3.914 977 | 3.915 04 | 3.977 453 | 3.977 52 |

4.169 249 | 4.122 196 | 4.122 27 | 4.155 786 | 4.155 86 |

4.530 27 | 4.669 34 | |||

4.797 780 | 4.627 694 | 4.627 77 | 4.774 484 | 4.774 56 |

4.942 362 | 4.947 461 | 4.947 54 | 4.929 546 | 4.929 63 |

5.086 565 | 5.045 944 | 5.046 03 | 5.097 846 | 5.097 93 |

5.161 420 | 5.139 197 | 5.139 28 | 5.208 598 | 5.208 69 |

5.257 03 | 5.275 40 | |||

5.305 400 | 5.435 047 | 5.435 14 | 5.310 701 | 5.310 79 |

5.323 344 | 5.475 603 | 5.475 70 | 5.347 387 | 5.346 48 |

5.501 625 | 5.479 522 | 5.479 61 | 5.512 956 | 5.513 05 |

5.561 940 | 5.684 228 | 5.684 32 | 5.576 348 | 5.576 44 |

5.739 862 | 5.957 796 | 5.957 90 | 5.734 950 | 5.735 05 |

5.965 31 | 5.806 21 | |||

5.866 725 | 6.153 921 | 6.154 03 | 5.842 906 | 5.843 00 |

6.235 88 | 6.058 46 | |||

6.370 460 | 6.460 665 | 6.460 77 | 6.370 622 | 6.370 73 |

6.488 406 | 6.486 284 | 6.486 39 | 6.462 084 | 6.462 19 |

### F. Monoclinic and triclinic

To our best knowledge, analysis software based on energy minimization methods that are publicly available has been limited to orthorhombic and higher crystal symmetry. Here, we demonstrate the calculation of eigenfrequencies for the monoclinic system with the data provided in Isaak and Ohno (2003). Their sample is a rectangular parallelepiped shape of chrome-diopside (93% diopside) with a mass of 0.3191 g. The side lengths are 0.4749, 0.4169 (parallel to *b* axis), and 0.4905 cm (parallel to *c* axis). The authors calculated the eigenfrequencies of this sample by using the RPR method with Legendre polynomials as the basis functions along with the reduction scheme by Ohno (1976) to invert the frequency spectrum to obtain the elastic constants. Using *N* = 12 and the 62 modes provided in the authors' Table V, we find excellent agreement with a rms error of ∼0.1% and results provided in Table VIII.

Elastic constant . | RUScal
. | Isaak and Ohno (2003) . | Levien et al. (1979)
. |
---|---|---|---|

c_{11} (GPa) | 228.040 ± 0.005 | 228.1 ± 1.0 | 223 |

c_{22} (GPa) | 181.110 ± 0.004 | 181.1 ± 0.6 | 171 |

c_{33} (GPa) | 245.350 ± 0.007 | 245.4 ± 1.3 | 235 |

c_{23} (GPa) | 61.131 ± 0.007 | 61.1 ± 0.7 | 57 |

c_{13} (GPa) | 70.271 ± 0.008 | 70.2 ± 0.7 | 81 |

c_{12} (GPa) | 78.811 ± 0.003 | 78.8 ± 0.5 | 77 |

c_{44} (GPa) | 78.869 ± 0.003 | 78.9 ± 0.3 | 74 |

c_{55} (GPa) | 68.145 ± 0.002 | 68.2 ± 0.2 | 67 |

c_{66} (GPa) | 78.072 ± 0.003 | 78.1 ± 0.2 | 66 |

c_{15} (GPa) | 7.877 ± 0.007 | 7.9 ± 0.5 | 17 |

c_{25} (GPa) | 5.807 ± 0.006 | 5.9 ± 0.5 | 7 |

c_{35} (GPa) | 38.756 ± 0.003 | 39.7 ± 0.4 | 43 |

c_{46} (GPa) | 6.421 ± 0.004 | 6.4 ± 0.2 | 7.3 |

K (GPa) | 116.45 ± 2.95 | 116.5 ± 0.9 | 113 |

G (GPa) | 72.72 ± 1.92 | 72.8 ± 0.4 | 67 |

Elastic constant . | RUScal
. | Isaak and Ohno (2003) . | Levien et al. (1979)
. |
---|---|---|---|

c_{11} (GPa) | 228.040 ± 0.005 | 228.1 ± 1.0 | 223 |

c_{22} (GPa) | 181.110 ± 0.004 | 181.1 ± 0.6 | 171 |

c_{33} (GPa) | 245.350 ± 0.007 | 245.4 ± 1.3 | 235 |

c_{23} (GPa) | 61.131 ± 0.007 | 61.1 ± 0.7 | 57 |

c_{13} (GPa) | 70.271 ± 0.008 | 70.2 ± 0.7 | 81 |

c_{12} (GPa) | 78.811 ± 0.003 | 78.8 ± 0.5 | 77 |

c_{44} (GPa) | 78.869 ± 0.003 | 78.9 ± 0.3 | 74 |

c_{55} (GPa) | 68.145 ± 0.002 | 68.2 ± 0.2 | 67 |

c_{66} (GPa) | 78.072 ± 0.003 | 78.1 ± 0.2 | 66 |

c_{15} (GPa) | 7.877 ± 0.007 | 7.9 ± 0.5 | 17 |

c_{25} (GPa) | 5.807 ± 0.006 | 5.9 ± 0.5 | 7 |

c_{35} (GPa) | 38.756 ± 0.003 | 39.7 ± 0.4 | 43 |

c_{46} (GPa) | 6.421 ± 0.004 | 6.4 ± 0.2 | 7.3 |

K (GPa) | 116.45 ± 2.95 | 116.5 ± 0.9 | 113 |

G (GPa) | 72.72 ± 1.92 | 72.8 ± 0.4 | 67 |

To our knowledge, there are no published RUS measurements of a crystal in the triclinic system for which the eigenfrequencies are provided publicly. Regardless, the triclinic symmetry has been included in *RUScal* and left to the community for testing.

### G. Mirror planes

Isotropic solids and singles crystals with orthorhombic or higher symmetry have translational invariance in the Cartesian coordinate system, i.e., *x* → −*x*, *y* → −*y*, and *z* → −*z* with three mutually perpendicular mirror planes that eliminate many off diagonal terms of the elastic stiffness matrix. This symmetry can be broken by, e.g., the introduction of a chiral axis for a single crystal or preferred orientation of crystallites (texture) in a polycrystal. The *RUScal* code allows zero, one, two, or all three orthogonal mirror planes to accommodate various scenarios. For example, a polycrystal with fiber-like texture that is aligned parallel to a sample axis can be described with a hexagonal stiffness matrix with a mirror plane normal to the sixfold axis. Mirror planes help to speed up the calculation of eigenfrequencies by dissecting the large *R* × *R* matrix into smaller blocks for faster diagonalization. Therefore, one should endeavor to find acoustical symmetry planes whenever possible. During fits of new or unknown materials, we advise using no mirror planes initially and then introducing them once acoustical symmetry is confirmed, e.g., by rotating the crystal to align the symmetry plane normal to a body axis. Acoustically equivalent orientations have identical resonance spectra. Warnings are given if both mirror planes and angular rotations are used, wherein the user should verify the spectra are identical with and without mirror planes.

## IV. CONCLUSION

Open-source software for the analysis of RUS measurements named *RUScal* has been developed in C language based on Lagrangian minimization and Rayleigh–Ritz methods. In comparison to other publicly available software, new features have been added that include numerical calculations for hollow cylinder and the irregular (but analytical) potato shapes as well as orthorhombic, trigonal, monoclinic, and triclinic crystal systems with and without mirror planes. Furthermore, multiple analysis tools are added, including Monte Carlo-based programs to help search for and evaluate elastic constants or Euler angles, two modes to assess the chi-squared error of a fit, and a calculation program that evaluates the effect of sample dimensions on its eigenfrequencies. Except for triclinic crystals and the octahedron shape, the analysis code was tested against published data for materials with various shapes, sizes, and crystal symmetries, thus confirming their results and validating the software. Details on practices and procedures were given to prepare new users for measurements using the *RUScal* code and the RUS experiment technique.

## ACKNOWLEDGMENTS

This work is funded by the Laboratory Directed Research and Development Program of Oak Ridge National Laboratory, managed by UT-Battelle, LLC, for the United States Department of Energy. We acknowledge Dr. Alaska Subedi for contributing Cholesky decomposition to the initial code development, Dr. A. Balodhi and Professor A. Zevalkink for providing data sets on LiYF_{4}, Professor V. Keppens and the D. Mandrus group for lending data to test the code, and Dr. D. Pierce for providing a 4140-steel sample. We thank Dr. R. Juneja and Dr. A. Biswas for helpful discussions on the implementation of the genetic algorithm and Dr. Watkins and Dr. Lara-Curzio for useful comments.

## References

_{4}

_{2}Si

_{2}uncovered by resonant ultrasound spectroscopy and machine learning

_{1.86}Sr

_{0.14}CuO

_{4}

_{2}Cu

_{3}O

_{6+δ}