Dirac Cone in Graphene

When a roller coaster slides down from a high place, potential energy is converted into kinetic energy, increasing its speed. The higher the energy, the faster it goes, this intuition is familiar from throwing a ball or pedaling a bicycle.

However, near the K point in graphene, electrons with different energies move at nearly constant speed. Why does this phenomenon occur, defying our everyday mechanical intuition? To unravel its essence, we need to treat electrons as "waves." Based on the band structure derived in the previous section, let us examine how electron waves propagate.

01

Looking at Electrons Near the Contact Point

In charge-neutral graphene, at absolute zero ($T=0$), the lower band (valence band) is completely filled with electrons, while the upper band (conduction band) is empty. The Fermi energy (Fermi level), which separates the two, is located at the contact point between the upper and lower bands at the K and K′ points. Therefore, only electrons located near this contact point can gain a small amount of energy and transition to empty states.

In the tight-binding approximation considering only interactions (hopping between $p_z$ orbitals) with the three nearest-neighbor carbon atoms, the energy dispersion of the two bands was expressed as follows:

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

where $t>0$ is the hopping parameter between adjacent orbitals, and $f(\mathbf k)$ is the sum of phase factors in three directions.

At the K point, these three phase factors cancel each other out, yielding $f(\mathbf K)=0$, so the energies of both bands become zero ($E_\pm(\mathbf K)=0$). Looking closely near the K point, we see that the energy dispersion extending from the contact point approaches a straight line with respect to changes in the wavevector.

×3.0
The colored lines in the hexagon indicate the cross-sectional directions. The lower panel expands the framed area in the upper diagram at equal magnification vertically and horizontally. The variable $a$ in $ak_x$ on the horizontal axis represents the carbon-carbon bond length.
02

Describing States via Displacement from the K Point

In §01, we looked at the band cross-section cut along the $\Gamma \to \text{K}$ direction. Now, let us broaden our view and consider displacing in an arbitrary direction from the K point within the wavevector plane. Writing the displacement from the K point as $\mathbf q = (q_x, q_y)$, the total wavevector $\mathbf k$ measured from the $\Gamma$ point is expressed as follows:

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

Here, the direction pointing from the $\Gamma$ point toward the K point of interest is defined as the $+q_x$ axis, and the direction rotated 90° counterclockwise from it within the wavevector plane is defined as the $+q_y$ axis. $q_x$ and $q_y$ are the components of the displacement from the K point along these two axes. The cross-section in §01 corresponds to the line $q_y=0$, and $\mathbf q=0$ represents the contact point of the upper and lower bands (the K point) itself.

03

Breakdown of Phase Factor Cancellation and Linear Approximation

Moving slightly away from the K point breaks the perfect cancellation of the phase factors from the three directions. Since the remaining components of these phase factors determine the energy dispersion, let us first calculate their sum $f(\mathbf k)$. Defining the bond vectors pointing from an A site to its three adjacent B sites as $\boldsymbol\delta_1, \boldsymbol\delta_2, \boldsymbol\delta_3$, the sum of the phase factors is defined as follows:

$$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}$$

Substituting $\mathbf k = \mathbf K + \mathbf q$ from Eq. (2), we decompose each exponential function into the "phase at the K point" and the "change from the K point."

$$\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}$$

Since we consider only the immediate vicinity of the K point, the displacement $\mathbf q$ is sufficiently small. The phase change $\mathbf q\cdot\boldsymbol\delta_j$ is also small, so we apply the linear approximation ($e^{ix} \simeq 1 + ix$), keeping terms up to first order.

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

Substituting this, we separate terms without and with $\mathbf q$.

$$\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 sum in the first term vanishes due to complete cancellation at the K point. Therefore, what remains is the sum of the changes in the three contributions caused by the displacement $\mathbf q$.

$$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 component calculation, let us use the coordinate system set up in the previous section (with carbon-carbon distance $a$):

$$\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, calculating the three phase factors at the K point yields:

$$\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, we write down the phase changes due to the displacement $\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}$$

Substituting these into Eq. (4) gives:

$$\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}$$

Grouping the terms with $q_x$, the difference between the two exponentials is $e^{i2\pi/3}-e^{-i2\pi/3} = 2i\sin(2\pi/3) = i\sqrt{3}$, so:

$$\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}$$

Similarly, all three bonds enter the terms for $q_y$. Grouping the sum of exponentials as $e^{i2\pi/3}+e^{-i2\pi/3} = 2\cos(2\pi/3) = -1$, we get:

$$\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 the $q_x$ and $q_y$ terms:

$$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 the $\Gamma \to \text{K}$ direction appears in the real part, while the orthogonal component $q_y$ appears in the imaginary part, with both coefficients having a magnitude of $3a/2$. Taking the magnitude of two orthogonal components in the complex plane:

$$\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}$$

