以第一原理計算分析砷化鎵與磷化銦紅外線雙聲子特徵譜線之研究

(1)

國立交通大學

電子工程學系電子研究所

碩士論文

以第一原理計算分析砷化鎵與磷化銦

紅外線雙聲子特徵譜線之研究

Analysis of two-phonon infrared spectral features

of gallium arsenide and indium phosphide by ab

initio calculation

研究生：林煌翔

Student: Huang-Hsiang Lin

指導教授：顏順通博士

Advisor: Dr. Shun-Tung Yen

指導教授：

Advisor:

(2)

以第一原理計算分析砷化鎵與磷化銦

紅外線雙聲子特徵譜線之研究

Analysis of two-phonon infrared spectral features of

gallium arsenide and indium phosphide by ab initio

calculation

研究生：林煌翔

Student: Huang-Hsiang Lin

指導教授：顏順通博士

Advisor: Dr. Shun-Tung Yen

指導教授：

國立交通大學

電子工程學系電子研究所

碩士論文

A Thesis

Submitted to Department of Electronics Engineering and

Institute of Electronic

College of Electrical and Computer Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of

Master of Science

in

Electronic Engineering

September 2014

Hsinchu, Taiwan, Republic of China

(3)

以第一原理計算分析砷化鎵與磷化銦紅外線特徵

譜線之研究

學生：林煌翔

指導教授：顏順通博士

國立交通大學

電子工程學系電子研究所碩士班

摘要

我們使用密度泛函微擾理論 (Density functional perturbation theory) 作自洽計算來分析雙聲子過程 (two-phonon process) 對於砷化鎵以及磷化銦在紅外線特頻譜的貢獻，這些特徵譜線可以藉由連結實驗特徵峰值與雙聲子聯合能態密度 (two-phonon joint density of states) 中的范霍夫奇異點 (van Hove singularity) 來完成標記，而我們也分析了在第一布里淵區 (first Brillouin zone) 中的相關臨界點 (critical points)。我們發現砷化鎵以及磷化 銦所有的臨界點都只會分布在對稱點 Γ, X, L, W ; 對稱軸 Σ, Λ, ∆, Q, S; 或是對稱面

ΓKL, ΓXWK, ΓXUL。

(4)

Analysis of two-phonon infrared spectral features of gallium

arsenide and indium phosphide by ab initio calculations

Student: Huang-Hsiang, Lin

Advisor: Shun-Tung, Yen

Department of Electronics Engineering

Institute of Electronics

National Chiao Tung University

ABSTRACT

We perform a self-consistent calculation based on the density functional perturbation theory to analyze the infrared spectral features of GaAs and InP contributed by two-phonon processes. Assignment of the features is made by connecting the experimental spectral peaks to the cal-culated Van Hove singularities of two-phonon joint density of states. The associated critical points in the first Brillouin zone are analyzed. We find that most critical points of GaAs and InP are located at symmetry points of Γ, X, L and W, symmetry lines of Σ, Λ, ∆, Q, and S, or symmetry faces of ΓKL, ΓXWK and ΓXUL.

(5)

Acknowledgments

I would like to thank for my advisor Dr.Shun-Tung Yen for his enthusiastic guidance, im-mense knowledge and inspiring comments. I learned a lot from his rigorous attitude on scientific research. Also, I would like to thank my friends in the Quantum and Optical Electronics Lab for their help. Finally, I would like show my sincere thankfulness to my parent for their kindness and supporting throughout my life.

(6)

List of Figures

1.1 The spectrum of electromagnetic wave. Terahertz radiation is ranged from 0.1-10 THz which is between the microwave and infrared. [1] . . . 2 1.2 A demonstration of terahertz image through sole and heel.[1] . . . 3 1.3 3D imaging of teeth reveals the cavity.[1] . . . 3 2.1 Two phonon processes. The sum absorption process involves absorbing a

pho-ton, virtually exciting an optically active TO phonon, and eventually decay-ing into two phonons, while in the difference process the virtually excited TO phonon merges with an existing phonon into a new one. . . 11 2.2 Illustration of types of critical points . . . 13 3.1 Single-phonon dispersions of (a) GaAs and (b) InP with the density-of-states

spectra displayed in the right panels. The curves are obtained by the first-principles calculation. The dots are for neutron scattering data[2, 3]. . . 20 3.2 The curve shows an individual two-phonon density of states of GaAs with the

branch combination ν1+ν1. . . 21 3.3 The first Brillouin zone of the face-centered cubic lattice with an irreducible

wedge highlighted. . . 22 3.4 Type P1 critical point behaviors with different scale cut-off parameter qc . . . . 23

(9)

3.5 The curves demonstrate how a DOS with only type P0 critical point is influ-enced by the background density of states with type P1 critical point. The blue/dashed curve is the DOS Dc(ω)/Db(ω) with a type P0/P1 critical point and the red curve represents the superposition of the both DOS. An inflection point has been generated during the superposition. . . 24 3.6 The curves demonstrate how a DOS with only type P1 critical point is

influ-enced by the background density of states with type P2 critical point. The blue/dashed curve is the DOS Dc(ω)/Db(ω) with a type P1/P2 critical point and the red curve represents the superposition of the both DOS. . . 25 3.7 The calculated total two-phonon density-of-states spectra of phonon sum

pro-cesses (in red) and phonon difference propro-cesses (in green) for GaAs. Some of the features are indicated with their spectral positions (in cm−1), corresponding phonon branch combinations and critical points (in 2π divided by the lattice constant), symmetry points (or lines, or surfaces) where the critical points are located, and the types of features. The straight red lines in the left panel indicate the spectral positions (in cm−1) of absorption features observed by Kotele et al.[4] 27 3.8 The calculated total two-phonon density-of-states spectra of phonon sum

pro-cesses (in red) and phonon difference propro-cesses (in green) for InP. Some of the features are indicated with their spectral positions (in cm−1), corresponding phonon branch combinations and critical points (in 2π divided by the lattice constant), symmetry points (or lines, or surfaces) where the critical points are located, and the types of features. The straight red lines in the left panel indicate the spectral positions (in cm−1) of Raman features observed by Bedel et al. [5] 29 3.9 One dimensional phonon dispersion relation of a linear chain in which alternate

ions have mass M1 and M2, and only nearest neighbors interact. Two cases of relative order between M1and M2have been shown respectively. . . 30 3.10 location of critical points projected onto the surface of the irreducible wedge . . 32

(10)

A.1 Individual two phonon density of states for GaAs with branches combinations 2ν1, ν1+ ν2, ν1+ ν3, ν1+ ν4, ν1+ ν5, ν1 + ν6, 2ν2, ν2+ ν3, ν2+ ν4. . . 36 A.2 Individual two phonon density of states for GaAs with branched combinations

ν2+ ν5, ν2+ ν6, 2ν3, ν3+ ν4, ν3+ ν5, ν3+ ν6, 2ν4, ν4+ ν5, ν4+ ν6, 2ν5, ν5+ ν6, 2ν6. . . 37 A.3 Individual two phonon density of states for GaAs with branched combinations

ν2− ν1, ν3− ν1, ν4− ν1, ν5− ν1, ν6− ν1, ν3− ν2, ν4− ν2, ν5− ν2, ν6− ν2. . 38 A.4 Individual two phonon density of states for GaAs with branched combinations

ν4− ν3, ν5− ν3, ν6− ν3, ν5− ν4, ν6− ν4, ν6− ν5. . . 39 A.5 Individual two phonon density of states for InP with branches combinations

2ν1, ν1+ ν2, ν1+ ν3, ν1+ ν4, ν1+ ν5, ν1 + ν6, 2ν2, ν2+ ν3, ν2+ ν4. . . 40 A.6 Individual two phonon density of states for InP with branched combinations

ν2+ ν5, ν2+ ν6, 2ν3, ν3+ ν4, ν3+ ν5, ν3+ ν6, 2ν4, ν4+ ν5, ν4+ ν6, 2ν5, ν5+ ν6, 2ν6. . . 41 A.7 Individual two phonon density of states for InP with branched combinations

ν2− ν1, ν3− ν1, ν4− ν1, ν5− ν1, ν6− ν1, ν3− ν2, ν4− ν2, ν5− ν2, ν6− ν2. . 42 A.8 Individual two phonon density of states for InP with branched combinations

ν4− ν3, ν5− ν3, ν6− ν3, ν5− ν4, ν6− ν4, ν6− ν5. . . 43

(11)

List of Tables

2.1 Pseudopotential files of the elements Ga, As, In and P. . . 17 3.1 The summarized characteristics for the DOS feature type. . . 26 3.2 The associated coordinates of symmetry points and symmetry axes. Note that

qx, qy, qz are reduced wavevector with unit 2π/a. . . 33 A.1 Two phonon spectral feature assignment for GaAs . . . 49 A.2 Two phonon spectral feature assignment for InP . . . 55

(12)

C H A P T E R

1

Introduction

(13)

2 Recently, terahertz/far-infrared radiation has growing applications in many areas such as medical imaging, security, scientific use and communication due to its excellent transmission properties. By definition, it refers to the electromagnetic wave with the frequency range from 0.1-10 THz (see Fig. 1.1) or wavelength from 30 µm-3mm.

Figure 1.1: The spectrum of electromagnetic wave. Terahertz radiation is ranged from 0.1-10 THz which is between the microwave and infrared. [1]

