Temperature-dependent dispersion model of float zone crystalline silicon

Applied Surface Science 421 (2017) 405–419

Contents lists available at ScienceDirect

Applied Surface Science journal homepage: www.elsevier.com/locate/apsusc

Temperature-dependent dispersion model of float zone crystalline silicon Daniel Franta a,∗ , Adam Dubroka b,c , Chennan Wang b,c , Angelo Giglia d , Jiˇrí Vohánka a , Pavel Franta a , Ivan Ohlídal a a

Department of Physical Electronic, Faculty of Science, Masaryk University, Kotláˇrská 2, 61137 Brno, Czech Republic Department of Condensed Matter Physics, Faculty of Science, Masaryk University, Kotláˇrská 2, 61137 Brno, Czech Republic c CEITEC – Central European Institute of Technology, Masaryk University, Brno, Czech Republic d IOM-CNR, s.s. 14, Km. 163.5 in AREA Science Park, 34149 Basovizza, Trieste, Italy b

a r t i c l e

i n f o

Article history: Received 29 July 2016 Received in revised form 2 February 2017 Accepted 5 February 2017 Available online 16 February 2017 Keywords: Crystalline silicon Optical constants Temperature dependence Ellipsometry Spectrophotometry Sum rule

a b s t r a c t In this paper, we present the temperature dependent dispersion model of float zone crystalline silicon. The theoretical background for valence electronic excitations is introduced in the theoretical part of this paper. This model is based on application of sum rules and parametrization of transition strength functions corresponding to the individual elementary phonon and electronic excitations. The parameters of the model are determined by fitting ellipsometric and spectrophotometric experimental data in the spectral range from far infrared (70 cm−1 ) to extreme ultraviolet (40 eV). The ellipsometric data were measured in the temperature range 5–700 K. The excitations of the valence electrons to the conduction band are divided into the indirect and direct electronic transitions. The indirect transitions are modeled by truncated Lorentzian terms, whereas the direct transitions are modeled using Gaussian broadened piecewise smooth functions representing 3D and 2D van Hove singularities modified by excitonic effects. Since the experimental data up to high energies (40 eV) are available, we are able to determine the value of the effective number of valence electrons. The Tauc–Lorentz dispersion model is used for modeling high energy electron excitations. Two slightly different values of the effective number of valence electrons are obtained for the Jellison–Modine (4.51) and Campi–Coriasso (4.37) parametrization. Our goal is to obtain the model of dielectric response of crystalline silicon which depends only on photon energy, temperature and small number of material parameters, e.g. the concentration of substituted carbon and interstitial oxygen. The model presented in this paper is accurate enough to replace tabulated values of c-Si optical constants used in the optical characterization of thin films diposited on silicon substrates. The spectral dependencies of the optical constants obtained in our work are compared to results obtained by other authors. © 2017 Elsevier B.V. All rights reserved.

1. Introduction Wafers of monocrystalline silicon (c-Si) are often used as a substrate for thin film growth in material science and microelectronics. The optical characterization of thin films is often performed in a wide temperature range from liquid helium temperature to high temperatures. In order to properly analyze the optical data of thin films, it is necessary to have a precise temperature-dependent optical constants of c-Si. Moreover, the optical constants of c-Si are dependent on dopants and impurities especially in the IR region

∗ Corresponding author. E-mail address: [email protected] (D. Franta). http://dx.doi.org/10.1016/j.apsusc.2017.02.021 0169-4332/© 2017 Elsevier B.V. All rights reserved.

[1–6]. The temperature dependence of the optical constants corresponding to various concentrations of dopants and impurities cannot be tabulated for practical purposes. It is necessary to point out that tabulation of the optical constants is impossible even when the wafers are taken from the same batch because the concentration of impurities and dopants varies along an ingot, i.e. for the individual wafers prepared from this ingot, but also in radial direction of this ingot. This means that individual pieces of silicon wafers will be slightly different. During the last four years we have dealt with the problem of parametrization of the optical constants of c-Si [7,8]. At present we have a temperature dependent dispersion model expressing the dependence of dielectric response on the concentration of oxygen, carbon and phosphorus. This model is implemented in newAD2 software [9].

406

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

In this paper, we present an advanced temperature-dependent dispersion model of float zone (FZ) c-Si. Wafers of FZ c-Si represent the purest c-Si available on the market at reasonable prices. Note that, the presented dispersion model also depends on the density of substituted carbon atoms. This is important because carbon contaminant is present even in FZ c-Si, and its concentration can be different for each sample. The model is developed on the basis of parametrizations of transition strength functions of individual elementary excitations normalized with respect to the sum-rule. The basic ideas were presented in our previous paper [10]. In that paper, only simple parametrizations describing individual elementary processes are discussed without details needed for construction of specific models. The current model is an improvement of the model of c-Si developed in [7] and [8]. In this paper, the theoretical background of the valence electron excitations will be presented in detail. The direct electronic transitions are modeled using Gaussian broadened piecewise smooth function representing 3D and 2D van Hove singularities (VHS) modified by many-body effects representing 3D Wannier and 2D Wannier (hyperbolic) excitons [11]. The indirect interband valence electron excitations are modeled using several truncated Lorentzian terms, the similar model was used for amorphous silicon [12,13]. Thermally induced free carrier contribution will be modeled by Drude-like term suggested in [10]. This contribution represents the indirect intraband transitions. The high energy electronic excitations will be modeled by two versions of Tauc–Lorentz model [14,15]. The theoretical background for phonon absorption will not be discussed in detail because it was already presented in our recent paper [8]. In that paper two-phonon absorption model suitable for describing absorption processes at the room temperature was presented. Above the room temperature the model anticipated lower absorption in the range 700–800 cm−1 . This discrepancy was explained by means of the absence of the higher order phonon absorption processes. This paper will show that the extension of the three-phonon absorption into this spectral range enables us to solve this problem. The model will be applied to temperature-dependent ellipsometric data (5–700 K) combined with ellipsometric and spectrophotometric data obtained at the room temperature (300 K). It will be shown that the presented model describes experimental data very well, thus it is accurate enough for practical purposes.

Ltd.) from the optical experiment so that no impact from vibrations was noticed in the spectra. The cryostat used KBr windows for mid-infrared range or silica windows for higher frequencies and was working under 8 × 10−8 mbar and 3 × 10−8 mbar pressure, respectively, at 300 K. We have tested that under these conditions, there were no signatures of freezing of residual atmosphere on the sample. The high temperature IR ellipsometric data were measured using Instec heating cell (HC) with ZnSe windows under reduced pressure, whereas UV-VIS data were measured under vacuum with UVISEL2 VUV Horiba Jobin Yvon ellipsometer with heating stage. The spectrophotometric data were measured at the room temperature with Bruker Vertex 80v in vacuum, Perkin Elmer Lambda 1050 in air and Elettra BEAR beamline in the UHV condition. Note that the ellipsometric quantities of sample #2 were also measured from the back side, i.e. from the side of aluminum film. The optical constants of alumina were determined from these measurements. The ellipsometric quantities Is , Ic and In used in this work correspond to the three independent elements of the normalized Mueller matrix M0 of isotropic systems:

2. Experiment

3. Data processing

Dispersion parameters of the advanced dispersion model were determined by processing experimental data obtained on two samples. The sample #1 was a piece of double side polished FZ c-Si wafer. The sample #2 was cut from the same wafer as the sample #1, but the back side was covered by aluminum film. Both samples were cut along (100) plane from 5 × 1013 cm−3 B-doped silicon with resistivity 80cm. The thickness of the samples was 0.58 mm. The instruments used for the experiment and the corresponding data sets are summarized in Table 1. It includes six table-top instruments covering the spectral range from the far IR to vacuum UV and the spectrophotometer at Elettra BEAR beamline [16,17] extending the spectral range of our measurements to the extreme UV. The ellipsometric data measured at low temperatures were obtained with the help of a closed-cycle helium cryostat mounted on VASE and IR-VASE Woollam ellipsometers. The cryostat was equipped with an ultra-low vibration interface (ColdEdge Technologies) that allowed to decouple the mechanical vibrations from the Gifford-MacMahon refrigerator (Sumitomo Heavy Industries,

All experimental data mentioned in Table 1 were fitted simultaneously in the newAD2 software [9]. This software can fit together an arbitrary set of optical data using a single consistent model. All parameters of the model are fitted simultaneously and the weights of individual data sets are set such that the residual sum of squares is equally distributed between these data sets [19–21]. This equalization permits us to combine individual data sets with different spectral and angular densities from different spectral ranges obtained using different instruments exhibiting different systematic errors while keeping the sensitivity of the method to effects manifested in the individual data sets. The influence of the native overlayer was taken into account as thin homogeneous film. The dielectric response of the overlayer was modeled by the universal dispersion model [21]. The thickness of the overlayer was determined separately for each measurement and each temperature. The dielectric function of the overlayer was assumed to be dependent on temperature, but it was assumed to be the same for different measurements. This is necessary because the native overlayer is an unstable system that depends on oxida-

