Spin-Orbit Dirac Fermions in 2D Systems

TL;DR

This paper introduces a new model for spin-orbit interactions in buckled 2D materials, revealing that intrinsic spin-orbit coupling can create Dirac cones, supported by DFT calculations.

Contribution

A novel theoretical model demonstrating how intrinsic spin-orbit coupling induces Dirac cones in buckled 2D systems, validated by first-principles calculations.

Findings

Intrinsic spin-orbit coupling leads to Dirac cone formation.

Model aligns with DFT calculations.

Provides microscopic understanding of spin-orbit effects.

Abstract

We propose a novel model for including spin-orbit interactions in buckled two dimensional systems. Our results show that in such systems, intrinsic spin-orbit coupling leads to a formation of Dirac cones, similar to Rashba model. We explore the microscopic origins of this behaviour and confirm our results using DFT calculations.

Equations

$$t=\left({\hat{o}}_1\cdot {\hat{o}}_2\right)V_{\pi }+\left({\hat{o}}_1\cdot \hat{d}\right)\left({\hat{o}}_2\cdot \hat{d}\right)\left(V_{\sigma }-V_{\pi }\right),$$

$$K$$

$$+4\cos \left(k_y\right)\sin \left(k_x\right)\left(\begin{array}{ccc}0& 0& -i\alpha \beta \Delta \\ 0& 0& 0\\ -i\alpha \beta \Delta & 0& 0\end{array}\right)+4\cos \left(k_x\right)\sin \left(k_y\right)\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& -i\alpha \beta \Delta \\ 0& -i\alpha \beta \Delta & 0\end{array}\right),$$

$$\frac{K}{\Delta }\to 4i{\alpha }^2\left(\begin{array}{ccc}0& -1& 0\\ 1& 0& 0\\ 0& 0& 0\end{array}\right)+2\sqrt{2}\alpha \beta k\left(\begin{array}{ccc}0& 0& -e^{i\varphi }\\ 0& 0& -e^{-i\varphi }\\ e^{-i\varphi }& e^{i\varphi }& 0\end{array}\right).$$

$$\frac{H_{\mathrm{SOI}}}{\Delta }=T\left(\frac{L_{+}\otimes s_{-}+L_{-}\otimes s_{+}}{2}+L_z\otimes s_z\right),$$

$$H_1=\left(\begin{array}{ccc}\frac{T}{2}& \mp 4i{\alpha }^2& 0\\ \pm 4i{\alpha }^2& -\frac{T}{2}& \frac{T}{\sqrt{2}}\\ 0& \frac{T}{\sqrt{2}}& 0\end{array}\right),$$

$$H_2=\left(\begin{array}{ccc}-\frac{T}{2}& \frac{T}{\sqrt{2}}& \mp 4i{\alpha }^2\\ \frac{T}{\sqrt{2}}& 0& 0\\ \pm 4i{\alpha }^2& 0& \frac{T}{2}\end{array}\right),$$

$$|{\Psi }_1^{\mathrm{I}}⟩$$

$$|{\Psi }_1^{\mathrm{II}}⟩$$

$$|{\Psi }_2^{\mathrm{I}}⟩$$

$$|{\Psi }_2^{\mathrm{II}}⟩$$

$$H_2^{\mathrm{eff}}={\left(H_1^{\mathrm{eff}}\right)}^{\ast }$$

Full text

Intrinsic SOI in two dimensions is limited to materials with heavy atoms, like some transition metal dichalcogenides, or graphene-like buckled Xenes. For both classes of materials, SOI results in band splitting, opening a gap Liu et al. (2011a, b); Xu et al. (2013); Kormányos et al. (2014, 2015); Molle et al. (2017).

In this work, we introduce a new model for two-dimensional systems with SOI. The system considered here has a fairly simple structure and can be constructed using the existing tehcnology. Gomes et al. (2012); Polini et al. (2013) By employing tight-binding formalism on a buckled homopolar square lattice with inversion symmetry, we show that instead of opening a gap, SOI can lead to a Rashba-like dispersion without an external field. After constructing a compact model, we use ab initio calculations to confirm our results.