The electromagnetic wave in terahertz regime can penetrate a wide variety of non-metallic materials like clothing, paper, wood, masonry, plastic, ceramic except for water [6]. Therefore, a lot of applications are visible. For example, it can be used in surveillance for uncovering the concealed weapons (see Fig. 1.2). On the other hand, it is a non-ionizing radiation which causes low-damage to human bodies. Hence, it can be applied in medical imaging such as epithelial cancer detection and 3D imaging of teeth (see Fig. 1.3). Moreover, it can be used in high-speed Wi-Fi (40-1200 times faster)[7] due to its large bandwidth. Miscellaneous applications in chemistry and astronomy are also emerging[6].

The optical properties, such as photon absorption/emission, of dielectric crystals in this regime are usually underlain by lattice vibration or, in quantum-mechanical terms, the single/multi-phonon processes.[8] To analyze the problems on single/multi-phonon, we use the famous ab initio methods based on the combination of density functional theory (DFT) and density functional perturba-tion theory(DFPT). Ab initio methods based on these methods have been common and well

(14)

3

Figure 1.2: A demonstration of terahertz image through sole and heel.[1]

Figure 1.3: 3D imaging of teeth reveals the cavity.[1]

established tools for studying vibrational properties of materials. In contrast, early calculations based on empirical modeling approaches, such as the shell model,[9] the bond charge-charge model, [10] and the rigid-ion model [11] have a serious limitation. For instance, the matrix elements in the associated interatomic force constant (electronic part) is usually set as fitting parameters, which are only available with comparing to the experimental dispersion data from neutron scattrering. As a result, they cannot be applied to crystals that are hygroscopic, or diffi-cult to prepare for neutron scattering, such as AlAs. Even for stable materials, exact calculation of phonon dispersion over the whole Brillouin zone is difficult because it requires adjustment of a lot of parameters for fitting to a full-zone set of experimental data of phonon dispersion and usually obtains a much larger deviation in some directions (X-Z-W and W-Q-L). Therefore, the ab initiomethods become much important in predicting the properties of unfamiliar materials. Ab initiomethod like frozen-phonon technique is also able to deal with the problem on phonon. However, it need to construct a supercell commensurate with wavevector q and will cost much

(15)

4 more computation resources. Hence, DFPT is still by far the most efficient way to solve the problem on phonon.

A lot of related works on phonon have already done recently. For example, Giannozzi et al [12] have first utilized the density functional theory (DFT) together with the linear response method to obtain the phonon dispersion for Si, Ge and many III-V semiconductor materials. Then Deinzer and Strauch[13] have further calculated the absorption spectra for Si and Ge by using DFPT. Moreover, Shirley and Lawler[14] compute the absorption spectra for Si and Ge again and futher assign the features by analyzing the associated critical points for the van Hove singularities. However, a first-principles calculation for the III-V semiconductor materials in-frared spectral featrues with critical points analysis is still missing.

In this study, we analyze two-phonon infrared spectral features of III-V semiconductors GaAs and InP. There has been much work on the infrared properties of GaAs and InP, [2, 4, 5, 15–18] but there is still lack of a detailed and convincing analysis of their two-phonon spectral features. We first calculate the phonon dispersions of GaAs and InP using DFPT and then the two-phonon densities of states (2PDOS). The 2PDOS spectra exhibit features strongly related to features in the absorption/emission infrared spectra.[13] The influence of the background DOS to the features are also discussed. As a byproduct, the formation of peak features can also be explained by the superposition of two features. We also perform a critical point analysis to further identify the 2PDOS spectral features. We find an interesting distribution of critical points; all points can only lie on symmetry planes of the reciprocal lattice.

We organize the thesis by first explaining the first-principles method and related theory about density of state in chapter 2. In chapter 3, we discuss the result of our calculation and give a detailed analysis and explanation. The miscellaneous materials such as fully assignment results are put in appendix.

(16)

C H A P T E R

2

Theoretical Approach

(17)

2.1 Density funtional theory 6

2.1 Density funtional theory

In many-body condense matter physics, the ground state of the electrons in the crystal lattice should be obtained by solving the coupled Schr¨odinger equation. However, It is technically impossible since the size of the linear system is much too large for the current computation ability. The density functional theory (DFT) provides an alternative approach. Instead of solving the eigenstate directly, it puts the emphasis on the charge density n(r) which still preserves a sufficient information of the system.

2.1.1 Kohn-Sham equations

In the density functional theory (DFT) framework, the ground state density can be obtained by solving the famous Kohn-Sham equations which read

(−¯h 2 2m∇

2 _{+ V}

ext(r) + VH(r) + ˆVxc)ψi(r) = iψi(r), i = 1, 2, ..., N (2.1) where ¯h is reduced Planck constant, m is electron mass, N is the number of Kohn-Sham or-bitals ψi(r), Vext(r) is external potential, VH(r) denotes Hartree potential and ˆVxcis non-local exchange-correlation potential. Vext(r) = − X I ZIe2 | r − RI | (2.2) VH(r) = e2 Z n(r0) | r − r0 _|dr 0 (2.3) ˆ Vxc = δExc[n(r)] δn(r) (2.4)

e is the electron charge and Exc denotes the exchange-correlation energy functional and n(r) represents the non-interating system ground state density of electrons. According to the for-mulation of Kohn-Sham equation, n(r) is identical to the ground state density of interacting many-body system.

(18)

2.2 Density functional perturbation theory 7

2.1.2 Local density approximation

In the local density approximation(LDA), the exchange-correlation energy functional is ap-proximated in the form

E_xcLDA[n(r)] = Z

xc(n(r))n(r)dr (2.5) where xc(n(r)) is the many-body exchange-correlation energy density of the interacting homo-geneous electron gas. We choose the form of xc(n(r)) given by Perdew and Zunger[19] who parametrized the numerical data obtained by Ceperly and Alder [20] with Monte Carlo calcula-tions. And exchange-correlation energy density is written in the sum of exchange contribution xand correlation contribution c. For the unpolarized homogeneous electron gas, we have:

xc(rs) = x(rs) + x(rs) (2.6) x(rs) = − 0.4582 rs (2.7) c(rs) =      −0.1423/(1 + 1.0529√rs+ 0.3334rs) rs ≥ 1 −0.0480 + 0.0311 ln rs− 0.0116rs+ 0.0020rsln rs rs ≤ 1 (2.8) In the above expression, engeries are in Hartrees (1 Hartree = 2 Rydbreg); rsis a dimensionless parameter defined by (4π/3)(rsaB)3 = 1/n where aB denotes the bohr radius. Now we plug Eq.(2.5) into Eq.(2.4), and Eq.(2.4) will become

Vxc(r) = xc(n(r)) + n(r)

dxc(n(r))

dn(r) (2.9)

2.2 Density functional perturbation theory

The charge density response ∆n(r) can be obtained self-consistently as a direct generaliza-tion of the DFT Kohn-Sham equageneraliza-tions. We briefly give the precedure which is usually referred to as the density functional perturbation theory (DFPT)[21].

(19)

2.2 Density functional perturbation theory 8 We begin with an initial charge density response and construct the linear correction of Kohn-Sham potential. Then we can solve the equation of order perturbation and obtain the first-order variation of Kohn-Sham orbitals. Finally, the updated charge density response can be calculated, which leads to the next step. From the first-order perturbation theory, we have:

(HKS − n) |∆ψni = −(∆VKS − ∆n) |ψni (2.10) where HKS = − ¯ h2 2m∇ 2_{+ V} ext(r) + VH(r) + Vxc(r) = − ¯ h2 2m∇ 2_{+ V} KS (2.11) is the unperturbed Kohn-Sham Hamiltonian,

∆VKS(r) = ∆Vext(r) + e2 Z ∆n(r0) | r − r0 _|dr 0 + dVxc(n) dn n=n(r) ∆n(r) (2.12) is the first-order correction of the Kohn-Sham potential, and ∆n = hψn| ∆VKS|ψni is the first-order variation of the Kohn-Sham eigenvalue n. Note that ∆ means the finite-difference operator which is essentially defined by

∆λF = X i

∂F ∂λi

∆λi (2.13)

with the corresponding variable λi.

The charge density can be constructed by the Kohn-Sham orbitals. n(r) =

N X n=1

| ψn(r) |2 (2.14) Direct linearizing of Eq.(2.14) gives

∆n(r) = 2Re( N X n=1

ψ_n∗(r)∆ψn(r)) (2.15) Then equations (2.10),(2.12) and (2.15) form a set of self-consistent equations for the per-turbed system, which is analogous to the unperper-turbed Kohn-Sham equations.

(20)

2.3 Phonons in DFPT 9

2.3 Phonons in DFPT

The most important quantity in the beginning of analyzing the properties of phonon is knowing the interatomic force constant or dynamical matrix. We discuss how to obtain these quantity in the following.

2.3.1 Interatomic force constant

In the harmonic approximation, the interatomic force constant is by definition Ciα,jβ(RI− RJ) ≡

∂2E

∂usα(RI)∂utβ(RJ)

(2.16) where the Greek subscripts indicate Cartesian components, s and t represent for the indices of basis atom, RI and RI mean the coordinates of atomic site I and J, and E is the Born-Oppenheimer energy surface, whose Hamiltonian HBOreads

