Engineering multifunctionality in gallium oxide: unveiling novel structural, electronic, and opto-mechanical attributes of gadolinium-doped β-Ga₂O₃ through advanced first-principles design

Abstract

The design of next-generation wide band gap semiconductors requires both intrinsic optimization and judicious doping strategies. In this work, we comprehensively examine pristine and gadolinium-doped β-Ga₂O₃ (Gd–Ga₂O₃) through a multi-scale computational approach combining HSE06 hybrid DFT, GGA + U and PBESol functionals, MLFF-accelerated ab initio molecular dynamics, and CALYPSO global structure prediction. Gd preferentially substitutes tetrahedral Ga sites in monoclinic β-Ga₂O₃, inducing local lattice distortions but preserving high thermal stability up to 700 K. Mechanical analysis reveals lattice softening without loss of ductility, while electronic calculations show Gd-induced impurity states that enable visible-light absorption and modest band-gap narrowing. Beyond the monoclinic phase, CALYPSO predicts a novel thermodynamically stable triclinic Gd–Ga₂O₃ structure where Gd adopts a distinctive icosahedral coordination. This phase exhibits a remarkably wide direct band gap of 7.0 eV alongside a ∼60% enhancement in visible-light absorption compared to the monoclinic counterpart. Furthermore, it combines high ductility with a greatly increased bulk modulus, reflecting enhanced incompressibility and mechanical robustness. Overall, these results demonstrate that Gd doping substantially tailors the structural, electronic, optical, and mechanical properties of Ga₂O₃. The discovery of a stable triclinic phase with a unique balance of deep-UV transparency, visible-light activity, and superior mechanical strength positions Gd–Ga₂O₃ as a promising multifunctional material for advanced optoelectronic and high-performance device applications.

---

1. Introduction

Because of its large bandgap (4.5–4.9 eV), high breakdown electric field strength, and excellent transparency in the ultraviolet range, monoclinic β-gallium oxide (β-Ga₂O₃) has emerged as a promising ultra-wide bandgap (UWBG) semiconductor with a wide range of applications in power electronics, deep-ultraviolet photodetectors, and high-voltage switching devices.^1–5 Its development for industrial applications has been driven by its low-cost availability and compatibility with a variety of single crystalline growth techniques, such as vertical gradient freeze (VGF), edge-defined film-fed growth (EFG), and Czochralski (CZ).^6–8

The doping chemistry and defect of β-Ga₂O₃ have a significant impact on its functioning. The search for efficient extrinsic doping techniques is still ongoing, despite the fact that native oxygen and gallium vacancies have been connected to deep donor and acceptor states that significantly affect electrical and optical properties.^9–11 The self-trapping of holes, their high effective mass, and the predicted formation of small hole polarons^12,13 make it difficult to achieve stable p-type conductivity, even though n-type conductivity can be easily obtained with shallow donor dopants such as Sn, Si, and Ge.^14–16

Different dopants, such as transition metals (Ti, Cr, Mn, Fe),^17–20 nonmetals (N),²¹ and rare-earth (RE) elements (Eu, Tb, Tm, Er, Nd),^22–25 have been studied in order to overcome the limitations caused by intrinsic defects and to fine-tune the optoelectronic properties of β-Ga₂O₃. In addition to their partially filled f-orbitals, which provide magnetic and optical tunability, rare-earth dopants like Gd are particularly interesting because of their potential for the use in radiation detection applications and as potential acceptors in the hunt for p-type doping.^23,24 Since both share a stable +3 oxidation state, substitutional Gd³⁺ at Ga³⁺ sites in β-Ga₂O₃ is chemically advantageous even though Gd has a much higher ionic radius than Ga.²⁶ The monoclinic lattice can tolerate the induced lattice deformation from Gd incorporation, but the strong electronic correlation from localized f-states requires precise theoretical treatments that go beyond the standard local or semi-local density functional theory (DFT).

Band gaps are often underestimated by standard approximations like the Perdew–Burke–Ernzerhof (PBE) functional because of self-interaction mistakes. On the one hand, this restriction can be somewhat overcome and better handling of localized electrons, such as the f-electrons in rare-earth dopants, can be provided by adding a Hubbard U correction (DFT + U). Hybrid functionals such as Heyd–Scuseria–Ernzerhof (HSE06), on the other hand, describes electronic structures and formation energetics in strongly correlated systems more accurately by incorporating some exact exchange.^27,28 HSE06 has become essential for accurately forecasting defect formation energies and structural stability in doped wide-bandgap semiconductors, despite its computational cost.²⁹

Few studies have used the HSE06 method to assess dopant effects in β-Ga₂O₃, but many have already used the GGA and GGA + U approaches. For instance, the effects of dopant atomic number and the selection of U parameter on the electrical structure and bandgap alteration have been emphasized by systematic studies of transition metal doping (Sc, Ti, V, Cr, Mn, Fe, Co, Ni, Cu, and Zn) utilizing GGA and GGA + U.^30–35 To the best of our knowledge, Gd doping in β-Ga₂O₃ at the HSE06 level has only been studied once before.²⁷ In that work, the optoelectronic characteristics of β-(Ga_1−xGd_x)₂O₃ alloys with Gd concentrations up to 50% were evaluated using hybrid density functional theory (HDFT). Additionally, the work suggested a feasible strategy for modeling f-orbital systems by combining hybrid functional exchange-correlation techniques with pseudopotentials. The authors assessed Gd substitution at both locations and determined that the octahedral site was energetically preferable in the monoclinic β-phase, where Ga atoms occupy both tetrahedrally and octahedrally coordinated sites.

Here, we performed a thorough first-principles study of Gd doping in β-Ga₂O₃ (Gd–Ga₂O₃), focusing on the consequences of Gd substitution at both tetrahedral and octahedral Ga sites. Our objective is to resolve persistent concerns about the thermodynamic site preference of rare-earth dopants and to comprehend the ensuing changes in structural, electrical, and optical characteristics. The HSE06 hybrid functional, which provides improved accuracy for systems with localized f-electrons, was used for structural optimization. The PBE + U method was used to further investigate the electronic and optical properties in order to attain computational efficiency without compromising important information.

The correct handling of Gd f-electrons, which are highly correlated and have a substantial impact on the defect energetics and electronic structure, is one of the primary issues in this field. Our study demonstrates the site-dependent behavior of Gd in β-Ga₂O₃ and shows how rare-earth doping can produce new optical transitions that are not present in the undoped system by methodically examining both coordination environments.

By using this integrated approach, which combines economical electronic structure analysis with precise hybrid-functional structural modeling, we help to resolve discrepancies in the literature and offer useful theoretical direction for the creation of RE-doped β-GaO₃ materials. The development of next-generation optoelectronic and transparent electronic devices, where accurate doping control is still a major design issue, is pertinent to these findings.

2. Computational methods

The structural and energetic features of pristine and doped β-Ga₂O₃ systems were investigated by first-principles calculations based on DFT utilizing the Vienna Ab Initio Simulation Package (VASP) version 6.3.2.^36–39 The projector augmented wave (PAW) approach was used to describe the interaction between valence and core electrons. The HSE06 hybrid functional was employed for structural optimization and formation energy calculations in order to directly capture the localized electronic interactions associated with the Gd f-orbitals.⁴⁰ Unlike a sequential GGA pre-relaxation, which can bias the equilibrium geometry through an incomplete treatment of strong correlations, direct optimization at the hybrid-functional level ensures a more accurate description of lattice distortions, site preference, and exchange–correlation effects.

As shown in Fig. 1, the monoclinic unit cell of β-Ga₂O₃, which consists of eight gallium and twelve oxygen atoms, served as the basis for the structural model. The initial structural parameters were taken from the experimental monoclinic β-Ga₂O₃ structure reported by Geller.⁴¹ In order to minimize interatomic forces below 10 meV Å⁻¹, geometry optimizations were carried out at constant volume, allowing atomic locations and cell shape to relax. Only atomic positions were further optimized while maintaining the fixed cell shape in a second static calculation that was performed to refine the total energy. With a constant plane-wave energy cutoff of 520 eV, the electronic configurations regarded as valence were Ga (3d¹⁰4s²4p¹), O (2s²2p⁴), and Gd (4f⁷5d¹6s²). In all self-consistent calculations, a 2 × 8 × 4 Gamma-centered Monkhorst–Pack k-point mesh⁴² was used to sample the Brillouin zone.

Fig. 1 First-principles calculated structural and electronic properties of pristine and Gd-doped β-Ga₂O₃. (a) The optimized monoclinic unit cell of pristine β-Ga₂O₃. (b) The optimized structure of Gd-doped β-Ga₂O₃ (with Gd at the tetrahedral site), both relaxed using the HSE06 hybrid functional within VASP. (c) Electronic band structure of pristine β-Ga₂O₃ and (d) Gd-doped β-Ga₂O₃, computed at the PBE + U level. (e) Partial Density of States (PDOS) for pristine β-Ga₂O₃ and (f) Gd-doped β-Ga₂O₃, also calculated using the PBE + U functional in CASTEP.

