Investigation of Zr and Si diffusion behaviors during reactive diffusion – a molecular dynamics study

Abstract

Molecular dynamics simulation was used to investigate the diffusion behaviors of Zr and Si atoms during a reactive diffusion which produces Zr silicide. The simulation results were compared with those in Roy's experimental results. The profiles of mean square displacements (MSDs) of Zr and Si atoms at different temperatures were first used to evaluate the melting point above which the significant inter-diffusions of Zr and Si atom occur. The diffusion coefficients near the melting point were derived by the Einstein equation from MSD profiles. On the basis of diffusion coefficients at different temperatures, the diffusion barriers of Zr and Si atoms can be calculated by the Arrhenius equation. Compared to the corresponding experimental values, the predicted diffusion barriers at the Zr–Si interface were 23 times lower than the measured values in Roy's study. The main reason for this is that the Zr and Si atoms within the inter-diffusion region form different local ZrSi crystal alloys in the experiment, resulting in the lower diffusion coefficients and higher diffusion barriers found in the experimental observation.

---

1 Introduction

Zirconium silicide has attracted much attention because it exhibits good thermal stability² and excellent mechanical properties at high temperatures over 1273 K. Besides these excellent mechanical properties at high temperature, Zr silicide also has a relatively low resistivity (13–34 μΩ cm⁻³), indicating its potential as a material comprising the gate structure in MOSFET fabrication.³

In the fabrication of Zr silicide, the inter-diffusion of Zr and Si atoms is the major mechanism in the production process. Roy's study¹ used 1 mm thick Zr film to contact the Si(100) surface, of 0.7 mm in thickness. This Zr–Si alloy was then annealed at a temperature range between 1100–1250 °C for 16 hours. The Zr silicide was formed between the pure Zr and Si thin films by the inter-diffusion of Zr and Si atoms into the Si and Zr matrices. In Bertolino's study,⁴ bulk diffusion alloys were prepared by using multiple Zr and Si layers. The diffusion experiments were conducted at temperatures between 1173 K and 1623 K for 2–30 hours. By the inter-diffusion of Zr and Si atoms, several Zr silicide crystal alloys formed between the bulk Zr and Si matrices. In another Zr–Si nanostructure fabrication process, Lu used physical vapor deposition (PVD) process to deposit Zr on the Si(111) surface at substrate temperature ranging from 300 K to 1073 K.² The growth of Zr silicide is also attributed to the Zr and Si inter-diffusion mechanism.

Since atomic inter-diffusion is the major mechanism to produce Zr silicide, it is essential to further investigate the dynamic behaviors of the Zr and Si atoms during the Zr silicide growth at the atomic scale. It is relatively difficult to directly observe these atomic behaviors by the empirical approach, so the numerical method is an alternative way to obtain detailed insight into the diffusion mechanism for Zr silicide growth. Molecular dynamics (MD) simulations is a powerful tool to overcome the limitations of traditional empirical approaches and enable detailed observations of material behaviors on the atomic scale. For instance, Albe used MD to observe the shear band formation process for CuZr BMG under compression loading.⁵ They found plasticity in the glass layers is realized via pronounced, stable shear banding. Wu et al.⁶ studied the dynamic properties of a liquid Cu₈₀Si₂₀ alloy, and results show that activation energy of Si is greater than that of Cu, indicating that the bonding of Si and its surrounding atoms is stronger than that of Cu. Molecular dynamics is also widely used for investigating phenomena which take place rapidly. Shen used molecular dynamics to understand the solidification mechanism of aluminum and found that FCC structures are gradually surrounded by HCP structures at the beginning of solidification. Then the FCC nucleus grows, with the unstable HCP ordered around the FCC structures transforming to FCC. Finally, perfect FCC structures form at 400 K.⁷ Consequently, in the current study, MD simulation was used to simulate the inter-diffusion of Zr and Si of bulk stacked layers during the annealing process. The diffusion behaviors of Zr and Si atoms at the interface were presented and the diffusion barriers of Zr and Si atoms were also derived by the Arrhenius equation on the basis of diffusion coefficients for different system temperatures obtained by the MD simulation.