HBO = He+ Vion−ion

2.3.2 Dynamical matrix

The dynamical matrix is simply given by Dsα,tβ(q) = 1 √ MsMt ˜ Csα,tβ(q) (2.18) where ˜Csα,tβ(q) is the Fourier transform of real space interatomic force constant.

˜

Csα,tβ(q) = X

I

(21)

2.4 Two phonon processes 10 In the monochromatic perturbation RI[us(q)] = Rl+ τs+ us(q)eiq·Rl(l means the l-th lattice site, τsis equilibrium position of the s-th atom in the unit cell), we have

˜ Csα,tβ(q) = 1 Nc ∂2_E ∂u∗ sα(q)∂utβ(q) = ˜C_sα,tβel (q) + ˜C_sα,tβion (q) = 1 Nc [ Z ( ∂n(r) ∂usα(q) )∗∂Vext(r) ∂utβ(q) dr + Z n(r) ∂ 2_V ext ∂u∗ sα(q)∂utβ(q) dr + ∂ 2_E ion−ion ∂u∗ sα(q)∂utβ(q) ] (2.20)

We find that the interatomic force constant in q-space is directly related to the charge density n(r) and charge density response ∂n(r)/∂uα(q), which have already been computed by DFT and DFPT calculation. Then we are able to obtained this quantity.

2.4 Two phonon processes

The two-phonon processes, including the sum and the difference processes(Fig.2.1), have been known as the underlying mechanism for most features in far-infrared absorption/Raman spectra. The sum absorption process involves absorbing a photon, virtually exciting an optically active TO phonon, and eventually decaying into two phonons, while in the difference process the virtually excited TO phonon merges with an existing phonon into a new one. Energy and mo-mentum conservation are preserved during the process. Since the speed of light is much faster than phonon, the photon momentum can be viewed as zero compared to phonon momentum.

2.4.1 Two phonon joint dispersion

The square roots of eigenvalues of the dynamical matrix form the dispersion relation ω(q). The two-phonon joint dispersion is calculated by:

Ω_λ±λ0(q) =| ω_λ(q) ± ω

λ0(q) | (2.21) where the subscript λ is the branch index.

(22)

2.4 Two phonon processes 11

hν=ħω(q

1

)+ ħω(q

2

)

ω

_TO

hν

ω