⎛

−In

0

1

0

⎝ 0

0

Ic

0

0

−Is

1

⎜ ⎜ −In

M0 = ⎜ ⎜

0

⎞

⎟ ⎟. ⎟ Is ⎠ 0⎟

(1)

Ic

These three ellipsometric quantities are suitable for the description of depolarization effects, which must be taken into account, especially if the back side reflection is enhanced by the aluminum film. The degree of polarization P is related to these quantities as follows: P=



Is2 + Ic2 + In2 .

(2)

Note that without the reflection from the back side of the sample the ellipsometric data would not carry enough information about slight absorption occurring in the slab. The depolarization caused by the back side reflections is crucial for the determination of temperature dependencies of phonon absorption and temperature dependent indirect absorption edge in c-Si. The multiple back side reflection were calculated using Mueller matrix formalism. The fact that the beam is slightly shifted by each back side reflection and thus part of it falls outside of the aperture of the detector was also taken into account.

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

407

Table 1 List of experimental and tabulated data used for construction of the advanced dispersion model. The superscripts f and b denote the front and back side measurements. Instrument/source

Sample

Spectral range

Measurement

Temperature (K)

Points

Figure

IR-VASE @ cryo. IR-VASE @ HC IR-VASE Vertex 80v Vertex 80v VASE @ cryo. UVISEL2 VUV UVISEL UVISEL Lambda 1050 Elettra BEAR CHARMS [18] CHARMS [18]

#2 #2 #2 #1 #2 #2 #2 #1 #2 #1 #2 Prism Prism

400–6200 cm−1 500–6200 cm−1 263–7800 cm−1 70–7500 cm−1 70–5000 cm−1 0.6–6.5 eV 0.6–8.3 eV 0.6–6.3 eV 0.6–6.5 eV 0.69–1.38 eV 6–40 eV 1.1–5.5 ␮m 2 ␮m

Isf , Icf , Inf @ 75◦ Isf , Icf , Inf @ 70◦ Isb , Icb , Inb @ 55-75◦ T @ 0◦ R @ 10◦ Isf , Icf , Inf @ 75◦ Isf , Icf , Inf @ 70◦ Isf , Icf , Inf @ 55–75◦ Isb , Icb , Inb @ 55–75◦ T @ 0◦ R @ 6◦ n n

5, 100, 200, 300, 320, 380 300, 400, 500, 600 300 300 300 10, 100, 200, 300, 330, 390 300, 325, 350, 400, 450, 500, 550, 600, 650, 700 300 300 300 300 30–295 5–700 extrapolated by (3)

18042 2716 1405 4007 1512 3218 1550 2855 1185 181 226 156 696

3,4 3, 4 – 5 5 6, 11 6 8, 10 – 7 9 12 13

5–700

37749

Total

0.0086–40 eV

tion, adsorption and desorption processes. Thus, the overlayer is influenced by temperature, pressure and history of the sample in general. The thickness of the aluminum film on the back side of the sample #2 was about 100 nm. Therefore this film was nontransparent in the whole measured spectral range. The dielectric constants were also modeled by the universal dispersion model with Drude-like term representing the free carrier contribution. Results concerning this film, will not be presented in this work. The spectrophotometric and ellipsometric measurements enable us to determine the optical constants with relative precision that is not better than approximately 10−3 . In order to improve the precision of optical constants the temperature and spectral dependencies of refractive index published by Frey et al. [18] in transparent region (1.1–5.5 ␮m) were fitted simultaneously with our experimental data. The data published by Frey were obtained by the minimal deviation method on silicon prisms with special instrument for twelve temperatures in the range 30–295 K. The relative precision of these data is in the range 10−5 . In order to suppress the oscillations in temperature dependency of optical constants outside of the twelve temperatures where the Frey’s data are available it was necessary to also consider data in temperature range 5–700 K calculated at wavelength 2 ␮m using the formula n(T ) = n0 +

1 2 + , exp(1 /T ) − 1 exp(2 /T ) − 1

(3)

where n0 = 3.42477, 1 = 0.01992, 1 = 210.72 K, 2 = 0.0737318, 2 = 636.501 K. The above formula assumes that the refractive index in the near infrared (NIR) region is inversely proportional to the position of absorption structures and, as we will show in Section 5.5, the position of these structures is governed by average Bose–Einstein statistical factor [22,23]. In order to reach sufficient accuracy two such factors were necessary. The constants in the formula were determined by fitting to Frey’s data.

transition density (joint density of states) function J based on correlation of density of initial and final states is presented. Do not confuse this function with the function J defined from the dielectric response (4). For the function F the sum rule takes the following simple form:



∞

F(E) dE = N,

where N is the total transition strength related to the density of atoms Na and mean number of electrons per atom Z as 2

N=

(eh) ZNa U = ZNa U, 8␲0 me

εˆ (E) = 1 +



Nt εˆ 0t (E),

where the unity can be understood as the vacuum contribution. The normalization of εˆ 0t (E) is performed with respect to sum rule, i.e.



∞

Eε0i,t (E) dE = 0

(4)

In [10] the quantum mechanical quantity called transition strength function F is introduced as integral form of the oscillator strength f in the Thomas–Reiche–Kuhn sum rule. In the frame of dipole approximation F is proportional to the real part of the complex conductivity  r or transition strength function F defined from the dielectric response (4). In [10] the proper definition of the

(7)

t



For the parametrization of the dielectric function εˆ it is helpful to define the transition strength function F and transition density (joint density of states) function J [24] as

(6)

where h, e, 0 and me are the universal physical constants, i.e. the Planck’s constant, electron charge, vacuum permittivity and electron mass, respectively. Number of electrons per atom Z and quantity U are also constants but specific for individual materials (Z = 14, U = 1.0002735 for silicon, for details see [10]). The quantity U is the correction factor taking into account the contribution of nuclei to the polarizibility of the material. The quantity Na is called the density parameter and it has the same unit as the total transition strength N, i.e. squared energy (eV2 ). In the formalism introduced in [10], the elementary excitations denoted t are described by the normalized dielectric functions εˆ 0t (E) and the parameters called transition strengths Nt (sometimes referred only as strengths). The dielectric function is then calculated as a sum of individual contributions

4. Parametrization of dielectric function

F (E) = Eεi (E) and J (E) = E 2 εi (E).

(5)

0

∞

Ft0 (E) dE = 1,

(8)

0

where Ft0 (E) are called the normalized transition strengths functions. The real and imaginary part of the normalized dielectric function of individual contributions can be calculated from Ft0 (E) as follows



ε0r,t (E)

∞

Ft0 (X) 2 = − dX,  0 X 2 − E2

ε0i,t (E) =

Ft0 (|E|) , E

(9)

where the first equation is the Kramers–Kronig relation [25,26] and the second follows from the definition in (4). The above definitions ensure the time-reversal symmetry [10].

408

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

The transition strength Nt can be expressed using dimensionless effective number of electrons per atom nt and density parameter Na as Nt = nt Na .

(10)

5. Excitation of valence electrons The effective number of valence electrons nve is assumed to be independent on temperature and it is a fitting parameter in our model. This means that the temperature dependence of the transition strength of valence electrons Nve (T) is only through thermal expansion [27], which is expressed by temperature dependence of density parameter Na (T) introduced in our previous paper [7]: Nve (T ) = nve Na (T ).

(11)

In our model the excitations of valence electrons to the conduction band ‘vc’ and higher energy states ‘vx’ are distinguished: Nvx (T ) = ˛x Nve (T ),

Nvc (T ) = (1 − ˛x )Nve (T ),

(12)

where ˛x is a temperature independent dimensionless fitting parameter. The corresponding ‘vc’ transition strength is divided between indirect and direct electronic transitions. The temperature dependence of strength of indirect transitions Nidt (T) is given by a phonon occupation number, i.e. by Bose–Einstein statistics [11,10]. The exact expression will be given by (19) in Section 5.1. The strength of direct transitions is then given as the complement to the total transition strength of valence electrons: Ndt (T ) = Nvc (T ) − Nidt (T ).

(13)

This represents the fundamental temperature dependence of valence electron excitations from the point of view of sum rules [28,10]. Of course, the temperature dependencies of the spectral distribution of individual contributions are not included into this fundamental dependence and they will be defined in Section 5.5.

processes. Thus the summation over creation branch, denoted by +p, and annihilation branch, denoted by −p, is performed as follows: Fp0 (E, T ) =

0 (E) + f BE (E˙ , T ) F 0 (E) [f BE (E˙ p , T ) + 1]F+p p −p

2f BE (E˙ p , 300 K) + 1

