Thermal Conductivity Analysis of Amorphous Silicon Formed by Natural Cooling: A Molecular-dynamics Study

Article information

J. Korean Ceram. Soc.. 2016;53(3):295-300

Publication date (electronic) : 2016 May 31

doi : https://doi.org/10.4191/kcers.2016.53.3.295

Byoung Min Lee

School of Energy, Materials & Chemical Engineering, Korea University of Technology and Education, Cheonan, 31253, Korea

^†Corresponding author: Byoung Min Lee, E-mail: bmlee@koreatech.ac.kr, Tel: +82-41-560-1365 Fax: +82-41-560-1360

Received 2016 March 25; Revised 2016 May 10; Accepted 2016 May 18.

Abstract

To investigate the thermal conductivity and the structural properties of naturally cooled excimer-laser annealed Si, molecular-dynamics (MD) simulations have been performed. The thermal conductivity of crystalline Si (c-Si) was measured by direct method at 1000 K. Steady-state heat flow was measured using a stationary temperature profile; significant deviations from Fourier’s law were not observed. Reliable processes for measuring the thermal conductivity of c-Si were presented. A natural cooling process to admit heat flow from molten Si (l-Si) to c-Si was performed using an MD cell with a size of 48.9 × 48.9 × 97.8 Å³. During the cooling process, the temperature of the bottom 10 Å of the MD cell was controlled at 300 K. The results suggest that the natural cooling system described the static structural property of amorphous Si (a-Si) well.

Keywords: MD simulations; Silicon; Thermal conductivity; Natural cooling; Static structure factors

1. Introduction

Polycrystalline silicon (poly-Si) has been an attractive material for thin film transistors (TFTs) in active matrix liquid crystal displays1) as well as thin film solar cells.2) Due to the technological significance of large-grained poly-Si thin films, a lot of research has been devoted to studying the kinetics and the transformation mechanisms of the excimer-laser annealing (ELA) process. For applications of TFTs and solar cells, because residual grain boundaries in recrystallized Si degrade the electronic characteristics, especially the carrier mobility of TFTs, large-grained poly-Si is desired to realize high performance of devices.3–6)

The ELA process involves rapid melting and solidification as well as lateral and vertical heat flow in a Si thin film when that film is irradiated with an incident beam, as shown in Fig. 1. It has been reported that heat flow into the substrate from molten Si is a predominant mechanism of heat removal from Si thin film.7) At present, the most widely used method to prepare poly-Si on SiO₂ glass substrate is excimer-laser crystallization. Even when the molten Si is heated above 900 °C, this high temperature can only be sustained for a short time (about 20 ns). Therefore, the thermal strain originating due to the high energy of the excimerlaser does not lead to severe damage of the SiO₂ glass substrate. Laser crystallization is a very fast process that can melt amorphous Si in a very short time and produce large grained poly-Si with low defect density.8)

Fig. 1

Schematics of nucleation and crystal growth of Si in TFTs during excimer-laser annealing process.

The thermal conductivity of the substrate materials is a key factor in determining the cooling rate of excimer-laser annealed Si because the nucleation and crystallization processes of Si are controlled by the unidirectional heat flow from molten Si to the substrate material. Many researchers have reported that MD simulations can be used to interpret the thermal and structural properties of Si.9) Therefore, the implementation of MD simulation to analyze the thermal conductivity under unidirectional heat flow will be a good technique to obtain the time-averaged atomic quantities and fluctuations that are needed to calculate the thermal properties of Si.

In order to obtain the thermal and structural properties of Si in atomic scale, we employed the natural cooling system introduced by Lee et al.10–15) They applied a combination NVT and NVE ensemble to reproduce the ELA process using MD simulation. Nucleation and crystallization processes of laser annealed Si, which result in heat flow from high energy regions that absorb the laser energy to low energy regions that act as heat sinks, were successfully reported. In the present study, MD simulations have been carried out to understand the thermal properties of c-Si and the effects of the microstructural characteristics of a-Si in the ELA process. The feasibility of performing thermal conductivity computing processes using MD simulations will be discussed.

2. Calculation Details