Model. We demonstrate that SOI can lead to Dirac dispersion in two-dimensional materials if certain geometric requirements are satisfied. We begin our discussion by considering a buckled square lattice, composed of a single atomic species. Because of the buckling, we can regard the lattice as being composed of two inequivalent shifted square sublattices A and B, see Fig. 1. We set the bond length to $a$ and the (buckling) angle that it makes with the horizontal to $\theta$. This yields the unit cell size of $2\alpha a\times 2\alpha a$, where $\alpha=\cos\theta/\sqrt{2}$. From this, the dimensions of the Brillouin zone are $-\frac{\pi}{2\alpha a}\leq q_{x},q_{y}\leq\frac{\pi}{2\alpha a}$.

To keep the model as simple as possible, we only consider $p$ orbitals for each atom. We will address the effects that including $s$ orbitals can have below. Given this approximation, we now construct a tight-binding Hamiltonian, allowing transport only between the nearest neighbors. Because all the atoms and, therefore, the orbitals are identical, we set the orbital on-site energy to 0. In addition, since we are limiting hopping to the nearest neighbors, there is no coupling between the atoms of the same sublattice. The hopping matrix connecting the different sublattices, recall that the coupling term between two $p$ orbitals is given by Slater and Koster (1954)

[TABLE]

where $\hat{o}_{i}$ is the orientation of the $i$th orbital, $\hat{d}$ is the unit vector pointing from atom 1 to atom 2, and $V_{\sigma}$ ($V_{\pi}$) is the hopping integral for the $\sigma$ ($\pi$) bonds. The vectors connecting sublattice A to its four nearest neighbors in sublattice B are $a\left(n\alpha,\,m\alpha,\,-\beta\right)$ for $m$ and $n=\pm 1$, and $\beta=\sin\theta$. Thus, the two sublattices are coupled by:

[TABLE]

where $\Delta=V_{pp\sigma}-V_{pp\pi}$, and the momentum $k_{x,y}\in\left[-\pi/2,\pi/2\right]$.

Even though it is convenient to use $p_{x}$ and $p_{y}$ orbitals to write down the hopping matrix, since we are interested in including SOI in our model, it is helpful to go to a basis which is more natural for the angular momentum operators: $|1,1\rangle=-(|p_{x}\rangle+i|p_{y}\rangle)/\sqrt{2}$ and $|1,-1\rangle=(|p_{x}\rangle-i|p_{y}\rangle)/\sqrt{2}$. In addition, we will focus our attention on the corner of the Brillouin zone. The reason for this is that, according to Eq. (2), a large number of the orbitals decouple at the M point. This both simplifies our analysis and further validates our dropping of $s$ orbitals. Applying the transformation described above and expanding the trigonometric functions around $(\pi/2,\,\pi/2)$ transforms $K$ into

[TABLE]

The first term couples in-plane orbitals with opposite angular momenta; the second term couples in-plane orbitals to $p_{z}$ at finite $k$. Note that this coupling is spin-preserving.

To add the spin effects to our model, we use the standard form describing the spin-orbit coupling:

[TABLE]

where $T$ is the dimensionless ratio between the spin-orbit coupling and the hopping energy. The last term modifies the diagonal elements of the self-energy for $|1,\pm 1\rangle$ by adding (subtracting) $T/2$ if $L_{z}$ and $s_{z}$ point in the same (opposite) direction. The first tem couples $|1,1\rangle\otimes|\downarrow\rangle$ with $|1,0\rangle\otimes|\uparrow\rangle$ and $|1,-1\rangle\otimes|\uparrow\rangle$ with $|1,0\rangle\otimes|\downarrow\rangle$ with the coupling strength $T/\sqrt{2}$.

Strictly at the M point, the system decouples into four block Hamiltonians:

[TABLE]

with the top (bottom) signs acting over $|1,1\rangle\otimes|\uparrow\rangle_{A}$, $|1,-1\rangle\otimes|\uparrow\rangle_{B}$, $|1,0\rangle\otimes|\downarrow\rangle_{B}$ ($|1,-1\rangle\otimes|\downarrow\rangle_{A}$, $|1,1\rangle\otimes|\downarrow\rangle_{B}$, $|1,0\rangle\otimes|\uparrow\rangle_{B}$) and

[TABLE]

where the top (bottom) signs operate on $|1,1\rangle\otimes|\downarrow\rangle_{A}$, $|1,0\rangle\otimes|\uparrow\rangle_{A}$, $|1,-1\rangle\otimes|\downarrow\rangle_{B}$ ($|1,-1\rangle\otimes|\uparrow\rangle_{A}$, $|1,0\rangle\otimes|\downarrow\rangle_{A}$, $|1,1\rangle\otimes|\uparrow\rangle_{B}$).

