Time-resolved dynamical Franz-Keldysh effect under elliptically
polarized laser
T. Otobe
arXiv:1609.07280v1 [physics.optics] 23 Sep 2016
Kansai Photon Science Institute, National Institutes for Quantum and
Radiological Science and Technology (QST), Kyoto 619-0215, Japan
Abstract
The analytical formula for the time-resolved dynamical Franz-Keldysh effect (Tr-DFKE) under
an elliptically polarized laser in sub-femtosecond time-scale is reported. The Houston function
is assumed as the time-dependent wave function of the parabolic two-band system. The resulting formula shows the sub-cycle change of the optical properties for elliptically polarization; the
modulation of the dielectric function becomes smaller than that of linear polarization. On the
other hand, the subcycle modulation of the dielectric function disappears for a circularly polarized
laser, which is a significant feature of the Tr-DFKE. This analytical formulas show good qualitative agreement with the first-principle calculation employing the time-dependent density functional
theory for diamond.
1
I.
INTRODUCTION
In the last two decades, advances in laser sciences and technologies have led to the
availability of intense coherent light with different characteristics. Ultra-short laser pulses
can be as short as a few tens of an attoseconds, leading to the new field of attosecond
science [1]. Intense laser pulses of mid-infrared (MIR) or terahertz (THz) frequencies have
also recently become available [2, 3]. By employing these extreme sources of coherent light,
it is possible to investigate the optical response of materials in real time with a resolution
much lower than an optical cycle [1, 4–8].
The dielectric function ε(ω) is the most fundamental quantity characterizing the optical
properties of matter. The dielectric function observed in an ultra-fast pump-probe experiment can be further considered as a probe time (Tp )-dependent function, ε(Tp , ω). The
modulation of the dielectric function ε(ω) in the presence of electromagnetic fields has been
a subject of investigation for many years. The change in the dielectric function under a static
electric field is known as the Franz-Keldysh effect (FKE) [9–16], and the change under an
alternating electric field is known as the dynamical FKE (DFKE) [17–22].
We determined the sub-cycle change in the optical properties, i.e., the time-resolved
DFKE (Tr-DFKE), which corresponds to the response of the dressed states and quantum
path interference of different dressed state [23]. In particular, this ultra-fast change exhibits
an interesting phase shift that depends on the field amplitude and probe frequency. By
utilizing this phenomenon, we can develop an ultra-fast modulator of light or an ultra-fast
optical switch. Recently, the Tr-DFKE was observed experimentally by a near-infrared
(NIR)-pump extreme ultraviolet (EUV)-probe with attosecond time-resolution for polycrystalline diamond [24]. A similar effect was also observed in an excitonic state in a quantum
well of GaAs by a THz-pump NIR-probe [25].
We previously showed that the field amplitude of the pump laser and the bandwidth of
the probe pulse are crucial parameters for controlling the Tr-DFKE [23]. Another possible
control parameter is the polarization of the laser. In the current work, we present analytical
formulas and first-principle calculations for the Tr-DFKE under elliptically and circularly
polarized pump lasers. The time-dependent oscillation of the optical properties decreases as
the ellipticity increases. We found that the time-dependent change in the optical properties
completely disappears under a circularly polarized pump laser for an isotropic band struc2
ture. We compare our formulas with first-principle calculations employing time-dependent
density-functional theory (TDDFT) for diamond.
The remainder of this paper is organized as follows. In Sec. II, we develop an analytical
formulation for the Tr- DFKE under a circularly polarized light using a parabolic two-band
model. In Sec. III, we present a real-time TDDFT approach for the Tr-DFKE. In Sec. IV,
we present the numerical results. In Sec. V, a summary is given.
II.
FORMULATION
To derive the time-dependent conductivity, we will revisit a simple model that we reported
in our previous work [23]. The probe’s electric field is assumed to be weak enough that is
can be treated using the linear response theory. We denote the electric current caused by the
probe field as Jp (t), which is assumed to be parallel to the direction of the probe’s electric
field. the relationship to the time-domain conductivity σ(t, t′ ) is given as:
Z t
Jp (t) =
dt′ σ(t, t′ )Ep (t′ ),
(1)
−∞
where Ep (t′ ) is the electric field of the probe pulse. We note that the conductivity σ(t, t′ )
depends on both times t and t′ rather than the just the time difference t − t′ , owing to the
presence of the pump pulse. If the probe laser is much shorter than the optical cycle of the
pump laser and has peaks at time t = Tp , we can define the time-dependent conductivity,
σ̃(Tp , ω), as:
σ̃(Tp , ω) =
R
dteiωt Jp (t)
R
.
dtEp (t)
(2)
In the following, we consider a simplified description: electron dynamics in the presence
of the pump and probe fields are assumed to be described by a time-dependent Schrödinger
equation for a single electron:
e ~ 2
1
∂
p~ + A(t) + V (~r) ψn (~r, t),
i ψn (~r, t) =
∂t
2me
c
(3)
~
where ψn (~r, t) is the time-dependent wave function of the n-th band, A(t)
is the vector
potential of the pump light field, and V (~r) is a time-independent, lattice periodic potential.
In this paper, we employ atomic units for all equations. We express the solution of this
P ~
equation using the time-dependent Bloch function vn~k (~r, t) as ψn (~r, t) = ~k eik~r vn~k (~r, t),
where ~k is the Bloch wavenumber vector.
3
We further assume that in the presence of the pump field described by a vector potential
~ P (t), the solution of Eq. (3) is well approximated by the Houston function [17, 26]. Using
A
static Bloch orbitals un~k (~r) and orbital energies ǫn~k that satisfy:
2
1
p~ + ~k + V (~r) un~k (~r) = ǫn~k un~k (~r),
2me
(4)
the Houston function can be expressed as:
wn~k (~r, t) = un ~kP (t) (~r) exp −i
Z
t
′
ǫn ~kP (t′ ) dt ,
(5)
~ P (t)/c.
where ~kP (t) is defined by ~kP (t) = ~k + eA
We now consider an elliptically polarized pump laser with the vector potential:
~ P (t) = A0 (η sin Ωt, 0, cos Ωt) (0 ≤ η ≤ 1),
A
(6)
and a two-band model, considering of only the two orbitals in the occupied valence (v) and
unoccupied conduction (c) bands. Here, η is the ellipticity of the pump laser. The excitation
energy from the valence band to the conduction band is assumed to have a parabolic form:
ǫc~k − ǫv~k ≃
k2
+ ǫg ≡ ǫk + ǫg ,
2µ
(7)
where ǫg is the band gap energy and µ is the reduced mass of the electron-hole pairs. We
assume the probe is a δ-functional linear polarized laser, fp δ(t − Tp ). Our formulas do not
depend on the polarization of probe laser.
Using the same procedure as in our previous work [23], we obtain the time-resolved
conductivity:
Z
Z 1
Z 2π
4πe2
2e2 |pcv |2 µ3/2 ∞ √
εE (Tp , ω) = 1 −
ne − √
εk dεk
d cos θ
dφ
mω 2
2m2e π
−1
0
0
X
×
(−1)m−n J˜l1 (α, β)J˜l1 −2m (α, β)Jl2 (γ)Jl +2(m−n) (γ)
2
l1 ,l2 ,m,n
ei2nΩTp
(ω + 2nΩ)(ω − (εg + εk + UE + (l1 + l2 − 2n)Ω))
e−i2nΩTp
,
−
(ω − 2nΩ)(ω + (εg + εk + UE + (l1 + l2 − 2n)Ω))
×
where:
ekA0
cos θ,
cµΩ
e2 A2
β = 2 0 (1 − η 2 ),
8c µΩ
ekA0
γ=η
sin θ cos φ,
cµΩ
α=
4
(8)
(9)
(10)
(11)
and:
e2 A20
(1 + η 2 ).
(12)
4c2 µ
Here Ue is the average kinetic energy in the elliptic polarized laser field. Here, Jl is the lth
Bessel function and J˜l is the generalized Bessel function [27] defined as:
X
J˜l (α, β) =
Jl−2i (α)Ji (β).
(13)
Ue =
i
We define the ~k by the polar angles, such that ~k = k(sin θ cos φ, sin θ sin φ, cos θ).
Equation (8) is connected to the Tr-DFKE under a linearly polarized laser at η = 0,
which we reported in our previous paper [23]:
Z 1
Z
X
√
4πe2
4e2 |pcv |2 µ3/2 ∞
√
ǫ
d(cos
θ)
J˜l (α, β)J˜l−2n (α, β)
εL (Tp , ω) = 1 −
n
−
dǫ
k
e
k
me ω 2
2πm2e ω
−1
0
l,n
i2nΩTp
e
×
(ω + 2nΩ){ω − (ǫg + ǫk + Up + (l − 2n)Ω)}
e−i2nΩTp
,
(14)
−
(ω − 2nΩ){ω + (ǫg + ǫk + Up + (l − 2n)Ω)}
where Up =
e2 A20
4µc2
is the well known ponderomotive energy.
The form of the circularly polarized laser (η = 1) becomes:
Z 1
X 4e2 |pcv |2 µ3/2 Z ∞ √
4πe2
√
εk dεk
d cos θJl2 (α1 )
εC (Tp , ω) = 1 −
ne −
2ω2
me ω 2
2πm
−1
0
e
l
1
1
,
−
×
ω − (εg + εk + Uc + lΩ) ω + (εg + εk + Uc + lΩ)
where:
ekA0 sin θ
α1 =
µcΩ
and:
e2 A20
.
Uc =
2µc2
In this case, we redefine the polarization as:
~ P (t) = A0 (cos Ωt, sin Ωt, 0),
A
(15)
(16)
(17)
(18)
to have a simpler formula. An interesting point is that the time-dependence of Tr-DFKE
disappears for a circularly polarized laser since the summation over ~k corresponds to the
time-average for one ~k-point isotropic system. We note that εC does not depend on the
parameter β that corresponds to the parameter γ = Up /Ω for the linear polarization [19].
Above formulas do not depend on the polarization direction of the probe laser since the
isotropic band structure is assumed.
5
III.
REAL-TIME TDDFT
In the above section, we derived the analytical forms for Tr-DFKE that can be used
for real materials. Phase-sensitive modulation attributed to the Tr-DFKE under a linearly
polarized laser has been reported by Lucchuni et al. [24] which agree with our previous work
[23]. This result indicates that real-time calculation based on the TDDFT is a the reliable
approach.
In this section, we revisit our first-principle approach. The detail of the computational
methods have been already reported elsewhere [23, 29, 30]. In real-time TDDFT, we describe
electron dynamics in a unit cell of a crystalline solid under a spatially uniform electric field
Rt
~
~ ′ ), the electron dynamics
E(t). Treating the field by a vector potential A(t)
= −c dt′ E(t
are described by the time-dependent Kohn-Sham (TDKS) equation [28]:
e ~ 2
1
∂
~p + A(t) + V (~r, t) ψi (~r, t),
i ψi (~r, t) =
∂t
2me
c
(19)
where V (~r, t) is composed of electron-ion, the Hartree, and exchange-correlation potentials. We use a norm-conserving pseudopotential for the electron-ion potential [31, 32].
For the exchange-correlation potential, we employ an adiabatic local density approximation
(LDA)[33].
We calculate electron dynamics in diamond, using a cubic unit cell containing eight carbon
atoms. The TDKS equation is solved in real time and real space. A total of 203 real-space
grids are used for the unit cell, and 323 grids are used for the k-points. The Taylor expansion
method is used for the time evolution [34] with a time step of ∆t = 0.02 in atomic unit.
The number of time steps is typically 75,000. An important output of the calculation is the
average electric current density as a function of time, J(t), given by:
Z
X
e
e~
J(t) = −
d~r
Reψi∗ ~p + A(t)
ψi + JN L (t),
me V V
c
i
(20)
where V is the volume of the unit cell and JN L (t) is the current caused by non-locality of
the pseudopotential.
In practice, we use the following electric fields. The pump field is of the form:
~ P (t) = E0,P fP (t)(η cos Ωt, 0, sin Ωt),
E
(21)
with a central angular frequency of Ω set to Ω = 0.4 eV. The field is turned on adiabatically
6
by fP (t), given by fP (t) = sin2
fs.
π
2TP
t for 0 < t < TP and fP (t) = 1 for t ≥ TP and TP =10
The probe field is of the form:
~ p (t) = ŷE0,p sin(ωp t) exp −(t − Tp )2 /ζ 2 .
E
(22)
To improve the numerical accuracy, we assume a linear polarized probe pulse laser that
is vertical to the pump laser field. The average frequency ωp is set to 5.5 eV, which is
equal to the calculated band gap energy of diamond in LDA. The field strength is set to
E0,p = 2.7×10−3 MV/cm, which is small enough to probe the linear response of the medium.
The pulse duration ζ is set to ζ = 0.7 fs.
The frequency-dependent conductivity is calculated as the ratio of the Fourier-transformed
current Jp (t) and field Ep (t), given by:
σ(Tp , ω) =
R
dteiωt g(t − Tp )Jp (t)
R
,
dteiωt Ep (t)
(23)
where g(t) is a filter function to damp the current and suppress spurious oscillations arising
from the numerical cutoff in the integration. The time-dependent conductivity is converted
to the dielectric function using the formula:
ǫ(Tp , ω) = 1 + i
IV.
4π
σ(Tp , ω).
ω
(24)
NUMERICAL RESULTS
In this section, we aim to demonstrate the quantitative difference between circular and
linear polarization by comparing the numerical results. Diamond is a typical dielectric used
in non-linear laser-matter interaction studies, and we selected it as an example with which
to illustrate the application of the above formalism.
Figure 1 shows the change in the imaginary part of the dielectric function, ∆Im[ε(Tp , ω)].
The maximum field amplitude is 10 MV/cm and the Ω is 0.4 eV. Figure 1 (b)-(d) show the
results from the real-time TDDFT, and (e)-(g) show the results of the analytical theory. We
assumed that the band gap, εg , is 5.5 eV, the reduced mass µ is 0.25me , and the transition
moment Pcv = 0.928 in atomic units for the analytical formulas. We adopt time components
of n = −1, 0, and 1 in Eqs. (8), (14), and (15). The vertical black solid lines in Fig. 1(b)-(g)
7
7654302*+/10,/.-,+*)
"!'"
"!"$
"!""
("!"$
("!'"
3?7
'&
AB:CDE
3@7
3,7
!"
'%
("!& ("!' "!"
'#
"!'
'8
"!&
3+7
#!$
#!"
$!$
$!"
$!"
%!$
%!$
'&
'#
&"
&%
#!$
317
!"
#!"
$!$
$!$
$!"
$!"
%!$
%!$
'&
'#
&"
&%
#!$
("!'
"!"
&%
"!'
"!&
'&
'#
&"
&%
'&
'#
&"
&%
'&
'#
&"
&%
#!$
#!"
!"
("!&
&&
#!$
$!$
74+30>=.+<)0
327
!"
#!"
!"
&"
9/:+031;7
3=7
!"
#!$
#!"
#!"
$!$
$!$
$!"
$!"
%!$
%!$
'&
'#
&"
9/:+031;7
&%
9/:+031;7
FIG. 1. Changes in Im[ε(Tp , ω)] ( ∆Im[ε(Tp , ω)]) under circularly ((b) and (e)), elliptically ((c)
and (f)), and linearly ((d) and (g)) polarized lasers. (b)-(d) show the numerical results from the
real-time TDDFT, and (e)-(g) show the results from analytical formulas Eqs. (8), (14), and (15),
respectively. Panel (a) shows the electric field of the pump laser. The horizontal axis presents the
probe time, whose range corresponds to the simulation by the real-time TDDFT.
represent the maximum of the absolute value of the electric field shown in Fig. 1 (a). In
numerical simulation by the real-time TDDFT, we assume the pump field rotating in x-z
plane, which corresponds to a- and c-axises of diamond, and the polarization of the probe
pulse is set to be parallel to b-axis.
The Tr-DFKE under circular polarization is presented in Fig. 1 (b) and (e) respectively. As we expected in Eq. (15), the time-dependence of DFKE almost disappears in
real-time TDDFT. However, for higher energy region around 6.5 − 7eV, small oscillation
of ∆Im[ε(Tp , ω)] can be seen. The ratio of the small oscillation between the peak and dip
is approximately 10 %. There are two possible explanations for this modulation. First, the
artificial error in real-time simulation is unavoidable. The Tr-DFKE requires high compu8
!"#
(9:3;<=>3(?@="
(0==:AB:;(?@="
(C:12>3(?@="
!"$
(DB>B:;
+*)('&%
!"!
!"$
!"#
/"!
/".
."!
.".
-"!
-".
,"!
012345(6278
FIG. 2. Time-average of the ∆Im[ε(Tp , ω)] under a laser field with frequency of 0.4 eV and a
field amplitude 10 MV/cm calculated from Eqs. (8), (14), and (15). For reference, the static FKE
calculated by the Tharmalingam formula [11] is represented by the black dashed line.
tational accuracy. Therefore, small errors in time evolution and the Fourier transformation
affect the result. The second is the anisotropy of the band structure, which changes the
definition of the reduced mass µ → µ(k, θ, φ). Although the band structure of diamond
is almost isotropic around the Γ-point, it has a four-fold symmetric structure far from the
Γ-point and higher energy bands. Since the DFKE depends on the reduced mass, timedependence in circular polarization should include the effect of the structure of the valence
and conduction bands.
Figures 1 (c) and (f) show the case of the elliptically polarized laser. The real-time
TDDFT shows 2Ω oscillation for all energies. The analytic formula Eq. (8) well reproduces the real-time TDDFT qualitatively for not only the frequency of the oscillation of
∆In[ε(Tp , ω)], but also the phase shift depending on the probe photon energy. The amount
of oscillation and the phase shift are weaker than that of the linear polarization shown in
Figs. 1 (d) and (g). However, the peaks of the modulation shows the same phase shift under
linear polarization.
Figure 2 shows the time-average of ∆Im[ε(Tp , ω)] of each polarization calculated by the
analytic formula. The time-average corresponds to the probe by continuous-wave lasers, and
depends only on the density of states [18].
The blue line denoted the case of the circularly polarized laser. The exponential tail
below the band gap, which is a feature of the FKE, can be seen. The oscillation above the
band gap is due to the response of the different dressed states and the blue shift of the band
9
gap by Uc . In this case, Uc is approximately 0.95 eV, which corresponds to the energy at
which ∆Im[ε(ω)] has a sharp peak. Since the response of the l = −1 dressed state is intense,
the apparent blue shift of the band gap becomes 0.5 eV.
√
The red line denotes the case of the elliptic polarization (η = 1/ 2). The most obvious
difference is the magnitude of the blue shift. The original band gap is blue-shifted by
Ue = 0.71 eV. Although the overall behavior is similar to that of circular polarization, the
contribution of the l1 + l2 = −1 dressed band becomes weaker than that of the circular
polarization.
The case of linearly polarized light is represented by the green line. Here, Up is 0.5 eV
which is consistent with the apparent blue shift of the band gap. There is a small shoulder
at 5.1 eV. This energy corresponds to the l1 + l2 = −3 dressed state.
For reference, the numerical result for the static FKE [11] is represented by the black
dashed line. Since the applied field amplitude is stable with circular polarization, the tunneling effect is expected to be similar to the static FKE. As expected, the behavior under the
band gap agrees with the case of circularly polarized light. However, the oscillation above
the band gap is different. This result indicates that the FKE under circularly polarized light
is not equivalent to the static FKE, even though the amplitude of the electric field is static.
Thus, using the dressed states is indispensable.
V.
SUMMARY
We presented an analytical formulation for the time-resolved dynamical Franz-Keldysh
effect under an elliptically polarized laser. This formula connects the cases of linear polarization and circular polarization though the ellipticity of the pump laser. We found that
the time-dependent change in the optical properties observed under linearly polarized light
disappears as the ellipticity increases. Therefore, linearly polarized light is suitable for ultrafast control of the material response. On the other hand, the response of the dressed
state is still important for understanding the change in the optical properties caused by a
circularly polarized light field.
10
ACKNOWLEDGEMENT
This work is supported by a JSPS KAKENHI (Grants No. 15H03674). Numerical calculations were performed on the supercomputer SGI ICE X at the Japan Atomic Energy
Agency (JAEA).
[1] M. Hentschel, R. Klenberger, Ch. Spielmann, G.A. Reider, N. Milosevic, T. Brabec, U. Heinzmann, M. Drescher, and F. Krausz, Nature 414, 509 (2001).
[2] H. Hirori, A. Doi, F. Blanchard, and K. Tanaka, Appli. Phys. Lett. 98, 091106 (2011).
[3] A.H. Chin, O.G. Calderon, and J. Kono, Phys.Rev. Lett. 86 , 3292 (2001).
[4] H.Hirori, K. Shinokita, M. Shirai, S. Tani, Y. Kadoya, and K. Tanaka, Nat. Comm. 2, 594
(2011).
[5] A. Schiffrin, T. Paasch-Colberg, N. Karpowicz, V. Apalkov, D. Gerster, S. Mühlbrandt, M.
Korbman, J. Reichert, M. Schultze, S. Holzner, J.V. Barth, R. Kienberger, R. Ernstorfer, V.S.
Yakovlev, M.I. Stockman, and F. Krausz, Nature 493, 70 (2013).
[6] M. Schultze, E.M. Botschafter, A. Sommer, S. Holzner, W. Schweinberger, M. Fless, M.
Hofstetter, R. Kienberger, V. Apalkov, V.S. Yakovlev, M.I. Stockman, F. Krausz, Nature
493, 75 (2013).
[7] F. Nobelli, D. Fausti, F. Giusti, F. Parmigiani, and M. Hoffmann, Scientific Reports 3 1227
(2013).
[8] M. Schultze, K. Ramasesha, C.D. Pemmaraju, S.A. Sato, D. Whitmore, A. Gandman, J.S.
Prell, L.J. Borja, D. Prendergast, K. Yabana, D.M. Neumark, S.R. Leone, Science 346, 1348
(2014).
[9] W.Franz, Z. Naturforsch. Teil A 13, 484 (1958).
[10] L. V. Keldysh, Sov. Phys. JETP 34, 788 (1958).
[11] K. Tharmalingam, Phys. Rev. 130, 2204 (1963).
[12] B. O. Seraphin and R. B. Hess, Phys. Rev. Lett. 14, 138 (1965).
[13] R. E. Nahory and J. L. Shay, Phys. Rev. Lett. 21, 1569 (1968).
[14] H.Shen, M. Dutta, J. Appl. Phys. 78, 2151 (1995).
[15] J.K. Wahlstrand, J.E. Sipe, Phys. Rev. B82, 075206 (2010).
11
[16] F. Duque-Gomez, J.E. Sipe, J. Phys. Chem. Solids 76, 138 (2015).
[17] Y. Yacoby, Phys. Rev. 169, 610 (1968).
[18] A. P. Jauho and K. Johnsen, Phys.Rev. Lett. 76, 4576 (1996).
[19] K. B. Nordstrom, K. Johnsen, S. J. Allen, A. P. Jauho, B. Birnir, J. Kono, T. Noda, H.
Akiyama, and H. Sakaki, Phys. Rev. Lett. 81, 457 (1998).
[20] A. Srivastava, R. Srivastava, J. Wang, and J. Kono, Phys.Rev.Lett. 93, 157401 (2004).
[21] Y. Mizumoto, Y. Kayanuma, A. Srivastava, J. Kono, and A. H. Chin, Phys. Rev. B 74, 045216
(2006).
[22] S. Ghimire, A.D. DiChiara, E. Sistrunk, U.B. Szafruga, P. Agostini, L. F. DiMauro, and D.
A. Reis, Phys. Rev. Lett. 107, 167407 (2011).
[23] T. Otobe, Y. Shinohara, S. A. Sato, and K. Yabana, Phys. Rev. B 93, 045124 (2016).
[24] M. Lucchini, S. Sato, J. Herrmann, A. Ludwig, M. Volkov, L. Kasmi, Y. Shinohara, K. Yabana,
L. Gallmann, and U. Keller, Science 353, 916 (2016).
[25] K. Uchida, T. Otobe, T. Mochizuki, C. Kim, M. Yoshita, H. Akiyama, L. N. Pfeiffer, K. W.
West, K. Tanaka, and H. Hirori, JSAP-OSA Joint Symposia 2015 Abstracts, (Optical Society
of America, 2015), paper 15p 2E 10. Paper is submitted.
[26] W. V. Houston, Phys. Rev. 51, 184 (1940).
[27] H. R. Reiss and V. P. Krainov, J. Phys. A: Math. Gen. 36, 5575 (2003).
[28] E. Runge and E. K. U. Gross, Phys.Rev. Lett. 52, 997 (1984).
[29] G.F. Bertsch, J.-I. Iwata, A. Rubio, and K. Yabana, Phys. Rev. B 62 , 7998 (2000).
[30] T. Otobe, M. Yamagiwa, J. -I. Iwata, K. Yabana, T. Nakatsukasa, and G. F. Bertsch, Phys.
Rev. B77, 165104 (2008).
[31] N. Troullier and J.L. Martins, Phys. Rev. B 43, 1993 (1991).
[32] L. Kleinman and D. M. Bylander, Phys. Rev. Lett. 48, 1425 (1982).
[33] J. P. Perdew and A. Zunger, Phys. Rev. B 23, 5048 (1981).
[34] K. Yabana and G.F. Bertsch, Phys. Rev. B 54, 4484 (1996).
[35] K. Yabana, T. Sugiyama, Y. Shinohara, T. Otobe, and G. F. Bertsch, Phys. Rev. B 85, 045134
(2012).
12