Previously, the infrared permittivity tensor of monoclinic β-Ga_{2}O_{3} crystals has been determined using ellipsometry reflection measurements from two differently oriented monoclinic β-Ga_{2}O_{3} crystals with surfaces parallel to (010) and (−201). The (010) surface places the crystallographic a-c plane in the table of the instrument. The permittivity tensor consists of four complex values, and in order to compute it, four or more combinations of measurements are required at selected table rotations and incidence angles. However, the (010) orientation also places the transverse optical (TO) modes with Au symmetry parallel to the z-axis of the instrument, and we find that these modes are not fully excited and, hence, not measurable due to underlying selection rules. This makes additional measurements on surfaces other than (010) necessary. The second orientation has been the (−201) crystal, which places the crystallographic b axis in the plane of the table to access the transverse Au phonons. In prior work, the overall tensor has been determined by combining measurements of the two crystal orientations [Schubert *et al.*, Phys. Rev. B **93**, 125209 (2016)]. The goal of the work here is to find single crystal orientations for which all TO modes can be determined from measurements. The use of a set of measurements employed for such a single crystal is inextricably linked to the choice of incidence angles and table rotations. Consequently, determining suitable angles for these is linked to the selection of a crystal orientation, which is, therefore, an integral part of the overall goal. The TO contribution to the permittivity strongly dominates at or near the TO mode wavenumber resonances and, therefore, are used in this work to identify suitable orientations for a single crystal. Any such crystal orientation will also provide measurements useful to compute permittivity across the entire measured wavenumber range. In principle, any crystal orientation that does not place the direction of any TO mode at or near the z-axis may be suitable due to the underlying physics and mathematics of the problem. We discuss which of these measurement parameters contain the most sensitivity for the (111) orientation. For accuracy, we seek the best or very good orientations. Our investigation follows a previously demonstrated approach where at a single wavelength, the full tensor of an orthorhombic absorbing crystal was obtained from a low-symmetry surface of stibnite [Schubert and Dollase, Opt. Lett. **27**, 2073 (2002)]. We discuss which of these measurement parameters contain the most sensitivity for the (111) orientation. The methods presented here will also be useful for other monoclinic materials as well as other materials of different crystal structures, including orthorhombic and triclinic materials.

## I. INTRODUCTION

Recent interest in crystalline materials with monoclinic symmetry is a result of their potential for application in electronic power devices, scintillators, high-power lasers, frequency stable laser local oscillators, light slowing and trapping devices, and optical quantum memory technologies.^{1–7} The material under investigation here is the monoclinic phase (β) of gallium oxide (β-Ga_{2}O_{3} or bGO), which is thermodynamically stable and can be grown by bulk and epitaxial methods.^{8} This particular material has potential for application in switches and transistors capable of operating at 10 000 or more volts and large current carrying capacities.^{2} Thus, it is a competitor for technologies based on silicon carbide and gallium nitride. It may be useful in electric vehicles as well as the transport of electric power. Consequently, measurement and understanding of the material properties is necessitated.

The topic of this work concerns optically nondepolarizing single crystals of monoclinic β-Ga_{2}O_{3}. The part of the spectrum under examination is in the infrared ranging from 300 to 1100 (1/cm) wave numbers, which corresponds to a wavelength range from 33 × 10^{3} to 9.1 × 10^{3} nm. Ellipsometry is the method used to determine the optical properties. In this work, we focus on the transverse optical (TO) modes at their respective wavenumbers at which the permittivity is strongly dominated by the mode contributions. The material is anisotropic; therefore, the permittivity is not a scalar but a 3 × 3 tensor. This permittivity tensor appears in Maxwell's equations, which are used to model the relationship between the electromagnetic field in the material and in the surrounding medium. The permittivity of the crystal may be solved by applying numerical methods to mathematical models using Maxwell's equations. The model relies upon formulating the equations for a reflection or transmission geometry, and this has been carried out by a number of authors.^{8–12} Berreman described a convenient formalism he evolved from early work in radar to treat the case of plane polarized light incident on a sample.^{8} We use this model, which is continuum mechanical in that it makes no explicit reference to the crystal microstructure on the electronic, ionic, or photonic levels. It is only concerned with the incident and reflected waves as they relate to the permittivity of the crystal.

Previously, Schubert used measurements for two separate crystals, (010) and (−201), at five azimuths (table rotations) and three incidence angles.^{13,14} In this seminal work, the values of the monoclinic permittivity tensor were determined using measurements from both crystals by simultaneously solving using the Levenberg–Marquardt numerical method. The process first solved one wavenumber at a time and then followed with additional processing across the wavenumber. The result was the determination of all four complex unknown components of the permittivity tensor in a least-squares sense. Details of this development are fully described elsewhere.^{13} It was shown previously that two types of phonons with Au and Bu symmetry exist in β-gallium oxide. The symbols refer to the irreducible presentation of the subgroups in C2h (the monoclinic symmetry class). There are four such groups: Au, Bu, Ag, and Bg, where the latter two are Raman active. Modes with Au symmetry (hereafter termed Au modes) are polarized parallel to the lattice b axis, and modes with Bu symmetry (hereafter Bu) are polarized within the monoclinic a-c plane; see Fig. 1. These groups of modes are identified by different symmetry operations under specific replacements of coordinates. Strain and stress relationships for optical phonon modes in monoclinic crystals of β-Ga_{2}O_{3} are an example.^{15} Features in the permittivity tensor elements, determined as a function of wavelength on a wavenumber-by-wavenumber approach, then permit lineshape analysis to quantify frequency, amplitude, broadening, and direction information within the monoclinic lattice.

