First-principles study of the phase competition, and mechanical and piezoelectric properties of pseudo-binary (SiC)_1−x(AlN)_x alloy

Abstract

The ongoing search for new piezoelectric materials offering adequate balance between piezoelectric response and other application-relevant properties has led to the investigation of various alloy systems. In this work we study the alloy of the widely used AlN with SiC for their relative abundance, current use in other electronic applications and expected phase competition between wurtzite and other polymorphs, the likes of which has led to some of the most interesting results notably between AlN and ScN. Here the pseudo-binary (SiC)_1−x(AlN)_x alloy is studied from first-principles over the entire composition range. Relevant crystalline phases are identified using the first-principles random structure sampling approach which, in accordance with previous bulk experiments, finds wurtzite, zincblende and rhombohedral phases to be the only statistically relevant phases of the alloy. Further study of these phases is done through special quasirandom structures (SQSs) and, in the case of the wurtzite phase, predictions of the stiffness, piezoelectric and dielectric tensors. Analysis of these tensors is done through the scope of a bulk acoustic wave (BAW) filter application, where trends and trade-offs between the c-axis acoustic velocity and piezoelectric response enable identification of relevant compositions.

---

1 Introduction

Alloying has been a common procedure for improving piezoelectric materials. Notable examples are the AlN-based alloys.¹ The aim has mostly been towards increasing the relatively low piezoelectric constant of the wurtzite AlN phase (d₃₃ ≈ 5 pC N⁻¹) achieved through alloying with scandium.^2,3 However, the scarcity of Sc has greatly increased the demand for alternatives, preferably more affordable and easily implemented with current fabrication techniques such as complementary metal-oxide semiconductors (CMOS).

In this work the choice of exploring alloys between AlN and SiC is motivated by several factors. First, this alloy has already been synthesized both in bulk and thin films.^4–9 SiC is also readily available, relatively inexpensive, already present in micro-electro-mechanical systems where it is used for its resistance to harsh conditions specifically at high temperatures. It is also known for its polymorphism and the phase competition between different crystalline structures. Phase competition in AlN-based alloys has been synonymous with significant piezoelectric property evolution and is thought to be at the origin of the surprising results of the AlScN alloy.¹⁰ Finally, while the piezoelectric response usually takes the focus of materials development, other properties such as the speed of sound are of high importance in applications such as resonators.^11,12 The choice of SiC is thus also driven by its high acoustic speed¹³ and its potentially beneficial transfer into the alloy.

This work presents a study of the pseudo-binary alloy (SiC)_1−x(AlN)_x through several steps in the following order: (i) identification of relevant phases and phase competition analysis, (ii) mechanical property calculations for the most relevant phase, (iii) comparison with experimental data, and finally, (iv) the analysis of the calculated electro-mechanical properties. To identify the statistically relevant phases of the alloy we perform first-principles random structure sampling (FPRSS), which consists of the generation of a large sample of random structures and their subsequent first-principles relaxations to the closest local minima on the potential energy surface (PES) of the system. It has been shown^14,15 that the frequencies of occurrence of various structure-types in the random samplings can be used to predict the likelihood for experimental realizations of identified phases. The FPRSS was successfully used to reproduce known polymorphism in crystalline systems,^14–16 create structural models for amorphous phases^17,18 and to elucidate polymorphism-linked phenomena in thin-films.^19,20

Using FPRSS in this work, for three different alloy compositions (x = 0.25, 0.5, 0.75) we have identified wurtzite, zincblende and rhombohedral as the three energetically and statistically relevant phases, more likely than others to be experimentally realized. The special quasirandom structures (SQS)²¹ used to approximate the random alloy provided the bases for the subsequent phase competition analysis, results of which reveal the wurtzite phase as the thermodynamically most stable phase in accordance with experimental results. Finally, the stiffness, piezoelectric and dielectric tensors are evaluated for the wurtzite phase across the entire composition range. From these tensors we obtain the chemical trends in the c-axis speed of sound and piezoelectric response as well as in the electromechanical coupling coefficient. These properties are of specific interest for BAW filters and allow the identification of the compositions having the optimal combination of the piezoelectric properties and the speed of sound.

2 Methods

2.1 First-principles random structure sampling