The most important element of MD simulations is the selection of the potential that describes the interactions among atoms. During calculations, the interatomic forces for Si were calculated using the Tersoff potential. The Tersoff potential consists of the many-body interactions included in the bond order term for Si. This potential has been used to study the lattice dynamics,16) thermo-mechanical properties,17) point defects,18) and the liquid and amorphous phases of Si.19,20) The Tersoff potential describes well the properties of both crystalline and liquid phases of Si.21) Although the Tersoff potential gives rise to a high melting temperature about 2500 K, it describes the thermal and structural properties of solid and liquid Si well.16)

Herein, the thermal conductivity of c-Si was obtained by means of a direct method. Fig. 2 shows the whole modeling procedure used to calculate the thermal conductivity of c-Si. The direct method is a non-equilibrium method that relies on imposing a temperature gradient across the MD cell. Therefore, it is similar to the standard experimental conditions, which involve measuring the stationary heat flux necessary to maintain a temperature gradient. This method is only valid when a stable, linear temperature profile is obtained in the system. In other words, it requires first the establishment of a steady state, and then that heat flow in the direction of the temperature gradient should be uniform. If a linear temperature profile is established, the thermal conductivity can be obtained using Fourier’s law of heat flow J,

Fig. 2

Schematic representation of periodic simulation cell. The cell was divided into slices to determine the thermal conductivity of c-Si. There is a layer of thickness δ at z = −Lz/4 into which energy Δɛ is added; in the layer at z = +Lz/4, energy Δɛ is removed at each MD step.

(1) J=ΔɛAΔt

where Δt is the time step used in the MD simulation, A is the cross section area of the system, and Δɛ is a constant energy that is added to and removed from the heat source and heat sink layer. The thermal conductivity can be expressed as:

(2) κ=-J∂T/∂x

Here, ∂T/∂x is the temperature gradient.

During the simulations, a steady-state heat flow was generated for a given sampling time by scaling the velocity of the atoms, which are located in the heat source and the heat sink layer, as shown in Fig. 2. The MD cells, which have a length of a few tens of nm in the direction of the heat flow, were used. The equations of motion were integrated using the velocity Verlet algorithm. All the atoms were placed in their initial positions and provided with initial velocities at 1000 K. We considered the heat source and heat sink layer as having a 2 Å thickness in the direction perpendicular to the heat current, which is located at z = Lz/4 and z = 3Lz/4 in the MD cell of length Lz. At each calculation, the atoms inside the slices were determined and constant energy was added to a thin slice centered at z = Lz/4 (heat source) and the same energy was removed from a layer centered at z = 3Lz/4 (heat sink) by rescaling the particle velocity at each MD time step of 0.5 fs. By imposing the heat in this manner, a constant heat flux per unit area J_z, which is equal to Δɛ / (2L²Δt), was obtained (the factor 2 comes from the fact that the heat flux coming from the hot layer divides equally into two parts to reach the cold layer). When the system achieved a steady state, the appropriate temperature gradient was obtained between the heat source and the heat sink layer. The calculations for the thermal conductivity of c-Si were performed after measuring the temperature gradients.

The structural property of a-Si prepared by natural cooling was investigated using the MD cell with the size of 48.9 × 48.9 × 97.8 Å³. MD simulations were performed under constant NVT conditions (i.e., constant number of particles N, constant volume V, and constant temperature T), using an MD cell with periodic boundary conditions. Periodic boundary conditions were employed in the [100] and [010] directions. The solid-liquid interface was prepared by setting the temperature of the atoms at 300 K and 3500 K for Z ≤ 35 Å and Z > 35 Å of the MD cell, respectively. Temperature was controlled using the ordinary Langevin equations; we employed the scheme developed by Gunsteren and Berendsen for the numerical integration of the Langevin equations. The Langevin equation for the motion of atom i is:

(3) mri(t)=Fi(t)-mγri(t)+Ri(t)

where F_i(t) is the sum of all forces exerted on atom i by other atoms and can be derived from the potential, which is dependent on the coordinates of all N particles. m_i is the mass of the atoms, and γ is the friction constant. The time step for the integration was set at 0.002 ps and the friction coefficient was set at 5 ps⁻¹. For the natural cooling processes employed in this study, the combination NVT and NVE ensemble (i.e., constant number of particles N, constant volume V, and constant energy E) was applied. During cooling, only the temperature of the bottom 10 Å of the MD cell was controlled at 300 K. The simulation procedures are described in detail elsewhere.10–13)

3. Results and Discussion