Preceding work by the authors has shown that solutions for the permittivity tensor from a single crystal, i.e., from measurements taken from a single surface of a monoclinic crystal, are possible but may not reveal all phonon modes due to the highly anisotropic optical nature of the material. There are spectral regions within which the experimental data contain no sensitivity to the existence of phonon modes. For example, the Au modes are not seen using measurements on the (010) crystal regardless of the angle of incidence or instrument table rotation. On the other hand, these Au modes are measurable in the case of the (−201) surface. This has led to the idea that two crystal orientations are required to measure all 12 TO modes. Indeed, if measurements from (010) and (−201) are combined in a common analysis, then all modes can be found.^{13}

The purpose of the work presented here is to determine how to obtain all TO modes using a single crystal sample and to determine which crystal orientations are more favorable for accuracy. We examine the usefulness of the crystal crystallographic plane at which the crystal is cut by considering the reflection intensity and then sensitivity of measured ellipsometric parameters (Ψ and Δ) to the tensor components. We select measurements and test for suitable mathematical condition of the Jacobian used in the solution process. Finally, it is possible to examine the expected optical response of various differently oriented crystals. The case of (111) surface orientation is shown here, for example, to determine if all of the TO phonon modes are selected and appear with sufficient intensity, sensitivity, and condition for measurement.^{16} Selection of a suitable crystal orientation for ellipsometry measurements for a low-symmetry crystal depends on a great many factors, such as the presence of TO phonon modes, their resonant frequencies, and their orientations within the crystal, for examples. Determining a suitable crystal orientation will, thus, depend upon the material of interest and all of the considerations included in this work. To accomplish this, it is necessary to determine useful incidence angles and table rotations, and these selections are also described here. Finally, it is true that any particular set of measurements have to have a useful solution of the Maxwell equations. Equations with good mathematical condition can be solved as opposed to a set of equations with poor condition. Thus, this topic is also evaluated and discussed. Therefore, there is no simple rule for the selection, but the steps in making such a selection are described in this work.

## II. THEORY

### A. Crystal orientation

Figure 1 shows the unit cell of the β-Ga_{2}O_{3} monoclinic crystal with a set of three axes (a, b, and c) separated by three angles (α, β, and γ). Of these angles, γ = 90° is the angle between a and b, α = 90° is the angle between b and c, and β ≠ 90° lies between a and c. Thus, the b axis is normal to the plane containing the a and c axes. The only axis of symmetry is b, and there is mirror symmetry in a–c.^{17} In the figure, the axes x, y, and z define the permittivity reference axes with respect to the crystal. The permittivity is calculated with reference to these axes.

The Miller indices of the plane under examination with reference to the crystal axes (hkl) correspond to the integer inverse of the intersections of the plane with the axes in the a, b, c order. A plane defined in this way must be rotated appropriately from the standard position given in Fig. 1 to a position that places the plane of interest parallel to the instrument table. This rotation must also be used to transform the permittivity tensor in the model. For consistency with prior work, we select the Euler angles to specify the rotations required. Thus, the rotation is ** Z_{3} X_{2} Z_{1}** (corresponding to φ, θ, ψ) in the sequence of the subscripts as applied to the standard position to obtain the position for measurement.

The values of ** Z** and

**are standard rotation matrices shown below. In this work, the convention for the sign of the rotation angle is defined as positive for a rotation that is counterclockwise when viewed by an observer looking along the axis of rotation and viewing toward the origin,**

*X*The rotation expression is not commutative, sequence matters. The 3D rotation is needed in order to transform the permittivity tensor to conform to the position of the crystal in measurement in which the cut face lies flat on the instrument table.

It can be seen in Fig. 2 that the (010) plane lies in the table (X–Y plane) and that no rotational transformation is required. The situation is different for the (−201) plane cut crystal. Figure 3 shows the crystal in its standard orientation with the (−201) plane in color and its normal as a black arrow.

This figure shows the crystal and the instrument axes in the same position as that shown in the immediately preceding Fig. 2. In order for the permittivity tensor of a crystal of this orientation to be solved, it must be transformed to be consistent with the crystal rotation, which places the plane of interest parallel to the instrument table. This is shown in Fig. 4.

Figure 4 shows the crystal rotated into the orientation, which places the (−201) crystal face parallel to the XY plane, which is required to make a measurement on that face. Taking care to perform the rotations in the proper sequence, the values required to transform the permittivity tensor are ψ = −49.92° followed by a rotation about the new x axis of θ = 90°. Similar rotations are used for other crystal orientations of interest.

### B. Reflection model