This result is of paramount importance. As long as the distance $|\mathbf q|$ from the K point is the same, the magnitude of $|f|$ remains constant regardless of the displacement direction. An isotropic form independent of direction emerges in the linear approximation from the anisotropic structure of three-directional bonds in the crystal lattice.

04

Appearance of Two Cones and Phase Relationship Between Sublattices

Let us substitute Eq. (6) back into the original energy dispersion equation (1):

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

First, consider the $q_y=0$ cross-section as in §01. The upper band energy becomes $E_+ = (3ta/2)|q_x|$. Moving away from the K point along this cross-section in either direction, the energy increases in proportion to the distance from the K point. This is why a V-shaped dispersion is formed with the K point at the bottom (apex). The lower band takes the inverted form of this V-shape.

Furthermore, points at equal distances from the K point have identical energies regardless of direction. Thus, rotating the V-shaped cross-section in all directions around the K point forms two symmetric cones. This is the Dirac cone. Below, we treat the electronic states near the K point within this linear approximation.

The base represents the displacement $\mathbf q$ from the K point ($q_x$: $\Gamma \to \text{K}$ direction, $q_y$: orthogonal direction), and the height represents energy. (Drag to rotate in 3D)

Note that an identical cone is formed near the K′ point located on the opposite side across the $\Gamma$ point. Inverting the wavevector from $\mathbf k$ to $-\mathbf k$ takes the complex conjugate of each of the three phase factors, so $f(-\mathbf k) = f^*(\mathbf k)$. Although the sign of the imaginary part is inverted, neither its absolute value $|f|$ nor the energy changes.

$$\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}$$

Returning to the K point, let us examine how shifting the wavevector in different directions changes the phases of the wavefunctions (coefficients) on the A and B sublattices. This relationship is preserved in the following two equations obtained from the tight-binding approximation (where $u_A, u_B$ are expansion coefficients of the $p_z$ orbitals at sublattices A and B after factoring out the phase factors for each atomic position):

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

Substituting the approximation $f(\mathbf K+\mathbf q) \simeq -\frac{3a}{2}(q_x-iq_y)$ from Eq. (5) yields:

$$\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}$$

In the upper band ($E>0$), let us first shift the wavevector in the direction extending $\Gamma \to \text{K}$ outward. This corresponds to the $+q_x$ direction ($q_x=q>0, q_y=0$). The energy is $E=(3ta/2)q$, so the first equation in Eq. (8) becomes:

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

indicating that the coefficients $u_A$ and $u_B$ match completely in both magnitude and phase.

What if we rotate the displacement direction by 90° counterclockwise while keeping the distance $|\mathbf q|=q$? Now $q_x=0, q_y=q>0$. Even though the energy is the same, $-i$ appears on 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}$$

Since $i=e^{i\pi/2}$, the phase advances by $\pi/2$ (a quarter cycle). That is, the coefficient for B now has a phase advanced by a quarter cycle compared to A. Even at the same energy, distinct characteristics remain in the phase relationship between A and B depending on the displacement direction.

Let us summarize for other directions. Take the direction pointing from the $\Gamma$ point toward the K point of interest and extending outward as the reference axis ($\theta=0$). In the wavevector plane, let $\theta$ be the angle measured counterclockwise from this reference to the displacement vector $\mathbf q$.

Points on the circle have equal $|\mathbf q|$ and identical energies within the same band. In the upper band, $u_B=u_A$ along the reference direction ($\theta=0$), and $u_B=i u_A$ at $\theta=\pi/2$.

Using polar coordinates $q_x=|\mathbf q|\cos\theta, q_y=|\mathbf q|\sin\theta$:

$$\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}$$

Substituting this into Eq. (8) and solving for $u_B$: For points away from the contact point ($|\mathbf q|>0$), $|\mathbf q|$ in the numerator cancels with $|\mathbf q|$ 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}$$

Selecting the upper band ($+$) gives $u_B = e^{i\theta}u_A$, meaning the B coefficient is equal to the A coefficient multiplied by $e^{i\theta}$. This operation leaves the magnitude untouched while advancing the phase by $\theta$. In the lower band ($-$), an additional factor of $-1=e^{i\pi}$ is multiplied, increasing the phase difference by half a cycle ($\pi$).

In this linear approximation, the ratio of coefficients (phase relationship) between sublattices A and B is determined solely by the direction $\theta$ of the displacement $\mathbf q$ from the K point. Holding that direction fixed and changing the distance $|\mathbf q|$ alters the energy, but does not change the phase relationship between A and B.

Taking the direction from the $\Gamma$ point toward the K point of interest as the reference, Eq. (9) holds commonly for the three K points of the hexagon. For the K′ point, taking $\Gamma \to \text{K}'$ as the reference yields $u_B=\pm e^{-i\theta}u_A$ from the corresponding equation to Eq. (8). Rotating $\mathbf q$ counterclockwise from the reference advances the phase of B at the K point, while delaying it at the K′ point, a symmetric inversion.

