ARTICLE pubs.acs.org/JPCA
Structural, Elastic Constant, and Vibrational Properties of Wurtzite Gallium Nitride: A First-Principles Approach Zahid Usman, Chuanbao Cao,* Waheed S. Khan, Tariq Mahmood, Sajad Hussain, and Ghulam Nabi Research Centre of Materials Science, School of Material Science and Engineering, Beijing Institute of Technology, Beijing 100081, People’s Republic of China. ABSTRACT: PerdewWang proposed generalized gradient approximation (GGA) is used in conjunction with ultrasoft pseudopotential to investigate the structural, elastic constant, and vibrational properties of wurtzite GaN. The equilibrium lattice parameters, axial ratio, internal parameter, bulk modulus, and its pressure derivative are calculated. The effect of pressure on equilibrium lattice parameters, axial ratio, internal parameter (u), relative volume, and bond lengths parallel and perpendicular to the c-axis are discussed. At 52 GPa, the relative volume change is observed to be 17.8%, with an abrupt change in bond length. The calculated elastic constants are used to calculate the shear wave speeds in the [100] and [001] planes. The finite displacement method is employed to calculate phonon frequencies and the phonon density of states. The first- and second-order pressure derivative and volume dependent Gruneisen parameter (γj) of zone-center phonon frequencies are discussed. These phonon calculations calculated at theoretical lattice constants agree well with existing literature.
I. INTRODUCTION GaN semiconductor inherits excellent physical properties like large ionicity, high thermal conductivity, greater mechanical hardness comparable to the diamond, high heat capacity, and a wide band gap of 3.5 eV. These novel physical properties make it a potential candidate for micro- and optoelectronics.1 Its structural properties are remarkable. It exists in two different structures: wurtzite and zinc blende. At ambient conditions, GaN structure is hexagonal, which is most stable structure in contrast to the cubic structure (zinc blende) grown on cubic substrates.2 Both these polytype structures have similar band gaps approximately. Their host atom’s radii have large mismatches (Ga, 1.26 Å; N, 0.75 Å), and their constituent atoms have quite different electronegativities (Ga, 1.8; N, 3.0). In addition to this, the X-ray absorption spectroscopy (XAS) measurements made by Perlin et al.3,4 also proposed its third structure called rock salt GaN, which is produced due to the phase transition of wurtzite GaN at a hydrostatic pressure of 47 GPa. Modern-day electroluminous devices5,6 are built through quantum size and multilayer heterostructures. The manufacturing of such structures is dependent upon a key parameter known as the potential energy step at heterojunctions, which is subsequently dependent upon the electrical and mechanical properties (elastic constants). To understand and improve the efficiency of these devices, knowledge of elastic constant properties is crucial. These elastic constants are usually calculated through ultrasonic measurements and Brillouin zone spectroscopy. The elastic continuum theory7 is employed to determine the deformation as a result of strain states related to pseudomorphic epitaxial film and the free superlattices of thin Ga and N layers, which not only provides a better understanding of changes in band structure but also leads to tailor the other physical properties for above-mentioned novel applications as well. Until now, elastic constant properties are r 2011 American Chemical Society
poorly known with scattered experimental810 and theoretical data.11,12 Therefore, it is much needed to make a systematic study of elastic constants and include its pressure effect to compare with experimentally obtained elastic constants. In spite of technological interest in GaN, still lacking is the knowledge of lattice dynamics theoretically as well as experimentally. The density functional theory (DFT) presents reliable phonon calculations for a wide variety of semiconductors using local density approximation (LDA) and plane wave pseudopotential methods.13 Recently, linear response techniques14,15 have been proposed to calculate the phonon frequencies throughout the Brillouin zone with good accuracy and speed. The results of this technique are comparable to neutron diffraction data, and other physical quantities such as specific heat, electron phonon interactions, thermal expansion, resistivity, and superconductivity can be extracted from these calculations. There exist few studies of optical phonons using Raman spectroscopy or infrared reflectivity under biaxial strains by exploring phonon deformation potential in GaN8,1619 or using the carrier concentration20 or by calculating stress and strains coefficient for high frequency mode E2 of GaN.21 In spite of massive research on GaN, the vibrational properties are still poorly known. For nonpolar zone-center phonons, first principles studies using mixed basis method22 or linear muffin-tin orbital method23,24 are performed along with some studies under hydrostatic pressure.4,25,26 In this paper, we have presented extensive first principles structural, elastic constant and lattice dynamical properties of wurtzite GaN using generalized gradient approximation, in comReceived: July 26, 2011 Revised: November 1, 2011 Published: November 01, 2011 14502
dx.doi.org/10.1021/jp207141k | J. Phys. Chem. A 2011, 115, 14502–14509
The Journal of Physical Chemistry A bination with ultrasoft pseudopotential. Unlike all previous calculations, we have carried out phonon calculations using finite displacement method. In the above-mentioned first principles studies,1327 the force constants are calculated through density functional perturbation theory (DFPT) along with LMTO or LDA potential or linear response calculations. In CASTEP, two important approaches are used to calculate the lattice dynamics calculations, the first one is density functional perturbation theory and the second one is the finite displacement method. DFPT in CASTEP have some restrictions, like the phonon calculations must be performed with fixed occupancies (insulators) and without spin polarization by using only norm conserving pseudopotentials. This approach is faster and accurate. The finite displacement method is based on numerical differentiation of forces when atoms are displaced by a small amount from their equilibrium positions. This may be used in cases where a DFPT calculation is not applicable, for example, in calculating phonon dispersion of metallic or magnetic systems, or when the use of ultrasoft pseudopotential is essential. With generalized gradient approximation and ultrasoft pseudopotential, this is first extensive study of structural, elastic constants, and phonons calculations. Section II presents the elaborate computational details of this study. The structural, elastic constant, and vibrational properties are discussed in section III. Finally, the conclusions are summarized in section IV.
II. METHOD The first principles calculations are subjected to the plane wave pseudopotential approach in context of density functional theory (DFT). The electron exchange correlation energy is taken into account by generalized gradient approximation proposed by PerdewWang (91).27 The Ga and N pseudopotentials are obtained by using ultrasoft pseudopotentials (USP) introduced by Vanderbilt.28 To keep the substantial weight of tight binding orbitals inside the atomic core region, usually a high cutoff energy with large plane-wave basis set is required. This cutoff energy can be dramatically reduced if USP is employed. Such potentials reduce the possible cutoff energy and hence the plane-wave basis set by dividing the electronic density of the Ga and N atoms into a softer part extending throughout the unit cell and a hard part localized inside the core region. Ga 3d electrons are considered as valence electrons in Ga due to their profound effect on band structure. The electronic wave function is expanded by using a plane wave basis set with a cutoff energy of 390 eV. At this energy cutoff, the structure is found to be well optimized with total energy convergence of smaller than 1 meV/Å using the BroydenFletcherGoldfarbSheno (BFGS) method and a residual force of 2 meV/ Å per atom. A MonkhorstPack grid29,30 of 3 3 3 is used to sum over all Brillouin zone. As the hexagonal GaN is a semiconductor, the convergence is obtained only for 6 special k-points. For the wurtzite GaN, total energy minimization is a two-step procedure because of three independent parameters, lattice constants a and c and internal parameter u. First, a set of c/a values is taken near equilibrium value, and internal parameter u is calculated for each c/a value. Then the equilibrium value of c/a is estimated leading to an equilibrium internal parameter u. At the end, the BirchMurnaghan equation of state31 is used to fit the total energy as a function of volume, to get the important ground-state properties like equilibrium volume (V0), bulk modulus B, and its pressure derivative B0 .
ARTICLE
Table 1. Comparison of Structural Parameters of Wurtzite GaN: Lattice Parameters a (Å), Ratio c/a, Internal Parameter (u), Bulk Modulus B0 (GPa), and B0 method
B0
B0
173
4.2
34
188
3.2
35
245
4.0
ref
present study experiment
33
other calculations
a
c/a
u
3.232847
1.617
0.3774
3.189
1.626
36 37
3.143 3.170
1.626 1.620
0.377 0.379
215 207
5.9 4.5
38
3.126
1.638
0.377
190
4.0
39
3.221
1.616
0.376
176.5
4.37
40
3.245
1.634
0.376
172
5.11
The finite displacement method32 is employed herein to calculate the phonon frequencies of wurtzite GaN. In this method, the atoms are displaced from their mean positions and the forces on atoms are calculated by using numerical central difference method. Depending upon the crystal symmetry, CASTEP calculates such an optimum number of displacements itself. These displacements with a very small magnitude give phonon frequencies very close to that of linear response theory. The cutoff energy and convergence parameters for phonon calculations are the same as prescribed above.
III. RESULTS AND DISCUSSION a. Structural Properties. The ground-state properties of wurtzite GaN are specified through the lattice constants a and c, and the internal parameter u. The total energy is minimized with respect to these parameters so that the force Fz acting on the atoms along z-axis and diagonal elements of stress tensors (σx, σy, σz) vanish. These data of total energy minimization against volume are fitted to the following third-order BirchMurnaghan equation of state.31
PðV Þ ¼ 1:5B0 ½ðV0 =V Þ7=3 ðV0 =V Þ5=3 f1 þ 0:75ðB0 4Þ½ðV0 =V Þ2=3 1g B0 is the first pressure derivative of bulk modulus B0, and V0 is the equilibrium volume at zero pressure. The equilibrium lattice parameters, bulk modulus, and its pressure derivative at zero pressure are presented in Table 1. Table 1 shows a comparison of our calculated results with experiments3335 and theoretical results.3640 The calculated lattice constant a is overestimated by 1.3% and axial ratio c/a is estimated to be 1.617, which is underestimated than the expeiment33 and LDA based studies3638,40 but in excellent agreement with ref 39 where GGA approximation is taken into account. Here, the value of internal parameter (u) is close to the ideal value of 0.375 and in excellent agreement with other results.3640 The value of bulk modulus calculated through the BirchMurnaghan equation of state in our case is 173 GPa, which agrees well with experimental results proposed by Xia et al.34 and other well converged first principles calculations.3840 But when compared to the bulk modulus determined experimentally by Ueno et al.,35 our results are far from even the reasonable agreement. With other LDA results,3638 our results are only a few percent underestimated. Because the experimental results show an inconsistency of more than 25%, it is very difficult 14503
dx.doi.org/10.1021/jp207141k |J. Phys. Chem. A 2011, 115, 14502–14509
The Journal of Physical Chemistry A
Figure 1. (a) Lattice parameters of wurtzite GaN under hydrostatic pressure:. The star (triangle) represents our calculations (experiment35). (b) Effect of hydrostatic pressure on axial ratio (c/a). The sphere (star, rectangle) represents the comparison of the present study with ref 4141 and experiment35, respectively.
Figure 2. (a) Relative volume as a function of pressure. The sphere and rectangle describe the comparison of our results with the experimental ones. (b) Change of internal parameter (u) with pressure. The sphere and star represent the comparison of this study with that of ref 41.
to make a clear statement about the accuracy of our results with that of both experimental values. But, as a whole, our results have an error of only a few percent from the other calculations34,3640 with close agreement with one of the experimental values, so we declare that our results agree well within the accuracy of the DFT calculations. The pressure derivative of the bulk modulus is overall in excellent agreement with experimental35 as well as with other theoretical results.3739 Figure 1a,b shows the lattice constant and axial ratio (c/a) as a function of hydrostatic pressure. It is evident from Figure 1a that lattice constants of GaN decrease monotonically with pressure up to 52 GPa. This monotonic decrease has the same trend as that of experimental35 lattice constants under pressure but is slightly overestimated because of the nature of potential applied herein. It is also apparent from Figure 1a that pressure is influencing our calculations significantly, approaching the experimental value at a pressure of 52 GPa. The axial ratio on the other hand is found to be unchanged virtually, in accord with ref 41 and experiment,35 but a slightly underestimated one. Figure 2 manifests the change of relative volume and internal parameter with pressure. The relative volume graph shows a monotonic decrease until 50 GPa, but an abrupt decrease is observed at 52 GPa. This sudden decrease verifies the fact that the structure is under phase transition from wurtzite to rock salt structure at this pressure. With respect to the equilibrium volume, this volume change is estimated to be 17.8%. This volume change is in excellent agreement with that of experimental value predicted by Ueno et al.35 This is further verified
ARTICLE
Figure 3. Effect of hydrostatic pressure on bond lengths parallel and perpendicular to the c-axis.
from the sudden bond length variation under pressure at 52 GPa, as shown in Figure 3. Figure 3 represents the bond lengths parallel and perpendicular to the c axis. The equilibrium bond length is 1.973 Å, in good comparison with the well converged first principles calculations proposed by Palummo et al.42 (1.93 Å). Figure 3 shows that bond lengths parallel and perpendicular to the c-axis decrease with the increase of pressure with a sudden change of bond lengths at pressure 52 GPA. Under external pressure, the bond angles are approximately constant, an effect that is consistent with that of Wagner’s41 theoretical studies, so we did not present its pictorial view here. b. Elastic Constants. The elastic constants are important to characterize the GaN crystal response to the external forces. The GaN is a hexagonal crystal at ambient conditions with five nonzero and independent elastic constants c11, c12, c13, c33, and c44, where the stress tensors σαβ are linearly related with the strain components as follows, 3 2 c11 σxx 7 6 6 σ 6 yy 7 6 c12 7 6 6 6 σzz 7 6c 7 ¼ 6 13 6 6 σyz 7 60 7 6 6 6σ 7 60 4 zx 5 4 σxy 0 2
c12 c11 c13 0 0 0
c13 c13 c33 0 0 0
0 0 0 c44 0 0
0 0 0 0 c44 0
3 2 3 εxx 0 7 6 7 0 7 6 εyy 7 7 6 7 0 7 6 εzz 7 76 7 6 7 0 7 7 6 εyz 7 6 7 0 7 5 4 εzx 5 εxy c66
ð2:1Þ
where 1, 2, 3, 4, 5, 6 denote xx, yy, zz, zx, xy, yz, and x, y, and z represent orthogonal directions, respectively. Elastic energy per unit cell volume is given by eq 2.2 c11 c33 ðεxx 2 þ εyy 2 Þ þ εzz 2 þ c12 εxx εyy 2 2
E¼
þ c13 ðεxx þ εyy Þεzz þ 2c44 ðεyz 2 þ εzx 2 Þ þ 2c66 εxy 2
ð2:2Þ
The elastic constants calculated through elastic constant theory are summarized in Table 2. And the bulk modulus can be calculated through the following expression B¼
c33 ðc11 þ c12 Þ 2ðc13 Þ2 c11 þ c12 4c13 þ 2c33
ð2:3Þ
We have also used elastic constant moduli of cubic GaN as we previously reported48 and calculated wurtzite elastic constants through Martin’s transformation49 using the following formulae: 1 ch11 ¼ ð3c11, c þ 3c12, c 6 þ 6c44, c Þ 3δ2 =ðc11, c c12, c þ c44, c Þ 14504
ð2:4Þ
dx.doi.org/10.1021/jp207141k |J. Phys. Chem. A 2011, 115, 14502–14509
The Journal of Physical Chemistry A
ARTICLE
Table 2. Elastic Constants for Wurtzite GaN method present study experiment
calcutaions
c11
c12
329 (341)
109 (125)
ref
c13
c33
c44
80.0 (56)
357 (410)
8
96
324
175
9
80.4
387
192
10
158
267
195
43
98
389
106 ( 20
398 ( 20
145 ( 20
91 (80)
B 176 (174)
201
12
390 ( 15
44
374
106
70
379
101
180
45 46
296 ( 18
130 ( 10
158 ( 5 104
267 ( 17 414
24 ( 2
195 207
47
367
135
103
405
95
202
37
(396)
(144)
(100)
(392)
(91)
(207)
1 ch12 ¼ ðc11, c þ 5c12, c 2c44, c Þ 6 þ 3δ2 =ðc11, c c12, c þ c44, c Þ
ð2:5Þ
1 ch13 ¼ ð2c11, c þ 4c12, c 4c44, c Þ 6
ð2:6Þ
1 ch33 ¼ ð2c11, c þ 4c12, c þ 8c44, c Þ 6
ð2:7Þ
1 ch44 ¼ ð2c11, c 2c12, c þ 2c44, c Þ 6δ2 =ðc11, c c12, c þ 4c44, c Þ 6
ð2:8Þ
210
of ch11, ch12, and ch33 are improved as compared to those calculated through stress strain theory. These results show a better agreement with other experimental and theoretical results summarized in Table 2. However, the values of c13 and c14 have been decreased as compared to its stressstrain counterparts. When compared with the values obtained through Martin’s transformations,49 our results lie very close to that of Kim et al.37 except for a smaller value of c13. But because of GGA potential, the results are slightly underestimated due to the overestimation of equilibrium lattice parameters. Here we calculate the acoustic wave speeds in the [100] and [001] plane using the elastic constants as presented in Table 2 as proposed by Truel et al.50 In the [100] plane, shear waves have one longitudinal and two transverse components calculated by the following formulas: Longitudinal wave speed (VL) in [100] plane:
where
pffiffiffi 2 2 ðc11, c c12, c 2c44, c Þ δ ¼ 6
105 ( 10
VL ¼ ð2:9Þ
Here the left side of the equations contain hexagonal elastic moduli (with superscript h) calculated through the GaN cubic elastic constants (with subscript c). The second terms in above equations are to compensate the differences between cubic and hexagonal systems. The calculated elastic constants (c13, c33) are very close to most of the experimentally calculated values.810 The elastic constant data are calculated through different experimental methods; for example, Shelag and Savastenko et al.45 estimated through temperature dependence broadening of GaN powder, Takagi et al.44 and Polian et al.12 through Brillouin zone scattering experiments. These sets of data have significant variation, suggesting that these results may be inaccurate. The first principles DFT calculations are mostly performed with local density approximation,37,46,47 so our results for elastic constants and especially c13 are underestimated by a few percent compared with these results, whereas the bulk modulus calculated through eq 2.3 is 176 GPa, close to that of the bulk modulus calculated in different experiments810,43 as well as with some theoretical results37,46,47 by a few percent. It is apparent from Table 2 that our results on a whole are underestimated as compared to the experiments and other ab initio studies. But because the elastic constant data show a significant variation, the elastic constant calculated by USP results seems to be acceptable and unique due to the method applied. The wurtzite elastic constants calculated through above Martin’s equations49 are given in parentheses in Table 2. Our calculated values
1=2 C11 F
ð3:0Þ
Transverse wave speed (VT,[001]) in the [100] plane with polarization along the [001] orientation: 1=2 C44 ð3:1Þ VT, ½001 ¼ F Transverse wave speed (VT,[010]) in the [100] plane with polarization along the [010] orientation: C11 C12 1=2 ð3:2Þ VT, ½010 ¼ 2F Longitudinal wave speed (Vl) in the [001] plane: 1=2 C33 Vl ¼ F Transverse wave speed (Vt) in the [001] plane: 1=2 C44 Vt ¼ F 14505
ð3:3Þ
ð3:4Þ
dx.doi.org/10.1021/jp207141k |J. Phys. Chem. A 2011, 115, 14502–14509
The Journal of Physical Chemistry A
ARTICLE
Table 3. Comparison of Acoustic Wave Speeds (105 cm/s)a plane
present study
other calculations50
VL
7.47 (7.6169)
7.96
VT,[001]
3.902 (3.680)
4.13
VT,[010]
4.322 (4.28)
6.31
Vl
7.776 (8.345)
8.04
Vt
3.90 (3.686)
mode
[100]
[001]
theoretical results quantity
4.13
The values of respective speeds in parentheses are those calculated by exploiting the elastic constants calculated through Martin’s transformation.
These wave speeds are calculated and compared with theoretically calculated wave speeds.50 Table 3 verifies that the wave speeds VL, VT,[001], and Vl calculated through Martin’s transformed elastic constants agree better with that calculated by Truel et al.50 than those calculated through elastic constant theory because of smaller values of C11 and C12. But, as a whole, these values are trustworthy as compared to their counterpart. To estimate the strain produced because of hydrostatic pressure, we have reduced elastic energy per unit cell with one of the strain components (εxx or εzz) for a specific fractional volumetric change (Δv/v0); we get the strain produced relative to the fractional volumetric change (Δv/v0) as c33 c13 Δv εxx ¼ ð2:10Þ c11 þ c12 4c13 þ 2c33 v0 εzz ¼
c11
ð2:11Þ
εxx/(Δv/v0) and εzz/(Δv/v0) for hexagonal GaN are calculated under hydrostatic pressure, and their values are in Table 4 for comparison with existing literature. It is apparent from the above table that both calculated results are in excellent agreement with that of experiment12 as well as with previous first principles results.37,44,45,47 c. Phonon Calculations. The phonon calculations of wurtzite GaN are computed mostly with density functional perturbation theory with local density approximation using frozen core approximation or through linear response calculations. Until now, no one has adopted a finite displacement method to calculate phonon dispersion and phonon density of states along high symmetry directions. We have employed a finite displacement method (supercell method) in combination with generalized gradient approximation using ultrasoft pseudopotential to calculate phonon frequencies. The finite displacement calculations are performed in a supercell, circumventing the sphere described by a cutoff radius of 3 Å, which directly results all the nonzero elements of the force constant matrix. Therefore, only the electron for which q = (0, 0, 0) (Γ point) will match this conditions and all other nonzero q values will be eliminated. The CASTEP shifts each atom slightly and calculates forces on the perturbed system by performing SCF calculations. The central difference method described below is used to calculate force constants in both positive and negative direction such as dFk, a dFk, a þ dFk, a ¼ du 2u
present study
12
experiment
ref 47 ref 37 ref 45 ref 44
εxx/(Δv/v0) 0.3329 (0.333)
0.322
0.336 0.316 0.332 0.323
εzz/(Δv/v0) 0.334 (0.3408)
0.356
0.329 0.368 0.335 0.355
a
a
c11 þ c12 2c13 Δv þ c12 4c13 þ 2c33 v0
Table 4. Comparison of εxx/(Δv/v0) and εzz/(Δv/v0) under Hydrostatic Pressurea
The values in parentheses are obtained from elastic constant measured through Martin’s transformation.
This single pair of calculations will lead to the calculation of an entire row of a dynamical matrix. Therefore, first we perform minimal perturbations and then extend these calculations to compute the entire dynamical matrix by assuming space charge symmetry. These calculations are efficient and very fast, like linear response calculations. But this method has limitation in calculation of nonanalytic terms responsible for LO-TO splitting in contrast to the linear response calculations only. In this case, only TO modes are reliably calculated here in the absence of the LO-TO term. The wurtzite GaN belongs to the C46v space group with four atoms per primitive cell. Therefore, its phonon calculations in wurtzite structure are more complicated than the cubic GaN due to the increased unit cell volume and uniaxial symmetry. The hexagonal GaN has 12 spaghetti-like phonon branches, out of which six are Raman active, and four are infrared. The calculated phonon frequencies are summarized in Table 5 along with experiments54,4,56,57 and other ab initio results.22,36,5154 Our computed value of E12 is in excellent conformity with the experiments4,5557 and theoretical results.22,36,5154 B11 is although not available experimentally, but its value agrees excellently with other theoretical results22,36,5153 except that of Shimada et al.,54 which apparently seems to be more overestimated. Our computed values for A1(TO), E1(TO), and Eh2 are underrated 44.3%, 55.5, and 55.4, respectively, from the experimental measurements4,5557 and show a good harmony with its other LDA results, whereas Bh1 differs from 11.4 to +8.4% than that from its theoretical counterparts.22,36,5154 This difference from DFPT based LDA findings might arise from the use of other potential (GGA) that overestimated the lattice constants and hence resulted in minor underestimation from its LDA counterpart. The angular dependence of TO modes is described by Bungaro et al.53 as [ωTO(E1) ωTO(A1)]/ωTO(E1). This angular dependence is a measure of anisotropy of wurtzite GaN, whose computed value is 0.0377, closely related with 0.04 of that proposed by Saib.51 The pressure dependence of phonon frequencies against hydrostatic pressure (GPa) is plotted in Figure 4. Here, all the phonon frequencies except E12 increase with the increase of pressure. Moreover the rate of change of ωTO(E1) is greater as compared to the ωTO(A1), consistent with Saib.51 This leads to the fact that, when TO modes are increased with pressure, the crystal anisotropy must increase. The volume dependence of phonon’s frequencies can be evaluated through a parameter called Gruneisen parameter (γj), which is defined as follows γj ¼
d ln ωj d ln ωj ¼ B0 d ln V dp
where ωj represents the phonon frequency relative to a specific mode j and B0 is the bulk modulus at zero pressure. 14506
dx.doi.org/10.1021/jp207141k |J. Phys. Chem. A 2011, 115, 14502–14509
The Journal of Physical Chemistry A
ARTICLE
Table 5. Comparison of Zone-Centered Phonon Frequencies calculations
experiment
present study
other51
ref 22
ref 36
ref 52
ref 53
ref 54
ref 55
ref 4
ref 56
ref 57
E12 B11
143
136
146
143
150
138
185
144
144
145
143
327
332
335
337
330
334
526
A1(TO)
510
534
534
541
537
550
544
533
531
533
533
E1(TO)
530
556
556
568
555
572
566
558
560
561
559
Eh2
539
565
560
579
558
574
557
568
568
570
568
Bh1
659
689
697
720
677
690
584
mode
Table 7. Comparison of First and Second-Order Pressure Coefficient, α (cm1 GPa1) and β (cm1 GPa2) for Optical Zone Centered Phonons experiment present study mode
Figure 4. Pressure dependence of phonon frequencies.
Table 6. Gruneisen Parameter (γj) for Wurtzite GaN experiment mode
theory
present study
ref 58
ref 59
ref 51
ref 26
ref 52
E12
0.48
0.4
0.43
0.4
B11
0.81
A1(TO)
1.42
1.51
E1(TO) Eh2
1.403 1.514
1.41 1.50
Bh1
1.35
0.35
0.20
0.94
0.99
1.04
1.18
1.51
1.21
1.52
1.61 1.80
1.47 1.56
1.19 1.28
1.48 1.60
1.30
1.29
Using above equation, we have approximated the mode Gruneisen parameter, and the results are listed in Table 6. To analyze our computed results, we have also included the existing experimental and theoretical results. This table finds that our calculated values of the Gruneisen parameter for the frequencies E12, E1(TO), and Eh2 are in excellent agreement with experiment.58 Because no experimental data for Gruneisen parameters B11 and Bh1 are available, we compare them with existing theoretical results and find that our results are consistent with that of other calculated results.51 Our calculated Gruneisen parameters for E1(TO) and Eh2 are underestimated whereas that of A1(TO) is overestimated with respective experimental Gruneisen parameters measured by Perlin et al.59 The Gruneisen parameters of this mode are overestimated with that of Wagner et al.26 and in reasonable agreement with ref 52. So, as a whole, our study based on GGA in conjunction with ultrasoft pseudopotential can be considered a good addition in phonon calculation too. Table 7 summarizes our calculated first- (α) and secondorder (β) derivatives of relative phonon frequencies. Some previous experimental and DFT based results are presented for
α
β
theory
ref 4 α
E12
0.398 0.0002
B11
1.45 0.005
A1(TO)
4.16 0.023 3.99
E1(TO)
4.2
Eh2
5.00 0.03
Bh1
5.21 0.0415
β
ref 22 α 0.015
0.035
0.02 4.63 0.01
β
ref 51 α
β
0.006 0.28 0.00387
1.72 0.007
1.60 0.01032
4.08
0.024
4.07 0.01848
4.10
0.013
4.13 0.01904
4.46
0.018
4.51 0.02004
4.36
0.033
4.54 0.01933
comparison too. These force constants are calculated by fitting the zone centered phonon frequencies with respect to pressure on the following relation: ωðpÞ ¼ ωð0Þ þ αp þ βp2 No experimental values of these force constants are available for the E12 mode; however, for this frequency mode, our results are slightly overestimated as compared to other first principles results.22,51 The phonon frequency for the E12 mode is negative, which represents a decrease of its value with respect to pressure, where the sum α + β is greater than that in ref 51, showing a relatively greater change. However, all other first derivative α = ∂ω/∂p are positive, and all are in good agreement with other results, except for Eh2 and Bh1, which are slightly overestimated as compared to its theoretical counterparts.
IV. CONCLUSION In the present study, PerdewWang27 proposed generalized gradient approximation is implemented along with ultrasoft pseudopotential to calculate the structural, elastic constant, and phonon properties of wurtzite GaN. This study reveals that our calculated lattice constants are in good agreement with others calculated through GGA potential but are 1.3% overestimated as compared to the experiment and LDA based calculations. The axial ratio is underestimated by 0.6% compared with the ideal value. The internal parameter (u) is estimated to be in excellent agreement with ideal value of 0.375 and with other first principles studies. The bulk modulus calculated through the Murnaghan equation of state is 173 GPa, which is very close to the experimental value of Xia et al.34 (188 GPa) but is underestimated significantly from those suggested 14507
dx.doi.org/10.1021/jp207141k |J. Phys. Chem. A 2011, 115, 14502–14509
The Journal of Physical Chemistry A by Ueno et al.35 (245 GPa), which is also inconsistent with its experimental counterpart. The calculated bulk modulus also agrees well with other theoretically calculated values of the bulk modulus using LDA and GGA methods. On the other hand, the pressure derivative is almost in excellent agreement with experimental and theoretical results. The lattice constant (a, c) exhibited a monotonic decrease with pressure. The axial ratio (c/a) and the internal parameter (u) remained constant. This behavior of lattice constants, axial ratio, and internal parameter (u) is in good agreement with other experimental and ab initio studies. The relative volume displayed a monotonic decrease with hydrostatic pressure. The relative volume change at this pressure was measured to be 17.8, consistent with the experimental value of 18. The equilibrium bond length was 1.973 Å near the experimental value of 1.93 Å . The bond lengths parallel and perpendicular to the c-axis decreased with the increase of pressure with a sudden change at 52 GPA. Because wurtzite GaN is a very hard material, the bond angles remained approximately constant with the increase of pressure in accordance with that of Wagner's theoretical studies. The elastic constants calculated through elastic constant theory and Martin’s transformation show a reasonable consistency with other experimental and theoretical results except for a smaller value of c13. The bulk modulus calculated using these two methods is better than that approximated through the Murnaghan31 equation of state. Using the concept of elastic constants, we have calculated shear wave speeds in [100] and [001], where we have found that the values of the wave speeds VL, VT,[001], and Vl calculated through Martin’s49 transformed elastic constants are credible with that calculated by Truel et al.50 than VT,[010] and Vt. The calculated strain relative to the fractional volumetric change εxx/(Δv/v0) and εzz/(Δv/v0) for hexagonal GaN under hydrostatic pressure are found to be in good agreement with experiment and other calculations. The phonon calculations are made through the “Supercell method”, where the phonon frequencies of various modes showed good agreement with experiments and LDA based calculations. These results showed a few percent variations from the other results, which makes our calculations reliable within DFT limits. This difference of phonon frequencies between our results and others may arise due to the fact that we have calculated phonon frequencies at geometrically optimized lattice constants rather than experimental ones. Because the lattice constants in our case are overestimated, the phonon frequencies are slightly underestimated. To determine the crystal anisotropy, the angular dependence is measured and found to be in accordance with that of Saib.51 The pressure dependence of phonon frequencies depicted an increase in all phonon frequencies with pressure modes except E12. The volume dependence of phonon frequencies is evaluated by finding the Gruneisen parameter (γj), which is harmonious with existing data. The first- and second-order derivatives of phonon frequencies are calculated and found to be comparable with experiment and theory.
’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected] (Chuanbao Cao) Tel: +86 10 6891 3792, Fax: +86 10 6891 2001.
’ ACKNOWLEDGMENT This work was supported by National Natural Science Foundation of China (20471007, 50972017) and the Research Fund
ARTICLE
for the Doctoral Program of Higher Education of China (20101101110026).
’ REFERENCES (1) Akasaki, I.; Amano, H. Jpn. J. Appl. Phys. 1989, 28, L 2112–L2114. (2) Eun-Ae, Choi; Chang, K. J. Physica B 2007, 401, 319–322. (3) Perlin, P.; Jauberthie Carillon, C.; Lti, J. P.; Miguel, A. S.; Gzegory, I. High Pressure Res. 1991, 7, 96–98. (4) Perlin, P.; Jauberthie Carillon, C.; Lti, J. P.; Miguel, A. S.; Gzegory, I. Phys. Rev. B 1992, 45, 83–89. (5) Vurgaftman, I.; Meyer, J. R. J. Appl. Phys. 2003, 94, 3675–3697. (6) Nakamura, S.; Fasol, G. The Blue Laser Diode—GaN Based Light Emitters and Lasers; Springer: Berlin, Germany, 1997. (7) Brantley, W. A. J. Appl. Phys. 1973, 44, 534–536. (8) Davydov, V. Yu.; Kitaev, Yu E.; Goncharuk, I. N.; Smirnov, A. N.; Graul, J.; Semchinova, O.; Uffmann, D.; Smirnov, M. B.; Mirgorodsky, A. P.; Evarestov, R. A. Phys. Rev. B 1998, 58, 12899–12907. (9) Deguchi, T.; Ichiryu, D.; Toshikawa, K.; Sekiguchi, K.; Sota, T.; Matsuo, R.; Azuhata, T.; Yamaguchi, M.; Yagi, T.; Chichibu, S.; Nakamura, S. J. Appl. Phys. 1999, 86, 1860–1866. (10) Savastenko, V. A.; Sheleg, A. U. Phys. Status Solidi A 1978, 48, K135–K139. (11) Wagner, J. M.; Bechstedt, F. Phys. Rev B 2002, 66, 115202– 115222. (12) Polian, A.; Grimsditch, M.; Grzegory, I. J. Appl. Phys. 1996, 79, 3343–3345. (13) Pickett, W. E. Comput. Phys. Rep. 1989, 9, 115–198. (14) Baroni, S.; Giannozzi, P.; Testa, A. Phys. Rev. Lett. 1987, 58, 1861–1864. (15) King Smith, D.; Needds, R. J. J. Phys.: Condens. Matter 1990, 2, 3431–3444. (16) Manchon, D.; Barker, A. S.; Dean, P. J.; Zetterstrom, R. B. Solid State Commun. 1970, 8, 1227–1231. (17) Giehler, M.; Ramsteiner, M.; Brandt, O.; Yang, H.; Ploog, K. H. Appl. Phys. Lett. 1995, 67, 733–736. (18) Barker, A. S.; Ilegems, M. Phys. Rev. B 1973, 7, 743–750. (19) Demangeot, F.; Frandon, J.; Renucci, M. A.; Briot, O.; Gil, B.; Aulombard, R. L. MRS Internet J. Nitride Semicond. Res. 1996, 1, U183–U189. (20) Wieser, N.; Klose, M.; Dassow, R.; Rohr, G. C.; Scholz, F. J. Mater. Sci. Eng., B 1997, 50, 88–92. (21) Klose, M.; Wieser, N.; Rohr, G. C.; Dassow, R.; Scholz, F. J. Cryst. Growth 1998, 189190, 634. (22) Miwa, K.; Fukumoto, A. Phys. Rev. B 1993, 48, 7897–7902. (23) Kim, K.; Lambrechet, W. R.; Segall, B. Phys. Rev. B 1994, 50, 1502–1505. (24) Kim, K.; Lambrechet, W. R.; Segall, B. Phys. Rev. B 1996, 53, 13310–13313. (25) Perlin, P.; Suski, T.; Ager, J. W.; Polian, G. A.; Christensen, N. E.; Gorczyca, I.; Grzegory, I.; Weber, E. R.; Haller, E. Phys. Rev. B 1999, 60, 1480–1483. (26) Wagner, J. M.; Bechstedt, F. Phys. Rev. B 2000, 62, 4526–4534. (27) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. Rev. B 1992, 46, 6671–6687. (28) Vanderbilt, D. Phys. Rev. B 1990, 41, 7892–7895. (29) Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1976, 13, 5188–5192. (30) Monkhorst, H. J.; Pack, J. D. Phys. Rev. B 1997, 16, 1748–1749. (31) Murnaghan, F. D. Proc. Natl. Acad. Sci. U. S. A. 1944, 30, 244. (32) Montanari, B.; Harrison, N. M. Chem. Phys. Lett. 2002, 364, 528–534. (33) Maruska, H. P.; Tietjen, J. J. Appl. Phys. Lett. 1969, 15, 327–330. (34) Xia, H.; Xia, Q.; Ruoff, A. L. Phys. Rev. B 1993, 47, 12925– 12928. (35) Ueno, M.; Yoshida, M.; Onodera, A.; Shimomura, O.; Takemura, K. Phys. Rev. B 1994, 49, 14–21. (36) Karch, K.; Wagner, J. M.; Bechstedt, F. Phys. Rev. B 1998, 57, 7043–7049. 14508
dx.doi.org/10.1021/jp207141k |J. Phys. Chem. A 2011, 115, 14502–14509
The Journal of Physical Chemistry A
ARTICLE
(37) Kim, K.; Lambrechet, W. R.; Segall, B. Phys. Rev. B 1996, 53 (16), 310–16326. (38) Van Camp, P. E.; Van Doren, V. E.; Devreese, J. T. Solid State Commun. 1992, 81, 23–26. (39) Arbouche, O.; Belgoumene, B.; Soudini, B.; Driz, M. Comput. Mater. Sci. 2009, 47, 432–438. (40) Stamp, C.; Van De Walls, C. G. Phys. Rev. B 1999, 59, 5521– 5535. (41) Wagner, J. M.; Bechstedt, F. Phys. Rev. B 2000-I, 62, 8003– 8011. (42) Palummo, M.; Bertoni Carlo, M.; Reining, L.; Finocchi, F. Physica B 1993, 185, 404–409. (43) Reeber, R.; Wang, K. MRS Internet J. Nitride Semicond. Res. 2001, 6, 1–5. (44) Takagi, Y.; Ahart, M.; Azuhata, T.; Sota, T.; Suzuki, K.; Nakamura, S. Physica B 1996, 219 & 220, 547–549. (45) Sheleg, A. U.; Savastenko, V. A. Izv. Akad. Nauk SSSR, Neorg. Mater. 1979, 15, 1598. (46) Wagner, J. M.; Bechstedt, F. Phys. Rev B 2002, 66, 115202– 115222. (47) Right, A. F. J. Appl. Phys. 1997, 82, 2833–2840. (48) Zahid, Usman; Chuabao, Cao; Ghulam, Nabi; Dou, Yan Kun; Waheed, S. Khan; Tariq, Mehmood; Sajad, Hussain. J. Phys. Chem. A 2011, 24, 6622–6628. (49) Martin, R. M. Phys. Rev. B 1972, 6, 4546–4553. (50) Truell, R.; Elbaum, C.; Chick, B. Ultrasonic Methods in Solid State Physics; Academic Press: New York, 1969. (51) Saib, S.; Bouarissa, N.; Rodriguez-Hernandez, P.; Munoz, A. Semicond. Sci. Technol. 2009, 24, 025007–025013. (52) Gorczyca, I.; Christensen, N. E.; Pelzer y Blanca, E. L.; Rodriguez, C. Phys. Rev. B 1995, 51, 11936–11939. (53) Bungaro, C.; Rapcewicz, K.; Bernholc, J. Phys. Rev. B 2000, 61, 6720–6725. (54) Shimada, K.; Sota, T.; Suzuki, K. J. Appl. Phys. 1998, 84, 4951– 4959. (55) Nakahara, J.; Kuroda, T.; Amano, H.; Akasaki, I.; Minomura, S.; Grzegory, I. Ninth Symposium Record of Alloy Semiconductor Physics and Electronics, Izunagaoka, Japan; 1990; p 391. (56) Filippidis, L.; Siegle, H.; Hoffmann, A.; Thomsen, C.; Karch, K.; Bechstedt, F. Phys. Status Solidi B 1996, 198, 621–627. (57) Cingolani, A.; Ferrara, M.; Lugara, M.; Scamarcio, G. Solid State Commun. 1986, 58, 823–824. (58) Siegle, H.; Go~ni, A. R.; Thomson, C.; Ulrich, C.; Syassen, K.; Sch€ottker, B.; As, D. J.; Schikora, D. In Gallium Nitride and Related Materials II; Abernathy, C. R., Amano, H., Zolper, J. C., Eds.; Materials Research Society: Pittsburgh, PA, 1997; pp 225. (59) Perlin, P.; Polian, A.; Itie, J. P.; Grzegory, I.; Litwin-Staszewska, E.; Suski, T. Physica B 1993, 185, 426–427.
14509
dx.doi.org/10.1021/jp207141k |J. Phys. Chem. A 2011, 115, 14502–14509
