Evolution and Intersection of Extended Defects and Stacking Faults in 3C‐SiC Layers on Si (001) Substrates by Molecular Dynamics Simulations: The Forest Dislocation Case

The crossing of stacking faults (SFs) (namely dislocation forests) in 3C‐SiC is experimentally probed to cause detrimental effects on the electronic properties of the layer. By means of classical molecular dynamics simulations, for the first time the evolution of two crossing partial dislocation loops is shown, revealing the mechanism that leads to the formation of a complex defect at their intersection. The obtained atomistic structure is then further refined by ab initio calculations, revealing that this defect does not cause defect states within the gap of the pristine material. The case of the intersection of two double‐dislocation loops, that is, involving double SFs, in the hypothesis that the evolution follows the same mechanism found for single‐dislocation loops is also investigated. Contrary to the single‐plane intersection, the simulations reveal that the junction formed by double‐dislocation loops does cause intragap states, thus suggesting detrimental effects on the electronic properties of the 3C‐SiC layers.

(called dislocation forests) can have a detrimental impact on the electronic properties of the 3C-SiC epilayer, [28] so that they could be considered as killer defects in 3C-SiC together with the dislocation complexes. [36] Although dislocation forests are naturally present in 3C-SiC layers on Si substrates because of the elevated SF densities, a specific theoretical analysis on these defects is still lacking, neither on their formation mechanism, nor on their structure, and nor on their electronic properties. Our work contributes to filling this gap.
In this article, a multiscale approach combining classical molecular dynamics (MD) simulations and ab initio calculations has been used to study the crossing of SFs in 3C-SiC epitaxial layers. The two methods resulted to be complementary in the description of the formation of defective complexes and the identification of the corresponding impact on the quality of the material in terms of electronic performances. We show that the imaging of two crossing SFs has to be interpreted as the result of the intersection of two adjacent partial dislocation loops that open along the same 〈110〉 direction in two different {111} planes. Indeed, under tensile strain, they can intersect and continue their motion, giving rise to a new line defect at their crossing point. Actually, the two SFs intersect and continue to expand, while the partial dislocations bounding them remain first pinned at the crossing point, change their line directions, and give rise to the formation of the new defect after depinning. The defect structure corresponds to the deposition of four 90 dislocation segments above and below the two SFs. Ab initio calculations demonstrate that defect states in the middle of the bandgap appear if the crossing SFs consist of multiple planes, as in the case of the dislocation complexes terminating multipledislocation loops.

