Dirac Cone in Graphene

As a roller coaster descends from a high point, potential energy turns into kinetic energy and it speeds up. More kinetic energy means greater speed. That connection is familiar from throwing a ball or pedalling a bicycle.

Near a K point in graphene, however, electrons with different energies travel at almost the same speed. How does this differ from familiar mechanics? To find out, we need to follow the electron as a wave. Using the bands derived on the previous page, let us examine how electron waves move.

01

Electrons near the contact points

In neutral graphene at zero temperature, the lower band is filled and the upper band is empty. The Fermi energy, at the boundary between them, lies at their contact points K and K′. Electrons close to these points can reach an empty state with only a small energy change.

With hopping only between nearest neighbours, the TB model gave the following two bands. Here $t>0$ is the magnitude of the hopping between neighbouring orbitals, and $f$ is the sum of three phase factors.

$$E_\pm(\mathbf k)=\pm t|f(\mathbf k)|\tag{1}$$

At K the three contributions cancel, bringing both energies to zero. Zooming in on the contact point shows the bands approaching straight lines.

×3.0
The highlighted line on the Brillouin zone gives the direction of the cut. The lower plot magnifies the boxed region of the full band by the same factor horizontally and vertically. In the horizontal coordinate $ak_x$, $a$ is the carbon–carbon bond length.
02

Measuring displacement from K

In §01, we viewed a band cross-section along Γ→K. Now consider displacements from K in different directions within the wavevector plane. Denote the displacement by $\mathbf q$; the wavevector $\mathbf k$ measured from Γ is then

$$\mathbf k=\mathbf K+\mathbf q,\qquad \mathbf q=(q_x,q_y)\tag{2}$$

Take $+q_x$ along Γ→the selected K, and $+q_y$ along the direction obtained by turning 90° counterclockwise in the wavevector plane. The components $q_x,q_y$ resolve the displacement from K along these two directions. The cross-section in §01 has $q_y=0$, and the contact itself is $\mathbf q=0$.

03

When the three contributions stop cancelling

Moving slightly away from K spoils the exact cancellation of the three contributions. Their remaining magnitude determines , so we first calculate their sum. The bond vectors $\boldsymbol\delta_1,\boldsymbol\delta_2,\boldsymbol\delta_3$ point from A to its three neighbouring B sites.

$$f(\mathbf k)=e^{i\mathbf k\cdot\boldsymbol\delta_1}+e^{i\mathbf k\cdot\boldsymbol\delta_2}+e^{i\mathbf k\cdot\boldsymbol\delta_3}\tag{3}$$

Insert and separate the phase at K from its change.

$$\begin{aligned}f(\mathbf K+\mathbf q) &=e^{i(\mathbf K+\mathbf q)\cdot\boldsymbol\delta_1}+e^{i(\mathbf K+\mathbf q)\cdot\boldsymbol\delta_2}+e^{i(\mathbf K+\mathbf q)\cdot\boldsymbol\delta_3}\\ &=e^{i\mathbf K\cdot\boldsymbol\delta_1}e^{i\mathbf q\cdot\boldsymbol\delta_1} +e^{i\mathbf K\cdot\boldsymbol\delta_2}e^{i\mathbf q\cdot\boldsymbol\delta_2} +e^{i\mathbf K\cdot\boldsymbol\delta_3}e^{i\mathbf q\cdot\boldsymbol\delta_3}. \end{aligned}$$

Consider a small neighbourhood of K, where the displacement $\mathbf q$ is small. The resulting phase change $\mathbf q\cdot\boldsymbol\delta_j$ is then small too. Keeping the constant and linear terms in $e^{ix}=1+ix+\cdots$ gives the first-order approximation:

$$e^{i\mathbf q\cdot\boldsymbol\delta_j}\simeq1+i\mathbf q\cdot\boldsymbol\delta_j$$

Apply this to each of and separate the constant and linear parts.