Monochromatic light that is incident on a crystal interacts with it and results in reflection and transmission. In the present work, we restrict ourselves to linear materials for which the optical response does not depend upon electric field strength of the incident wave. The ellipsometry method and mathematics for the light reflection and transmission are well described elsewhere.^{8–12} For the present discussion, the derivation of Berreman is used as previously noted.^{8} The continuum mechanical model based on the Maxwell equations and a selected geometry for the reflection results in expressions that relate the permittivity of the crystal to ellipsometry measurements. Measurements may be the Jones matrix elements or the Mueller matrix elements. In this simple model, a two-phase, anisotropic substrate-isotropic ambient is included, where no surface overlayers are considered. Surface overlayers are often required to account for effects of small-scale roughness, e.g., due to polishing processes and/or surfactants but are considered irrelevant at infrared wavelengths due to the smallness of their effective optical thickness. For this particular dataset, we selected the data identified by the ellipsometer as having the least estimated experimental error, which is reported in the same file as the measurements. In other appropriate cases, the Mueller matrix data could be used.

The incident light causes a reflection from the crystal by which permittivity can be solved considering boundary conditions between the surrounding medium and the crystal. This model provides the necessary mathematical relationship between known parameters, including the light incidence angle, light wavelength, light reflection measurements, crystal orientation, and the sought-for unknown parameters consisting of the values of the permittivity tensor. More detail is provided by Berreman.^{8}

It is important that we review a point concerning the permittivity reference axes (x, y, z) fixed in the crystal with respect to the crystal axes (a, b, c); see Fig. 1. First, let us consider the representation of the ɛ tensor with respect to the (x, y, z) axes. The general permittivity tensor for a monoclinic crystal for which we solve is expressed relative to the permittivity crystal axes by

The tensor is symmetric so that ɛ_{xy} = ɛ_{yx}. We refer to this representation as the standard reference tensor consisting of diagonal values and an off-diagonal value. Here, the b crystal axis corresponds to the z direction, and the x direction is aligned with the a crystal direction. Thus, the remaining axis (y) is not aligned with the c crystal axis but with an axis c* instead, as the angle between a and c is not 90° but rather 103.7°.^{5} Physical rotation of the crystal with respect to the measurement axes requires that the standard reference tensor be transformed appropriately, mathematically, for the model to correspond to the physical orientation.

Because each TO mode corresponds to a distinct direction within the crystal structure, it is clear that the exciting electric field at the appropriate frequency must have a component (or generate a component) in a direction capable of influencing the corresponding lattice vibration, which is measurable and useful. This is the reason that not all combinations of light incidence angles, crystal orientations (with respect to characteristic measurement directions), and sample rotations (azimuths) are equally capable of elucidating a given mode. This is more extensively described in our prior work.^{16} The background in this area is extensive and far too lengthy to be repeated here in its totality but instead can be accessed through the citations.

### C. TO modes

In this work, we examine the TO modes of which there are four of the Au type and eight of the Bu type. The Au type mode directions are parallel to the b axis within the crystal, and the Bu type modes all have directions lying in the crystal a–c plane.^{13} Figure 5 shows the mode normals overlaid upon the crystal in the standard orientation. In prior work, these modes were observed experimentally at specific wavenumbers that were found to be in reasonable agreement with quantum mechanical modeling and theoretical calculations of long-wavelength active Γ-point phonon frequencies performed by plane wave density functional theory (DFT) using Quantum ESPRESSO (QE).^{13}

Note the difference between the concepts of normal directions of phonon displacement (a.k.a. mode normal) and the concepts of normal electromagnetic propagation modes (a.k.a. normal mode). The latter applies to electromagnetic waves, which follow from the solutions of the wave equation as derived from Maxwell's postulates. On the other hand, the set of eigenvalue/eigenfunction equations does not deliver the direction under which a certain phonon mode is polarized. This direction of a mode normal is the result of the physical arrangement of matter in local and global atomic and molecular structures. First principles methods, such as density functional theory, can predict such properties from computational methods. The rationale following our previous work is to identify whether or not permittivity values can be determined from given sets of measurements taken on different surface cuts of monoclinic symmetry β-gallium oxide. It was found that permittivity solutions using only (010) measurements did reveal all tensor components within the monoclinic a-c plane, and these did not reveal ɛ_{zz} values. Thus, while all Bu modes were found, none of the Au modes could be seen. The wavenumber-by-wavenumber, numerically determined values for ɛ_{zz} carried infinite error bars in the vicinity of TO modes since there was no sensitivity to their physical behaviors within the measured data. In other words, one can insert practically any random value for ɛ_{zz} during the model calculations and that will have no effect on the other solved permittivity parameters. The physical cause of this phenomenon, sometimes also considered analog to selection rules for TO modes in anisotropic materials, can be found simply in Snell's law.^{18} As the frequency approaches the TO resonance, the index of refraction for electric field components perpendicular to the surface reaches very large values. The boundary conditions require the normal component of the displacement to be continuous across the interface. With the index of refraction being unity at the ambient side, the electric field component approaches zero within the anisotropic material when the frequency reaches toward TO. Thus, the reflected wave contains no information in the limit toward the TO mode about the TO mode itself. Hence, one cannot identify via model calculations the amplitude and center resonance of a TO phonon whose displacement direction is perpendicular to a given surface. Accordingly, analysis of measurements from the (−201) surface revealed almost all Au and Bu modes except for those now oriented close to the surface normal, e.g., Bu(3) at 572 (1/cm) and Bu(5) at 357 (1/cm). Some of the other Bu modes that did appear were quite noisy. It is clear that due to these effects, ellipsometry measurements using either of these crystal orientations will not provide permittivity spectra, which completely permit identification of all existing phonon modes (Figure 6).