Using the formula E_f = E_Gd–Ga2O3 – E_Ga2O3 – μGd + μGa (1), the thermodynamic viability of Gd substitution in β-Ga₂O₃ was evaluated. Here, E_f, E_Gd–Ga2O3, and E_Ga2O3 stand for the formation energy, the total energy of the Gd-doped system, and the total energy of pristine β-Ga₂O₃, respectively. The chemical potentials of the Gd and Ga atoms are denoted by the terms μGd and μGa, respectively. Additional calculations were conducted using the Perdew–Burke–Ernzerhof functional for solids (PBESol), which has been shown to provide improved accuracy for lattice constants, elastic properties, and mechanical stability of crystalline solids compared to the standard PBE functional.^43–45 In the case of Ga₂O₃ and related oxides, PBESol reproduces experimental elastic constants and bulk moduli more reliably, whereas PBE tends to overestimate lattice parameters and underestimate stiffness. Therefore, PBESol was selected here to assess the mechanical stability and elastic characteristics of pristine and Gd-doped Ga₂O₃, ensuring quantitative consistency with experimental benchmarks. The VASP framework was also used to run these simulations. The PBESol functional within VASP was used to evaluate mechanical properties for the structures that had previously been optimized using HSE06 and the 2 × 8 × 4 k-point mesh.

Although explicit phonon calculations were not performed, the stability of Gd–β-Ga₂O₃ structures was evaluated using complementary approaches. Ab initio molecular dynamics (AIMD) simulations were greatly accelerated by machine-learned force fields (MLFFs), allowing for the investigation of thermal stability over long timescales. The VASP code was used for all MLFFs-accelerated AIMD computations. The core-valence electron interactions were described using PAW pseudopotentials. The HSE06 hybrid functional was used to perform the structural optimizations that were used as the training data for the MLFFs. A plane-wave energy cutoff of 500 eV was employed, guaranteeing electronic convergence for SCF iterations within 10⁻⁷ eV. A Γ-centered k-point mesh was used for structural relaxations in the MLFF training procedure. The on-the-fly learning features included into VASP were used to create the MLFFs. With this method, the force field may be iteratively refined as the simulation goes on, guaranteeing great accuracy while reducing the amount of computationally costly ab initio stages. The reliability of the machine-learned force fields (MLFFs) used in the AIMD simulations was systematically verified against ab initio reference data. The Bayesian error estimation module in VASP was employed to monitor model uncertainty; configurations exceeding the adaptive threshold were recalculated at the DFT (HSE06) level and added to the training dataset until convergence was achieved. After iterative retraining, the Bayesian error estimations stabilized at ≈3 × 10⁻³ eV atom⁻¹ for energies and ≈0.13 eV Å⁻¹ for atomic forces. The stress error remained below 3 kB. This iterative refinement ensured that the MLFF maintained high fidelity to DFT throughout the production runs. As a result, the MLFF-enabled AIMD trajectories extended well beyond the time and length scales typically accessible to direct DFT, allowing thousands of MD steps and structural optimizations to be sampled efficiently without compromising accuracy. The canonical (NVT) ensemble was used for AIMD simulations, with the number of atoms (N), volume (V), and temperature (T) held constant. A Nose–Hoover thermostat was used to regulate the temperature. For MD simulations, a timestep of 0.01 fs was employed. In order to guarantee numerical stability and precise integration of the equations of motion, this timestep was used.

The CASTEP software⁴⁶ was used to calculate optical properties and electronic structures, taking advantage of its efficient linear-response routines for dielectric functions. To ensure full internal consistency with CASTEP's ultrasoft pseudopotentials, all structures were reoptimized using the Perdew–Burke–Ernzerhof (PBE) functional within the generalized gradient approximation (GGA).⁴⁷ The CASTEP–PBE geometries deviated by less than 1% from the VASP–HSE06 optimized structures, confirming that no significant structural bias was introduced. The valence configurations were O (2s²2p⁴), Gd (4f⁷5d¹6s²), and Ga (3d¹⁰4s²4p¹). The GGA + U approach^48,49 was further applied to mitigate band gap underestimation and to appropriately treat the localized d- and f-electrons. The Hubbard U correction within the GGA + U framework was employed specifically to treat localized states rather than as a generic band gap correction. The optimal values were determined as 2.0 eV for Ga 3d, 7.5 eV for O 2p, and 6.0 eV for Gd 4f, consistent with CASTEP's default settings for rare-earth elements. Although less common, the U applied to O-2p states is justified because the valence band maximum of β-Ga₂O₃ is dominated by O-2p orbitals, which tend to be excessively delocalized in GGA. This over-delocalization leads to an underestimated band gap and inaccurate Ga–O hybridization. As shown in prior reports,^43–45 introducing a moderate U (∼7–8 eV) on O-2p significantly improves the reproduction of the experimental electronic density of states and optical band gap. By calibrating these values against experimental benchmarks, we ensured that the chosen U parameters provide a physically meaningful and accurate electronic description. It should be noted that applying U on top of a hybrid functional (HSE06 + U) would be inconsistent, as HSE06 already mitigates self-interaction and partially corrects d/f electron localization via exact exchange. Therefore, U corrections were only combined with GGA, while HSE06 was used independently for structural optimization and reference-level energetics. To preserve clarity and prevent repetition in the computational methods section, a thorough explanation and rationale of these U values—including a comparison with experimental findings—are given in the results and discussion section. With Monkhorst–Pack k-point grids of 2 × 8 × 4 for structural optimizations and 5 × 20 × 11 for electronic structure and optical property calculations, a plane-wave cutoff energy of 600 eV was utilized. The Broyden–Fletcher–Goldfarb–Shanno (BFGS) optimization technique⁵⁰ was used to perform structural relaxations, and the convergence criteria were set at 10 meV Å⁻¹ for forces and 1 × 10⁻⁶ eV per atom for energy. The difference in plane-wave cutoff energies (520 eV in VASP vs. 600 eV in CASTEP) stems from the use of distinct pseudopotential libraries: in VASP, the PAW pseudopotentials are well converged at ∼500 eV (<1 meV/atom), while CASTEP's ultrasoft pseudopotentials require a higher cutoff (∼600 eV) to ensure reliable convergence of dielectric response calculations.^39,46 The choice of k-point meshes in all calculations was guided by convergence testing (convergence tests confirmed that the total energies calculated with 2 × 4 × 3 and 2 × 8 × 4 differed by only ∼0.03 eV) and validated against available experimental benchmarks, with the 2 × 8 × 4 grid ensuring reliable structural parameters at the HSE06 level, the PBESol-based meshes yielding mechanical properties in close agreement with experimental data, and the 5 × 20 × 11 grid accurately reproducing band gap and optical absorbance trends at the PBE + U level.

Using a multi-tiered first-principles approach, the crystal structure prediction and energetic analysis of Gd–Ga₂O₃ were carried out. In order to conduct global structural searches without symmetry constraints and objectively identify low-energy doped configurations, the CALYPSO (Crystal structure AnaLYsis by Particle Swarm Optimization) code⁵¹ was utilized. The CALYPSO framework uses a particle swarm optimization (PSO) technique to repeatedly develop a population of structures. A specific proportion of new structures are produced using updated velocity vectors in each generation. These velocity vectors take into account the worldwide lowest enthalpy structure found in the current generation as well as data from earlier optimized structures. Stability was assessed through enthalpy calculations from CALYPSO global structure prediction, which identified it as a thermodynamically low-energy configuration. Until a predetermined maximum number of generations (GenNumb) is achieved, this iterative process keeps on. We set the GenNumb to 70 and the population size (PopSize) at 30 for our study. It is generally understood that in order to guarantee a comprehensive investigation of the potential energy surface and to validate the stability of the lowest-energy structure throughout several generations, systems with a larger number of atoms require bigger PopSize and GenNumb values.⁵¹

The VASP was used to apply a sequential four-stage relaxation process to every structure produced by CALYPSO. The PAW method was used to perform the DFT computations. The high correlations in O-p, Ga-d, and Gd-f electrons were taken into consideration by using the PBESol with an additional on-site Hubbard U correction (PBESol + U). In order to obtain precise total energies and structural parameters, the four-stage relaxation scheme consisted of: (I) initial pre-relaxation at moderate cutoff energy and k-point density; (II) intermediate relaxation with increased accuracy; (III) high-accuracy optimization using PBESol + U with tight convergence thresholds; and (IV) final refinement using the hybrid HSE06 functional.