$$\begin{aligned} f(\mathbf K+\mathbf q) &\simeq\sum_{j=1}^3e^{i\mathbf K\cdot\boldsymbol\delta_j}(1+i\mathbf q\cdot\boldsymbol\delta_j)\\ &=\underbrace{\sum_{j=1}^3e^{i\mathbf K\cdot\boldsymbol\delta_j}}_{f(\mathbf K)=0} +i\sum_{j=1}^3 e^{i\mathbf K\cdot\boldsymbol\delta_j}(\mathbf q\cdot\boldsymbol\delta_j). \end{aligned}$$

The first sum is the cancellation at K and equals zero. What remains is the sum of the changes in the three contributions.

$$f(\mathbf K+\mathbf q)\simeq i\sum_{j=1}^3e^{i\mathbf K\cdot\boldsymbol\delta_j}(\mathbf q\cdot\boldsymbol\delta_j)\tag{4}$$

For the component calculation, use the same coordinate choice as in the geometry and TB pages.

$$\begin{aligned} \boldsymbol\delta_1&=a(0,1),\\ \boldsymbol\delta_2&=a(\sqrt3/2,-1/2),\\ \boldsymbol\delta_3&=a(-\sqrt3/2,-1/2),\\ \mathbf K&=(4\pi/(3\sqrt3a),0). \end{aligned}$$

First calculate all three phases at K.

$$\begin{aligned} \mathbf K\cdot\boldsymbol\delta_1&=0,\\ \mathbf K\cdot\boldsymbol\delta_2&=\frac{4\pi}{3\sqrt3a}\frac{\sqrt3a}{2}=\frac{2\pi}{3},\\ \mathbf K\cdot\boldsymbol\delta_3&=\frac{4\pi}{3\sqrt3a}\left(-\frac{\sqrt3a}{2}\right)=-\frac{2\pi}{3}. \end{aligned}$$
$$e^{i\mathbf K\cdot\boldsymbol\delta_1}=1,\qquad e^{i\mathbf K\cdot\boldsymbol\delta_2}=e^{i2\pi/3},\qquad e^{i\mathbf K\cdot\boldsymbol\delta_3}=e^{-i2\pi/3}$$

Next write the phase changes produced by $\mathbf q$.

$$\begin{aligned} \mathbf q\cdot\boldsymbol\delta_1&=aq_y,\\ \mathbf q\cdot\boldsymbol\delta_2&=a(\sqrt3q_x/2-q_y/2),\\ \mathbf q\cdot\boldsymbol\delta_3&=a(-\sqrt3q_x/2-q_y/2). \end{aligned}$$

Insert and into :

$$\begin{aligned}f(\mathbf K+\mathbf q)\simeq i\Big[&aq;_y\\ &+e^{i2\pi/3}a(\sqrt3q_x/2-q_y/2)\\ &+e^{-i2\pi/3}a(-\sqrt3q_x/2-q_y/2)\Big].\end{aligned}$$

Collect the $q_x$ terms. The difference of the two exponentials is $2i\sin(2\pi/3)$:

$$\begin{aligned} i\frac{\sqrt3a}{2}(e^{i2\pi/3}-e^{-i2\pi/3})q_x &=i\frac{\sqrt3a}{2}\,2i\sin(2\pi/3)q_x\\ &=i\frac{\sqrt3a}{2}(i\sqrt3)q_x\\ &=-\frac{3a}{2}q_x. \end{aligned}$$

All three bonds contribute to the $q_y$ terms. This time the sum of exponentials becomes $2\cos(2\pi/3)$.

$$\begin{aligned} ia\left[1-\frac12e^{i2\pi/3}-\frac12e^{-i2\pi/3}\right]q_y &=ia[1-\cos(2\pi/3)]q_y\\ &=ia(1+1/2)q_y\\ &=i\frac{3a}{2}q_y. \end{aligned}$$

Combining and gives