_(q

₁

₎

ω

_(q

₂

₎

(a) Sum process

hν=ħω(q

2

) - ħω(q

1

)

ω

_TO

hν

ω

_(q

₁

₎

ω

_(q

₂

₎

(b) Difference process

Figure 2.1: Two phonon processes. The sum absorption process involves absorbing a photon, virtually exciting an optically active TO phonon, and eventually decaying into two phonons, while in the difference process the virtually excited TO phonon merges with an existing phonon into a new one.

(23)

2.5 Critical points 12

2.4.2 Two phonon joint density of states

The 2PDOS can be calculated with the following equaiton: D_λ±λ0(ω) = V

(2π)3 Z

BZ

δ(ω − Ω_λ±λ0(q))dq (2.22)

where V is the volume of crystal and the sign + (-) represents the sum (difference) process and the integration runs over the whole first Brillouin zone.

2.5 Critical points

Recall that the formula of density of states can be written as : V (2π)3 Z S(¯hω) dS | ∇¯hω(q) | (2.23) where S is the isoenergy surface with energy ¯hω. A critical point of a two-phonon joint dis-persion Ω_λ±λ0(q) is the one in the reciprocal space at which the gradient of Ω

λ±λ0(q) vanishes. The critical points correspond to van Hove singularities of the 2PDOS spectrum D_λ±λ0(ω).

The dispersion ω(k) can be expanded about the critical point as ω(k) = ω(k0) ± 1 2cxk 2 x± 1 2cyk 2 y ± 1 2czk 2 z (2.24)

cx, cy, czare corresponding to the curvature of the phonon dispersion in the x,y and z directions. The single-phonon DOS spectra have five types of singular features denoted by P , P0, P1, P2, P3, respectively. The shapes of all types of critical points are shown in Fig.2.2

2.5.1 Type P

0

critical points

- P0 denotes a critical points with zero negative sign in the Eq.(2.24), the critical point is a minimum. For type P0 critical points about the frequency ω0

(24)

2.5 Critical points 13

(a)Type P0 (b)Type P1

(c)Type P2 (d)Type P3

(e)Type P (Peak)

(25)

2.5 Critical points 14 D(ω) = V (2π)3 Z δ(ω0+ 1 2cxk 2 x+ 1 2cyk 2 y + 1 2czk 2 z − ω)dk = V (2π)3 s 8 cxcycz Z δ(ω0+ q2x+ q 2 y + q 2 z − ω)dqxdqydqz = V π2 s 2 cxcycz Z δ(ω0+ q2− ω)q2dq = V π2_p2c xcycz Z δ(ω0+ u − ω) √ udu =      V π2√_2c xcycz √ ω − ω0 ω > ω0 0 ω < ω0 (2.25)

2.5.2 Type P

1

critical points

- P1 denotes a critical point with a negative sign in the Eq.(2.24), the critical point is a saddle point. For type P1 critical points about the frequency ω1 :

D(ω) = V (2π)3 Z δ(ω1+ 1 2cxk 2 x+ 1 2cyk 2 y − 1 2czk 2 z− ω)dk = V (2π)3 s 8 cxcycz Z δ(ω1+ qx2+ q 2 y− q 2 z − ω)dqxdqydqz = V π2_p2c xcycz Z δ(ω1+ q2− q2z− ω)qdqdqz = V π2_p2c xcycz Z δ(qz−pω1+ q2 − ω) + δ(qz+pω1+ q2− ω) 2pω1 + q2− ω qdqdqz , q2 > ω − ω1 = V π2_p2c xcycz Z qdq pω1+ q2− ω , q2 > ω − ω1 = V π2_p2c xcycz        R_√qc ω−ω1 qdq √ ω1+q2−ω ω > ω1 Rqc 0 qdq √ ω1+q2−ω ω < ω1

(26)

2.5 Critical points 15 = V π2_p2c xcycz      pω1+ qc2− ω ω > ω1 pω1+ qc2− ω − √ ω1− ω ω < ω1 =      C1 ω > ω1 C1− _π2√_2cV xcycz √ ω1− ω ω < ω1 (2.26) where qcis a cut-off parameter, and C1 is a smoothly frequency-dependent quantity.

2.5.3 Type P

2

critical points

- P2denotes a critical point with two negative signs in the Eq.(2.24), the critical point is a saddle point.

For type P2critical points about the frequency ω2

D(ω) = V (2π)3 Z δ(ω2− 1 2cxk 2 x− 1 2cyk 2 y+ 1 2czk 2 z − ω)dk = V (2π)3 s 8 cxcycz Z δ(ω2− q2x− q 2 y + q 2 z − ω)dqxdqydqz = V π2_p2c xcycz Z δ(ω2− q2+ q2z− ω)qdqdqz = V π2_p2c xcycz Z δ(qz−pω + q2− ω2) + δ(qz+pω + q2− ω2) 2pω + q2_{− ω} 2 qdqdqz , q2 > ω2− ω = V π2_p2c xcycz Z qdq pω + q2_{− ω} 2 , q2 > ω2− ω = V π2_p2c xcycz        Rqc 0 qdq √ ω+q2_−ω 2 ω > ω2 R_√qc ω2−ω qdq √ ω+q2_−ω 2 ω < ω2 = V π2_p2c xcycz      pω + q2 c − ω2− √ ω − ω2 ω > ω2 pω + q2 c − ω2 ω < ω2

(27)

2.5 Critical points 16 =      C2− _π2√_2cV xcycz √ ω − ω2 ω > ω2 C2 ω < ω2 (2.27)

2.5.4 Type P

3

critical points

- P3 denotes a critical point with three negative signs in the Eq.(2.24), the critical point is a maximum .

For type P3critical points about the frequency ω3

D(ω) = V (2π)3 Z δ(ω3− 1 2cxk 2 x− 1 2cyk 2 y − 1 2czk 2 z − ω)dk = V (2π)3 s 8 cxcycz Z δ(ω3− q2x− qy2− q2z− ω)dqxdqydqz = V π2 s 2 cxcycz Z δ(ω3− q2− ω)q2dq = V π2_p2c xcycz Z δ(ω3− u − ω) √ udu =      0 ω > ω3 V π2√_2c xcycz √ ω3− ω ω < ω0 (2.28)

2.5.5 Type P critical points

- P denotes a critical point which can result in ”peak” type of singularity. It can be essentially type P1/P2critical points. This will be discussed in detail in chapter 3.

(28)

2.6 Technical details 17 pseudopotential file Ga Ga.pz-bhs.UPF As As.pz-bhs.UPF In In.pz-hgh.UPF P P.pz-bhs.UPF

Table 2.1: Pseudopotential files of the elements Ga, As, In and P.

2.6 Technical details

Both DFT and DFPT calculation is implemented in the QUANTUM ESPRESSO[22], which is an integrated suite of Open-Source codes for electronic-structure calculations. In the DFT cal-culations (self-consistent charge density and Kohn-Sham orbitals), we sample the Brillouin zone by the 5×5×5 Monkhorst-Pack special k-point method[23] and we use the norm-conserving type pseudopotential with the LDA exchange-correlation funcitonal proposed by von Barth and Car[24] with 25 Ry cut-off energy. The pseudopotential files can be obtained on the website of QUANTUM ESPRESSO. The files used in our calculation are listed in Table. 2.1. The asso-ciated lattice constant we use is 10.49 bohr for GaAs and 11.07 bohr for InP. In the DFPT cal-culation (especially in the charge density response), we first compute the dynamical matrix on 5×5×5 q-point grids. We then perform the Fourier transform from the k-space to the real space and do interpolation. Finally, we compute the dynamical matrix again on the 100×100×100 sampling points. In implementing the integral in Eq.(2.22) to obtain the DOS, the improved tetrahedron method (Bl¨ochl version)[25] is being used. The procedure mentioned above can be achieved by following the guide [26]. However, the two-phonon joint DOS is not implemented. It is mandatory to modify the source code by ourself.

(29)

C H A P T E R

3

Results and Discussion

(30)

3.1 Single Phonon Dispersion Relation and Density of States 19

3.1 Single Phonon Dispersion Relation and Density of States

The single-phonon dispersion relations of GaAs and InP are shown Fig. 3.1. Since the full dis-persion relation is a four dimension manifold, we can only plot the disdis-persion curve in particular lines. In order to compare to the experimental data given by neutron scattering [2, 3], we choose the six direction Γ-∆-X, X-Σ-Γ, Γ-Λ-L, L-X, X-Z-W and W-Q-L which are the same to that in their research. The curves are obtained by the first-principles calculation. The dots are for neutron scattering data. The acoustic phonon branches (ν1,ν2 and ν3) for both materials always give an excellent agreement with the experimental data. For optical phonon branches (ν4,ν5and ν6), calculation results match the experimental data in most directions for GaAs. Even though we use the first-principles calculations, we still have some small discrepancies in X-Z-W and W-Q-L directions. However, it exhibits a great improvement than any other empirical model in this region.[2] For InP, we can only compare the calculation results with experimental data in much fewer directions (Γ-∆-X, X-Σ-Γ and Γ-Λ-L) due to the limitation of number of the neu-tron scattering data. As the same case to the GaAs, the curves match well for acoustic mode. For optical mode, the discrepancies are much obvious than that of GaAs. This problem may be due to the choices of the combinations of the proposed the psuedopotentials and exchange-correlation functionals of indium and phosphorus or essentially owing to the reliability of the such limited experimental data. Of course, the same argument can be applied to GaAs to explain the deviation in X-Z-W and W-Q-L directions. We also compute the associated single phonon density of states (DOS) for both materials. The results are shown in the right-hand side of the dispersion curves. As we can see, there is a number of peak corresponding to the flat or less dispersive region in the dispersion relation. This is an expected consequence by the definition of the DOS. From the Eq.(2.23), the vanishing gradient regions give rise to a large amount of contribution to the DOS. On the other hand, We also check our result by comparing to the oth-ers’ works[21, 27] and we have excellent agreements with them. As a result, we have arrived a steady situation, which is ready for the further analysis of two-phonon joint density of states (2PDOS).

(31)

3.1 Single Phonon Dispersion Relation and Density of States 20

DOS (a. u.) 0 50 100 150 200 250 300 F requency (cm −1 ) Γ ∆ _X Σ Γ Λ _L _{X Z W Q L} ν 1 ν 2 ν 3 ν 4 ν 5 ν 6 (a) GaAs

DOS (a. u.) 0 100 200 300 400 F requency (cm −1 ) Γ ∆ _X Σ Γ Λ _L _{X Z W Q L} ν 1 ν 2 ν 3 ν 4 ν 5 ν 6 (b) InP

Figure 3.1: Single-phonon dispersions of (a) GaAs and (b) InP with the density-of-states spectra dis-played in the right panels. The curves are obtained by the first-principles calculation. The dots are for neutron scattering data[2, 3].

(32)

3.2 Assignment of 2PDOS’s and critical points analysis 21

3.2 Assignment of 2PDOS’s and critical points analysis

The assignment scheme is given by the following:

1. Find the van-Hove singularities in the indiviual 2PDOS spectra.

2. Identify the critical point coordinates in two-phonon dispersion relation.

3. Check the whether the critical point frequencies are identical to the singularity positions. 4. Determine the types of the critical points by comparing the nearby behavior to the stan-dard five types of behaviors in Fig. 2.2.

0 50 100 150 200 250 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035 Frequency (cm−1)

Density of states (a. u.)

128.1 P

1 ,L(0.5,0.5,0.5)

165.4 P ,S(0.14,1,0.14)

207.1 P₃ ,W(0.5,1,0)

Figure 3.2: The curve shows an individual two-phonon density of states of GaAs with the branch com-bination ν1+ν1.

(33)

3.2 Assignment of 2PDOS’s and critical points analysis 22 K W X q_x q_z q_y L U ī ȁ Ȉ ǻ

FIG. 1

Q Z S

Figure 3.3: The first Brillouin zone of the face-centered cubic lattice with an irreducible wedge high-lighted.

The example is demonstrated in Fig. 3.2. First, we find out three van-Hove singularities in the curve. The corresponding position of them are 128.1, 165.4 and 207.1 cm−1. Secondly, we search the vanishing gradient points in the two-phonon dispersion relation with the branch combination ν1+ν1 in the irreducible wedge (the truncated octahedron in Fig. 3.3). We find that the critical points are located at L(0.5,0.5,0.5), S(0.14,1,0.14) and W(0.5,1,0). In the third step, we check the corresponding frequencies of these critical points whether they are exactly or approximately equal to the singularity positions in step 1. If the deviation is acceptable (within 5% tolerance), then we will assign this feature. In this case, the corresponding frequencies at L(0.5,0.5,0.5), S(0.14,1,0.14) and W(0.5,1,0) are 126.78, 163.35 and 208.00 cm−1 respectively. With this observation, we can easily attribute the van-Hove singularities at 128.1, 165.4 and 207.1 cm−1 to the corresponding critical points L(0.5,0.5,0.5), S(0.14,1,0.14) and W(0.5,1,0) respectively. In the final step, we have to determine the types of the critical points by com-paring the nearby behaviors to the standard type P0 (minimum), P1 (saddle point), P2 (saddle point), P3 (maximum) and P (peak) critical points (see Fig. 2.2). For instance, the behavior of critical point at 128.1 cm−1 manifests an similar shape to that in Fig. 2.2(b), which is concave upward on the left-hand side (behavior on the right-hand side is less important and this point

(34)

3.2 Assignment of 2PDOS’s and critical points analysis 23 (a) large qc Type P2 Type P3 (Peak) (b) medium qc Type P2 Type P3 (Peak) (c) small qc

Figure 3.4: Type P1critical point behaviors with different scale cut-off parameter qc

will be discussed later). The same procedure can be applied to the other two van-Hove singu-larities. For critical point at 165.4 cm−1, it shows an prominently sharp behavior similiar to that in Fig. 2.2(e), which is a nondifferentiably local maximum. For critical point at 207.1 cm−1, it exhibits an similar behavior to that in Fig. 2.2(d), which is vanishing on the right-hand side and increasing (concave downward) on the left-hand side.

For a type P1 critical point, the behavior on the right-hand side is less important. This is because that the smoothly frequency-dependent quantity C1 in Eq.(2.26) contains the cut-off parameter qc. C1 behaves like a constant, when qc is large enough which means the nearby dispersion relation: ω(k) = ω(k0) + 1 2cxk 2 x+ 1 2cyk 2 y − 1 2czk 2 z (3.1)

holds the validity in the much wider region. In contrast, C1 becomes vanishing when qcis small which means the dispersion relation in the vicinity holds its validity much narrower domain. Based on this fact, the right-hand side behavior of a type P1 critical point is dependent on the

(35)

3.2 Assignment of 2PDOS’s and critical points analysis 24 region of validityof Eq.(3.1). However, the the left-hand side behavior is always concave up-ward regardless of the region of validity. Therefore, we can conclude that the left-hand side behavior is always concave upward in any situation. This argument can apply to the type P2 critical point as well. We can claim that the behavior on the left-hand side is less important. On the contrary, we face no difficulty upon assigning the type P0, P3 and P critical points since they have well-defined behaviors on the both side.

There is one relevant effect for us to consider which is the influence of the background density of states. The background density of states will smoothly influence the behavior of critical points. The reason for the cause of background density of state can be simply due to the other critical points or essentially by the dispersion behavior beyond the standard quadratic form in Eq.(2.24). The latter one is much difficult to discuss due to the arbitrariness of the function form of dispersion. On the contrary, the former one can be illustrated in the Fig. 3.5.

Ω DHΩL inflection point DcHΩL with type P0 D_bHΩL_{with type P}1 DcHΩL DbHΩL

Figure 3.5: The curves demonstrate how a DOS with only type P0 critical point is influenced by the

background density of states with type P1critical point. The blue/dashed curve is the DOS Dc(ω)/Db(ω)

with a type P0/P1 critical point and the red curve represents the superposition of the both DOS. An

(36)

3.2 Assignment of 2PDOS’s and critical points analysis 25 As we can see, behavior of superposition of the background DOS and the small region DOS with the type P0 critical point has slightly different to the original DOS but it will not affect too much to our assigning criterion. Note that an inflection point may be potentially generated due to the competition of the concavity of the both DOS. If the background concavity can flip the direction of the original concavity, then an inflection point will be generated and will not affect the DOS too much.

Let’s see the next example. A DOS Dc(ω) with type P1 critical point on the left-hand side is superimposed by another DOS Db with type P2 critical point on the right-hand side. The process is demonstrated in the Fig. 3.6. We can find that there are two features on the total DOS Dc(ω) + Db and observe that the peak at the feature 1 is essentially caused by the type P1 critical point. Based on this concept, the peak behavior can be explained. Notice that the total DOS can not be ploted in full region since the individual DOS’s are only simultaneously well-defined in some region.

Ω DHΩL DcHΩLwith type P1 DbHΩLwith type P2 DcHΩL DbHΩL IHDWXUH1SHDN IHDWXUH23

Figure 3.6: The curves demonstrate how a DOS with only type P1 critical point is influenced by the

background density of states with type P2critical point. The blue/dashed curve is the DOS Dc(ω)/Db(ω)

(37)

3.2 Assignment of 2PDOS’s and critical points analysis 26 Dispersion behavior Left-hand side DOS Right-hand side DOS Type P0 minimum vanishing/varying concave downward Type P1 saddle point concave upward varying Type P2 saddle point varying concave upward Type P3 minimum concave downward vanishing/varying Type P (peak) saddle point varying varying

Table 3.1: The summarized characteristics for the DOS feature type.

Now we can fairly say that the feature shape of a critical point would be influenced slightly but still preserves part of the original behaviors. We summarize the criterion of type identifica-tion in the Table. 3.1. Note that even though the feature shape of a peak are varying on the both side, a nondifferentiably local maximum is still a clear characteristic.

Now we are able to specify the features for all of the individual 2PDOS for GaAs and InP. After a very tedious work, we find out 124 features for GaAs and 117 features for InP. The assignment is all listed in the appendix. On the other hand, the relatively prominent features of the total 2PDOS is also identified and shown in Fig. 3.7 and Fig. 3.8. As we can see, there are two types of curves. One (in red) corresponds to the sum process, and the other (in green) to the difference process. There are plenty of experimental peaks in the spectra which are connected to the Van Hove singularities.

3.2.1 Assignment for GaAs

For GaAs, we can see a plenty of solid red straight lines on the left-hand side in Fig. 3.7 which are the experimental feature positions obtaining by Kotele et al [4] via infrared absorption mea-surement. The data are ranging from 108 cm−1to 590 cm−1and there is no observation of peak between 235 cm−1 to 310 cm−1. This is due to that the photons have a strong resonance with

(38)

3.2 Assignment of 2PDOS’s and critical points analysis 27 0 100 200 300 400 500 600 0 0.5 1 1.5 Frequency (cm − 1 )

108 150 161 172.6 195 214 216 230 235 310 323 332 338 353 376 385 410 414 446.5 459.5 475 481 490 498 507 516 531 538 566 590 13.5 ν₆−ν₅ (0.5 0.5 0.5) L P 19.41 ν₂−ν₁ (0.5 0.77 0.23) Q P 69.9 ν₅−ν₃ (0.5 0.76 0.24) Q P 83.6 ν₆−ν₃ (0.5 0.78 0.22) Q P 93.8 ν₃−ν₁ (0.5 1 0) W P 138.2 ν₅−ν₂ (0.5 1 0) W P 151.8 ν₆−ν₂ (0.72 0.72 0) Σ P 165.4 ν₁+ν₁ (0.14 1 0.14) S P 172.6 ν₄−ν₁ (0.5 0.5 0.5) L P 201.1 ν₆−ν₁ (0.5 0.5 0.5) L P 219.5 ν₁+ν₂ (0.5 1 0) W P 3 237.2 ν₂+ν₂ (0.67 0.67 0) Σ P₃ 265.8 ν₁+ν₃ (0.62 0.62 0.25) ΓKL P 313.2 ν₂+ν₃ (0 1 0) X P 319.1 ν₁+ν₄ (0.38 0.38 0) Σ P 330.5 ν₁+ν_5,6 (0.5 0.5 0.5) L P 1 337.1 ν₁+ν₅ (0 1 0) X P 358.8 ν₂+ν₆ (0 0.52 0) ∆ P 372.2 ν₃+ν₃ (0.5 0.78 0.22) Q P 384.0 ν₂+ν₆ (0.5 1 0) W P 406.8 ν₃+ν₄ (0.5 1 0) W P 415.6 ν₃+ν₄ (0.27 0.72 0.27) ΓXUL P₂ 442.5 ν₃+ν₅ (0.5 0.8 0.2) Q P 455.6 ν₃+ν₆ (0.5 0.78 0.22) Q P 496.8 ν₄+ν₆ (0 1 0) X P 1 509.2 ν₄+ν₄ (0 0.73 0) ∆ P 530 ν₅+ν₆ (0.5 0.5 0.5) L P 538.8 ν₆+ν₆ (0.75 0.5 0.25) Q P

Figure 3.7: The calculated total two-phonon density-of-states spectra of phonon sum processes (in red) and phonon difference processes (in green) for GaAs. Some of the features are indicated with their spec-tral positions (in cm−1), corresponding phonon branch combinations and critical points (in 2π divided by the lattice constant), symmetry points (or lines, or surfaces) where the critical points are located, and the types of features. The straight red lines in the left panel indicate the spectral positions (in cm−1) of absorption features observed by Kotele et al.[4]

(39)

3.2 Assignment of 2PDOS’s and critical points analysis 28 the TO phonons in this region which is so-called reststrahlen band. Therefore, the light cannot penetrate through the material and cause the transmittance become approximately zero. On the other hand, the material will also reflect the light strongly and give a roughly unity reflectance. Since the absorption coefficient depends on the transmittance, it will be very sensitive to the background noise and cause a large deviation and unreliable data. That’s the reason why they report empty peak in this region. For the other observed features, we can see that many of the experimental data have been assigned. For example, a feature assigned at 150 cm−1 in the left-hand side of the Fig. 3.7 and the phonon peak position carried out by our calculation is 151.8 cm−1which somewhat indicates a successfulness of our analysis. Besides, the associated branch combination and critical point are ν6−ν2and (0.72,0.72,0) (on the Σ axis). Now, we can clearly explain how the peak here is formed. In the beginning, a incoming photon carried with the frequency (energy) about 150 cm−1 scatters with the GaAs lattice, and transforming into an virtually excited TO phonon. Then the intermediate phonon scatters another existed acous-tic phonon mode ν2 with wavevector or momentum at (0.72,0.72,0) (in 2π/a cm−1) and they combine together become an new optical phonon mode ν2 with the same wavevector or mo-mentum at (0.72,0.72,0) due to the momo-mentum conservation (the photon momo-mentum is roughly zero or negligible compared to that of phonon). There is one point for us to know that the low frequency data below 108 cm−1 are missing due to the experimental limit but we can still give a clear result which is dominated by difference process. For the region inside the reststrahlen band, we have also observed some features. As a result, the features in low frequency region or inside reststrahlen band can be predicted and may be observed by the improvement of the experiment.

3.2.2 Assignment for InP

As for InP, most of the mentioned descriptions can read from the Fig. 3.8 and the procedure is all most the same. As this time, the experimental data is obtained from Bedel et al [5] via

(40)

Ra-3.2 Assignment of 2PDOS’s and critical points analysis 29 0 100 200 300 400 500 600 700 0 0.5 1 1.5 Frequency (cm − 1 )

30 69 154 191 227 253 276 311 385 401 410 441 465 498 624 656 110 195 612 662 55 136 306 348 165 370 525 633 688 697 16.66 ν₂−ν₁ (0.5 1 0) W P 27.01 ν₆−ν₅ (0.33 0.33 0) Σ P 67.41 ν₃−ν₂ (0.27 0.73 0.27) ΓXUL P 79.09 ν₃−ν₁ (0.5 1 0) W P 111.9 ν_1,2+ν_1,2 (0.5 0.5 0.5) L P₁ 142.0 ν₁+ν₁ (0 1 0) X P 155 ν₄−ν₃ (0.5 0.72 0.28) Q P 166.1 ν₅−ν₃ (0.5 1 0) W P 188.6 ν₆−ν₃ (0.5 1 0) W P 198.6 ν₂+ν₂ (0.5 1 0) W P₃ 224.1 ν₁+ν₃ (0.5 0.5 0.5) L P 253.3 ν₅−ν₁ (0 1 0) X P 278.6 ν₆−ν_1,2 (0.5 0.5 0.5) L P 314.3 ν₃+ν₃ (0.5 0.78 0.22) Q P 381.2 ν₁+ν₄ (0.65 0.65 0) Σ P 386.7 ν₁+ν₄ (0.5 1 0) W P 402.9 ν₁+ν₄ (0.5 1 0) W P 410.4 ν₁+ν₆ (0 0.56 0) ∆ P 443.7 ν₂+ν₆ (0.65 0.65 0) Σ P 465.3 ν₃+ν₄ (0.5 1 0) W P 479.0 ν₃+ν₅ (0.5 0.75 0.25) Q P 503.1 ν₃+ν₆ (0.5 0.5 0.5) L P 608.2 ν₄+ν₄ (0.5 1 0) W P 0 631.6 ν₄+ν₅ (0.5 1 0) W P 638.3 ν₄+ν_4,5 (0 0 0) Γ P 658.1 ν₄+ν₆ (0.5 0.75 0.25) Q P 666.5 ν₄+ν₆ (0.35 0.35 0) Σ P 697.9 ν₆+ν₆ (0.5 1 0) W P

Figure 3.8: The calculated total two-phonon density-of-states spectra of phonon sum processes (in red) and phonon difference processes (in green) for InP. Some of the features are indicated with their spectral positions (in cm−1), corresponding phonon branch combinations and critical points (in 2π divided by the lattice constant), symmetry points (or lines, or surfaces) where the critical points are located, and the types of features. The straight red lines in the left panel indicate the spectral positions (in cm−1) of Raman features observed by Bedel et al. [5]

(41)

3.2 Assignment of 2PDOS’s and critical points analysis 30 k Z (a) M1≈ M2 k Z (b) M1 M2

Figure 3.9: One dimensional phonon dispersion relation of a linear chain in which alternate ions have mass M1 and M2, and only nearest neighbors interact. Two cases of relative order between M1and M2

have been shown respectively.

man spectroscopy instead of infrared absorption. Since the wavelength of Raman experiment is around visible light which is completely different order compared to infrared and the measure-ment is based on the scattering of light instead of absorption, it will not face with the problem of reststrahlen effect. The data are ranging from 30 cm−1 to 697 cm−1. It is noteworthy that they still observed a gap between 525 and 612 cm−1and we also obtain a very similiar result. This is a called the fingerprint of semiconductor. The basic reason of this phenomenon is owning to the large mass difference (e.g. MIn = 114.818 amu MP = 30.947 amu). Such semiconductors have much less dispersive optical phonon branches and a wider gap than other semiconductors, such as GaAs (MGa = 69.72 amu ≈ MAs = 74.92 amu). This can easily be explained in one dimensional case. Consider a linear chain in which alternate ions have mass M1 and M2, and only nearest neighbors interact with string constant K. One can use the classical model to show that the dispersion relation for the normal modes is :

ω(k) = r K M1M2 (M1+ M2 ± q M2 1 + M22+ 2M1M2cos ka) (3.2) where a represents the lattice constant. Two extreme cases for relative order between M1 and M2 have been shown in Fig. 3.9. The fact that the gap is much wider than each of the optical phonon branches, as shown in Fig 3.1, causes the fingerprints in the two-phonon DOS spectrum

(42)

3.2 Assignment of 2PDOS’s and critical points analysis 31 of InP. Be more specific, the fingerprint will appear when the condition Ων4+ν4(k) > Ων3+ν6(k)

holds for every k-point in the first Brillouin zone, which means we cannot form a possible linear combination of the branches and give a contribution to the 2PDOS in the fingerprint region. In appendix, this condition is also held.

3.2.3 Critical Point Analysis

The cubic symmetry of the crystals allows us to consider only the critical points in the irredici-ble wedge of the first Brillouin zone as shown in Fig. 3.3. The wedge is one of 48 equivalent parts of the first Brillouin zone and enclosed by five symmetry planes with respective to which reflection leaves the reciprocal lattice invariant.

Specifically, we choose the wedge with the following condition:

qz ≥ 0, qx+ qy + qz ≤ 3/2, qy ≤ 1, qx ≤ qy, qz ≤ qx

The associated coordinates of symmetry points and symmetry axes are listed in Table 3.2. Inter-estingly, a very important phenomenon has been observed: All critical points are only located at symmetry points, on symmetry axes or on symmetry faces. No critical point inside the irre-ducible wedge of the first Brillouin zone(Fig. 3.3) is being found, which means the extrema and saddle points will only appear in these high symmetry domain and this is a reasonable discov-ery. Although the phenomenon may be elegantly explained by the space group, however, it is beyond our discussion.

Since all the critical points are located on the surface of the irreducible wedge of fcc first Brillouin zone, it is straightforward to make a scatter plot of the critical points which is shown in Fig.3.10. As we can see, critical points manifest an aggregation behavior, which means a critical point would be usually accompanied by other critical points in the vicinity. This is due to that the vanishing gradient criterion is replaced by a small tolerance value in practical, which depends on the critical points locations. Some critical points simply need small tolerance and

(43)

3.2 Assignment of 2PDOS’s and critical points analysis 32

Γ

L

U

X

W

K

Q

Λ

Σ

∆

S

Z

(a) GaAs

Γ

L

U

X

W

K

Q

Λ

Σ

∆

S

Z

(b) InP

(44)

3.3 Conclusion 33 symmetry points symmetry axes symmetry faces Γ (0,0,0) Z (X-W) (ξ,1,0) ΓXUL qx = qz X (0,1,0) ∆ (Γ-X) (0,ξ,0) ΓXWK qz = 0 K (0.75,0.75,0) Σ (Γ-K) (ξ,ξ,0) ΓKL qx = qy L (0.5,0.5,0.5) Λ (Γ-L) (ξ,ξ,ξ) LUWK qx+ qy + qz = 3/2 W (0.5,1,0) Q (L-W) (0.5,1-ξ,ξ) XWU qy = 1 U (0.25,1,0.25) S (U-X) (ξ,1,ξ)

Table 3.2: The associated coordinates of symmetry points and symmetry axes. Note that qx, qy, qz are

reduced wavevector with unit 2π/a.

the other way around. This implies the curvature depends on the location of the critical point. Therefore, it will cause the numbers of points in aggregation are much more in somewhere. This problem can be simply resolved by dynamically adjusting the tolerance value to limit the number of points in aggregation with a acceptable threshold value (in our analysis the value is 10).

The distribution of the critical points are also interesting. In Fig.3.10(a), the critical point responsible for the features of GaAs are mainly located around symmetry point L,W,X and Γ. Some of them lie on Σ, ∆, Q, Λ and S axis. The others reside on symmetry faces ΓKL, ΓXWK and ΓXUL. As for InP, the distribution is roughly similiar to that of GaAs since they are both zincblende structure and with similar profile of dispersion manifold.

3.3 Conclusion

We have performed a convincing analysis of far-infrared features resulting from two-phonon processes for GaAs and InP. The phonon dispersions have been calculated exactly by DFPT. The resulting 2PDOS spectra can successfully explained many observed features from either

(45)

3.3 Conclusion 34 infrared absorption or Raman scattering measurement, and further give a prediction beyond the limit of experiment. We have also identified the features by assigning them the corresponding two-branch combinations, the critical points, and the types of critical points. The analysis indicates that phonons of wave vectors around symmetry points and along symmetry lines are collectively responsible for strong infrared features in two-phonon processes.

(46)

Appendix

(47)

A.1 Individual two phonon density of states for GaAs 36

A.1 Individual two phonon density of states for GaAs

A.1.1 Sum process

0 100 200 300 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

2ν1 0 100 200 300 0 0.005 0.01 0.015 0.02 Frequency (cm−1)

ν1+ ν2 0 200 400 0 0.01 0.02 0.03 Frequency (cm−1)

ν1+ ν3 0 200 400 0 0.02 0.04 0.06 Frequency (cm−1)

ν1+ ν4 0 200 400 0 0.05 0.1 Frequency (cm−1)

ν1+ ν5 0 200 400 0 0.02 0.04 0.06 Frequency (cm−1)

ν1+ ν6 0 100 200 300 0 0.005 0.01 0.015 0.02 Frequency (cm−1)

2ν2 0 200 400 0 0.005 0.01 0.015 Frequency (cm−1)

ν2+ ν3 0 200 400 0 0.05 0.1 Frequency (cm−1)

ν2+ ν4

Figure A.1: Individual two phonon density of states for GaAs with branches combinations 2ν1, ν1+ ν2,

(48)

A.1 Individual two phonon density of states for GaAs 37 0 200 400 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν2+ ν5 0 200 400 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν2+ ν6 0 200 400 600 0 0.005 0.01 0.015 0.02 Frequency (cm−1)

2ν3 0 200 400 600 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν3+ ν4 0 200 400 600 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν3+ ν5 0 200 400 600 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν3+ ν6 0 200 400 600 0 0.005 0.01 0.015 0.02 Frequency (cm−1)

2ν4 0 200 400 600 0 0.01 0.02 0.03 Frequency (cm−1)

ν4+ ν5 0 200 400 600 0 0.01 0.02 0.03 Frequency (cm−1)

ν4+ ν6 0 200 400 600 0 0.1 0.2 0.3 0.4 Frequency (cm−1)

2ν5 0 200 400 600 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

ν5+ ν6 0 200 400 600 0 0.05 0.1 Frequency (cm−1)

2ν6

Figure A.2: Individual two phonon density of states for GaAs with branched combinations ν2 + ν5,

(49)

A.1 Individual two phonon density of states for GaAs 38

A.1.2 Difference process

0 20 40 0 0.02 0.04 0.06 Frequency (cm−1)

ν2− ν1 0 50 100 150 0 0.01 0.02 0.03 Frequency (cm−1)

ν3− ν1 0 100 200 300 0 0.005 0.01 0.015 0.02 Frequency (cm−1)

ν4− ν1 0 100 200 300 0 0.02 0.04 0.06 Frequency (cm−1)

ν5− ν1 0 100 200 300 0 0.01 0.02 0.03 Frequency (cm−1)

ν6− ν1 0 50 100 150 0 0.02 0.04 0.06 Frequency (cm−1)

ν3− ν2 0 100 200 300 0 0.005 0.01 0.015 Frequency (cm−1)

ν4− ν2 0 100 200 300 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν5− ν2 0 100 200 300 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν6− ν2

Figure A.3: Individual two phonon density of states for GaAs with branched combinations ν2 − ν1,

(50)

A.1 Individual two phonon density of states for GaAs 39 0 100 200 300 0 0.01 0.02 0.03 Frequency (cm−1)

ν4− ν3 0 100 200 300 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν5− ν3 0 100 200 300 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν6− ν3 0 20 40 60 0 0.05 0.1 Frequency (cm−1)

ν5− ν4 0 10 20 30 40 50 60 0 0.01 0.02 0.03 0.04 0.05 0.06 Frequency (cm−1₎

ν6− ν4 0 10 20 30 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

ν6− ν5

Figure A.4: Individual two phonon density of states for GaAs with branched combinations ν4 − ν3,

(51)

A.2 Individual two phonon density of states for InP 40

A.2 Individual two phonon density of states for InP

A.2.1 Sum process

0 100 200 0 0.02 0.04 0.06 Frequency (cm−1)

2ν1 0 100 200 0 0.01 0.02 0.03 Frequency (cm−1)

ν1+ ν2 0 100 200 300 0 0.02 0.04 0.06 Frequency (cm−1)

ν1+ ν3 0 200 400 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

ν1+ ν4 0 200 400 600 0 0.02 0.04 0.06 Frequency (cm−1)

ν1+ ν5 0 200 400 600 0 0.02 0.04 0.06 Frequency (cm−1)

ν1+ ν6 0 100 200 0 0.01 0.02 0.03 Frequency (cm−1)

2ν2 0 100 200 300 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν2+ ν3 0 200 400 600 0 0.05 0.1 Frequency (cm−1)

ν2+ ν4

Figure A.5: Individual two phonon density of states for InP with branches combinations 2ν1, ν1+ ν2,

(52)

A.2 Individual two phonon density of states for InP 41 0 200 400 600 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν2+ ν5 0 200 400 600 0 0.01 0.02 0.03 Frequency (cm−1)

ν2+ ν6 0 200 400 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

2ν3 0 200 400 600 0 0.02 0.04 0.06 Frequency (cm−1)

ν3+ ν4 0 200 400 600 0 0.02 0.04 0.06 Frequency (cm−1)

ν3+ ν5 0 200 400 600 0 0.05 0.1 Frequency (cm−1)

ν3+ ν6 0 500 1000 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

2ν4 0 500 1000 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

ν4+ ν5 0 500 1000 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

ν4+ ν6 0 500 1000 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

2ν5 0 500 1000 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

ν5+ ν6 0 500 1000 0 0.05 0.1 Frequency (cm−1)

2ν6

Figure A.6: Individual two phonon density of states for InP with branched combinations ν2+ν5, ν2+ν6,

(53)

A.2 Individual two phonon density of states for InP 42

A.2.2 Difference process

0 10 20 30

0 0.05 0.1

Frequency (cm−1)

ν2− ν1 0 50 100 150 0 0.02 0.04 0.06 Frequency (cm−1)

ν3− ν1 0 200 400 0 0.01 0.02 0.03 0.04 Frequency (cm−1)

ν4− ν1 0 200 400 0 0.05 0.1 Frequency (cm−1)

ν5− ν1 0 200 400 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

ν6− ν1 0 50 100 150 0 0.02 0.04 0.06 Frequency (cm−1)

ν3− ν2 0 200 400 0 0.01 0.02 0.03 Frequency (cm−1)

ν4− ν2 0 200 400 0 0.02 0.04 0.06 Frequency (cm−1)

ν5− ν2 0 200 400 0 0.05 0.1 Frequency (cm−1)

ν6− ν2

Figure A.7: Individual two phonon density of states for InP with branched combinations ν2−ν1, ν3−ν1,

(54)

A.2 Individual two phonon density of states for InP 43 0 200 400 0 0.02 0.04 0.06 Frequency (cm−1)

ν4− ν3 0 200 400 0 0.05 0.1 Frequency (cm−1)

ν5− ν3 0 200 400 0 0.02 0.04 0.06 Frequency (cm−1)

ν6− ν3 0 10 20 30 0 0.1 0.2 0.3 0.4 Frequency (cm−1)

ν5− ν4 0 20 40 60 0 0.05 0.1 Frequency (cm−1)

ν6− ν4 0 20 40 0 0.05 0.1 0.15 0.2 Frequency (cm−1)

ν6− ν5

Figure A.8: Individual two phonon density of states for InP with branched combinations ν4−ν3, ν5−ν3,

(55)

A.3 Spectral feature assignment for GaAs 44

A.3 Spectral feature assignment for GaAs

Feature type Position (cm−1) Branch combination Critical points P0 9.8 ν4− ν3 W(0.5,1,0) P 12.4 ν2− ν1 W(0.5,1,0) P 13.5 ν6− ν5 L(0.5,0.5,0.5) P2 15.2 ν6− ν5 W(0.5,1,0) P 19.41 ν2− ν1 Q(0.5,0.77,0.23) P 20.12 ν6− ν5 Σ(0.59,0.59,0) P 20.4 ν6− ν4 Λ(0.05,0.05,0.05) P2 23.2 ν6− ν4 ΓXWK(0.22,0.48,0) P 26.3 ν5− ν4 L(0.5,0.5,0.5) P3 27 ν6− ν4 L(0.5,0.5,0.5) P0 27.1 ν4− ν3 L(0.5,0.5,0.5) P3 27.4 ν6− ν5 ∆(0,0.43,0) P0 30 ν5,6− ν3 X(0,1,0) P3 35.7 ν2− ν1 Σ(0.67,0.67,0) P2 37.9 ν5− ν4 Σ(0.71,0.71,0) P 42.7 ν4− ν3 Q(0.5,0.73,0.27) P3 46.2 ν5− ν4 W(0.5,1,0) P0 51.5 ν5,6− ν3 L(0.5,0.5,0.5) P2 53.6 ν6− ν4 Σ(0.69,0.69,0) P1 57.2 ν5− ν3 W(0.5,1,0) P3 59.2 ν6− ν4 W(0.5,1,0)

(56)

A.3 Spectral feature assignment for GaAs 45 continued from previous page

Feature type Position (cm−1) Branch combination Critical points P1 68.2 ν6− ν3 W(0.5,1,0) P 69.9 ν5− ν3 Q(0.5,0.76,0.24) P 74.5 ν3− ν2 Q(0.5,0.84,0.16) P2 81.4 ν3− ν2 W(0.5,1,0) P 83.6 ν6− ν3 Q(0.5,0.78,0.22) P0 92.1 ν4− ν2 W(0.5,1,0) P 93.8 ν3− ν1 W(0.5,1,0) P0 98.1 ν4− ν2 Σ(0.7,0.7,0) P0 104.5 ν4− ν1 W(0.5,1,0) P2 106.1 ν3− ν1 ΓKL(0.66,0.66,0.18) P1 128.1 2ν1 L(0.5,0.5,0.5) P1 130.2 ν1+ ν2 L(0.5,0.5,0.5) P1 132.8 2ν2 L(0.5,0.5,0.5) P1 133 ν4− ν1 Σ(0.72,0.72,0) P0 136 ν5− ν2 Σ(0.68,0.68,0) P 138.2 ν5− ν2 W(0.5,1,0) P3 148.5 ν3− ν1,2 L(0.5,0.5,0.5) P0 150.8 ν5− ν1 W(0.5,1,0) P0 151 ν6− ν2 W(0.5,1,0) P 151.8 ν6− ν2 Σ(0.72,0.72,0) P 161.5 ν4− ν1,2 X(0,1,0) P0 163.8 ν6− ν1 W(0.5,1,0) P 165.4 2ν1 S(0.14,1,0.14) P1 165.9 ν1,2+ ν2 X(0,1,0)

(57)

Feature type Position (cm−1) Branch combination Critical points P2 171 ν5,6− ν2 X(0,1,0) P 171.6 ν5− ν1 S(0.16,1,0.16) P2 171.6 ν4− ν2 ∆(0,0.8,0) P 172.6 ν4− ν1 L(0.5,0.5,0.5) P3 172.6 ν4− ν2 L(0.5,0.5,0.5) P2 198 ν5,6− ν2 L(0.5,0.5,0.5) P 201.1 ν6− ν1 L(0.5,0.5,0.5) P2 201.1 ν5− ν1 L(0.5,0.5,0.5) P 201.1 ν1+ ν2 Σ(0.67,0.67,0) P3 207.1 2ν1 W(0.5,1,0) P3 219.5 ν1+ ν2 W(0.5,1,0) P 232 2ν2 W(0.5,1,0) P3 237.2 2ν2 Σ(0.67,0.67,0) P 265.8 ν1+ ν3 ΓKL(0.62,0.62,0.25) P3 278 ν1+ ν3 L(0.5,0.5,0.5) P1 278 ν2+ ν3 L(0.5,0.5,0.5) P0 299 ν1+ ν4 L(0.5,0.5,0.5) P2 301.2 ν1+ ν3 W(0.5,1,0) P1 302.4 ν1+ ν4 L(0.5,0.5,0.5) P0 302.4 ν2+ ν4 L(0.5,0.5,0.5) P3 306.6 ν1+ ν3 X(0,1,0) P 306.6 ν2+ ν3 X(0,1,0) P1 307.5 ν1+ ν4 ΓKL(0.61,0.61,0.21) P2 312.9 ν2+ ν3 W(0.5,1,0)

(58)

Feature type Position (cm−1) Branch combination Critical points P 313.2 ν1+ ν4 W(0.5,1,0) P 319.1 ν1+ ν4 Σ(0.38,0.38,0) P3 319.5 ν2+ ν3 K(0.75,0.75,0) P2 322.8 ν1+ ν4 ΓXUL(0.26,0.36,0.26) P0 326.6 ν2+ ν4 ΓXUL(0.31,0.35,0.31) P1 330.5 ν1+ ν5,6, ν2+ ν6 L(0.5,0.5,0.5) P1 335.5 ν1+ ν6 Λ(0.24,0.24,0.24) P3 335.9 ν1+ ν4 ∆(0,0.83,0) P1 335.9 ν2+ ν4 ∆(0,0.83,0) P1 337 ν2+ ν6 Λ(0.23,0.23,0.23) P 337.1 ν1+ ν5 X(0,1,0) P1 337.6 ν2+ ν5 ∆(0,0.82,0) P 338.4 ν2+ ν4 ΓKL(0.66,0.66,0.13) P2 342.4 ν2+ ν4 ΓXWK(0.33,0.64,0) P3 346.4 ν2+ ν4 ΓXUL(0.2,0.69,0.2) P1 349.2 ν1+ ν6 ΓXUL(0.19,0.7,0.19) P 352 ν1+ ν6 Σ(0.57,0.57,0) P 354 ν1+ ν6 ∆(0,0.52,0) P3 357.7 ν1+ ν5 W(0.5,1,0) P 358.8 ν2+ ν6 ∆(0,0.52,0) P 370.4 ν2+ ν5 W(0.5,1,0) P3 370.7 ν1+ ν6 W(0.5,1,0) P 372.2 ν3+ ν3 Q(0.5,0.78,0.22) P3 373.4 ν2+ ν5 Σ(0.66,0.66,0)

(59)

Feature type Position (cm−1) Branch combination Critical points

P 384 ν2+ ν6 W(0.5,1,0) P3 389.7 ν2+ ν6 Σ(0.64,0.64,0) P2 397.5 3ν3 W(0.5,1,0) P 406.8 ν3+ ν4 W(0.5,1,0) P2 415.6 ν3+ ν4 ΓXUL(0.27,0.72,0.72) P0 416.1 2ν4 W(0.5,1,0) P3 428 2ν3 L(0.5,0.5,0.5) P1 433 2ν4 Σ(0.71,0.71,0) P 442.5 ν3+ ν5 Q(0.5,0.8,0.2) P3 449 ν3+ ν4 L(0.5,0.5,0.5) P2 451.6 ν3+ ν5 W(0.5,1,0) P 455.6 ν3+ ν6 Q(0.5,0.78,0.22) P2 457.6 ν3+ ν6 ΓKL(0.64,0.64,0.19) P0 462 ν4+ ν5 W(0.5,1,0) P2 465.5 ν3+ ν6 W(0.5,1,0) P3 467.8 ν3+ ν4 X(0,1,0) P0 470 ν4+ ν5 Σ(0.72,0.72,0) P0 475.3 ν4+ ν6 W(0.5,1,0) P 475.9 2ν4 L(0.5,0.5,0.5) P3 477.8 ν3+ ν5,6 L(0.5,0.5,0.5) P1 484 ν4+ ν6 Σ(0.73,0.73,0) P 496.1 ν4+ ν5 X(0,1,0) P1 496.8 ν4+ ν6 X(0,1,0) P 501.5 ν4+ ν5,6 L(0.5,0.5,0.5)

(60)

Feature type Position (cm−1) Branch combination Critical points P0 507.1 2ν5 S(0.09,1,0.09) P0 507.5 ν5,6+ ν6 X(0,1,0) P 508.5 ν4,5+ ν5 ∆(0,0.8,0) P 509.2 2ν4 ∆(0,0.73,0) P1 521.5 ν5 + ν6 W(0.5,1,0) P 530 ν5 + ν6 L(0.5,0.5,0.5) P2 530 2ν5 L(0.5,0.5,0.5) P0 530 2ν6 L(0.5,0.5,0.5) P2 532 ν4+ ν4,5 Λ(0.31,0.31,0.31) P1 534.5 2ν6 W(0.5,1,0) P 538.8 ν6 + ν6 Q(0.5,0.75,0.25) P3 549.9 ν4,5+ ν4,5 Γ(0,0,0) P3 570 ν4,5+ ν6 Γ(0,0,0) P3 590.1 ν6 + ν6 Γ(0,0,0)

(61)

A.4 Spectral feature assignment for InP 50

A.4 Spectral feature assignment for InP

Feature type Position (cm−1) Branch combination Critical points P 0.52 ν5− ν4 Q(0.5,0.54,0.46) P3 2.79 ν5− ν4 Σ(0.32,0.32,0) P2 10.4 ν5− ν4 U(0.25,1,0.25) P0 13 ν6− ν4,5 L(0.5,0.5,0.5) P 16.66 ν2− ν1 W(0.5,1,0) P1 18.08 ν2− ν1 Q(0.5,0.78,0.22) P1 22.3 ν6− ν5 W(0.5,1,0) P3 22.6 ν5− ν4 W(0.5,1,0) P1 24.23 ν6− ν5 Q(0.5,0.79,0.21) P3 24.43 ν2− ν1 Σ(0.65,0.65,0) P1 26.2 ν6− ν4 ΓXUL(0.18,0.52,0.18) P 27.01 ν6− ν5 Σ(0.33,0.33,0) P 29.2 ν6− ν5 Σ(0.33,0.33,0) P3 31.3 ν6− ν5 Σ(0.53,0.53,0) P3 36.1 ν6− ν4,5 Γ(0,0,0) P 38.3 ν6− ν4 Σ(0.65,0.65,0) P3 44.9 ν6− ν4 W(0.5,1,0) P 62.7 ν3− ν2 W(0.5,1,0) P 67.41 ν3− ν2 ΓXUL(0.27,0.73,0.27) P 79.09 ν3− ν1 W(0.5,1,0) P1 89.92 ν3− ν1 ΓXUL(0.31,0.78,0.31)

以第一原理計算分析砷化鎵與磷化銦 紅外線雙聲子特徵譜線之研究

國 立 交 通 大 學

電子工程學系 電子研究所

碩 士 論 文

以第一原理計算分析砷化鎵與磷化銦

紅外線雙聲子特徵譜線之研究

Analysis of two-phonon infrared spectral features

of gallium arsenide and indium phosphide by ab

initio calculation

研 究 生：林煌翔

Student: Huang-Hsiang Lin

指導教授：顏順通 博士

Advisor: Dr. Shun-Tung Yen

指導教授：

Advisor:

以第一原理計算分析砷化鎵與磷化銦

紅外線雙聲子特徵譜線之研究

Analysis of two-phonon infrared spectral features of

gallium arsenide and indium phosphide by ab initio

calculation

研 究 生：林煌翔

Student: Huang-Hsiang Lin

指導教授：顏順通 博士

Advisor: Dr. Shun-Tung Yen

指導教授：

國立交通大學

電子工程學系 電子研究所

碩士論文

A Thesis

Submitted to Department of Electronics Engineering and

Institute of Electronic

College of Electrical and Computer Engineering

National Chiao Tung University

in Partial Fulfillment of the Requirements

for the Degree of

Master of Science

in

Electronic Engineering

September 2014

Hsinchu, Taiwan, Republic of China

以第一原理計算分析砷化鎵與磷化銦紅外線特徵

譜線之研究

學生：林煌翔

指導教授： 顏順通博士

國立交通大學

電子工程學系 電子研究所 碩士班

摘要

Analysis of two-phonon infrared spectral features of gallium

arsenide and indium phosphide by ab initio calculations

Student: Huang-Hsiang, Lin

Advisor: Shun-Tung, Yen

Department of Electronics Engineering

Institute of Electronics

National Chiao Tung University

ABSTRACT

Acknowledgments

Contents

List of Figures

List of Tables

C H A P T E R

1

Introduction

C H A P T E R

2

Theoretical Approach

2.1

Density funtional theory

2.1.1

Kohn-Sham equations

2.1.2

Local density approximation

2.2

Density functional perturbation theory

2.3

Phonons in DFPT

2.3.1

Interatomic force constant

2.3.2

Dynamical matrix

2.4

以第一原理計算分析砷化鎵與磷化銦紅外線雙聲子特徵譜線之研究

國立交通大學

電子工程學系電子研究所

碩士論文

研究生：林煌翔

指導教授：顏順通博士

研究生：林煌翔

指導教授：顏順通博士

電子工程學系電子研究所

指導教授：顏順通博士

電子工程學系電子研究所碩士班

_(q

₎

_(q

₎

_(q

₎

_(q

₎