Why is there this relationship between quaternions and Pauli matrices?

1. At the level of formulas, the three quaternionic units $i_a$, $a\in~\{1,2,3\}$, in $\mathbb{H}\cong \mathbb{R}^4$ satisfy $$i_a i_b ~=~ -\delta_{ab} + \sum_{c=1}^3\varepsilon_{abc} i_c, \qquad\qquad a,b~\in~\{1,2,3\}, \tag{1}$$ while the three Pauli matrices $\sigma_a \in {\rm Mat}_{2\times 2}(\mathbb{C})$, $a\in~\{1,2,3\}$, $\mathbb{C}=\mathbb{R}+\mathrm{i}\mathbb{R}$, satisfy $$\sigma_a \sigma_b ~=~ \delta_{ab} {\bf 1}_{2\times 2} + \mathrm{i}\sum_{c=1}^3\varepsilon_{abc} \sigma_c\quad\Leftrightarrow \quad \sigma_{4-a} \sigma_{4-b} ~=~ \delta_{ab} {\bf 1}_{2\times 2} - \mathrm{i}\sum_{c=1}^3\varepsilon_{abc} \sigma_{4-c}, $$ $$ \qquad\qquad a,b~\in~\{1,2,3\},\tag{2}$$ with complex unit $\mathrm{i}\in\mathbb{C}.$ In other words, we evidently have an $\mathbb{R}$-algebra monomorphism $$\Phi:~~\mathbb{H}~~\longrightarrow ~~{\rm Mat}_{2\times 2}(\mathbb{C}).\tag{3}$$ by extending the definition $$\Phi(1)~=~{\bf 1}_{2\times 2},\qquad \Phi(i_a)~=~\mathrm{i}\sigma_{4-a}, \qquad\qquad a~\in~\{1,2,3\},\tag{4}$$ via $\mathbb{R}$-linearity. This observation essentially answers OP title question (v2).

2. However OP's question touches upon many beautiful and useful mathematical facts about Lie groups and Lie algebras, some of which we would like to mention. The image of the $\mathbb{R}$-algebra monomorphism (3) is $$\Phi(\mathbb{H}) ~=~ \left\{\left. \begin{pmatrix} \alpha & \beta \cr -\bar{\beta} & \bar{\alpha} \end{pmatrix}\in {\rm Mat}_{2\times 2}(\mathbb{C}) \right| \alpha,\beta \in\mathbb{C}\right\}$$ $$~=~ \left\{ M\in {\rm Mat}_{2\times 2}(\mathbb{C}) \left| \overline{M} \sigma_2=\sigma_2 M\right. \right\}.\tag{5}$$ Let us for the rest of this answer identify $\mathrm{i}=i_1$. Then the $\mathbb{R}$-algebra monomorphism (3) becomes $$ \mathbb{C}+\mathbb{C}i_2~=~\mathbb{H}~\ni~x=x^0+\sum_{a=1}^3 i_a x^a ~=~\alpha+\beta i_2$$ $$~~\stackrel{\Phi}{\mapsto}~~ \begin{pmatrix} \alpha & \beta \cr -\bar{\beta} & \bar{\alpha} \end{pmatrix} ~=~ x^0{\bf 1}_{2\times 2}+\mathrm{i}\sum_{a=1}^3 x^a \sigma_{4-a}~\in~ {\rm Mat}_{2\times 2}(\mathbb{C}),$$ $$ \alpha~=~x^0+\mathrm{i}x^1~\in~\mathbb{C},\qquad \beta~=~x^2+\mathrm{i}x^3~\in~\mathbb{C},\qquad x^0, x^1, x^2, x^3~\in~\mathbb{R}.\tag{6}$$

3. One may show that $\Phi$ is a star algebra monomorphism, i.e. the Hermitian conjugated matrix satisfies $$ \Phi(x)^{\dagger}~=~\Phi(\bar{x}), \qquad x~\in~\mathbb{H}. \tag{7}$$ Moreover, the determinant becomes the quaternionic norm square $$\det \Phi(x)~=~ |\alpha|^2+|\beta|^2~=~\sum_{\mu=0}^3 (x^{\mu})^2 ~=~|x|^2, \qquad x~\in~\mathbb{H}.\tag{8}$$ Let us for completeness mention that the transposed matrix satisfies $$\Phi(x)^t~=~\Phi(x|_{x^2\to-x^2})~=~ \Phi(-j\bar{x}j), \qquad x~\in~\mathbb{H}. \tag{9} $$

