Magnetic Multilayer Edges in Bernal-Stacked Hexagonal Boron Nitride

In \FigrefMonolayer, we present the N-edge of a monolayer h-BN sheet. The ferromagnetic (FM) and antiferromagnetic (AFM) states are presented in \FigrefMonolayer(a,b), respectively. The isosurface plots for magnetization with the isovalue $n_{\uparrow}-n_{\downarrow}=\pm 0.02\ \left|e\right|/a_{0}^{3}$ show a spatially localized spin polarization on the dangling $sp^{2}$ orbitals on the edge N atoms. Both states are lower in energy than the nonmagnetic state, and the ferromagnetic state is the ground state [\TabrefEdges]. In both the FM and AFM states, the magnetization per edge N atom is $1$ $\mu_{\text{B}}$. There is negligible relaxation ($<0.01\ \mathring{\text{A}}$) at these edges arising from the change in the spin configuration. The spin-resolved densities-of-states (DOS) plots show that the ferromagnetic edge is a half-metal, i.e. the majority spin channel is insulating and the minority spin channel is a metal, which should give rise to a fully spin-polarized edge current and potentially enable spintronic devices such as tunnel junctions, spin filters, diodes and transistors (Prinz, 1998; Wolf et al., 2001; Žutić et al., 2004; Felser et al., 2007). The antiferromagnetic state is a semiconductor with a gap of 0.51 eV. Because the spin configuration and the electronic structure are closely linked, these edges may be useful in metal spintronics as well as semiconductor spintronics applications (Prinz, 1998; Wolf et al., 2001; Žutić et al., 2004; Awschalom and Flatté, 2007; Felser et al., 2007; Yazyev, 2010; Awschalom et al., 2013). Our findings on the monolayer h-BN agree well with the existing computational studies (Barone and Peralta, 2008; Du et al., 2009; Lai et al., 2009; Si et al., 2014; Deng et al., 2017).

When the bilayer h-BN is formed in the AB stacking, the N-edge remains “open” upon relaxation, i.e. no interlayer bonds form. This results in essentially two copies of the monolayer h-BN sheet including the edge. This bilayer configuration can be repeated along the out-of-plane direction to form a bulk N-edge. The total energies of the FM state and the AFM state (denoted AFM$\uparrow\downarrow$) for the bilayer and bulk N-edge are listed in \TabrefEdges. The FM state remains the ground state and is preferred to the AFM$\uparrow\downarrow$ state by a similar difference in energy, compared to that of the monolayer h-BN. If the consecutive layers are each ferromagnetic but spin-polarized in the opposite direction with respect to each other (denoted AFM${}_{\downarrow}^{\uparrow}$), the total energy is slightly higher than the FM case. All of these results reinforce our finding that the neighboring layers are electronically and magnetically largely independent from each other.

In \FigrefBilayer_mag, the details of the atomic and magnetic configuration of the bilayer and bulk AB1-h-BN N-edge are presented. In \FigrefBilayer_mag(a), the isosurface plots for magnetization (isovalue $n_{\uparrow}-n_{\downarrow}=\pm 0.02\ \left|e\right|/a_{0}^{3}$ ) for the two AFM configurations of the bilayer edge are displayed (AFM${}_{\ \uparrow\ \downarrow}^{\uparrow\ \downarrow}$=AFM$\uparrow\downarrow$ and AFM${}_{\ \downarrow\ \downarrow}^{\uparrow\ \uparrow}$=AFM${}_{\downarrow}^{\uparrow}$). In \FigrefBilayer_mag(b,c,d), the DOS plots show that the bilayer FM edge is a half-metal, the bilayer AFM$\uparrow\downarrow$ state is a semiconductor, the bilayer AFM${}_{\downarrow}^{\uparrow}$ state is a half-semiconductor, and the bulk AFM${}_{\downarrow}^{\uparrow}$ state is a metal. The bilayer AFM${}_{\downarrow}^{\uparrow}$ state is similar to a bipolar magnetic semiconductor in which the valence band maximum (VBM) and the conduction band minimum (CBM) of one spin channel are at a higher energy than the those of the other spin channel (Li et al., 2012; Li and Yang, 2016). Here, the VBMs of the two spin channels are aligned but the CBMs are not aligned, leading to unequal gaps in the two spin channels, which, in the case of slight electron doping, would give rise to a spin-polarized edge current that would also be spatially confined to the top layer. Manipulating spins in semiconductors (semiconductor spintronics) is a rapidly developing field with promising applications in quantum computing (Awschalom and Flatté, 2007; Awschalom et al., 2013).

