Four Hermitian matrices, denoted by $ \alpha_{1} $, $ \alpha_{2} $, $ \alpha_{3} $ and $ \beta $, of dimension $ 4 \times 4 $ that satisfy the relations \begin{gather} \alpha_{k} \alpha_{j} + \alpha_{j} \alpha_{k} = 2 \delta_{k j} \mathsf{I}_{4}, \\ \alpha_{k} \beta + \beta \alpha_{k} = \mathbf{0}_{4}, \\ \alpha_{k} \alpha_{k} = \beta^{2} = \mathsf{I}_{4}, \end{gather} where $ \mathsf{I}_{4} $ is the $ (4 \times 4) $ identity matrix. The matrices $ \alpha_{1} $, $ \alpha_{2} $, $ \alpha_{3} $ and $ \beta $ may also be replaced by the Hermitian matrices $ \gamma^{k} = - i \beta \alpha_{k} $, where $ k \in \{ 1,2,3 \} $, and by the anti-Hermitian matrix $ \gamma^{0} = i \beta $. These then satisfy the relation $$ \gamma^{\alpha} \gamma^{\beta} + \gamma^{\beta} \gamma^{\alpha} = - 2 \eta^{\alpha \beta} \mathsf{I}_{4}, \qquad \forall \alpha,\beta \in \{ 0,1,2,3 \}. $$ Here, $ \eta^{\alpha \beta} \stackrel{\text{df}}{=} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & -1 & 0 & 0 \\ 0 & 0 & -1 & 0 \\ 0 & 0 & 0 & -1 \end{bmatrix} $. It is therefore possible to write the Dirac equation in a form that is covariant with respect to the Lorentz group of transformations. The matrices $ \alpha_{k} $, $ \beta $ and $ \gamma^{k} $, where $ k \in \{ 0,1,2,3 \} $, are defined up to an arbitrary unitary transformation and may be represented in various ways. One such representation is $$ \gamma^{0} = - i \begin{bmatrix} \mathsf{I}_{2} & \mathbf{0}_{2} \\ \mathbf{0}_{2} & - \mathsf{I}_{2} \end{bmatrix}; \qquad \gamma^{k} = - i \begin{bmatrix} \mathbf{0}_{2} & \boldsymbol{\sigma}_{k} \\ - \boldsymbol{\sigma}_{k} & \mathbf{0}_{2} \end{bmatrix}, $$ where the $ \boldsymbol{\sigma}_{k} $’s are the $ (2 \times 2) $ Pauli matrices, while $ \mathsf{I}_{2} $ and $ \mathbf{0}_{2} $ are the $ (2 \times 2) $ identity and zero matrices respectively. Dirac matrices may be used to factorize the Klein–Gordon equation in the following manner: $$ (\Box - m^{2}) E \psi = \left( \sum_{k = 0}^{3} \gamma^{k} \frac{\partial}{\partial x^{k}} - m E \right) \! \left( \sum_{l = 0}^{3} \gamma^{l} \frac{\partial}{\partial x^{l}} + m E \right) \psi = 0, $$ where $ \Box $ denotes the d’Alembert operator.
The Dirac matrices were Introduced by P. Dirac in 1928, in his derivation of the Dirac equation.
For references, see the article on the Dirac equation.