Correct identification of local configurations of dopants and point defects in random alloys poses a challenge to both computational modeling and experimental characterization methods. In this paper, we propose and implement a computationally efficient approach to address this problem. Combining special quasirandom structures, virtual crystal approximation, and real-space lattice static Green’s functions, we are able to calculate, at moderate computational cost, the local phonon density of states (LPDOSs) of impurities in a random alloy crystal for system sizes, surpassing the capabilities of a conventional, cubic-scaling, density functional theory. We validate this method by showing that our LPDOS predictions of substitutional silicon in GaAs and InAs are in excellent agreement with the experimental data. For the case study, we investigate a variety of local configurations of Si and Se substitutional dopants and cation vacancies in quasirandom $In0.5Ga0.5As$ alloys. In all cases, the impurity LPDOS in a random alloy exhibits qualitatively different signatures from those in the pure binary compounds GaAs and InAs. Specifically, they are characterized by a wide continuous band (rather than narrow discrete peaks) of vibrational modes at frequencies typically higher than the bulk modes, a sign of coupling between localized vibrations of the impurity and those of its random neighboring host atoms. The accuracy and computational cost of this approach open a way to the simulation of impurities in random structures on a large scale and the prediction of vibrational signatures of alloys with defects.

## I. INTRODUCTION

Many technologically important semiconductor alloys, such as $Si1\u2212xGex$ and III–V compounds, exhibit a variety of different compositional orderings, depending on growth and annealing conditions.^{1,2} This is important since the spatial arrangements of atoms on their respective lattices often critically determine the alloy’s electronic and optical properties due to the preference of local impurities (dopants or point defects) to adopt specific local environments. For example, the alloy-induced electronic level splitting of an interstitial iron–boron pair is significantly enhanced by an increase in local Ge distribution in a $Si1\u2212xGex$ alloy, leading to significant deterioration of electrical conductivity in SiGe.^{3}

As a result, it would be extremely valuable if it were possible to obtain a local energy spectrum of atomic-level disorder (due to dopants and point defects) corresponding to specific spatially *local* atomic environments. For non-alloys and *ordered* alloys, phonon spectra can be retrieved either directly from experimental spectroscopic methods (such as infrared and Raman spectroscopy^{4,5}) or from *ab initio* computational methods such as density functional theory (DFT). The DFT results typically produce near-experimental accuracy,^{6–8} as local disorder typically lowers the point group symmetry in its immediate neighborhood compared to the rest of the crystal, leading to distinguishable features (peaks) in the measured spectra.

In contrast, for random and near-random alloys, such detailed information is very difficult to obtain from an experiment alone. This is because the point group symmetry of the parent lattice is broken almost everywhere in the random alloy, leading to significant background “noise” in the measured spectrum. This noise makes it very difficult to unambiguously resolve individual peaks due to the impurities alone. A suitable simulation route could provide a detailed map allowing the extraction of relative concentrations of local environments from the experimentally measured energy spectrum and a spatial distribution of disorder within the alloy lattice.^{3} Such a procedure could serve as an important step in predicting how to achieve effective doping of semiconductors, which is a critical task in creating electronic devices.

There have been a limited number of computational studies on dopant- or defect-induced vibrational modes in *near-random* semiconductor alloys. One prominent example of the difficulties in not considering truly random alloys can be illustrated in the work on a *ternary* III–V system that highlighted the resolution of preferred local bonding for nitrogen impurities in high In/Ga content InGaAs alloys. Experimental studies using Fourier transform infrared (FTIR) spectroscopy^{9–11} and Raman spectroscopy^{12–14} measured the local vibrational modes of nitrogen in InGaAs. However, due to random variations in the local environment, these experimental studies provide conflicting and ambiguous results. On the computational side, numerical simulations attempted to address the same challenge^{15} using approximate schemes such as empirical lattice dynamical models. These models require *a priori* calibration with experimental data and thus are system-dependent and lack generality and transferability. For example, the use of an Abell–Tersoff semi-empirical model to make predictions of the extrinsic “substitution energy” of a Si dopant on a cationic lattice site in InGaAs was found to be strongly dependent on the strain induced by a local arrangement of In and Ga cationic atoms.^{16} Moreover, current state-of-the-art *ab initio* density functional theory (DFT) codes often make extensive use of point group and space group symmetries for treating phonons in crystals,^{17} but which are not applicable to highly disordered random alloys. Several *ab initio* computational studies on defects in alloys have been performed.^{19,20} However, these studies only considered variations among chemical species of nearest-neighbor atoms, while ignoring the randomness of atomic arrangements in the larger host crystal. They also ignored any influence of dopants and defects on the lattice structure that extends over several shells of neighbors, which can be important depending on the lattice strain created by the dopant or defect. A more rigorous theoretical treatment must take these crucial aspects into account, while maintaining a feasible computational cost.

