Substitutional Si impurities in monolayer hexagonal boron nitride

TL;DR

This study reports the first direct observation of silicon atoms substituting boron vacancies in monolayer hexagonal boron nitride, combining experimental STEM imaging with DFT calculations to analyze stability and electronic properties.

Contribution

It provides the first experimental evidence and theoretical analysis of substitutional silicon impurities in monolayer h-BN, revealing stable configurations and electronic effects.

Findings

Silicon atoms fill boron vacancies in h-BN.

Si_B^{+1} is the most stable substitutional configuration.

Substitutional Si introduces defect levels in the band gap.

Abstract

We report the first observation of substitutional silicon atoms in single-layer hexagonal boron nitride (h-BN) using aberration corrected scanning transmission electron microscopy (STEM). The medium angle annular dark field (MAADF) images reveal silicon atoms exclusively filling boron vacancies. This structure is stable enough under electron beam for repeated imaging. Density functional theory (DFT) is used to study the energetics, structure and properties of the experimentally observed structure. The formation energies of all possible charge states of the different silicon substitutions (Si$_\mathrm{B}$, Si$_\mathrm{N}$ and Si$_\mathrm{{BN}}$) are calculated. The results reveal Si$_\mathrm{B}^{+1}$ as the most stable substitutional configuration. In this case, silicon atom elevates by 0.66{\AA} out of the lattice with unoccupied defect levels in the electronic band gap above the Fermi level. The formation energy shows a slightly exothermic process. Our results unequivocally show that heteroatoms can be incorporated into the h-BN lattice opening way for applications ranging from single-atom catalysis to atomically precise magnetic structures.

Equations

$$\begin{gathered} E^f[X^q]=E_{tot}[X^q]-E_{tot}[host] \\ -\sum n^i{\mu }_i+q[E_F+E_{VBM}]+E_{corr}, \end{gathered}$$

Full text

††preprint: AIP/123-QED

The study of two-dimensional (2D) materials has since the introduction of graphene Novoselov et al. (2004) opened an active research field in material science. Graphene was quickly followed by other 2D material such as hexagonal boron nitride (h-BN) Pacilé et al. (2008) and transition metal dichalcogenides Wang et al. (2012), that in contrast to graphene exhibit an electronic band gap. Among 2D materials, h-BN has attracted attention due to its high thermal and chemical stability, high thermal conductivity and low dielectric constant, besides its wide electronic band gap Cassabois, Valvin, and Gil (2016), although due to its honeycomb-like arrangement of $sp^{2}$-hybridized structure that is almost a perfect match to graphene, it is most often regarded just a suitable substrate for graphene-based applications Geim and Grigorieva (2013).

Defects play a crucial role in semiconductors in determining the applicability of the material. For example, vacancies and impurities change the electronic and optoelectronic properties by adding localized defect levels into their band gaps Oba and Kumagai (2018); Kapetanakis et al. (2016); Kotakoski et al. (2011) providing charge carrier traps and/or combination centers Hu et al. (2018). Although silicon is the most often encountered impurity atom in graphene samples Zhou et al. (2012); Ramasse et al. (2013a); Jalili et al. (2018), it has not been until now observed in h-BN. In contrast, both oxygen and carbon have been found in single-layer h-BN, probably due to electron beam damage, during a STEM experiment Krivanek et al. (2010). On the one hand, Si impurities in h-BN could be interesting from material engineering point of view, especially for applications in electronics Meyer et al. (2009); Kim et al. (2018); Murata et al. (2013); Majety et al. (2013), quantum computing Sajid, Reimers, and Ford (2018) and spintronic devices Tang and Cao (2010); Asshoff et al. (2018); Wang et al. (2018). On the other hand, they could be detrimental when h-BN is used as a gate dielectric in field-effect transistors. Previous theoretical studies have shown that silicon substitution in boron vacancy is more stable than substitution in nitrogen vacancy Liu et al. (2014); Mapasha, Igumbor, and Chetty (2016); Mapasha et al. (2016). Both $+1$ and $-1$ have also been considered in addition to the neutral state Mapasha et al. (2016). However, the number of valence electrons and the position of the defect levels have been ignored.

In this Letter, we show the atomic resolution STEM images of substitutional silicon impurities in h-BN with DFT calculations and image simulations revealing the details of the atomic configuration. The formation energies for each of the possible silicon substitutions (in boron, nitrogen and double vacancy) are calculated with different charge states. In accordance with the experiments that show exclusively impurity atoms in the B lattice site, our simulations reveal Si${}_{\mathrm{B}}^{+1}$ as the configuration with the lowest formation energy. Our results demonstrate the possibility of incorporating heteroatoms into h-BN opening way for atomic-scale engineering of the material for applications.

