Snell's law in vector form
Let a unit vector $\;\mathbf{n}=(\rm n_1,n_2,n_3)\,, \Vert\mathbf{n}\Vert=1$. Any vector $\;\mathbf{r}\;$ could be decomposed in two components with respect to $\;\mathbf{n}\;$, see Figure-01 in the bottom \begin{equation} \mathbf{r}=\mathbf{r}_\|+\mathbf{r}_\bot \tag{01}\label{01} \end{equation} one parallel and the other normal to axis $\mathbf{n}$ respectively \begin{align} \mathbf{r}_\| &=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{r}\right)\mathbf{n} \tag{02.1}\label{02.1}\\ \mathbf{r}_\bot &=\left(\mathbf{n}\boldsymbol{\times}\mathbf{r}\right)\boldsymbol{\times}\mathbf{n}= \mathbf{r}-(\mathbf{n}\boldsymbol{\cdot}\mathbf{r})\mathbf{n} \tag{02.2}\label{02.2} \end{align} that is \begin{equation} \mathbf{r}=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{r}\right)\mathbf{n}+\left(\mathbf{n}\boldsymbol{\times}\mathbf{r}\right)\boldsymbol{\times} \mathbf{n} \tag{03}\label{03} \end{equation} The vectors $\;\mathbf{t},\mathbf{i}\;$ are decomposed as follows \begin{equation} \mathbf{t}=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)\mathbf{n}+\left(\mathbf{n}\boldsymbol{\times}\mathbf{t}\right)\boldsymbol{\times} \mathbf{n} \tag{04}\label{04} \end{equation} \begin{equation} \mathbf{i}=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}+\left(\mathbf{n}\boldsymbol{\times}\mathbf{i}\right)\boldsymbol{\times}\mathbf{n} \tag{05}\label{05} \end{equation} Now, Snell's Law is expressed as \begin{equation} \left(\mathbf{n}\boldsymbol{\times}\mathbf{t}\right)=\mu\left(\mathbf{n}\boldsymbol{\times}\mathbf{i}\right) \tag{06}\label{06} \end{equation} see Figure-02 in the bottom.
Equation \eqref{04} combined with \eqref{05} and \eqref{06} yields \begin{equation} \mathbf{t}=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)\mathbf{n}+\mu\left[\mathbf{i}-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}\right] \tag{07}\label{07} \end{equation} Taking norms in \eqref{07} and since the vector $\;\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)\mathbf{n}\;$ is normal to the vector $\;\left[\mathbf{i}-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}\right]\;$ \begin{equation} \Vert\mathbf{t}\Vert^2=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)^2+\mu^2\Vert\mathbf{i}-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}\Vert^2 \tag{08}\label{08} \end{equation} or \begin{equation} 1=\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)^2+\mu^2\left[1-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)^2\right] \tag{09}\label{09} \end{equation} so \begin{equation} \left(\mathbf{n}\boldsymbol{\cdot}\mathbf{t}\right)=\pm\sqrt{1-\mu^2\left[1-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)^2\right]} \tag{10}\label{10} \end{equation} Since the angle between $\;\mathbf{n},\mathbf{t}\;$ is less than $\;\pi/2\;$ we keep the plus sign in \eqref{10} and \eqref{07} yields finally \begin{equation} \mathbf{t}=\sqrt{1-\mu^2\left[1-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)^2\right]}\mathbf{n}+\mu\left[\mathbf{i}-\left(\mathbf{n}\boldsymbol{\cdot}\mathbf{i}\right)\mathbf{n}\right] \tag{11}\label{11} \end{equation}
Two parts of this answer: (1) a quick sketch of how one gets the result cited and (2) why I believe one can't get a less awkward, more succinct result.
Proof Sketch
The one thing that the numerical version of Snell's law does not give us, but which we always unconsciously use, is the equally important fact that the incident $\vec{i}$ and refracted rays $\vec{r}$ and the interface normal $\vec{n}$ all lie in the same plane. In other words, the refracted ray lies in the span $\left<\left\{\vec{i},\,\vec{n}\right\}\right>$ of $\vec{i}$ and $\vec{n}$ and so there are constants $\alpha$ and $\beta$ such that:
$$\vec{r} = \alpha\,\vec{i}+\beta\,\vec{n}\tag{1}$$
Another implication of this assertion is as follows: the vector $\vec{n}\times\vec{i}$ is the normal to this span plane. So is $\vec{n}\times\vec{r}$. So therefore, $\vec{n}\times\vec{r}\propto \vec{n}\times\vec{i}$. Now use Snell's numerical law to find the proportionality constant and you get the first equation you cite.
Take the cross product of both sides of (1) with $\vec{n}$:
$$\vec{n}\times\vec{r}=\alpha\,\vec{n}\times\vec{i}\tag{2}$$
and apply the form of Snell's law you have cited to find:
$$\alpha=\mu\tag{3}$$
Let's agree to work with normalized vectors, in which case taking the inner product of (1) with itself yields:
$$0 = 1-\mu^2 + \beta^2 + 2\,\mu\,\beta\, \langle \vec{i},\,\vec{n}\rangle\tag{4}$$
whence, on finding $\beta$ from (3), the result follows with a bit of messing about. You also need to use the fact that the components of $\vec{r}$ and $\vec{i}$ normal to the interface are in the same, not the opposite, direction to decide the sign of the root of the equation for $\beta$.
A Simpler Version?
This is a complicated and awkward expression, hard to work with, so one wonders whether there might be an eleganter expression for something so trivial as Snell's law. I don't believe so, and this is because:
(1) Snells law results from a boundary condition involving medium optical densities, and is thus not simply a geometric result. If it were, simple vector or tensor expressions would result;
(2) As a result of (1), Snell's law is a nonlinear function of $\vec{i}$ and the ratio $\mu$. It is this nonlinearity that begets optical wavefront aberration in lens systems. You can write the equation derived above as:
$$\vec{r} = \vec{n} + \mu\,(\mathrm{id}-\vec{n}\otimes\vec{n})\,\vec{i} +\left(-\frac{\mu^2}{2}\left(1-\langle\vec{n},\,\vec{i}\rangle^2\right)+\frac{\mu^4}{8}\left(1-\langle\vec{n},\,\vec{i}\rangle^2\right)^2 + \cdots\right)\vec{n}$$
The first two terms above are the normal component ($\vec{n}$) and $\mu$ times the projection $(\mathrm{id}-\vec{n}\otimes\vec{n})\,\vec{i}$ projection of $\vec{i}$ onto the interface. This second term expresses the conservation of the transverse optical momentum that upholds the Hamiltonian description of ray optics and conservation of Étendue (the optical Liouville Theorem or second law of thermodynamics applied to light) across discontinuous interfaces. See the latter part of my answer here for a discussion of this point. The other terms are those responsible for optical aberrations; consideration of the linear terms together with the two written explicitly out above alone gives the Seidel theory of aberration. If Snell's law were simply expressible in vector notation, lens design would be trivial.
(3) If we wanted a purely geometric theory, one could perhaps begin with defining the metric tensor in the medium as $n^2(x,\,y,\,z)$ times the identity matrix; the distance function is then the optical path difference, and, as discussed in my answer, ray optics becomes the purely geometric theory of finding geodesics in this conformally flat Riemannian geometry, with the geodesics found from the Euler-Lagrange equation $\ddot{X}^k + \Gamma^k_{ij}\,\dot{X}^i\,\dot{X}^j = 0$ parameterized by the affine parameter equal to the optical pathlength. So we might hope for simple vector expressions expressing purely geometrical relationships in this case. Unfortunately, as discussed in the answer, the six dimensional Hamiltonian / Riemann-geometric approach breaks down at discontinuous interfaces and we are forced back to the four dimensional Hamiltonian approach involving only transverse components, because Snell's law only conserves the transverse optical momentums but not the normal and so we are again forced to deal with awkward boundary conditions. So again I think we are forced to the same conclusion that no eleganter geometric expression exists for Snell's law.