## III. EXPERIMENT

### A. Crystals

Bulk single crystals of β-Ga_{2}O_{3} were grown at Novel Crystal Technologies, Tokyo, Japan (formerly part of the Tamura Company) using the edge-defined film-fed growth method.^{13} Substrates of different crystallographic orientations, (010) and (−201), were cut from the bulk crystals and then polished on one side and oriented within a few tenths of a degree. The resulting samples were 10 × 10 mm^{2} with a thickness of 0.65 mm. The substrates were doped with Sn to an estimated activated electron density of N_{d}−N_{a} ≈ (2–9) × 10^{18} cm^{−3}. Cut crystals are analyzed by XRD prior to shipping. The orientation angles obtained in the ellipsometry analysis for axis b are in full agreement with the XRD data information of the two samples. In this work, the angle ψ refers to the first rotation about the initial Z axis as shown, for example, in Fig. 2, which rotates the crystal from its initial standard position; see Eq. (1). The angle θ refers to the second rotation about the new X axis (resulting from the first rotation) to bring the plane of interest to be parallel to the table, X_Y, plane; see Eq. (2).

### B. Measurements

The original measurements by Schubert *et al.* were made across two different spectral ranges using two different ellipsometers. Details of these measurements are given in a prior publication.^{13} All measurements were taken directly from data files output from the ellipsometer, including Ψ and Δ for p-p, p-s, and s-p. Mueller matrix measurements are also available in the instrument-generated data files but were not used in this work. The generalized ellipsometry parameters are simply defined in terms of Jones matrix elements in reflection,

Errors in the measurements due to instrumental factors arise and are reported along with the measurements as estimated experimental errors.

### C. Solutions

The method for solving for the permittivity tensor is briefly described here and fully described elsewhere.^{11} Computations are carried out at each wavenumber using the trust-region dogleg algorithm employing measured Ψ_{pp} and Δ_{pp}. This algorithm finds a numerical solution while seeking a match to 10^{−6} between the four complex measured and computed data at each wavenumber. A set of nine complex measurements at a selected wavenumber are picked from which four are taken at a time for a total of 126 combinations. The algorithm finds values of the permittivity tensor of Eq. (4), which consists of four complex numbers. If between 20 and 40 solutions are found, these are recorded and the algorithm moves on to the next wavenumber. If fewer or more than between 20 and 40 are found, the stopping condition of function tolerance is adjusted up or down accordingly by a factor of 1.1. Adjustments in the function tolerance continue until the target number of found solutions falls between 20 and 40.

## IV. RESULTS AND DISCUSSION

### A. Determination of a single crystal orientation choice

The first step in selecting a single crystal orientation, which potentially would reveal all the modes of interest, is to examine the angle between the mode normals and the Z axis of the measurement table. If the angle is 0, they are parallel, the sine, thus, is 0 and that is “bad,” and if perpendicular, the angle is 90°, so the sine is 1, which is “good.” This step is necessary but not sufficient. The first two angles of an Euler rotation (ψ and θ in that order) can bring any crystal plane from the initial standard orientation shown in Fig. 1 so that the specified crystal face is parallel to the measurement table. Thus, the mode normal to the Z axis angle as a function of crystal cut can be found beginning with the standard crystal orientation in Fig. 1, rotating for a range of values of ψ and θ, and computing the angles. Each mode is computed individually, and an example is shown for a Bu 3 TO mode in Fig. 7. The orientation of (010) is shown for reference. The magnitude of the sine of the angle between instrument Z and mode normal is plotted according to the colorbar reference to the right in the figure. Thus, the center of the dark region corresponds to parallel (sine = 0) and the light to perpendicular (sine = 1). We shall term this kind of dataset a “mode rotation array.”

Figure 7 is periodic due to the periodicity of angle space and is plotted thusly to allow convenient display, which may include negative angles. Note the circle labeled “(010),” which indicates the rotational combination, which places an (010) plane parallel to the instrument table, X-Y plane. A number of other rotations would also accomplish that, for example, ψ = 0 and θ = 0, and there are others also not noted in the figure. Because the BU modes for this orientation all lie in the X-Y plane as can be seen in Fig. 5, this kind of pattern is repeated for every BU TO mode for a (010) cut.

By contrast, Fig. 8 shows the mode rotation array for all Au modes, which lie at the same orientation, as seen in Fig. 5, to be aligned with the Z axis. Interestingly, the plot of the magnitude of the sine of the angle between Au mode normal and Z axis is not a function of table rotation (equivalent to φ) as can be readily understood from Fig. 5. This alignment has been shown by experiment to result in the Au modes not affecting experimental data measured from the (010) orientation where measurement is insensitive to ɛ_{zz}. The object here is to select a point in the θ and ψ plane of all mode rotation arrays, which will provide useful experimental data.

