Abstract

Tight-binding molecular-dynamics (TBMDs) simulations are performed to study atomic and electronic structures during high-temperature consolidation processes of nanocrystalline silicon carbide under external pressure. We employ a linear-scaling method (the Fermi-operator expansion method) with a scalable parallel algorithm for efficient calculations of the long time-scale phenomena. The results show that microscopic processes of the consolidation depend strongly on initial orientations of the nanocrystals. It is observed that an orientational rearrangement of the nanocrystals initially misaligned is induced by an instantaneous shearing force between nanocrystals, whereas the aligned system undergoes densification without shearing. Analysis on an effective-charge distribution and an average bond-order distribution reveals electronic-structure evolutions during these processes.

1. Introduction

Nanocrystalline ceramics have been one of the most fascinating research subjects because of their promising properties for industrial applications, such as larger fracture toughness, higher sinterability than conventional ceramics, and superplastic behavior at elevated temperatures [1, 2]. Especially, Silicon-Carbide-based micro/nanostructures have been also viewed as a building block for new optoelectronic devices and nuclear-radiation detectors [3, 4]. Polycrystalline SiC has been known to have limitation of its usefulness due to the low sinterability [5]. As in other nanocrystalline ceramics, reduction of the size of SiC powder leads to lower sintering temperature and higher degree of densification. Ohyanagi et al. have reported [6] that nanocrystalline SiC (nc-SiC) without additives has an onset temperature of sintering, which is significantly lower than that of coarse-grained powders. Such a high sinterability of pure nc-SiC and the better radiation resistance possessed by highly densified nc-SiC have also been indicated by classical molecular-dynamics simulations [7, 8]. Further understanding of consolidation processes of nc-SiC from an electronic-structure level leads to an essential improvement in designing and optimizing materials performance especially in optoelectronic applications. The density-functional-theory (DFT-) based electronic-structure calculation is known to have an advantage in its accuracy and transferability compared with other semiempirical methods including the tight-binding approach and, thus, has been widely adopted for various types of materials simulation. However, due to rather heavy demands for computer resources in the DFT calculations, high-temperature and/or long time-scale phenomena such as consolidation processes of nanocrystals are far beyond being approached by the method. Semi-empirical tight-binding molecular dynamics method [9] with highly optimized parallel algorithms is a key to approaching alternatively to this demand.

In this paper, we report on parallel tight-binding molecular-dynamics (TBMDs) simulations of consolidation of nanocrystalline SiC with a diameter 2.4 nm. The Fermi-operator expansion method (FOEM) [10] has been employed to calculate efficiently the electronic part of the energy and forces in the TBMD simulations. We have also implemented a constant-pressure MD algorithm to virtually simulate the consolidation processes and it has been run on our parallel PC cluster. Using the parallel TBMD method, we investigate characteristics in the initial stage of consolidation, that is, interparticle interactions just before and after the initial contact between SiC nanocrystallites, by comparing two different mutual orientations of the crystallites. Analysis on an effective ionic charge and a bond-order distribution is used to characterize evolutions of local electronic states during the consolidation.

2. Methodology

2.1. Tight-Binding Total Energy Model of Silicon Carbide

In the tight-binding (TB) model [9], the total energy of a system consisting of N atoms is defined as where stands for the kinetic energy of ions, is the band-structure energy of valence electrons, and represents the repulsive term that takes into account the core-core interactions between ions and neglected contributions in to the true electronic energy, such as a correction for double counting of electron-electron interactions. The repulsive energy is modeled by the sum of short-range 2-body interaction. The band-structure energy is calculated by diagonalizing the effective one-electron Hamiltonian matrix. Each off-diagonal element of the TB Hamiltonian involves interactions between valence electrons within the two-center hopping approximation. The electronic contribution to interatomic forces, that is, derivatives of with respect to the atomic coordinates, can be obtained through the Helmann-Feynman (HF) theorem.