where fBE (E, T) = 1/(exp[E/(kB T)] − 1) is the Bose–Einstein statis0 (E) are tics and kB is the Boltzmann constant. The functions F±p normalized transition strength functions describing distribution of creation or annihilation branches at 300 K (room temperature) [11,10]. Apart from the temperature dependence given by the Bose–Einstein statistics, there is also weaker temperature dependence given by shifts in electron and phonon structures and by changes in values of matrix elements. For this reason some dispersion parameters in our model has to be considered weakly temperature dependent. The parameter Ep in the above equation has such weak temperature dependence, which causes small changes in the normalized transition strength functions and almost negligible change in the value of f BE (E˙ p , T ). In the whole article the temperature dependent parameters will be dotted if this temperature dependence influences only the distribution of transition strength. The temperature dependencies influencing the integral values of transition strengths will be denoted by symbol T in brackets. The parametrization of temperature dependence of dotted quantities will be discussed in Section 5.5. 0 (E) are modThe normalized transition strength functions F±p eled by a sum of the analytical truncated Lorentzian terms, indexed by j. The truncated Lorentzian terms, which are nonzero only between the bandgap energy and the high energy limit of indirect interband transitions, are usually used to describe response of amorphous or disordered solids [12,13,21]. Therefore, the normalized transition strength function of one term must be indexed by j and ±p: F 0j,±p (E)

=

2 2 (E˙ g ± E˙ p − |E|) (E˙ h ± E˙ p − |E|) E˙ g ±E˙ p ,E˙

εi,idt (E, T ) =

Nvc (T )

[ ˛p F 0p (E, T ) + ˛fc F 0fc (E, T )], E

(14)

p

where the index p represents TO (transversal optical), L (longitudinal), and TA (transversal acoustical) branches of interband indirect transitions. The intraband indirect transitions are denoted ‘fc’ (free carrier contribution). The dimensionless parameters ˛p and ˛fc are relative transition strengths related to the total strength of electron transitions between valence and conduction bands Nvc (T) at 300 K. In principle the summation over p should include summation over all six phonon branches and integration over initial and final electron states in the Brillouin zone. From the practical point of view it is not necessary to do this integration and perform the summation over all phonon branches. In our case it was found that it is enough to sum over three phonon branches as indicated in (14). The phonon energies Ep of these three branches correspond to those in X direction that are available as common fitting parameters employed in phonon part of the presented dispersion model (parametrized using the frequencies fTO (X), fL (X) and fTA (X)). The X point of symmetry was chosen because minimum energy of interband indirect electronic transitions is realized between absolute maximum of valence bands in the  point and absolute minimum of conduction band located at the line close to the X point [11]. The temperature dependence of indirect interband transitions is given mostly by changes of transition probability that are proportional to the occupation numbers of phonons involved in absorption

˙

h ±Ep

2

