#### Abstract

Structural and elastic properties of , a novel semiconductor alloy, are studied from the first principles in both zinc-blende and wurtzite structures. Performances of the* finite difference* (FD) method and the* density functional perturbation theory* (DFPT) are tested and compared. Both of these methods are applied to two different approaches of alloy simulation, a supercell of 16 and 32 atoms (for zinc-blende and wurtzite structures, resp.) and the* alchemical mixing* (AM) method, where the pseudopotentials are mixed in an appropriate way to form an alloy. All elastic properties, including the elastic tensors, elastic moduli, Poisson’s ratio, B/G, and relaxation coefficient, as well as lattice parameters are calculated using all said methods. Conclusions about the use of the approaches investigated in this paper and about their performance are drawn. In addition, in both crystal structures, the band gap is studied in the whole composition range using the MBJLDA functional. The band gap bowings are unusually high, which confirms earlier reports.

#### 1. Introduction

Semiconductor technology plays a crucial role in optoelectronics. Nitride semiconductor alloys, for example, allow constructing optoelectronic devices for wavelengths varying from ultraviolet (AlN [1]) through the visible spectrum (GaN [2]) to infrared (InN [3]). Since the application and performance of such devices depend mostly on the band gap of used alloy, the main focus of researchers is placed on studying their electronic structure, while other properties such as elastic constants are usually neglected. However, designing a semiconductor based optoelectronic device requires knowledge of not only the electronic band structure of the material but also its lattice parameters and elastic properties. Because the layers are usually grown using epitaxial methods, the substrates crystal structure induces strains on the grown layer, modifying its electronic properties. The knowledge of elastic constants and other related elastic properties allows the calculation of these changes, thus making the design much more accurate.

This paper reports on the structural parameters, elastic constants, and other elastic properties such as elastic moduli, Poisson’s ratio, B/G, and relaxation coefficient for a novel semiconductor alloy in both zinc-blende and wurtzite structures. All mentioned parameters are calculated in a number of different approaches including supercell (SC) calculations, the alchemical mixing of pseudopotentials (AM), the finite difference (FD) method of calculating the elastic constants, and the density functional perturbation theory (DFPT). All the different approaches are then compared and conclusions are drawn. In addition, the band gap of the alloy in the whole composition range is calculated using the MBJLDA functional [4], which is known for excellent performance in III-V semiconductors for obtaining the band structure [5].

#### 2. Computational Methods

In the beginning, the structure was relaxed in order to obtain the equilibrium lattice constants and atomic positions. This was performed via energy/forces relaxation with the Perdew and Zunger [6] parametrization of the LDA exchange-correlation energy functional. Two approaches of alloy simulation have been employed: a supercell of 16 and 32 atoms (for zinc-blende and wurtzite structure, resp.) and the alchemical mixing of pseudopotentials [7, 8].

Then, the elastic constants were calculated. Since this paper is focused on the analysis of computational methods, there were a few different approaches to calculating the elastic properties of the alloy. First, we used the so-called* finite difference* (FD) method, where the elastic constants are calculated directly from the stress-strain relations. For the ZB structure, two strains (presented as vectors in Voigt notation) have been applied: a tensile strain and a shear strain . For the WZ structure, one more tensile strain was needed. To obtain the most accurate results, we used six values of for each kind of strain, three positive and three negative. The obtained values of Hellmann-Feynman stress were then interpolated to to produce the values of elastic constants at equilibrium state, corresponding to infinitesimal strains.

The second approach to calculating the elastic constants was the use of the density functional perturbation theory (DFPT) as implemented in the ABINIT code [9]. The calculations started from previously optimized cell geometry with the atoms relaxed to their equilibrium positions. Then, the perturbations were introduced and first-order wave functions were calculated. These were used to obtain the second derivatives of the total energy with respect to those perturbations. The perturbations needed for elastic constant calculations were strains and atomic displacement, both in all nonequivalent configurations. This allowed obtaining the so-called* clamped-ion* elastic tensor, which comes from the second derivatives of the total energy with respect to strain components where atoms are displaced proportionally, and* relaxed-ion* elastic tensor which comes from the second derivatives with respect to one strain component and one component of each internal atomic displacement, the* internal strain* [9].

There were also two different approaches to alloy simulation. The first one was the use of a multiplication of a primitive cell and a supercell (SC) containing 16 or 32 atoms (for zinc-blende and wurtzite structure, resp.). The calculations in the SC approach were performed for 8 compositions, in all of which the metal atoms were distributed as uniformly as possible (the average distance to a respective number of nearest neighbors of the same kind was the biggest). The second one was the alchemical mixing (AM) which was, due to its similarity to the virtual crystal approximation (VCA), called the VCA in our previous papers [21]. A detailed description of the alchemical mixing method can be found in the computational details section of [7].