05

Motion of Electron Wavepackets and Constant Velocity

Let us track the propagation of an electron as the movement of the position where the electron is most likely to be found (peaks of probability density). Superimposing states with slightly different wavevectors produces regions of constructive and destructive interference. Let us first superimpose two states and calculate how the peak of probability density moves.

In the upper band, select two nearby wavevectors along the direction extending $\Gamma \to \text{K}$ outward. Set these reference components to $K_x+q_1, K_x+q_2$ ($q_1,q_2>0$) and orthogonal components to 0. Writing time as $\tau$ (to avoid confusion with the hopping parameter $t$), a state with energy $E_j$ evolves in time with a phase factor $e^{-iE_j\tau/\hbar}$. Setting $\omega_j=E_j/\hbar$, the phase of each wave becomes $(K_x+q_j)x-\omega_j\tau$.

Factoring out the common phase factor $e^{iK_x x}$ of unit magnitude, we focus on the part affected by the displacement from the K point. Superimposing the two waves with equal amplitude gives the composite wave $F(x,\tau)$:

$$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 fine atomic-scale variations, is proportional to $|F|^2$. Multiplying by its complex conjugate yields:

$$\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}$$

Peaks of probability density occur where the cosine argument equals $2\pi n$ ($n$ is an integer). Tracking the same peak requires fixing the integer $n$:

$$\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 multiplying time, $\delta\omega/\delta q$, represents the propagation velocity of the peak. Taking the limit as the two wavevectors approach each other, this ratio approaches the derivative (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. While superimposing two waves produces periodically arranged peaks, superimposing many close wavevectors forms a localized wavepacket. The propagation velocity of the entire wavepacket is also determined by this slope.

Returning to the upper cone, along the $+q_x$ direction chosen above, the energy dispersion from Eq. (7) is $E=(3ta/2)q$. Dividing its slope $3ta/2$ by $\hbar$ gives the wavepacket speed. Let us define this speed as the Fermi velocity $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}$$

The slope of the straight line is constant regardless of the distance $q$ from the K point, so the speed remains constant regardless of energy. Although Fermi velocity strictly refers to the group velocity at the Fermi energy, in neutral graphene where the Fermi level lies at the contact point (Dirac point), the constant speed $v_F$ as one approaches the contact point is referred to as such.

For instance, substituting $t=2.7\,\mathrm{eV}$ and $a=0.142\,\mathrm{nm}$ into Eq. (13) yields:

$$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}$$

which is an extraordinarily high speed, approximately 1/300th the speed of light.

A state with wavevector $\mathbf k$ carries a crystal momentum of $\hbar\mathbf k$. The shift in crystal momentum corresponding to the displacement $\mathbf q$ from the K point is $\mathbf p=\hbar\mathbf q$. For a free particle of mass $m$ with parabolic dispersion $E=p^2/(2m)$, the group velocity is $v_g = dE/dp=p/m$, which increases with energy. In contrast, on the upper cone where $E_+=v_Fp$, changing the height does not change the slope, maintaining a constant speed $v_F$.

Waves in each group
Dispersion
0 ℏ/t
Select central wavevector $\bar q$ for ① and ②, moving them from the same starting position over the same duration. Dashed lines indicate the envelope $\pm|F|$; the lower panel shows probability density. When central wavevectors differ, peaks advance in unison on the cone, but drift apart on the parabola. Display modes toggle the interference wave structure.

Differentiating the cone dispersion equation (13) with respect to the 2D wavevector yields the velocity vector ($\hat{\mathbf q}=\mathbf q/|\mathbf q|$ is the unit vector in the direction of $\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}$$

Electrons move in the direction of $\mathbf q$ in the upper band and in the opposite direction in the lower band, but the speed is equal to $v_F$ in both cases. Furthermore, around the K point, the phase difference between sublattices derived in Eq. (9) was $\theta$ for the upper band and $\theta+\pi$ for the lower band. Measured from $\Gamma \to \text{K}$, the phase difference of the wavefunction perfectly aligns with the direction angle of the velocity vector.

The property that speed does not change even as energy increases is shared with photons in vacuum. The relationship between energy and momentum for photons is $E=cp$, and Eq. (13) derived above takes the exact same form with the speed of light $c$ replaced by $v_F$. That is, near the K point, electrons in the crystal behave like "massless particles (massless Dirac particles)".

Looking at the shape of the cone, one can intuitively read the electron energy from the height and the propagation speed from the slope. The geometrical shape derived from the bonds with three adjacent carbon atoms also reflects the microscopic motion of electrons.

Furthermore, examining the difference between two states on the cone allows a quantitative and visual understanding of how much energy an electron gains or loses during transitions. This Dirac cone plays a fundamental role in unraveling energy exchanges with light and lattice vibrations in processes such as Raman scattering and phonon emission in graphene.

← Back to Learn