Because the energy differences between the FM and AFM configurations are low, thermodynamic considerations are needed to determine the expected configurations at a given temperature. Because spin polarization is localized to each edge N atom, we can use a discrete lattice model with interactions between neighbors instead of an itinerant spin model. In order to extract the interaction strengths, we have computed the energies of 10 single-layer, 5 bilayer and 5 bulk spin configurations. Because of the high precision required from these calculations, only the supercells compatible with the $12\times 1$ $k$-point grid are considered, which are $1\times 1$, $2\times 1$, $3\times 1$, $4\times 1$ and $6\times 1$. The spin configurations, total energies and the electronic gaps of these systems are listed in Table S1 of the Supplemental Material. These calculations have allowed us to extract the in-plane interaction energies for nearest, next-nearest and next-next-nearest neighbors ($J_{1}$, $J_{2}$ and $J_{3}$, respectively), as well as the out-of-plane interaction energy ($J_{4}$). The schematic representations of these interactions are presented in Figure S1, and the fitted values with the 90% confidence intervals are listed in Table S2 of the Supplemental Material. We find that zero is within the confidence interval for $J_{3}$ and $J_{4}$, meaning that the open N-edge of each layer can be thought of as an independent one-dimensional spin chain with nearest and next-nearest neighbor interactions with the Hamiltonian:

$$H=J_{1}\sum_{i}\boldsymbol{S_{i}}\cdot\boldsymbol{S_{i+1}}+J_{2}\sum_{i}\boldsymbol{S_{i}}\cdot\boldsymbol{S_{i+2}}-\mu_{\text{B}}H\sum_{i}S_{i}^{z},$$

(1)

where $\boldsymbol{S_{i}}$ is a unit vector with $x,y,z$ components that represents the spin at a lattice site $i$, $J_{1}=-17$ meV is the nearest neighbor interaction energy, $J_{2}=-12$ meV is the next-nearest neighbor interaction energy, $\mu_{\text{B}}$ is the Bohr magneton, and $H$ is the external magnetic field. Because spin–orbit interaction in h-BN is extremely weak (Zollner et al., 2019) (two or three orders of magnitude smaller than $J_{1}$ and $J_{2}$), we do not expect magnetic anisotropy in this system above $\sim 10$ K (Yazyev and Katsnelson, 2008). As a result of this lack of anisotropy, there is no easy axis present in this system which leads to a continuous spin model (Heisenberg) instead of a discrete spin model (Ising); furthermore, the external magnetic field can be assumed to apply in the $z$-direction without loss of generality.

This model corresponds to classical one-dimensional isotropic Heisenberg model with nearest and next-nearest neighbor interactions. For such a model, according to the Mermin–Wagner theorem, in the absence of an external magnetic field, the ferroelectric ground state is not stable for any temperature greater than absolute zero (Mermin and Wagner, 1966). Although an exact analytical solution to this model is not known, when both $J_{1}$ and $J_{2}$ are negative, i.e. both the interactions are ferromagnetic, its behavior is analogous to the Heisenberg model with only nearest-neighbor interactions (Majumdar and Ghosh, 1969a, b) whose exact solution is known (Fisher, 1964). In order to determine the behavior of the expected value of spin ($\left\langle S_{z}\right\rangle$) and the spin pair correlation function ($\xi_{z}$) in our system where $J_{1}$ and $J_{2}$ have specific numerical values, we have conducted a Monte Carlo study a modified version of the Metropolis algorithm (Metropolis et al., 1953). In the usual Metropolis algorithm for the Heisenberg model, at each simulation step, an attempt is made to randomize one of the spins (Slanic et al., 1991). However, this leads to a critical slowing down where the decorrelation time diverges near the critical temperature (Barkema and Newman, 1997; Dogan and Ismail-Beigi, 2019), which in this case is 0 K. To overcome this problem to some degree, instead of choosing the proposed direction of the spin at random, we choose it according to the formula