Pseudopotentials have been generated with the use of the Opium package [22]. The Perdew and Zunger form of the LDA [6] for the exchange-correlation functional was employed in the scalar relativistic mode. The valence states have been chosen as for the N atom and for the Al and P atoms. All the calculations were performed with the use of the ABINIT package [23]. The 8 × 8 × 8 and 4 × 4 × 4 -point grids [24] (for primitive cell and SC, resp.) have been applied with standard shifts for ZB and WZ structures. A tolerance of Ha in energy difference was used for sufficient convergence in FD calculations and geometry optimization, and a tolerance of Ha for potential residual was used in DFPT calculations. The pseudopotentials employed in the calculations required 60 Ha energy cut-off for the plane wave basis set.

The quantities related to elastic constants, bulk (B), shear (G), and tensile (E) moduli, Poisson’s ratio, B/G, and relaxation coefficients have been calculated within the Voigt-Reuss-Hill approximation [25].

For the band structure, we used the LDA optimized geometries with the MBJLDA [4] functional. This functional has been proven numerous times to be one of the most accurate ones for band structure calculations, matching even the state-of-the-art hybrid functionals [5, 26]. However, while hybrid functionals are hundreds of times more computationally expensive, the MBJLDA does not require much more than the standard LDA/GGA approach. In these calculations, also both approaches of alloy simulations have been tested: the supercell and alchemical mixing.

#### 3. Results and Discussion

The calculated values of equilibrium lattice parameters (ZB), ratio, and internal parameter (WZ) for parent AlN and AlP compounds are presented in Table 1. The overall agreement with the reference data suggests a satisfying accuracy of pseudopotentials used here so all subsequent predictions of elastic constants and related values should be correct. Since the main goal of this study was to investigate the performance of various calculation methods, further discussion of the results for parent compounds, for example, the elastic constants, has been omitted. The tables containing all calculated parameters would be huge and confusing so all results are presented in a graphical form, where they can be easily read (Figures 2, 3, and 4). Reference data for parent compounds are well known and have been reported many times, for example, in [27, 28] for AlN and [29, 30] for AlP.

As has been stated in our previous studies [21], an alloy simulated by the alchemical mixing method exhibits a bowing of lattice constants, as presented in Figure 1 (see also Figure in [21]), which is due to the lack of the local relaxation of atomic positions, since the mixing of pseudopotentials is applied to all of the atoms as if they were identical and they always occupy the characteristic of particular lattice high symmetry positions. This bowing in lattice constant manifests itself on all further structural and elastic properties where forces and stresses acting on atoms have to be taken into account. It can be seen in Figures 2 and 3, where the elastic constants are presented, that the calculations that employ the alchemical mixing show an abnormal bowing, much bigger than that appearing when supercells are used. As was said earlier, this excess of bowing in elastic constants is caused mainly by the improper geometry optimization and lack of local relaxation of atomic positions, the same as in the case of lattice constants.

The supercell approach to simulating an alloy should in principle give much better results for both lattice parameters and elastic constants. The calculations reported in this paper confirm this fact, which can be easily seen in Figures 1, 2, and 3. The lattice constants are nearly perfectly linear, thus as expected, obeying linear Vegard’s law. The elastic constants exhibit much smaller bowing than that in the AM. However, the values obtained with the supercells exhibit some irregularities and deviations from their expected, parabolic character. This is justified and easily understandable, since total relaxation and proper convergence in a bigger cell with lower symmetry are much more difficult to reach, and even after a successful calculation, there may be some residual stresses or forces acting on some of the atoms. The deviations from regular behavior are negligible when calculating elastic constants or lattice parameters but are enhanced when calculating the related quantities like bulk and shear moduli or Poisson’s ratio, especially while using the Voigt-Reuss-Hill approximation.

The finite difference (FD) method of calculating the elastic constants was used in our previous work [7] and gave excellent results. In principle, it is relatively easy and is expected to produce accurate results. Each strain requires an independent calculation, and the use of this method even on big supercells (32 atoms for wurtzite structure) is rather straightforward. The detailed procedure used in this approach is described in previous section, and the results obtained are consistent and reliable.