4. Consider the Lie group of quaternionic units, which is also the Lie group $$U(1,\mathbb{H})~:=~\{x\in\mathbb{H}\mid |x|=1 \} \tag{10}$$ of unitary $1\times 1$ matrices with quaternionic entries. Eqs. (7) and (8) imply that the restriction $$\Phi_|:~U(1,\mathbb{H})~~\stackrel{\cong}{\longrightarrow}~~ SU(2)~:=~\{g\in {\rm Mat}_{2\times 2}(\mathbb{C})\mid g^{\dagger}g={\bf 1}_{2\times 2},~\det g = 1 \} $$ $$~=~\left\{\left. \begin{pmatrix} \alpha & \beta \cr -\bar{\beta} & \bar{\alpha} \end{pmatrix} \in {\rm Mat}_{2\times 2}(\mathbb{C}) \right| \alpha, \beta\in\mathbb{C}, |\alpha|^2+|\beta|^2=1\right\}\tag{11}$$ of the monomorphism (3) is a Lie group isomorphism. In other words, we have shown that

   $$ U(1,\mathbb{H})~\cong~SU(2).\tag{12}$$

5. Consider the corresponding Lie algebra of imaginary quaternionic number $$ {\rm Im}\mathbb{H}~:=~\{x\in\mathbb{H}\mid x^0=0 \}~\cong~\mathbb{R}^3 \tag{13}$$ endowed with the commutator Lie bracket. [This is (twice) the usual 3D vector cross product in disguise.] The corresponding Lie algebra isomorphism is $$\begin{align}\Phi_|:~{\rm Im}\mathbb{H}~~\stackrel{\cong}{\longrightarrow}~~ su(2)~:=~&\{m\in {\rm Mat}_{2\times 2}(\mathbb{C})\mid m^{\dagger}=-m \}\cr ~=~&\mathrm{i}~{\rm span}_{\mathbb{R}}(\sigma_1,\sigma_2,\sigma_3),\end{align}\tag{14}$$ which brings us back to the Pauli matrices. In other words, we have shown that

   $$ {\rm Im}\mathbb{H}~\cong~su(2).\tag{15}$$

6. It is now also easy to make contact to the left and right Weyl spinor representations in 4D spacetime $\mathbb{H}\cong \mathbb{R}^4$ endowed with the quaternionic norm $|\cdot|$, which has positive definite Euclidean (as opposed to Minkowski) signature, although we shall only be sketchy here. See also e.g. this Phys.SE post.

   Firstly, $U(1,\mathbb{H})\times U(1,\mathbb{H})$ is (the double cover of) the special orthogonal group $SO(4,\mathbb{R})$.

   The group representation $$\rho: U(1,\mathbb{H}) \times U(1,\mathbb{H}) \quad\to\quad SO(\mathbb{H},\mathbb{R})~\cong~ SO(4,\mathbb{R}) \tag{16}$$ is given by $$\rho(q_L,q_R)x~=~q_Lx\bar{q}_R, \qquad q_L,q_R~\in~U(1,\mathbb{H}), \qquad x~\in~\mathbb{H}. \tag{17}$$ The crucial point is that the group action (17) preserves the norm, and hence represents orthogonal transformations. See also this math.SE question.

   Secondly, $U(1,\mathbb{H})\cong SU(2)$ is (the double cover of) the special orthogonal group $SO({\rm Im}\mathbb{H},\mathbb{R})\cong SO(3,\mathbb{R})$.

   This follows via a diagonal restriction $q_L=q_R$ in eq. (17).

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