Computational Methods
We conducted classical MD simulations using the software large-scale atomic/molecular massively parallel simulator (LAMMPS) [40] and Vashishta interatomic potential. [41] Vashishta potential has some advantages over other empirical potentials for silicon carbide [42] : this potential considers the atom interactions beyond the first-neighbor shells; therefore, different total energies for the hexagonal and the cubic configurations are provided, thus enabling for an estimate of nonzero SF energies. The simulations are performed in the canonical (i.e., the number of particles, the volume, and the temperature are constant, NVT) ensemble, using a Nose-Hoover thermostat. [43] The atomic trajectories were analyzed using the open visualization tool (OVITO) software, [44] which identifies extended defects in virtual crystals, highlighting SFs and dislocation lines along with the associated Burgers vector. Orthogonal simulation cells, oriented along the directions determined by the vectors u ¼ a=2½110, v ¼ a=2½110, and w ¼ a½100, were used, and periodic boundary conditions (PBCs) were applied in the u and v directions. The bottom three layers in the simulation cell were fixed at the atomic position of the bulk, to model the effect of the Si substrate below the 3C-SiC layer. A large box (38.8 Â 38.8 Â 26.0 nm 3 corresponding to 3'456'000 atoms) was needed to analyze the evolution and interaction of multiple dislocations.
Dislocation loops were inserted in the simulation cell by shifting each atom according to the displacement field vectors calculated by the dislocation theory in the framework of the linear elasticity theory. [45] The obtained configurations were geometrically optimized using a conjugated gradient minimization algorithm. For an arbitrary dislocation segment, the field components can be calculated by rotating the reference system and decomposing the dislocation Burgers vector into a screw (b k ) and an edge component (b ⊥ ). At this, b z axis has to be oriented along the dislocation line ξ, b x axis perpendicular to ξ in the glide plane, and b y axis oriented following the right-hand rule, respectively. In this case, the components of the displacement field vectors turn to be the following: where ν is the Poisson ratio of 3C-SiC. [46] To insert a hexagonal loop, one has to divide the space in six regions. In each region, the displacement field due to the presence of a dislocation segment can be calculated by rotating the reference frame according to the segment direction and glide plane and exploiting Equation (1)- (3). Then, the total displacement can be obtained as a sum of the displacement fields of the six segments. Notice that the displacement field components (1)-(3) have to be corrected to have the stress normal to the free surface equal to zero after system equilibration. This has been done by adding an image dislocation loop in the vacuum above the free surface. [45] The annealing of the systems was carried out at the thermostat temperature of 1400 K, corresponding to the typical experimental temperature for 3C-SiC deposition. The integration time step was 1 fs. Moreover, we scaled the simulation box dimensions by a factor calculated using the Nose-Hoover barostat, to avoid the presence of residual strain due to the thermal dilation of the layer. The electronic structure calculations are based on density functional theory (DFT) and exploit the recently developed nonempirical strongly constrained and appropriately normed (SCAN) meta-generalized-gradient approximation (GGA) functionals, [47] as implemented in the Quantum Espresso package. [48] This allows to improve the prediction of the bulk SiC bandgap (%1.9 eV) as compared with the value of about 1.35 eV obtained by GGA. [36] The atomistic models obtained by the classical MD calculations were used as inputs for the DFT calculations, after scaling the supercells' size to match the lattice parameters of bulk 3C-SiC. [49] These were obtained by using GGA exchangecorrelation potentials [50] and accounting the van der Waals interaction by the semiempirical method of Grimme (DFT-D2). [51] Within the projector-augmented wave (PAW) method, [52] a planewave cutoff energy of 80 Ry was used for ionic relaxations. The cutoff for the atomic forces was set to 10 À4 Ry/Bohr, exploiting a conjugate gradient method and a 1 Â 4 Â 1 Monkhorst-Pack grid for sampling the Brillouin of the supercells.
The supercells have been constructed by replicating 16 Â 2 Â 4 times the unitary slab of the zinc-blende structure with the x axis aligned along the ð11 2Þ crystallographic direction, y along the ð110Þ direction, and the z along the ð111Þ direction, resulting in 768 atoms. The final structure has been fit into an orthogonal 4.26 Â 0.62 Â 3.2 nm 3 box. Then, multiple straight partial dislocations were inserted by shifting the atoms in the cell and considering the displacement fields (Equation (1)

Results and Discussion
Looking at the interaction between SFs, several structural configurations can be probed. [27,37,38] In particular, it has been suggested that forest dislocations lead to deleterious effects in the electronic properties of the layer, [27] but a study of the atomistic structure of these defects and its relative electronic properties is still lacking.
In this section, we present novel results on the dislocation forest divided in two sections: in the first one, we show MD simulations of the evolution of crossing loops, leading to the formation of a dislocation forest, and in the second, atomistic details of the complex defect thus formed are presented, and their relative electronic states obtained by ab initio calculations are discussed.

Evolution of Crossing SFs and Formation of Dislocation Forest
Forest dislocations are complex defects that form at the intersection of the two expanding dislocation loops whose threading arms cross each other. Due to the mismatch with the substrate, the cubic SiC layer is in a tensile state in the first deposition stage. In this case, 90 partial dislocations parallel to the interface are known to take part in strain relaxation, [31] creating SFs within the layer. Notice that partial dislocations are formed by the dissociation of perfect dislocations with 1 2 h110i Burgers vector. The dissociation lowers the dislocation energy, replacing the perfect dislocation with a pair of glissile Shockley partial dislocations, with Burgers vectors 1 6 h112i and an SF between them. Such a process is favored when the dislocations are placed in narrowly spaced {111} planes (i.e., in the glide set). Indeed for dislocation in the shuffle set (in widely spaced {111} planes), the shear needed for dissociation produces SF of very high energy. [53] Following the procedure described in the previous section, we construct in our cell two dislocation loops in the glide set with , lying on ð111Þ and ð111Þ planes, respectively. Loops are placed consecutively, one close to the other with the segment parallel to the interface along the same ½110 direction. Being the loop constituted by partial dislocation segments in its internal part, SF is present. Such an initial configuration is shown in (a) of Figure 1. Loops have initially a hexagonal shape due to the fact that the preferred dislocation lines in the zinc-blende structure are 〈110〉 in type, so that six favored dislocation segments are present on each {111} plane. Consequently, dislocation segments that are not parallel to the interface turn out to be 30 partial dislocations. Immediately the evolution computed by MD changes the loop shape, making them more rounded (see also the study by Barbisan [37] ). We applied a tensile strain of 4% in the cell. The value of such strain is set to be sufficient to observe the loops' evolution in a reasonable computational time but not too high so as to prevent disorder in the cell (as shown in ref. [37]). Indeed, the main part of the relaxation of the initial 20% misfit strain in the layer is known to be done by an array of Lomer dislocations formed in the first deposited layers as a structural reconstruction directly at the interface. [32] So, it is impossible to correctly estimate the real strain acting in the loop opening being dependent on the local formation of Lomer dislocations. The two loops expand, as shown in (b), driving two threading arms one close to the other at a simulation time of 15 ps. When they reach, a junction is developed, as shown in Figure 1d, by an enlarged image of the crossing configuration. Upon further evolution, we can observe that the dislocations do not cross the SF plane of the encountered loop. Instead, each threading dislocation remains pinned at the junction point, changes locally the dislocation line direction, and deposits two segments parallel to the interface, as shown in Figure 1e, where the threading arms are colored differently (light blue and red) to clarify the mechanism. Then, threading dislocations further move apart via a depinning process and leave a complex line defect along the junction. The visualization program OVITO cannot resolve it; therefore, a white cylinder is drawn in the simulation. However, the resolved defect structure can be easily analyzed by cutting a slice along the intersection line (white cylinder) and looking along those lines shown in Figure 2a. Clearly, forming the junction, the threading arms change their direction, also changing the angle between Burgers vector and the dislocation lines, that is, the dislocation type. Each segment becomes a 90 dislocation, as all the segments parallel to the interface. As a consequence, the new complex defect developed, namely, the forest dislocations, turns out to be formed by four adjacent 90 dislocations and two SFs, as shown in Figure 2b. The four composing 90 dislocation cores are highlighted in red in (c-f ) of Figure 2.
In the study by Nagasawa et al., [28] it is clearly shown that the majority of SFs near the interface of 3C-SiC layers on Si(001) substrates turns out to be multiple-plane SFs. In the study by Scalise et al., [36] a mechanism for the formation of double and triple SFs via the gliding of adjacent Shockley partials is proposed. It is validated by MD simulations of the evolution of only the segments parallel to the interface. It is almost impossible to follow the evolution of the entire multiple loops in three dimensions in reasonable computational time. However, it can be supposed that also multiple SFs can intersect each other, www.advancedsciencenews.com www.pss-b.com as illustrated earlier (Figure 1) for the single loop. In particular, double SFs are composed of adjacent 30 and 90 partial loops, forming a defect with the 30 total Burgers vector, and therefore it can further evolve. In contrast, triple SFs have a null total Burgers vector, so they do not move and their intersection is statistically less probable. [36] In analogy with the mechanism found by simulation for single loops, we constructed in a periodic cell the complex defect that should appear at the intersection of double SFs. It is shown in Figure 3a. Notice that the cell is tilted to recover full periodicity. Eight partial dislocation segments parallel to the interface are inserted in the cell, 30 and 90 in type. The procedure was Figure 2. Analysis of the junction structure. a) View of the cut plane considered (light blue atoms are in the perfect cubic lattice, orange atoms belong to SFs, in green are the dislocation lines); b) the projected view in 〈110〉 direction of the junction structure (white and colored atoms, except for orange and light blue, do not belong to any recognized crystal structure); c-f ) zoomed images of the junction core and atoms highlighted in red can be associated with the four 90 dislocation segment cores. www.advancedsciencenews.com www.pss-b.com illustrated by showing the four different double-partial dislocations inserted in the four panels on the right. In each one of the panels (b-d), we labeled the ð111Þ and ð111Þ planes in terms of α, A, β, B, γ, and C monolayers, following the Ramsdell notation. [54] Such a notation is useful to analyze the stacking order change in the formation process of the junction structure shown in Figure 3. Red atoms belong to 90 partial dislocation cores, while green and pink atoms to 30 partial dislocation and are antiparallel in the screw component. Each double dislocation shown in each panel corresponds to one of the two segments deposited in the junction by one of the crossing threading arms of the double loop after depinning. (b) and (c) show the segments relative to the threading arms traveling in a ð111Þ plane (equivalent to the red one in Figure 2) and (d) and (e) to the threadings traveling in the ð111Þ plane (equivalent to the light-blue one in Figure 2). The former threadings change the conventional cycling stacking order ( : : : αAβBγC : : : ) into the anticyclic one ( : : : γCβBαA : : : ) for four consecutive ð111Þ monolayers ( (b) and (c)), and the latter threading does the same for four consecutive ð111Þ monolayers ((d) and (e)). The cell was then relaxed by a conjugated gradient procedure using Vashista potential. The realistic atomistic structures of forest dislocations are established; we used them in a superior ab initio calculation to further relax the structure and analyze their electronic properties.