The GGA + U functional (PBEsol + U) is used in stages I–III with ever-increasing precision: Low (ENCUT = 450 eV, force tolerance = 0.1 eV Å⁻¹, KSPACING = 2 × 2 × 2), Normal (ENCUT = 500 eV, force tolerance = 0.05 eV Å⁻¹, KSPACING = 2 × 2 × 2), and High (ENCUT = 600 eV, force tolerance = 0.02 eV Å⁻¹, KSPACING = 3 × 3 × 3) are the three types. For high-accuracy electronic and energetic evaluation, Stage IV uses the hybrid HSE06 functional with dense k-point mesh (KSPACING = 11 × 11 × 11), ENCUT = 600 eV, and force tolerance = 0.02 eV Å⁻¹. For accurate structure prediction and validation, this pipeline strikes a balance between computing efficiency and accuracy. Through the comparison of total energies and the search for consistent coordination environments, the final structures were verified to be both dynamically and thermodynamically stable.

3. Results and discussion

3.1 Monoclinic β-Ga₂O₃ and Gd-doped β-Ga₂O₃

This section provides a thorough and in-depth examination of our results, arranged into four main subsections to help readers grasp how Gd doping affects monoclinic β-Ga₂O₃. We start by carefully analyzing the energy stability and optimal structural geometries of both pure and Gd–Ga₂O₃, paying special attention to the favored integration sites and the ensuing local structural distortions. We then analyze the electrical structures and intrinsic properties of these systems, such as their band structure and density of states, to elucidate the fundamental changes induced by Gd doping. The optical characteristics are then thoroughly covered in the third subsection, with particular attention paid to the analysis of the dielectric functions and absorption spectra in order to comprehend their potential for optoelectronic uses. In order to evaluate the practical utility of both the pristine and Gd–Ga₂O₃ systems, we lastly examine their mechanical properties, focusing on their stability and elasticity under mechanical stress.

3.1.1 Structural geometries and energetic stability of β-Ga₂O₃ and Gd–Ga₂O₃. The HSE06 hybrid functional with a 2 × 8 × 4 k-point sample grid was used to optimize the structure of the monoclinic β-Ga₂O₃ crystal. In Fig. 1(a), the optimized structure is shown. Table 1 displays the computed lattice parameters for comparison with previously published theoretical and experimental values. For β-Ga₂O₃, we have determined the lattice parameters a = 12.29 Å, b = 3.06 Å, c = 5.84 Å, and β = 103.75°. As shown in Table 1, these values are in a good accordance with experimental results⁴¹ and previously published theoretical calculations using hybrid functionals (e.g., HSE06 (ref. 52) and B3LYP⁵³). With variations of 0.06 Å, 0.02 Å, and 0.04 Å, respectively, our results for the a, b, and c lattice parameters are particularly near to the experimental values. Additionally, the computed monoclinic angle β of 103.75° falls between 103.71° and 103.9°, closely matching the observed value of 103.7° and other theoretical predictions. The strong correlation between our HSE06 results and experimental data highlights how well the selected computational approach and k-point sample density describe the structural characteristics of β-Ga₂O₃.

Table 1 Calculated structural parameters of β-Ga₂O₃ with HSE06 hybrid functional compared with previous theoretical and experimental results

Parameters	This work	HSE06 (ref. 52)	B3LYP^53	Expt.^51
a (Å)	12.29	12.27	12.34	12.23
b (Å)	3.06	3.05	3.08	3.04
c (Å)	5.84	5.82	5.87	5.80
β (°)	103.75	103.82	103.9	103.7

---

Substituting a Ga atom at both the octahedral (Ga-oct) and tetrahedral (Ga-tet) locations allowed us to assess Gd incorporation. According to our HSE06 hybrid functional calculations, it is energetically more advantageous to incorporate Gd into the β-Ga₂O₃ lattice when Gd replaces a tetrahedrally coordinated Ga site (Ga-tet). The computed formation energies provide substantial support for this conclusion. Under normal growth conditions, this structure is exceedingly unlikely and thermodynamically unstable due to the large positive value of +13.82 eV for the formation energy for Gd substituting an octahedrally coordinated Ga site (Ga-oct). In sharp contrast, a negative value of −1.29 eV was obtained for the production energy of Gd substituting a tetrahedrally coordinated Ga site (Ga-tet). The negative formation energy indicates that Gd will spontaneously absorb into these sites because the substitution of Gd at the Ga-tet site is energetically advantageous and consequently thermodynamically stable.

We noticed alterations in the local bond lengths surrounding the pristine and replaced sites upon Gd doping (see Fig. 1(b)). The determined Ga–O bond lengths for the tetrahedral locations in the pristine β-Ga₂O₃ were 1.841 Å, 1.839 Å, and 1.875 Å. Gd–O bond lengths at the tetrahedrally coordinated site were found to be longer in the Gd-doped system, measuring 2.067 Å, 2.132 Å, and 2.131 Å. This increase in bond length causes a local lattice distortion around the dopant site and is consistent with Gd having a larger ionic radius (0.94 Å) than Ga (0.47 Å).⁵⁴ The energetically advantageous substitution of Gd at the tetrahedral sites causes these structural alterations, which are essential for comprehending the electrical and optical characteristics of Gd-doped β-Ga₂O₃.

Our calculated formation energies for octahedral substitution show a large positive value, indicating strong repulsive interactions caused by Gd at this site. This suggests that such a configuration is thermodynamically unstable. This finding contrasts with an earlier study,²⁷ a discrepancy we attribute to differences in computational methodology. The previous work involved a two-step approach: initial structural relaxation used pseudopotentials where Gd f-electrons were treated as core electrons for computational efficiency, followed by electronic structure calculations with f-electrons treated as valence electrons, based on the pre-relaxed structures. In contrast, our study employed full HSE06-level calculations, consistently treating Gd f-electrons as valence electrons throughout both the structural relaxation and subsequent electronic structure analysis. This consistent treatment provides a more accurate representation of the electronic interactions, leading to the observed instability.

