Why is there this relationship between quaternions and Pauli matrices?

QMechanics has given you the straight answer in terms of group isomosphisms. So please go with that, but in case you go further into mathematics of quarterions and their applications in physics you will find many twists and turns on this subject.

I personally found the book On Quaternions and Octonions : Conway, Smith (2003) finally gave me some clarity on this whole subject. I will summarise a few key points. Apologies if this goes a bit further than your original question.

Qauternions are part of a series of division algebras used by mathematicians. They only occur in dimensions which are powers of 2 but only up to 8, namely:
1. Real numbers
2. Complex numbers
4. Quaterions
8. Octonians

As you should know, unit complex numbers are related to rotations in 2 dimensions - you might expect this to be part of a pattern and indeed it is roughly stated (from page 89):

• multiplication by unit complex numbers generate rotations in 2 dimensions
• multiplication by unit quaternions generate rotations in 4 dimensions (not 3 dimensions - see below!)
• multiplication by unit octonions generate rotations in 8 dimensions

The subtlety (which links back to your question) is that in 4 dimensions there are 6 rotations (you will know this if you have studied special relativity) so you actually need 2 copies of the quaternions. If you restrict to just one copy you get back to 3 dimension rotations.

In summary:
3d rotations: one copy of unit quaternions relate to pauli matrices
4d rotations: two copy of unit quaternions relate to 2 copies of pauli matrices

In group language:
Spin(3) = SU(2) (3 dimensions)
Spin(4) = SU(2) x SU(2) (6 dimensions)

As 3 and 4 dimensions are the two most important to physics this comes up in many guises in Quantum Physics.

1. Pauli matrices-Rotations-Special Unitary matrices $:mathrm{SU}(2):$