Parameterization of matrix elements in the semi-empirical TB Hamiltonian is an essential ingredient for TBMD simulations. For SiC systems, we have chosen Mercer’s parameterization of the TB Hamiltonian [11], which is based on an orthogonal basis. It also includes environment-dependent contributions to the onsite energies through intra-atomic terms. These terms give important contributions to variation of the electron transfer between Si and C atoms, especially for inhomogeneous systems, although it does not include explicitly the Coulomb interaction between ions. The original parameterization by Mercer, however, leads to some discrepancies with experimental data such as the lattice constant and the interfacial energies. Minor modifications of the parameters, have therefore, been made so that magnitudes of these discrepancies are reduced [12].

2.2. Linear Scaling Parallel Algorithm for TBMD with NPT Ensemble

In ordinary electronic-structure calculations, the Schrödinger equation with a one-electron Hamiltonian is solved through numerical diagonalization of the Hamiltonian matrix which requires the computational cost proportional to , where is the rank for valence band in the matrix. Although the semi-empirical TB method reduces significantly the cost for constructing the Hamiltonian by the a priori parameterization of each matrix element as functions of atomic positions, this N dependence of the CPU time becomes overwhelming when the simulation deals over one thousand of atoms and/or for more than thousands of time steps, which are essential for realistic simulations of nanostructured materials. Several new methods have already been proposed to resolve this problem [13]. Among them, we have adopted the Fermi-operator expansion method, proposed originally by Goedecker and Colombo [10]. This algorithm is based on a moment expansion of a pseudodensity matrix (the Fermi operator) of the TB Hamiltonian by Chebyshev polynomials. Combining an appropriate truncation at a finite order of the expansion and a truncation in the multiplication between each element of the matrix at a physical cut-off distance, the computational complexity in evaluating the band-structure energy and the HF forces is reduced to that proportional to (linear scaling). Moreover, a second-order approximation in updating the Chebyshev polynomials is used in order to gain further efficiency. This technique accelerates the calculations about twice as fast as those with ordinary recursion algorithm for calculating the polynomials.

The truncation at the physical distance between atomic sites can make the data structures localized, hence a scalable parallelization can easily be achieved in the present algorithm. We employed a spatial domain-decomposition technique for calculating matrix-matrix multiplication. In this algorithm, data transmission between neighboring nodes is performed only for near-boundary atoms and matrix elements because of the spatially localized nature of the density matrix. Internode data transmission is performed using the standard MPI library. We have been developing and testing the fortran code based on this algorithm on various computational platforms, including massively (over 1,000 CPU) parallel architectures [13]. Also, this algorithm has been used successfully for other studies on nanostructured SiC systems [12, 14]. In the present study, we have performed the parallel TBMD simulations on our 160 cores parallel PC cluster.

In the TBMD simulations, we used a velocity-Verlet algorithm [15], with the time step of 0.2 fs (femtosecond). Temperature of the system and size/shape of the MD cell were controlled using the explicit reversible integrator algorithm [16] so as to realize the system to be an NPT ensemble at a thermal equilibrium. In the calculation of internal stresses for NPT dynamics, we employed a symmetric decomposition method of the TB Hamiltonian [17]. The criterion of convergence in the HF force calculation by the Chebyshev expansion was set with 1.00−2 eV/Å per atom, and the pseudoelectronic temperature was chosen to be 0.1 eV. These values gave a reasonable accuracy and numerical stability in the present simulations at a high temperature (1500 K).

3. Results

3.1. Preparation of SiC Nanocrystals

A spherical nanocrystal of 3C-SiC was initially prepared by removing atoms within a spherical region of radius about 12 Å in the prefect crystal. In addition, some atoms less bonded (i.e., having only one or two nearest neighbors within a cutoff distance 2.3 Å) on the surface of the cutout sphere were removed so that the overall stoichiometry in the nanocrystal was maintained. The nanocrystal, thus, prepared contains 504 Si atoms and 504 C atoms. Four copies of the nanocrystal were then placed at the positions of face-centered cubic (FCC) cell with a periodicity 42 Å for three directions, respectively. The FCC configuration was chosen because we assumed that before consolidation the system had been well thermalized so that it could be modeled as an ideal close-packed mutual positions. In setting the crystallographic orientations of each nanocrystal, we prepared two sets of the system; one with four nanocrystals aligned mutually so that all the axes of the crystallites coincide with those of the cubic periodic cell, and the other with randomly oriented four nanocrystals. Figure 1 illustrates the initial configuration of the system with (a) aligned and (b) misaligned nanocrystals. Each shows one crystallite at center plus three crystallites and their images appeared due to the periodic-boundary condition in the simulation cell. The central nanocrystal is thus surrounded virtually by twelve nanocrystals. As shown in the figure, the initial cell size (42 Å) was given so that the nanocrystals do not touch each other and, hence, can be regarded initially as independent systems.