2 Simulation model

In order to model the Zr–Si alloy by molecular dynamics simulation, two potential functions, tight-binding⁸ and Tersoff,⁹ were used to describe the interactions between different atomic pairs. The tight-binding potential is responsible for modeling the interaction between the Zr–Zr pair, and the Si–Si and Zr–Si pairs are described by Tersoff potential. The parameters for these two potentials from our previous force-matching method¹⁰ fitting process are listed in Tables 1 and 2.

Table 1 Tight-binding potential parameters for Zr–Zr

Type	A (eV)	ζ (eV)	p	q	r0 (Å)
Zr–Zr	0.162	2.095	10.727	2.257	3.138

---

Table 2 Tersoff potential parameters for Si–Si and Si–Zr

Type	Si–Si	Zr–Si
A (eV)	7835.380	2251.660
B (eV)	45.087	175.073
λ	3.851	2.603
μ	1.079	1.474
β	0.429	0.468 × 10^−05
n	21.161	39.960
c	27 340.700	4061.980
d	119.344	3.252
h	−0.330	−0.062
R (Å)	2.783	3.216
S (Å)	2.986	3.562

---

The schematic diagram of the Zr–Si inter-diffusion model is shown in Fig. 1. To create the bulk diffusion blocks used in the experiment,¹ the many layers of Zr and Si thin films were contacted by the (001) surfaces. There are 56 000 Zr atoms and 74 088 Si atoms in our simulation model, and the MD simulation was conducted by LAMMPS.¹¹

Fig. 1 Stacked layers of Zr and Si blocks (56 000 Zr atoms, 74 088 Si atoms).

The inter-diffusion behaviors of Zr and Si atoms at the bulk Zr–Si interface were investigated by the MD temperature elevating process, starting from an initial temperature of 300 K. During this process, the isothermal–isobaric ensemble (NPT) was applied to the model shown in Fig. 1. The Nose–Hoover thermostat and barostat were adopted to adjust the system temperature and maintain a constant pressure of 1 atm during the temperature elevating simulation. The heating rate is 10 K per 5 ps and the average heating rate is about 2 × 10¹² K s⁻¹.

3 Results and discussion

The square displacement (SD) at different temperatures was used to indicate the specific temperature at which the inter-diffusions of Zr and Si atoms at the interface occur. The SD at temperature T is defined by a function of T as shown in eqn (1):¹²where r_i(T) represents the position of the atom i at system temperature T, and r_i(T₀) indicates the referenced position of the corresponding atom at the referenced system temperature T₀; and N represents the total atom number in the investigated system.

Fig. 2 shows the SD profiles of Zr and Si atoms at the interface from 300 K to 1500 K. It is obvious that the SD value of the Zr atom linearly increases with temperature increasing from 300 K to 650 K, while the Si SD profile remains almost constant throughout this temperature range, indicating the slow mobility of Si atoms at the interface at lower temperatures. Once the system temperature is higher than 650 K, abrupt increases in both the Zr and Si SD values can be clearly seen in Fig. 2, indicating the interfacial Zr and Si atoms possess enough kinetic energies to leave their equilibrium positions and conduct the inter-diffusion.

Fig. 2 SD as function of temperature for Zr (red curve) and Si (blue curve).

The morphologies of the Zr–Si interface at different temperatures are illustrated in Fig. 3(a)–(d). At 300 K, the mismatch between Zr and Si lattices causes the rearrangement of Zr and Si atoms at the interface. As the temperature continuously increases, the thickness of the inter-diffusion layer becomes thicker, growing from 3.92 to 14.28 Å as the temperature increases from 300 to 1500 K.

