This answer started off good but stumbled pretty badly near the end. Instead of talking about rotations about an axis, talk about rotations in a plane. In 2D, which we can identify with the complex plane, there is obviously one plane in which to rotate. Going up a dimension to 3D doesn't give us just one more plane but two more mutually orthogonal planes. The generators of the quaternion algebra, i.e. the $i$, $j$, $k$, directly correspond to the generators of rotations in each of these planes. This is most clearly captured by identifying the quaternion algebra with the even subalgebra of 3D geometric algebra in which case $i$, $j$, and $k$ are identified with the unit bivectors which directly and intuitively correspond to planes. Introducing a fourth spatial dimension confuses things immensely and doesn't mesh with your presentation. Moving to 3D immediately produces 3 square roots of -1. Your presentation makes it sound like adding a dimension adds a single "axis" to rotate about.

Let me start by saying that I came at this as someone who learnt complex numbers “properly”, but is entirely self-taught in quaternions out of sheer amateur interest. So I’m far from being an expert, and intended simply to share with the asker a mental model that worked for me.

With that proviso, I will dare to disagree with you. My axis of rotation and your plane of rotation would seem to be equivalent, given that the former is the normal of the latter.

If you think describing it in terms of bivectors is more intuitive, more power to you! Go ahead and add that as an answer. I never quite wrapped my head around them, myself.

Derek Elkins‭ And to a couple of your specific points that I didn’t have room for in the above…

“Moving to 3D immediately produces 3 square roots of -1.” Does it? I had thought it produced infinitely many roots lying on the circle $x=0,y^2+z^2=1$. And I worried that I’d glossed over that part too quickly by insisting on roots lying only on the orthogonal axes.

“Your presentation makes it sound like adding a dimension adds a single "axis" to rotate about.” Doesn’t it? A single orthogonal axis, yes… or so I thought!

Yes, I was implying orthogonal when I said it produces 3 as I said explicitly earlier in the comment. Even if you don't limit to orthogonal, infinitely many is still more than 1. For your latter question, in 4D with an orthogonal basis labelled $x$, $y$, $z$, $w$, what is the "axis" orthogonal to the $x$-$y$ plane? There are infinitely many lines orthogonal to that plane in 4D. You could say by "axis" you mean all the lines orthogonal to the plane of rotation, which would produce a 2-dimensional thing in 4D and different things in higher dimensions and nothing in 2D. This sort of works (again, in a way that can be naturally explained by geometric algebra) but is definitely not what people think of when they hear the term "axis". The benefit of the plane approach (as represented by a bivector) is that embedding it into higher dimensions changes nothing. The rotation in a plane represented by a bivector is exactly the same no matter the dimension of those vectors.

Even in 3D, the notion of axis is problematic. Physicists have a distinction between the notion of an "axial vector" and a "position vector". Reflecting a "position vector" in a plane orthogonal to it negates it, while reflecting an "axial vector" does nothing. When stated as a fact lacking context as I just have, this is not obvious. It is a bit more obvious when the "axial vector" is a derived result, e.g. a torque $\vec tau=\vec r \times \vec F$. If you recompute the torque after reflecting $\vec r$ and $\vec F$ in the plane orthogonal to $\vec tau$, i.e. the plane spanned by $\vec r$ and $\vec F$, then clearly $\vec r$ and $\vec F$ will be unchanged and thus the result will be unchanged. Again, everything is clear if we use rotation planes and bivectors instead of rotation axes. An "axial vector" is a dual bivector. The dualization process (and thus the cross product) is sensitive to the handedness of the space, but the plane of rotation, i.e. the bivector, is not.

Derek Elkins‭ Ah, I see my mistake now. I hadn’t even considered that the normal to a plane would no longer be unique in 4D. But just as there’s only one line orthogonal to a line in 2D, but infinitely many (a whole plane of them) in 3D, it makes sense that the same would happen to the normal(s) to a plane going from 3D to 4D. And my choice not to even try drawing the 4D quaternion space, but to instead stick to the 3D non-real sub-space, meant I never ran into the limitations of my assumption.

I think my answer still has value as a way to visualise, and maybe even to remember, the rules of quaternion multiplication. But now I’m worried I’ve missed some other error that makes it incorrect even for that purpose!

‘I was implying orthogonal when I said it produces 3….’ You’re still beyond me here. What three roots are they? I only count two new ones, plus the two that we had in 2D: $(0,\pm 1,0)$ and $(0,0,\pm 1)$.

‘…an "axial vector" and a "position vector".’ Ah, now we’re back on ground I know—or once knew, anyway. As far as I can tell, the distinction isn’t essential to understanding the nature of quaternion multiplication. But I could be wrong there too!