2] CN,j,±p |E| [(E˙ cj ± E˙ p − |E|) + B˙ cj

5.1. Indirect electronic transitions The complete model of indirect electronic transitions must include both interband and intraband transitions:

(15)

,

(|E|) ,

(16)

where the expression for the complex dielectric function and normalization constant CN,j,±p are given in Appendix A in [12]. The function a,b (x) is unity for x ∈ [a, b] and zero outside of this interval. The truncated Lorentzian terms are described by parameters Ecj (peak position) and Bcj (peak width) and strength parameter Aidt,j . In the presented model five terms j = 1, . . ., 5 were chosen. There is no substantial physical reason for choosing this number of terms, but it ensures sufficiently free distribution function describing interband indirect transitions with correct course above the bandgap energy Eg , i.e. the quadratic Tauc form. This function also has quadratic course below the maximum transition energy Eh . The summation over the j-index is performed as follows: 0 F±p (E) =

5

Aidt,j

5

A j=1 idt,j

j=1

0 Fj,±p (E)

(17)

which ensures the normalization of the resulting functions. One of the Aidt,j parameters must be fixed in an arbitrary positive value, to avoid correlation between the parameters. Temperature dependence of transition strength of intraband indirect transitions is given as [7,29] 0 Ffc (E, T ) =

T 300 K

3/2



exp −

E˙ g 2kB

1 T

−

1 300 K



0 Ffc (E) ,

(18)

where F 0fc (E) is a normalized Drude-like transition strength function at 300 K. The temperature dependent factor in front of F 0fc (E), which will be denoted by ffc (T), is derived using the parabolic band approximation from Fermi–Dirac statistics of electrons and holes. This approximation is sufficient for temperatures studied in this

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

409

Table 2 List of critical point of energies used as parameters in the presented dispersion model. Note that all critical point energies in the table are identified by the corresponding points of symmetry. In the standard notation the identification by points of highest symmetry , X, L is omitted, e.g. E0 , E1 , and points with lower symmetry are placed in brackets, e.g. E2 () [23,11]. Energy

@ 0 K (eV)

@ 300 K (eV)

 (K)

Description

1.087

248.6 –

Bandgap Maximum energy

Direct transitions (critical points)  3.398 3.367 E0 E1L 3.457 3.393 4.068 4.149 E0 E1K 4.239 4.195 X E2 4.333 4.302  E2 4.517 4.427 L E1 5.442 5.398 7.159 7.119 E3 X E2 10.72 10.51 E3X 12.18  E3X 19.33

– 115.7 – – 349.2 – 403.6 – – – –

T dep. by E1L

High energy excitations Egx 8.748

–

Indirect transitions 1.14 Eg 24.6 (fixed) Eh



T dep. by E1L T dep. by E2X

Fig. 1. Schematic diagram of unbroadened transition strength function corresponding to 3D van Hove singularities representing one integration path. The Arabic numbers 0–3 in areas correspond to individual J contributions in (22)–(25), (28) and (30). The Roman numbers I–III define intervals between singularities, see (26) and (27).

T dep. by E2X T dep. by E2X T dep. by E1L

paper. The temperature dependence of the strength of indirect transitions is then given by slight dependence caused by thermal expansion and temperature factor corresponding to Bose–Einstein statistics of interband transitions:



Nidt (T ) = Nvc (T )





˛p

p

2f BE (E˙ p , T ) + 1 + ˛fc ffc (T ) 2f BE (E˙ p , 300 K) + 1

singularities approximating anisotropic critical points [30]. Thus, the transition strength function can be divided into the 3D and 2D parts. Within the parabolic band approximation the 3D part can be modeled by the sum of four transition density (joint density of states) functions as follows: J0,k (E) = A˙ 0,k



.

(19)

XI,k (E) + YII,k (E) ,

J1,k (E) = A˙ 1,k 1 −



YI,k (E) + YII,k (E) ,



Although the term ˛fc ffc (T) corresponding to intraband electronic transitions has negligible contribution to the value of Nidt (T) in (19), it cannot be neglected in (14), since it is important for the dielectric response in the IR region. The intraband electronic transitions will not be further discussed in this paper, because this was already done in [7], where the value of ˛fc = nfc /nvc = 8 ×10−13 was found.

J2,k (E) = A˙ 2,k XII,k (E) + 1 − and

(22)



J3,k (E) = A˙ 3,k XII,k (E) +



(23)



XIII,k (E)

(24)

YIII,k (E) ,

(25)

5.2. Direct electronic transitions

where the functions X and Y are nonzero only in intervals I, II, III (see Fig. 1). For the interval I, i.e. for E ∈ (E˙ 0,k , E˙ 1,k ), they are defined as:

A contribution of the direct electronic transitions is modeled 0 (E) as follows: using normalized transition strength function Fdt

XI,k (E) =

εi,dt (E, T ) = Ndt (T )

0 (E) Fdt

E

,

(20)

where Ndt (T) was defined in (13). The normalized transition strength function is parametrized using broadened piecewise smooth function. The discontinuities of this function are in 11 critical points introduced in Table 2. The positions of the critical points depend on temperature which is usually investigated temperature dependence and it will be discussed in Section 5.5. The presented model replaces 3D integration over the Brillouin zone by several contributions between sets of four M0 , M1 , M2 and M3 critical points called integration paths (see Table 3). In our model, five integration paths are sufficient for describing complete transition strength function: 0 Fdt (E)

=

5

k=1

Adt,k

5

A k=1 dt,k

0 Fdt,k (E),

(21)

where Adt,k are fitting parameters controlling the relative strengths of individual integration paths and, therefore, one parameter must be fixed in arbitrary value, otherwise the parameters would be 0 (E) among the four criticorrelated. The course of functions Fdt,k cal points is based on combination of isotropic 3D and 2D van Hove

YI,k (E) =

⎧ ⎨ E − E˙ 0,k

for E˙ 0,k < E < E˙ 1,k

E˙ 1,k − E˙ 0,k

⎩

0

⎧ ⎨ E˙ 0,k − E ⎩

(26)

otherwise, for E˙ 0,k < E < E˙ 1,k

E˙ 1,k − E˙ 0,k 0

(27)

otherwise.

The functions XII,k (E), YII,k (E), XIII,k (E) and YIII,k (E) are defined in the same way with energies E˙ 0,k and E˙ 1,k in the above definitions replaced by energies at boundaries of the corresponding interval (see Fig. 1). The sum of functions in Eqs. (22)–(25) represents transition density without excitonic effect. If the electron–hole interactions represented by Rydberg energies R˙ k are taken into account, the square root functions in Eqs. (22) and (25) must be modified as



J0,k (E) = A˙ 0,k

∞

2R˙

k n3

RI,k (E˙ 1,k )

ı

|E| − E˙ 0,k

R˙ k + 2 n



+

RI,k (E˙ 1,k ) RI,k (E)

+ YII,k (E)

,(28)

n=1

where RI,k (E) = 1 − exp

 −2

 R˙ k E − E˙ 0,k

 (29)

410

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

Table 3 Schema of integration paths and related parameters. Broadening parameters B and Rydberg constants R are in meV units. Path k 1 2 3 4 5

E0,k 

E0  E0  E0  E0 E0

and



J3,k (E) = A˙ 3,k where

E1,k

E2,k

E3,k

Adt,k

0K B1,k

100K B1,k

300K B1,k

E1L E1L E1K E1K  E1L

E2X  E2X E2X E2  E2X

E3 E3X E3X E3X  E3X

1 (fixed) 0.683 0.765 0.914 0.739

29.2 70.6 216 97.3 89.7

50.7 295 277 149 127

57.3 360 265 155 145

RIII,k (E˙ 2,k ) XII,k (E) + RIII,k (E)





RIII,k (E) = exp

57.8 170 501 134 83

0K B2,k

100K B2,k

300K B2,k

700K B2,k

Rk0K

Rk300K

26 456 567 106 201

41 456 492 112 201

59.5 456 424 120 201

98.4 456 700 190 201

7.93 120 3.68 0.0113 1.07

6.12 132 2.56 0.368 1.73

 (30)

,



R˙ k

2

700K B1,k

− 1.

E˙ 3,k − E

(31)

Note that in the limit for small R the Eqs. (28) and (30) are identical with Eqs. (22) and (25). Similarly as in the 3D part the 2D part can be modeled by a sum of eight transition density functions given as follows: J4,k (E) = A˙ 4,k YII,k (E),

(32)

J5,k (E) = A˙ 5,k XII,k (E),

(33)

J6,k (E) = −A˙ 6,k ln YI,k (E),

(34)

J7,k (E) = −A˙ 7,k ln XII,k (E),

(35)

J8,k (E) = −A˙ 8,k ln XIII,k (E) ,

(36)

J9,k (E) = −A˙ 9,k ln YII,k (E),

(37)

J10,k (E) = A˙ 10,k

(38) +YI,k (E)

and J11,k (E) = A˙ 11,k XIII,k (E),

 ⎧  ˙ ⎨ ln E1,k − E for E˙ 0,k < E < E˙ 1,k E˙ 1,k − E˙ 0,k ln YI,k (E) = ⎩ 0

 32R˙ k



∞

(2n − 1)

3

ı

|E| − E˙ 1,k +

R˙ k E − E˙ 1,k

RIIL,k (E) = 1 + exp

−2

J5,k (E) = A˙ 5,k XII,k (E)

RIIR,k (E˙ 1,k ) , RIIR,k (E)

RIIR,k (E) = 1 + exp

4R˙ k (2n − 1)

2

(41)





2

|E| − E˙ 0,k +

4R˙ k

(2n − 1)2 (45)



R˙ k E˙ 2,k − E

 ,

(42)

(43)

 ,



 RI,k (E) = 1 + exp

−2

J11,k (E) = A˙ 11,k XIII,k (E)



RIIL,k (E˙ 2,k ) , RIIL,k (E)





 RI,k (E˙ 1,k ) , RI,k (E)

(40)

otherwise.

In the case of the 2D part the excitonic effects must be applied to contributions 4, 5, 10 and 11:

n=1

(2n − 1)3

␦

(39)

where ln must be understood as an operator in the corresponding interval, i.e.

J4,k (E) = A˙ 4,k

∞ 

32R˙ k n=1

J10,k (E) = A˙ 10,k YI,k (E)

+ YII,k (E)

Fig. 2. Schematic diagram of unbroadened transition strength function corresponding to 2D van Hove singularities representing one integration path. The Arabic numbers 4–11 in areas correspond to individual J contributions in (32)–(39), (41), (43), (45) and (47). The Roman numbers I–III define intervals between singularities, see (26) and (27).

(44)

 ,

RIII,k (E˙ 2,k ) , RIII,k (E)

 RIII,k (E) = 1 + exp

R˙ k E − E˙ 0,k

2

 R˙ k E˙ 3,k − E

(46)

(47)

 .

(48)

Note that in both 3D and 2D parts the excitonic effects cause redistribution of transition strength from higher energies to lower energies (see Figs. 1 and 2 ). Moreover, excitonic effects cause existence of discrete states below minimum energy of transitions, which are represented by delta functions in Eqs. (28), (41) and (45). The values of Rydberg constants R˙ k control magnitudes of these excitonic effects. In the presented model there is just one for each integration path and their temperature dependencies express the change of lifetime of electron–hole pairs. The quantities A˙ 0,k –A˙ 11,k control relative strength of the individual contributions and their temperature dependencies express changes in the matrix elements (probability of transitions) in consequence of thermal changes of the band structure. It should be emphasized that neither the parameters R˙ k nor A˙ 0,k –A˙ 11,k influence to the integral value of transition strength of direct transitions which is given by Eq. (13). The temperature dependencies will be discussed in Section 5.5.

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

The normalization constants CN,k are calculated by means of integrating overall unbroadened transition strength function as follows:



Ji,k (E)

∞ 11

CN,k = 0

i=0

E

(49)

dE.

For the calculation of the dielectric response of the direct transitions it is necessary to carry out ε-broadening conserving the transition strength:



Ji,k (E)

ˆ∗ ␤

␧ˆ i,k (E) =

 ,

E2

(50)

ˆ is the complex where the symbol * denotes convolution and ˇ Dawson–Gaussian function. The broadening procedure is described in detail in our previous paper [24]. The broadening parameter depends on temperature and also on energy. However, it was sufficient to divide broadening in two parts: contributions 0, 1, 4, 6, 7, 10 were broadened by parameter B˙ 1,k and 2, 3, 5, 8, 9, 11 by B˙ 2,k . The temperature dependencies will be discussed in Section 5.5. Taking into account (21) the complex dielectric function corresponding to the direct transitions is then given as: 5

εˆ dt (E) = Ndt (T )

k=1

Adt,k

CN,k

5

11

A l=1 dt,l i=0

εˆ i,k (E).

(51)

5.4. Core electrons excitations The core electron excitations were modeled using the same simple model as we used for amorphous hydrogenated silicon in [12]. 5.5. Temperature dependencies In the literature the temperature dependencies of the critical points are modeled by average Bose–Einstein statistical factor with average phonon energy [22,23]. For example, the critical point E1L is parametrized as:



E1L (T ) = E1L,0K + E1L,300K − E1L,0K



T 300



(55) 



Bx (|E| − Egx )2 Egx ,∞ (|E|) 2

[(|E| − Egx )2 − (Ecx − Egx )2 ] + Bx2 (|E| − Egx )2

,

(52)

where Ecx and Bx are parameters describing position and width of the truncated Lorentzian peak. The imaginary part of the dielectric function εi,vx (E, T) is calculated as follows 0 (E) Fvx E

(54)

,



where L0K and L300K are given by parameters E0,0K and E0,300K in the following logarithmic substitutions:

Transition strength of the excitations of valence electrons into the states above the conduction band is modeled by a four parameter function. The shape of this distribution function corresponds to the classical damped harmonic oscillator model truncated in minimum energy of transitions denoted Egx . In the vicinity of Egx the function has quadratic Tauc course F ∼ (E − Egx )2 and for high energies the decay is classical F ∼ 1/E2 , i.e. corresponding to the Lorentz model. This function can be parametrized in two ways. The first of them was suggested by Campi and Coriasso [14] and the second one called Tauc–Lorentz (TL) was suggested by Jellison and Modine [15]. Both the parametrization give similar results so that they can be called two versions of the Tauc–Lorentz model. The difference between these models is only for higher energies. For these models the decrease of transition strength in the tail is classical but tail is higher for Jellison–Modine (JM) version than for Campi–Coriasso (CC) version for the same values in the peak region at the lower energies. This difference implies the advantage of the CC version in comparison with JM version of TL model in practice if the parameter ε(∞) is not used. If the TL model is used in limited spectral range of low energies together with the parameter ε(∞), both version are practically equivalent. Within the CC version of TL model the normalized transition strength function is written as:

εi,vx (E, T ) = Nvx (T )

exp(1 /T ) − 1

E0 (T ) = E1L (T ) + exp L0K + (L300K − L0K )

0K

2 

 exp(1 /300) − 1

where the parameters E1L,0K and E1L,300K are critical point energies at 0 and 300 K. The parameter 1 represents average phonon energy for this critical point in Kelvin units. From the practical point of view this three-parameter equation can not be used to parametrize the critical points that are close to each other. For example, the  inequality E0 < E1L must be ensured for all temperature values for  which the model is calculated. This problem E0 < E1L was solved so  L that E0 (T ) was derived from E1 (T ) by means of the following two parameter formula:

5.3. High energy excitations

0 Fvx (E) =

411

(53)

and the real part must be expressed by the Kramers–Kronig integral.

L

= ln

E0,0K E1L,0K





and

300K

L

= ln

E0,300K E1L,300K



.

(56)

The notes “T dep. by E1L ”, etc. in Table 2 distinguish this type of temperature dependency from the dependency given by (54).  The critical point energies E3X and E3X were assumed to be independent on temperature because no temperature dependent data were available. The bandgap energy Eg was parametrized by average Bose–Einstein statistical factor. The same trick with logarithmic substitutions was utilized for temperature dependencies of parameters A, B and R introduced in Sections 5.1 and 5.2. These parameters were parametrized at several selected temperatures and the dependencies among the logarithmic substitutions were interpolated by cubic splines. The different parameters were modeled in different number of points. The parameters R were modeled at 0 and 300 K, while parameters B were modeled at 0, 100, 300 and 700 K (see Table 3). In order to accurately express the temperature dependence of the change of the matrix elements the parameters A had to be modeled in eight points, i.e. at 0, 100, 200, . . ., 700 K. 6. Results and discussion 6.1. Phonon absorption The model for infrared region was taken from [8], where the model of dielectric response based on room temperature data was presented. In that paper, the two- and three-phonon absorptions could not be distinguished, since only data at one temperature were processed. In the older model which was used in that paper, the three-phonon absorptions below 2 O () ≈ 1000 cm−1 , which is the theoretical upper limit for the two-phonon absorptions, were neglected. It was shown that this assumption does not lead to the correct temperature dependency of the dielectric response. In the current paper, the dispersion model is extended to the temperature range from 0 to 600 K (see Figs. 3–5), thus it is possible to separate one-, two- and three-phonon absorptions in the range below

412

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

Fig. 5. Spectral dependencies of near-normal reflectance R and normal transmittance T in IR region and corresponding separation of transition strength to one-, twoand three-phonon contributions. Note that reflectance was measured for sample #2 whereas transmittance was measured for sample #1.

Fig. 3. Spectral dependencies of ellipsometric quantities Is , Ic , In and P in IR for chosen temperatures. Some ellipsometric quantities are shifted by offsets labeled in graphs.

Table 4 Phonon frequencies in points of symmetry at 0 and 300 K. The values are in THz units and dagger means fixed values. The degenerated phonon frequencies are not repeated in the table. The values of phonon frequencies in the empty fields are identical with the values on the left. fTA1

Fig. 4. Spectral and temperature dependencies of ellipsometric quantity In in IR region and corresponding transition strengths. Every ellipsometric quantity measured above 5 K is shifted by offset +0.05 relative to the previous curve.

T=0K  X L W K

0† 4.74 3.13 6.37 5.31

T = 300 K  X L W K

0† 4.72 3.11 6.31 4.82

fTA2

fLA

6.66

12.13 11.08 10.70 10.77

6.74

12.00 11.04 10.65 10.84

fLO

fTO1

fTO2

14.17 15.05 14.38 14.31

14.18

14.08 14.98 14.28 13.89

14.87

15.78 13.02 11.53 15.59 12.94 11.34

2 O (). In the presented model the three-phonon absorption was extended bellow 1000 cm−1 and an absorption peak representing one-phonon absorption on substituted carbon atoms at wavenumbers C0K = 612.7cm−1 and C300K = 610.1cm−1 was added. From Figs. 3–5 it is evident that the presented model describes the phonon response perfectly for spectrophotometric purposes at the room temperature and sufficiently outside the room temperature. The parameter values of the phonon absorption are not presented here with the exception of phonon frequencies (see Table 4). The temperature dependencies of the phonon frequencies (energies) was modeled using the three parameter formula (54) with the common average phonon energy ph = 426.1 K.

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

413

Fig. 6. Spectral and temperature dependencies of ellipsometric quantities Ic and P in VUV–NIR region and corresponding transition strengths. Ellipsometric quantities measured above 10 K are shifted by offsets in brackets.

Fig. 7. Spectral dependence of transmittance T in NIR region and corresponding transition strengths. The effective phonon energy ±Ep corresponds to phonon frequency of L branch in X point of symmetry at 300 K, i.e. fL300K (X) = 12 THz. Transmittance calculated using Herzinger optical constants [31] is added for comparison.

6.2. Electron excitations The experimental data and their fits corresponding to the excitations of the valence electrons to the conduction band and above are introduced in Figs. 6–9 . The temperature dependence of the bandgap energy Eg is seen as shift of the fundamental absorption edge located around 1 eV in Fig. 6. The corresponding parameters are summarized in Table 2.

Fig. 8. Spectral and angular dependencies of ellipsometric quantities Is , Ic , In and P in UV–NIR region at 300 K and corresponding separation of transition strength to individual contributions forming the dispersion model. The vertical bars represent positions of critical point energies labeled in graph and minimum energy of indirect transitions Eg (bandgap). The denotation of critical point energies by the points of highest symmetry , X, L is omitted in this figure.

414

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

Fig. 9. Room temperature reflectance in EUV region and corresponding transition strengths. The vertical bars represent positions of critical point energies labeled in graph, maximum energy of indirect transitions Eh and minimum energy of high energy excitations Ex . The denotation of critical point energies by the points of highest symmetry X is omitted in this figure.

In this table the value of maximum energy of indirect transition Eh is also introduced. This energy cannot be determined from the reflectance spectra (see Fig. 9), since it is not associated with any visible structure at these energies. Therefore, the fit of the experimental data was insensitive to this parameter and it was fixed in value of 24.6 eV estimated from ab initio calculation which will be presented in the forthcoming paper. In the present dispersion model the summation over three effective phonons corresponding to the three phonon branches in X point was performed. The relative transition strengths at 300 K are: ˛TA = 7 ×10−5 , ˛L = 0.0014 and ˛TO = 0.00072. The phonon frequencies fL0K (X) = 12.13THz and fL300K (X) = 12THz of the strongest L branch are in a good agreement with effective phonon frequency fp ≈ 11.6THz (Ep ≈ 48meV) determined in our earlier work [7] or with frequency fp ≈ 12.5THz (Tp ≈ 600 K) determined in classical Macfarlane and Roberts work [32,11]. The course of transmittance and transition strength function in the vicinity of the bandgap is plotted in Fig. 7. The vertical bars at Eg + Ep and Eg − Ep indicate the absorption edges corresponding to the minimal energy of excitation of electrons from the valence to the conduction band with simultaneous creation and annihilation of phonons. The total relative transition strength of indirect transitions is ˛idt = ˛TO + ˛L + ˛TA =0.00214, which is much less than the value ˛idt = nidt /nvc = 0.13 presented in our earlier paper [7]. In that paper optical characterization was based on ellipsometric data in the temperature range 300–600 K and simpler model, i.e. with smaller number of parameters. The discrepancy between the current value and the value obtained in this older work can be explained by different separation of the electron transition strength into parts corresponding to direct and indirect transitions. We believe that some of indirect transitions are incorrectly included into direct transitions in the current model. This is supported by the fact that the current model has maximum of the transition strength for indirect electron excitations around 2.5 eV (see Fig. 8) whereas the older model has maximum around 3.5 eV, i.e. around E1L critical point. The maximum of indirect transitions near E1L exists also in the dis-

ordered structure of amorphous silicon. The lower portion of path 2 in the current model of direct transitions probably corresponds to indirect transitions (see Fig. 8). Although the separation that was used in the older work may be better from the physical point of view, more accurate description of the temperature dependent dielectric function is achieved with the current model. As was already described in the theoretical part, five integration paths in the first Brillouin zone were used for calculating contributions from direct electronic transitions. The separation of the transition strengths into parts corresponding to different integration paths is in Fig. 8, the corresponding parameters are in Table 3. The energies of critical points and their temperature dependencies are summarized in Table 2. The identification of critical points is based on ab initio calculations. At highest points of symmetry (, X, L) it coincides with the standard identification [23,11]. The critical point energies E0 , E1K and points above 6 eV are not usually mentioned in the literature or they are not mentioned at all. We included them on the basis of ab initio calculations. Note that the identification of these critical points forming weak structures in the spectra may not be completely correct and could be revised in the future. For example, our model uses three critical points for the structure around 4.3 eV, similar model was used by Kondo and Moritani [33]. It is not clear whether E1K point is real one or just artifact. These ab initio calculations and a comparison with the current model will be presented in the forthcoming paper. The contributions from path 1, path 4 and path 5, which have small broadening factors, describe sharp structures in the spectra corresponding to critical points. The large values of the broadening factors for both path 2 and path 3 suggest that they describe processes that take place further from the critical points. Moreover, as was already discussed in the previous paragraph, at least part of the path 2 seems to model indirect transitions instead of the direct ones. This is supported by unjustifiably large value of the Rydberg constant for this path. The transition strength for path 2 looks almost like Gaussian broadened discrete spectrum at energy E1L . In order to verify that our model does not introduce any artifacts to dielectric function we can compare the calculated dielectric function with the dielectric function obtained by single-wavelength method. In this method the optical constants of c-Si are determined point-by-point from the ellipsometric data measured for one wavelength and one or several angles of incidence. The influence of the overlayer is taken into account in the single-wavelength method. The dielectric functions and their second derivatives obtained at 10 and 300 K are plotted in Figs. 10 and 11 . It can be seen that the values obtained by the single-wavelength method are in a good agreement with the values calculated by our model. So far, the exact positions of critical points were determined from the positions of peaks in the second derivative of the dielectric function. This is not necessary in this work since the positions of the critical points are found by fitting the experimental data. The positions of the critical points are also displayed in Figs. 10 and 11 . Since the experimental data in the range above 6.5 eV were available from VUV ellipsometer (Fig. 6) and EUV spectrophotometer at Elettra synchrotron (Fig. 9), we were able to consider also   high energy critical points E3 , E2X , E3X and E3X . The presence of these high energy critical points can be seen in Fig. 9. Unfortunately, the VUV data were measured only above the room temperature and EUV data were measured only at the room temperature, thus the high energy part of the model will be accurate only at the room temperature. Since the data up to 40 eV were available, it was possible to determine the effective number of valence electrons. This approach is different form the one used in our previous work [7], where the number of valence electrons was fixed to the estimated value 4.12. The effective number of valence electrons that we found was strongly dependent on the behavior of the dielectric function in the

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

Fig. 10. Spectral dependencies of the dielectric function and its second derivative determined by the single-wavelength method at 300 K. The data calculated by the presented model are displayed for comparison. The vertical bars represent positions of critical point energies labeled in graph. The denotation of critical point energies by the points of highest symmetry , X, L is omitted in this figure.

Fig. 11. The same as Fig. 10 but for 10 K.

limit of high energies, i.e. on the form of high energy excitations contribution. Two variants of the Tauc–Lorentz model were tried. The first variant, which used the Jellison–Modine parametrization, gave us the value nve = 4.51. The second variant, which used the Campi–Coriasso parametrization, gave us equally good fit with the lower value nve = 4.37. We believe that both the values obtained are too large, which is a result of 1/E2 (classical) decay of the transition strength function at high energies in the TL model. 7. Comparison with previous works The monocrystalline silicon (c-Si) is one of the most studied materials because of its use in microelectronics and also because samples of high purity become easily available. The optical constants of silicon have been studied many times in the literature. The works investigating spectral dependencies of optical constants can be divided into two groups, the first group uses single-wavelength method and the second group uses dispersion models.

415

In the single-wavelength method the optical constants are determined point-by-point from measurements (transmittance, reflectance, ellipsometry or minimal deviation method) performed at fixed wavelength. The transmittance is used mostly in the region of phonon absorption where the absorption coefficient can be easily determined if the thickness of the c-Si slab is known [34]. The reflectance can be used in the region of interband electronic transitions. Since the optical constants in this method has to be calculated using Kramers–Kronig analysis it is necessary to perform the measurement up to sufficiently high photon energy or extrapolate the data [35]. Advent of ellipsometry allowed to determine the optical constants directly from ellipsometric measurements at one wavelength [36,37,23,38,39]. The minimal deviation method uses the prism to measure the refractive index with high accuracy, but is limited only to the region in which the c-Si is transparent [18]. The disadvantage of single-wavelength methods is that it is difficult to compensate for the influence of the native overlayer, which is the main source of systematic errors in optical constants. The collection of optical constants obtained by the single-wavelength methods published by Palik in [40] will be used as a reference to which we will compare our results. The extinction coefficient in the region of partial transparency can be obtained from attenuation of interference fringes on thin slabs or films. If the refractive index and thickness of c-Si are known with high accuracy then this method could be very precise. The data recently obtained by ˇ [41] on silicon-on-insulator and homoepitaxial Humlíˇcek and Sik samples will be compared with our data in the region of indirect interband transitions above the bandgap. The methods that use the dispersion model are based on the interpretation of experimental data using the least squares method. In these methods the data measured over the entire spectral range are processed together. These methods have several advantages over single-wavelength methods: the influence of the noise is suppressed, the influence of the overlayer can be compensated more accurately, and if the dispersion model is Kramers–Kronig consistent then some of systematic errors can be eliminated. The easiest way to model c-Si optical constants in the regions of interband electronic transitions and phonon absorption is to use model of several damped harmonic oscillators (Lorentz model). The damped harmonic oscillator can be interpreted as a Lorentzian broadened discrete transition (delta function) and therefore it is not appropriate for the description of critical points where electronic and two-phonon JDOS is described by van Hove singularities. The delta function is a good approximation only in the case of 1D M1 critical point, discrete excitonic states or for the description of one-phonon absorption. The overview of several models that can be used to model the c-Si optical constants are in Table 5. The first practically useful model was introduced by Adachi [42,43] who based his approach on Lorentzian broadened JDOS functions corresponding to van Hove singularities. This approach provides Kramers–Kronig consistent description of direct electronic transitions, but cannot be used for the description of the indirect transition between the bandgap and critical point E 0 . In order to solve this problem Adachi introduced contribution to the imaginary part of the dielectric function, but left the real part of the dielectric function unchanged thus the resulting model was not Kramers–Kronig consistent. The correct solution to this problem can be based on models combining Tauc law and classical decay F ∝ 1/E2 for high energies, i.e. Jellison–Modine model [15], Campi–Coriasso [14] or similar models suggested in our work [10]. In principle it should be possible to use the Adachi model for modeling of two-phonon absorption bands but we are not aware of any work that would do this. Although the Lorentz broadening leads to analytic expressions, it is not suitable for broadening of direct electronic transitions. As was noted by Kim et al. [30] the Gaussian broadening is more

416

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

Table 5 Overview of dispersion models of crystalline solids. Acronyms: VHS, van Hove singularities; Ex, excitonic effects; P, Polynomials; L, Lorentzian; G, Gaussian, PG, pseudo Gaussian; ı, delta function. Model

Lorentz Adachi Kim et al. Tanguy Herzinger–Johs Our works

Refs.

e.g. [11] [42,43] [30] [44,45] [31,46] [10,24,7,8]

Direct electronic transitions

Indirect electronic transitions

Phonon absorption

JDOS

Broadening

Interband

Intraband

JDOS

Broadening

ı VHS VHS VHS, Ex P VHS, Ex

L L PG L G G

– Yes – – – Yes

– – – – – Yes

ı – – – – VHS, ı

L – – – – Gaussian

suitable in this case. The model that they proposed used pseudoGaussian broadening that mimics the Gaussian broadening and allows writing the broadened functions in an analytical form. The Adachi model was further extended by Tanguy [44,45] who added corrections for many-body effects. The corrections consist in adding the contribution representing the discrete electronic transitions corresponding to discrete excitonic bound states and modification of JDOS functions between critical points M0 and M1 corresponding to unbound states. As in the Adachi model, the Lorentzian broadening was utilized. As was already mentioned, the Lorentzian broadening is inappropriate for direct electronic transitions and Gaussian broadening should be used if high accuracy of the model is required. This was done in the model proposed by Herzinger and Johs [31,46] in which the JDOS function between the critical points is not modeled by functions corresponding to van Hove singularities but it is approximated by polynomial functions. The use of polynomial functions has the advantage that the Gaussian broadened functions can be written in an analytical form using special functions which can be efficiently calculated with arbitrary precision. This model is used to model the direct electronic transitions although it could be used also for two-phonon absorption. The excitonic effects are not incorporated in this model. The optical constants of c-Si published in [31] will be used as a third reference to which we will compare our results. Our model is based on similar ideas as the previous models, and it has several advantages: (i) The model distinguishes direct and indirect electronic transitions. The direct electronic transitions and two-phonon absorption use Gaussian broadened JDOS functions corresponding to the combination of 2D and 3D van Hove singularities. Note that the Gaussian broadening and Kramers–Kronig integration must be done numerically [24]. The indirect electronic transitions are modeled by two absorption edges representing electronic excitations with simultaneous creation or annihilation of phonons. (ii) Although all of the above models can be made temperature dependent by means of temperature dependent dispersion parameters, none of them incorporates temperature dependency based on physical principles. For example, in our model the phonon occupation number controls the transition strength (probability of transitions) of multi-phonon absorptions and also phonon assisted electronic excitations (indirect transitions). This feature of our model is ensured by the formalism based on sum-rule normalization of individual contributions. (iii) The direct electronic excitations contain the same corrections as used by Tanguy for many-body effects around the critical point M0 but with Gaussian broadening. Moreover, the influence of many-body effects in the vicinity of the critical points M3 in 3D and M2 in 2D is also taken into account [47]. (iv) At each critical energy all possible 3D and 2D van Hove singularities are considered. (v) Our model also contains the phonon and free carrier parts [7,8]. The comparison of our model with Palik’s, Herzinger’s and Frey’s data in the region of transparency is plotted in Fig. 12. Because the the Frey’s data were included in the fitted data, their difference from

Fig. 12. Comparison of the refractive index (top) and difference of Frey’s data [18] from our model (bottom) in the region of transparency. The Frey’s data for 300 K were extrapolated by (3). Herzinger’s [31] and Palik’s [40] data are added for comparison.

Fig. 13. Temperature dependence of refractive index at wavelengths 2 ␮m (top) and difference between extrapolated Frey’s data [18] and our model (bottom). Extrapolation was performed by (3).

our model, especially at 2 ␮m (0.6199 eV), is very small. The temperature dependencies of extrapolated Frey’s data and our model at wavelength 2 ␮m is plotted in Fig. 13. Note that the undulation in differences between extrapolated Frey’s data and our model is given by systematic error of our model, whereas the noise is given by the precision of numerical integration. Our model exhibits the

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

Fig. 14. Comparison of spectral dependencies of optical constants calculated using our model at 300 K and Palik’s [40], Frey’s [18] and Herzinger’s [31] data in the region of phonon absorption.

systematic error which is approximately five times larger than the accuracy of Frey’s data (see Fig. 12). This indicates that the distribution of transition strength in our model is not perfect. The Herzinger’s data are approximately ın = 2.0 × 10−3 above the Frey’s data at 295 K and ın = 1.2 × 10−3 above extrapolated Frey’s data at 300 K. Note that the exact temperature is not stated in the work [31], but the relative error better than 1 × 10−3 is a very good result for ellipsometric measurements. The Palik’s data show larger differences (ın =−4.5 × 10−3 at wavelength 2 ␮m and temperature of Frey’s data 295 K). This should be expected since the data are older and the single-wavelength method was used. The comparison of our model with Palik’s data in the region of phonon absorption is plotted in Fig. 14. The Herzinger’s and Frey’s data are available only in the region of transparency (the right part of the plot) and were already discussed in the previous paragraph. The Palik’s data for the extinction coefficient contain all two-phonon and three-phonon structures that we have in our model except for LO−TA structure at 160 cm−1 . The higher value of the peak at 610 cm−1 in Palik’s data could be a systematic error or it could be caused by one-phonon absorption corresponding to higher carbon content in c-Si. The comparison of our model with Palik’s, Herzinger’s and Humlíˇcek’s data in the region of indirect electronic transitions is plotted in Fig. 15. The values of the extinction coefficient calculated by our model lies between Palik’s and Herzinger’s data. The comparison of our model and Herzinger’s data in the vicinity of bandgap is in Fig. 7. The two-edge absorption structure is clearly visible in our model while this structure is absent in Herzinger’s data. This is because the model used by Herzinger does not include contribution from the indirect electronic transitions, thus the data in this region are modeled as direct electronic transitions. The comparison of our model with Palik’s and Herzinger’s data in the region of direct electronic transitions is plotted in Fig. 16. The maximal value 5.355 of the extinction coefficient at critical energy E2X is lower in our model than the values 5.435 (Palik) and 5.465 (Herzinger). This lower value is supported by the results obtained by single-wavelength method (see Fig. 10) and also by the fact that if the weight of Frey’s data is increased then this value of extinction coefficient decreases. On the other hand the higher value of peak value of the extinction coefficient was obtained by Herzinger et al. using multi-sample multi-wavelength method based on simultaneous processing of data measured on sample covered by native

417

Fig. 15. Comparison of spectral dependencies of optical constants calculated using our model at 300 K and Palik’s [40], Frey’s [18], Herzinger’s [31] and Humlíˇcek’s [41] data in the region of indirect electronic transitions.

Fig. 16. Comparison of spectral dependencies of optical constants calculated using our model at 300 K and Palik’s [40] and Herzinger’s [31] data in the region of direct electronic transitions.

overlayer as well as on samples covered by thermal oxide layers. The higher value 5.448 was also obtained in our previous work [48] where multi-sample single-wavelength method was utilized for similar sample set as in Herzinger et al. work. The comparison of our model with Palik’s data in the region of high energy excitations of valence electrons and core electron excitations is plotted in Fig. 17. We do not have experimental data above 40 eV thus our model in this region is purely theoretical. The positions of edges of core electron excitations in this region are fixed at values taken from [49], the ratios of their strengths are taken from Shiles work for aluminum [50] with absolute values of strengths determined from the sum rule. Moreover, our model assumes classical asymptotic behavior, i.e. k ∝ 1/E3 . The discrepancy between our model and Palik’s data in the region 20–50 eV can be explained by low sensitivity of our near-normal reflectance measurements (see Fig. 9). This discrepancy could explain higher than expected value of effective number of valence electrons and since it contributes to Kramers–Kronig integral it could also be

418

D. Franta et al. / Applied Surface Science 421 (2017) 405–419

c-Si commonly used in microelectronics industry is significantly less expensive than the pure float zone silicon. At this moment we performed these extensions for slightly phosphorus doped CZ cSi. So far, only the part corresponding to interstitial oxygen was published in [7]. This paper could be understood as a part of broader work leading to complete presentation of temperature dependent model of c-Si accurate enough for use in spectroscopic ellipsometry and spectrophotometry. Although the main goal of this paper is the presentation of the theoretical background of our model it is clear that our model is already sufficient for spectroscopic purpose from far IR to vacuum UV. We are planning to further improve the precision of our model in the future by processing additional experimental data, for example, transmittance measurements around bandgap for different wafer (slab) thicknesses or ellipsometric data with different thicknesses of oxide overlayers. Acknowledgements

Fig. 17. Comparison of spectral dependencies of optical constants calculated using our model at 300 K and Palik’s [40] data in the region of high energy excitations of valence electrons and core electron excitations.

a reason for worse then ideal agreement of our model with dispersion of Frey’s data in the transparency region (see Fig. 12). In order to solve this problem it is necessary to add relevant experimental data and change the model of high energy excitations of valence electrons. The discrepancy between 1.0–1.8 keV could indicate deviation from the classical behavior. In order to explain it the relativistic corrections should be considered.

This work was supported by the projects ‘R&D center for low-cost plasma and nanotechnology surface modifications’ (CZ.1.05/2.1.00/03.0086) and ‘CEITEC – Central European Institute of Technology’ (CZ.1.05/1.1.00/02.0068) from European Regional Development Fund, projects LO1411 and LQ1601 (NPU I) funded by Ministry of Education, Youth and Sports of the Czech Republic and project TA02010784 of the Technology Agency of the Czech Republic. The experiment at BEAR was carried out in proposal 20155284 at Elettra. The research leading to these results has received funding from the European Community’s Seventh Framework Programme (FP7/2007–2013) under grant agreement no. 312284.

8. Conclusion

References

The temperature dependent dispersion model of float zone crystalline silicon was presented in this paper. The model is based on data in the temperature range 5–600 K in the infrared region and 10–700 K in the visible and ultraviolet region. The current model in IR region is an extension of the model presented in [21], which worked well only at the room temperature. The extended IR model correctly distinguishes between one-, two- and threephonon absorption processes and it should give accurate dielectric response in the range 0–600 K. The theoretical background for valence electronic excitations is presented in detail in the theoretical part of this paper. Although the separation of direct and indirect electronic transitions may not be accurate from physical point of view, the agreement between modeled and experimental ellipsometric and spectrophotometric quantities was excellent. We would like to address the problem of incorrect separation of direct and indirect electronic transitions in our future work. Since the experimental data up to 40 eV were available, we were able to determine the value of effective number of valence electrons. Two slightly different values were obtained for Jellison–Modine (nve = 4.51) and Campi–Coriasso (nve = 4.373) parametrization of Tauc–Lorentz model that were utilized for the description of high energy excitations. Both of these values are higher than expected, which is probably caused by classical asymptotic behavior at high energies. In reality the decay should be faster, but in order to solve this problem experimental data for energies above 40 eV are needed. The current model is concerned with float zone c-Si, it can be extended for Czochralski (CZ) silicon if the influence of interstitial and precipitated oxygen is taken into account. Also the model can be relatively easily extend to include the influence of dopants. These extensions of our model are very important because the doped CZ

[1] S.M. Hu, Infrared absorption spectra of SiO2 precipitates of various shapes in silicon: calculated and experimental, J. Appl. Phys. 51 (1980) 5945–5948. [2] R.C. Newman, Defects in silicon, Rep. Prog. Phys. 45 (1982) 1163–1210. [3] R. Jones, A. Umerski, S. Öberg, Ab initio calculation of the local vibratory modes of interstitial oxygen in silicon, Phys. Rev. B 45 (1992) 11321–11323. [4] B. Pajot, O-related IR absorption in c-Si, in: R. Hull (Ed.), Properties of Crystalline Silicon, The Institution of Engineering and Technology, 1998, pp. 492–498. [5] R.C. Newman, Oxygen diffusion and precipitation in Czochralski silicon, J. Phys.: Condens. Matter 12 (2000) R335–R365. [6] B. Pajot, B. Clerjaud, Optical Absorption of Impurities and Defects in Semiconducting Crystals, Springer Series in Solid-State Sciences 169, Springer-Verlag, 2013. [7] D. Franta, D. Neˇcas, L. Zajíˇcková, I. Ohlídal, Utilization of the sum rule for construction of advanced dispersion model of crystalline silicon containing interstitial oxygen, Thin Solid Films 571 (2014) 490–495. [8] D. Franta, D. Neˇcas, L. Zajíˇcková, I. Ohlídal, Dispersion model of two-phonon absorption: application to c-Si, Opt. Mater. Express 4 (8) (2014) 1641–1656. [9] D. Franta, et al. Software for optical characterization newAD2, http://newad. physics.muni.cz. [10] D. Franta, D. Neˇcas, L. Zajíˇcková, Application of Thomas-Reiche-Kuhn sum rule to construction of advanced dispersion models, Thin Solid Films 534 (2013) 432–441. [11] P.Y. Yu, M. Cardona, Fundamentals of Semiconductors, Springer, 2001. [12] D. Franta, D. Neˇcas, L. Zajíˇcková, I. Ohlídal, J. Stuchlík, D. Chvostová, Application of sum rule to the dispersion model of hydrogenated amorphous silicon, Thin Solid Films 539 (2013) 233–244. [13] D. Franta, D. Neˇcas, L. Zajíˇcková, I. Ohlídal, J. Stuchlík, Advanced modeling for optical characterization of amorphous hydrogenated silicon films, Thin Solid Films 541 (2013) 12–16. [14] D. Campi, C. Coriasso, Prediction of optical properties of amorphous tetrahedrally bounded materials, J. Appl. Phys. 64 (1988) 4128–4134. [15] G.E. Jellison Jr., F.A. Modine, Parameterization of the optical functions of amorphous materials in the interband region, Appl. Phys. Lett. 69 (1996) 371–373. [16] S. Nannarone, F. Borgatti, A. DeLuisa, B.P. Doyle, G.C. Gazzadi, A. Giglia, P. Finetti, N. Mahne, L. Pasquali, M. Pedio, G. Selvaggi, G. Naletto, M.G. Pelizzo, G. Tondello, The BEAR beamline at Elettra, AIP Conf. Proc. 705 (2004) 450–453. [17] L. Pasquali, A.D. Luisa, S. Nannarone, The UHV experimental chamber for optical measurements (reflectivity and absorption) and angle resolved photoemission of the BEAR beamline at Elettra, AIP Conf. Proc. 705 (2004) 1142–1149.

D. Franta et al. / Applied Surface Science 421 (2017) 405–419 [18] B.J. Frey, D.B. Leviton, T.J. Madison, Temperature-dependent refractive index of silicon and germanium, in: Optomechanical Technologies for Astronomy, Vol. 6237 of Proc. SPIE, 2006, p. 62732J. [19] D. Franta, I. Ohlídal, P. Klapetek, Analysis of slightly rough thin films by optical methods and AFM, Mikrochim. Acta 132 (2000) 443–447. ˇ P. Pokorny, ´ Optical [20] D. Franta, I. Ohlídal, D. Neˇcas, F. Viˇzd’a, O. Caha, M. Hason, characterization of HfO2 thin films, Thin Solid Films 519 (2011) 6085–6091. [21] D. Franta, D. Neˇcas, I. Ohlídal, Universal dispersion model for characterization of optical thin films over wide spectral range: application to hafnia, Appl. Opt. 54 (2015) 9108–9119. [22] P. Lautenschlager, P.B. Allen, M. Cardona, Temperature dependence of band gaps in Si and Ge, Phys. Rev. B 31 (1985) 2163–2171. [23] P. Lautenschlager, M. Garriga, L. Vi na, M. Cardona, Temperature dependence of the dielectric function and interband critical points in silicon, Phys. Rev. B 36 (1987) 4821–4830. [24] D. Franta, D. Neˇcas, L. Zajíˇcková, I. Ohlídal, Broadening of dielectric response and sum rule conservation, Thin Solid Films 571 (2014) 496–501. [25] R. de Laer Kronig, On the theory of dispersion of X-ray, J. Opt. Soc. Am. Rev. Sci. Instrum. 12 (1926) 547–557. [26] F. Wooten, Optical Properties of Solids, Academic Press, New York, 1972. [27] H. Watanabe, N. Yamada, M. Okaji, Linear thermal expansion coefficient of silicon from 293 to 1000 K, Int. J. Thermophys. 25 (2004) 221–236. [28] D.Y. Smith, Dispersion theory, sum rules, and their application to the analysis of optical data, in: E.D. Palik (Ed.), Handbook of Optical Constants of Solids, Vol. 1, Academic Press, 1985, pp. 35–68. [29] C. Kittel, Introduction to Solid State Physics, 5th edition, Wiley, New York, 1976. [30] C.C. Kim, J.W. Garland, H. Abad, P.M. Raccah, Modeling the optical dielectric function of semiconductors: extension of the critical-point parabolic-band approximation, Phys. Rev. B 45 (1992) 11749–11767. [31] C.M. Herzinger, B.D. Johs, W.A. McGahan, J.A. Woollam, W. Paulson, Ellipsometric determination of optical constants for silicon and thermally grown silicon dioxide via a multi-sample, multi-wavelength, multi-angle investigation, J. Appl. Phys. 83 (1998) 3323–3336. [32] G.G. Macfarlane, V. Roberts, Infrared absorption of silicon near the lattice edge, Phys. Rev. 98 (1955) 1865–1866. [33] K. Kondo, A. Moritani, Symmetry analysis of the E2 structures in Si by low-field electroreflectance, Phys. Rev. B 15 (1977) 812–815. [34] F.A. Johnson, Lattice absorption bands in silicon, Proc. Phys. Soc. 73 (1959) 265.

419

[35] H.R. Philipp, E.A. Taft, Optical constants of silicon in the region 1 to 10 eV, Phys. Rev. 120 (1960) 37–38. [36] D.E. Aspnes, A.A. Studna, Dielectric functions and optical parameters of Si, Ge, GaP, GaAs, GaSb, InP, InAs, and InSb from 1.5 to 6.0 eV, Phys. Rev. B 27 (1983) 985–1009. [37] G.E. Jellison Jr., F.A. Modine, Optical functions of silicon between 1.7 and 4.7 eV at elevated temperatures, Phys. Rev. B 27 (1983) 7466–7472. [38] G. Vuye, S. Fisson, V. Nguyen Van, Y. Wang, J. Rivory, F. Abelès, Temperature dependence of the dielectric function of silicon using in situ spectroscopic ellipsometry, Thin Solid Films 233 (1993) 166–170. ˇ J. Hora, J. Humlíˇcek, Optical functions of silicon at high temperatures, J. [39] J. Sik, Appl. Phys. 84 (1998) 6291–6298. [40] D.F. Edwards, Silicon (Si), in: E.D. Palik (Ed.), Handbook of Optical Constants of Solids, Vol. 1, Academic Press, New York, 1985, pp. 547–569. ˇ Optical functions of silicon from reflectance and [41] J. Humlíˇcek, J. Sik, ellipsometry on silicon-on-insulator and homoepitaxial samples, J. Appl. Phys. 118 (2015) 195706. [42] S. Adachi, Model dielectric constants of GaP, GaAs, GaSb, InP, InAs, and InSb, Phys. Rev. B 35 (1987) 7454–7463. [43] S. Adachi, Model dielectric constants of Si and Ge, Phys. Rev. B 38 (1988) 12966–12976. [44] C. Tanguy, Optical dispersion by Wannier excitons, Phys. Rev. Lett. 75 (1995) 4090–4093. [45] C. Tanguy, Complex dielectric constant of two-dimensional Wannier excitons, Solid State Commun. 98 (1996) 65–68. [46] B.D. Johs, C.M. Herzinger, J.H. Dinan, A. Cornfeld, J.D. Benson, Development of a parametric optical constant model for Hg1−x Cdx Te for control of composition by spectroscopic ellipsometry during MBE growth, Thin Solid Films 313–314 (1998) 137–142. ´ J. Sak, Excitonic effects in the interband absorption of [47] B. Velicky, semiconductors, Phys. Status Solidi 16 (1966) 147–157. [48] I. Ohlídal, D. Franta, E. Pinˇcík, M. Ohlídal, Complete optical characterization of the SiO2 /Si system by spectroscopic ellipsometry, spectroscopic reflectometry and atomic force microscopy, Surf. Interface Anal. 28 (1999) 240–244. [49] G.P. Williams, X-Ray Data Booklet, http://xdb.lbl.gov/Section1/Sec 1-1.html. [50] E. Shiles, T. Sasaki, M. Inokuti, D.Y. Smith, Self-consistency and sum-rule tests in the Kramers–Kronig analysis of optical data: applications to aluminum, Phys. Rev. B 22 (1980) 1612–1628.