In this work, we tackle this challenge by presenting a systematic approach to modeling the vibrational modes of impurity atoms in a random alloy. For our case study, we choose two prototypical $n$-type dopants, silicon and selenium, and one prevalent compensating defect, the cation vacancy, in the framework of a random InGaAs alloy. Specifically, we use special quasirandom structures (SQSs) to model a random alloy within the framework of DFT and a real-space Green’s function to calculate local phonon density of states of a dopant or defect within an extended supercell that is beyond the typical capabilities of DFT calculations. This approach allows us to accurately determine the local vibrational modes of defects in a random alloy directly from DFT calculations, without any empirical parameters and, importantly, without extensive computational costs. Specifically, our method reveals the key qualitative difference between dopant and defect modes in *random* semiconductor alloys and in pure semiconductor compounds, namely, the formation of satellite peaks across a continuous range of frequencies.^{15} We will discuss the methodology in detail in Sec. II and provide results in Sec. III.

## II. METHODOLOGY

### A. Structural model of random In_{0.5}Ga_{0.5}As

For density functional theory (DFT) calculations, we use 216-atom cubic supercells of $In0.5Ga0.5As$ containing one impurity [a substitutional silicon on cation site ($SiIII$) (III denotes an In or Ga site), a substitutional selenium atom on an anion site ($SeAs$), or a cation vacancy ($VIII$)]. This particular size is chosen such that the impurity strain field, created by placing the impurity at the center of the supercell, is almost fully contained within the cell boundaries (above a small force threshold of 0.03 eV/Å) for all defect/dopant species. (More details are given in Sec. III A). This is a crucial test needed to minimize finite-size effects of elastic interactions between neighboring images of impurities. Compared with the typical supercell choices for defect in GaAs in DFT phonon calculations,^{18,19} our choice has larger dimensions, which means a larger range of defect–host interactions can be captured in our model.

The random arrangement of the cation species is modeled within the formalism of “special quasirandom structures.”^{21} We used an in-house code to generate candidate structures, minimizing the sum of cluster functions for six clusters: pairs with first-, second-, third-, and fourth-nearest-neighbor separations on the cation sublattice, as well as triplets and quadruplets with the first nearest-neighbor separation on the cation sublattice. We used simulated annealing^{22} to find the optimal 216-atom SQS model such that the sum of the cluster functions lies below a chosen threshold of $10\u22122$. Figure 1 shows the SQS model of random $In0.5Ga0.5As$ we used in all subsequent calculations and analysis.

We use the plane wave density functional code Quantum Espresso for all calculations in this work. We apply PBEsol functional with ultrasoft pseudopotentials, with an energy cutoff of 50 Ry and a density cutoff of 400 Ry. Due to the large number of expensive calculations that will be involved, we use a single $k$-point $2\pi a(14,14,14)$. This minimizes the computational cost while retaining the sufficient computational accuracy, based on past work which showed it to be the most accurate single $k$-point sampling for simple cubic systems.^{23} We verified that, upon atomic relaxation, this setting yields convergence of atomic forces to within 5 meV/$\xc5$. We use the Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm^{24} for atomic and cell relaxation while keeping the cubic symmetry of the cell, with a force criterion of $10\u22123$ Ry/bohr and a stress criterion of 0.5 kbar. The relaxed supercell of our $In0.5Ga0.5As$ SQS assumes an equal side length of 16.96 $\xc5$.