$$f(\mathbf K+\mathbf q)\simeq-\frac{3a}{2}q_x+i\frac{3a}{2}q_y=-\frac{3a}{2}(q_x-iq_y)\tag{5}$$

The component $q_x$ along Γ→K enters the real part, while the perpendicular component $q_y$ enters the imaginary part, with equal coefficient magnitudes $3a/2$. Taking the length of these perpendicular components in the complex plane gives

$$\begin{aligned}|f|^2&\simeq\left(\frac{3a}{2}\right)^2(q_x^2+q_y^2),\\ |f|&\simeq\frac{3a}{2}\sqrt{q_x^2+q_y^2}=\frac{3a}{2}|\mathbf q|. \end{aligned}\tag{6}$$

At a given distance from K, the magnitude is the same in every direction. The three bonds have produced this isotropic form at first order.

04

Two cones emerge

Substitute into .

$$E_\pm(\mathbf K+\mathbf q)\simeq\pm\frac{3ta}{2}|\mathbf q|\tag{7}$$

First take the same $q_y=0$ cut as in the opening figure. The upper energy is $(3ta/2)|q_x|$: moving away from K in either direction along this cut raises it in proportion to the distance from K. That gives a V with its bottom at K. The lower band is the upside-down V.

Points at the same distance from K have the same energy in every direction. Rotating the V-shaped cut around K produces two cones, the Dirac cones. We now work within this first-order description near K.

The base plane shows displacement from K: $q_x$ along Γ→K, and $q_y$ perpendicular to it. Height is energy. Drag to rotate.

The opposite corner K′ has the same cone. Replacing $\mathbf k$ by $-\mathbf k$ conjugates all three phase factors. This reverses the imaginary part of their sum $f$, while leaving its magnitude and the energy unchanged.

$$\begin{aligned}f(-\mathbf k)&=\sum_{j=1}^3 e^{-i\mathbf k\cdot\boldsymbol\delta_j}=f^*(\mathbf k),\\ E_\pm(-\mathbf k)&=E_\pm(\mathbf k).\end{aligned}$$

Return to K and compare displacements in different directions in the wavevector plane. The two TB equations below give the relative phase of the A/B coefficients. Here $u_A,u_B$ are the A and B $p_z$ orbital coefficients after separating out the position-dependent phase factors.

$$\begin{aligned}Eu_A&=-tf(\mathbf k)u_B,\\Eu_B&=-tf^*(\mathbf k)u_A.\end{aligned}$$

Insert :

$$\begin{aligned}Eu_A&=\frac{3ta}{2}(q_x-iq_y)u_B,\\Eu_B&=\frac{3ta}{2}(q_x+iq_y)u_A.\end{aligned}\tag{8}$$

On the upper band, first displace the wavevector outwards along the extension of Γ→K, the $+q_x$ direction. With $q_x=q>0,\ q_y=0$, the energy is $E=(3ta/2)q$, so reads

$$\begin{aligned}\frac{3ta}{2}q\,u_A&=\frac{3ta}{2}q\,u_B,\\u_B&=u_A.\end{aligned}$$

The A and B coefficients have equal magnitude and phase. Now keep $|\mathbf q|$ fixed and turn the displacement direction 90° counterclockwise to $+q_y$, so $q_x=0,\ q_y=q>0$. The energy stays the same, but a factor $-i$ enters the right-hand side.

$$\begin{aligned}\frac{3ta}{2}q\,u_A&=-i\frac{3ta}{2}q\,u_B,\\u_B&=i\,u_A.\end{aligned}$$

Multiplication by $i=e^{i\pi/2}$ advances the phase by $\pi/2$. The B coefficient now has a phase one quarter of a cycle ahead of A. The energy is unchanged, but the A/B phase relation is different.

To include other directions, use the outward extension of Γ→the selected K as the reference direction. Define $\theta$ as the angle through which $\mathbf q$ turns counterclockwise from this reference in the wavevector plane.