The creation of the random structures of a given size (number of atoms) and a given chemical composition is done through the following three steps: (i) a random cell is created with random lattice parameters (a, b, c, α, β, γ), (ii) atoms are distributed within the cell with an algorithm designed for distribution homogeneity, and also favoring cation–anion coordination in relevant cases, and (iii) first-principles relaxation of the structures to the closest local minimum.²²

In the relaxations, a standard first-principles setup is employed. The electron–electron interactions are treated using the PBE²³ exchange–correlation functional and the projector augmented wave (PAW)²⁴ method as implemented in the VASP computer code.²⁵ The automatic generation of the Γ-centered k-point grid is employed with an R_k value of 20. All degrees of freedom, including volume, cell shape and atomic positions, are relaxed using the conjugate gradient algorithm.²⁶ Volume and cell-shape relaxations are restarted at least four times for the real-space grid to be re-created. Structural relaxations are considered converged when the maximal force on any atom gets below 0.02 eV Å⁻¹, the total energy falls below 10⁻⁶ eV and the hydrostatic pressure drops below 0.5 kbar. The obtained structures are understood to be dynamically stable in the majority of cases, as per previous work by Jankousky et al.,¹⁶ with additionally the high symmetry phases of interest in the later part of this work shown to be proper local minima with high transformation barriers.

Following this procedure, the potential energy surface of the system is probed at three different compositions (x = 0.25, 0.5, 0.75) where ∼5000 random structures are generated and relaxed at each composition for a total of ∼15 000 structures. The structures used here contain 24 atoms, per previous FPRSS size dependence studies^15,16,20 which offers a good balance between capturing statistically and energetically relevant phases and the computational cost. The three chosen compositions offer integer number of atoms within each structure while allowing for a symmetric sampling of the compositional range. To obtain accurate statistics (frequencies of occurrence), one must sort relaxed random structures into groups with the same structure-type.

Because of the possible Al–Si disorder present in different structures, the grouping is done based on the symmetry (the space-group number) of the underlying structure and the local first-shell coordination. The symmetry of the underlying (parent) crystal structure is obtained by labeling all cations with the same label and using the spglib library²⁷ with a tolerance of 0.6 Å. The first-shell coordination is obtained using the same spatial tolerance and a count tolerance of 0.1. For the more frequently occurring groups, further confirmation of the symmetry of the structures was done by analyses of simulated powder XRD spectra as well as a visual inspection. Any structures that would indicate a possible phase separation would most often be classified as low symmetry space groups (<10).

2.2 Special quasirandom structures (SQS)

In order to approximate the occupation disorder expected in random alloys, atomic distributions over the cation and anion sublattices are obtained using the Monte Carlo special quasirandom structures (MC-SQS) program as implemented in the ATAT package.²⁸ These structures are created with 128 atoms for the wurtzite or rocksalt phases and 144 atoms for the rhombohedral phase (see Fig. 1). A step in composition of Δx = 0.125 is used for each phase to sample the whole composition range of the alloy. For the Monte Carlo procedure, clusters are considered up to the 5th neighbor for the pair interactions and up to the 3rd neighbor for the three-body interactions.

Fig. 1 An example of a 128 atom (SiC)_0.125(AlN)_0.875 SQS relaxed using the first-principles methods described in Section 2. Light blue, blue, white and brown spheres represent Al, Si, N and C atoms, respectively.

The procedure is carried out four times for each instance with only the best three structures retained. These three structures are then relaxed through a similar DFT procedure as that used for the FPRSS structures with a stricter total energy convergence criterion of 10⁻¹⁰ eV, to facilitate the subsequent mechanical properties calculations.

While the SQSs are created with a specific symmetry, nothing guarantees that it will be conserved throughout their relaxation. In that regard, the symmetry of the SQSs is determined again after relaxation following a procedure similar to that used for FPRSS structures.

2.3 Mechanical properties calculations

Relevant mechanical properties are calculated using VASP with similar parameters than the ones used in the previous Section 2.2. The stiffness tensor is obtained through strain finite differences with the default step of 0.015 Å while the piezoelectric tensor is calculated from density functional perturbation theory (DFPT).

The isotropic Young's modulus and the isotropic Poisson's ratio are used for comparison with experimental measurements on bulk. These isotropic properties are obtained from the stiffness tensor through the Voigt and Reuss averages (upper and lower bounds) of the bulk modulus B_V and B_R and the shear modulus G_V and G_R^29,30 with:

9B_V = C₁₁ + C₂₂ + C₃₃ + 2(C₁₂ + C₂₃ + C₃₁), (1)