### B. Dynamical matrix

To calculate the interatomic force constants required for analyzing the phonon density of states, we apply the “frozen phonon” method^{25} by displacing each atom within the impurity strain field by $\xb10.05\xc5$ relative to its equilibrium position in each of the Cartesian directions, with the impurity placed at the center of the supercell. The resulting force constants are then used to construct corresponding entries of the much larger $3000\xd73000$ force constant matrix for the 1000-atom effective supercell of defective random $In0.5Ga0.5As$, which, in turn, is used to calculate an atom-projected local phonon density of states using a real-space Green’s function formalism (details are given in Sec. II C). This approach should in principle be more accurate than the mass disorder approximation,^{26} where only the mass of the defect, but not the interatomic force constants, is modified in the construction of the dynamical matrix. For materials such as InGaAs, incorporating force constant disorder in the dynamical matrix is especially important, as it dominates the contributions to phonon-related properties such as lattice thermal conductivity.^{27}

A simplified picture that illustrates our formulation is shown in Fig. 2. We divide all the atoms in the 216-atom defective supercell into two regions (1 and 2). Region 1 consists of atoms close enough to the impurity that their interatomic interactions are strongly influenced by the impurity. Region 2 consists of the rest of the atoms, whose interatomic interactions remain largely unaffected by the impurity. In our calculations, region 1 consists of all atoms within the strain field of the impurity, as determined by DFT (for details, see Sec. III A). We then embed this 216-atom supercell at the center of a larger, 1000-atom, cubic *effective* supercell of random $In0.5Ga0.5As$ to construct the total dynamical matrix. We denote all the atoms outside of the 216-atom supercell as region 3. The reason to use a larger effective supercell is that the size of the 216-atom supercell imposes a cutoff on the range of the interatomic force constants that is insufficient for a real-space Green’s function to yield a converged local phonon density of states. The total dynamical matrix, $D$, of the 1000-atom effective supercell can then be formally expressed in the following form:

In this representation, $Hi$ are *intra*-region dynamical submatrices, describing interactions between atom pairs both belonging to the same region. $Vij=VjiT$ are the *inter*-region dynamical submatrices, describing interactions between atom pairs belonging separately to two regions $i$ and $j$. Our previous 216-atom DFT calculations have already determined force constants for all pairs of atoms with at least one atom lying within the strain field of the impurity (region 1, designated as the “inner region”); thus, we have determined all elements of $H1$, $V12$, and $V21$.

To determine the rest of the total dynamical matrix, we use the “virtual crystal approximation,” (VCA) treating the interatomic interactions in regions less impacted by the impurity (regions $2+3$, designated as the “outer region”) as the corresponding interatomic interactions within a homogeneous medium of random $In0.5Ga0.5As$. The VCA force constants are calculated by averaging the force constants of bulk GaAs and of bulk InAs and then mapped onto the corresponding lattice positions of the 1000-atom effective supercell using space group symmetry operations. Note that since the impurity is sufficiently far from any atom in region 3, we can effectively approximate the interatomic force constant between the impurity and any atom in region 3 as zero. This treatment has minimal impact on the local vibrational modes of the impurity but is necessary in order to construct the larger dynamical matrix needed for proper convergence in the real-space Green’s function method. Finally, we symmetrize $D$ by iteratively applying the force constant sum rules together with the diagonal symmetrization procedure.^{25} Our formalism effectively reduces the heavy computational cost of phonon calculations, while retaining the required accuracy for the local vibrational spectrum of the impurities. Our approach is similar to the “combined dynamical matrix” method proposed by Shi and Wang.^{28} In all of our calculations, we use charge-neutral dopants and defects, since Coulombic interactions between the supercell images are long-ranged and require special treatment in phonon calculations.^{29} Experimentally determined vibrational modes of charged shallow dopants (such as Si and Se) are only slightly different from those of neutral ones, as the dopant charge is delocalized throughout the space.^{30} The vibrational modes of deep-level defects (such as cation vacancies) can shift by as much as several tens of cm$\u22121$ under different charge states,^{30} but this effect does not modify our qualitative conclusions (see Secs. III C 7 and IV).