Single-layer h-BN suffers from significant electron beam damage in transmission electron microscopy experiments Meyer et al. (2009); Jin et al. (2009); Kotakoski et al. (2010), caused by a combination of knock-on processes and ionization damage. The vacancies created during electron irradiation tend to grow fast into triangular holes, but can also be filled by atoms from the ubiquitous hydrocarbon-based contamination covering most samples, as presumably happened in the study presented in Ref. Krivanek et al. (2010). In contrast to such non-intrinsic defects created during the experiment, we discuss here intrinsic impurity atoms that were found in the samples after preparation with no additional processing and only a minimal electron dose.

We point out that unlike most h-BN samples, that are prepared via mechanical exfoliation, ours were grown via chemical vapor deposition (Graphene Laboratories, Inc.), which may explain why Si impurities in h-BN have not been reported until now. The samples were directly transferred onto golden transmission electron microscopy grids with perforated amorphous carbon membrane (Quantifoil®) without a polymer. The copper was etched in FeCl over night and the samples were baked in vacuum at 150∘C overnight before being inserted into the microscope. Fig. 1a shows an atomically resolved medium angle annular dark field (MAADF) image of the suspended h-BN membrane with four silicon impurity atoms in the h-BN lattice. All images were acquired with the Nion UltraSTEM 100 microscope Krivanek et al. (2008) in Vienna at 60 keV with near-ultrahigh vacuum conditions at the objective area (pressure $<10^{-9}$ mbar) Leuthner et al. (2018). The beam convergence semiangle was 30 mrad and the MAADF detector angular range was 60-200 mrad. Typical beam current of the device is on the order of 30 pA.

The silicon atom can be easily distinguished by its higher contrast as compared to boron and nitrogen Krivanek et al. (2010). In all images, Si impurities are on boron sites (see Fig. 1b). The image intensities are shown in Fig. 1c along the path marked by arrows in Fig. 1b, normalized to the intensity corresponding to a single boron atom. Intensity ratio between Si and B is expected to be $(14/5)^{1.64}=5.41$ Krivanek et al. (2010), close to our experimental value of 5.25. Although the substituted silicon atom is stable enough to allow its repeated scanning, it can not sustain the electron dose required for electron energy loss spectroscopy Ramasse et al. (2013b) or energy-dispersive X-ray spectroscopy to independently confirm the chemical identity of the individual impurity atoms. Since the hydrocarbon-based contamination covering much of the sample also contains silicon, spectra recorded over larger areas would remain inconclusive. Fig. 1d shows a simulated MAADF image of the relaxed structure (+1 charge state) with a good agreement between the projected Si-N distances between the nearest nitrogen atoms and the impurity in calibrated experimental and simulated images. For both cases the projected Si-N distance is around 1.55 Å. The distance between neighboring nitrogen atoms shows an increase from 2.51 Å in pristine h-BN to 2.74 Å around the impurity, consistent with the optimized model structure.

The typical hydrocarbon-based contamination that covers practically all graphene and h-BN samples contains a large number of silicon atoms. In fact, it is possible to dope such structures with Si simply by creating vacancies into them at elevated temperatures Inani et al. (2019). Hence, it is natural to assume that this is also the source of those Si atoms found in our h-BN samples.

We turn to density functional theory (DFT) calculations (as implemented in the Vienna ab initio simulation package (VASP) Kresse and Hafner (1993)) to try to understand why they are exclusively found in the B lattice sites. The electron exchange and correlation was treated by Perdew-Burke-Ernzerdorf (PBE) functional Perdew, Burke, and Ernzerhof (1996). The total energy of the system was calculated via the pseudopotential-momentum-space formalism using projector-augmented-wave (PAW) method Kresse and Joubert (1999). The Kohn-Sham wavefunctions are expanded over plane-wave basis sets with the kinetic energy cut off set to 525 eV. Converged locally optimized configurations and formation energies were found for a supercell of $8\times 8\times 1$. The interlayer vacuum space of 43.46 Å was selected according to "special vacuum" proposed in Refs. Komsa et al. (2014, 2018). The Brillouin-Zone integration was done over a $\Gamma$-centered $5\times 5\times 1$ k-point mesh. The damped molecular dynamics method was used to optimize the ionic degrees of freedom until residual forces were below 0.01 eV/Å. We point out that although it is known that the band gaps calculated using PBE underestimate the true band gap of semiconductors, it yields formation energies similar to those calculated with the HSE (Heyd-Scuseria-Ernzerhof) Heyd, Scuseria, and Ernzerhof (2003) formalism Berseneva, Krasheninnikov, and Nieminen (2011). This allows rescaling of the electron chemical potential of PBE calculations using the difference in the band gap obtained from the two methods leading to a significant saving in the computational cost. Here, HSE (HSE06 functional with 0.25 fraction of exchange Krukau et al. (2006)) was used to calculate the band gap of bulk h-BN (5.72 eV as compared to 4.48 eV calculated with PBE) for this purpose.