Having constructed a simplified Hamiltonian around the M point of the square Brillouin zone, we now demonstrate that the presence of the atomic spin-orbit coupling leads to a formation of the Dirac cones. From the form of $H_{1}$ and $H_{2}$ it is clear that they result in a four-fold degenerate states at the M point. The most general form for these eigenstates is

[TABLE]

where $a$, $b$, and $c$ are real.

Using Eq. (3), we get that the coupling between the degenerate pairs for $H_{1}$ and $H_{2}$ is

[TABLE]

which has the form of the Rashba perturbation. This results in the lifting of the degeneracy for finite values of $k$. Since the coupling between the two states is linear in momentum, the splitting forms a Dirac cone. Superimposed on the curvature of a spinless band, this linear coupling from Eq. (8) leads to a Rashba-like dispersion. On the other hand, if the spinless band is flat, the result is a well-defined Dirac cone. Figure 2 shows two middle bands over the whole Brillouin zone for $T=1/2$, $V_{\pi}=0$, and $\theta=\pi/12$ with the linear bands at the M point.

An important feature that sets our system apart from others Dirac systems is the fact that we have only one valley. In other materials (like graphene and TMDC’s), the cones appear at the inequivalent valleys at K and K’ points of the Brillouin zone. Here, on the other hand, the cone appears at the unique $M$ point.

In the case of Rashba splitting, the Dirac cones describe helical states, where in-plane spin is linked to momentum. At first glance, it might appear that our system will behave in a similar fashion. Indeed, the Rashba-like Hamiltonian in Eq. (8) does give rise to helical states. However, the two Hamiltonians in Eq. (8) have opposite helicities, resulting in a cancellation of the spin texture. In order to avoid this cancellation one can get rid of the equivalence of $H_{1}$ and $H_{2}$ by, for example, replacing one of the sublattices by a different atomic species.

Finally, we look at the pseudo-spin, as defined by Eq. (8). The eigenstates for $H_{2}^{\mathrm{eff}}$ are $(1,\,\pm ie^{i\phi})/\sqrt{2}$. Treating the first component of the spinor as $|\uparrow\rangle$ for the pseudo-spin and the second component as $|\downarrow\rangle$, we see that the pseudo-spin forms a vortex around the origin, pointing at $90^{\circ}$ to the direction of the momentum. The direction of the pseudo-spin for the top and bottom cones are opposite to each other. In addition, the direction for $H_{1}^{\mathrm{eff}}$ is opposite to that of $H_{1}^{\mathrm{eff}}$, yielding no net pseudo-spin texture for the degenerate system. Interestingly, the vortex-like texture is different from what one observes in graphene, where the pseudo-spin is parallel to the momentum.

Discussion. Having demonstrated that our simplified system does, in fact, possess spin-split bands, let us now look at the mechanism that leads to their formation. We start by turning our attention to the composition of the matrix element in Eq. (8). We have already stated that its linear $k$ dependence results in Dirac-like bands. The second factor, $\sin 2\theta$, shows that for a flat lattice ($\theta=0$), the linear coupling vanishes. The following factor, $ac$, indicates that the band splitting takes place only if the degenerate states contain orbitals of opposite spins. The only mechanism that we have which mixes spin is SOI. Treating $T$ as a perturbation for the blocks in the first term of $H_{1}$ and $H_{2}$ reveals that the $ac$ term is proportional to $T$ for small $T$. Combined with the $\Delta$ factor, this means that in its leading-order behavior, the matrix element in Eq. (8) is proportional to $T\Delta$. Unlike Rashba effect, where the coupling element is always linear in the spin-orbit term, here, increasing $T$ yields a non-linear dependence of Eq. (8) on the ratio between the spin-orbit and hopping terms.