Electronic Properties of Dislocation Forests
The electronic density of states (DOS) calculated for an intersection between the single loop is shown in Figure 4. The simulation cell has been built periodically, inserting the structure of the junction obtained by classical MD simulations, developed by two crossing double SFs shown in Figure 3, to facilitate the DFT calculations. An ionic relaxation has been performed by DFT. Then, the DOS was calculated and integrated inside 20 slices of the cell and plotted as a function of the energy (with respect to the calculated Fermi energy) and its intensity, whose color scale is also shown in Figure 4. The integrated DOS in each slice is plotted along one of the directions of the cell, and specifically the horizontal one is shown in the corresponding atomistic model of Figure 4. This approach may help to identify possible defective electronic states in the bandgap of the material and localize them within a real space region in the crystal cell model. The calculated DOS shows quite clearly that the single SF does not introduce any relevant electronic state in the bandgap of the 3C-SiC. In fact, far from the central region plotted in Figure 4,  Figure 1e). On each side of the dislocation complex that binds the SF, the ð111Þ and ð111Þ planes are labeled using the Ramsdell notation (e.g., : : : αAβBγC : : : ). the bandgap is very close to our predicted gap for bulk 3C-SiC (of about 1.9 eV), in agreement with our previous findings. [36] In the central region of the cell modeled, corresponding to the intersection zone between the two SFs, a few states appear, close to the valence band edge. However, these states may be identified as very shallow states, similar to the states found for the dislocation complex terminating the single SFs. [36] Thus, they should not be very detrimental for the devices. The same approach has been used to investigate the electronic properties of the intersection between two double SFs. Contrary to the single-loop case, the DOS of the intersection between double SFs, illustrated in Figure 5, shows several defects states with the bandgap of the pristine material. These intragap states are in energy just a few hundreds of meV one from each other and spatially very localized at the intersection point. Thus, their impact on the electronic transport and particularly on the leakage current into the device may be marked, as also already shown for the dislocation complexes terminating the double SFs. [36] A deeper analysis on the physical origin of these intragap states needs a further and specific study, which is out of the scope of this article, but few important points may be highlighted here. The two states above the Fermi level and well below the conduction band edge are induced by the dangling bonds of C atoms close to the intersection point. In our previous work, [36] we have found that typically C─C bonds involved in the reconstruction of the dislocation core are not responsible for deep intragap states, but we see here that the C dangling bonds act differently. In fact, by pushing closer the two C atoms, each with a dangling bond in the atomic structure of Figure 5, a new reconstruction is formed after ionic relaxation. A comparison between the two structures revealed that the induced C─C bond is very stretched and thus less stable (proved by an higher total energy of the simulated cell) than the structure shown in Figure 5 and including the dangling bonds. However, the two intragap states above the Fermi level disappear when the C─C bond is formed, but the other intragap states below the Fermi level remain unaffected. Note that the discrete states that we are commenting here are very likely single-dispersive states if plotted along the k-space, as already shown for the dislocation complexes. [36] Note also that our analysis does not include charged defects; thus, further analysis is necessary to study exhaustively the stability of these structures and understand how this may impact on the defect states. Nevertheless, the appearance of the intragap states below the Fermi level in both the different recontractions simulated suggests that they are unaffected by the specific reconstructions that may be formed at the intersection of the two double SFs. Therefore, they are very likely to cause detrimental effects on the performance of the SiC devices, in analogy with the dislocation complexes terminating double and triple SFs. [36]

Conclusion
MD simulations were performed to reveal the atomistic behavior at the crossing point of two adjacent partial dislocation loops that open along the same 〈110〉 direction in two different {111} planes. A pinning-depinning process is activated, leading to the formation of a complex junction defect formed by four close 90 dislocation segments parallel to each other. Such a configuration has been analyzed and associated with the so-called dislocation forest. Moreover, superior ab initio calculations showed that the resulting structure of the junction does not cause defect states in the electronic bandgap of the material. In contrast, supposing that the same mechanism is activated also at the intersection of double-plane dislocation loops (hence double-plane SFs), we constructed another complex defect formed by 30 and 90 parallel and close-dislocation segments. The DOS calculated for this junction of double-dislocation loops reveals several intragap defect states, which are very likely detrimental to the properties of the devices. To our knowledge, this work is the first attempt to model the atomistic-level forest dislocations and predict their impact on the electronic properties of the material.