The density functional perturbation theory (DFPT), although much more complicated mathematically, should give the most precise and accurate results [9]. Unfortunately, it requires more computational resources, since all the perturbations, the values, and configuration of strains and atomic displacements are included in a single run of the program. This fact prevented the calculation of elastic constants with the DFPT on 32-atom wurtzite supercell on the resources that were available to us. However, it may be, at the same time, a benefit. For structures with a symmetry that contains more independent elastic constants and would require more configurations of strains in the FD method, this approach offers the opportunity to calculate them all at once.

The DFPT also allows obtaining two variants of elastic constants: the so-called* clamped ion* and* relaxed ion*. The clamped-ion results correspond to the results from the finite difference method perfectly. The relaxed-ion results, however, hold a significant discrepancy in some cases. This should be easily explained by the fact that, to reproduce the relaxed-ion results using the FD method, an extreme accuracy in relaxing the atomic positions is needed [9] which is not obtainable in practice. That is also why the relaxed-ion results in all cases exhibit stronger irregularities and deviation from their expected parabolic character. These deviations in both DFPT and AM (SC) result in a very erratic behavior of most of the other related values (Poisson’s ratio, B/G, and relaxation coefficient), when they are calculated from the Voigt-Reuss-Hill approximation. Surprisingly, those quantities, when calculated using the alchemical mixing, behave very regularly and their values are close to those obtained in SC calculations. This can be explained by the fact that those values are calculated by various multiplications and divisions of the elastic constants [25]; thus, the most important is the smoothness of the curves and their consistency between each other. The fact that the elastic constants exhibit unnatural bowing does not play a significant role here, since all of them exhibit similarly deviated character. If the results for parent compounds are consistent for both alloy simulation methods, SC and AM (as they should be if the calculations were performed correctly), the most efficient one to calculate Poisson’s ratio, B/G, and relaxation coefficients seems to be the alchemical mixing, giving the smoothest curve and correct values and being very efficient in calculation.

Figures 5 and 6 show the composition dependence of different band gaps for ZB and WZ structures, respectively. The calculated band gaps agree very well with our previous studies [21] for the ZB structure. However, in this paper, we expanded our research to include both the direct and indirect band gaps and studied both crystal structures, ZB and WZ. The bowing parameters for the ZB are 7.8 eV and 6.7 eV for the direct and indirect -X gap, respectively. The WZ structure shows bowings of 5.64 eV for direct and 4.7 (6.4) eV for indirect -K (-M) gap. The bowings are given for the supercell approach and seem to be unusually large, which we have also found in our previous research [21]. As can be clearly seen, the alchemical mixing significantly deviates from the supercell calculations. As we have shown before, in these calculations, the SC results should present more trustworthy values as well. In our previous research, we have shown that the local geometry optimization and bond lengths between atoms of different species crucially impact the value of the band gaps [7]. The alchemical mixing does not account for that properly which explains the deviation.

#### 4. Conclusions

This paper presents the elastic properties of a novel semiconductor alloy . It also tests the performance of different methods of calculating the elastic constants of an alloy: the finite difference method and the density functional perturbation theory, using two different approaches to simulating an alloy, a supercell and the alchemical mixing of pseudopotentials. Both the FD and DFPT give accurate results when using supercells, and the alchemical mixing gives an abnormal bowing, as expected from our previous works. However, the alchemical mixing approximation turned out to give excellent and consistent results of Poisson’s ratio, B/G, and relaxation coefficient when calculated within the Voigt-Reuss-Hill approximation.

The DFPT enabled the calculations of the so-called clamped-ion and relaxed-ion cases, while the finite difference method allowed calculating only one kind of elastic constants, which reproduced the clamped-ion DFPT values perfectly. It has been expected, since the relaxation of the atomic positions to the accuracy needed in the relaxed-ion case is virtually impossible in FD calculations.

The DFPT calculations on a 32-atom supercell appeared to exceed the efficiency of computational resources we had an access to, while the FD calculations converged and finished without any problems. That may point at the FD as a method of choice, when dealing with larger systems and bigger atoms. On the other hand, the DFPT method provides the possibility of calculating all the elastic constants in one run, which is impossible in the FD method. This feature is very beneficial in investigations of low symmetry systems where there are more independent elastic constants.

The band gap behavior with the composition exhibits a very big band gap bowing within the supercell approach. It is worth noticing that in the WZ structure the system starts as an indirect band gap material and with the increasing phosphorus composition gradually changes to a direct one. The alchemical mixing, as expected, performs very poorly for both structures.

#### Competing Interests

The authors declare that they have no competing interests.

#### Acknowledgments

The calculations were performed in Wrocław Centre for Networking and Supercomputing. M. P. Polak acknowledges the support within “Diamond Grant” (DI2013 006143) from the MNiSzW.