### C. Local phonon density of states from a real-space Green’s function

To calculate the vibrational density of states in a crystalline solid, the conventional approach is to uniformly sample a sufficient number of $k$-points in the Brillouin zone. The vibrational density of states can then be constructed either from the sum of states over eigenvalues and **k**-points or from the imaginary part of the trace of the static lattice Green’s function. However, for defective crystals, this approach is not practical, since the real-space “unit cell” must be large enough to minimize interactions between periodic images of the defect, which renders such method prohibitively expensive. On the other hand, random crystals are characterized by their local environments, instead of periodicity. This suggests that the optimal method should take advantage of the locality of the atomic interactions in the real space, rather than the periodicity of the supercell.

Based on these considerations, we choose to use the real-space Green’s method^{31} that is particularly suitable for treating large disordered systems. With the total dynamical matrix $D$ calculated following the procedure in Sec. II B, we can construct the real-space Green’s function, $G$, defined as^{32}

where $\omega $ is the phonon frequency and $I$ is the identity matrix with the same dimensions as $D$. This expression is of little practical use, as inversion operation on a large matrix such as $D$ is infeasible to compute. On the other hand, since $D$ is symmetric by construction, we can use Haydock’s recursion method^{33} to transform $D$ to a tridiagonal matrix $D~$. From there, we can easily construct the diagonal elements of Green’s function by a continued fraction expansion using entries of tridiagonal $D~$^{34} and obtain the local phonon density of states (LPDOS) of the defect using the spectral expression^{35}

where $G\alpha j,\alpha j(\omega 2+i\epsilon )$ is the real-space Green’s function for atom $\alpha $ in direction $j$. According to our convergence tests, 300 iterations are adopted in the Lanczos recursions to construct the diagonal elements of Green’s function.

## III. RESULTS

### A. Local atomic environments

The vibrational signature of a substitutional dopant or defect (henceforth called an “impurity”) in a random crystal is strongly affected by its local atomic environment, as its interatomic interaction with different species of atoms generally have different strengths. An impurity can sit on either group III lattice sites or group V lattice sites in InGaAs. We searched and recorded the local environments throughout the 216-atom quasirandom supercell of random $In0.5Ga0.5As$, described by the number of In and Ga atoms in the nearest group III sites of the unrelaxed lattice. The results of the local environments for groups III and V sites are provided in the supplementary material.

Based on these results, we made the assumption that we can ignore any local environment with a $<2%$ probability of being substituted by a defect. With this assumption, we were able to mark six local environments for a group III site defect (Si and vacancy) and seven local environments for a group V site defect (Se), which are considered subsequently for vibrational calculations.

### B. Strain field induced by dopants and defects

Introducing an impurity into a pristine crystal reshapes the electrostatic interactions and atomic force constants in the neighborhood of the impurity. This creates a strain field, defined as the force response of each atom on the $relaxed$ $pristine$ lattice due to the presence of the defect. In DFT simulations of periodic crystals, if the strain field of the impurity is not fully contained within the supercell, elastic interactions between image impurities could occur, which introduces a finite-size error to the calculated phonon spectrum. Furthermore, the interatomic force constant between any two atoms within the strain field of the defect may be different from the bulk value due to the change in local lattice geometry. As a result, in order to accurately study the vibrational properties of an impurity, we must ensure that the supercell size is sufficiently large to fully contain the strain field of the impurity. The impact of the strain field surrounding the impurity on the interatomic force constants must also be fully taken into account in the DFT calculations.

In this work, we study three kinds of substitutional defects: $SiIII$, $VIII$ (cation vacancy), and $SeAs$ [subscript III denotes group III atom (In/Ga)]. We use a force criterion of 0.03 eV/A to determine if a particular atom lies within the strain field of the impurity. We verified that for all three impurities and all the local atomic environments considered in this work, the atoms within the strain field are fully contained within the 216-atom supercell boundary.