$$\Delta\theta=\cos^{-1}\left(K^{-1}\log(\exp(-K)+2r\sinh(K))\right),$$

(2)

where $\Delta\theta$ is the angle between the new spin and the old spin, $K=\nicefrac{{\left(J_{1}+J_{2}\right)}}{{\left(k_{\text{B}}T\right)}}$ and $r$ is a random number between 0 and 1 (Lurie et al., 1974). This attempt is then accepted or rejected based on the Boltzmann factor of the energy difference between the proposed and original configurations, which generates a sampling of the configuration space that obeys the detailed balance condition (Dogan and Ismail-Beigi, 2019).

We have obtained a collection of results by running these simulations with free boundary conditions where the number of spins is chosen between $N=20$ and $N=200$ to be at least a few times the spin pair correlation length, the temperature varies between $T=12$ K and $T=232$ K and the external magnetic field varies between $H=0$ T and $H=100$ T. As an example, the correlation functions and the evolution of the average site energy for a simulation for $H=0.1$ T and $T=70$ K are presented in Figures S2, S3 and S4 of the Supplemental Material. $\left\langle S_{z}\right\rangle$ and $\xi_{z}$ for all temperatures and magnetic fields are plotted in \FigrefHeisenberg_res. The error bars are obtained from the standard deviations over ten simulations conducted for each $H,T$ pair. According to our results on $\left\langle S_{z}\right\rangle$ [\FigrefHeisenberg_res(a)], it is possible to magnetize significantly these edges with external fields of 4 T and above. This is higher than the field for magnetic saturation (about 1 T) found in an experiment (Si et al., 2014). However, in this experiment, the observed magnetization may have originated from not only edges but also other defects such as vacancies and pores, which are also known to have magnetic properties (Dogan et al., 2020). Although the average spin for $H=0$ T is zero for all nonzero temperatures, the spin pair correlation length is not zero [\FigrefHeisenberg_res(b)] and increases exponentially as $T$ approaches zero. As a result, the expectation value of the site energy decreases as $T\rightarrow 0$ [Figure S5 of the Supplemental Material]. At low-temperature experiments (liquid N₂ temperatures or lower), correlation lengths comparable to the edge length may produce an enhanced magnetization in the presence of small and intermediate-sized nanopores. The heat capacity ($C_{v}$) and magnetic susceptibility ($\chi$) of this system for all $H,T$ pairs are also presented in Figures S6 and S7 of the Supplemental Material, whose behaviors are analogous to similar models (Bonner, 1978). Experiments that can isolate the magnetic behavior of edges in h-BN sheets, such as electron energy loss spectroscopy (EELS) (Alem et al., 2012) and scanning tunneling microscopy (STM) (Wong et al., 2015) are needed to test these findings. The ability to induce the ferromagnetic state at the edge via an external field, which is also the state with less electrical resistance (\FigrefMonolayer), indicates the potential for magnetoresistive effect (Prinz, 1998; Wolf et al., 2001; Felser et al., 2007). Furthermore, if edges could be stacked so that FM and AFM edges alternate, giant magnetoresistance may be observed (Prinz, 1998; Wolf et al., 2001; Felser et al., 2007).

We now turn to other possibilities of bilayer edge reconstructions in h-BN which result in “closed” edges. First, as reference, we compute the relaxation of the N-edge in AA- and $\text{AA}^{\prime}$-h-BN. In AA-h-BN, edge N atoms in neighboring layers are directly aligned, and relax out-of-plane to form an $sp^{2}$ type bond [\FigrefBilayer_AA(a,b)]. The degree of out-of-plane relaxation is quite large, where the highest interlayer distance becomes $6.66$ Å. In $\text{AA}^{\prime}$-h-BN, due to the inequivalence between the neighboring layers, an N-edge in the bottom layer must align with a B-edge in the top layer. These edges form an interlayer bond which creates a structure that can be understood as a half BN nanotube [\FigrefBilayer_AA(d,e)]. Our findings on the closed $\text{AA}^{\prime}$-h-BN edge agrees with previous calculations and observations (Alem et al., 2012). In both of these edge reconstructions, each atom makes three bonds, which eliminates unpaired electrons and results in nonmagnetic insulators [\FigrefBilayer_AA(c,f)].