It has been observed that modes Bu 3 and Bu 5 are also insensitive to the permittivity components in solution for the (−201) crystal orientation. Examination of the mode rotation array in Fig. 9 shows that the Bu 5 mode normal is nearly parallel to the Z axis as can be seen clearly in the figure.

Finding the optimum crystal cut for mode solutions can be accomplished using the following approach. First, select the modes of interest by a chosen criterion. For example, the criterion could be that all the modes fall within the measured wavenumber range. Then, compute mode rotation arrays for all modes of interest, as shown, for example, for three mode rotation arrays in Figs. 7–9. Since mode amplitudes differ, examination of all modes in the aggregate takes this difference into account by weighting each by its amplitude to produce a set of mode rotation array “footprints.” These new footprint arrays are then normalized and expressed as a percentage of individual footprint maximum to give each equal standing in determining the usefulness of crystal orientations. A different weighting could be chosen on a different basis than equal standing.

Each mode rotation array footprint has a unique value at each orientation expressed as the intersection of θ and ψ. Thus, if all 12 TO modes are considered, at each point in the ranges of θ and ψ, one of these modes will have the smallest footprint. The importance of this minimum is that at each corresponding orientation, all other modes will have a greater footprint. Thus, if the minimum footprint is sufficiently visible, all the other modes will be more visible at that orientation. Figure 10 shows a plot of the overall minimum of the normalized footprints of all Bu and Au modes. We shall call the values of this set of minima the orientation “score” in percent arising due to normalization. It is possible but not necessary to optimize this choice as it is clearly evident that many choices will be sufficient. An optimal choice by this numerical method is (101); however, it does not guarantee sensitivity to be evaluated following. The overall unsuitable (010) and (−201) are seen for which the footprints are zero or near zero. Note that (111) has a score of about 40%; yet, there are better choices with greater values. It now becomes necessary to examine the sensitivities.

The choice of orientation can be made on other considerations, such as availability of the crystal cut. Some possible nonoptimal but potentially sufficient examples are shown in Table I, which include (111), (212), and (210) as well as the not useful (010) and (−201). These are chosen for illustration purposes only. The table lists the sine of the angle between each mode and the instrument z-axis. A sine of zero indicates that the measurement is not sensitive to that mode. A sine that is close to zero indicates that a signal-to-noise ratio problem can be expected, which would result in a noisy solution. Note that for Au modes, values are summed over all modes.

Mode . | Wavenumber (1/cm) . | (010) . | (−201) . | (111) . | (212) . | (210) . |
---|---|---|---|---|---|---|

sin(κ) . | sin(κ) . | sin(κ) . | sin(κ) . | sin(κ) . | ||

Au all | 663, 449, 297, 155 | 0 | 1 | 0.5447 | 0.7924 | 0.4555 |

Bu(1) | 743 | 1 | 0.9993 | 0.8585 | 0.6659 | 0.9261 |

Bu(2) | 692 | 1 | 0.7093 | 0.9673 | 0.9302 | 0.8928 |

Bu(3) | 572 | 1 | 0.5581 | 0.9046 | 0.7845 | 0.9998 |

Bu(4) | 432 | 1 | 0.8753 | 0.927 | 0.8381 | 0.8921 |

Bu(5) | 357 | 1 | 0.07111 | 0.9919 | 0.9827 | 0.9556 |

Bu(6) | 279 | 1 | 0.6956 | 0.97 | 0.9354 | 0.8935 |

Bu(7) | 262 | 1 | 0.3186 | 0.9999 | 0.9999 | 0.9281 |

Bu(8) | 214 | 1 | 0.8574 | 0.9999 | 0.637 | 0.9843 |

Mode . | Wavenumber (1/cm) . | (010) . | (−201) . | (111) . | (212) . | (210) . |
---|---|---|---|---|---|---|

sin(κ) . | sin(κ) . | sin(κ) . | sin(κ) . | sin(κ) . | ||

Au all | 663, 449, 297, 155 | 0 | 1 | 0.5447 | 0.7924 | 0.4555 |

Bu(1) | 743 | 1 | 0.9993 | 0.8585 | 0.6659 | 0.9261 |

Bu(2) | 692 | 1 | 0.7093 | 0.9673 | 0.9302 | 0.8928 |

Bu(3) | 572 | 1 | 0.5581 | 0.9046 | 0.7845 | 0.9998 |

Bu(4) | 432 | 1 | 0.8753 | 0.927 | 0.8381 | 0.8921 |

Bu(5) | 357 | 1 | 0.07111 | 0.9919 | 0.9827 | 0.9556 |

Bu(6) | 279 | 1 | 0.6956 | 0.97 | 0.9354 | 0.8935 |

Bu(7) | 262 | 1 | 0.3186 | 0.9999 | 0.9999 | 0.9281 |

Bu(8) | 214 | 1 | 0.8574 | 0.9999 | 0.637 | 0.9843 |