Figure 3 shows the strain fields of selected local atomic environments for all three impurities considered projected onto the XY-plane. Specifically, Si$III$ exhibits an ellipsoidally symmetric compressive field. This field induces strain on 30 surrounding atoms (reaching out to its fifth nearest neighbors), as shown in the leftmost picture of Fig. 3. Three of the four first nearest neighbors of the defect move *toward* the Si dopant atom, while the remaining atom stays in its original position for geometrical reasons. The presence of Si creates an off-center ion in which its equilibrium position is shifted away from the original lattice site. This is because Si is much smaller than the regular group III ions that it replaces (the atomic radius of Si is $rSi=1.11$ Å compared to $rIn=1.42$ Å and $rGa=1.22$ Å^{36}). Consequently, this significantly weakens the repulsive forces between the impurity ion and its nearest neighbors that stabilize the ion at a regular site.

The strain field around a $VIII$ defect also exhibits ellipsoidal symmetry. On the other hand, $VIII$ exerts a stronger and more isotropic compressive field with less symmetry reduction than that created by $SiIII$. The vacancy compressed around 55 atoms, extending its influence out to its sixth nearest neighbors. Similar to the Si-induced case, three of the vacancy’s four first nearest neighbors moved significantly toward the unoccupied site, while the remaining atoms remain in their original positions.

Unlike the case for the two defects mentioned above, the strain field around a $SeAs$ defect is asymmetric and more dispersive. Se creates a long-range tensile strain field, which distorts the zincblende structure and reduces the symmetry of tetrahedral bonds to a greater extent than either Si or vacancy. Atoms in a Se-doped crystal are more likely to shift toward one particular direction (+y). This can be understood by the fact that, in binary compounds such as InSe, selenium atoms tends to form a wurtzite structure, which is non-centrosymmetric, unlike the tetrahedral bonds in a zincblende lattice.^{37}

### C. Local phonon density of states

#### 1. Method validation

As mentioned in Sec. II B, we determine the *intra*-region force constants by displacing every atom within the defect’s strain field. Empirically, all the atoms within the defect’s first 3 nearest-neighbor shells may have non-negligible effects on the local density of states of the defect. But we noticed that the insertion of the impurity does not stress all the atoms inside these three shells. This is caused by the fact that the strains imposed on the host atoms may balance some forces by bringing the presence of the defect. Thus, some atoms that are not stressed by the defect may still have non-negligible interatomic interactions with the defect. To ensure we have considered all non-negligible interatomic interactions in our calculations, we compute a few more displacement jobs for those “non-overlapping” atoms (that is, atoms within the defect’s third-nearest-neighbor distance but not inside the defect’s strain fields). Since Si has the most non-overlapping atoms, compared to the other cases we studied, we conduct a test calculation to examine the contribution of those non-overlapping atoms to the LPDOS of the defect.

The results, shown in the supplementary material, confirm our assumption that it is unnecessary to calculate additional displacement jobs for the non-overlapping atoms since it will not change the local vibrational modes of the defect. The two curves shown have very similar shape with local modes (peaks) occurring at nearly the same frequencies. These results validate our approach of displacing only the atoms within the defect's strain field as an accurate and efficient method.

#### 2. Sensitivity analysis

To study how uncertainty in the force constants affects the vibrational density of states, we conduct several sensitivity tests based on a Gaussian process in which random variables have a multivariate distribution. Specifically, this distribution is a joint distribution of all those random variables and, as such, it depicts a distribution over functions with space. Thus, in our case, this process is a homogeneous (stationary and isotropic) process with zero mean, and the covariance function depends only on the Euclidean distance between each atom and the defect. We implement the Gaussian process module from scikit-learn (v0.20.1)^{38} to add noise to the force constants of each pair of atoms.

The results, shown in Fig. 4, provide evidence that our model is robust and that local vibrations remain unchanged even when the introduced noise on the force constants is as high as 0.005 in Rydberg atomic units. In particular, if the uncertainty of force constants is 0.001, the vibrational density of states perfectly matches the results without noise. Local vibrations start to exhibit minor inconsistencies with the original output when the introduced noise is 0.005, shown in the middle plot of Fig. 4, while the peaks still occur at the same frequencies. These results confirm the robustness of our approach.

