What is the general form of the Pauli matrices for a $3\times 3$ matrix (spin $s=1$)? Do these generalized Pauli matrices satisfy all the properties exhibited by the Pauli $2\times 2$ matrices, such as being Hermitian? Furthermore, is it possible to represent all Hermitian $3\times 3$ matrices using these generalized Pauli matrices?
As @Charlie asks, the Pauli matrices have several properties, which generalize in different manners. They certainly, together with the identity, provide a complete basis for 2×2 matrices, but they are also hermitian. If Hermiticity is important to you, you generalize them as in the link provided, that is along the Gell-Mann matrices' route for 3×3 matrices.
However, a far more tasteful and systematic basis is J J Sylvester's 1882 one of clock and shift matrices for d×d unitary matrices which you should know about, anyway. They are not hermitean in general, but they are more systematic (some would say "analytic in d").
For $ω= \exp(2iπ/d)$, a root of unity not equal to 1. The sum of all roots annuls, $1 + \omega + \cdots + \omega ^{d-1} = 0 $, so integer indices may be cyclically identified mod d.
The shift matrix is defined as $$ \Sigma _1 = \begin{bmatrix} 0 & 0 & 0 & \cdots & 0 & 1\\ 1 & 0 & 0 & \cdots & 0 & 0\\ 0 & 1 & 0 & \cdots & 0 & 0\\ 0 & 0 & 1 & \cdots & 0 & 0\\ \vdots & \vdots & \vdots & \ddots &\vdots &\vdots\\ 0 & 0 & 0 & \cdots & 1 & 0\\ \end{bmatrix} $$ and the clock matrix as $$ \Sigma _3 = \begin{bmatrix} 1 & 0 & 0 & \cdots & 0\\ 0 & \omega & 0 & \cdots & 0\\ 0 & 0 & \omega^2 & \cdots & 0\\ \vdots & \vdots & \vdots & \ddots & \vdots\\ 0 & 0 & 0 & \cdots & \omega^{d-1} \end{bmatrix}. $$ These matrices generalize $σ_1$ and the diagonal $σ_3$, respectively. Since Pauli matrices describe quaternions, Sylvester dubbed the higher-dimensional analogs "nonions", "sedenions", etc.
These two matrices are also the cornerstone of Weyl's celebrated quantum mechanical dynamics in finite-dimensional vector spaces The clock matrix amounts to the exponential of position in a "clock" of d hours, and the shift matrix is just the translation operator in that cyclic vector space, so the exponential of the momentum. They are (finite-dimensional) representations of the corresponding elements of the Heisenberg group on a d -dimensional Hilbert space.
The following relations echo and generalize those of the Pauli matrices:
$\Sigma_1^d = \Sigma_3^d = I$,
and the braiding relation,
$\Sigma_3 \Sigma_1 = \omega \Sigma_1 \Sigma_3 = e^{2\pi i / d} \Sigma_1 \Sigma_3$, and can be rewritten as
$\Sigma_3 \Sigma_1 \Sigma_3^{d-1} \Sigma_1^{d-1} = \omega ~$.
The complete family of $d^2$ unitary (but non-Hermitian) independent matrices $$ \left(\Sigma_1\right)^k \left(\Sigma_3\right)^j = \sum_{m=0}^{d-1} |m+k\rangle \omega^{jm} \langle m|, $$ then provides Sylvester's well-known trace-orthogonal basis for $\mathfrak{gl} (d,ℂ)$, known as "nonions" $\mathfrak{gl} (3,ℂ)$, "sedenions" $\mathfrak{gl} (4,ℂ)$, etc... Since all indices are defined cyclically mod d, $\mathrm{tr}\Sigma_1^j \Sigma_3^k \Sigma_1^m \Sigma_3^n = \omega^{km} d ~\delta_{j+m,0} \delta_{k+n,0}$
Pauli matrices (plus the identity matrix) are just a choice of matrices that allow decomposition of an arbitrary 2-by-2 matrix - i.e. a matrix with 4 independent parameters. One could choose them differently, so this particular choice is more due to the tradition and the fact that all the three matrices are already Hermitian. It is not uncommon to use $\hat{\sigma}_\pm$ instead of $\sigma_{x,y}$.
For 3-by-3 case one needs in principle 9 matrices, one of which can be the identity matrix. They also can be all chosen Hermitian. However, one will have more freedom than in the 2-by-2 case. This freedom might be further restricted by a particular application - e.g., describing certain type of rotations.