For each possible impurity site in h-BN, we first structurally optimize the structure, and then calculate its electronic properties and formation energy. The calculations are repeated for different charge states. Fig. 2a-b show the top and side view of the structures for neutral state of SiB. The N-Si bond length is found to be $\sim{1.72}$ Å, which is significantly longer than $\sim{1.45}$ Å between B and N in the pristine structure and the 1.55 Å measured from the experimental images. The silicon atom rises 1.241 Å above the h-BN plane. The side view of the structures corresponding to different charge states are shown in Fig. 2c-d. As for the +1 charge state, the N-Si bond length drops to 1.63 Å. For -1 charge state the bond length becomes 1.82 Å. Simply taking into account the projected distance between silicon and neighboring nitrogen atoms is not conclusive enough to estimate the charge state of the defect. The simulated MAADF images show the projected Si-N distance within $1.53-1.55$ Å, which are very close to the calibrated distances of experimental data. Negative charge elevates Si further away from the negatively charged N neighbors, whereas positive charge has the opposite effect.

The band structures of silicon substitution in boron vacancies at different charge states are plotted in Fig. 2e-g. The Fermi level is set to zero. The red (blue) lines correspond to spin up (down) band structure. For the neutral state, Si substitution adds two defect states within the band gap, where only one of the defect levels (spin down) is empty. Therefore, SiB has only two expected charge states (-1 and +1); adding further electrons or holes to the structure leads to electrons in the nearly-free-electron (NFE) state of conduction band minimum (CBM) or holes at the valence band maximum (VBM) Berseneva, Krasheninnikov, and Nieminen (2011); Huang and Wei (2011). This NFE state is estimated to be 2 Å away from the h-BN plane Huang and Lee (2012). In the +1 and -1 charge states the system becomes spin-unpolarized. An interesting case is the band structure of Si${}_{\mathrm{B}}^{-1}$ (Fig. 2f). Here, the defect level is so close to the Fermi level that added electron is almost free and should be easy to extract.

We calculate the defect formation energy using the supercell method Van de Walle et al. (1989). Here, the formation energy at a charge state $q$ is defined as

[TABLE]

where $E_{tot}[X^{q}]$ is the total energy of the supercell containing a defect or impurity $X$, $E_{tot}[host]$ is the total energy for equivalent supercell of perfect crystal. $n_{i}$ are the number of atoms which are added ($n_{i}>0$) or removed ($n_{i}<0$) from the supercell and $\mu_{i}$ are the chemical potentials of the constituent atoms $i$. The formation energy is expressed as a function of electron chemical potential (i.e., Fermi energy $E_{F}$ with respect to the valence band maximum of the pristine structure). $E_{corr}$ corresponds to all spurious electrostatic corrections due to employing the supercell method. We calculate the chemical potential for boron and nitrogen atoms as the total energy of $\beta$-rhombohedral boron (per atom), containing 106 atoms per unit cell van Setten et al. (2007), and half of the chemical potentials of the nitrogen molecule (N2). The chemical potential of silicon has been calculated from total energy of bulk silicon (per atom).

The calculated SiB formation energies in N-rich condition are shown in Fig. 3. The formation energy for the neutral defect is around 0.46 eV, which is higher than the previously reported value of -0.29 eV, possibly due to the small vacuum size (15 Å) used in the earlier calculations Mapasha et al. (2016); Mapasha, Igumbor, and Chetty (2016). However, the most stable charge state for Si${}_{\textrm{B}}$ substitution is +1 which undergoes a slightly exothermic process. The defect energy transition states for $e$(+1/0) and $e$(0/-1) are expected to be at 0.61 eV and 2.96 eV with respect to VBM as calculated with the Perdew-Burke-Ernzerdorf (PBE) Perdew, Burke, and Ernzerhof (1996) functional. Taking into account the difference in the band gaps calculated with PBE and HSE, as described above, we would expect the actual transition levels to be 1.23 eV and 3.58 eV above VBM. Both defect levels are considered deep and expected to influence the optical properties of h-BN.