#### 3. Single component substrate: Local phonon density of defects in Si

To validate our method, we calculate the local phonon density of states for point defects in a single component material, here Si. Specifically, we consider two defect species, germanium (Ge) and carbon (C); the defect atom is placed at the center of a cubic, 216-atom Si supercell, substituting one Si atom at that lattice site. Instead of calculating force constants for all the atoms within the defect’s strain field, we simply displace the defect atom and determine the force constants between the defect and other atoms in the supercell. Our results, shown in Fig. 5, are in reasonable agreement with the existing DFT work^{26} in terms of reproducing both the frequency and the localization character of the vibrational mode peak. We find that, as expected, the phonon mode of a lighter dopant (C) compared to the host atom species (Si) is localized at high frequencies with a relatively high intensity, whereas the phonon mode of a heavier dopant (Ge) appear at lower frequencies with a relatively low intensity. Results for InGaAs are also consistent with this finding (details are given in Secs. III C 5–III C 7).

#### 4. Binary cases: Local phonon density of states of defects in GaAs and InAs

To further validate our method, we apply our method to calculate the local phonon density of states for defects in binary alloys: GaAs and InAs. Our results, shown in Figs. 6 and 7, are in good agreement with existing experimental data. Specifically, the calculated vibrational modes of Si occur at 12.124 and 10.762 THz (404 and 359 cm$\u22121$) in GaAs and InAs, while the experimentally observed absorption peaks occur at very similar locations, namely, 399^{39} and 359 cm$\u22121$,^{40,41} at the current doping level. For $SeV$, the peaks of local phonon density of states in GaAs and InAs occur at 9.181 and 8.010 THz. The vibrational modes of the vacancy’s first-NNs are all resonant states where the peaks occur at frequencies ranging from 3.503 to 5.627 THz and 3.073 to 3.912 THz, respectively, in GaAs and InAs alloys.

#### 5. Local phonon density of states of Si_{III}

The local phonon density of states of $SiIII$ in the quasirandom $In0.5Ga0.5As$ alloy is shown in Fig. 8 for all the local environments we considered. Since the first nearest neighbors of Si, which have the most significant interactions with the defect, are all occupied by As atoms, the peaks in different local environments occur at similar frequencies. Specifically, the most intense vibrations occur at 11.492, 12.075, 10.964, 10.860, 11.389, and 11.544 THz, respectively, in each local environment. As expected, all the vibrational modes of $SiIII$ in a quasirandom $In0.5Ga0.5As$ alloy occur at frequencies between 10.762 THz and 12.123 THz, corresponding, respectively, to the frequency of $SiIn$ in InAs and that of $SiGa$ in GaAs based on results from our control group. We observe that when the number of Ga atoms equals, or is close to, the number of In atoms in the Si atom’s neighborhood (within a 3rd-NN distance), the peak occurs at a lower frequency. We do not see a specific pattern of vibrational modes with an increasing number of In/Ga atoms in the Si’s neighborhood. These indicate that the phonon modes of Si in a random III–V alloy do not depend very sensitively on the nearest-neighbor atom species.

The local vibrational modes of Si’s first nearest neighbors and a “bulk” atom (an atom which sits in the virtual crystal and is far from the defect) are shown in Fig. 9. The phonon modes of the bulk atom exhibit the global vibrations of the quasirandom $In0.5Ga0.5As$ crystal since this atom has negligible interactions with the defect. We can clearly see that the insertion of Si breaks the symmetry and disrupts the global vibrations. All four nearest neighbors react strongly to this interruption and show similar phonon modes with those of $SiIII$. The local modes in all scenarios occur at frequencies ranging from 10.3 to 11.8 THz. One of the four first nearest neighbors in each scenario, such as the one shown in the upper right plot in Fig. 9, displays weak vibrations (reflected as low intensity). It does not respond as strongly as the other three first nearest neighbors to the introduction of Si, which is consistent with our results for the strain fields (one atom is not strained by the defect), as discussed in Sec. III A.

#### 6. Local phonon density of states of Se_{As}