The circle joins displacements with equal $|\mathbf q|$, hence equal energy on one band. On the upper band, the reference direction ($\theta=0$) gives $u_B=u_A$; a quarter-turn counterclockwise ($\theta=\pi/2$) gives $u_B=i u_A$.
$$\begin{aligned}q_x&=|\mathbf q|\cos\theta,\qquad q_y=|\mathbf q|\sin\theta,\\ q_x-iq_y&=|\mathbf q|(\cos\theta-i\sin\theta)=|\mathbf q|e^{-i\theta}.\end{aligned}$$

Use in and solve for $u_B$. Away from the contact, where $|\mathbf q|>0$, the factor $|\mathbf q|$ in cancels the one in the denominator.

$$\begin{aligned}u_B&=\frac{E}{(3ta/2)|\mathbf q|}\,e^{i\theta}u_A\\&=\frac{\pm(3ta/2)|\mathbf q|}{(3ta/2)|\mathbf q|}\,e^{i\theta}u_A\\&=\pm e^{i\theta}u_A.\end{aligned}\tag{9}$$

The upper band takes the $+$ sign: multiplying $u_A$ by $e^{i\theta}$ preserves its magnitude and advances its phase by $\theta$. The lower band adds a factor $-1=e^{i\pi}$, increasing the phase difference by half a cycle. In this first-order approximation, the A/B coefficient ratio depends on the direction of $\mathbf q$. Holding that direction fixed while changing $|\mathbf q|$ changes the energy, but leaves the coefficient ratio unchanged.

Using Γ→the selected corner as the angular reference makes common to all three K corners of the hexagon. At K′, using Γ→K′ as the reference and applying gives $u_B=\pm e^{-i\theta}u_A$. As $\mathbf q$ turns counterclockwise from the reference, the B phase advances at K and retreats at K′.

05

Following the motion of the probability density

Let us follow the electron's motion by tracking where it is most likely to be found. Superposing states with nearby wavevectors creates regions of constructive and destructive interference. Start with two states and calculate how the peaks of their probability density move.

On the upper band, choose two nearby wavevectors along the outward extension of Γ→K. Their components along this reference direction are $K_x+q_1,K_x+q_2$ ($q_1,q_2>0$), and their perpendicular components are both zero. Write time as $\tau$. A state of energy $E_j$ carries the time factor $e^{-iE_j\tau/\hbar}$, so with $\omega_j=E_j/\hbar$ its phase is $(K_x+q_j)x-\omega_j\tau$.

is common to both states, as is the unit-magnitude phase factor $e^{iK_xx}$. Factoring these out leaves the variation due to displacement from K. For two states with equal weights, this is the sum $F$:

$$F(x,\tau)=e^{i\phi_1}+e^{i\phi_2},\qquad \phi_j=q_jx-\omega_j\tau\tag{10}$$

The probability density, averaged over the fine structure within a unit cell, is proportional to $|F|^2$. Multiplying by the complex conjugate gives

$$\begin{aligned}|F|^2 &=(e^{i\phi_1}+e^{i\phi_2})(e^{-i\phi_1}+e^{-i\phi_2})\\ &=2+e^{i(\phi_2-\phi_1)}+e^{-i(\phi_2-\phi_1)}\\ &=2+2\cos\!\big[(q_2-q_1)x-(\omega_2-\omega_1)\tau\big]\\ &=2+2\cos(\delta q\,x-\delta\omega\,\tau),\\ \delta q&=q_2-q_1,\qquad \delta\omega=\omega_2-\omega_1. \end{aligned}\tag{11}$$

A density maximum occurs wherever the cosine argument is $2\pi n$. Keeping the same integer $n$ follows the same maximum.

$$\delta q\,x_n-\delta\omega\,\tau=2\pi n \quad\Longrightarrow\quad x_n(\tau)=\frac{\delta\omega}{\delta q}\tau+\frac{2\pi n}{\delta q}$$

