In calculating the electron dispersion you probably obtained the diagonalized Hamiltonian in the momentum space
$$
H=\sum_\mathbf{k}\left[c^{\dagger}_{\mathbf{k}A},c^{\dagger}_{\mathbf{k}B}\right]\left[\begin{array}{cc}0 & \Delta(\mathbf{k})\\ \Delta^{\dagger}(\mathbf{k}) &0\end{array}\right]\left[\begin{array}{c}c_{\mathbf{k}A} \\ c_{\mathbf{k}B}\end{array}\right].
$$
If you you chose your $x$ axis along the zigzag direction (arXiv:1004.3396), the two nonequivalent Dirac valleys are $\mathbf{K}_\kappa=\left(\kappa\frac{4\pi}{3\sqrt{3}a},0\right)$, $\kappa=\pm1$ and $\mathbf{K}_{-1}=\mathbf{K}^{\prime}$, where $a$ is the C-C distance. Then $\Delta(\mathbf{k})=-t\left(1+e^{-i\mathbf{k}\cdot\mathbf{a}_1}+e^{-i\mathbf{k}\cdot\mathbf{a}_2}\right)$, where $t$ is the hopping term, and $\mathbf{a}_1=\left(\sqrt{3}a/2,3a/2\right)$ and $\mathbf{a}_2=\left(-\sqrt{3}a/2,3a/2\right)$ are the lattice vectors.
Taylor expanding $\Delta(\mathbf{k})$ up to linear terms around those two points you obtain
$$
\Delta(\mathbf{k})=\kappa\frac{3ta}{2}q_x-i\frac{3ta}{2}q_y
$$
where $\mathbf{q}$ is the displacement momenta from the $\mathbf{K}_\kappa$ point. Promoting these displacement momenta to operators you obtain the Hamiltonian
$$
H=\hbar v_F\left[\begin{array}{cccc}0 & q_x-iq_y & 0 & 0\\q_x+iq_y & 0 & 0 & 0\\0 & 0 & 0 & -q_x-iq_y\\0 & 0 & -q_x+iq_y & 0\end{array}\right]
$$
where $v_F=\frac{3ta}{2\hbar}$ is the Fermi velocity. This is in $\left[\Psi_{A\mathbf{K}},\Psi_{B\mathbf{K}},\Psi_{A\mathbf{K}^{\prime}},\Psi_{B\mathbf{K}^{\prime}}\right]^T$ basis, if you rearrange your basis as $\left[\Psi_{A\mathbf{K}},\Psi_{B\mathbf{K}},\Psi_{B\mathbf{K}^{\prime}},\Psi_{A\mathbf{K}^{\prime}}\right]^T$ you get the compact form
$$
H=\hbar v_F\tau_z\otimes\boldsymbol{\sigma}\cdot\mathbf{k}
$$
where $\tau_z$ acts in the valley space. This is similar to the Dirac-Weyl equation for relativistic massless particles, where instead of $v_F$ you get the speed of light
$$
H=\pm\hbar c\boldsymbol{\sigma}\cdot\mathbf{k}
$$
where $+$ denotes right-handed antineutrions, and $-$ denotes left-handed neutrions. The differences are that $\boldsymbol{\sigma}=\left(\sigma_x,\sigma_y\right)$ for graphene acts in pseudospin space and $\boldsymbol{\sigma}=\left(\sigma_x,\sigma_y,\sigma_z\right)$ for neutrinos acts in real spin space.
Best Answer
Something seems to be off with the energy term here... Your last term goes as $\cos\theta$.
Let's try again :)
$$ H = \begin{pmatrix} 0&-t f_\mathbf{q} \\ -tf_\mathbf{q}^*&0 \end{pmatrix}\,, $$ where $f_\mathbf{q} = 1 + e^{i\mathbf{q}\cdot\mathbf{d}_1}+ e^{i\mathbf{q}\cdot\mathbf{d}_2}$ and $\mathbf{d}_1$ and $\mathbf{d}_2$ are the basis vectors. If the Dirac point is at $\mathbf{q} = \mathbf{K}$, for momenta close to it, we write
$$ f_\mathbf{q} = 1 + e^{i\left(\mathbf{q}+\mathbf{K}\right)\cdot\mathbf{d}_1}+ e^{i\left(\mathbf{q}+\mathbf{K}\right)\cdot\mathbf{d}_2} = 1 + e^{i\mathbf{k}\cdot\mathbf{d}_1} e^{i\mathbf{K}\cdot\mathbf{d}_1} + e^{i\mathbf{k}\cdot\mathbf{d}_2} e^{i\mathbf{K}\cdot\mathbf{d}_2}\,. $$
For small $k$, we get $$ f_\mathbf{q} \approx 1 + (1+i\mathbf{k}\cdot\mathbf{d}_1 - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_1\right)^2) e^{i\mathbf{K}\cdot\mathbf{d}_1} + (1+i\mathbf{k}\cdot\mathbf{d}_2 - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_2\right)^2) e^{i\mathbf{K}\cdot\mathbf{d}_2} \\ = 1 + e^{i\mathbf{K}\cdot\mathbf{d}_1} + e^{i\mathbf{K}\cdot\mathbf{d}_2} + (i\mathbf{k}\cdot\mathbf{d}_1 ) e^{i\mathbf{K}\cdot\mathbf{d}_1} + (i\mathbf{k}\cdot\mathbf{d}_2) e^{i\mathbf{K}\cdot\mathbf{d}_2} - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_1\right)^2 e^{i\mathbf{K}\cdot\mathbf{d}_1} - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_2\right)^2 e^{i\mathbf{K}\cdot\mathbf{d}_2} \,. $$
Using $\mathbf{K} = 4\pi(0, 1)/(3\sqrt{3}a)$ and $\mathbf{d}_{1/2}=a(3,\pm\sqrt{3})/2$, we have $e^{i\mathbf{K}\cdot\mathbf{d}_1} = e^{ 2\pi i/3}$ and $e^{i\mathbf{K}\cdot\mathbf{d}_2} = e^{- 2\pi i/3}$:
$$ f_\mathbf{q} = i\mathbf{k}\cdot(\mathbf{d}_1 e^{ 2\pi i/3} + \mathbf{d}_2 e^{- 2\pi i/3}) - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_1\right)^2 e^{2\pi i/3} - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_2\right)^2 e^{- 2\pi i/3} \\ = -\frac{3a}{2}ik e^{-i \theta} - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_1\right)^2 e^{ 2\pi i/3} - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_2\right)^2 e^{-2\pi i/3} \,, $$ where $\mathbf{k} = k(\cos\theta,\sin\theta)$.
We know that the dispersion is $\pm t^2\sqrt{f_\mathbf{q}f_\mathbf{q}^*}$, so let's deal with the square root. For small $k$ (keeping to 3d order), we have
$$ f_\mathbf{q}f_\mathbf{q}^* \approx \frac{9a^2}{4}k^2 -\frac{3a}{2}ik e^{-i \theta} \left[ - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_1\right)^2 e^{ -2\pi i/3} - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_2\right)^2 e^{2\pi i/3} \right] +\frac{3a}{2}ik e^{i \theta} \left[ - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_1\right)^2 e^{ 2\pi i/3} - \frac{1}{2}\left(\mathbf{k}\cdot\mathbf{d}_2\right)^2 e^{-2\pi i/3} \right] \\ =\frac{9}{4}a^2k^2+\frac{9}{8}a^3k^3\sin(3\theta)\,. $$ You can check this in, say, Mathematica.
From this, $$ \sqrt{f_\mathbf{q}f_\mathbf{q}^*} = \sqrt{\frac{9}{4}a^2k^2+\frac{9}{8}a^3k^3\sin(3\theta)} = \frac{3}{2}ak\sqrt{1+\frac{1}{2}a^2k^2\sin(3\theta)} \approx \frac{3}{2}ak\left(1+\frac{1}{4}a^2k^2\sin(3\theta)\right)\,, $$ as required.