The local phonon density of states of $SeAs$ in the quasirandom $In0.5Ga0.5As$ alloy are shown in Fig. 10 for all the local environments considered. Compared to the vibrational modes of $SiIII$, those of $SeAs$ depend more sensitively on the local environment in both frequency and intensity due to variations of its first nearest neighbors in different scenarios. When the first-NN shell sites are all occupied by In atoms and most of the third-NN shell sites are occupied by Ga atoms (such as LE1 and LE2), the most intense peaks occur inside the range of the vibrational frequencies of the host lattice, making them resonant states. This may be because the difference in covalent radius of Se and As is much smaller compared to In vs Ga ($rGa=1.26$ Å, $rIn=1.44$ Å^{36}). Hence, the lattice structure around the Se dopant is similar to a central As atom with In atoms as all first-NNs, As atoms as all second-NNs, and Ga atoms as all third-NNs, which makes the local structure ordered and $SeAs$ showing a stronger global vibration. Under this circumstance, with an increasing number of In as third-NNs, the local structure becomes more disordered and Se shows relatively stronger local vibrations, which is reflected in much weaker global vibrations and more evident peaks in the local frequency zone in LE2 compared to LE1. In particular, the peaks occur at 4.487 and 4.176 THz (stronger global vibrations) under LE1 and LE2; and 7.606, 6.680, 7.814, 6.041, and 7.358 THz (stronger local vibrations) under LE3–LE7, respectively. Comparing to our results for the binaries, these peaks occur at slightly lower frequencies than in an ordered host lattice. Since more energy (a higher frequency) is required to vibrate a stronger bond, our results suggest that the bond strengths between Se atom and host atoms in the random lattice are weaker than those in an ordered lattice.

Since local atomic environments involve random atomic arrangements, multiple vibrational modes are observed even for the same type of neighboring atoms as the defect. The intensity of $SeAs$’s vibrations is much lower than that of $SiIII$ due to its larger atomic mass ($mSi=28.08$ amu, $mSe=79.92$ amu^{42}) and its similar covalent radius to As ($rSe=1.16$ Å, $rAs=1.19$ Å^{36}). These peaks occur at lower frequencies than those in the Si case, which is consistent with empirical research in which the vibrations of heavier defects result in lower frequencies.

The local vibrational modes of the first nearest neighbors of Se and a bulk atom are shown in Fig. 11. The most intense peaks in different local environments occur at frequencies ranging from 9.2 to 10.7 THz. The peaks of these first nearest neighbors shift to higher frequencies compared to that of the defect. In contrast to the results of the Si case, three of the four nearest neighbors of Se still have some global oscillations reflected as smaller peaks in the low-frequency zone. As expected, bulk-like atoms shows similar phonon density of states curves to those for the Si case.

#### 7. Local phonon density of states of V_{III}

The local phonon density of states of a cation vacancy’s first nearest neighbors in the quasirandom $In0.5Ga0.5As$ alloy are shown in Fig. 12 for all the local environments considered.

All the first nearest neighbors show a complex distribution of local vibrations, to different extents, in addition to global vibrations, as shown by the varying intensities of the bulk-like peaks in the LPDOS. Whereas the first-NNs of vacancies in GaAs and InAs show resonant modes in the phonon spectra, the first-NNs of vacancies in quasi-random InGaAs give rise to local vibration modes due to the disordered atomic arrangements. Unlike cases involving dopants, where localized modes often show higher intensities compared to bulk-like modes, here, global vibrations are more pronounced relative to local vibrational modes in the first nearest-neighbor projected modes. This is mainly because, for a vacancy, there is no impurity $atom$ to which the nearest-neighbor atoms are bonded. Hence, the local vibrational modes related to the impurity atom are absent. A vacancy lowers the local lattice symmetry less in comparison to Si and Se dopants, and this is also reflected in a more symmetrical strain field, as shown in Sec. III A. The vibrational signatures of the cationic vacancies suggest that the presence of such defects in a random InGaAs would introduce a noisy “background” of signals on top of the more distinguishable peaks of the dopants, such as Si and Se.

#### 8. Ensemble-averaged local phonon density of states