Fig. 3 Diffusion process of Zr–Si layers: (a) 300 K (b) 700 K (c) 1100 K (d) 1500 K.

Fig. 4 shows the diffusion distance from the reference interfacial surface for Zr and Si atoms versus system temperature. The average height of Zr and Si layer contact area at 300 K was used as the reference interfacial surface. The distances from the z coordinates of the bottommost Zr and the uppermost Si to the reference interfacial surface were defined as the diffusion distances. At temperatures under 650 K, the diffusion distances of both Zr and Si vary at constant values close to zero, indicating the inter-diffusion mechanism at this temperature range is not significant. At temperatures above 650 K, the diffusion distance of Si atoms displays a continuous and significant increase, while that of Zr atoms increases slightly with increasing temperature. Because the number density of Si is about 1.3 times higher than that of Zr, the Si atoms more easily diffuse into the Zr matrix. This is the main reason why the diffusion distance of Si is larger than that of Zr.

Fig. 4 Diffusion distance of (a) Zr and (b) Si during heating process.

The mean-square displacement (MSD) profiles at temperatures ranging from 300 to 1500 K for Zr, Si and interfacial atoms were used to investigate their dynamical properties. The MSD is defined by a function of time as shown in eqn (2):

where r_i(t) represents the position of atom i at delay time t, and r_i(t₀) indicates the reference position of the corresponding atom at reference time t₀; N represents the total atom number of the investigated system. The brackets are interpreted as average over time origins and the total atom number of the system. From Fig. 5(a)–(c), it is clear that the slopes of the MSD profile are generally larger at higher temperatures for all cases. It should be noted that the MSD profiles for 300 K remain constant, indicating that the Zr and Si atoms undergo thermal vibrations around their equilibrium positions. It should be noted that although the SD and MSD were calculated on the basis of the coordinates of atoms in the bulk region and at the Zr–Si interfaces of our simulation models, the variations of these two functions are very sensitive to the changes of atomic coordinates. Since the atoms in the bulk region only conduct the thermal vibrations at their equilibrium positions and almost have no contribution to the SD and MSD values, the variations of these two functions can be attributed to the atomic diffusion at the interface.

Fig. 5 The mean-square displacement plots (MSD) of Zr–Si interlayer: (a) total (b) Zr (c) Si.

It is known that the MSD profile is linear to the delay time in the long-time limit, and thus the diffusion coefficients of Zr and Si atoms at the interface can be derived from the slopes of MSD profiles after a longer delay time by the Einstein equation:¹³

where D is the self-diffusion coefficient and N is the number of atoms. According to the MSD profiles shown in Fig. 5 for different temperatures and the Einstein equation, the total, and average Zr and Si diffusion coefficients at different temperature can be obtained. Further, in Roy's study, they used the Arrhenius equation to derive the diffusion barrier on the basis of the calculated diffusion coefficient near the melting points of the crystal ZrSi₂.¹ The formula of the Arrhenius equation for describing the diffusion coefficient at different temperatures is:where Q is the activation energy, T is the temperature, D₀ is the preexponential factor, and R is the Boltzmann constant. To calculate the diffusion barrier, the profiles of ln(D) versus 1/T for total, Zr, and Si atoms are shown in Fig. 6. It can be inferred that the diffusion coefficients of both Zr and Si atoms significantly increase with the increasing temperature when the system temperature exceeds the melting temperature. The values of these diffusion coefficients range between 10⁻¹⁰ and 10⁻¹¹, which are two orders larger than those found in the experimental measurement for the bulk inter-diffusion region.¹ In Roy's study, they investigated the atomic diffusivity within the inter-diffusion region at which a certain amount of local denser C49 crystal structure forms. These denser local structures lead to the considerable cage effect, which causes the backflow of atoms when an atom interacts with atoms of local C49 crystal structure, resulting in it jumping back to its initial position and lowering diffusivity.

Fig. 6 Diffusion coefficient for the Zr–Si layer.