The coefficient of time, $\delta\omega/\delta q$, is the speed of the maximum. As the wavevectors approach each other, this ratio approaches the slope of the dispersion curve.

$$v_g=\lim_{\delta q\to0}\frac{\delta\omega}{\delta q} =\frac{d\omega}{dq}=\frac{1}{\hbar}\frac{dE}{dq}\tag{12}$$

This is the group velocity. Two waves give repeated maxima; many nearby wavevectors can form a localized wave packet. The same slope determines its motion.

Return to the upper cone. Along $+q_x$, where we selected the two states, is $E=(3ta/2)q$. Dividing its slope $3ta/2$ by $\hbar$ gives the wave-packet speed. Call this speed $v_F$.

$$\begin{aligned}v_F&=\frac{1}{\hbar}\frac{dE}{dq}=\frac{3ta}{2\hbar},\\E_\pm&=\pm\hbar v_F|\mathbf q|.\end{aligned}\tag{13}$$

Moving farther from K leaves the straight-line slope unchanged: different energies give the same speed. The Fermi velocity is the group velocity at the Fermi energy. In neutral graphene that energy lies at the contact, so $v_F$ denotes the limiting speed on approaching it.

For example, inserting $t=2.7\,\mathrm{eV}$ and $a=0.142\,\mathrm{nm}$ into gives

$$v_F=\frac{3(2.7\,\mathrm{eV})(0.142\times10^{-9}\,\mathrm m)}{2(6.582\times10^{-16}\,\mathrm{eV\,s})} \simeq8.7\times10^5\,\mathrm{m/s}$$

A state with wavevector $\mathbf k$ has crystal momentum $\hbar\mathbf k$. The therefore corresponds to a crystal-momentum displacement $\mathbf p=\hbar\mathbf q$. For a particle of mass $m$ with a parabolic dispersion $E=p^2/(2m)$, is $dE/dp=p/m$: a higher energy means a higher speed. On the upper cone, $E_+=v_Fp$ has the same slope at every energy, so the speed stays at $v_F$. The following figure lets us compare these two shapes.

Waves in each group
Dispersion
0 ℏ/t
Choose the central wavevectors $\bar q$ of groups ① and ②, and compare their motion from the same position after the same elapsed time. Dashed envelopes trace $\pm|F|$; the lower curves compare their probability densities. With different central wavevectors, the peaks stay together on the cone and separate on the parabola. Choosing gives the repeated maxima derived above.

Differentiating in two dimensions gives the following velocity. Here $\hat{\mathbf q}=\mathbf q/|\mathbf q|$ is a unit vector along $\mathbf q$.

$$\begin{aligned}(v_x,v_y) &=\frac1\hbar\left(\frac{\partial E_\pm}{\partial q_x},\frac{\partial E_\pm}{\partial q_y}\right)\\ &=\pm v_F\left(\frac{q_x}{|\mathbf q|},\frac{q_y}{|\mathbf q|}\right) =\pm v_F\hat{\mathbf q}\qquad(|\mathbf q|>0). \end{aligned}\tag{14}$$

The upper-band velocity points along $\mathbf q$, and the lower-band velocity points oppositely. Both have speed $v_F$. Around K, was also $\theta$ on the upper band and $\theta+\pi$ on the lower band. Measured from Γ→K, the velocity direction has this same angle.

The speed stays constant as the energy increases. Light in a vacuum shares this property: a photon's energy and momentum obey $E=cp$. has the same form, with the speed of light $c$ replaced by $v_F$. Near K, electrons in the crystal behave like massless particles.

We can read an electron's energy from the cone's height, and its velocity from the slope. The shape that emerged from three neighbouring bonds also tells us how the electron moves.

The height difference between two states tells us how much energy an electron must receive or give up to move between them. The same cone helps describe the exchanges of energy with light and lattice vibrations involved in graphene Raman scattering and phonon emission.

← Back to Learn