Subsequently, we ran TBMD simulations with an NVT ensemble in order to thermalize the systems before applying an external pressure. Starting at the room temperature, the system was gradually heated to 1500 K and was then thermalized for 5000 MD steps (1 ps). During this heating/thermalization procedure, the centers-of-mass of the four nanocrystals were kept at the original FCC positions so that the nanocrystals do not interact with each other. Also, angular momentum of each nanocrystal was kept zero during the simulation in order to eliminate artificial rotations of the crystallites.

3.2. Consolidation Processes: Evolution of Internal Stress

After the initial thermalization processes were completed, we switched to constant-pressure MD runs for consolidation simulation. The applied external pressure is 1 GPa and the temperature is kept at 1500 K. The external pressure 1 GPa is chosen so that fluctuation of the internal stresses is suppressed by the pressure during the densification. Also, in the previous reports [6, 7, 14], an onset of sintering has been found around the 1500–1700 K. Thus, we focus on monitoring the time evolution at the temperature of 1500 K. Figures 2 and 3 depict evolutions of the system volume and hydrostatic and shear component of the internal stress during the MD runs for the aligned (Figure 2) and the misaligned (Figure 3) systems. Starting from the original cell size, the system volumes in both cases are decreased linearly with time. At around 0.5 ps in Figure 2, the hydrostatic component of the internal stress is shifted from zero to a negative value in average, indicating the system being under tension. Since the surface atoms begins to interact with others on neighboring crystallites at the time step, this behavior of the internal stress can be attributed to an attractive force exerted between nanocrystals in the aligned case. At around 2 ps, the aligned system is changed to be under compression and the hydrostatic stress of the system is then increased to 1 GPa so as to balance with the external pressure.

On the other hand, figure 3 shows that the hydrostatic pressure is essentially zero upto about 3 ps in the case of misaligned nanocrystals. Since the average surface distances between neighboring crystallites are approximately the same in both cases of Figures 2 and 3 at around the same elapse time, this difference indicates that an effective intercluster interaction depends sensitively on the mutual orientations of the nanocrystals. At around 3 ps in Figure 3, the hydrostatic pressure in the system with the misaligned nanocrystals begins to increase until the stress reaches around the equilibrium value of 1 GPa. Also, the shear component of the internal stress in the misaligned case exhibits a pulse-like increase at around 3.5 ps (highlighted by a circle in the figure). At this time step, some of the nanocrystal pairs begin to touch with each other and rotate their mutual orientations as if they spontaneously optimize their grain-boundary structure. This observation is discussed again in the next section.

Sintering is, in general, driven by free-energy reduction, and the free energy one should consider in an NPT ensemble is the Gibbs free energy. It is not easy to estimate the absolute value of the free energy directly from MD simulations, here we only consider the enthalpy reduction, assuming the entropy contribution to the energy is relatively small. In the present NPT-MD simulations, the enthalpy has also been monitored. In both the aligned and the misaligned cases, it is found that the decrease of the enthalpy, which is about 10% reduction, after the first contact between nanocrystals is attributed to the decrease of both the internal energy and the pressure-volume term (PV), but the main contribution comes from that of the internal energy. Thus, it implies that the free-energy change which drives the initial stage of consolidation of SiC nanocrystals is mainly due to the surface-energy reduction by the formation of inter-crystallite bonds. We will see in the next section that this interpretation is also consistent with the analysis on time evolutions of effective-charge distributions.

