So, why do we care about the general outer product $\mathbf u\otimes\mathbf v=\mathbf u\mathbf v^{\sf T}$?
I understand what it is, but why do we need it? It’s related to the cross product I understand, but why doesn’t that suffice?
The outer product $u v^T$ more often than not appears as an operator; for instance, as Botond pointed out, if $u$ is a unit vector, then $P_u = uu^T$ is the projector onto the $u$.
There are, however, many other uses. Some of which are:
We can express a matrix as the sum of outer products of its eigenvectors according to $$A = \lambda_1 u_1 v_1^T + \cdots \lambda_n u_n v_n^T,$$ where $\lambda_i$ is the $i$th eigenvalue and $u_i$ and $v_i$ are the $i$th right and left eigenvectors of $A$. This is sometimes called the spectral theorem.
This is idea from linear algebra that can be extended to general spectral theory where they are used to represent rank-1 operators. Rank one operators are used a lot in quantum mechanics: for instance, one central operator in Grover's algorithm, known as Grover's diffusion operator is given by $$ U_s = 2|s \rangle \langle s| - I,$$ where $|s \rangle \langle s|$ is the outer product of the vector $|s\rangle$ in bra-ket notation.
In numerical algebra, we can have matrices that are updated partially. For instance, given $A^{-1}$ and two vector $u$ and $v$, we may need to compute $$(A + uv^T)^{-1}.$$ This can be done with the Sherman-Morrison formula which gives us a numerically cheap way to compute the inverse of a matrix when it is perturbed by a rank-1 perturbation: $$ (A + uv^T)^{-1} = A^{-1} - \frac{A^{-1}u v^T A^{-1}}{1 + v^T A^{-1} u}.$$
A similar result holds for the determinant: $$ \det(A +u v^T) = (1 + v^T A^{-1} u)\det(A).$$ This is called the matrix determinant lemma.
A personal favorite of mine is to use of the duality between the inner product and the trace of the outer product: $$ v^Tu = \mathrm{trace}\{uv^T\}.$$
For instance, it allows for an easy (and fun) computation of the expected value of a $\chi^2$ random variable. If $X$ is $\chi^2$ distributed with $n$ degrees of freedom, then it can be seen as the sum of the squares of $n$ standard normal random variables: $$ X = E_1 ^2 + \cdots E_n^2.$$
Let $E$ be the vector whose $i$th component is $E_i$, then $E$ is an $n$ dimensional standard normal random vector, $$E \sim N(0,I).$$
Hence, $$ \mathbf E [X] = \mathbf E[E^TE] = \mathrm{trace}\{\mathbf E[ E E^T]\} = \mathrm{trace}\{I\} = n.$$