In the case of AB-h-BN, the starting configuration of the bilayer N-edge can be created in two different ways by terminating the layers at different N rows. The first is presented in \FigrefBilayer_AB1(a). When the atoms are allowed to relax, no interlayer bonds form, leading to the open edge structure [\FigrefBilayer_AB1(b,c)] which we have discussed above. This bilayer edge is the one we were able to identify in the HRTEM study (Dogan et al., 2020). However, it is conceivable that while the edges form, the B atom bonded to the edge N atom in the top layer finds the edge N atom in the bottom layer and forms a bridge between the two edges. This reconstruction is shown in \FigrefBilayer_AB1(d). This highly distorted $2\times 1$ edge structure has not been observed, but may arise under different experimental conditions. Because each atom makes three bonds, eliminating unpaired electrons, the DOS plot shows a nonmagnetic insulator [\FigrefBilayer_AB1(e)].

The second way of creating a starting configuration for bilayer N-edge in AB-h-BN is presented in \FigrefBilayer_AB2(a). When the atoms are allowed to relax, large out-of-plane relaxations occur, leading to the closed edge structure shown in \FigrefBilayer_AB2(b). The large out-of-plane relaxations are accompanied by a significant in-plane contraction perpendicular to the edge, which was not observed in the HRTEM study (Dogan et al., 2020). It is possible that while the edges form, the B atom bonded to the edge N atom in the bottom layer finds the edge N atom in the top layer and forms a bridge between the two edges. This reconstruction is shown in \FigrefBilayer_AB2(d). Although this forms a closed edge with a $2\times 1$ reconstruction, the in-plane relaxations are strictly confined to the edge atoms and no significant contraction perpendicular to the edge is observed. As a result, this structure would be very difficult to identify in an HRTEM study such as Ref. Dogan et al., 2020, and therefore is not ruled out. In both of these reconstructions, because each atom makes three bonds, eliminating unpaired electrons, the DOS plots reveal nonmagnetic insulators [\FigrefBilayer_AB2(c,e)].

In real h-BN sheets, it is difficult to manufacture bilayer edges that are aligned to form a collection of stacked open N-edges as in \FigrefBilayer_AB1(b) although a preference for such alignment has been observed in our study where edges were fabricated through a high-energy electron beam (Dogan et al., 2020). Thus, it is important to investigate the behavior of a monolayer edge when its neighbor does not also have an edge nearby. To this end, we have investigated the behavior of step edges in h-BN for AA, $\text{AA}^{\prime}$ and AB stacking sequences. We have found that the monolayer edges remain intact and can be in an FM or AFM state. In \TabrefSteps, we summarize these findings. We observe that the energy difference between the FM and AFM states depends on the stacking sequence, and hence the interaction energies in 1 for each edge would be different. However, the FM state remains the ground state and the AFM state is higher in energy by similar amounts, thus the resulting models are not expected to be qualitatively different. Therefore, the N-edge in a single layer should remain magnetically active when it is sandwiched by other sheets that do not have an edge in the immediate vicinity, regardless of the stacking sequence.

In \FigrefBilayer_steps, we present the atomic and electronic structures of the four step edges we have investigated. We find that the ground state electronic structure of all four step edges is a FM half-semiconductor with band gaps in the infrared region ($0.12-0.18$ eV), and the AFM state is a nonmagnetic semiconductor with band gaps also in the infrared region ($0.31-0.47$ eV). Our findings indicate that step edges can be utilized to stabilize a magnetic edge in $\text{AA}^{\prime}$-h-BN. We also find that the amount of bending of the step edge toward the full sheet is significantly different between the AB1 and AB2 stacking sequences (0.23 Å and 0.06 Å, respectively), due to the alignment of the edge N atom in the top layer with a B atom vs. a hollow site in the bottom layer.