Simulation of bonding effects in HRTEM images of light element materials

The accuracy of multislice high-resolution transmission electron microscopy (HRTEM) simulation can be improved by calculating the scattering potential using density functional theory (DFT) [1–2]. This approach accounts for the fact that electrons in the specimen are redistributed according to their local chemical environment. This influences the scattering process and alters the absolute and relative contrast in the final image. For light element materials with well defined geometry, such as graphene and hexagonal boron nitride monolayers, the DFT based simulation scheme turned out to be necessary to prevent misinterpretation of weak signals, such as the identification of nitrogen substitutions in a graphene network. Furthermore, this implies that the HRTEM image does not only contain structural information (atom positions and atomic numbers). Instead, information on the electron charge distribution can be gained in addition. In order to produce meaningful results, the new input parameters need to be chosen carefully. Here we present details of the simulation process and discuss the influence of the main parameters on the final result. Furthermore we apply the simulation scheme to three model systems: A single atom boron and a single atom oxygen substitution in graphene and an oxygen adatom on graphene.

: Graphene k convergence test. On the left we see the dependence of the total energy and the electric field gradient on the number of k-points. On the right the corresponding computation time is visualized.
As we are mainly interested in the projected potential we also checked how this quantity is influ-23 enced by the number of k-points. The projected potential was calculated from a single slice parallel 24 to the beam direction. Each 2d slice was normalized to the smallest value and no cutoff was used 25 during the projection process. Four calculations (50k, 100k, 1000k and 10000k) were compared 26 to the calculation with the highest number of k-points that was performed (20000k). The range 27 of the absolute difference (left side of Figure 2) for each calculation is relatively small while the 28 absolute high of the difference plot is probably very sensitive to the normalization. On the right 29 side of Figure 2 we see the relative difference between the projected potentials V z obtained for 30 difference numbers of k-points compared to the 20000k calculation in % that was obtained by: . Surprisingly, the projected potential does not converge smoothly with 32 the number of k-points: The 100k calculation is in much better accordance to the 20000k calcu-33 lation than the 1000k calculation. One reason might be the low DFT convergence conditions that 34 were used for these calculations (charge convergence: -cc 0.0001 C and energy convergence: -ee 35 0.0001 Ry).

36
In conclusion it can be said that 100k are enough and that the projected potential is not very sen-37 sitive to the number of k-points. As the calculation time scales only linear with the number of k-38 points this parameter is not critical for the quantity we are interested in.

Size of the basis set 40
The second parameter we tested was the size of the basis set that is determined by the RKMAX 41 value. This value was increased starting from 5.5 to 9.0 with a step width of 0.5. In contrast to 42 the k-convergence test the DFT convergence conditions were increased to: -cc 0.  should be necessary to use a value of RKMAX = 9.0. For bigger systems this high value may not 47 be affordable due to strongly increasing computation time. We again calculated the projected po-48 tential for several RKMAX values and compared them to the calculation with the highest accuracy.

49
The results are visualized in Figure 4 and we see that luckily the projected potential is not very sen-50 sitive to RKMAX: The relative difference between RKMAX = 7 and RKMAX = 9 (green line on 51 the right side of Figure 4) is only 0.06% at its maximum.

52
RKMAX=9.00 8.00−9.00 7.00−9.00 6.00−9.00 5.50−9.00 RKMAX=9.00 8.00−9.00 7.00−9.00 6.00−9.00 5.50−9.00    Figure 6 where the time consumed to create one 66 slice is printed against the GMAX value. We also checked the influence of the GMAX parameter on the projected potential and found that 68 for GMAX = 10 the relative error is up to 4% while for the default value of 12 the error is smaller 69 than 0.25%.

71
Graphene is a true 2d material in 3d space while our DFT calculations were always using 3d 72 unit cells and 3d periodic boundary conditions. How large do we have to make the unit cell in z-73 direction to 'isolate' the graphene layers? The graphene layer separation was increased from 5 Å 74 to 35 Å with a stepsize of 5 Å and the total energy and electric field gradient were plotted against 75 the layer separation in Figure 8. One problem that complicates the interpretation of this study is the 76 fact that each calculation was using a different k-grid. We always were using 500k for the WIEN2k 77 input but the software constructed different k-grids with variating number of k-points for each case 78 so the results are not directly comparable. Nevertheless, it appears that a vacuum separation of 79 20 − 25 Å is sufficient to suppress the interaction between two neighboring graphene layers.