In the case of crystalline materials, heat is carried by phonons, which can influence the transport of heat. Usually, long simulation times are needed to achieve a smooth temperature profile because phonon-phonon scattering events in a solid are far less frequent than particle collisions in a fluid and are amorphous when phonon-phonon scattering is the dominant mechanism for thermal resistance.22) Therefore, the simulation times of the 2 × 10⁶ MD steps with a time step of 0.5 fs were applied to achieve a smooth temperature profile.

As mentioned above, the steady-state heat flow is an important criterion for the reasonability of Fourier’s law when calculating the thermal conductivity. The steady-state heat flow can be measured using a stationary temperature profile as a function of time. Fig. 3 displays the time evolution of the time-averaged temperature at z = 59.9 nm and z = 99.9 nm at 1000 K. At the specific position z, the time-averaged temperature in a thin slice is given by

Fig. 3

Time evolution of time-averaged temperature for slices at 20 nm away from heat sink. Gray and black curves are the temperatures at Z = 59.9 nm and Z = 99.9 nm for a 4 × 4 × 196 system at 1000 K.

<T(z)>M=1M∑m=0M-1TN-m(z)

where <T(z)>_M is the time-averaged temperature at z, M is the time step, and T_N_-_m(z) is the instantaneous temperature at z for the time step N-m. Here, N is the total number of time steps. To obtain steady-state heat flow, the temperatures in both slices were first averaged over 1000 MD steps to eliminate the large temperature fluctuations. The temperatures at both positions fluctuate significantly for a few ps and stabilize after 150 ps. The steady state was achieved in 1 ns. We found a difference of about 2 K between the two layers.

The temperature gradient was calculated by dividing the MD cell into slices with thickness of 2 Å along the heat flow direction. The average temperature in a thin slice was computed at each MD step. As shown in Fig. 2, the presence of the heat source and heat sink and the application of periodic boundary conditions in all directions generated a heat current in two opposite directions. The thermal conductivity (κ) is calculated by linear fitting of the temperature gradient in the intermediate region of the average temperature profile.

To determine whether Fourier’s law is obeyed, it is essential to check whether the value of κ does in fact depend on the choice of Δɛ. Because the heat current is proportional to Δɛ/A, this is performed by computing for several different energy values from 0.0003 to 0.01 eV. Fig. 4 shows the effect of changing Δɛ for a 4 × 4 × 96 system at 1000 K. Although these results show a bit variation of with Δɛ, the values of κ are not appreciably affected by the changes of Δɛ. No significant deviations from Fourier’s law are apparent in Fig. 4. Extrapolation of the system size L_z to infinity yields an accurate κ of 58 W/mK at 1000 K.

Fig. 4

Effect of changing Δɛ for a 9 × 9 × 96 system at 1000 K.

Because the thermal conductivity can be calculated well using MD simulations, we performed a natural cooling computation to prepare a-Si. In this way, the heat is moved from l-Si to the temperature-controlled region. The heat and temperature changes during natural cooling are shown in Fig. 5 and 6. Rapid heat release by heat diffusion from l-Si to the temperature-controlled region is observed and, as a result, the heat flow through c-Si induces a considerable temperature gradient in the l-Si region. A temperature transition layer, with a thickness of about 10 Å at the solid-liquid interface after melting, can be observed in Fig. 6, but this soon disappears during natural cooling. This phenomenon is attributed to the temperature control of 10 Å of c-Si region and the heat flow through the c-Si region, both of which lead to the linear temperature gradient being at the interface, unlike previous results.10–12) Considering the fact that the heat release rate from l-Si into the temperature-controlled regions depends not only on the heat flow within the bulk l-Si but also on the thermal conductivity of c-Si, the relatively rapid cooling rate of ~10¹³ K/sec can be attributed to the high thermal conductivity of c-Si and the short distance between the heat removal regions and the solid-liquid interface. Because we consider that the MD cell used in this simulation is relatively too small to allow a correct interpretation of the cooling rate, it can be expected that significant improvements will be possible in the future.

Fig. 5

Temperature change (a) and variation of heat flow (b) of system during natural cooling. The temperature of bottom 10 Å of MD cell was set at 300 K.

Fig. 6

Local temperature changes obtained by mean kinetic energy during natural cooling with time.

