{'fraction_non_alphanumeric': 0.031310923376853546, 'fraction_numerical': 0.015064689549240859, 'mean_word_length': 4.858560221504384, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Band-to-band transitions, selection rules, effective mass and exciton binding energy parameters in monoclinic \beta-Ga2O3

We employ an eigen polarization model including the description of direction dependent excitonic effects for rendering critical point structures within the dielectric function tensor of monoclinic \beta-Ga2O3 yielding a comprehensive analysis of generalized ellipsometry data obtained from 0.75 eV--9 eV. The eigen polarization model permits complete description of the dielectric response, and we obtain single-electron and excitonic band-to-band transition anisotropic critical point structure model parameters including their polarization eigenvectors within the monoclinic lattice. We compare our experimental analysis with results from density functional theory calculations performed using a recently proposed Gaussian-attenuation-Perdue-Burke-Ernzerhof hybrid density functional, and we present and discuss the order of the fundamental direct band-to-band transitions and their polarization selection rules, the electron and hole effective mass parameters for the three lowest band-to-band transitions, and their exciton binding energy parameters, in excellent agreement with our experimental results. We find that the effective masses for holes are highly anisotropic and correlate with the selection rules for the fundamental band-to-band transitions, where the observed transitions are polarized closely in the direction of the lowest hole effective mass for the valence band participating in the transition.

We employ an eigen polarization model including the description of direction dependent excitonic effects for rendering critical point structures within the dielectric function tensor of monoclinic β-Ga2O3 yielding a comprehensive analysis of generalized ellipsometry data obtained from 0.75 eV-9 eV. The eigen polarization model permits complete description of the dielectric response, and we obtain single-electron and excitonic band-to-band transition anisotropic critical point structure model parameters including their polarization eigenvectors within the monoclinic lattice. We compare our experimental analysis with results from density functional theory calculations performed using a recently proposed Gaussian-attenuation-Perdue-Burke-Ernzerhof hybrid density functional, and we present and discuss the order of the fundamental direct band-to-band transitions and their polarization selection rules, the electron and hole effective mass parameters for the three lowest band-to-band transitions, and their exciton binding energy parameters, in excellent agreement with our experimental results. We find that the effective masses for holes are highly anisotropic and correlate with the selection rules for the fundamental band-to-band transitions, where the observed transitions are polarized closely in the direction of the lowest hole effective mass for the valence band participating in the transition.

I. INTRODUCTION
Single crystalline group-III sesquioxides are currently at the forefront of research for applications in electronic and optoelectronic devices due to unique physical properties. Such conductive oxides, including tin doped In 2 O 3 or Ga 2 O 3 , can be utilized as transparent thin film electrodes for various devices such as photovoltaic cells 1 , flat panel displays 2 , smart windows 1,3 , and sensors 4 . The highly anisotropic monoclinic β-gallia crystal structure (β phase) is the most stable crystal structure among the five phases (α, β, γ, , and δ) of Ga 2 O 3 ( Fig. 1) 5,6 . It belongs to the space group 12 and has base center monoclinic lattice. Ga 2 O 3 shows potential for use in transparent electronics and high energy photonic applications due to its large band gap of 4.7-4.9 eV 7-10 .
Precise and accurate knowledge of the band gap energies, their polarization selection rules and energetic order, and the resulting anisotropy in the dielectric function is a prerequisite to correct assignment of basic physical properties in low symmetry materials. Electronic band-to-band transitions cause critical point (CP) features in the joint density of states, which result in CP structures in the dielectric function 11 . Accurate parametric models yield access to CP parameters such as transition energies and polarization selection rules, which allow for direct comparison with results both from exper- iment, e.g., optical absorption and reflectance measurements, as well as from theory, e.g., density functional theory (DFT) band structure calculations. Additionally, the accurate and explicit knowledge of the spectral dependence of the dielectric function tensor of monoclinic β-Ga 2 O 3 is equally important as it will allow for accurate and precise prediction of optical and optoelectronic arXiv:1704.06711v1 [cond-mat.mtrl-sci] 21 Apr 2017 device performances. A suitable and precise technique to determine the complex dielectric function of a material without a-priori implication of a physical model is ellipsometry 12 . An extension of this approach permitting for the appropriate concept to acquire experimental data from arbitrarily anisotropic materials has been developed recently by generalized spectroscopic ellipsometry (GSE) [12][13][14] . This concept has been used for accurate determination of the spectral dependencies of dielectric function tensor elements for a wide variety of materials with all known crystal symmetries 13,[15][16][17][18][19][20][21][22][23][24] .
The CP model approach was used frequently in the past to model dielectric function spectra as model dielectric function (MDF) approach. Within the MDF, individual contributions corresponding to, for example, band-to-band transitions and associated excitonic contributions, and contributions from spectral regions outside a given spectrum under consideration, sum to the overall dielectric response. The CP-MDF approach was also used successfully to explain anisotropy for many different types of materials, e.g., with tetragonal, hexagonal, and orthorhombic crystal symmetries in single crystal bulk and thin film heterostructures 13,[15][16][17][18][19]21,25 . Different model functions exist for different types of electronic excitations, e.g., single electronic transitions characterizable by the transition's Van-Hoofe singularity, and displacement and losses caused by many-body interactions (excitons) 11,26,27 . Functions for electronic transitions in the parabolic band approximation, for example, are described by Adachi, Kimura, and Suzuki 28 . Due to the complexity of the potential distribution in many body problems, excitonic contributions are often considered by Lorentzian or Gaussian broadened oscillators. Tangyu provided explicit solutions for Wannier-type excitonic contributions including all bound and unbound states in the parabolic band approximation [29][30][31] .
Dielectric anisotropic materials can be described by a symmetric dielectric tensor consisting of nine generally independent complex valued elements -complex valued wavelength dependent functions. In general, these dielectric function tensor elements will display variations with wavelength reflecting the underlying physical mechanisms responsible for the dielectric response of a given material. For example, at long wavelengths the dielectric response is governed by long wavelength active phonon modes and free charge carrier excitations. For nearinfrared (NIR), visible (VIS), ultra violet (UV) to vacuum ultra violet (VUV) photon energies, the dielectric response is caused by electronic band-to-band transitions and associated excitonic absorption processes. For cubic, tetragonal, trigonal, hexagonal, and orthorhombic crystal symmetries, physically meaningful transformations (spatial rotations with real valued angular arguments) can be found independent on wavelength (photon energy) which diagonalize the dielectric function tensor. Jellison et al., investigating the tensor elements of single crystals of monoclinic CdWO 4 by GSE in the NIR-UV spectral region, pointed out the necessity to include four inde-pendent dielectric tensor elements in order to describe a material with monoclinic crystal symmetry within spectral regions of absorption, equivalent to the fact that no physically meaningful rotations can be found for diagonalization independent of wavelength 32 . This fact has long been known for long wavelength analysis of materials with monoclinic and triclinic crystal symmetries. Born and Huang developed a microscopic description of the long wavelength lattice vibrations in the harmonic approximation, where the inter-atomic forces are considered constant and the equations of motion are determined by harmonic potentials 33  Fundamental band-to-band transitions in β-Ga 2 O 3 have been investigated very recently using a variety of techniques including density functional theory (DFT) calculations 9,10,44-46 , optical absorption 44 , reflection 47,48 , and ellipsometry 8,9 . Due to monoclinic symmetry, all band-to-band transitions are expected to depend on polarization. The eigen polarization of the transitions may not necessarily align with any of the high symmetry crystal axes. In addition, the formation of excitons, upon the optical excitation of a band-to-band transition, strongly modifies the frequency dependence of the dielectric response in semiconductors 11 . Thus, in order to accurately determine the transition energies, the excitonic contribution must be accounted for. For β-Ga 2 O 3 , there has been significant discrepancies in reported properties of the fundamental band-to-band transitions. Ricci et al., ignoring excitonic effects, recently showed optical absorption anisotropy in β-Ga 2 O 3 with the lowest onset of absorption occurring with polarization in the a-c plane at 4.5-4.6 eV 44 . For polarization along the crystal axis b the absorption onset was unambiguously shifted by 0.2 eV towards shorter wavelength. Onuma et al. investigated polarized transmittance and reflectance spectra 48 . As a result of their investigations, an indirect gap band-to-band transition around 4.43 eV and a direct gap transition around 4.48 eV parallel to the c axis were proposed, without considering excitonic effects. Sturm et al. considered contributions from both bound and unbound Wannier-type excitons 29,31,49 and reported the lowest direct gap band-to-band transition at approximately 4.88 eV, polarized within the a-c plane nearly parallel to c. Furthmüller and Bechstedt presented quasiparticle band structures and density of states of β-Ga 2 O 3 obtained from DFT combined with Hedin's GW approximation for single-particle excitations 10 . The lowest transition energy was determined by this approach to be around 5.04 eV, with polarization mainly along the c axis in the monoclinic a-c plane.
Setyawan and Curtarolo 50 demonstrated that materials with base center monoclinic lattices exhibit five distinct first Brillouin zones. Recently, Peelaers and Van de Walle 51 discussed the proper Brillouin zone and band structure of β-Ga 2 O 3 following the recommendations of Setyawan and Curtarolo. 50 High-symmetry points and lines connecting points in the Brillouin zone are of general interest for identifying physical properties of materials. Band structure computations for β-Ga 2 O 3 were presented in Refs. 9, 10, 44-46, 52-55, however, few of the authors specify the coordinates of high symmetry points and lines. Lines between high symmetry points are not necessarily high symmetry lines. In particular, the coordinates specified in previous theory works define highsymmetry points whose connecting lines do not present high symmetry lines. As pointed out by Peelaers and Van de Walle 51 most band structures published to date for β-Ga 2 O 3 do not consider all high-symmetry points, and include very few high symmetry lines, thus missing significant physics of the material. In the present manuscript we are expanding the results presented by Peelaers and Van de Walle 51 , and we study the effect of the choice of the lattice vectors onto the Brillouin zone and the resulting band structure.
The dielectric tensor of β-Ga 2 O 3 has been studied computationally by many authors. The most recent studies on the subject 9,10,44 use the following procedure: (1) ground state calculations at the generalized gradient approximation density functional theory (GGA-DFT) level, followed by (2) calculation of the wavefunction using a hybrid functional mixing a fraction of the Hartree-Fock exchange (HF) into the GGA-DFT; (3) the hybrid wavefunction is then used for quasiparticle correction in the GW approximation 56 , and (4) the dielectric tensor components including excitonic effects are calculated by solving the Bethe-Salpeter equation (BSE) 57,58 . The quality of the wavefunction obtained in step (2) has critical impact onto the quality of the resulting band structure. The first simple test of that quality is whether the chosen hybrid functional satisfactorily reproduces the experimental band gap. A recent comprehensive paper by Furthmüller and Bechstedt 10 as well as papers by other authors 9,10,44-46 use hybrid functionals of Heyd, Scuseria, and Ernzerhof (HSE) 59 for β-Ga 2 O 3 . However, we observe that the HSE band gap energy converges to a value of about 4 eV, significantly less than the experimental value. We tested all popular hybrid functionals: B3LYP 60,61 (published results for β-Ga 2 O 3 using B3LYP functional in Ref. 54), PBE0 62 , HSE, and Gau-PBE 63,64 and found that all of them except HSE produce large band values around 4.7 eV while only HSE produces band gap value around 4 eV. Varley et al. 45 and Peelaers and Van de Walle 51 addressed this problem by increasing the fraction of direct exchange in HSE functional from 25% to 35% to reproduce the experimental band gap of β-Ga 2 O 3 . However, a most recent study by Deák et al. 65 shows that increasing the amount of HF exchange compensates the effect of the screening parameter, which is the reason of the failure of the HSE functional in the case of β-Ga 2 O 3 . Removing the screening, on the other hand reduces the HSE functional to the PBE0 functional. In this work, we use Gau-PBE which to the best of our knowledge has not yet been used for band structure calculation of β-Ga 2 O 3 . The primary purpose of our DFT calculations is to identify band-to-band transitions and to calculate parameters of CP transitions, and compare these with contributions to the experimental dielectric function of β-Ga 2 O 3 . For this purpose, the band structure calculations at the hybrid HF-DFT level are sufficient. 65,66 Hence, we only perform steps (1) and (2) above in this present work.
Sturm et al. 9 recently investigated the NIR-UV dielectric function tensor elements of β-Ga 2 O 3 and noted the same observation made earlier by Jellison et al., where a fourth, independent tensor element was needed to correctly model calculate measured GSE data from multiple samples and from multiple sample azimuth orientations. Following up, Sturm et al. 8 then provided the first CP analysis of monoclinic dielectric function spectra by adapting the concept of dielectric eigen displacement polarizations in the Born and Huang model towards electronic band-to-band transitions and associated exciton contributions for β-Ga 2 O 3 . In their analysis, Sturm et al. described a dipole analysis of the dielectric function by including Tanguy's model into their direction dependent MDF approach and by assuming a band-to-band transition independent exciton binding energy parameter. In this present work, we also adapt the Born and Huang model which we have described for use in the Far-IR-IR spectral region for accurate analysis of GSE data of β-Ga 2 O 3 prior to Sturm et al.'s work 23 . In the eigen polarization approach, transitions are associated with a displacement direction dependent transition amplitude in analogy to the quantum mechanical concept where the transition matrix element is proportional to the transition probability depending on the choice of coordinates. The coordinates describe a given electronic system relative to the polarization direction of the incident photon. Instead of using Tanguy's expression for bound and unbound Wannier exciton contributions, we use functions described by Higginbotham, Cardona and Pollak (HCP) 67 for contributions from 3-dimensional Van-Hoofe singularities, and we augment asymmetrically broadened Lorentz oscillators to account for excitonic contributions. As a result, we find direction dependent (band-to-band transition dependent) exciton binding energies as well as band-to-band transition energies for β-Ga 2 O 3 which differ in detail from those reported by Sturm et al., and we provide additional information on band-to-band transitions into the VUV range not previously reported. We compare our observations with effective mass parameters determined in our DFT analysis for the topmost valence and lowest conduction bands, and we discuss recent findings reported in the literature.

A. Mueller matrix generalized ellipsometry
Ellipsometry is an indirect measurement technique that utilizes changes in the polarization state of light transmitted through or reflected off a sample surface. In general, optically anisotropic materials necessitate the application of generalized ellipsometry 13 . We have previously provided detailed discussion on a variety of single crystalline materials and thin films with uniaxial and biaxial anisotropy using this technique 13,18,19,[68][69][70] . Here, we investigate off-axis cut surfaces of β-Ga 2 O 3 single crystals. Multiple samples cut at different angles from the same crystal are investigated using Mueller matrix generalized ellipsometry (MMGE) at multiple angles of incidence and multiple sample azimuthal angles, and all data are then analyzed simultaneously.
The Mueller matrix connects incident and emergent real-valued Stokes vector components by: where the Stokes vectors are described as: with I p , I s , I +45 , I −45 , I σ+ , and I σ− refer to intensities with p, s, +45 • , -45 • , right-hand, and left-hand polarization, respectively 12 .
In order to extract physical parameters from a GSE measurement (MMGE data), appropriate physical boundary conditions must be applied. For this investigation we use the so-called substrate-ambient approximation, where the single crystalline β-Ga 2 O 3 samples correspond to the half-infinite substrate 12 . For the substrate, the only unknown physical parameters are the dielectric function tensor elements 71 . For monoclinic materials, in order to obtain physically meaningful dielectric function tensor elements, it is necessary to assign coordinate relations between laboratory coordinate axes (x,ŷ,ẑ) and crystallographic axes (a, b, c). 23 We choose theẑ axis of the laboratory coordinate system to be normal to the sample surface, thereby defining the sample surface as the laboratoryx-ŷ plane. By our choice, the (x,y,z ) system is described in Fig.1 with respect to the crystal structure. Euler angles (φ, θ, and ψ) are then determined to describe angular rotations necessary to relate (x,y,z ) with (x,ŷ,ẑ). Effects of finite roughness always present on the nanoscale when measuring polished crystal surfaces must be accounted for. An effective medium approximation (EMA) approach is commonly used to mimic the effect of a very thin effective layer with thickness much smaller than all wavelengths for data analysis 72 . Rigorous treatment of the combination of roughness and anisotropy has not been investigated yet, hence, an isotropic averaging approach was employed here. Thus, our roughness layer model was calculated by averaging all four dielectric tensor elements and then added together in the EMA approach assuming 50% void.

B. Monoclinic dielectric tensor description
With the coordinate choices in Fig. 1 we select the dielectric tensor cross-term element ε xy as the fourth independent complex dielectric function tensor element to describe the monoclinic properties in the corresponding a-c plane, and the remaining elements we set to zero

C. Eigen polarizability critical point model
We adopt the concept of the Born and Huang model and consider electronic contributions to the dielectric response of monoclinic β-Ga 2 O 3 as the result of eigen dielectric displacement processes. Each individual contribution l is characterized by a CP model function, l (ω) and its eigen dielectric polarizability unit vector,ê l . The same approach was adopted by Schubert et al. 23 and Sturm et al. 8 for analysis of GSE data for FIR-IR and NIR-VUV spectral regions, respectively, as follows: When l =0, the quasi-static (ω → 0) dielectric tensor dyadicε DC is determined. For monoclinic materials the tensor can be described according to: with α j equal to the angle of the shear projection into the a-c plane and m,n equal to the number of CP contributions in the a-c plane and b direction, respectively. a. Fundamental band-to-band transitions: We use photon energy (hω) dependent functions described by Higginbotham, Cardona and Pollak 67 for rendering electronic contributions at 3-dimensional Van-Hoofe singularities ("M 0 "-type CP in Adachi's CP composite approach 28 ): with χ = (hω + iB)/E, and A, E, and B are, respectively, amplitude, transition energy, and broadening parameters, and i 2 = −1. Our choice is directed by inspection of the symmetry and band curvatures for the lowest band-to-band transitions observed in our DFT calculations. b. Exciton contributions at fundamental band-to-band transitions: The contributions to the dielectric function due to exciton absorption arise from two parts, one from bound states and another from continuum states 11,26,27 . For Wannier-type excitons, Tanguy developed model functions for parabolic bands taking into account bound and unbound states [29][30][31] . These functions, strictly valid for parabolic bands and isotropic materials only, were used by Sturm et al. 8 for analysis of GSE data from β-Ga 2 O 3 . In their work, Sturm et al. did not determine the exciton binding energy parameter from using the Tangyu model approach. In our present work, and because the dominant contribution to exciton absorption processes in direct-band gap semiconductors is the ground state (n=1) transition 11,26 , we employ a single Lorentz oscillator with non-symmetric broadening to account for, and to spectrally locate the ground state excitonic contribution, thereby further following Adachi's CP composite approach 28 with A, E, B, and b are, respectively, amplitude, energy, broadening, and asymmetric broadening parameter, respectively. The exciton binding energy parameter, R [eV] can then be approximated by where E t and E x are the band-to-band transition energy and the energy of the respective exciton CP contribution, respectively, ε −1 ∞ jj is the element jj of the inverse dielectric function tensor at photon energies above the phonon mode spectral region, and µ −1 jj is the element jj of the inverse effective reduced mass tensor, for any given unit vector directionĵ. The inverse effective reduced mass tensor element is obtained from summation over the inverse effective mass tensor element of conduction and valence bands participating in the transition: vb jj , i.e., electron and hole effective mass parameters along polarization directionĵ, respectively.
c. Above-band gap band-to-band transitions: At photon energies far above the band gap, multiple transitions originating at multiple points in the Brillouin zone often overlap, and CP features due to individual transitions cannot be differentiated by experiment. Hence, broadened Lorentzian or Gaussian oscillators are often used to account for broad CP features typical for aboveband gap spectra. Here, we use the same anharmonic broadened functions as in Eq. 7.
d. Higher energy band-to-band transitions: Transitions above the spectral range investigated here contribute to the overall lineshape of the dielectric functions at wavelengths within the investigated spectral range. Such higher energy contributions are usually accounted for by Gaussian broadened oscillator functions: where amplitude A, center energy E, and broadening B are adjustable parameters. The real part, ε 1 , is obtained by the Kramers-Kronig integration 73 : Note that each non-trivial sum in Eqs. 5 satisfies the Kramers-Kronig integral condition 23,43,74 , and which can be set as additional side condition during the CP-MDF analysis.

D. Density Functional Theory
Theoretical calculations were performed using plane wave density functional theory (DFT) code Quantum ESPRESSO (QE) 75 . The primitive cell of β-Ga 2 O 3 consisting of six oxygen and four gallium atoms was used, and the initial atomic positions and parameters of the unit cell were taken from Ref. 76. The atoms were represented by norm-conserving pseudopotentials from the QE library. Following our previous calculation of phonon modes in β-Ga 2 O 3 23 , we do not include the gallium 3d electrons in the valence configuration. Wright and Nelson 77 in their computational studies for GaN observed a significant overlap and hybridization of gallium 3d states with 2s energy levels of nitrogen, and suggested explicit treatment of these shallow states for gallium compounds. However, band structures of β-Ga 2 O 3 published to date (for example Refs. 10 and 53) clearly show that the gallium 3d levels are positioned not only far (>10 eV) below the valence band maximum, but also well separated from other states, and little hybridization can be observed within the energy window considered in the current study. Therefore we opted for pseudopotentials for gallium that leave the 3d electrons in the core.
First, the structure was carefully relaxed to force levels less than 10 −5 Ry/bohr using the exchange-correlation functional of Perdew, Burke and Ernzerhof (PBE). 78 A 4×4×4 regular shifted Monkhorst-Pack grid was used for sampling of the Brillouin zone 79 . A convergence threshold of 10 −12 was used to reach self-consistency with a large electronic wavefunction cutoff of 100 Ry. The resulting lattice parameters obtained are shown in Tab. I in comparison to results from previously reported studies using GGA-DFT methods. We find very good agreement between our values and those reported previously. Second, the structure fully relaxed at the PBE level was used for electronic structure calculations employing the hybrid Gau-PBE functional. This calculation was performed with a 6×6×6 Γ-centered Monkhorst-Pack grid (and after testing the convergence with respect to the grid of k-points up to 8×8×8), and with otherwise the same parameters as for the preceding PBE calculations. The converged Gau-PBE wavefunction was used to analyze the band structure.
The visualization program XCrySDen 82 produces the proper Brillouin zones for all known sets of lattice vectors, including base-centered monoclinic, using the Voronoi decomposition. The orientation of the Brillouin zones for base center monoclinic lattices, however, depends on the definition of the primitive cell used. While commonly the symmetry axes is assigned with 'unique axis b' and the monoclinic angle β > 90 • , Setyawan and Curtarolo 50 used an unconventional definition using 'unique axis a' and monoclinic angle α < 90 • . In the standard convention of using 'unique axis b' and the monoclinic angle β > 90 • , a primitive cell with vectors p 1 = (a + b)/2 and p 2 = (−a + b)/2 produces a Brillouin zone like in Fig. 2 in Ref. 51, whereas a primitive cell with vectors p 1 = (a−b)/2 and p 2 = (a+b)/2 produces a Brillouin zone identical with Fig. 17 in Ref. 50 if β < 90 • , and one other when flipped vertically for β > 90 • , shown in Fig. 2(a) and Fig. 2(b), respectively. In our present study, we use the structure shown in Fig. 2(b). If high symmetry points are correctly identified, then physical properties descending from such high symmetry points should not depend on the choice of coordinates or convention for the Brillouin zone description. For convenience, we define these points in Fig. 2 for the Brillouin zone oriented exactly like in the paper of Setyawan and Curtarolo. Our coordinates of the high symmetry points in both orientations are given in Tab. II.
The band structure along a high symmetry path was plotted using the band interpolation method based on the maximally localized Wannier functions 83,84 as implemented in the software package WANNIER90 85 . We used s and p orbitals on both Ga and O atoms and performed disentanglement of the bands in a frozen energy window from -5 eV to 22 eV. The disentangled bands were also used for calculating the effective masses of the carriers. The bands were sampled in the range ±0.005Å −1 from the Γ point parallel to the crystal directions a, b, and c. II. Example coordinates of the high symmetry points in the Brillouin zone. Note that one can draw four symmetryequivalent paths through each of the Brillouin zones, i.e., one for each possible irreducible BZ.

Label
Coordinates Coordinates for BZ in Fig. 2(a) for BZ in Fig. 2  Parabolic curves were used to fit the dispersions of the respective energy bands and the quadratic terms of the parabolas were converted to inverse effective mass tensor parameters as follows: where derivatives are taken along directions k = k jĵ with unit vectorĵ, for example, parallel to a, b, or c. Finally, the significant band-to-band transitions contributing to the experimental dielectric tensor measured by ellipsometry were identified by analyzing the matrix elements |M cv | 2 of the momentum operator between conduction and valence bands at the Γ point (Tab. V). The signatures (parallel or anti-parallel) of the projections of |M cv | 2 along the crystal directions a and c * were obtained from inspecting the complex argument of M cv . Transition eigenvectors with parallel arguments (antiparallel) were plotted in the first (second) quadrant of the Cartesian (a − c ) plane.

III. EXPERIMENT
Bulk single crystalline β-Ga 2 O 3 was grown by Tamura Corp., Japan by the edge-defined film fed growth process as described in Refs. 86-88. Samples were then cut at different orientations to the dimensions of 650 µm × 10 mm × 10 mm, and then polished on one side. In this paper we investigate a (010) and a (201) surface.
Mueller matrix generalized spectroscopic ellipsometry data were collected from 133 nm to 1690 nm. The vacuum ultraviolet (VUV) measurements were obtained using a rotating-analyzer ellipsometer with an automated compensator function (VUV-VASE, J.A. Woollam Co., Inc.). Data were acquired at three angles of incidence (Φ a =50 • , 60 • , 70 • ), and at several azimuthal angles by manually rotating the sample about the sample normal in steps of ≈45 • . Note that in the VUV range, due to limitations of the instrument, no elements in the 4 th row of the Mueller matrix are available. Measurements from the near infrared to near ultraviolet (NIR-NUV) were performed using a dual-rotating compensator ellipsometer (RC2, J.A. Woollam Co., Inc.) allowing for the determination of the complete 4×4 Mueller matrix. Measurements were taken at three angles of incidence (Φ a =50 • , 60 • , 70 • ), and at different orientations by auto-rotating the sample by steps of 15 • beginning at the same azimuthal orientation as in the VUV measurements. All model calculations were conducted using WVASE32 TM (J. A. Woollam Co., Inc.). Each data set (sample, azimuthal orientation, angle of incidence) is unique, however, characteristic features are shared between them all at energies indicated by vertical lines. While we do not show all data in Figs. 3 and 4, we note that all data sets are identical when samples are measured at 180 • rotated azimuth orientation. Most important to note in the experimental Mueller matrix data is the clear anisotropy shown by the nonzero off-diagonal block elements (M 13 , M 14 , M 23 , M 24 ) and strong dependence on sample azimuthal orientation in all Mueller matrix elements. All data gathered by the measurement of multiple samples, with multiple orientations, and at multiple angles of incidence were analyzed simultaneously using a best-match model data regression procedure (polyfit). For each energy, up to 144 indepen-dent data points were included from 2 samples, 3 angles of incidence, and as many as 24 different azimuthal orientations. Only 8 independent model parameters for real and imaginary parts of ε xx , ε yy , ε zz , ε xy as well as 2 sets of energy-independent Euler angles describing the sample orientation and crystallographic structure and 2 roughness layer thickness parameters were fit for. The thickness parameters for the roughness layer of the (010) and (201) samples were determined to be 1.78(1) nm and 1.61(1) nm, respectively. The best match model calculated Mueller matrix elements from the polyfit procedure are shown in Figs. 3 and 4 as red solid lines. We obtain an excellent agreement between model calculated and experimental Mueller matrix data. Euler angle parameters noted in the captions of Figs. 3 and 4 are in agreement with anticipated orientations of the crystallographic axes of each of the samples. The dielectric function tensor elements, ε xx , ε yy , ε xy , and ε zz determined from the wavelength-by-wavelength polyfit procedure are shown in Fig. 5, Fig. 6, Fig. 7 and Fig. 8, respectively, as dotted green lines.

B. CP model analysis
We identify 11 differentiable contributions in ε xx , ε yy , and ε xy , and 5 in ε zz . Distinct features can be seen, e.g., in the imaginary part of each tensor element in Figs. 5-8. Vertical lines are drawn corresponding to the results from our CP analysis at the respective CP transition energy model parameters. Note that while vertical lines are identical for ε xx , ε yy , and ε xy , a different set is seen for ε zz which corresponds to the difference between the monoclinic a-c plane and the axis parallel to b.
a. a-c plane Eleven CP features are needed to match the tensor elements. Functions described in Sect. II C are used to model individual CP contributions as projections in the x-y plane with angular parameters α CP,j . The lowest band-to-band transition, CP ac 0 was modeled with the HCP CP function (Eq. 6) with an excitonic contribution determined by an asymmetrically broadened Lorentzian oscillator (Eq. 7). It was assumed that excitonic and band-to-band transition contributions share the same unit eigenvector (angular parameter α CP,0 ). The exciton CP contribution is labeled CP ac 0x . We identify a second pair of CP contributions (CP ac 1 , CP ac 1x , α CP,1 ) using the same functions. A Gaussian oscillator was used to model a very small CP contribution at 5.64 eV which could not be further differentiated (CP ac 2 , α CP,2 ). Above-band gap CP contributions (CP ac 3−6 , α CP,3−6 ) were identified at higher photon energies, which were modeled by asymmetrically broadened Lorentzian oscillators (Eq. 7). We were unable to differentiate between excitonic and band-to-band transition contributions associated with these higher energy CPs. Contributions due to higher-energy transitions, outside the investigated spectral region, were accounted for by a Gaussian function CP with projection along x (CP ac 7 ), and y (CP ac 8 ). The resulting best-match CP-MDF parameters are listed in Tab. III, and are shown as solid red lines in Figs. 5-7. An excellent agreement between our   tions (CP b 1−2 ) were identified and modeled by functions in Eq. 7. Here again, we were unable to differentiate between excitonic and band-to-band transition contribution. Contributions due to higher-energy transitions, outside the investigated spectral region, were accounted for by a Gaussian function CP with projection along z (CP b 3−4 ). The resulting best-match CP-MDF parameters are listed in Tab. IV, and are shown as solid red lines in Fig. 8. Again, an excellent agreement between our GSE wavelength-by-wavelength obtained and CP-MDF calculated data is noted. We also note close agreement with the GSE wavelength-by-wavelength obtained data reported by Sturm et al.   III. CP-MDF parameters for polarization within the a-c plane of β-Ga2O3 obtained in this work from MMGE wavelength-by-wavelength data analysis of (010) and (201) surfaces of single crystalline bulk β-Ga2O3.   two plots is the expected opening of the energy gap between the valence and conduction bands by about 2 eV. The most important feature, common to both plots, is the lowest conduction band which dominates the Brillouin zone center. A direct comparison between our band structures and previously published data is rather difficult as the majority of the authors only considered some of the high symmetry points in the first Brillouin zone. The only complete band structure to date was published by Peelaers and Van de Walle 51 . They used the HSE density functional with the fraction of HF exchange adjusted to reproduce an assumed value of the band gap. Their valence bands are very similar to ours, whereas their conduction bands are slightly shifted to higher energies and steeper than ours, most likely due to the effect of the higher amount of HF exchange included into the calculations (35% vs 24% in Gau-PBE). Otherwise the agreement is rather good. The character of the band gap in β-Ga 2 O 3 can be obtained from the band structure. The broad valence band maximum (VBM) has been reported previously along the However, the energy difference between the actual VBM and, for example, the top valence band energy at the L-point only amounts to few meV. Thus studying the band properties at high-symmetry points L and at the Γ points is accurate enough. Local density approximation DFT methods typically predict the band gap to be indirect, and render the valence band at the L-point about 100 meV higher than the direct gap at Γ 52 . At the GGA-DFT level this difference is reduced to about 20-50 meV 10,53 , and this usually holds for hybrid HF-DFT as well 10,45,51,54 . Our results at the Gau-PBE level show the VBM near the L point about 50 meV higher than the top valence band at Γ. Interestingly, at the GW level (quasiparticle bands) the band gaps are completely degenerate, or even the direct gap appears marginally higher 10 . However, Ratnaparkhe and Lambrecht 89 used the quasiparticle self-consistent version of GW, QSGW 90 , and obtained the indirect band gap energy smaller by nearly 100 meV than the direct band gap energy.

Band-to-band transitions
We analyze band-to-band transitions by identifying all allowed transitions, i.e., transitions with non-zero matrix elements of the momentum operator between conduction and valence bands, and whose transition energies are less than 10 eV. We find eight such transitions, summarized in Tab. V and Tab. VI, respectively, with polarization within the a-c plane and along axis b, respectively. The transition labels are according to the numbers of the bands involved, where numberings start at the top of the valence and at the bottom of the conduction bands. The matrix elements of the momentum operator are obtained from the overlap of the wavefunctions for the respective energy bands. Hence, their values represent the probabilities of the transitions, i.e., transition amplitudes, which can be compared to experimental ones. In case of a-c plane transitions, which do not have any pre-defined orientation within the plane, the transition probabilities along the crystallographic directions a and c * constitute Cartesian components of the corresponding transition vectors, thereby defining the polarization direction of the band-to-band transition in space. Their orientations are shown in Fig. 12(b), and can be compared to the eigen dielectric displacement unit vectors obtained from GSE analysis. The fundamental (lowest energy) band-to-band transition is polarized nearly parallel to the crystallographic axis c. It is closely followed by a second transition polarized at a small angle from the crystallographic axis a. The lowest transition along the crystallographic axis b occurs about 0.6 eV above the fundamental transition in the a-c plane. This tendency agrees well with GW results shown in Tab. VIII of Ref. 10.  and within the a-c plane (black solid arrows). Note that the reciprocal space direction Γ-X corresponds to the real space directions parallel to the crystallographic vector b, and the reciprocal space direction Γ-Y corresponds to the real space direction lying within the a-c plane, inclined at a small angle from the crystallographic vector a. Fig. 10 shows the vicinity of the Γ point of the Brillouin zone, with top valence bands and the first conduction band. The four lowest transitions are schematically shown as vertical arrows. The first conduction band is clearly parabolic, and rather symmetric, indicating a nearly isotropic electron effective mass, which is consistent with many previous studies 10, 52 . It has been assumed previously 8 that due to the valence bands being generally flat, the hole effective masses are expected to be large, and that the electron effective mass parameter hence dominates the carrier reduced masses for the zone center band-to-band transitions. As can be inferred from Fig. 10, however, the valence band structure is far from being isotropic. For example, the curvature of the second valence band in the direction Γ − Y , and the curvature of the fourth valence band in the direction X −Γ are similar to the curvature of the first conduction band, and reveal strong anisotropy as well.

Conduction and valence band effective mass parameters
a. Electron effective mass: The electron effective mass for β-Ga 2 O 3 has been studied previously, both by computation and experiment. Computational results consistently predict a very small anisotropy, but span a relatively wide range of values: from (0.12...0.13)m e (GGA-DFT) 53  with the HSE functional, and 0.34m e for the B3LYP functional 54 . Our results (Gau-PBE) are presented in Tab. VII, which fall within this broad range of reported values, but exhibit a slightly higher anisotropy than found in previous studies.
b. Hole effective mass parameters: We have analyzed the effective mass parameters for the three valence bands involved in the lowest band-to-band transitions, and data are presented in Tab. VII. It is obvious that the hole effective mass anisotropy cannot be neglected for these bands. To our best knowledge, hole effective mass parameters at the Γ point have not been reported for β-Ga 2 O 3 thus far. Yamaguchi 52 presented values of the top valence band effective mass parameter at point labeled "E" away from the zone center and thus not relevant for zone center transitions.
We note an interesting observation from our analysis here: The lowest values of the hole effective mass for each valence band occurs in the approximate polarization direction of the transition that connects this particular valence band and the lowest conduction band, and which we observe and identify from our GSE and DFT analyses. For the first and topmost valence band, the lowest value of the effective mass occurs along axis c, and the transition Γ 1−1 is polarized nearly along axis c as well. For the second valence band, the lowest effective mass is along axis a and the transition Γ 1−2 is polarized near axis a. For the fourth valence, band the lowest hole effective mass is along b and the transition Γ 1−4 is polarized along b. We thus observe here a clear correlation between the transition selection rules for electronic band-to-band transitions and the values of the carrier effective masses for these transitions. In contrast to previous studies, we find that not only do the hole effective mass parameters matter, but due to their very large anisotropy these parameters may play a decisive role for the polarization of the band-to-band transitions and the associated exciton CP contribution transition parameters.

Exciton binding energy parameters
Table VII also lists the exciton binding energies for the three lowest band-to-band transitions, for which the excitonic CP model parameters can be resolved in our experimental data. Exciton binding energy parameters were calculated according to Eq. 8 from reduced mass parameters of each electron and hole pair for the lowest three transitions. Because of the large values of R observed in our experiment, which exceed the energy values of all infrared-active phonon modes 23 , we chose the static dielectric function parameters from experimental GSE wavelength-by-wavelength determined epsilon tensor element values at photon energy well above the reststrahlen range, at 0.75 eV (Tab. VII). We observe a good agreement between calculated and experimental values for R . Values for R are found different for each transition in contrast to those calculated by Sturm et al.   (E=9.714 eV). The CP contributions may contain further transitions from different Brillouin zone regions not further investigated. In Fig. 11 we also indicate a small contribution (CP ac 2 at 6.53 eV) for which we do not observe an equivalent transition in our DFT results. This energy is close to a strong contribution identified along axis b, and it appearance in the a-c plane may originate from lattice defects or from slight experimental misalignment. Figure 12 (a), (b) depict CP transition eigenvectors multiplied with their respective CP transition amplitude parameters, or transition matrix element, obtained from GSE and DFT, respectively. Colors and linestyles are as in Fig. 11, for convenient guiding of the eye. As one can see, the agreement between eigen dielectric displacement unit vectors in our CP-MDF approach and the polarization selection conditions obtained from DFT is re-TABLE VII. Effective mass parameters for conduction (c) and valence (v) bands as indexed, and reduced effective mass parameter for lowest transition along directionsĵ in units of me. The exciton binding energy R calc is then obtained for the lowest transition according to Eq. 8. R exp is obtained from GSE analysis in this work. Data for ε∞,jj were obtained from GSE data measured in this work at 0.75 eV.  markably good, in particular for the first two band-toband transitions. None of the identified contributions are purely polarized along either axis a, c, or c*. At higher energies we see a considerable shift between GSE and DFT. We attribute this to an increase in error associated with both the experimental results and the calculations at higher energies. The higher energy transitions predicted by DFT calculations which cannot be resolved from our GSE investigation are shown all in dark blue, the remaining colors correspond to the associated transitions identified by our GSE analysis. Previous work assumed transitions were independently polarized along crystallographic axes. We find here that the lowest two transitions are indeed polarized close to crystal axes c and a, respectively. We find that the energy levels of the lowest exciton contributions agree very well with those found in our work. However, because Sturm et al. imposed the constraint of fixed and isotropic exciton binding energy parameters, we find in detail different band-to-band transition energy parameters is our work. Table VIII summarizes Table VII lists exciton binding energy parameters for the 2 lowest transitions within the a-c plane and the lowest transition for axis b. GSE data did not provide sufficient sensitivity to differentiate exciton contributions from higher energy transitions. Using the effective mass parameters obtained from our DFT calculations and DFT calculated transition energies, we derived binding energies for Wannier-type excitons, and assuming absence of phonon interaction. We observe that both data sets reflect direction dependent binding energy parameters, and both sets show similar magnitude. Our calculated data predict the largest [smallest] binding energy parallel axis b (0.205 eV) [a (0.171 eV)], while our GSE results results indicate axis a (0.23 eV) [c (0.12 eV)]. The deviations may originate from the simplification made by assuming an isotropic exciton picture separately for each transition. Nonetheless, our results clearly show that the exciton contributions should not be assumed direction independent, and further theoretical and experimental analysis for semiconductor materials with monoclinic crystal symmetry may be needed. Our finding further suggests that the assumed binding energies of 270 meV by Sturm et al. overestimated the exciton shift of the onset of absorption in β-Ga 2 O 3 .

V. CONCLUSIONS
The eigen dielectric displacement model, previously described for analysis of long wavelength generalized ellipsometry data as well as for near-infrared to ultra violet spectral regions for for β-Ga 2 O 3 , was applied for an analysis expanding into the vacuum ultra violet spectral region. In our analysis we permitted for direction dependent exciton binding parameters. We differentiated 9 critical point contributions in the a-c plane, the lowest 2 of which were modeled with excitonic contributions. Additionally, we observed five critical point contributions in the b direction, with an excitonic contribution associated with the lowest transition. The binding energy parameter of these excitons associated with the lowest band-toband transitions along each of the experimental coordinate system axes were shown to present with distinct values ranging from 0.12-0.23 eV. Additionally, transitions in the monoclinic plane, which does not contain any symmetry operation, were found to be distributed within the plane and none of them would align with major crystal directions a or c. Our experimental analysis compares well with results from density functional theory calculations performed using a recently proposed Gaussianattenuation-Perdue-Burke-Ernzerhof density functional, We presented and discussed the order of the fundamental direct band-to-band transitions and their polarization selection rules, the electron and hole effective mass parameters for the three lowest band-to-band transitions, and their exciton binding energy parameters, in excellent agreement with our experimental results. We find that the effective masses for holes are highly anisotropic and correlate with the selection rules for the fundamental band-to-band transitions, where the observed transitions are polarized closely in the direction of the lowest hole effective mass for the valence band participating in the transition. The MDF approach and parameter set for β-Ga 2 O 3 presented here will become useful for ellipsometry analysis of heterostructures, and may be expanded for description of alloys with monoclinic crystal symmetry.