Table I lists the angle between mode normals and the instrument z-axis for five separate crystal orientations. All four Au modes are at the same angle, and values given are sums for all modes. The Bu modes are listed individually. Mode wave numbers are shown to three digits of accuracy. Note that two of the modes that are within the data range, Bu(3) and Bu(5), do not appear in (−201) data solutions; hence, we suggest that the threshold for measurements is somewhere between 0.5 and 0.6. The other modes appear in measurements on (−201). The other three orientations are all well above this threshold with the possible exception of (212) Bu(1). Table I column (010) lists the sine of the angle from each TO mode to the instrument Z axis. The Au mode angles are 0. This does not mean that the values of ɛ_{zz} are zero, but that they are not determinable, in the vicinity of the TO mode. The sines of all of the Bu mode angles are 1, and examination of (010) solutions shows all Bu modes.^{17} Column 2 contains the sines of the angles for the (−201) crystal plane. The sine of the angle for the Au modes is maximum at 1. Examination of (−201) solutions shows the appearance of the three Au modes, which lie in the range of wave numbers in the measurement. At the same time, the values of the sines of the angles for the Bu modes have changed significantly. Examination of (−201) for these modes reveals that not all of the Bu modes appear.^{17} Several mode contributions do not give rise to features in the ellipsometry measured data above the noise floor, and thus, computational methods do not detect the presence of a phonon mode. Modes Bu(3) and Bu(5) are absent in the solutions as they correspond to small factors in the table. We note that in principle, these modes with small sine values (meaning, modes that are almost parallel to the surface normal) could be measured if the measurement accuracy were improved so that their response would appear distinctly above the measurement error shown in part in Fig. 11.

Figure 11 shows the experimental error to be the best for wavenumbers greater than about 400 with the exception of a 100 wavenumber band around 700 wavenumbers. Measurements for wavenumbers less than 400 show rather significant estimated experimental errors on the order of the measurement itself. This explains the challenge of getting accurate measurements for wave numbers other than in the low-error ranges. A consequence is that outside of the best ranges, the solution algorithm is computing with noise.

### B. Simulation of (111) measurements

Finally, we arrive at an examination of a newly selected crystal orientation, which is not optimal and is chosen because the Au modes appear to be above the margin of noise. Also, we select the (111) orientation knowing beforehand that the crystal surface will be oriented such that none of the crystals a, b, and c axes are either within or perpendicular to (111); see Fig. 13.

The effect on the location in space of the TO phonon mode dipole orientation due to rotation as shown in Fig. 12 is shown in Fig. 13 for visualization. It can be seen in this figure that none of the mode normals are parallel to the table Z axis. Consequently, all modes should be detectable in data measured by ellipsometry under various angles of incidence and table rotations as was predicted in Fig. 10, which also shows orientations with better scores.

To determine to what degree permittivity solutions would reveal all modes from measurements on a (111) crystal, we first need simulated measurements obtained considering the tensor and measurement geometry. In prior work, we used table rotations of 0°, 45°, 90°, and 135° along with incidence angles of 50°, 60°, and 70°. The initial purpose here is to discover if Au and Bu modes can be determined from theoretical measurements at these table rotations and incidence angles. The issue of improving the choice of initial rotation on the table will be dealt with the following.

It is then a simple matter to compute theoretical measurements using this tensor transformed to correspond to the crystal rotation of (111) and the Maxwell equations as expressed by Berreman.^{8} At each wave number, 12 different measurements were computed for different incidence angles and table rotations. Finally, we perform solutions using the computed measurements, four at a time from a group of nine as already described. The only difference in the solution process is due to the absence of experimental error. For confirmation, we sought and found all 126 solutions rather than a number in a range for 20–40 solutions.

### C. Sensitivity to table rotations and incidence angles

Once the crystal cut is determined, the next question is what incidence angles and table rotations to select for the set of measurements. This is far less urgent because they can be adjusted in the ellipsometer and do not require obtaining a new crystal. The first consideration is the intensity of the reflection. As has been shown in our prior work, intensity of the reflection is sufficient for practically all wave numbers for the pp and ss reflections and insufficient for the ps and sp reflections. Consequently, we select the Ψ_{pp} and Δ_{pp} measurements, which rely on these pp and ss reflections. In this notation, p refers to the component of the E field, which is polarized parallel to the plane of incidence and s refers to the E field component perpendicular to the plane. The order of the subscripts refers to the polarization of the incident and reflected light, respectively. The opposite sequence does appear occasionally in the literature.

The next consideration is the sensitivity of the measurements to the underlying components of the permittivity tensor (ɛ_{xx}, ɛ_{yy}, ɛ_{zz}, and ɛ_{xy}). Computing sensitivities with respect to the real and imaginary parts of the four permittivity components shown in Eq. (4) results in eight three-dimensional plots of sensitivity as functions of table rotation, φt, and incidence angle, α, at each wavenumber. For this, we consider just the wavenumbers for the TO modes. Sensitivity is obtained by computing the derivatives of Ψ_{pp} and Δ_{pp} with respect to each of the eight parameters to be solved. Given the choices of subdivisions of the angles, this gives rise to 2 × 8 × 12 = 196 plots (Ψ_{pp} Δ_{pp}, tensor components, modes). Each plot consists of 80 × 180 = 14 400 points for a total of >2.8 × 10^{6} points.