That part is talking about orthogonal planes. But where are you getting triplets in 2D? For the complex plane, we usually identify $(1,0)$ with 1 and $(0,1)$ with $i$. With complex multiplication, we then get $i^2 = -1$ and, modulo sign (the negation is not orthogonal), that's the only square root of $-1$. In 3D, there are infinitely many square roots of $-1$, one for each plane through the origin. However, these square roots of $-1$ can all be expressed as a linear combination of 3 orthogonal vectors (quaternions in this case) which themselves represent square roots of $-1$. These three orthogonal planes are each naturally represented by a bivector as is $i$ in the complex plane.

The digression on "position vectors" and "axial vectors" was just to illustrate why the notion of "axis" is problematic. The distinction doesn't matter for quaternions alone as vector quaternions are essentially all "axial vectors". (More clearly, vector quaternions correspond to bivectors.)

Derek Elkins‭, I'm not sure I get it either. If we have a rotation in a plane around the origin which, applied twice, takes us from $(1,0,0)$ to $(-1,0,0)$ then the plane must contain the origin, the start point, and the end point. By symmetry the square roots of $-1$ end up all being in the plane through the origin perpendicular to $(1,0,0) - (-1,0,0)$, so expressable as a linear combination of two vectors.

In quaternion algebra, $1$ is not identified with the vector $(1,0,0)$. If we used quaternion terminology for complex numbers, we'd say complex numbers have a scalar part (the real part) and a 1D vector part (the imaginary part). These spaces are (double covers of) spaces of rotation operators that act upon a vector space (2D for complex, 3D for quaternions). Complex numbers let us conflate these because a vector can be conflated with the complex number that rotates $(1,0)$ to that vector. Solving $q^2=-1$ in quaternions means find the operator represented by $q$ that when composed with itself gives the operator represented by $-1$. Quaternions representing rotations in the y-z plane solve this equation even though they won't send $(1,0,0)$ to $(-1,0,0)$.

Your constraints amount to asking for the quaternions whose plane of rotation is represented by a bivector with $0$ $y$-$z$ component. This will restrict to a 2D subspace of bivectors, but there's nothing special about this. The space of rotation operators with $0$ $y$-$z$ component, however, is not closed under composition (quaternion multiplication) so this doesn't identify a subalgebra.

Derek Elkins‭ “But where are you getting triplets in 2D?” From the trivial $(x,y)\mapsto (x,y,0)$.

“…modulo sign (the negation is not orthogonal)…” Ah, see, I was counting the negations, so your three I would’ve called six. But I still don’t see how we get that many! In two dimensions, there’s one (pair of) root(s), called $\pm i$ (as either a complex number or a quaternion). In four dimensions, there are three (pairs of) roots, called $\pm i$, $\pm j$, and $\pm k$. But surely in three dimensions, there are only two (pairs), because the $k$ dimension is not in consideration? Where am I going wrong?

“In quaternion algebra, $1$ is not identified with the vector $(1,0,0)$.” I think I’m starting to get it. Or rather, I think I understand it, but I still don’t get it. My lack of formal quaternion education is probably betraying me here. Why doesn’t the real unit 1 correspond to the 4D vector $(1,0,0,0)$? And $i=(0,1,0,0)$, $j=(0,0,1,0)$ and $k=(0,0,0,1)$?

Because that’s what I was assuming, in my answer and in my comments claiming to only find two (pairs of) roots in three dimensions. Now, if your three dimensions are $i$, $j$ and $k$, and it’s the real dimension that’s not in consideration, then yes, in that 3D space, all three (pairs of) roots exist. But how is that helpful, when $-1$ doesn’t exist?

If you want to talk about square roots of -1, you need to define what "-1" means and what multiplication means so that you can solve the equation $x \cdot x=-1$. If you want to use the embedding $(x,y)\mapsto(x,y,0)$, that's fine, but you haven't defined a multiplication on triplets nor what -1 means in this case, so it's not meaningful to talk about "square roots of -1" in this context. The algebra of quaternions is built on a 4-dimensional vector space and we certainly can identify the generators of the quaternion algebra with the 4D vectors you mention.

To be clear, when I say "in 3D", I mean in the (even subalgebra of the) geometric algebra on a 3-dimensional space. The geometric algebra on an $n$-dimensional space is itself a $2^n$ dimensional vector space. For $n=2$, we have a 4-dimensional space and its even subalgebra is 2D and corresponds to the complex numbers. For $n=3$, we have an 8-dimensional space and its even subalgebra is 4D and corresponds to the quaternions.