As we said in the introduction, the goal of this paper is to gain understanding of the effects of spin-orbit coupling on microscopic level. To do so, we focus on the composition of the states in Eq. (7). We can see that both of them include $|1,0\rangle_{B}$, but with the opposite spin. The fact that these states are coupled means that there is a finite amplitude for spin-flipping processes for $|1,0\rangle$. Figure 3 shows sets of steps for changing the spin for one of the sublattices. There are several important features in the diagram, so let us address them carefully. First, the illustration makes it clear why the buckling in the system is crucial. The path involves a transition between in-plane and out-of-plane orbitals. For a flat lattice, $p_{z}$ orbitals decouple from $p_{x}$ and $p_{y}$, making these hops impossible. This is in agreement with the $\sin 2\theta$ term in Eq. (8). Next, the path involves only one SOI-mediated transition. In other words, the spin-flipping process is first-order in $T$. Of course, as more complex paths are added, higher powers of $T$ will be included. This is equivalent to going to higher order perturbation terms in the expansion of $ac$ in Eq. (8). The atoms on the right of the vertical dashed line all belong to the same state in Eq. (7). The first hop is the transition between the states and the amplitude of this transition is proportional to the product of the relevant orbital amplitudes $ac$, as expected.

From our cartoon illustration in Fig. 3, it might appear that we have done unnecessary work by considering a square lattice because one only needs a zigzag 1D chain, similar to polyacetylene. It is possible to show that 1D chains are not sufficient by performing a full tight-binding analysis on such a chain. However, an easier way to see that 1D chain is insufficient involves rotating it $90^{\circ}$ around the longitudinal axis so that the zigzags are in the $xy$ plane. In this orientation, $p_{z}$ orbitals are decoupled and, as we determined above, one needs coupling between in-plane and out-of-plane orbitals to observe the SOI-induced band splitting.

For the sake of completeness, let us analyze the effects that including $s$ orbitals would have. If one were to construct a spinful Hamiltonian around $\Gamma$ point, similarly to how we did it around M point, one could see that the end result would not include the degenerate states, coupled at finite $k$. The two equivalent Hamiltonian blocks, while degenerate, would still be uncoupled. This means that our model would not have Dirac cones at the $\Gamma$ point. Now, we include $s$ orbitals. Because of their symmetry, $s$ orbitals couple to $p_{z}$ orbitals of the other sublattice if the lattice is buckled. This would enable the coupling between the Hamiltonian blocks and lift the degeneracy. It is important to keep in mind that because of the energy difference between $s$ and $p$ orbitals, the eigenstates will be predominantly one or the other. Since the magnitude of the band splitting depends on the presence of both orbital types, it will generally be weaker than what’s expected at the M point.

First Principles Calculations. We performed density functional theory (DFT) calculations implemented in Quantum ESPRESSO package Giannozzi et al. (2009) to simulate heavy elements Pb and Sn in a square lattice geometry. We employed Projector Augmeneted-Wave (PAW) pseudopotential with Perdew-Burke-Ernzerhof (PBE) for the exchange and correlation functional within the generalized gradient approximation (GGA) Perdew et al. (1996). The Kohn-Sham orbitals were expanded in a plane-wave basis with a cutoff energy of 70 Ry, and for the charge density a cutoff of 280 Ry was used. A $k$-point grid sampling grid was generated using the Monkhorst-Pack scheme with 16$\times$16$\times$1 points Monkhorst and Pack (1976), and a finer regular grid of 40$\times$40$\times$1 was used for spin texture calculations. For electronic band structure calculations, the spin orbit interaction was included using noncollinear calculations with fully relativistic pseudopotentials.

The optimized lattice constant $2\alpha a$ and buckling angle $\theta$ for Pb (Sn) are 3.44 Å(3.83 Å) and 44.3*∘* (27.4*∘*), respectively. We found Dirac cone or Rashba-like dispersion near $M$-point which are consistent with our TB prediction, shown in Fig. 4. As expected from Eq. (8), the lighter element Sn has a smaller Dirac slope due to smaller value of $T$ and buckling angle $\theta$.

The top and bottom Dirac cone are colored blue and red, respectively. Each band is doubly degenerate with opposite spin texture resulting a zero net spin texture. Recently, it has been proposed that centrosymmetric materials might have spin polarization Zhang et al. (2014). It was found that the top and the bottom sector of LaOBiS2 have opposite spin texture Zhang et al. (2014). Our works support such claims as Pb and Sn monolayers are centrosymmetric materials.

Conclusions We have demonstrated that intrinsic spin-orbit coupling in buckled 2D materials can lead to the appearance of the Dirac cones. Even though the effective Hamiltonian and the band structure is reminiscent of the well-known Rashba effect, the underlying physics is fundamentally different. In our study, we focused on a square lattice, but the results hold for any 2D material where in-plane orbitals couple to out-of-plane ones.