What is needed is to select at least nine combinations of incidence angle and table rotation from which to select four at a time. In prior work, we actually had 12 combinations (φt = 0°, 45°, 90°, and 135° by α = 50°, 60°, and 70°). The goal is to choose a new set of 9 or more, which could be an improvement over those listed just above. The approach used is as follows. First, a sensitivity limit is selected. Observation of the individual plots suggested 1% as a sensitivity limit from which to start. It can be seen that most modes and parameters result in sensitivities ranging from 0.1% or so to greater than 1%. It stands out that the sensitivities of the permittivity components associated with the Au modes (ɛ_{zz}) are many orders of magnitude lower, on the order of 10^{−4}%, which correlates with the previous finding that these components were not able to be computed for measurements on (010). Figure 14 is a plot of sensitivity for the Bu1 mode with the real part of ɛ_{xx} shown as contours. There are 97 other plots like this to include the 12 modes for each of 8 permittivity components.

In order to examine all 196 figures, a method of consolidation is required. First of all, we will consider one mode at a time, thus requiring the analysis of eight figures like that in Fig. 14. The method employed is to first choose a cut-off sensitivity and set all points in the sensitivity plot so that all points below are zero. The result is a figure that indicates a region that does not meet the criteria and that is being recorded as zero, and the remaining region that meets or exceeds the criteria is set to unity. The figure of this kind is shown in Fig. 15, which displays Fig. 14 with a minimum sensitivity set to 3% as an example.

For illustrative purposes only, Fig. 15 displays the result of limiting sensitivity to 1%. This is restrictive as it is only an example, but it can be seen that greater than 1% can be obtained for a range of table rotations from approximately 50°–100° at a 50° angle of incidence and for greater than 60° angle of incidence sensitivity is less than 1%.

It is useful to consider all relevant sensitivities one mode at a time. To accomplish this, a new data set is formed by choosing a cut-off sensitivity and then setting all points, which fall below it to zero and all points above it to 1. Thus, the range of sensitivity applicable to any selected group, such as the one considered here, which is all parameters in mode Bu1, can be found by multiplying these arrays together. The result is Fig. 16, which shows the region within the plot for which all parameter sensitivities for mode Bu 1 exceed 1% for a (111) crystal. The narrow stripes are regions for which the sensitivity for the remaining six parameters is less than 1% and regions greater than 1%. It is concluded that a wide range of incidence angles and table rotations will acquire useful data for mode Bu 1 from a (111) crystal.

### D. Condition of measurements in a solution

Finally, it becomes necessary to determine if Au and Bu modes can be solved from sets of four measurements at particular table rotations and incidence angles by determining the mathematical condition of sets of eight real number equations and eight unknowns resulting from the four complex measurements. Thus, the usefulness of sets of four measurements can be determined. It is already understood that placing all of these four measurements at the same incidence angle or at the same table rotation very likely provides a poor mathematical condition. Each individual mode responds differently to these angles; in fact, the response of the material at every wavenumber is different. To avoid the necessity of examining over a thousand sets of combinations of table rotation and incidence angle, we examine only at the previously used table rotations and incidence angles as an example. In practice, the angles determined in the sensitivity analysis above would be used to produce a similar study.

There are 126 combinations of sets of 4 measurements taken from 9 measurements. To check if a given set of four can potentially be solved, we compute the condition at solution for which there exists a Jacobian matrix, which is an 8 × 8 matrix of all of the first-order partial derivatives of the equation parameters with respect to the tensor components. A simple test of the Jacobian reveals the probability of achieving accurate solutions because it must be well conditioned to facilitate its use in the numerical process. The “condition number for inversion” or the “condition number” is equal to the ratio of the largest to the smallest singular values of the Jacobian. A condition number near 1 indicates a well-conditioned matrix, and larger values indicate increasingly ill-conditioned matrices. A list of the combinations is provided in Fig. 17, which shows the 126 combinations across the horizontal axis and the selected measurements on the vertical axis. Table II lists the measurement conditions for the nine measurements from which four are repeatedly selected. As can be seen in Fig. 18, the condition varies significantly for each combination of four measurements.

Measurement number . | Incidence angle (°) . | Table rotation (°) . |
---|---|---|

1 | 70 | 0 |

2 | 70 | 45 |

3 | 70 | 90 |

4 | 70 | 135 |

5 | 60 | 0 |

6 | 60 | 45 |

7 | 60 | 90 |

8 | 60 | 135 |

9 | 50 | 0 |

Measurement number . | Incidence angle (°) . | Table rotation (°) . |
---|---|---|

1 | 70 | 0 |

2 | 70 | 45 |

3 | 70 | 90 |

4 | 70 | 135 |

5 | 60 | 0 |

6 | 60 | 45 |

7 | 60 | 90 |

8 | 60 | 135 |

9 | 50 | 0 |

The figure shows a visual of each of the 126 combinations of 4 measurements out of 9, which are used in Fig. 18. Each column displays four of the nine measurements listed from one to nine as shown in Table II. The first number of the horizontal axis is the “selection index,” which extends to 126, as shown in Fig. 18.