For theoretical understanding of the thermal equilibrium during natural cooling, the velocity distribution of atoms (VDA) in the local area was investigated. Fig. 7 shows the changes of VDA in the region of 90 ~ 95 Å at two different times of 160 and 170 fs, respectively, during natural cooling. The dashed lines indicate the VDA obtained theoretically using the Maxwell-Boltzmann function with the averaged local temperatures. The Maxwell-Boltzmann distribution of VDA indicates that the system is in thermal equilibrium during natural cooling. The cooling rate of ~10¹³ K/sec leads to the formation of a-Si.

Fig. 7

Velocity distribution changes during natural cooling in the region of 90 ~ 95 Å at (a) t = 160 fs and (b) t = 170 fs, respectively.

In order to determine the dynamical properties of naturally cooled a-Si, the static structure factors S(k) are calculated by Fourier transforming of the pair-correlation functions g(r),10)

S(k)=1+4πρ0∫0∞r2{g(r)-1}sin(kr)krdr

Here, ρ₀ is the average number density of atoms and k is the wave number of the diffracted wave. The atoms lying over 50 Å of the MD cell were sampled to investigate the structural properties. It should be noted that the setting temperature for preparing l-Si is about 3500 K, with a density of 2.33 g/cm³. The k < 1.5 Å⁻¹ region contains truncation error and should be ignored. In the case of l-Si, it can be clearly seen that the first peak is asymmetric. This shows the existence of some covalent bonds and the reminiscent of the crystalline characteristics in the liquid.23) The peak shape resembles that of S(k) for a supercooled metallic liquid Si.24) This result supports the metallic behavior of the melt, as expected. With decreasing temperature, the first peak splits into two and these two peaks develop and appear at low temperature. The appearance of a new shoulder of the third peak is in good agreement with the experimental results.25) This spectrum of the peak indicates the short-range order of a-Si.

4. Conclusions

MD simulations have been performed to calculate the thermal conductivity of c-Si via the direct method. Calculation procedures for thermal conductivity were proposed and reliable results were obtained using Fourier’s law at 1000 K. The results for the steady-state heat flow and the effect of changing Δɛ confirmed that Fourier’s law was satisfied during the thermal conductivity calculation of c-Si. To prepare a-Si, natural cooling simulation was performed with a combination NVT and NVE ensemble. A temperature transition layer was observed at the interface after melting faded away during natural cooling, and a linear temperature gradient was found at the interface instead. The static structure of a-Si quenched with a cooling rate of ~10¹³ K/sec showed the appearance of a new shoulder on the third peak, suggesting that a-Si generated by natural cooling described the characteristic of experimentally prepared a-Si well. MD simulations using a large system, which could be compared with the film thickness used in the ELA process, are under study.

Fig. 8

Static structure factors S(k) of Si. Closed and open circles denote the static structure factors of a-Si naturally cooled at 300 K and l-Si prepared at 3500 K, respectively.

References

1. Samb ML, Jacques E, Belarbi K, Coulon N, Mohammed-Brahim T. μc-Si Thin Film Transistors with very Thin Active Layer. Solid-State Electron 89(1):128–33. 2013;

2. Orpella A, Martin I, Lopez-Gonzalez JM, Ortega P, Munoz J, Sinde DC, Voz C, Puigdollers J, Alcubilla R. Experimental Determination of Base Resistance Contribution for Point-like Contacted c-Si Solar Cells Using Impedance Spectroscopy Analysis. Sol Energy Mater Sol Cells 141(5):350–55. 2015;

3. Sposili RS, Im JS. Sequential Lateral Solidification of Thin Silicon Films on SiO₂. Appl Phys Lett 69(19):2864–66. 1996;

4. Hatano M, Moon S, Lee M. Excimer Laser-Induced Temperature Field in Melting and Resolidification of Silicon Thin Films. J Appl Phys 87(1):36–43. 2000;

5. Nerding M, Dasson R, Christiansen S, Kohler JR, Krinke J, Werner JH, Strunk H-P. Microstructure of Laser-Crystallized Silicon Thin Films on Glass Substrate. J Appl Phys 91(8):4125–30. 2002;

6. Lee M, Suzuki K, Grigoropoulos CP. Relationship between Fluence Gradient and Lateral Grain Growth in Spatially Controlled Excimer Laser Crystallization of Amorphous Silicon Films. J Appl Phys 88(9):4994–99. 2000;

