Can the hidden state of an RNN be a matrix?
Yes, I would say more, that hidden state can be a tensor of arbitrary dimensionality. For vanilla RNN the update rules of the hidden state and the output are:
$$
h_t = \sigma_h(W_h x_t + U_h h_{t-1} + b_h) 
$$
$$
y_t = \sigma_y(W_y h_t + b_y)
$$
Here $W_h$ is input to hidden state matrix, and $U_h$ is hidden state to hidden state, $W_y$ is hidden state to output. One can simply upgrade these matrix products to more general case of tensor products, such that they could handle multidimensional input data and hidden state.
The contraction of the $m$-dimensional tensor $X$ with $n$-dimensional Y (we assume $m \geq n$) by maximal amount of indices will be $m-n$ dimensional tensor $Z$:
$$
X_{i_1 \ldots i_m} Y_{i_1 \ldots i_n} = Z_{i_{n+1} \ldots i_m}
$$
Namely, imagine, that the $x_t = (x_t)_{i_1 \ldots i_X}$ is an $d_x$-dimensional tensor. Then take $W_h$ to be  $d_{wh}$ dimensional tensor, hidden state $h$ to be $d_h = d_{wh} - d_{x}$ -dimensional, $U_h$ is $d_{uh}$ - dimensional, $W_y$ is $d_{wy}$ -dimensional. The output $y_t$ of the RNN will be $d_{wy} - d_{h}$ dimensional tensor.
Here note, that in order for the sum to be make sense in the first equation, dimensions need to match:
$$
d_{uh} - d_{h} = d_{h}
$$