Suppose we have a point $x_0\in{\bf H}$ (where by $\bf H$ I denote the ring of quaternions).
What I'm curious about is what can the set of solutions of $x^2=x_0$ look like?
From what I've checked, for $x_0$ a positive real there are only the two real solutions, as we can write $(a+bi+cj+dij)^2=(a^2-b^2-c^2-d^2)+2abi+2acj+2adij$, which implies $b=c=d=0$ if the result is to be a positive real. On the other hand, for real $x_0<0$, all solutions must have zero real part and the solutions form a $3$-sphere, while for $x_0=r\cdot i$ ($r$ real) there are again only two solutions. Is that right, or did I make some stupid mistake?
The general question I'm interested in is, what do solutions of $x^n=x_0$ look like (geometrically) for arbitrary $n\in{\bf N},x_0\in{\bf H}$?
(I'm tagging this with algebraic geometry because I'm actually asking about some specific varieties, though I'm aware that the context is rather unusual; feel free to remove the tag if you find it inappropriate.)
Edit: Using the argument in Lutzl's answer, we can see that there are exactly $n$ solutions whenever $x_0$ is nonreal. Furthermore, it's rather obvious we can assume $\lvert x_0\rvert=1$ (as scaling $x_0$ by a positive real factor will only scale the set of solutions by an appropriate positive real root, and for $x_0=0$ the only solution is zero, since quaternions are a domain), so that leaves us with the case $x_0=1$ and $x_0=-1$, and I'm interested in what do solutions look like exactly, like in case of $n=2$ we know that they are two distinct points or a single $3$-sphere, respectively.
Assuming that $x_0$ is not purely real, it defines a unique embedding of the complex plane into the quaternions and any solution of a polynomial equation with real coefficients and $x_0$ as constant term will only have solutions in that embedded complex plane.
Split $x_0=r_0+s_0e_0$ into its real and imaginary part where $e_0=Im(x_0)/|Im(x_0)|$ is a complex unit with $e_0^2=-1$. Consider any polynomial $p\in \mathbb R[x]$. Then this polynomial evaluated at some quaternion $x=r+se$ with imaginary unit $e^2=-1$ will result in a value $p(x)=p(r+se)\in \mathbb R[e]$. So that $p(x)=x_0$ is only possible if $e$ and $e_0$ are essentially the same imaginary quaternion unit, up to a sign. That is, any solution of $p(x)=x_0$ is contained in the embedding $\mathbb R[e_0]$ of the complex plane, i.e., $x=r+se_0$.
This is true as long as no other non-commuting complex units enter the picture.
Update: The $n$-th degree roots: if $x_0=\cos(ϕ)+\sin(ϕ)e$ is a unit quaternion with imaginary unit $e$, then the solution set of $x^n=x_0$ is as in the complex plane $x_m=\cos(\tfrac{ϕ+2πm}{n})+\sin(\tfrac{ϕ+2πm}{n})e$, $k=1,...,n$.
However, in the case $\sin(ϕ)=0$, i.e., $ϕ=ℓπ$, there is no imaginary unit determined by $x_0$ and so any imaginary unit will produce a solution. Instead of single points the solution set thus consists of 2-spheres $$\left\{\cos(\tfrac{ℓπ+2πk}{n})+\sin(\tfrac{ℓπ+2πk}{n})e\middle|Re(e)=0,\ |e|=1\right\}.$$
In the case $n=2$ and $x_0=1$, $ℓ=0$, the radii of the spheres are $\sin(kπ)=0$, thus the solution set is the two single points $x=\pm 1$. In the case $x_0=-1$, $ℓ=1$, the radius of the single sphere is $\sin(\tfrac{π}{2}+πk)=(-1)^k$. All conjugate pairs $\pm e$ are on the same sphere of purely imaginary unit quaternions, so here the solution set has only one component.
This remains true in general, the solution sets for $x^n=\pm1$ are 2-spheres inside the unit 3-sphere perpendicular to the real line and with real parts corresponding to the real parts of the usual complex $n$-th unit roots. Since complex conjugate solution pairs have the same real part, there are about $n/2$ component spheres in the solution set.
Loosely related comment: In general, it is a rather convoluted task to say what a polynomial over the quaternions with quaternion coefficients should be, a quadratic polynomial could have the general form
$$f(x)=a_0xa_1xa_2+b_0xb_1+c_0.$$
In the context of Julia sets and Mandelbrot fractals, one explores the properties of the quadratic iterations $x_+=f(x)$. Since the quality of the iteration remains unchanged under linear transformations $w=u_0xu_1+v_0$ one can try to reduce the number of free real parameters in those iterations.
In the classical complex case, the same reduction leads to the normal form of the complex one-parametric (real two-parametric) Mandelbrot iteration $z_+=z^2+c$. In the quaternion case there remain more parameters, I'd say (10+7+4)-(7+4)=10 real parameters for the quaternionic Mandelbrot fractal.