7. Ishikawa K, Oazwa M, Oh CH, Matsumura M. Excimer-Laser-Induced Lateral-Growth of Silicon Thin-Films. Jpn J Appl Phys, 37(3A):731–36. 1998;

8. Boyce JB, Mei P, Fulks RT, Ho J. Laser Processing of Polysilicon Thin-Film Transistors: Grain Growth and Device Fabrication. Phys Stat Sol A 166(2):729–41. 1998;

9. Schelling PK, Phillpot SR, Keblinski P. Comparison of Atomic-Level Simulation Methods for Computing Thermal Conductivity. Phys Rev B 65(14):144306. 2002;

10. Lee BM, Seong BS, Munetoh S, Motooka T. Heat Flow and Structural Properties of Naturally Cooled a-Si: A Molecular Dynamics Study. J Korean Phys Soc 49(6):2353–61. 2006;

11. Lee BM, Munetoh S, Motooka T. Molecular Dynamics Study of Velocity Distribution and Local Temperature Change during Rapid Cooling Processes in Excimer-laser Annealed Silicon. Comp Mater Sci 37(3):198–202. 2006;

12. Lee BM, Baik HK, Seong BS, Munetoh S, Motooka T. Molecular-Dynamics Analysis of the Nucleation and Crystallization Process of Si. Phys B 392(2):266–71. 2007;

13. Lee BM. Crystal Growth and Nucleation in Liquid Silicon During Natural Cooling: A Molecular-Dynamics Simulation. Korean J Met Mater 45(12):681–87. 2007;

14. Lee BM, Kuranaga S, Munetoh T, Motooka T. Surface Nucleation of the (111) Plane of Excimer Laser Annealed Si on SiO₂ Substrates: A Molecular Dynamics Study. J Appl Phys 101(5):0543161–67. 2007;

15. Lee BM, Motooka T, Munetoh S. Molecular-Dynamics Simulations of Nucleation and Crystallization Processes of Laser Crystallized Poly-Si. J Phys: Condens Matter 20(5):055205–15. 2008;

16. Bazant MZ, Kaziras E. Modeling of Covalent Bonding in Solids by Inversion of Cohesive Energy Curves. Phys Rev Lett 77(21):4370–73. 1996;

17. Broughton JQ, Li XP. Phase Diagram of Silicon by Molecular Dynamics. Phys Rev B 35(17):9120–27. 1987;

18. Tersoff J. New Empirical Model for the Structural Properties of Silicon. Phys Rev Lett 56(6):632–35. 1986;

19. Luedtke WD, Lendman U. Preparation and Melting of Amorphous Silicon by Molecular-Dynamics Simulations. Phys Rev B 37(9):4656–63. 1988;

20. Porter LJ, Yip S, Yamaguchi M, Kaburaki H, Tang M. Empirical Bond-Order Potential Description of Thermodynamic Properties of Crystalline Silicon. J Appl Phys 81:96–106. 1997;

21. Ishimaru M, Yoshida K, Motooka T. Potential Model for Silicon Clusters. Phys Rev B 39(2):1212–18. 1989;

22. Mountain RD, MacDonald RA. Thermal Conductivity of Crystals: A Molecular-Dynamics Study of Heat Flow in a Two-Dimensional Crystal. Phys Rev B 28(6):3022–25. 1983;

23. Takeuchi N, Garzon IL. Ab initio Molecular-Dynamics Study of Structural, Dynamical, and Electronic Properties of Liquid Ge. Phys Rev B 50(12):8342–47. 1994;

24. Stich I, Car R, Parrinello M. Structural, Bonding, Dynamical, and Electronic Properties of Liquid Silicon: An ab initio Molecular-Dynamics Study. Phys Rev B 44(9):4262–74. 1991;

25. Ishmaru M, Munetho S, Motooka R. Generation of Amorphous Silicon Structures by Rapid Quenching: A Molecular-Dynamics Study. Phys Rev B 56(23):15133–38. 1997;

Article information Continued

This is an Open Access article distributed under the terms of the Creative Commons Attribution Non-Commercial License (http://creativecommons.org/licenses/by-nc/4.0/) which permits unrestricted non-commercial use, distribution, and reproduction in any medium, provided the original work is properly cited.