If *y* is a quaternion, how many solutions does *x*² = *y* have? That is, does every quaternion have a square root, and if so, how many square roots does it have?

A quaternion can have more than two roots. There is a example right in the definition of quaternions:

*i*² = *j*² = *k*² = -1

and so -1 has at least three square roots. It follows that every negative real number has at least three square roots.

## Negative reals

In fact, every square root of -1 corresponds to a square root of any other negative number. To see this, let *p* be a positive real number and let *r* = √*p* > 0. If *x* is a square root of -1, then *rx* is a square root of –*p*. Conversely, if *x* is a square root of –*p*, then *x*/*r* is a square root of -1.

So all negative real numbers have the same number of square roots in the quaternions, and that number is at least 3. In fact, it turns out that all negative real numbers have infinitely many quaternion roots. A calculation shows that

(*a* + *bi* + *cj* + *dk*)² = -1

if and only if *b*²* + c*²* + d*² = 1. So there are as many square roots of -1 as there are points on a sphere.

## Not negative reals

If a quaternion is not a negative real number, then the story is much more familiar: the equation *x*² = *y* has one solution if *y* = 0 and two solutions otherwise.

To summarize, *x*² = *y* has

- one solution if
*y*= 0, - infinitely many solutions if
*y*is negative and real, and - two solutions otherwise.

It’s easy to show that 0 has only one root: |*x*²| = |*x*|², and so if *x**² *= 0, then |*x*| = 0 and so *x* = 0. Here the norm of a quaternion is defined by

|*a* + *bi* + *cj* + *dk| = *√(*a²* + *b²* + *c²* + *d²*).

If *p* > 0, the solutions to *x*² = –*p* are *bi* + *cj* + *dk* where *b*²* + c*²* + d*² = *p*.

What’s left is to show that if *y* is not zero and not a negative real, then *x*² = *y* has two solutions.

## Polar form and calculating roots

Quaternions have a polar form analogous to the polar form for complex numbers. That is, we can write a quaternion *y* in the form

*y* = *r*(cos θ + *u* sin θ)

where *r* is real and positive, θ is real, and *u* is a unit quaternion with zero real part. We can find *r*, θ, and *u* as follows. Write *y* as

*a* + *bi* + *cj* + *dk* = a + *q*

so that *a* is the real part and *q* is the pure quaternion part. Then we can find *r*, θ, and *u* by

*r* = |*y*|

cos θ = a/|*y*|

*u* = *q* / |*q*|

Then the two roots are

*x* = ±√*r *(cos θ/2 + *u* sin θ/2).

## Python code

First we’ll implement the algorithm above, then show that it produces the same results as the Python package `quaternionic`

.

First, our code to compute *x* satisfying *x*² = *y*.

import numpy as np np.random.seed(20210106) def norm(q): return sum(t*t for t in q)**0.5 y = np.random.random(size=4) r = norm(y) theta = np.arccos(y[0]/r) u = y + 0 # make a copy u[0] = 0 u /= norm(u) x = np.sin(theta/2)*u x[0] = np.cos(theta/2) x *= r**0.5

And now we check our results using `quaternionic`

.

import quaternionic qx = quaternionic.array(x) qy = quaternionic.array(y) print(f"y = {y}") print(f"x*x = {qx*qx}") # should equal y print(f"x = {x}") # our root print(f"x = {np.sqrt(qy)}") # root via package

This produces

y = [0.61615367 0.07612092 0.09606777 0.11150865] x*x = [0.61615367 0.07612092 0.09606777 0.11150865] x = [0.79189641 0.04806243 0.06065678 0.07040609] x = [0.79189641 0.04806243 0.06065678 0.07040609]

## More quaternion posts

- Finnegans Wake and quaternions
- How far is
*xy*from*yx*on average for quaternions? - Probability of commuting

For an exploration of the roots of general polynomial equations over polynomials, see Equations in Quaternions by Ivan Niven. The American Mathematical Monthly, Dec., 1941, Vol. 48, pp. 654-661.

Seems like there is a geometric interpretation: Quaternions are 3D rotations. So each square root is half of a “flip movement”, and there are sphere-surface-many points defining the flip plane.

An interesting thing about the polar form is that it tells you how to extend analytic functions of a single variable to quaternions. If a is the real part and q is the purely quaternionic part (I’ll assume non-zero here) and r = |q|, u = q/r then q = ru, q^2 = -r^2, q^3 = -r^3 u, etc. In algebraically, u always behaves like the complex number “i”.

So if z = x + i y and f(z) = U(x,y) + i V(x,y) an analytic function, a quaternionic extension for f is f(a+q) = U(a, r) + V(a,r) u.

If you can compute a function as a function of a complex variable, you can modify the method to compute quaternionic results for you, as long as it’s a single variable.

The thing I like about this is that it gives you the point of leverage to characterize all the field automorphisms of H.

You’ve probably seen the proof that id is the only field automorphism of R, where you start by showing that only positive numbers have square roots? And of course “# of square roots” is preserved by an automorphism, so positives are mapped to positives –> order preserving + fixed on Q –> must be the identity.

For H, “having more than 2 square roots” is automorphism preserved, so negative reals are mapped to negative reals –> R is mapped to R –> must be identity on R –> R-linear. So your range of possibilities is reduced to 4×4 real matrices, and a bit of computation shows that conjugations are the only field automorphisms. It’s one of my favorite proofs.