The condition also depends upon wave number, and therefore, mode wavenumbers only are presented here. Thus, Fig. 18 shows all 126 results for Bu modes for the example crystal orientation (111). In this figure, the logarithm to the base 10 is represented by as shading as shown in the bar on the right. Dark shading corresponds to usable condition, which fades to light shading as the highest condition number displayed.

A number of interesting features are immediately evident in the figure. First note that each horizontal line corresponds to a selection of four measurements from the nine. It is evident that there are many combinations for which the condition is usable. Thus, the figure identifies the measurements, which may be excluded in the solution process if desired, although at present, as 126 are attempted, so not all of the combinations are needed. The point here is that many combinations have a good condition.

Figure 19 shows all 126 results for Au modes for the example crystal orientation (111). In this figure, the logarithm to the base 10 is represented by a shading as shown in the shading bar. Dark blue corresponds to usable condition, which fades to light shading as the highest condition number displayed.

It is evident that there are many combinations for which the condition is usable. Thus, the figure identifies the measurements, which may be excluded in the solution process if desired.

These figures indicate that a large number of measurement combinations will be solvable in the (111) case. Considering the speed of parallel computing used here, the overload of attempting to solve using poor condition combinations did not present a problem. The time to identify which combinations not to use was greater than just trying each and letting the algorithm stopping conditions rule them out.

## V. SUMMARY AND CONCLUSIONS

The desired result of identifying potentially useful single crystal orientations that allow measurement and computation of all phonon modes of interest has been achieved. In the past, infrared permittivity tensors of monoclinic β-Ga_{2}O_{3} crystals have been determined using ellipsometry reflection measurements from two differently oriented crystals, (010) and (−201). The (010) orientation does not allow measurement of the Au TO modes even as it clearly places the Bu modes in their optimum orientations. The (−201) orientation places the Au modes in the optimum orientation; however, the Bu modes are not well oriented.

This work shows the steps to identify crystal orientations, which reveal all TO modes. At first, each mode normal must not be parallel with the ellipsometer Z axis. Examination of the rotation space, Euler ψ and θ starting from the standard orientation, is demonstrated for particular modes. A minimax solution is shown, which permits identification of suitable crystal orientations considering all TO modes. Next, for this crystal orientation, suitable choices of incidence angle, α, and table rotation, φ, are evaluated. This is shown as sensitivity of measurements (Ψ_{pp} and Δ_{pp}) to the real and imaginary parts of permittivity tensor components ɛ_{xx}, ɛ_{yy}, ɛ_{zz}, and ɛ_{xy} as a function of incidence angle, α, and table rotation, φt. Sensitivity is taken to be the first derivatives, for example, *d*Ψ_{pp}/*d* Real(ɛ_{xx}). The sensitivities can be examined individually for improved insight, examined one mode at a time, or examined for all modes under consideration. It is seen that useful values that can work in combination in the next step can be selected. Either way, these numbers can be easily changed once the crystal is in the instrument. Finally, we examine the computation feasibility using the condition of 126 combinations of four measurements taken from a subset of nine. For this, we examine all mode wave numbers in order to make a selection. Further experimental work is called for in order to confirm these findings. These methods will also be useful for other monoclinic materials as well as other materials of different crystal structures, including orthorhombic and triclinic materials.

## ACKNOWLEDGMENTS

We acknowledge support by the National Science Foundation (Award Nos. DMR 1808715 and OIA-2044049), by the Air Force Office of Scientific Research (Award Nos. FA9550-18-1-0360, FA9550-19-S-0003, and FA9550-21-1-0259), by the Knut and Alice Wallenbergs Foundation, by the University of Nebraska Foundation and the J. A. Woollam Foundation, by the Swedish Research Council VR (Award No. 2016-00889), by the Swedish Foundation for Strategic Research (Grant Nos. RIF14-055 and EM16-0024), by the Swedish Governmental Agency for Innovation Systems VINNOVA under the Competence Center Program (Grant No. 2016-05190), and by the Swedish Government Strategic Research Area in Materials Science on Functional Materials at Linköping University (Faculty Grant SFO Mat LiU No. 2009-00971). Special thanks to M. Stokey for Fig. 17.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

All three of the authors are university professors who worked closely together over a period of years to complete the work. The expertise and decades-long publication record of the authors is easily available by Internet search. It was a true team with each member contributing in his areas of expertise and in overlapping areas. These areas cover a broad range encompassing mathematics, physics, electromagnetics, materials science, crystallography, and numerical methods to name a few. In addition, each author contributed approximately equally to the ideas and into the writing. Thus, a simple rubric cannot give an accurate overview of the contributions. Hence, F. K. Urban III, D. Barton, and M. Schubert contributed to and took leadership in specific areas at different times in all aspects of the present work principally, including conceptualization, formal analysis, software, writing/original draft preparation, and writing/review and editing. F. K. Urban III, D. Barton, and M. Schubert also contributed equally at different times in all other areas.

**F. K. Urban, III:** Investigation (equal); Writing – original draft (lead); Writing – review & editing (lead). **D. Barton:** Investigation (equal); Writing – review & editing (supporting). **M. Schubert:** Writing – review & editing (supporting).

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.