3.3. Consolidation Processes: Effective Ionic Charge and Bond Order

In order to characterize evolution of local electronic states during the initial stage of consolidation of SiC nanocrystals, we perform the electron-population analysis based on the Mulliken’s scheme [12, 18]. An effective ionic charge for each atom is then calculated by subtracting the formal valence (4) of Si and C atom from the calculated electron population at th atomic site calculated by where is the orbital index, and represents the diagonal matrix element of the Fermi-operator as a function of the energy level , obtained approximately by the Chebyshev polynomials explained in Chap. 2. In our TB model, these values for cations and anions in the perfect crystal are ±1.83, respectively.

In Figure 4, we plot a time evolution of spatial distribution of the atoms in the misaligned nanocrystals and their effective charge on a cross section at the (111) plane of the MD cell. At the initial step, depicted in Figure 4(a), the effective charge varies from −2 to −1 on anions and from 1 to 2 on cations, and larger variations can be found mainly on the surface of each nanocrystal due to the existence of dangling bonds on the surfaces. After starting the densification, the system is shrinking and at 3.5 ps nanocrystal pairs (A-C and B-C pair in Figure 4) touch with each other forming grain boundaries, as shown in Figure 4(b). It should be noted that the crystallographic orientations of the central (C) and the touched (A, and B) crystallites are changed slightly from the original ones. This indicates that the shearing force exerted between these nanocrystals results in a spontaneous adjustment of mutual orientations between them, as discussed in the previous section. After about 5 ps elapsed, all the nanocrystals begin to form grain boundaries, as shown in Figure 4(c). At this stage, the effective-charge distribution becomes more uniform compared to that in Figures 4(a) and 4(b) due to a significant reduction of dangling bonds on the interfaces. Thus the system evolves to the next consolidation stage, that is, the neck growth via diffusion processes on the interface [14].

We also characterize the local electronic states, which are related particularly to covalent bonding, via an average bond order [18], calculated as where “n.n.” stands for the nearest neighbor atoms within 2.3 Å around the ith atomic site with being the number of its n.n. atoms. This quantity is a measure of strength of covalent bonding in the system, whereas the variation of the electron population given by (2), which is an energy integral of the local density of states, represents the electron transfer between cation and anion.

Figure 5 shows time evolution of the average bond order during the consolidation. The color indicates the value of scaled between 0 (blue) and 1 (white). Unlike the inhomogeneity shown in Figure 4, the bond-order distribution is rather uniform in Figure 5(a) and does not change much even when the nanocrystal pairs begin to touch with each other, as in Figure 5(b). It is, however, remarkable in Figure 5(c) that the bond order decreases rapidly once after the process is evolved to the stage where the internal stress balances with the external pressure. This phenomenon may be indicating a delocalization of electrons at the bonds as an effect of the compressive stress.

4. Conclusion

Using a linear-scaling parallel algorithm, we have performed large-scale/long time-scale TBMD simulations with NPT ensembles for analyses of the initial stage of consolidation of nanocrystalline SiC with a diameter of 2.4 nm. It is revealed that an effective intercluster interaction is strongly dependent on the mutual orientations of nanocrystals at high temperature. The densification under an external pressure at 1500 K proceeds without large change of crystallographic orientation between aligned nanocrystals, whereas initially misaligned nanocrystals exhibits an instantaneous shearing motion of the crystallites on the contact of their surfaces prior to the grain-boundary formation. The time evolution of local electronic structures in the consolidation process has been studied by monitoring the effective ionic charge and the average bond order on each atomic site. It is found that the inhomogeneity of the charge distribution at the surfaces of crystallites is observed to be reduced during the consolidation due to the bond formations and reconstruction at the interfaces. The evolution of the bond-order distribution indicates a tendency of electron delocalization on covalent bonds due to the high internal stress induced by the external pressure.

Acknowledgment

This work was supported by Grant-in-Aid for Scientific Research on Priority Areas “Nanomaterials Science for Atomic Scale-Modification 474” from MEXT of Japan.