Any vector in $mathbb{R}^3$ can be represented by a $2times2$ hermitian traceless matrix and vice versa. So, there exists a bijection (one-to-one and onto correspondence) between $mathbb{R}^3$ and the space of $2times2$ hermitian traceless matrices, let it be $mathbb{H}$ : begin{equation} mathbf{x}=(x_1,x_2,x_3)in mathbb{R}^3;longleftrightarrow ; X= begin{bmatrix} & x_3 & x_1-ix_2 & x_1+ix_2 & -x_3 end{bmatrix} in mathbb{H} tag{001} end{equation} From the usual basis of $mathbb{R}^3$ begin{equation} mathbf{e}_{1}=left(1,0,0right),quad mathbf{e}_{2}=left(0,1,0right),quad mathbf{e}_{3}=left(0,0,1right) tag{002} end{equation} we construct a basis for $mathbb{H}$ begin{eqnarray} mathbf{e}_1 &=&(1,0,0)qquad longleftrightarrow qquad sigma_1= begin{bmatrix} &0&1& &1&0& end{bmatrix} tag{003a} mathbf{e}_2 &=&(0,1,0)qquad longleftrightarrow qquad sigma_2= begin{bmatrix} &0&-i &i&0 end{bmatrix} tag{003b} mathbf{e}_3 &=&(0,0,1)qquad longleftrightarrow qquad sigma_3= begin{bmatrix} &1&0 &0&-1 end{bmatrix} tag{003c} end{eqnarray} where $:boldsymbol{sigma}equiv(sigma_{1},sigma_{2},sigma_{3}):$ the Pauli matrices⁽¹⁾, essentially the components of the spin $:s=1/2:$ angular momentum by a factor $:1/2:$ begin{equation} S_1=dfrac{1}{2}sigma_{1};, quad S_2=dfrac{1}{2}sigma_{2};, quad S_3=dfrac{1}{2}sigma_{3}, quad text{or} quad mathbf{S}=dfrac{1}{2}boldsymbol{sigma} tag{004} end{equation} Suppose now that the vector $:mathbf{x}=(x_1,x_2,x_3):$ is rotated around an axis with unit vector $:mathbf{n}=(n_1,n_2,n_3)$ through an angle $theta$⁽²⁾ begin{equation} mathbf{x}^{prime}= costheta ;mathbf{x}+(1-costheta);(mathbf{n}boldsymbol{cdot}mathbf{x});mathbf{n}+sintheta;(mathbf{n}boldsymbol{times}mathbf{x}) tag{005} end{equation} and let to the vectors $:mathbf{x},mathbf{x}^{prime}:$ correspond the matrices begin{eqnarray} X & equiv & mathbf{x}boldsymbol{cdot} boldsymbol{sigma} = x_1sigma_1+x_2sigma_2+x_3sigma_3= begin{bmatrix} x_3&x_1-ix_2 x_1+ix_2&-x_3 end{bmatrix} tag{006a} X{'} & equiv & mathbf{x}{'}boldsymbol{cdot} boldsymbol{sigma} = x_1^{'}sigma_1+x_2^{'}sigma_2+x_3^{'}sigma_3= begin{bmatrix} x^{'}_3&x^{'}_1-ix^{'}_2 x^{'}_1+ix^{'}_2&-x^{'}_3 end{bmatrix} tag{006b} end{eqnarray}

Taking the inner product of equation (005) with $boldsymbol{sigma}$ begin{equation} (mathbf{x}{'}boldsymbol{cdot}boldsymbol{sigma}) = costheta(mathbf{x}boldsymbol{cdot}boldsymbol{sigma})+(1-costheta)(mathbf{n}boldsymbol{cdot}mathbf{x})(mathbf{n}boldsymbol{cdot}boldsymbol{sigma})+sintheta[(mathbf{n}boldsymbol{times}mathbf{x})boldsymbol{cdot}boldsymbol{sigma})] tag{007} end{equation}
we have begin{equation} X{'} = costheta ;X+(1-costheta)(mathbf{n}boldsymbol{cdot}mathbf{x})N+sintheta[(mathbf{n}boldsymbol{times}mathbf{x})boldsymbol{cdot}boldsymbol{sigma})] tag{008} end{equation} where begin{equation} N equiv mathbf{n}boldsymbol{cdot}boldsymbol{sigma}= begin{bmatrix} n_3&n_1-in_2 n_1+in_2&-n_3 end{bmatrix} tag{009} end{equation}

After a not so easy elaboration equation (008) turns to be begin{equation} X{'}=left[Icosfrac{theta}{2}-i(mathbf{n} boldsymbol{cdot} boldsymbol{sigma})sinfrac{theta}{2} right];X;left[Icosfrac{theta}{2}+i(mathbf{n}boldsymbol{cdot}boldsymbol{sigma})sinfrac{theta}{2} right] tag{010} end{equation} and in compact form begin{equation} X{'}=U;X;U^{boldsymbol{*}} tag{011} end{equation} where
begin{equation} Uequiv cosfrac{theta}{2}-i(mathbf{n} boldsymbol{cdot} boldsymbol{sigma})sinfrac{theta}{2} tag{012} end{equation} with hermitian conjugate begin{equation} U^{boldsymbol{*}}=Icosfrac{theta}{2}+i(mathbf{n} boldsymbol{cdot} boldsymbol{sigma})sinfrac{theta}{2} tag{013} end{equation} We choose the $2 times 2$ complex matrix $U$ to represent the rotation (005).

Now, because of the identity begin{equation} (mathbf{n} boldsymbol{cdot} boldsymbol{sigma})^2=left|mathbf{n}right|^{2} I=I tag{014} end{equation} we have begin{equation} UU^{boldsymbol{*}}=I=U^{boldsymbol{*}}U tag{015} end{equation} Operators with this property are called unitary operators, symbol $:mathrm{U}(2):$ for our case, and in general $:mathrm{U}(n):$ for $n times n$ complex matrices. Any unitary matrix $:U:$ has as determinant a unit complex number $:det(U)=e^{iphi}, phi in mathbb{R}:$.

An explicit expression of $U$ in (012) is begin{equation} U= begin{bmatrix} cosfrac{theta}{2}-isinfrac{theta}{2}n_{3} & & -sinfrac{theta}{2}left( n_{2}+in_{1}right) sinfrac{theta}{2}left( n_{2}-in_{1}right) & & cosfrac{theta}{2}+isinfrac{theta}{2}n_{3} end{bmatrix} = begin{bmatrix} alpha & beta -beta^{boldsymbol{*}} & alpha^{boldsymbol{*}} end{bmatrix} tag{016} end{equation}
where here begin{equation} alpha =cosfrac{theta}{2}-isinfrac{theta}{2}n_{3} qquad beta=-sinfrac{theta}{2}left( n_{2}+in_{1}right) tag{017} end{equation} but more generally $left(alpha,beta right)$ any pair of complex numbers satisfying the condition begin{equation} alpha alpha^{boldsymbol{*}}+betabeta^{boldsymbol{*}}=left|alpharight|^2 + left|betaright|^2=1 tag{018} end{equation} So, the unitary matrix $:U:$ in (012) has as determinant the real positive unit $:det(U)=+1:$. Unitary matrices with $:det(U)=+1:$ are called special unitary and the set symbol is $:mathrm{SU}(n):$ in general. So for the unitary matrix $:U:$ in (012) we have $:U in mathrm{SU}(2):$.

2. Quaternions-Rotations

The unitary matrix representation (016) is simplified if we define the following quantities begin{align} mathbf{1} & equiv I = begin{bmatrix} 1&0 0&1 end{bmatrix} tag{019a} mathbf{i} & equiv -isigma_{1} = begin{bmatrix} 0&-i -i&0 end{bmatrix} tag{019b} mathbf{j} & equiv -isigma_{2} = begin{bmatrix} 0&-1 1&0 end{bmatrix} tag{019c} mathbf{k} & equiv -isigma_{3} = begin{bmatrix} -i&0 0&i end{bmatrix} tag{019d} end{align}

with properties begin{equation} mathbf{i}^{2}=mathbf{j}^{2}=mathbf{k}^{2}=-mathbf{1} tag{020} end{equation} begin{equation} mathbf{i} cdot mathbf{j}=mathbf{k}=-mathbf{j}cdot mathbf{i} quad , quad mathbf{j} cdot mathbf{k}=mathbf{i}=-mathbf{k}cdot mathbf{j} quad , quad mathbf{k} cdot mathbf{i}=mathbf{j}=-mathbf{i}cdot mathbf{k} tag{021} end{equation} begin{equation} mathbf{i} cdot mathbf{j}cdot mathbf{k}= -mathbf{1} tag{022} end{equation}

Then begin{equation} U= left(cosfrac{theta}{2}right)mathbf{1}+left(n_{1}sinfrac{theta}{2}right)mathbf{i}+left(n_{2}sinfrac{theta}{2}right)mathbf{j}+left(n_{3}sinfrac{theta}{2}right)mathbf{k} tag{023} end{equation} and setting begin{equation} cosfrac{theta}{2}equiv q_{0}quad , quad n_{1}sinfrac{theta}{2} equiv q_{1} quad , quad n_{2}sinfrac{theta}{2} equiv q_{2} quad , quad n_{3}sinfrac{theta}{3} equiv q_{3} tag{024} end{equation} we have begin{equation} U= q_{0}mathbf{1}+ q_{1}mathbf{i}+q_{2}mathbf{j}+q_{3}mathbf{k} quad , quad q_{kappa}in mathbb{R}quad , quad q_{0}^{2}+q_{1}^{2}+q_{2}^{2}+q_{3}^{2}=1 tag{025} end{equation} Inversely, an expression $ U $ defined by (025) represents a rotation with parameters
$ mathbf{n},theta $ determined by equations (024).

If in equation (012) we replace $theta$ by $-theta$ or exclusively $mathbf{n}$ by $-mathbf{n}$, then we have the inverse rotation begin{equation} U^{-1}= Icosfrac{theta}{2}+i(mathbf{n} boldsymbol{cdot}boldsymbol{sigma})sinfrac{theta}{2}equiv U^{boldsymbol{*}} tag{026} end{equation} and so begin{equation} U^{-1}=U^{boldsymbol{*}}= q_{0}mathbf{1}-q_{1}mathbf{i}-q_{2}mathbf{j}-q_{3}mathbf{k} quad , quad q_{kappa}in mathbb{R}quad , quad q_{0}^{2}+q_{1}^{2}+q_{2}^{2}+q_{3}^{2}=1 tag{027} end{equation} Ignoring the condition begin{equation} q_{0}^{2}+q_{1}^{2}+q_{2}^{2}+q_{3}^{2}=1 tag{028} end{equation} we define the so called quaternions by begin{equation} boldsymbol{mathsf{Q}}= q_{0}mathbf{1}+ q_{1}mathbf{i}+q_{2}mathbf{j}+q_{3}mathbf{k} quad , quad q_{kappa}in mathbb{R} tag{029} end{equation} In analogy to the properties of complex numbers begin{equation} z=a+ib quad , quad z^{boldsymbol{*}}=text{conjugate of } z =a-ib quad , quad Vert z Vert ^{2}=zz^{boldsymbol{*}}=a^{2}+b^{2} tag{030} end{equation}
we define the conjugate of quaternion $boldsymbol{mathsf{Q}}$ to be begin{equation} boldsymbol{mathsf{Q}}^{boldsymbol{*}}= q_{0}mathbf{1}- q_{1}mathbf{i}-q_{2}mathbf{j}-q_{3}mathbf{k} tag{031} end{equation} but since, making use of properties (020) and (021), the expression $boldsymbol{mathsf{Q}}boldsymbol{mathsf{Q}}^{boldsymbol{*}}$ in not a number but a scalar multiple of the identity quaternion begin{equation} boldsymbol{mathsf{Q}}boldsymbol{mathsf{Q}}^{boldsymbol{*}}= left( q_{0}mathbf{1}+q_{1}mathbf{i}+q_{2}mathbf{j}+q_{3}mathbf{k}right) left( q_{0}mathbf{1}- q_{1}mathbf{i}-q_{2}mathbf{j}-q_{3}mathbf{k}right)=left( q_{0}^{2}+q_{1}^{2}+q_{2}^{2}+q_{3}^{2}right) mathbf{1} tag{032} end{equation} we define the norm of quaternion $boldsymbol{mathsf{Q}}$ of (029) to be
begin{equation} Vert boldsymbol{mathsf{Q}} Vert ^{2}=q_{0}^{2}+q_{1}^{2}+q_{2}^{2}+q_{3}^{2} tag{033} end{equation} As the space of complex numbers begin{equation} mathbb{C} equiv lbrace z: z=a+ib quad a,b in mathbb{R}rbrace tag{034} end{equation} is in many respects identical to the 2-dimensional real space $mathbb{R}^{boldsymbol{2}}$, so the space of quaternions begin{equation} mathcal{Q} equiv lbrace boldsymbol{mathsf{Q}}:boldsymbol{mathsf{Q}}= q_{0}mathbf{1}+ q_{1}mathbf{i}+q_{2}mathbf{j}+q_{3}mathbf{k} ; , ; q_{kappa}in mathbb{R}rbrace tag{035} end{equation} is identical to the 4-dimensional real space $mathbb{R}^{boldsymbol{4}}$.

A quaternion of unit norm begin{equation} boldsymbol{mathsf{Q}}= q_{0}mathbf{1}+ q_{1}mathbf{i}+q_{2}mathbf{j}+q_{3}mathbf{k} ; , ;q_{kappa}in mathbb{R} ; ,; Vert boldsymbol{mathsf{Q}} Vert ^{2}=q_{0}^{2}+q_{1}^{2}+q_{2}^{2}+q_{3}^{2}=1 tag{036} end{equation} or any quaternion normalized,$;boldsymbol{mathsf{Q}}/Vert boldsymbol{mathsf{Q}} Vert;$, represents a unique rotation in the 3-dimensional real space $mathbb{R}^{boldsymbol{3}}$, but inversely to any rotation corresponds a pair $; lbraceboldsymbol{mathsf{Q}},-boldsymbol{mathsf{Q}}rbrace; $, where $;boldsymbol{mathsf{Q}};$ is a unit norm quaternion.

Let the quaternions $;boldsymbol{mathsf{Q}},boldsymbol{mathsf{P}} in mathcal{Q}$ begin{equation} boldsymbol{mathsf{Q}}= q_{0}mathbf{1}+ q_{1}mathbf{i}+q_{2}mathbf{j}+q_{3}mathbf{k} quad , quad boldsymbol{mathsf{P}}= p_{0}mathbf{1}+ p_{1}mathbf{i}+p_{2}mathbf{j}+p_{3}mathbf{k} tag{037} end{equation} Using properties (020) and 021) their product is ^{begin{equation} boldsymbol{mathsf{P}}boldsymbol{mathsf{Q}}= left( p_{0}mathbf{1}+ p_{1}mathbf{i}+p_{2}mathbf{j}+p_{3}mathbf{k}right)left( q_{0}mathbf{1}+q_{1}mathbf{i}+q_{2}mathbf{j}+q_{3}mathbf{k}right) = h_{0}mathbf{1}+h_{1}mathbf{i}+h_{2}mathbf{j}+h_{3}mathbf{k}=boldsymbol{mathsf{H}} tag{038} end{equation}} where begin{align} h_{0} & = q_{0}p_{0}-left(mathbf{q} boldsymbol{cdot} mathbf{p}right) tag{039a} mathbf{h} & = p_{0}mathbf{q} +q_{0}mathbf{p}- left(mathbf{q} boldsymbol{times} mathbf{p}right) tag{039b} end{align}
and $;mathbf{q},mathbf{p},mathbf{h} in mathbb{R}^{boldsymbol{3}};$ the 3-dimensional real vectors begin{equation} mathbf{q}= left[q_{1},q_{2},q_{3}right] quad , quad mathbf{p}= left[p_{1},p_{2},p_{3}right] quad , quad mathbf{h}= left[h_{1},h_{2},h_{3}right] tag{040} end{equation}

Note that begin{equation} boldsymbol{mathsf{H}}=boldsymbol{mathsf{P}}boldsymbol{mathsf{Q}}Longrightarrow Vertboldsymbol{mathsf{H}}Vert ^{2}=Vertboldsymbol{mathsf{P}}Vert ^{2}Vertboldsymbol{mathsf{Q}}Vert ^{2} tag{041} end{equation}

If both quaternions $;boldsymbol{mathsf{Q}},boldsymbol{mathsf{P}};$ are of unit norm, $;Vertboldsymbol{mathsf{Q}}Vert ^{2}=1=Vert boldsymbol{mathsf{P}}Vert^{2};$, then they represent rotations in $;mathbb{R}^{boldsymbol{3}}$ and $;boldsymbol{mathsf{H}};$ is of unit norm also,$;Vertboldsymbol{mathsf{H}}Vert ^{2}=1;$, representing their composed rotation. In this case equations (039a) and (039b) are identical to (043a) and (043b) respectively, see 3. Addendum, under the following substitutions
begin{align} q_{0} & = cosfrac{alpha}{2} & mathbf{q}& = sinfrac{alpha}{2}mathbf{n}_alpha tag{42a} p_{0} & = cosfrac{beta}{2} & mathbf{p}& = sinfrac{beta}{2}mathbf{n}_beta tag{42b} h_{0} & = cosfrac{phi}{2} & mathbf{h}& = sinfrac{phi}{2}mathbf{n} tag{42c} end{align}

---

3. Addendum

In above Figure it's shown the rotation $U(mathbf{n}_phi,phi)$, composition of two rotations $U(mathbf{n}_alpha,alpha)$ and $U(mathbf{n}_beta,beta)$ applied in this sequence. Note that this composed rotation is determined by the following equations ^{begin{equation} cosfrac{phi}{2}=cosfrac{alpha}{2}cosfrac{beta}{2}-left(mathbf{n}_alpha boldsymbol{cdot} mathbf{n}_betaright)sinfrac{alpha}{2}sinfrac{beta}{2}=cosfrac{alpha}{2}cosfrac{beta}{2}-cosomegasinfrac{alpha}{2}sinfrac{beta}{2} tag{043a} end{equation} begin{equation} sinfrac{phi}{2} mathbf{n}_{phi}= sinfrac{alpha}{2}cosfrac{beta}{2} mathbf{n}_alpha+sinfrac{beta}{2}cosfrac{alpha}{2} mathbf{n}_beta-sinfrac{alpha}{2}sinfrac{beta}{2}left(mathbf{n}_alpha boldsymbol{times} mathbf{n}_betaright) tag{043b} end{equation}}

---

^{(1) See my answer here as user82794 Construction of Pauli Matrices}

^{(2) See my answer here Rotation of a vector}

Summary: The Pauli matrices span the vector space of $2times2$ traceless Hermitian matrices, the unit quaternions span the Lie algebra $2times2$ traceless skew-Hermitian matrices (whence the $i$ factor), the latter being the Lie algebra of the Lie group of rotations and that group's universal cover.

As a particular basis for the lie algebra, the unit quaternions' defining relationships follow (1) from their being orthonormal with respect to the unique (up to a scaling constant) billinear form (i.e. here inner product) that is invariant under the action of $SO(3)$ (this gives the commutator relationships) and (2) requiring a $2times2$ matrix representation of the Lie algebra. The Pauli matrices are then these elements divided by $i$, and you can think of the factor $i$ as being what is needed to make the Hermitian matrices close under a Lie bracket.

---

The Pauli matrices span the vector space of traceless, $2times 2$ Hermitian matrices and the quaternion units span the vector space of traceless, skew-Hermitian matrices, when we think of the faithful matrix representation of the vector space of so-called pure quaternions of the form $alpha,mathbf{i}+beta,mathbf{j}+gamma,mathbf{k}$ where $alpha,,beta,,gammainmathbb{R}$. Therefore, the Pauli matrices must be linear superpositions of the scalar $iinmathbb{C}$ times the quaternion units $mathbf{i},, mathbf{j},,mathbf{k}$ and indeed are related by a simple scaling by $i$ as you have found and as QMechanic's Answer elaborates with much detail.

The question is then: what has all this to do with rotation and spin?

To answer this, let's work in the other direction from the other answers by beginning with the group $SO(3)$ of rotations and seeing how the quaternions and Pauis follow from the geometry.

$SO(3)$, the group of proper, homogeneous isometries of Euclidean 3-space obviously captures the geometry of rotations. From this starting point, we need the following Lie group theoretic facts to go forward:

1. $SU(2)$ is the double cover, indeed the universal cover of $SO(3)$;

2. $SO(3)$ is the image of $SU(2)$ under the Adjoint representation $mathrm{Ad}:SU(2)to SO(3)$ whereby $SU(2)$ acts on its own Lie algebra $mathfrak{su}(2)congmathfrak{so}(3)$ by the spinor map. Thus $SU(2)$ acts on $mathfrak{su}(2)$ by rigid rotations;

3. $SO(3)$ and $SU(2)$ are simple Lie groups, containing no normal proper Lie subgroups of dimension 1 or greater and thus $mathfrak{su}(2)congmathfrak{so}(3)$ is a simple Lie algebra;

4. The Killing form for the Lie algebra $mathfrak{g}$ of a simple Lie group $mathfrak{G}$ is the unique billinear form, to within a real scaling constant, that is $mathrm{Ad}$-invariant. That is, if $X,,Y inmathfrak{g}$ and $F(X,,Y)$ is the value of a billinear form defined on the Lie algebra, and if further $F(X,,Y)=F(mathrm{Ad}(gamma),X, , mathrm{Ad}(gamma),Y);forall X,,Yinmathfrak{g},,forall gammainmathfrak{G}$, then $F$ is the Killing form, to within a real scaling constant.

So now the question is, what's the simplest basis for the Lie algebra $mathfrak{su}(2)$ that exponentiates to our rotation group and its double cover?

$SU(2)$ acts, through the Adjoint representation, on $mathfrak{su}(2)$ by rigid rotations, so let's define an inner product on $mathfrak{su}(2)$ that is left invariant by the group action on the algebra. This inner product will then yield the correct angles between vectors in $mathfrak{su}(2)$ since this is what is left invariant by rotations.

But we know $SU(2),,SO(3)$ are simple, so there is only one billinear form that is invariant in this way to within a scaling constant: the Killing form.

So now, it is obviously convenient to construct a basis that is orthonormal with respect to the Killing form.

With a little bit of work, one can show that such an orthonormal basis for $mathfrak{su}(2)$ must fulfill $[mathbf{i},,mathbf{j}]=2,mathbf{k}$, $[mathbf{k},,mathbf{i}]=2,mathbf{j}$, $[mathbf{j},,mathbf{k}]=2,mathbf{i}$. Furthermore, if we use $2times2$ unitary matrices to represent $SU(2)$, then such an orthonormal basis for $mathfrak{su}(2)$ must also fulfill the famous "Broughham bridge" relationships (equation (1) in QMechanic's Answer).

So now our quaternion units are defined, to within a similarity transformation. One obvious choice is the one that gives the relationship between the Paulis and the quaternion units that you cite.

For a very easy and quick explanation, use a search engine to find the paper titled "Imaginary Numbers are not Real". Eight pages in you get your answer.

Pauli matrices are mathematically and physically obtuse, a poor way to teach anything, let alone quantum mechanics. In a nutshell, quaternions are the proper way to think about spin in 3D. Imagine two objects spinning in two separate planes perpendicular to each other in that higher dimensional space. Let one "multiply" the other. Literally the first element in the operation, acts on the second, moving it into the 3rd plane that defines the space while maintaining its spin orientation. The operation is anti-commutative: reverse the order of the elements being multiplied, and you end up in that 3rd plane, but with spin in the opposite direction.

That is how spinning particles interact. It's a function of the geometry of space, not a property of the particles.

All the answers above are very good and illuminating mathematically. I like to add a simple physical picture. See this article and quote below; http://wwwf.imperial.ac.uk/~jdg/nonlinquat.pdf ‘’A quaternion can be constructed from a scalar s and a 3-vector r by forming the tetrad q = [s,r] that is defined by q = [s, r] = sI − r • σ, where r• σ = Sum(ri σi); i=1-3, and I is the 2 × 2 unit matrix. {σ1, σ2, σ3} are the Pauli spin matrices σ1 = [0 1; i 0] , σ2 = [0 1; −1 0] , σ3 =[ i 0; 0 –i] , that obey the relations; σiσj=−δij I – εijk σk . A multiplication rule between two tetrads q1 = [s1, r1] and q2 = [s2, r2] can easily be determined from these properties q1 * q2 = [s1s2 − r1 • r2, s1r2 + s2r1 + r1 × r2].’’

The article is in fluids, so what is the physics connecting the two. Radiation or fields in general follow hyperbolic type equations- the wave equation. This equation represents advancing in space and time and doesn't allow any back flow/vortices/curvature/spin. To be able to have flow in the reverse direction we need an elliptic type equation. This can represent vortices and spin. Dirac equation is elliptic while Klein-Gordon is hyperbolic. This is where we need the rotation matrices/the Pauli matrices or quaternions. In the flow of fluids, the wall does the job of deflecting the flow and creating vortices and turbulence. In the present case, it is done by the imaginary 'i', since multiplying by i rotates a quantity by 90 degrees as we know.

There is a simple answer to your question. First, in $mathrm{Cl3}$ (geometric algebra of 3D Euclidean vector space), the elements of the even part of the algebra are just quaternions.

Second, the simplest representation of orthonormal vectors in $mathrm{Cl3}$ are the Pauli matrices. In addition, we have a geometric meaning of both the Pauli matrices and quaternions.

The $σ$ matrices generate a Clifford algebra, like the Dirac $γ$ matrices. In this case it is the three-dimensional algebra of signature $+{+}+$.

In odd dimensions, the volume element/pseudoscalar of a Clifford algebra (in this case $mathbf{xyz}$ or $mathbf e_1 mathbf e_2 mathbf e_3$ or $σ_xσ_yσ_z$) commutes with everything. If in addition it squares to $-1$, which in this signature it does, then it can be identified with the complex $i$. Note that $σ_xσ_yσ_z = mathrm{diag}(i,i)$. The $σ$ matrices therefore generate the real Clifford algebra only, and $i$ represents the Clifford pseudoscalar and not the separate $i$ that would appear in a complexified algebra. This is different from the Dirac case: in 3+1 dimensions, while the pseudoscalar still squares to $-1$, it anticommutes with the odd part of the algebra, and so can't be identified with $i$.

The even subalgebra of any Clifford algebra contains a representation of the spin group (and the full algebra contains a representation of the pin group). In signature $+{+}+$ and $-{-}-$, the even subalgebra is isomorphic to the quaternions, and elements of norm $1$ represent elements of $mathrm{Spin}(3)$ in the same way that unit quaternions do.

The basis quaternions can be written as pairwise products of $σ$ matrices or as products of single $σ$ matrices with the pseudoscalar. But note that you can't take $mathbf i, mathbf j, mathbf k = iσ_x, iσ_y, iσ_z$, because then you get $mathbf{ijk}=+1$, which is fine in the abstract but isn't what Hamilton chose. You need $mathbf i, mathbf j, mathbf k = -iσ_x, -iσ_y, -iσ_z$ instead. (You could also take any odd permutation of $iσ_x, iσ_y, iσ_z$, but that makes no geometric sense if you're associating the $σ$ matrices and basis quaternions with coordinate axes, as you should.)

---

You can alternately identify the $σ$ matrices with the generating set $pmmathbf{tx}, pmmathbf{ty}, pmmathbf{tz}$ of the even, real Clifford algebra of 3+1 dimensions, where the sign should be chosen so that $σ_xσ_yσ_z$ is the pseudoscalar. Because this is only the even subalgebra, the pseudoscalar is central as before and it can be (and is) identified with $i$. The quaternions appear as the subalgebra of rotations in the hypersurface normal to $mathbf t$. The 6-dimensional subalgebra isomorphic to $mathrm{Spin}(3,1)$ turns out to be $mathrm{SL}(2,mathbb C)subset mathrm M_2(mathbb C)$.