For h-BN layer which is on top of Cu(111), we expect that the copper (as during growth) Fermi level lies around $1.4-1.5$ eV above the VBM of h-BN Laskowski, Blaha, and Schwarz (2008); Zhang et al. (2018) and therefore the SiB should be in neutral state. However it has been shown that under electron beam, BN-nanotubes are positively charged, which is attributed to emission of secondary electrons that shift the Fermi level to VBM and the SiB switches to +1 state Wei et al. (2013).

Fig. 4a-b show the locally optimized structure of silicon substitution in nitrogen site in h-BN in its neutral state. Compared to SiB, in SiN the silicon atom is extruded much higher out from h-BN. The distance between the silicon and the h-BN plane is $\sim{2.0}$ Å and the bond length between silicon and neighboring boron atom is 1.956 Å. This buckling could be attributed to the electrostatic repulsion between silicon and neighboring boron atoms. Bader analysis Tang and Cao (2010) of charge density shows that the silicon atom has lost all valence electrons while the neighboring boron atoms have lost a fraction of their electrons to the neighboring nitrogen atoms. Adding an electron (Fig. 4) would again slightly elevate the silicon atom. However, most of the charge is again transferred to the nitrogen atoms.

The formation energy for SiN in boron rich environment is shown in Fig. 4e. The formation energy for neutral defect is around 4.86 eV, which high compared to silicon in boron vacancy. Based on band structure calculation (Supplementary Material), possible charge states range from -1 to +2. The most stable case is when the silicon atom is in neutral state. The transition level $e$(0/-1) is at 1.90 eV as calculated with PBE. By rescaling the Fermi energy based on the calculated HSE band gap, we expected that this transition energy would rise to 2.51 eV.

We also calculated formation energy for silicon substitution in a double vacancy, where adjacent boron and nitrogen atoms are missing. The optimized structures are shown in Fig. 5. In this case the h-BN monolayer stays almost flat with the silicon-boron bond length of $\sim{2.09}$ Å and silicon-nitrogen bond length of $\sim{1.78}$ Å. One added electron changes the bond length for Si-B and Si-N to $\sim{1.96}$ Å and $\sim{1.88}$ Å, respectively. From the band structure calculation (Supplementary Material) it is evident that both +1 and +2 states are possible. However, the silicon is pushed even more toward the nitrogen atoms due to electrostatic repulsion between boron atoms and the impurity. The bond lengths in this case are $\sim{2.20}$ Å and $\sim{1.73}$ Å for Si-B and Si-N bonds, respectively.

The calculated formation energies are shown in Fig. 5e. Interestingly, for the neutral structure of SiBN, the formation energy is around 5.25 eV, which is slightly higher than SiN substitution despite the flat structure. The +2 charge state has the highest formation energy, and the PBE calculated transition levels are 0.51 eV for (+1/0) and 3.11 eV with respect to VBM. For HSE method of calculation, we would expect transition levels for $e$(+1/0) and $e$(0/-1) to be 1.12 eV form VBM and 2.00 eV from CBM.

We further calculated the migration barrier for SiB to pass from one side of the h-BN plane to the other, a process that has been recently observed in graphene Hofer et al. (2019). The barrier was calcalated with the nudged elastic band method Henkelman and Jónsson (2000). The estimated energy barrier of 2.42 eV is much higher than 1.08 eV for silicon substitution in graphene Susi et al. (2014), which could be attributed to the ionic bonding between silicon and h-BN. The same transition can be achieved by providing the Si atom 3.28 eV of kinetic energy, as confirmed with molecular dynamics calculation with DFT (0.5 fs time step). This amount of energy can be transferred to a static Si atom by an electron with a kinetic energy of 40.5 keV and easily achieved under our experimental conditions (acceleration voltage 60 kV). Thus most Si atoms are expected to face away from the electron beam.

In this work we present the direct experimental observation of Si impurities in a free standing monolayer of hexagonal boron nitride (h-BN) using scanning transmission electron microscopy. Our density functional theory calculations show that Si atoms have the lowest formation energy in a boron vacancy, which is also the only experimentally observed configuration. Our simulations indicate +1 as the most likely charge state. Although the origin of Si atoms in our samples remains unknown, the results expands the number of observed impurity atoms in h-BN from the previously reported C and O to an element from the third row of the periodic table. This shows that heteroatom doping of h-BN with heavier elements is possible similar to graphene Susi et al. (2017); Tripathi et al. (2018), opening the way towards applications ranging from single-atom catalysis to atomically small magnetic structures.