15G_V = C₁₁ + C₂₂ + C₃₃ − C₁₂ − C₂₃ − C₃₁ + 3(C₄₄ + C₅₅ + C₆₆), (3)

where, the S_αβ are components of the compliance tensor S = C⁻¹. In the above equations the stiffness tensor is represented by a 6 × 6 matrix through the Voigt notation such that xx → 1, yy → 2, zz → 3, yz → 4, xz → 5, xy → 6. In this representation the stiffness tensor C has components C_αβ(α,β = 1,…,6). The isotropic Young's modulus E and the isotropic Poisson's ratio ν are then obtained from:where and are averages between the high and low bounds named the Voigt–Reuss–Hill averages.

The speed of sound along any unit vector n = (n₁,n₂,n₃) is obtained by solving the Christoffel eigenvalue equation,³¹

where i, j, k, l indices run over the three Cartesian directions, C_iklj is a 4-index element of the stiffness tensor, v = |v| is the norm of the phase velocity of the monochromatic plane wave polarized along s = (s_x, s_y, s_z) direction. Solving this eigenvalue equation offers three solutions two of them corresponding to transversal modes while the third corresponds to the longitudinal one; the eigenvectors corresponding to the polarization direction and the eigenvalue to the associated velocities.

The piezoelectric stress tensor e is a 3 × 6 matrix that follows the same Voigt notation as the stiffness tensor.

While calculations performed here return the piezoelectric stress tensor e, the experimentally reported piezoelectric strain tensor d can be obtained through the compliance tensor S.³²

3 Results

3.1 FPRSS identification of relevant phases

Identification of the relevant phases of the (SiC)_1−x(AlN)_x alloy is done by performing FPRSS at three specific compositions x = 0.25, 0.50, 0.75. The resulting space-group-resolved thermodynamic density of states (TDOS) are presented in Fig. 2. Similar structures are grouped together and referred to by the space group number of the underlying/parent structure as described in Section 2.1.