To directly assess the impact of the exchange–correlation treatment on dopant site energetics, we performed additional structural optimization and formation-energy calculations for Gd substituting at the octahedral and tetrahedral Ga sites using the PBE + U. The computed formation energies are −0.05 and +4.38 eV for Gd–oct and Gd–tet, respectively (under HSE06 the tetrahedral site is preferred by 15.18 eV (tet–oct = −15.18 eV), i.e., the relative site stability is reversed compared with PBE + U. These results demonstrate a qualitatively important difference: the choice of functional not only shifts absolute formation energies but changes the preferred substitutional site for Gd. The reversal and the large energetic differences reflect the strong sensitivity of formation energies to the treatment of exchange and self-interaction for localized 4f (and correlated d) electrons. Hybrid functionals such as HSE06 include a fraction of exact exchange that corrects self-interaction errors present in semilocal functionals and can strongly modify the localization, energetic position and hybridization of Gd 4f states; this in turn affects total energies and defect formation energies.

3.1.2 Electronic structure analysis. The electrical characteristics of pristine and Gd-doped β-Ga₂O₃ are described in this subsection, with a particular emphasis on the band structure, total density of states (DOS), and partial density of states (PDOS), as determined by our PBE + U calculations.

For pristine β-Ga₂O₃, our PBE + U calculations show that it is a semiconductor with a direct band gap. Fig. 1(c) shows that the computed band gap is 4.87 eV, which is astonishingly close to experimental data for β-Ga₂O₃ band gaps, which are generally about 4.9 eV.^55,56 Additionally, this result is consistent with earlier theoretical studies conducted at the DFT + U and HSE06 hybrid levels.^57,58 Crucially, the valence band maximum (VBM) and conduction band minimum (CBM) are both situated near the Brillouin zone's Γ point, clearly demonstrating its direct band gap nature. For many optoelectronic applications, this straight band gap is ideal.

As illustrated in Fig. 1(d), the overall band structure of Gd-doped β-Ga₂O₃ retains its semiconducting nature, despite a somewhat smaller band gap of 4.77 eV. The addition of Gd produces localized states inside the band gap or close to the band boundaries, but the material still has a large band gap. These Gd orbital-derived impurity states have the ability to greatly affect the electronic transport characteristics by serving as recombination or trapping sites or even opening up new channels for electronic transitions.

Visual confirmations of these changes would be provided by a thorough analysis of the band structure plots for both pristine and doped systems, which would reveal the development of Gd-induced states in the doped counterpart and the absence of mid-gap states in pristine β-Ga₂O₃. We carefully examined the total DOS and projected DOS (PDOS) for both pristine and Gd–Ga₂O₃ in order to obtain a better understanding of the contributions of specific atomic orbitals to the electronic states. In β-Ga₂O₃, O-p states mostly form the valence band, especially close to the Fermi level as can be found from the plot (see Fig. 1(e)). Additionally, Ga-p states make a noteworthy contribution, while Ga-s and Ga-d states make a lower but still important contribution. The strong covalent character of the Ga–O bonds inside the β-Ga₂O₃ lattice is highlighted by this complex hybridization. On the other hand, Ga-s states contribute most to the CBM, with Ga-p states and O-p orbitals playing a minor role. These results are essential for comprehending the carrier dynamics and electrical transitions in undoped gallium oxide.

The presence of the Gd impurity at the tetrahedral Ga site in Gd–Ga₂O₃ causes a significant change in the electronic structure. The valence band is still mostly dominated by O-p states, but the addition of Gd creates new, noticeable features. As can be seen from Fig. 1(f), our PDOS analysis clearly shows that the Gd-d and Gd-f states are important in forming the electronic structure, especially close to the Fermi level and, crucially, within the conduction band. In particular, Gd-d states contribute the most to the conduction band, followed by Gd-f states, while Gd-s and Gd-p states contribute quite slightly. The O-p states remain the primary contributors to the valence band in the Gd-doped system.

Numerous investigations have documented that the creation of defect levels or shallow impurity states within the band gap is often caused by the extremely restricted nature of f states in metal oxides.^59–61 Because these states open up new pathways for charge carriers and light absorption, they can have a substantial effect on a material's electrical and optical properties. In accordance with the theoretically expected and experimentally observed bond length variations, the substantial hybridization between Gd orbitals (particularly Gd-d and Gd-f) and nearby O-p states suggests the development of strong covalent bonds between Gd and O atoms. Since they can affect carrier dynamics and possibly open up new channels for optical transitions—which will be covered in more detail in the section on optical characteristics that follows—the appearance of these Gd-derived states close to the band boundaries is crucial.

3.1.3 Optical properties analysis. The optical characteristics of both pristine and Gd-doped β-Ga₂O₃ systems are thoroughly examined in this subsection and are computed using the PBE + U functional. In order to clarify the significant influence of Gd doping on the material's interaction with electromagnetic radiation, especially in the visible spectrum, we concentrate on the dielectric functions, energy loss functions, and absorption spectra.

The frequency-dependent dielectric function, ϵ(ω) = Re[ϵ(ω)] + iIm[ϵ(ω)], offers essential information on the mechanisms of electronic polarizability and energy absorption. Both the real and imaginary components of the dielectric function for pristine β-Ga₂O₃ (Fig. 2(a)) show distinctive little peaks, mostly in the ultraviolet (UV) region (corresponding to energies above the band gap). Im[ϵ(ω)], the imaginary part that is directly related to light absorption, exhibits first peaks at higher photon energies, which is consistent with the material's broad band gap and indicates very little absorption in the visible spectrum. The transparent character of β-Ga₂O₃ is indicated by the real portion, Re[ϵ(ω)], which is typically positive and rather constant in the visible range. Experimental studies^62–64 using UV–vis spectroscopy and spectroscopic ellipsometry report that optical absorption in β-Ga₂O₃ begins in the deep UV, consistent with the energy range where our computed features appear.

Fig. 2 Comparative optical properties of pristine and Gd-doped β-Ga₂O₃ obtained via PBE + U calculations (CASTEP). Panels (a) and (b) show the real and imaginary components of the dielectric function for pristine and Gd-doped β-Ga₂O₃, respectively. Panels (c) and (d) present the energy loss functions and absorbance spectra for both systems, highlighting the impact of Gd doping.

On the other hand, as Fig. 2(b) illustrates, the dielectric response of the Gd–Ga₂O₃ system is drastically changed, especially in the visible area (wavelengths greater than 400 nm). In the visible range, the Im[ϵ(ω)] and Re[ϵ(ω)] both exhibit wide, noticeable peaks. Notably, in some regions of the visible spectrum, the real part of the dielectric function, Re[ϵ(ω)], shows negative values. At those particular frequencies, when the material acts as a plasma and reflects received light instead of transmitting or absorbing it through typical interband transitions, this phenomenon is a clear sign of metallic or plasmonic-like behavior. Significant reflection or absorption may result from the material's electric field oscillating out of phase with the external field, as indicated by the negative real part. The additional Gd-d and Gd-f states that are introduced within the band gap and close to the conduction band boundaries are directly responsible for these broad characteristics in the dielectric function of Gd–Ga₂O₃, which enable novel low-energy optical transitions.

The energy loss function, L(ω) = Im[−1/ϵ(ω)], describes the energy lost by fast electrons moving across a medium and provides important information on the collective excitations within a material. Usually, plasma oscillations correlate to peaks in this function. At low wavelengths, about 100 nm, the loss function for β-Ga₂O₃ shows a strong peak (see Fig. 2(c)). The strong interband transitions seen in the UV absorption spectrum and the corresponding high-energy features in the imaginary part of the dielectric function, which represent collective excitations of the valence electrons, are closely correlated with this prominent characteristic. It's interesting to note that this peak decreases throughout the visible spectrum (400–700 nm), suggesting a markedly lower energy loss. At longer wavelengths (more than 1000 nm), it exhibits a modest rise, which might be explained by the material's other low-energy excitations or by residual free-carrier absorption.

The behavior of the energy loss function is noticeably more complex and enriched in the Gd–Ga₂O₃ system. According to the host material's natural UV response, the sharp peak at low wavelengths (about 100 nm) remains present, but as seen in Fig. 2(c), new, broad peaks clearly emerge throughout the visible area. The alterations in the dielectric response, especially the existence of negative values in the real part of the dielectric function (Re[ϵ(ω)]) within the observable range, are strongly supported by these novel features in the loss function. In contrast to the pristine material, this suggests a greater tendency for plasma oscillations or intense collective excitations at lower energy. The changed electronic structure brought about by Gd doping, where the introduction of Gd-d and Gd-f states produces more accessible electronic pathways and highly polarizable states, is directly responsible for these events. Together with the modifications to the system's dielectric response and absorption properties, the prominent features in the visible spectrum of the Gd-doped system's loss function highlight the substantial functional tunability made possible by Gd doping and its potential for new optoelectronic uses.

The material's capacity to absorb photons at various energies (wavelengths) is directly measured by the absorption coefficient, α(ω). In accordance with its broad band gap, pristine β-Ga₂O₃ displays an extremely sharp absorption edge (see Fig. 2(d)). In accordance with its huge band gap and the usual absorption properties of wide band gap semiconductors, a prominent and abrupt absorption peak is seen at around 100 nm (deep UV area) beyond this edge. Its transparency is confirmed by the fact that there is almost little absorption in the visible spectrum.

The absorption spectra change significantly and critically upon Gd-doping. As illustrated in Fig. 2(d), in addition to the intrinsic UV absorption from the host β-Ga₂O₃, a second, broad absorption peak stands out in the visible range, spanning from roughly 200 nm to 2000 nm with a minor maximum at 700 nm. The Gd incorporation is the direct cause of this recently observed wide absorption in the visible spectrum. Lower-energy photons (i.e., visible light) can be absorbed thanks to the Gd-d and Gd-f impurity states, which are covered in the section on electronic structure. They either serve as intermediate energy levels or offer new accessible final states for electronic transitions. Gd–Ga₂O₃ is a viable option for applications that need visible light interaction, including photocatalysis or visible-light photodetectors, because of its increased absorption in the visible spectrum.

3.1.4 Mechanical properties. According to calculations of elastic constants using the PBESol functional on HSE06-optimized structures, the mechanical characteristics of pure β-Ga₂O₃ and Gd–Ga₂O₃ (containing Gd at the tetrahedral site) are described in detail in this subsection. Determining a material's strength, stiffness, and appropriateness for a certain application requires an understanding of these characteristics.

β-Ga₂O₃ is a monoclinic crystal structure with 13 distinct elastic constants (space group C2/m). Below are the determined stiffness tensors (C_ij in GPa) for the pristine and Gd-doped systems. All independent elastic constants must meet the Born stability requirements in order for a material to be considered mechanically stable. The mechanical stability of both systems is confirmed by the fact that all estimated elastic constants for both pristine and Gd–Ga₂O₃ meet these requirements.

Stiffness Tensor C_ij (in GPa) for pristine β-Ga₂O₃:

Stiffness Tensor C_ij (in GPa) for Gd–Ga₂O₃:

The overall magnitudes of the diagonal elements (C₁₁, C₂₂, C₃₃, C₄₄, C₅₅, and C₆₆) for Gd–Ga₂O₃ are generally smaller than those for pristine β-Ga₂O₃, according to a comparison of the elastic constants. Given that the higher ionic radius of Gd replaces the smaller Ga atom, causing local lattice distortions and a softening of the material, this implies a decrease in stiffness and resistance to deformation following integration of Gd.

Using the Voigt–Reuss–Hill (VRH) averaging approximation, we obtained a number of macroscopic mechanical properties from the computed elastic constants, which are shown in Tables 2–4.^65–68 The bulk response of the material to stress can be better understood thanks to these characteristics.

Table 2 Comparison of elastic constants and derived mechanical properties with DFT-PBESol and Voigt–Reuss–Hill approximation for pristine β-Ga₂O₃

Property	This work	AM05 (ref. 64)	GGA^65	PBESol^66	Exp.^63
C11 (GPa)	237.08	223	199	208	243
C12 (GPa)	120.60	116	112	118	128
C13 (GPa)	126.34	125	125	146	160
C15 (GPa)	−14.86	−17	−2	0	−1.6
C22 (GPa)	328.83	333	312	335	344
C23 (GPa)	74.31	75	62	83	71
C25 (GPa)	13.51	12	1	0	0.4
C33 (GPa)	329.39	330	298	318	347
C35 (GPa)	8.91	7	17	19	1
C44 (GPa)	49.53	50	39	50	48
C46 (GPa)	16.86	17	3	9	5.6
C55 (GPa)	69.53	69	77	77	89
C66 (GPa)	93.52	94	95	96	104
Derived properties
BH (GPa)	170.33	167	155	171	183
EH (GPa)	195.85	194	182	192	210
GH (GPa)	74.88	74	70	73	80
νH	0.31	0.31	0.31	0.31	0.31

---

The mechanical parameters we computed for pristine β-Ga₂O₃ using the PBESol functional correspond well with experimental data and previous theoretical and experimental research, particularly those that use PBESol or related GGA functionals (see Table 2). This gives confidence in the computed trends for the doped system and verifies our computational approach. Our determined mechanical moduli and elastic constants for pure β-Ga₂O₃ closely match experimental results and are in good agreement with earlier theoretical investigations, especially those that use the PBESol functional.^69–72 Young's Modulus (E) and Shear Modulus (G) show comparable agreement. Variations in computational settings, pseudopotentials, and, in particular, our use of HSE06-optimized structures as input for the PBESol mechanical computations, are likely to cause minor variations. These variations may result in slightly different equilibrium volumes, which in turn may affect elastic characteristics. The accuracy of our computational methods for β-Ga₂O₃ is validated by the consistency observed throughout this research.

It is clear from Tables 3 and 4 that Gd doping causes a considerable drop in the average Bulk Modulus (B), Young's Modulus, Shear Modulus, and P-wave Modulus. For example, under hydrostatic pressure, the average Bulk Modulus drops from roughly 170.33 GPa (pristine) to 135.67 GPa (Gd-doped), suggesting that the Gd-doped material is noticeably less resistant to volume change. A significant reduction in stiffness and elastic response to uniaxial stress is also suggested by the decreased average Young's modulus (from about 195.85 GPa to approximately 142.48 GPa).

Table 3 Mechanical properties of pristine β-Ga₂O₃ with PBESol functional and Voigt–Reuss–Hill averages

Property	Voigt (GPa)	Reuss (GPa)	Hill (GPa)	Average (GPa)
Bulk modulus B	170.87	169.785	170.326	170.33
Young's modulus E	209.37	182.246	195.931	195.85
Shear modulus G	80.79	68.975	74.881	74.88
Poisson's ratio ν	0.30	0.321	0.308	0.31
P-wave modulus	278.58	261.752	270.168	270.17
Pugh's ratio (B/G)	2.12	2.462	2.275	2.29

---

Table 4 Mechanical properties of Gd-doped β-Ga₂O₃ (PBESol/Voigt–Reuss–Hill averages)

Property	Voigt (GPa)	Reuss (GPa)	Hill (GPa)	Average (GPa)
Bulk modulus B	136.42	134.928	135.674	135.67
Young's modulus E	166.64	117.970	142.821	142.48
Shear modulus G	64.27	43.554	53.913	53.91
Poisson's ratio ν	0.30	0.354	0.325	0.33
P-wave modulus	222.11	193.001	207.558	207.56
Pugh's ratio (B/G)	2.12	3.098	2.517	2.58

---

Additionally, the average Shear Modulus drops (from around 74.88 GPa to about 53.91 GPa), suggesting a decreased ability to withstand shear deformation. The idea that adding Gd makes the material mechanically less stiff is further supported by these reductions, which are in line with the softening effect shown in the elastic constants.

For β-Ga₂O₃, the average Poisson's Ratio (v) is roughly 0.31; for the Gd-doped system, it slightly increases to 0.33. Even though this shift is tiny, it indicates that doping may somewhat increase the material's propensity to deform perpendicular to the applied force. In addition, pure β-Ga₂O₃ has a Pugh's Ratio (B/G) of around 2.29, which rises to about 2.58 for the Gd-doped system. B/G > 1.75 denotes ductility, whereas B/G < 1.75 implies brittleness. This implies that β-Ga₂O₃ remains ductile in both its pristine and Gd-doped states, with the doped system possibly showing somewhat increased ductility.

One important observation is the reported decrease in hardness and stiffness with Gd doping. Because Gd has a larger ionic radius than Ga, the lattice local to the dopant site expands and weakens, resulting in this softening effect while retaining ductility. Applications needing particular mechanical resilience and material processing may be affected by this characteristic.

3.1.5 MLFF-Accelerated AIMD simulation of Gd-doped β-Ga₂O₃. We used machine-learned force fields (MLFFs) based on HSE06 hybrid functional level to speed ab initio molecular dynamics (AIMD) simulations in order to assess the thermal stability of Gd–Ga₂O₃. Targeting both realistic working circumstances and thermal endurance limits for optical and electrical applications, the simulations were carried out at two critical temperatures: 300 K (ambient) and 700 K (elevated).

Our MLFFs-AIMD simulations reveal that the structure maintained its crystallographic framework and local coordination contexts at 300 K, remaining stable throughout the AIMD trajectory. Initially replaced at a tetrahedral Ga site, the Gd dopant stayed stable in its local surroundings, showing no signs of atomic migration or lattice drift. The system's dynamical stability was further supported by the temperature and total energy changes over time.

During the simulation window, the overall lattice symmetry and bonding network remained intact, despite an increase in thermal agitation at 700 K, which was caused by amplified atomic vibrations. Notably, there was no site diffusion or bond breakdown, demonstrating exceptional heat resistance. The feasibility of Gd–Ga₂O₃ in high-temperature optoelectronic devices depends on this observation.

Temperature-dependent elongation was discovered by a thorough examination of the Gd–O bond lengths. As displayed in Fig. 3(a), the average Gd–O distances at 300 K were around 2.00 and 2.10 Å for Gd–O(I) and Gd–O(II), respectively, which were in good agreement with the results of static HSE06 optimizations. At 700 K, the Gd–O(II) bonding distances exhibited a slight increase, reaching 2.70 Å, as presented in Fig. 3(b). The average lengths of Gd–O(I) bonds remained constant despite an increase in instability. Initial bond length fluctuations, such as those that reach ∼3.2 Å at the beginning of the 700 K simulation (Gd–O(II)), are typical of MD simulations, particularly at higher temperatures. These are the initial relaxation and thermalization of the system. The initial few picoseconds, usually about 5 ps, are usually regarded as an equilibration phase and are not indicative of the final stable state. The data following this equilibration is the main focus of the analysis. This behavior is consistent with normal thermal expansion and does not suggest a destabilizing structural change. At 700 K (see Fig. 3(b)), a slight increase in Gd–O(II) bond lengths to 2.70 Å was noted. Gd–O(I) bond lengths increased in volatility but stayed stable. This behavior does not point to a destabilizing structural transition and is compatible with normal thermal expansion. Crucially, no coordination shift or bond rupture was seen, demonstrating that Gd substitution at the tetrahedral location maintains its structural stability even in the face of heat stress.

Fig. 3 Evolution of Gd–O bond distances in Gd-doped β-Ga₂O₃ during MLFF-accelerated AIMD-HSE06 simulations. (a) Gd–O(I) and Gd–O(II) bond distances at 300 K, demonstrating local structural stability. (b) Gd–O(I) and Gd–O(II) bond distances at 700 K, indicating thermal expansion without bond rupture or coordination change.

These findings support previous static predictions that the energetically advantageous Gd substitution site is the tetrahedral Ga site. The thermodynamic and kinetic stability of this configuration is demonstrated by the doped structure's capacity to retain its geometry at high temperatures. Additionally, MLFF-accelerated AIMD based on hybrid functional data (HSE06) is a useful method for evaluating finite-temperature behavior in doped complex oxides since it drastically lowers computational costs without sacrificing accuracy.

3.2 Crystallographic features of the Gd-doped structure

To accurately determine the substitutional site preference of Gd in the Ga₂O₃ lattice, we employed a multi-stage computational strategy combining high-accuracy DFT relaxations with VASP and global structure prediction using the CALYPSO framework. CALYPSO, based on particle swarm optimization (PSO), explores the potential energy surface (PES) in an unbiased manner, making it well suited for systems with significant local distortions and complex dopant environments. In the previous section (3.1), Gd substitution in monoclinic β-Ga₂O₃ was investigated to identify site preference within the experimentally known lattice, considering both tetrahedral and octahedral Ga sites. Complementarily, CALYPSO global searches were performed to uncover alternative low-energy phases, leading to the discovery of a thermodynamically stable triclinic phase with high-coordination Gd sites, structurally and energetically distinct from the monoclinic phase. These combined analyses provide a comprehensive picture: (i) the local site preference and property modifications in the established monoclinic lattice, and (ii) the identification of a novel global minimum structure exhibiting unique electronic and optical characteristics.

3.2.1 Structural geometries. Our CALYPSO-guided simulations demonstrate that Gd preferentially occupies a high-coordination icosahedral (20-fold coordinated) position in the Ga₂O₃ host lattice, as illustrated in Fig. 4(a), in contrast to original expectations preferring tetrahedral locations. In contrast, the host structure's Ga atoms are stable in both square-planar and octahedral environments, which is in line with the β-Ga₂O₃ coordination motifs that are known to exist. Comparative total energy assessments after a stringent four-stage relaxing technique validated this updated dopant placement, which appeared consistently across the lowest-enthalpy structures produced during the global search.

Fig. 4 Structural and electronic properties of the predicted triclinic Gd-doped Ga₂O₃ phase. (a) Optimized triclinic unit cell of Gd-doped Ga₂O₃ (left: polyhedral and right: ball-and-stick style representation), showcasing Gd at its preferred icosahedral coordination site, obtained from CALYPSO-VASP global structure prediction and subsequent PBESol + U/HSE06 hybrid functional relaxation. (b) Electronic band structure of the triclinic Gd-doped system, highlighting its wide band gap, computed at the PBE + U level. (c) Corresponding PDOS illustrating orbital contributions to the electronic structure. (d) Calculated absorbance spectrum, showing prominent absorption in both the DUV and visible regions, determined using the PBE + U functional in CASTEP.

Physically speaking, the preference for a high-coordination site makes sense because Gd³⁺ (0.94 Å) is substantially larger than Ga³⁺ (0.47 Å) and is hence more stable in prolonged coordination settings with reduced steric strain. The structural adaptability of the Ga₂O₃ lattice and its ability to locally accept large trivalent ions without suffering a substantial loss of crystallographic integrity are demonstrated by the observed icosahedral coordination.

The doped structure takes on a triclinic distortion with lattice parameters after hybrid functional relaxation:

α = β = γ ≈ 109.47°

a = b = c ≈ 6.52 Å

Strong local strain fields brought on by Gd incorporation; the lack of symmetry constraints in CALYPSO's P1 search, which permits complete relaxation of lattice degrees of freedom; and the soft lattice nature of Ga₂O₃, which can withstand rare-earth-induced distortion without sacrificing thermodynamic stability, are the reasons for this departure from the monoclinic C2/m symmetry of pristine β-Ga₂O₃.

This result highlights how important the coordination environment is for stabilizing dopants and shows how incorrect it is to assume straightforward replacement into low-coordination sites (such tetrahedral) based only on ionic size. Instead, the choice of dopant location is determined by local bonding topology, electrostatic balance, and spatial accommodation. Additionally, hybrid functional is crucial in this situation because it correctly depicts the subtle electronic interactions involving localized f-electrons of rare-earth dopants, which ordinary DFT or even DFT + U might overlook.

In brief, our results demonstrate the CALYPSO-DFT approach's resilience in revealing dopant site preferences that are unexpected yet physically significant. A more accurate and updated structural model for Gd-doped Ga₂O₃ is presented by assigning Gd to an icosahedral environment as opposed to Ga occupying octahedral and square-planar regions. This highlights the significance of correct electrical treatments and global structural search in the study of dopant action in complex hosts, and it has direct consequences for comprehending the electronic, optical, and mechanical characteristics of such doped oxides.

3.2.2 Electronic structure. For optoelectronic applications, the band gap and other electrical characteristics of Gd–Ga₂O₃ systems are essential. Our PBE + U functional analysis sheds light on the effects of structural phase transitions and Gd doping on the electronic band structure and density of states.

The newly anticipated triclinic Gd-doped Ga₂O₃ phase, in contrast, shows a much higher band gap of 7.02 eV as shown in Fig. 4(b). This obvious increase is remarkable—about 2.13 eV higher than pristine β-Ga₂O₃ and 2.25 eV greater than the monoclinic Gd-doped phase. This triclinic phase may be a suitable choice for next-generation DUV optoelectronic devices, potentially outperforming existing β-Ga₂O₃ materials due to its wide band gap. A number of mechanisms, including structural reconfigurations that alter the orbital interactions and possibly a more efficient suppression of defect states within the band gap because of the special atomic arrangement, are likely responsible for the triclinic phase's larger band gap.

A thorough examination of the PDOS for the triclinic Gd–Ga₂O₃ is necessary to comprehend the cause of this noticeably expanded band gap. We can identify the precise electronic interactions controlling the band boundaries thanks to the PDOS, which clarifies the contributions of individual atomic orbitals to the total density of states (see Fig. 4(c)).

Valence Band (VB) contributions. O-p orbitals are the main contributors to the valence band states, which generally dictate the material's bonding properties and the upper edge of the band gap. Because of their greater energy and widespread hybridization, oxygen p-states typically compose the main portion of the valence band in oxides. The Ga-d orbitals contribute significantly to the lower energy areas of the valence band after the O-p contributions. Significant contributions are also seen from Ga-p orbitals when they approach the Fermi level. More interestingly, the valence band is also influenced by the Gd-d, Gd-p, and Gd-f orbitals. The presence of Gd-f states, which are usually deeper in energy and confined, indicates hybridization effects with the host lattice orbitals, which could affect the stability and overall electronic structure of this new triclinic phase. Gd states' active involvement in the bonding network is indicated by their presence in the valence band.
Conduction Band (CB) contributions. Ga-p and Ga-s orbitals mostly dominate the conduction band states, which define the lower edge of the band gap and are in charge of electron transport. The main constituents of the conduction band minimum are these empty states of gallium. Furthermore, Gd-f orbitals contribute significantly, especially at energies near the conduction band minimum. Lesser contributions from Gd-d orbitals come next. A complex electronic interaction inside the Gd-doped triclinic structure is suggested by the presence of Gd-f orbitals in the valence band and their involvement in the conduction band.

The final band gap value depends on the hybridization between Ga orbitals and Gd f-states in both the valence and conduction bands, which is regulated by the special triclinic symmetry. The observed 7.02 eV band gap is supported by the separation of these contributing states, specifically the upward movement of the conduction band minimum or downward shift of the valence band maximum caused by these particular orbital interactions.

More thorough examination of the electronic band structure and projected DOS is necessary to comprehend how these particular orbital interactions cause the notable band gap opening. This could reveal areas of decreased overlap or greater energy separation between important bonding and anti-bonding states.

3.2.3 Optical properties. Important information about a material's electronic transitions and optoelectronic application possibilities can be found in the optical absorbance spectrum. By contrasting it with pristine β-Ga₂O₃ and monoclinic Gd–Ga₂O₃, we analyze the expected absorbance properties of the Gd–Ga₂O₃ triclinic phase.

The absorbance spectrum of the predicted triclinic Gd–Ga₂O₃ phase reveals two well-defined peaks and a markedly richer spectral profile (Fig. 4(d)), highlighting its distinct optical response compared to both pristine and monoclinic phases.

3.2.3.1 First peak (UV region). The triclinic Gd–Ga₂O₃ has a prominent and sharp absorption peak in the deep UV region, just like its monoclinic counterparts. The broad direct band gap of 7.02 eV, which marks the beginning of basic interband electronic transitions, is intimately linked to this peak. This peak's remarkable strength and location highlight the material's potential for cutting-edge DUV optoelectronic devices, providing substantial absorption outside of the visible range and great transparency within.
3.2.3.2 Second peak (visible region). The existence of a broad peak in the visible light region is a noteworthy and distinctive characteristic of the triclinic Gd–Ga₂O₃ absorbance spectra. The visible absorption seen in monoclinic Gd–Ga₂O₃, which also exhibits some absorption in this region, is evocative of this peak. In contrast to the monoclinic Gd–Ga₂O₃, the triclinic phase exhibits a substantially higher intensity (height) of this visible peak, reaching about 60% stronger.

This apparent absorption peak's presence and increased intensity in the triclinic phase are especially noteworthy. Absorption in the visible spectrum is usually ascribed to sub-bandgap electronic states, in contrast to inherent band-to-band transitions in broad bandgap semiconductors. These conditions may result from a number of mechanisms.

3.2.3.3 Gd-f intra-configurational transitions. Due to their high localization, Gd's f electrons have the ability to change within the 4f shell, creating distinctive absorption bands in the visible or near-UV spectrum. The stronger oscillator or more conducive environment for these transitions in the triclinic lattice is suggested by the greater intensity.
3.2.3.4 Defect states. Although the PBE + U functional often offers acceptable band gaps, observable absorption may be caused by impurity levels or residual defect states within the band gap. The triclinic phase's distinct structural configuration may produce new defect states or alter preexisting ones, improving absorption.
3.2.3.5 Charge transfer transitions. The visible spectrum may also include transitions involving the transfer of charge between Gd ions and oxygen ions or between distinct Gd sites. Such transitions may be more successfully facilitated by the triclinic phase's changed local environment and bonding.

A higher effectiveness in absorbing visible light is indicated by the triclinic phase's increased absorbance in the visible region when compared to the monoclinic Gd–Ga₂O₃. This feature could be used in photocatalysis, where visible light harvesting is essential, or in applications where regulated absorption of visible light is desirable, such photodetectors with customized spectral response. A bifunctional optical material with promise for a variety of applications is suggested by the coexistence of a strong DUV absorption and a sizable visible absorption peak.

3.2.4 Mechanical properties. The mechanical properties of the Gd–Ga₂O₃ triclinic phase predicted by CALYPSO-VASP are thoroughly analyzed in this section and compared to those of β-Ga₂O₃ and Gd–Ga₂O₃ (monoclinic phase with Gd in tetrahedral site). The impact of structural phase and Gd doping on ductility, hardness, and elastic moduli will be discussed. The PBESol functional calculations accompanied with Voigt, Reuss, and Hill approximations are used to derive the mechanical properties, such as B, E, G, v, P-wave Modulus, and Pugh's Ratio (B/G). For polycrystalline aggregates, the Hill average yields a good estimate (see Table 5). Compared to both pristine and monoclinic Gd–Ga₂O₃, the mechanical properties of the Gd–Ga₂O₃ triclinic phase are very different. With a Hill average of 394.89 GPa, the Bulk Modulus has dramatically increased, which is the most noticeable aspect. This value is over three times that of the monoclinic Gd-doped material and more than double that of pristine β-Ga₂O₃. The anticipated triclinic phase appears to have a far higher incompressibility due to its very high Bulk Modulus, which indicates that it is very resistant to volume change under hydrostatic pressure. This might be explained by an atomic arrangement that is more structurally stiff or tightly packed and stable in this new phase.

Table 5 Mechanical properties of triclinic Gd-doped Ga₂O₃ (PBESol/Voigt–Reuss–Hill averages)

Property	Voigt (GPa)	Reuss (GPa)	Hill (GPa)	Average (GPa)
Bulk modulus B	394.89	394.886	394.886	394.887
Young's modulus E	235.25	100.992	169.45	168.564
Shear modulus G	83.97	34.649	59.311	59.31
Poisson's ratio ν	0.4	0.457	0.428	0.428
P-wave modulus	506.85	441.084	473.968	473.967
Pugh's ratio (B/G)	4.7	11.397	6.658	7.585

---

On the other hand, the behavior of the Young's modulus and shear modulus varies considerably in the Voigt, Reuss, and Hill approximations. For Young's modulus and shear modulus, the Hill average values are 169.450 GPa and 59.311 GPa, respectively. The Young's modulus is higher than the monoclinic Gd-doped phase, indicating adequate stiffness, but it is lower than that of pristine β-Ga₂O₃. The Shear Modulus is much lower than that of pristine β-Ga₂O₃ and is similar to the monoclinic Gd-doped phase. With a Hill average of 0.428, the triclinic phase's Poisson's Ratio is exceptionally high. The very high Bulk Modulus is further supported by the fact that this value is near the theoretical limit of 0.5 for incompressible materials. When crushed longitudinally, a material with a high Poisson's ratio usually experiences significant lateral expansion, which is frequently linked to ductile behavior. The unusually high Pugh's Ratio of 6.66 supports this ductility finding even more. This number is much higher than 1.75, indicating that the triclinic Gd–Ga₂O₃ may be ductile rather than brittle, which is a very uncommon and perhaps advantageous oxide feature.

4. Conclusion

Using a variety of sophisticated computational methods, such as DFT with HSE06 hybrid functional, GGA + U, PBESol functional, MLFF-accelerated AIMD, and CALYPSO global structure prediction, our thorough first-principles investigations have yielded key insights into the structural, energetic, electronic, optical, mechanical, and thermal stability characteristics of pristine and Gd-doped β-Ga₂O₃.

Our calculations validated our methodology by precisely reproducing the known direct band gap semiconducting nature (4.87 eV) and structural properties of pure β-Ga₂O₃. Gd preferentially replaces the tetrahedral Ga site in monoclinic Gd-doped β-Ga₂O₃, causing local lattice distortions because of its greater ionic radius, according to preliminary investigations. MLFF-accelerated AIMD simulations verified its exceptional thermal durability in spite of these distortions, with the doped structure retaining its integrity up to 700 K. In terms of electronics, the addition of Gd created impurity states that allowed for broad absorption in the visible spectrum and a modest reduction of the band gap to 4.77 eV. The material's natural ductility was maintained despite considerable mechanical softening brought on by this initial doping.

Importantly, we discovered a new and extremely stable triclinic phase for Gd–Ga₂O₃ using our CALYPSO-guided global structural search. In this phase, Gd preferentially occupies a high-coordination icosahedral (20-fold coordinated) site, which contradicts more straightforward substitution assumptions and emphasizes the significance of thorough structural investigation. Compared to both pure and monoclinic Gd–Ga₂O₃, this recently identified phase shows a substantially wider direct band gap of 7.02 eV, a remarkable increase of over 2 eV. The triclinic phase is a strong contender for next-generation deep-ultraviolet (DUV) optoelectronic devices due to its high band gap. It exhibits enhanced and pronounced absorption of visible light optically, with a noticeable broad peak that is about 60% more intense than that of the monoclinic doped phase. Because of this, the material is ideal for uses like photocatalysis or sophisticated photodetectors that need to gather visible light. Mechanically, the triclinic phase exhibits a remarkable degree of ductility (Pugh's Ratio of 6.7) and a huge increase in Bulk Modulus (∼395 GPa), indicating unusually high incompressibility. This is an uncommon and highly desired combination for an oxide.

To sum up, Gd doping significantly alters the characteristics of β-Ga₂O₃. The finding of the stable triclinic phase offers an even more alluring prospect, even though the monoclinic doped phase provides controlled tunability and thermal resistance. Superior DUV capabilities, greatly improved visible light interaction, and strong mechanical integrity—including exceptional incompressibility—are all provided by this phase. Together, these thorough results highlight the significant functional tunability and robustness that Gd doping in β-Ga₂O₃ offers, making it a very promising material system for a variety of demanding optoelectronic applications, including those that operate at high temperatures.

Conflicts of interest

There are no conflicts to declare.

Data availability

All data supporting the findings of this study are included within the article.

Supplementary information: optimized atomic geometries with HSE06 hybrid functional for: 1- pristine β-Ga₂O₃, 2- monoclinic Gd-doped β-Ga₂O₃ (tetrahedral substitution), 3- monoclinic Gd-doped β-Ga₂O₃ (octahedral substitution), and 4- predicted triclinic Gd–Ga₂O₃ phase (icosahedral coordination). See DOI: https://doi.org/10.1039/d5ra05763a.

Acknowledgements

The work was supported by the Virtual Engineering Platform Project (Grant No. RS-2025-15373147), funded by the Ministry of Trade, Industry & Energy (MoTIE, South Korea).

References

1. M. Higashiwaki and G. H. Jessen, Appl. Phys. Lett., 2018, 112, 060401.
2. J. Y. Tsao, et al., Adv. Electron. Mater., 2018, 4, 1600501.
3. D. Guo, et al., Mater. Today Phys., 2019, 11, 100157.
4. B. Fu, Z. Jia, W. Mu, Y. Yin, J. Zhang and X. Tao, J. Semicond., 2019, 40, 011804.
5. A. J. Green, et al., APL Mater., 2022, 10, 061101.
6. A. Kuramata, K. Koshi, S. Watanabe, Y. Yamaoka, T. Masui and S. Yamakoshi, Jpn. J. Appl. Phys., 2016, 55, 1202A2.
7. M. Saleh, A. Bhattacharyya, J. B. Varley, S. Swain, J. Jesenovec, S. Krishnamoorthy and K. Lynn, Appl. Phys. Express, 2019, 12, 085502.
8. J. Jesenovec, J. Varley, S. E. Karcher and J. S. McCloy, J. Appl. Phys., 2021, 129, 225702.
9. E. Korhonen, F. Tuomisto, D. Gogova, G. Wagner, M. Baldini, Z. Galazka, R. Schewski and M. Albrecht, Appl. Phys. Lett., 2015, 106, 242103.
10. T. Onuma, S. Fujioka, T. Yamaguchi, M. Higashiwaki, K. Sasaki, T. Masui and T. Honda, Appl. Phys. Lett., 2013, 103, 041910.
11. E. Chikoidze, A. Fellous, A. Perez-Tomas, G. Sauthier, T. Tchelidze, C. Ton-That, T. T. Huynh, M. Phillips, S. Russel, M. Jennings, B. Berini, F. Jomard and Y. Dumont, Mater. Today Phys., 2017, 3, 118.
12. H. He, R. Orlando, M. A. Blanco, R. Pandey, E. Amzallag, I. Baraille and M. Rerat, Phys. Rev. B:Condens. Matter Mater. Phys., 2006, 74, 195123.
13. J. L. Lyons, Semicond. Sci. Technol., 2018, 33, 05LT02.
14. N. Ueda, H. Hosono, R. Waseda and H. Kawazoe, Appl. Phys. Lett., 1997, 70, 3561.
15. K. Sasaki, M. Higashiwaki, A. Kuramata, T. Masui and S. Yamakoshi, Appl. Phys. Express, 2013, 6, 086502.
16. E. Ahmadi, O. S. Koksaldi, S. W. Kaun, Y. Oshima, D. B. Short, U. K. Mishra and J. S. Speck, Appl. Phys. Express, 2017, 10, 041102.
17. B. Wang, D. Look and G. Farlow, J. Phys. D: Appl. Phys., 2020, 53, 444001.
18. A. Luchechko, et al., ECS J. Solid State Sci. Technol., 2020, 9, 045008.
19. X. Wang, et al., Superlattices Microstruct., 2019, 125, 330–337.
20. I. Hany, et al., Mater. Lett., 2019, 257, 126744.
21. S. Luan, et al., J. Alloys Compd., 2020, 812, 152026.
22. K. Sawada and T. Nakamura, J. Lumin., 2019, 215, 116616.
23. L. Chen, et al., J. Alloys Compd., 2020, 846, 156425.
24. S. Y. Lee, K. H. Choi and H. C. Kang, Mater. Lett., 2016, 176, 213–218.
25. Z. Wu, et al., Appl. Phys. Lett., 2015, 106, 171910.
26. P. Gribisch, A. Roy Chaudhuri and A. Fissel, ECS Trans., 2019, 93, 57–60.
27. E. Welch, P. Borges and L. Scolfaro, Physica B, 2023, 651, 414558.
28. J. Heyd, G. E. Scuseria and M. Ernzerhof, J. Chem. Phys., 2003, 118, 8207–8215.
29. C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti and C. G. Van de Walle, Rev. Mod. Phys., 2014, 86, 253.
30. Z. Suo, et al., Results Phys., 2020, 17, 103058.
31. M. D. Ganji, H. Song, J.-H. Park, D.-W. Jeon and H. Ko, Phys. B, 2025, 712, 417334.
32. I. Bouziani, et al., Physica A, 2018, 512, 1249–1259.
33. I. Bouziani, et al., J. Magn. Magn. Mater., 2019, 477, 220–225.
34. I. Bouziani, et al., Mater. Sci. Eng., B, 2020, 253, 114484.
35. S. Gao, W. Li, J. Dai, Q. Wang and Z. Suo, Mater. Res. Express, 2021, 8, 025904.
36. G. Kresse and J. Hafner, Phys. Rev. B:Condens. Matter Mater. Phys., 1994, 49, 14251–14269.
37. G. Kresse and J. Furthmuller, Comput. Mater. Sci., 1996, 6, 15–50.
38. G. Kresse and J. Furthmuller, Phys. Rev. B:Condens. Matter Mater. Phys., 1996, 54, 11169–11186.
39. G. Kresse and D. Joubert, Phys. Rev. B:Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
40. J. Paier, R. Hirschl, M. Marsman and G. Kresse, J. Chem. Phys., 2005, 122, 234102.
41. S. Geller, J. Chem. Phys., 1960, 33, 676–684.
42. H. J. Monkhorst and J. D. Pack, Phys. Rev. B:Condens. Matter Mater. Phys., 1976, 13, 5188–5192.
43. S. Poncé and F. Giustino, Phys. Rev. Res., 2020, 2, 033102.
44. A. S. Grashchenko, S. A. Kukushkin, V. I. Nikolaev, A. V. Osipov, E. V. Osipova and I. P. Soshnikov, Phys. Solid State, 2018, 60, 852.
45. A. V. Osipov, A. S. Grashchenko, S. A. Kukushkin, V. I. Nikolaev, E. V. Osipova, A. I. Pechnikov and I. P. Soshnikov, Continuum Mech. Thermodyn., 2018, 30, 1059.
46. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. I. J. Probert, K. Refson and M. C. Payne, Z. Kristallogr. Cryst. Mater., 2005, 220, 567–573.
47. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865.
48. A. I. Liechtenstein, V. I. Anisimov and J. Zaanen, Phys. Rev. B:Condens. Matter Mater. Phys., 1995, 52, R5467.
49. S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys and A. P. Sutton, Phys. Rev. B:Condens. Matter Mater. Phys., 1998, 57, 1505–1509.
50. B. G. Pfrommer, M. Côté, S. G. Louie and M. L. Cohen, J. Comput. Phys., 1997, 131, 233–240.
51. Y. Wang, J. Lv, L. Zhu and Y. Ma, Comput. Phys. Commun., 2012, 183, 2063–2070.
52. J. Varley, J. Weber, A. Janotti and C. Van de Walle, Appl. Phys. Lett., 2010, 97, 142106.
53. J. B. Varley, J. R. Weber, A. Janotti and C. G. Van de Walle, Appl. Phys. Lett., 2010, 97, 14.
54. S. Luan, L. Dong and R. Jia, J. Cryst. Growth, 2019, 505, 74–81.
55. M. Orita, H. Ohta, M. Hirano and H. Hosono, Appl. Phys. Lett., 2000, 77, 4166–4168.
56. S. Kumar, C. Tessarek, G. Sarau, S. Christiansen and R. Singh, Adv. Eng. Mater., 2014, 17, 709–715.
57. X. Ma, Y. Zhang, L. Dong and R. Jia, Results Phys., 2017, 7, 1582–1589.
58. H. Peelaers and C. G. Van de Walle, Phys. Status Solidi B, 2015, 252, 828.
59. A. Polman, J. Appl. Phys., 1997, 82, 1.
60. S. M. Sze and K. K. Ng, Physics of Semiconductor Devices, Wiley, 2007.
61. P. Wellenius, Rare-Earth Doped Wide Bandgap Oxide Semiconductor Materials and Devices, PhD thesis, North Carolina State University, 2009, https://repository.lib.ncsu.edu/items/e023437e-484a-4903-ab80-cd68565aa156.
62. F. Ricci, F. Boschi, A. Baraldi, A. Filippetti, M. Higashiwaki, A. Kuramata, V. Fiorentini and R. Fornari, J. Phys.: Condens. Matter, 2016, 28, 224005.
63. F. Ricci, F. Boschi, A. Baraldi, A. Filippetti, M. Higashiwaki, A. Kuramata, V. Fiorentini and R. Fornari, Adv. Mater. Sci. Eng., 2018, 2018, 9450157.
64. J. Zhang, J. Shi, D.-C. Qi, L. Chen and K. H. L. Zhang, APL Mater., 2020, 8, 020906.
65. V. Wang, N. Xu, J.-C. Liu, G. Tang and W.-T. Geng, Comput. Phys. Commun., 2021, 267, 108033.
66. W. Voigt, Lehrbuch der Kristallphysik, Teubner, Leipzig, 1928.
67. A. Reuss and Z. Angew, Math. Mech., 1929, 9, 49–58.
68. R. Hill, Proc. Phys. Soc., Sect. A, 1952, 65, 349–354.
69. K. Adachi, H. Ogi, N. Takeuchi, N. Nakamura, H. Watanabe, T. Ito and Y. Ozaki, J. Appl. Phys., 2018, 124, 085102.
70. J. Furthmüller and F. Bechstedt, Phys. Rev. B:Condens. Matter Mater. Phys., 2016, 93, 115204.
71. M. de Jong, W. Chen, T. Angsten, A. Jain, R. Notestine, A. Gamst, M. Sluiter, C. K. Ande, S. van der Zwaag, J. J. Plata, C. Toher, S. Curtarolo, G. Ceder, K. A. Persson and M. Asta, Sci. Data, 2015, 2, 150009.
72. A. S. Grashchenko, S. A. Kukushkin, V. I. Nikolaev, A. V. Osipov, E. V. Osipova and I. P. Soshnikov, Phys. Solid State, 2018, 60, 852.

---