Note that the ln(D) profile is proportional to the inverse of the temperature, and the diffusion barriers of Zr and Si atoms at the interface can be derived by the slopes of the ln(D) profiles. The D₀ values and the diffusion barriers are listed in Table 3. The diffusion barriers of Zr and Si are 12.6 and 14.9 kJ mol⁻¹, values which are about 23 times smaller than those obtained in Roy's measurement.¹ The main reason is that the Zr and Si atoms within the inter-diffusion region form different local ZrSi crystal alloys, resulting in the lower diffusion coefficients and higher diffusion barriers in the experimental observation.

Table 3 The estimated pre-exponential factor (D₀) and activation energy (Q) of the Zr–Si layer

Temperature interval	Type	D0 (m^2 s^−1)	Q (kJ mol^−1)
700–1500 K	Zr	8.5 × 10^−11	12.6
	Si	1.0 × 10^−10	14.9
	Total	9.6 × 10^−11	13.5
1000–500 K	Zr	4.4 × 10^−10	31.0
	Si	3.1 × 10^−10	28.1
	Total	4.1 × 10^−10	30.0

---

4 Conclusion

MD simulation was used to investigate the inter-diffusion behaviors of Zr and Si atoms at the Zr–Si interface during a reactive diffusion of Zr/Si bulk matrices. The temperature at which the inter-diffusion at the Zr–Si interface begins was first predicted by the MSD profiles. From the simulation result, the inter-diffusion begins at 700 K, which is about 400–500 K lower than the annealing temperatures used in the related experiments for bulk diffusion measurement.² The diffusion coefficients near the melting point derived by the Einstein equation from MSD profiles were two orders higher than the experimental measurement, while the predicted diffusion barriers at the Zr–Si interface were 23 times lower than the values measured in Roy's study.¹ The main reason is that the Zr and Si atoms within the inter-diffusion region form different local ZrSi crystal alloys in the experiment, resulting in the lower diffusion coefficients and higher diffusion barriers in the experimental observation.

References

1. S. Roy and A. Paul, Mater. Chem. Phys., 2014, 143, 1309–1314.
2. M. Le Flem, J. Canel and S. Urvoy, J. Alloys Compd., 2008, 465, 269–273.
3. B. Lu, W. Xiao, S. S. Kushvaha, P. He and X.-S. Wang, J. Phys.: Condens. Matter, 2008, 20, 225015.
4. F. Maglia, U. Anselmi-Tamburini, C. Milanese, N. Bertolino and Z. A. Munir, J. Alloys Compd., 2001, 319, 108–118.
5. D. Sopu, Y. Ritter, H. Gleiter and K. Albe, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 83, 134208.
6. S. Wu, M. Kramer, X. Fang, S. Wang, C. Wang, K. Ho, Z. Ding and L. Chen, Phys. Rev. B: Condens. Matter Mater. Phys., 2011, 84, 134208.
7. B. Shen, C. Y. Liu, Y. Jia, G. Q. Yue, F. S. Ke, H. B. Zhao, L. Y. Chen, S. Y. Wang, C. Z. Wang and K. M. Ho, J. Non-Cryst. Solids, 2014, 383, 13–20.
8. F. Cleri and V. Rosato, Phys. Rev. B: Condens. Matter Mater. Phys., 1993, 48, 22–33.
9. J. Tersoff, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 39, 5566–5568.
10. F. Ercolessi and J. B. Adams, Europhys. Lett., 1994, 26, 583–588.
11. S. Plimpton, J. Comput. Phys., 1995, 117, 1–19.
12. C. Dietrich, L. A. Bagatolli, Z. N. Volovyk, N. L. Thompson, M. Levi, K. Jacobson and E. Gratton, Biophys. J., 2001, 80, 1417–1428.
13. M. Meunier, J. Chem. Phys., 2005, 123, 134906–134907.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c4ra16962j

---