Fig. 2 Space group resolved thermodynamic density of states of the (SiC)_1−x(AlN)_x alloy for x = 0.25, 0.50 and 0.75. Each plot is derived from ∼5000 structures obtained from random structure sampling. Groups of structures with the same underlying/parent structure (see Section 2.1 for details) are designated by their parent-structure space group numbers and are ordered in ascending (minimal) total energy in each legend. For clarity low symmetry structures (space group # < 10) have been grouped together. For better visualization, a Gaussian broadening of 0.005 eV has been applied to the data.

When analyzing the TDOSs in the three different compositions, one can observe two distinct regions. A low energy region, from 0.0 eV to ∼0.3 eV per atom, where the high symmetry structures dominate, and a high energy region, centered around ∼0.5–0.6 eV per atom, dominated by the low symmetry structures. One can associate the high symmetry structures with different crystalline phases of the two parent compounds, while the low symmetry structures (s.g. # < 10) collectively can be thought of as representing the amorphous phase of the system.^17,18

The noticeable energy difference between the high-symmetry structures and the amorphous phase is an indication of the system's preference to crystallize rather than grow amorphous. This conclusion is supported by previous studies of both crystalline phases²⁰ and amorphous systems,^18,33 where the presence of relatively large energy separation between low-energy, high-symmetry (ordered) phases and low-symmetry states has been associated with a system favoring high symmetry crystalline phases. This can be further supported by estimating the thermal energy per atom at the upper bound of the melting point of the alloy, T_M(x) = (1 − x) × T_M(SiC) + x × T_M(AlN) denoted in Fig. 2 by the vertical dashed lines. Considering the system under some processing conditions at temperatures close to T_M, structures far above the dashed lines (the majority of the low-symmetry ones) would not be easily accessible. This is in contrast to the high-symmetry structures presenting energies per atom lower than the expected kinetic energy of the system at that point, thus further demonstrating the preference of the system to crystallize.

Upon closer inspection of the low energy region of each plot in Fig. 2 the lowest energy structure groups are consistently those with space group numbers #186, #216 and #160, that is, those with the underlying wurtzite, zincblende, and rhombohedral symmetry, respectively.

In addition to their low energies, as shown in Table 1, these three groups are always the most represented (the most frequently occurring) in the FPRSS, except for the #186 structures at 25% AlN. The combination of these two arguments leads to the conclusion that experimental growth of the three different compositions studied here will most likely crystallize in one of the three phases identified (hexagonal #186, cubic #216, rhombohedral #160). This conclusion is in accordance with previous experimental results in growing both bulk and films of (SiC)_1−x(AlN)_x alloys.^4,5,9

Table 1 Frequency of occurrence of high symmetry structures in the random structure sampling over different %AlN stoichiometry. High-symmetry structure groups not present in this table appear significantly less than the 3 groups presented here

	25% AlN	50% AlN	75% AlN
#186 (Wurtzite)	1	10	35
#216 (Zincblende)	11	10	20
#160 (Rhombo.)	9	11	33

---

3.2 Mixing thermodynamics from SQS

Once the relevant phases are identified, the special quasirandom structures (SQSs) are generated for all three phases over a spectrum of compositions. Using SQS to model properties of a random alloy is more practical than the averaging over the structures obtained by the FPRSS. While the two methods are shown to produce very similar results,³⁴ converging ensemble averages is much more computationally demanding in comparison to the phase identification. Hence in this work the FPRSS is used to identify the relevant phases and the tendency for crystallization, whereas the SQSs are employed to evaluate relevant properties of a random alloy.

Large supercells reproducing the occupation-disorder of the (SiC)_1−x(AlN)_x alloys are created through the SQS procedure as detailed in Section 2.2. The results presented in this section are averages over the three SQSs obtained at each instance. Unless indicated otherwise, the deviation between the three SQSs was found to be negligible, showcasing a good overall approximation of a random alloy.

Fig. 3 shows density (upper) and mixing enthalpy (lower) results. In the former, one can note a similar bowing of all three phases. While there is no physical reason to expect a linear evolution of the density or the lattice parameters, it is always interesting to notice departure from the empirical Vegard's law. In experiments, downward bowing of the density or lattice parameters appears in some instances^6,9 and not in others.⁸ The values of the density however, when compared to our results, do not go beyond 2%, which when brought back to linear dimensions (1%) is within the expected error from the PBE functional (1–2%).³⁵

Fig. 3 Density and mixing enthalpy evolution over the entire range of the pseudo-binary (SiC)_1−x(AlN)_x alloy using 128 atoms SQSs in both wurtzite (orange) and zincblende (blue) phases and 144 atoms SQSs for the 9R phase.

In the lower plot, the mixing enthalpies are calculated relative to the wurtzite and the zincblende ground states of AlN and SiC, respectively. Here the wurtzite phase is favored at most compositions, with significant competition between phases only appearing for SiC-rich compositions x ≤ 0.25. This corroborates experimental results where the wurtzite phase of the (SiC)_1−x(AlN)_x alloy is attainable at virtually any composition and dominates at most of them (x > 0.2).^4–9 These experimental results also show the importance of the growth parameters on the resulting alloy phase, supporting the predicted relatively tight thermodynamic competition between phases.

3.3 Mechanical properties of the wurtzite phase

Compared to the other two phases of the alloy, the wurtzite phase not only prevails thermodynamically but also presents important piezoelectric properties, the likes of which have been of great interest in multiple applications for both pure compounds (AlN,^36–38 SiC^39,40). In that regard, the rest of this work will focus solely on the wurtzite phase and its mechanical properties.

X-ray diffraction (XRD) data are simulated for all wurtzite SQS structures. The (100) and (002) peak positions are used to derive the wurtzite lattice parameters following the procedure from the experimental data^6,8 used for comparison in the left panels of Fig. 4. The lattice parameters of the SQSs behave similarly to experimental results throughout the different compositions. The apparent absolute value discrepancy is still within the expected error from the PBE functional and is consistent with the linear dimension uncertainty obtained from the density (1%). Similarly, it is interesting to correlate the bowing of the density from Fig. 3 to the uncoupled bowing of the lattice parameters seen here.

Fig. 4 Evolution of lattice parameters a and c as well as the isotropic Young's modulus and Poisson's ratio of the wurtzite phase of the pseudo-binary (SiC)_1−x(AlN)_x alloy over compositions. Experimental data taken from Lubis et al.⁶ and Rafaniello et al.⁸

The isotropic Young's modulus and Poisson's ratio are presented in the right panels of Fig. 4 alongside experimental results from ref. 6 and 8. These quantities are obtained from the stiffness tensor following the procedure presented in Section 2. As observed in the plots, both the Poisson's ratio and Young's modulus show very good agreement between the computational results and experiments, where once again, the non-linear trends are reproduced between the two. Quantitatively the Poisson's ratio data are in excellent concordance, while the Young's modulus results show more variations with a relative difference of at most 8% which is still within the expected error of 11% for similar mechanical material moduli obtained from the PBE functional.³⁵

Overall, the agreement between the computed and experimental data shown in Fig. 4 demonstrates the reliability of the theoretical methodology used to model the pseudo-binary wurtzite alloy (SiC)_1−x(AlN)_x and warrants confidence in the subsequent quantities presented in this work.

As a platform for discussion, further analyses and calculations are aimed towards properties relevant to the use of piezoelectric materials in bulk acoustic wave (BAW) filters. In that context we present in Fig. 5 results for the longitudinal speed of sound v_l,z (green), the piezoelectric strain constant d₃₃ (blue) and the electromechanical coupling coefficient k₃₃², all of which are taken along the c direction of the wurtzite lattice, see Section 2 for definitions. The electromechanical coupling coefficient presented here is defined as^2,41

where ε₃₃ is a component of the total dielectric tensor ε while C₃₃ and e₃₃ are components respectively of the stiffness and piezoelectric stress tensor. In the plot each point corresponds to an average between the calculated values of 3 different SQSs generated at that chemical composition while the error bars show maximum and minimum values of the sample. For pure compounds presenting no occupational disorder, a single structure is used in those cases and one can already notice the adequate reproduction of the presented properties for AlN (v_l,33 ≃ 10 500 m s⁻¹ (ref. 3) and d₃₃ ≃ 5.5 pC N⁻¹ (ref. 42)).

Fig. 5 Composition dependence of the longitudinal speed of sound along the c-axis (in green) and of the piezoelectric strain constant d₃₃ as defined in Section 2.3 (in blue) for the wurtzite pseudo-binary alloy (SiC)_1−x(AlN)_x. The inset shows the evolution of the electromechanical coupling coefficient k₃₃². Error bars show maximal and minimal values of the sample while markers are averages. Complete stiffness, piezoelectric stress and dielectric tensors can be found in the SI.

4 Discussion

The properties presented in Fig. 5 are of particular interest for applications in bulk acoustic wave (BAW) filters, which initially gained interest in the early 2000s,⁴³ where the acoustic resonators were based on an AlN piezoelectric layer. Improvements to these filters were achieved in part through modifications of the piezoelectric material, a domain that saw a significant push from the 2010s onward with the introduction of the AlScN alloy⁴² and its large increase in the piezoelectric response.

Improvement of the piezoelectric layer has generally two routes, the first one is to improve the piezoelectric response of the material which will improve several parameters such as broadening the bandwidth or increasing the efficiency of the filter, as can be seen from eqn (9). On the other hand one can improve the speed of sound in the piezoelectric material which will lead to higher working frequencies for the filter, allow for thicker films and lower the acoustic loss. These two approaches are usually understood as mutually exclusive, since the piezoelectricity is understood as benefiting from a more “malleable” material where atoms can more easily move and create dipoles⁴⁴ while on the other hand acoustic velocities are understood to take advantage of more “stiff” materials with more correlated interatomic displacements. The equations presented in Section 2 can further be approximated to also reveal this inverse relationship where the speed of sound can be and the piezoelectric strain response d₃₃ ≈ e₃₃/C₃₃. Here the presence of the stiffness tensor coefficient C₃₃ on different sides of the fraction clearly shows the opposition between these two properties.

The curves in Fig. 5 support this understanding, where d₃₃ and v_l,33 clearly have an inverse relationship. The choice of the composition will then be dependent on the desired improvement to AlN. A small increase in the piezoelectric response can be obtained by introducing small amounts of SiC in AlN (x = 0.875). Conversely, going to x = 0.500, the speed of sound increases by about 5% but that lowers the piezoelectric response by about 40%. Simultaneously, the efficiency of the material measured by k₃₃² diminishes for x < 0.875 compositions. As indicated the variations in these different properties will open choices in optimizations of the material stoichiometry for applications but will largely be limited by the evolution of the k₃₃² parameter.

Trying to elucidate a microscopic origin to these trends, Fig. 6 presents the evolution of the wurtzite internal parameter throughout the composition range. From the slight increase in u at x = 0.875 we can understand that small incorporation of SiC in AlN initially distorts the wurtzite. Similar distortions have been reported in other AlN-alloys^10,45 and related to greater atomic displacement and increase of the piezoelectricity. However the distortion is here on a smaller scale and further incorporation of SiC lowers the u value bringing the wurtzite structure closer to its pristine state. This lowering of the distortion accompanied with the lower bond polarity of SiC can explain the fall in the piezoelectric properties of the alloy.

Fig. 6 Evolution of the wurtzite internal parameter u of the pseudo-binary (SiC)_1−x(AlN)_x alloy over compositions.

5 Conclusion

In this work we presented a first-principles study of the pseudo-binary (SiC)_1−x(AlN)_x alloy. First we performed first-principles random structure sampling (FPRSS) at select compositions in order to get an understanding of the relevant phases of the alloy. Wurtzite, zincblende and rhombohedral (9R) were found to be the three phases of importance as in previous bulk material results. We continued by creating large cells (128 or 144 atoms) through the use of SQSs which after DFT-PBE relaxation offered an analysis of the enthalpy of mixing competition between them. The wurtzite phase was shown to either be prevailing (x > 0.25) or not significantly higher than the other phases. Stiffness, piezoelectric response and dielectric tensors of all the wurtzite structures were then calculated in order to obtain properties relevant to the BAW filter application. The data show the known inverse relation between the piezoelectric strain d₃₃ and the acoustic velocity v_l,z with an interesting peak of the former for x = 0.875 and an overall downward bowing of the latter, allowing for different optimization of the piezoelectric material with regards to the BAW filter application. The complete tensor data is also made available, enabling further property calculation for other MEMS applications or as parameters for lower theory level simulations.

Conflicts of interest

There are no conflicts to declare.

Data availability

The data that support the findings of this study are available in the published manuscript and its supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d6tc00220j.

Moreover, additional information and results are available from the authors upon reasonable request.

Acknowledgements

This work was supported by the NSF PFI program under Award No. 2234617.

Notes and references

1. J. Startt, M. Quazi, P. Sharma, I. Vazquez, A. Poudyal, N. Jackson and R. Dingreville, Adv. Electron. Mater., 2023, 9, 2201187.
2. G. Wingqvist, F. Tasnádi, A. Zukauskaite, J. Birch, H. Arwin and L. Hultman, Appl. Phys. Lett., 2010, 97, 112902.
3. O. Ambacher, S. Mihalic, M. Yassine, A. Yassine, N. Afshar and B. Christian, J. Appl. Phys., 2023, 134, 160702.
4. M. Landon and F. Thevenot, Ceram. Int., 1991, 17, 97–110.
5. A. Zangvil and R. Ruh, J. Am. Ceram. Soc., 1988, 71, 884–890.
6. A. H. Lubis, N. L. Hecht, G. A. Graves and R. Ruh, J. Am. Ceram. Soc., 1999, 82, 2481–2489.
7. J.-F. Li and R. Watanabe, J. Mater. Sci., 1991, 26, 4813–4817.
8. W. Rafaniello, K. Cho and A. V. Virkar, J. Mater. Sci., 1981, 16, 3479–3488.
9. Z. Gu, L. Du, J. H. Edgar, E. A. Payzant, L. Walker, R. Liu and M. H. Engelhard, MRS Internet J. Nitride Semicond. Res., 2005, 10, e5.
10. F. Tasnádi, B. Alling, C. Höglund, G. Wingqvist, J. Birch, L. Hultman and I. A. Abrikosov, Phys. Rev. Lett., 2010, 104, 137601.
11. Y. Liu, Y. Cai, Y. Zhang, A. Tovstopyat, S. Liu and C. Sun, Micromachines, 2020, 11, 630.
12. T. Pensala, Thin film bulk acoustic wave devices: Performance optimization and modeling, PhD thesis, VTT, 2011.
13. K. Sarasamak, S. Limpijumnong and W. R. L. Lambrecht, Phys. Rev. B: Condens. Matter Mater. Phys., 2010, 82, 035201.
14. E. B. Jones and V. Stevanović, Phys. Rev. B, 2017, 96, 184101.
15. V. Stevanović, Phys. Rev. Lett., 2016, 116, 075503.
16. M. Jankousky, E. M. Garrity and V. Stevanović, Phys. Rev. Mater., 2023, 7, 053606.
17. E. Jones and V. Stevanović, npj Comput. Mater., 2020, 6, 56.
18. L. Wolf, A. Novick and V. Stevanović, J. Appl. Phys., 2025, 137, 095101.
19. R. Woods-Robinson, V. Stevanović, S. Lany, K. N. Heinselman, M. K. Horton, K. A. Persson and A. Zakutayev, Phys. Rev. Mater., 2022, 6, 043804.
20. A. Zakutayev, M. Jankousky, L. Wolf, Y. Feng, C. L. Rom, S. R. Bauers, O. Borkiewicz, D. A. LaVan, R. W. Smaha and V. Stevanović, Nat. Synth., 2024, 3, 1471–1480.
21. A. Zunger, S.-H. Wei, L. G. Ferreira and J. E. Bernard, Phys. Rev. Lett., 1990, 65, 353–356.
22. V. Stevanović, S. Lany, X. Zhang and A. Zunger, Phys. Rev. B: Condens. Matter Mater. Phys., 2012, 85, 115104.
23. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
24. P. E. Blöchl, Phys. Rev. B: Condens. Matter Mater. Phys., 1994, 50, 17953–17979.
25. G. Kresse and D. Joubert, Phys. Rev. B: Condens. Matter Mater. Phys., 1999, 59, 1758–1775.
26. M. P. Teter, M. C. Payne and D. C. Allan, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 40, 12255–12263.
27. A. Togo, K. Shinohara and I. Tanaka, Sci. Technol. Adv. Mater.: Methods, 2024, 4, 2384822.
28. A. van de Walle, P. Tiwary, M. de Jong, D. Olmsted, M. Asta, A. Dick, D. Shin, Y. Wang, L.-Q. Chen and Z.-K. Liu, Calphad, 2013, 42, 13–18.
29. R. Hill, Proc. Phys. Soc., London, Sect. A, 1952, 65, 349.
30. M. de Jong, W. Chen, T. Angsten, A. Jain, R. Notestine, A. Gamst, M. Sluiter, C. Krishna Ande, S. van der Zwaag and J. J. Plata, et al., Sci. Data, 2015, 2, 150009.
31. J. W. Jaeken and S. Cottenier, Comput. Phys. Commun., 2016, 207, 445–451.
32. J. F. Nye, Physical properties of crystals: their representation by tensors and matrices, Oxford University Press, 1985.
33. O. V. Pshyk, S. Zhuk, J. Patidar, A. Wieczorek, A. Sharma, J. Michler, C. Cancellieri, V. Stevanović and S. Siol, Adv. Mater., 2025, 37, 2501074.
34. A. Novick, Q. Nguyen, R. Garnett, E. Toberer and V. Stevanović, Phys. Rev. Mater., 2023, 7, 063801.
35. G.-X. Zhang, A. M. Reilly, A. Tkatchenko and M. Scheffler, New J. Phys., 2018, 20, 063020.
36. H. Yang, J. Sun, H. Wang, H. Li and B. Yang, J. Alloys Compd., 2024, 989, 174330.
37. H. P. Nguyen, Light: Sci. Appl., 2022, 11, 164.
38. R. Ding, W. Xuan, S. Dong, B. Zhang, F. Gao, G. Liu, Z. Zhang, H. Jin and J. Luo, Nanomaterials, 2022, 12, 3082.
39. H. Morkoç, S. Strite, G. B. Gao, M. E. Lin, B. Sverdlov and M. Burns, J. Appl. Phys., 1994, 76, 1363–1398.
40. Y. Zhou, World J. Text. Eng. Technol., 2024, 12, 1034–1045.
41. K. Hirata, Y. Mori, H. Yamada, M. Uehara, S. A. Anggraini and M. Akiyama, Materials, 2021, 14, 309.
42. M. Akiyama, T. Kamohara, K. Kano, A. Teshigahara, Y. Takeuchi and N. Kawahara, Adv. Mater., 2009, 21, 593–596.
43. R. Ruby, P. Bradley, J. Larson and Y. Oshmyansky, Electron. Lett., 1999, 35, 794–795.
44. R. M. Martin, Phys. Rev. B, 1972, 5, 1607–1613.
45. S. Manna, K. R. Talley, P. Gorai, J. Mangum, A. Zakutayev, G. L. Brennecka, V. Stevanović and C. V. Ciobanu, Phys. Rev. Appl., 2018, 9, 034026.

---