'Higher order' complex numbers I recently learned that complex numbers can also be represented as matrices of the form: $$\begin{pmatrix} x & -y \\ y & x\end{pmatrix}$$ where complex multiplication corresponds to matrix multiplication and conjugation corresponds to transposition. I thought it was interesting how this representation corresponds to the $2 \times 2$ unit rotation matrix:
$$\begin{pmatrix} \cos \theta & -\sin\theta \\ \sin \theta & \cos \theta \end{pmatrix}$$
I was wondering if, then, there is a higher order version of the complex numbers? e.g.  numbers that correspond to the $3 \times 3$ unit rotation matrix:
$$\begin{pmatrix} \sin \theta \cos \phi & \sin \theta \sin \phi & \cos \theta \\ \cos \theta \cos \phi & \cos \theta \sin \phi & - \sin \theta \\ -\sin \phi & \cos \phi & 0\end{pmatrix}$$ etc...
 A: It's a good question, albeit one already long since pondered.
In short: there isn't, at least not one that maintains a lot of the properties we would likely want from such a system.

To clarify, let's first note that the next step up is quaternions, four-dimensional instead of three- or two-dimensional. They have the natural representation as numbers of the form
$$a + b i + cj + dk$$
where
$$i^2 = j^2 = k^2 = ijk = -1$$
We can represent these as matrices (either in $M_{2\times 2}(\mathbb{C})$ or $M_{4 \times 4}(\mathbb{R})$) (Wikipedia reading):
$$\begin{bmatrix}a+bi & c+di \\ -c + d i & a - b i \end{bmatrix}$$
for complex matrices, and one possible (non-unique, among $48$ total) representation over $\mathbb{R}$ by
$$\begin{align*}
&\begin{bmatrix}
 a & -b & -c & -d \\ 
 b & a & -d & c \\
 c & d & a & -b \\
 d & -c & b & a 
\end{bmatrix}\\
&= a
\begin{bmatrix}
 1 & 0 & 0 & 0 \\ 
 0 & 1 & 0 & 0 \\
 0 & 0 & 1 & 0 \\
 0 & 0 & 0 & 1 
\end{bmatrix}
+ b
\begin{bmatrix}
 0 & -1 & 0 & 0 \\ 
 1 & 0 & 0 & 0 \\
 0 & 0 & 0 & -1 \\
 0 & 0 & 1 & 0 
\end{bmatrix}
+ c
\begin{bmatrix}
 0 & 0 & -1 & 0 \\ 
 0 & 0 & 0 & 1 \\
 1 & 0 & 0 & 0 \\
 0 & -1 & 0 & 0 
\end{bmatrix}
+ d
\begin{bmatrix}
 0 & 0 & 0 & -1 \\ 
 0 & 0 & -1 & 0 \\
 0 & 1 & 0 & 0 \\
 1 & 0 & 0 & 0 
\end{bmatrix}
\end{align*}$$