Any vector in $\mathbb{R}^3$ can be represented by a $2\times2$ 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 $2\times2$ 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,0\right),\quad \mathbf{e}_{2}=\left(0,1,0\right),\quad \mathbf{e}_{3}=\left(0,0,1\right) \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}= \cos\theta \;\mathbf{x}+(1-\cos\theta)\;(\mathbf{n}\boldsymbol{\cdot}\mathbf{x})\;\mathbf{n}+\sin\theta\;(\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_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{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}) = \cos\theta(\mathbf{x}\boldsymbol{\cdot}\boldsymbol{\sigma})+(1-\cos\theta)(\mathbf{n}\boldsymbol{\cdot}\mathbf{x})(\mathbf{n}\boldsymbol{\cdot}\boldsymbol{\sigma})+\sin\theta[(\mathbf{n}\boldsymbol{\times}\mathbf{x})\boldsymbol{\cdot}\boldsymbol{\sigma})] \tag{007} \end{equation}
we have \begin{equation} X{'} = \cos\theta \;X+(1-\cos\theta)(\mathbf{n}\boldsymbol{\cdot}\mathbf{x})N+\sin\theta[(\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[I\cos\frac{\theta}{2}-i(\mathbf{n} \boldsymbol{\cdot} \boldsymbol{\sigma})\sin\frac{\theta}{2} \right]\;X\;\left[I\cos\frac{\theta}{2}+i(\mathbf{n}\boldsymbol{\cdot}\boldsymbol{\sigma})\sin\frac{\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} U\equiv \cos\frac{\theta}{2}-i(\mathbf{n} \boldsymbol{\cdot} \boldsymbol{\sigma})\sin\frac{\theta}{2} \tag{012} \end{equation} with hermitian conjugate \begin{equation} U^{\boldsymbol{*}}=I\cos\frac{\theta}{2}+i(\mathbf{n} \boldsymbol{\cdot} \boldsymbol{\sigma})\sin\frac{\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^{i\phi}, \phi \in \mathbb{R}\:$.

An explicit expression of $U$ in (012) is \begin{equation} U= \begin{bmatrix} \cos\frac{\theta}{2}-i\sin\frac{\theta}{2}n_{3} & & -\sin\frac{\theta}{2}\left( n_{2}+in_{1}\right) \\ \sin\frac{\theta}{2}\left( n_{2}-in_{1}\right) & & \cos\frac{\theta}{2}+i\sin\frac{\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 =\cos\frac{\theta}{2}-i\sin\frac{\theta}{2}n_{3} \qquad \beta=-\sin\frac{\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{*}}+\beta\beta^{\boldsymbol{*}}=\left\|\alpha\right\|^2 + \left\|\beta\right\|^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 -i\sigma_{1} = \begin{bmatrix} 0&-i\\ -i&0 \end{bmatrix} \tag{019b}\\ \mathbf{j} & \equiv -i\sigma_{2} = \begin{bmatrix} 0&-1\\ 1&0 \end{bmatrix} \tag{019c}\\ \mathbf{k} & \equiv -i\sigma_{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(\cos\frac{\theta}{2}\right)\mathbf{1}+\left(n_{1}\sin\frac{\theta}{2}\right)\mathbf{i}+\left(n_{2}\sin\frac{\theta}{2}\right)\mathbf{j}+\left(n_{3}\sin\frac{\theta}{2}\right)\mathbf{k} \tag{023} \end{equation} and setting \begin{equation} \cos\frac{\theta}{2}\equiv q_{0}\quad , \quad n_{1}\sin\frac{\theta}{2} \equiv q_{1} \quad , \quad n_{2}\sin\frac{\theta}{2} \equiv q_{2} \quad , \quad n_{3}\sin\frac{\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}= I\cos\frac{\theta}{2}+i(\mathbf{n} \boldsymbol{\cdot}\boldsymbol{\sigma})\sin\frac{\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 $\; \lbrace\boldsymbol{\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 \Vert\boldsymbol{\mathsf{H}}\Vert ^{2}=\Vert\boldsymbol{\mathsf{P}}\Vert ^{2}\Vert\boldsymbol{\mathsf{Q}}\Vert ^{2} \tag{041} \end{equation}

If both quaternions $\;\boldsymbol{\mathsf{Q}},\boldsymbol{\mathsf{P}}\;$ are of unit norm, $\;\Vert\boldsymbol{\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,$\;\Vert\boldsymbol{\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} & = \cos\frac{\alpha}{2} & \mathbf{q}& = \sin\frac{\alpha}{2}\mathbf{n}_\alpha \tag{42a}\\ p_{0} & = \cos\frac{\beta}{2} & \mathbf{p}& = \sin\frac{\beta}{2}\mathbf{n}_\beta \tag{42b}\\ h_{0} & = \cos\frac{\phi}{2} & \mathbf{h}& = \sin\frac{\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} \cos\frac{\phi}{2}=\cos\frac{\alpha}{2}\cos\frac{\beta}{2}-\left(\mathbf{n}_\alpha \boldsymbol{\cdot} \mathbf{n}_\beta\right)\sin\frac{\alpha}{2}\sin\frac{\beta}{2}=\cos\frac{\alpha}{2}\cos\frac{\beta}{2}-\cos\omega\sin\frac{\alpha}{2}\sin\frac{\beta}{2} \tag{043a} \end{equation} \begin{equation} \sin\frac{\phi}{2}\ \mathbf{n}_{\phi}= \sin\frac{\alpha}{2}\cos\frac{\beta}{2}\ \mathbf{n}_\alpha+\sin\frac{\beta}{2}\cos\frac{\alpha}{2}\ \mathbf{n}_\beta-\sin\frac{\alpha}{2}\sin\frac{\beta}{2}\left(\mathbf{n}_\alpha \boldsymbol{\times} \mathbf{n}_\beta\right) \tag{043b} \end{equation}}

---

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

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

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.