We acknowledge the Austrian Science Fund (FWF) for funding through projects I3181 and P31605 as well as the Vienna Science and Technology Fund through project MA14-009. HPK acknowledges Academy of Finland for the support under Project No. 311058. We further acknowledge the Vienna Scientific Cluster for generous grants of computational resources. ReferencesKrivanek et al. (2008); Leuthner et al. (2018); Kresse and Hafner (1993); Perdew, Burke, and Ernzerhof (1996); Kresse and Joubert (1999); Komsa et al. (2014, 2018); Heyd, Scuseria, and Ernzerhof (2003); Berseneva, Krasheninnikov, and Nieminen (2011); Krukau et al. (2006); Koch (2002); Henkelman and Jónsson (2000) are cited in the supplement.

Supplementary Material

Substitutional Si impurities in monolayer hexagonal boron nitride

—

Supplementary Material

Mohammad Reza Ahmadpour Monazam

Ursula Ludacka

Hannu-Pekka Komsa

Jani Kotakoski

Methods

The samples were commercially available single-layer h-BN grown via chemical vapor deposition on copper by Graphene Laboratories, Inc. They were directly transferred onto golden transmission electron microscopy grids with perforated amorphous carbon membrane (QUANTIFOIL®) without the use of a polymer, which decreases the amount of contamination on the samples. The copper was etched in a bath of FeCl over night and the samples were cleaned with deionized water and isopropyl alcohol. Samples were baked in vacuum at 150∘C for at least eight hours before being inserted into the microscope.

We acquired the experimental data using a Nion UltraSTEM 100 microscope Krivanek et al. (2008) with a cold field-emission electron gun operated at 60 keV. The near-ultrahigh vacuum conditions at the objective area (pressure below $10^{-9}$ mbar) around the sample ensure a minimum influence of chemical reactions Leuthner et al. (2018) on the sample during observation. The beam convergence semiangle was 30 mrad and the used medium angle annular dark field detector angular range 60-200 mrad. Typical beam current of the device is on the order of 30 pA.

We used density functional theory as implemented in the Vienna ab initio simulation package (VASP) Kresse and Hafner (1993). The electron exchange and correlation was treated by Perdew-Burke-Ernzerdorf (PBE) functional Perdew, Burke, and Ernzerhof (1996). The total energy of the system was calculated via the pseudopotential-momentum-space formalism using projector-augmented-wave (PAW) method Kresse and Joubert (1999). The Kohn-Sham wavefunctions are expanded over plane-wave basis sets with the kinetic energy cut off set to 525 eV. A supercell of $8\times 8\times 1$ was employed to study different defect states in the membrane with the assumption of minimizing the lateral interaction of the defect with its periodic images. The interlayer vacuum space of 43.46 Å was selected according to "special vacuum" proposed in Refs. Komsa et al. (2014, 2018). The results were compared to those calculated with a supercell of $6\times 6\times 1$ and vacuum size of 30.27 Å. The locally optimized configurations and formation energies were in good agreement for the two different system sizes. The Brillouin-Zone integration was done over a $\Gamma$-centered $5\times 5\times 1$ k-point mesh. The damped molecular dynamics method was used to optimize the ionic degrees of freedom until the residual forces were below 0.01 eV/Å. Although it is known that the band gaps calculated using PBE underestimate the true band gap of semiconductors, we restricted our calculation to the level of PBE due to the agreement between PBE formation energies and those calculated with the HSE formalism Heyd, Scuseria, and Ernzerhof (2003); Berseneva, Krasheninnikov, and Nieminen (2011). Therefore one only needs to re-scale the electron chemical potential using the difference in the band gap obtained from the two methods. Due to the computational cost, we carried out only one HSE calculation for bulk h-BN for estimating the band gap. The size of the band gap in this case is 5.72 eV as compared to 4.48 eV calculated with PBE. The HSE calculation is performed using HSE06 functional with 0.25 fraction of exact exchange Krukau et al. (2006). For STEM image simulations, we used the QSTEM package Koch (2002), where all parameters were set to correspond to our experimental setup. The energy barrier estimation is based on the nudged elastic band (NEB) method implemented in VASP Henkelman and Jónsson (2000). A set of calculations with five images between the initial and final configurations were performed. The standard dynamic calculation was performed using in DFT-based molecular dynamics with Nosé-thermostat ; an increasing initial vertical velocity toward the h-BN plane is applied to the silicon atom until it passed through the membrane. The time step is set to 0.5 fs.