The ensemble average of the vibrational density of states of $SiIII$, $VIII$, and $SeV$ in a random InGaAs matrix are obtained by adding the weighted LPDOS of all possible local environments for 20 atoms around the defect (including the defect). The results are shown in Fig. 13. Experimental techniques cannot probe the phonon density of states on an atomic scale, so the “local” phonon spectrum is generated based on the nearest neighbors around the defect. The DOS summation of 20 neighboring atoms can preserve the local properties of the defect without introducing too much noise, while this size is detectable in experiments. Note that this number can be modified based on the specific precision in experiments.

## IV. CONCLUSIONS

In this work, we proposed a computational method that efficiently predicts the vibrational signatures of single-atom dopants/defects in a random alloy with fully *ab initio* accuracy. These predicted vibrational signatures can then serve as a guide for experimentalists to assess the dopant activation efficiency of novel doping schemes for semiconducting random alloys such as SiGe or InGaAs. Our calculations take into account a number of local environments for a defect site in a random alloy, which have not been considered in previous studies. We demonstrated that our method provides a quantitatively accurate description of phonon modes as evidenced by the excellent agreement we obtained between our simulated binary alloy results and experimental data.

Using this approach, we studied two prototypical $n$-type dopants, silicon and selenium, and one prevalent compensating defect, the cation vacancy, in a random $In0.5Ga0.5As$ alloy. We found that compared to the non-random crystal, random crystalline alloys introduce variations for the defect/dopant vibrational mode, whose frequencies and localization character now depend sensitively on the local environment. Specifically, the phonon modes of dopant/defect in a ternary alloy, here InGaAs, tend to occur at an intermediate frequency between those in the binary compounds (GaAs and InAs), and the local phonon density of states typically has more satellite peaks. In addition, the local vibrational modes of $SiIII$ under all the local environments studied occur at frequencies higher than those of host atoms, ranging from 10.9 to 12.1 THz. In contrast, $SeAs$ exhibits resonant vibrational modes for some local configurations, and the most intense peaks occur at very low frequencies ranging from 4.2 to 4.5 THz. For other local configurations, the vibrational modes occur at higher frequencies (6.0 to 7.8 THz) than those of the host atoms, assuming a more localized character. Compared to dopants like Si and Se, the phonon modes of the nearest neighbors of a cation vacancy are much less affected by the defect center. More than half of the vacancy’s first nearest neighbors show a stronger global vibrational character than local ones, while all the first nearest neighbors of Si and Se strongly respond to the presence of the defect and show dominant local vibrations. Our study reveals the significant influence of local environments on the vibrational signature of impurities in random alloys, and our method opens the door to investigate such phenomena in multicomponent crystalline materials.

## SUPPLEMENTARY MATERIAL

See the supplementary material for local atomic environments of group III and V sites in a 216-atom quasirandom supercell of $In0.5Ga0.5As$ and the method validation results of neglecting the interactions between “non-overlapping” atoms and the defect.

## ACKNOWLEDGMENTS

The authors thank Dr. Matthias Polozcek, Dr. Marco Arrigoni, and Dr. Xiaokun Gu for fruitful discussions. Part of the code used to generate the defect force constant matrix was inspired by the code by Dr. Robin Stern and the LAMMPS phonon code written by Dr. Ling-Ti Kong. The calculations were performed with the resources provided by the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by the National Science Foundation (NSF) Grant No. TG-DMR160062 and by the Maryland Advanced Research Computing Center (MARCC) funded by the State of Maryland.

## AUTHORS' CONTRIBUTIONS

H.J. and J.W. contributed equally to this work.

## REFERENCES

*Semiconductors and Semimetals*, edited by E. R. Weber (Academic, New York, 1993), Vol. 38, p. 59.

*Cliques, Coloring and Satisfiability: 2nd DIMACS Implementation Challenge*, edited by D. S. Johnson and M. A. Trick (American Mathematical Society, 1993), Vol. 26, pp. 533–558.

*Semiconductors and Semimetals*, edited by M. Stavola (Academic, New York, 1999), Vol. 51B, p. 153.