A first order linear system is any PDE of the form
$$q_t(x, t) + A q_x(x, t) = 0$$
Where $q(x,t):\mathbb{R}\times\mathbb{R} \to \mathbb{R}^n$ is a vector function, and $A$ is an $n\times n$ matrix satisfying the hyperbolicity condition, i.e. diagonalizable with real eigenvalues.
$q_t, q_x$ designate the derivatives with respect to time $t$ and space coordinate $x$, respectively.
We can easily see that the second order -classical- wave equation in one dimension $$u_{tt}= c^2u_{xx}$$ can be written as a first order linear system of two equations, simply by taking $v_t=cu_x\ ,\ v_x= \frac{1}{c}u_t$ ; \begin{align*} &v_{tx}= cu_{xx} \\ &v_{tx}= \frac{1}{c}u_{tt} \end{align*} To obtain $u_{tt}= c^2u_{xx}$ , so here the first order system is \begin{align*} &v_t- cu_x= 0 \\ &u_t- cv_x= 0 \end{align*} Meaning we take $q= (u,v):\ \mathbb{R}\times\mathbb{R} \to \mathbb{R}^2$ and $A= \Big(\begin{matrix} 0 & -c \\ -c & 0 \end{matrix}\Big)$
Now for the $3$-dimensional wave equation : $$u_{tt}= c^2\Delta u= c^2(u_{xx}+u_{yy}+u_{zz})$$ My question is : Can we turn this equation into a first order linear system or no ?
Maybe let's try first the 2D case: $u(x,y,t)$ and $u_{tt} = c^2\Delta (u_{xx} + u_{yy})$.
Then define $v^{(x)}_t = cu_x, v^{(y)}_t = cu_y$. For the spatial derivatives of $\boldsymbol v = \begin{pmatrix} v^{(x)} \\ v^{(y)} \end{pmatrix}$ require that $v^{(x)}_x + v^{(y)}_y = \frac 1c u_t$. Then for continuously differentiable $v^{(x)}, v^{(y)}$: \begin{align} \frac 1c u_{tt} &= \partial_t \Big(v^{(x)}_x + v^{(y)}_y \Big) \\ &= v^{(x)}_{xt} + v^{(y)}_{yt} \\ &= v^{(x)}_{tx} + v^{(y)}_{ty} \\ &= cu_{xx} + c u_{yy}\end{align} so that the original PDE (Wave equation) remains preserved. So the first order system reads $$ \partial_t\begin{pmatrix} v^{(x)} \\v^{(y)} \\ u \end{pmatrix} + \nabla \cdot \begin{pmatrix} -cu & 0 \\ 0 & -cu \\ -cv^{(x)} & -cv^{(y)} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0\end{pmatrix} $$ where the $\nabla \cdot$ acts row-wise.
This is the general form of a conservation law in multiple dimensions in divergence form: $$ \partial_t \boldsymbol u + \nabla \cdot \boldsymbol f(\boldsymbol u) = \boldsymbol{0}.$$
I am not sure how you would write a linear system in divergence form, since for a row-wise acting divergence $\nabla \cdot$ you need that $A \boldsymbol u \in \mathbb R^{m\times d}$, with $m$ the number of variables $| \boldsymbol u |$ and $d$ the spatial dimension. While you can ensure that $A$ has $m$ rows, there is no way for $A$ having $d$ columns if you multiply it with the column vector $\boldsymbol u \in \mathbb R^{m \times 1}$.
The extension to 3D is then straightforward, here you have $$ \partial_t\begin{pmatrix} v^{(x)} \\v^{(y)} \\ v^{(z)} \\ u \end{pmatrix} + \nabla \cdot \begin{pmatrix} -cu & 0 & 0\\ 0 & -cu &0 \\ 0 & 0 & -cu \\ -cv^{(x)} & -cv^{(y)} & -cv^{(z)} \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \\ 0 \\ 0\end{pmatrix} $$