Dirac equation from scratch
These are some notes I took while reading Modern Particle Physics by Mark Thomson.
Probability current and density
The change in probability density $\rho=\psi^*\psi$ of a volume $V$ over time equals to the total probability leaving the surface enclosing the volume, giving us equation
\[\begin{equation} \dfrac{\partial}{\partial t}\int_V\rho dV=-\int_S\mathbf{j}\cdot d\mathbf{S}, \end{equation}\]where $\mathbf{j}$ is the probability flux density and $d\mathbf{S}$ is an elemental surface. By using the divergence theorem we obtain
\[\begin{equation} \int_V\dfrac{\partial\rho}{\partial t} dV=-\int_V\nabla\cdot\mathbf{j} dV. \end{equation}\]The relationship between the probability density and the probability flux density is given by equation
\[\begin{equation} \nabla\cdot\mathbf{j}+\dfrac{\partial\rho}{\partial t}=0. \end{equation}\]The Klein-Gordon equation
The Klein-Gordon equation ($\mathrm{KGE}$) relates the energy of a free particle with its momentum according to equation
\[\begin{equation} E^2=\mathbf{p}^2+m^2, \end{equation}\]which can also be expressed in terms of operators with equation
\[\begin{equation} \hat{E}^2\psi=(\hat{p}^2+m^2)\psi \end{equation}\]For wavefunctions of form $Ne^{i(\mathbf{p}\cdot\mathbf{x}-Et)}$ we have energy and momentum operators $\hat{E}=i\dfrac{\partial}{\partial t}$ and $\hat{p}=-i\nabla$. The Klein-Gordon equation expressed with these operations is
\[\begin{equation} \dfrac{\partial^2\psi}{\partial t^2}=\nabla^2\psi-m^2\psi. \end{equation}\]By taking $\psi^*\times\mathrm{KGE}-\psi\times\mathrm{KGE}^*$ we obtain
\[\begin{equation} \begin{split} & \psi^*\dfrac{\partial^2\psi}{\partial t^2}-\psi\dfrac{\partial^2\psi^*}{\partial t^2}=\psi^*\left(\nabla^2\psi-m^2\psi\right)-\psi\left(\nabla^2\psi^*-m^2\psi^*\right) \newline \Rightarrow & \dfrac{\partial}{\partial t}\left(\psi^*\dfrac{\partial\psi}{\partial t}-\psi\dfrac{\partial\psi^*}{\partial t}\right)=\nabla\cdot\left(\psi^*\nabla\psi-\psi\nabla\psi^*\right) \end{split}, \end{equation}\]which in turn gives us probability density $\rho=i\left(\psi^*\dfrac{\partial\psi}{\partial t}-\psi\dfrac{\partial\psi^*}{\partial t}\right)$ and probability flux density $\mathbf{j}=-i\left(\psi^*\nabla\psi-\psi\nabla\psi^*\right)$. The factor of $i$ is included to ensure that the probability density is real. For wavefunctions of form $Ne^{i(\mathbf{p}\cdot\mathbf{x}-Et)}$ we get
\[\begin{equation} \rho=2|N|^2E~\mathrm{and}~\mathbf{j}=2|N|^2\mathbf{p}. \end{equation}\]Notice that according to this equation negative energy solutions have negative probability densities, which is unphysical. A more sophisticated equation is needed.
The Dirac equation
To tackle the problems with the Klein-Gordon equation, Dirac formulated the now-famous
Dirac equation ($\mathrm{DE}$)
Expressed with operators $\hat{E}=i\dfrac{\partial}{\partial t}$ and $\mathbf{\hat{p}}=-i\nabla$, the equation takes form
\[\begin{equation} i\dfrac{\partial\psi}{\partial t}=-i\alpha_x\dfrac{\partial \psi}{\partial x}-i\alpha_y\dfrac{\partial \psi}{\partial y}-i\alpha_z\dfrac{\partial \psi}{\partial z}+\beta m\psi. \end{equation}\]From the equation above we form equation
\[\begin{equation} \begin{split} -\dfrac{\partial^2\psi}{\partial t^2} = &~\left(i\alpha_x\dfrac{\partial}{\partial x}+i\alpha_y\dfrac{\partial}{\partial y}+i\alpha_z\dfrac{\partial}{\partial z}-\beta m\right)\left(i\alpha_x\dfrac{\partial}{\partial x}+i\alpha_y\dfrac{\partial}{\partial y}+i\alpha_z\dfrac{\partial}{\partial z}-\beta m\right)\psi \newline \Rightarrow\dfrac{\partial^2\psi}{\partial t^2} = &~\alpha_x^2\dfrac{\partial^2\psi}{\partial x^2}+\alpha_y^2\dfrac{\partial^2\psi}{\partial y^2}+\alpha_z^2\dfrac{\partial^2\psi}{\partial z^2}-\beta^2m\psi \newline &+(\alpha_x\alpha_y-\alpha_y\alpha_x)\dfrac{\partial^2\psi}{\partial x\partial y}+(\alpha_y\alpha_z-\alpha_z\alpha_y)\dfrac{\partial^2\psi}{\partial y\partial z}+(\alpha_x\alpha_z-\alpha_z\alpha_x)\dfrac{\partial^2\psi}{\partial x\partial z} \newline &+i\left(\alpha_x\beta+\beta\alpha_x\right)m\dfrac{\partial\psi}{\partial x}+i\left(\alpha_y\beta+\beta\alpha_y\right)m\dfrac{\partial\psi}{\partial y}+i\left(\alpha_z\beta+\beta\alpha_z\right)m\dfrac{\partial\psi}{\partial z} \end{split}. \end{equation}\]In order for this equation to satisfy the Klein-Gordon equation
\[\begin{equation} \dfrac{\partial^2\psi}{\partial t^2}=\nabla^2\psi-m^2\psi, \end{equation}\]the following relations
\[\begin{equation} \begin{split} \alpha_x^2=\alpha_y^2=\alpha_z^2=\beta^2 & =I\newline \alpha_j\beta+\beta\alpha_j & =0\newline \alpha_j\alpha_k+\alpha_k\alpha_j & =0,~j\neq k \end{split} \end{equation}\]must be satisfied. The second and the third relations are anticommutative and can not be satisfied if $\alpha_i$ and $\beta$ are normal numbers. The simplest mathematical objects satisfying the relations are matrices. From the cyclic property of traces, $\mathrm{Tr}\left(ABC\right)=\mathrm{Tr}\left(BCA\right)$, we obtain
\[\begin{equation} \mathrm{Tr}(\alpha_i)=\mathrm{Tr}(\alpha_i\beta\beta)=\mathrm{Tr}(\beta\alpha_i\beta)=-\mathrm{Tr}(\alpha_i\beta\beta)=-\mathrm{Tr}(\alpha_i), \end{equation}\]which implies that the $\alpha_i$ and $\beta$ matrices must have trace zero. The eigenvalues $\lambda$ of the matrices must satisfy
\[\begin{equation} \alpha_iX=\lambda X\Rightarrow X=\alpha_i^2X=\lambda^2X\Rightarrow\lambda^2=\pm1. \end{equation}\]Since trace of a matrix is the sum of the eiganvalues and in the case of matrices $\alpha_i$ and $\beta$ has to equal to zero, the $\alpha_i$ and $\beta$ matrices must be of even dimension.
In order for the Dirac Hamiltonian operator $\hat{H}_D$ to have only real eigenvalues, the $\alpha_i$ and $\beta$ matrices must be Hermitian \(\begin{equation} \alpha_i=\alpha_i^\dagger~\mathrm{and}~\beta=\beta^\dagger. \end{equation}\)
Because there can only be three mutually anticommuting $2\times2$ traceless matrices, the lowest dimension object that can represent $\alpha_x$, $\alpha_y$, $\alpha_z$ and $\beta$ are $4\times4$ matrices and thus $\hat{H}_D$ is a $4\times4$ matrix of operators acting on a four-component wavefunction, known as Dirac spinor,
\[\begin{equation} \psi=\begin{pmatrix}\psi_1 \\ \psi_2 \\ \psi_3 \\ \psi_4 \end{pmatrix}. \end{equation}\]One of the commonly used representations for the matrices is the Dirac-Pauli representation, where
\[\begin{equation} \beta=\begin{pmatrix}I & 0 \\ 0 & I\end{pmatrix}~\mathrm{and}~\alpha_i=\begin{pmatrix}0 & \sigma_i \\ \sigma_i & 0\end{pmatrix}, \end{equation}\]with
\[\begin{equation} I=\begin{pmatrix}1 & 0 \\ 0 & 1\end{pmatrix},~\sigma_x=\begin{pmatrix}0 & 1 \\ 1 & 0\end{pmatrix},~\sigma_y=\begin{pmatrix}0 & -i \\ i & 0\end{pmatrix}~\mathrm{and}~\sigma_z=\begin{pmatrix}1 & 0 \\ 0 & -1\end{pmatrix}. \end{equation}\]To obtain an equations for probability density we take $\psi^\dagger\times\mathrm{DE}-\psi\times\mathrm{DE}^\dagger$ and get
\[\begin{equation} \begin{split} i\psi^\dagger\dfrac{\partial\psi}{\partial t}+i\psi\dfrac{\partial\psi^\dagger}{\partial t}=&~\psi^\dagger\left(-i\alpha_x\dfrac{\partial\psi}{\partial x}-i\alpha_y\dfrac{\partial\psi}{\partial y}-i\alpha_z\dfrac{\partial\psi}{\partial z}+\beta m\psi\right) \newline & -\left(i\alpha_x\dfrac{\partial\psi}{\partial x}+i\alpha_y\dfrac{\partial\psi}{\partial y}+i\alpha_z\dfrac{\partial\psi}{\partial z}+\beta m\psi\right)\psi^\dagger \end{split} \end{equation}\]By writing $\psi^\dagger\alpha_x\dfrac{\partial\psi}{\partial x}+\dfrac{\partial\psi^\dagger}{\partial x}\alpha_x\psi\equiv\dfrac{\partial\left(\psi^\dagger\alpha\dagger\right)}{\partial x}$ and $\psi^\dagger\dfrac{\partial\psi}{\partial t}+\psi\dfrac{\partial\psi^\dagger}{\partial t}=\dfrac{\partial\left(\psi^\dagger\psi\right)}{\partial t}$ the equation becomes
\[\begin{equation} \nabla\cdot\left(\psi^\dagger\mathbf{\alpha}\psi\right)+\dfrac{\partial\left(\psi^\dagger\psi\right)}{\partial t}=0, \end{equation}\]where $\psi^\dagger=\left(\psi_1^\dagger,\psi_2^\dagger,\psi_3^\dagger,\psi_4^\dagger\right)$. The probability density is
\[\rho=|\psi_1|^2+|\psi_2|^2+|\psi_3|^2+|\psi_4|^2,\]which indicates that all the solutions of the Dirac equation have positive probability density, giving physical solutions for the Klein-Gordon equation.