Loosely, Hamilton (the discoverer of quaternions) sought a way to multiply and divide three-dimensional numbers, in the same sense we can do this for two-dimensional numbers.
More specifically: in $\mathbb{C}$, we can naturally add and multiply numbers together in the usual way, and even subtract and divide them. The complex numbers form what we call a "field" owing to its properties; another, and more relevant, term is that of "associative division algebra".
Hamilton wondered if we could extend these notions to a three-dimensional version.
In his attempts, issues arose however in trying to make division (or multiplying by a multiplicative inverse) work out for all elements of the set. As it turns out, this is provably impossible, through the fact that there is no three-dimensional associative division algebra over $\mathbb{R}$ (Frobenius' theorem & proof).

Granted, your question seems a little more grounded in the language of matrices. Let's notice something: if we let $\sin \theta = x, \cos \theta = y$, then your matrix representations in the two-dimensional case notably correspond identically (in a loose sense).
Performing this for the three dimensional matrix, you should notice we already need four variables, not three, in a similar sense.
And, of course, if we wish to keep this grounded towards matrices and not necessarily in the language of division algebras... Well, I don't have a conclusive answer to offer you there.
What I can tell you is that the system of numbers such matrices would naturally correspond to would be, in some way or another ... well, jank. They would be missing properties we typically like to have.
Whether you like to have those properties is, of course, entirely up to you. There's nothing stopping such a thing from being useful or interesting after all, in the right context.
A: Usually with even number of dimensions you have better results than with odd. Thus, you do not have anything useful or interesting in 3 dimensions. But in 4 dimensions you have multiple interesting algebras, some of which are more or less numbers-like.
Eevee Trainer has already mentioned quaternions, which are a division algebra but not commutative, but there are also other approaches which give algebras with different properties in 4D.

*

*First of all, possibly the simplest approach is to consider the $2\times2$ real matrices themselves, which is a 4D space. Matrix operations on $2\times2$ matrices are isomorphic to so-called split-quaternions:

$\begin{align}
\boldsymbol{1} =\begin{pmatrix}1&0\\0&1\end{pmatrix},\qquad&\boldsymbol{i} =\begin{pmatrix}0&1\\-1&0\end{pmatrix},\\
\boldsymbol{j} =\begin{pmatrix}0&1\\1&0\end{pmatrix},\qquad&\boldsymbol{k} =\begin{pmatrix}1&0\\0&-1\end{pmatrix}.\end{align}$
They are neither commutative, nor a division-algebra, but allow to apply any functions that can be generalized to square matrices. They have zero divisors, nilpotents and idempotents.


*Dual complex numbers. Formed by adding the dual unity $\varepsilon$ such that $\varepsilon^2=0$ to complex numbers. They can be represented by matrices of the form
$\left(
\begin{array}{cc}
 u & w \\
 0 & v \\
\end{array}
\right)$, where $u,v,w$ are complex numbers or matrices of the form $\left(
\begin{array}{cccc}
 a & b & c & d \\
 -b & a & -d & c \\
 0 & 0 & a & b \\
 0 & 0 & -b & a \\
\end{array}
\right)$ or $\left(
\begin{array}{cccc}
 a & c & b & d \\
 0 & a & 0 & b \\
 -b & -d & a & c \\
 0 & -b & 0 & a \\
\end{array}
\right)$ where $a,b,c,d$ are real numbers ($a$ is real part, $b$ is imaginary part, $c$ is dual part, $d$ is imaginary dual part). This system is commutative but has zero divisors and nilpotents.


*Tessarines. Formed by adding split-complex unity $j$ such as $j^2=1$ to complex numbers. Can be represented as matrices of the form $\left(
\begin{array}{cc}
 u & w \\
 w & u \\
\end{array}
\right)$, where $u,w$ are complex numbers or matrices of the form $\left(
\begin{array}{cccc}
 a & b & c & d \\
 -b & a & -d & c \\
 c & d & a & b \\
 -d & c & -b & a \\
\end{array}
\right)$ or $\left(
\begin{array}{cccc}
 a & c & b & d \\
 c & a & d & b \\
 -b & -d & a & c \\
 -d & -b & c & a \\
\end{array}
\right)$ where $a,b,c,d$ are real numbers ($a$ is real part, $b$ is imaginary part, $c$ is split-imaginary part, $d$ is the coefficient of $ij$ or $-ij=k$). This system is commutative algebra, but it has zero divisors and idempotents.
If one adds another complex unity $i_1$ such that $i_1^2=-1$ instead of split-complex unity $j$, one will get bicomplex numbers, which are isomorphic to this system (the split-complex unity will arise automatically as $j=-i i_1$).


*Dual split-complex numbers. They are formed by adding split-complex unity $j$ and dual unity $\varepsilon$ to real numbers. They can be represented as real matrices of the form $\left(
\begin{array}{cccc}
 a & b & c & d \\
 b & a & d & c \\
 0 & 0 & a & b \\
 0 & 0 & b & a \\
\end{array}
\right)$ or $\left(
\begin{array}{cccc}
 a & c & b & d \\
 0 & a & 0 & b \\
 b & d & a & c \\
 0 & b & 0 & a \\
\end{array}
\right)$ where $a,b,c,d$ are real numbers ($a$ is real part, $b$ is split-imaginary part, $c$ is dual part, $d$ is split-imaginary dual part).
This is a commutative algebra but with zero divisors, nilpotents and idempotents.


*4D split-complex numbers. These are formed by adding two split-complex unities $j$ and $j_1$ such that $j^2=j_1^2=1$. They are isomorphic to the matrices of the form $\left(
\begin{array}{cccc}
 a & b & c & d \\
 b & a & d & c \\
 c & d & a & b \\
 d & c & b & a \\
\end{array}
\right)$ or $\left(
\begin{array}{cccc}
 a & c & b & d \\
 c & a & d & b \\
 b & d & a & c \\
 d & b & c & a \\
\end{array}
\right).$ This is a commutative system, with zero divisors and idempotents.


*4D dual numbers. These are formed by adding two dual unities to real numbers: $\varepsilon$ and $\varepsilon_1$, such that $\varepsilon^2=\varepsilon_1^2=\varepsilon\varepsilon_1=0$. They are isomorphic to the real matrices of the following form:
$\left(
\begin{array}{cccc}
 a & b & c & d \\
 0 & a & 0 & c \\
 0 & 0 & a & b \\
 0 & 0 & 0 & a \\
\end{array}
\right)$. This system is commutative but has zero divisors and nilpotents.
A: This is discussed a lot in physics. You're representing your complex numbers in the form $m\cdot e^{i\cdot\theta}$ here, where the direction-part, the unit-modulus $e^{i\cdot\theta}$, can be seen as an element of the special orthogonal group $\mathrm{SO}(2)$, aka the circle group, which is what you can represent with those $2\times2$ rotation matrices. Then, $\mathbb{C}$ as a whole is the group product of the (multiplicative) group $\mathbb{R}^+$ with $\mathrm{SO}(2)$.
...Well, almost. $\mathbb{R}^+$ is only a group if you exclude zero, i.e. $\mathbb{R}_{>0}$. As a consequence you also only get $\mathbb{C}\setminus\{0\}$ from this construction.
By considering $3\times3$ rotation matrices instead, you're swapping out the group for $\mathrm{SO}(3)$. That is, well, a group just as good as $\mathrm{SO}(2)$. In that sense, all works much the same – we get the product group $\mathbb{R}_{>0}\oplus\mathrm{SO}(3)$.
The complex numbers are more than just a group under multiplication though. They are also a group under addition. That might seem to be the simpler of the two operations, but actually it doesn't in any way fit in the abstract product-of-groups formulation. The way to get an addition on our $\mathbb{C}$ construction is to embed the $\mathrm{SO}(2)$ mappings into the space of arbitrary linear endomorphisms on $\mathbb{R}^2$, aka the set of arbitrary $2\times2$ matrices. Specifically, we have
$$\begin{align}
  E_n:& (\mathbb{R}_{>0}\oplus\mathrm{SO}(n)) \to \mathrm{GL}(n)
\\ & E_n(m, r) = mr
\end{align}$$
where $mr$ just denotes scaling of the linear operator given by $r\in \mathrm{SO}(n)\subset \mathrm{GL}(n)$.
On matrices, addition is simply component addition. However, a priori there is no reason why the sum of two matrices should again correspond to an element of our group product. This happens to work out in the 2D case, but not in 3D. For example,
$$
  \begin{pmatrix}0 & 1 & 0 \\ 0 & 0 & 1 \\ 1 & 0 & 0 \end{pmatrix}
  + \begin{pmatrix}0 & 1 & 0 \\ 1 & 0 & 0 \\ 0 & 0 & -1 \end{pmatrix}
 = \begin{pmatrix}0 & 2 & 0 \\ 1 & 0 & 1 \\ 1 & 0 & -1 \end{pmatrix}
$$
which maps $\mathbf{e}_\mathrm{x}$ to a vector of length $\sqrt2$, but $\mathbf{e}_\mathrm{y}$ to a vector of length $2$, something that is impossible for $\mathbb{R}_{>0}\oplus\mathrm{SO}(3)$.
