$\mathbb{H}$ is a division ring that is not commutative, hence not a field: $q^{-1} = \bar q / N(q)$ for $q \ne 0$, while $ij = k$ and $ji = -k$

The additive group. Addition on $\mathbb{H}$ is defined coordinatewise from addition on $\mathbb{R}$, so it is associative and commutative, $0_{\mathbb{H}}$ is a two-sided identity, and $(-x_0,-x_1,-x_2,-x_3)$ is a two-sided additive inverse of $x$. Hence $(\mathbb{H},+,0_{\mathbb{H}})$ is an abelian group and $-x = (-x_0,-x_1,-x_2,-x_3)$.

Coefficient form of the product. For $m, p, q \in 4$ let $\varepsilon_m(p,q)$ be the real coefficient of the monomial $x_p y_q$ in the $m$-th coordinate of the defining product formula, and $\varepsilon_m(p,q) = 0$ when that monomial does not occur there; every such coefficient is $1$, $-1$ or $0$. Reading the four coordinates of the formula off one at a time, $[xy]_m = \sum_{p<4}\sum_{q<4} \varepsilon_m(p,q)\, x_p\, y_q$ for all $x, y \in \mathbb{H}$ and all $m \in 4$, the right-hand side being a sum of sixteen real numbers.

The cyclic symmetry of the table. Let $\gamma : 4 \to 4$ fix $0$ and send $1 \mapsto 2 \mapsto 3 \mapsto 1$, and let $\Gamma : \mathbb{H} \to \mathbb{H}$ be the coordinate permutation determined by $[\Gamma x]_{\gamma(p)} = x_p$, that is $\Gamma(x_0,x_1,x_2,x_3) = (x_0,x_3,x_1,x_2)$. Then $\Gamma$ is a bijection, $\Gamma(e_p) = e_{\gamma(p)}$, $\Gamma$ fixes each $\hat{\lambda}$, and $\Gamma(\lambda x) = \lambda\Gamma(x)$, all immediately from the definition of $\Gamma$.

Claim 2: from the table, $ij = k = (0,0,0,1)$ and $ji = -k = (0,0,0,-1)$. By [L1], $1 \ne 0$, so [L3] gives $1=1^2>0$; the ordered-field definition in [L2] then gives $-1<0$, hence $1 \ne -1$. Thus the two products differ, and multiplication on $\mathbb{H}$ is not commutative.

Claim 3, the norm identity. Evaluating the product formula at $y = \bar{x}$ gives coordinates $x_0x_0 - x_1(-x_1) - x_2(-x_2) - x_3(-x_3) = N(x)$, then $x_0(-x_1) + x_1x_0 + x_2(-x_3) - x_3(-x_2) = 0$, then $x_0(-x_2) + x_2x_0 + x_3(-x_1) - x_1(-x_3) = 0$, then $x_0(-x_3) + x_3x_0 + x_1(-x_2) - x_2(-x_1) = 0$; so $x\bar{x} = \widehat{N(x)}$. Evaluating it at $x := \bar{x}$, $y := x$ gives $x_0x_0 - (-x_1)x_1 - (-x_2)x_2 - (-x_3)x_3 = N(x)$ and, in the same way, $0$ in each of the other three coordinates; so $\bar{x}x = \widehat{N(x)}$ as well.

Claim 3, positivity. Let $x \ne 0_{\mathbb{H}}$; then $x_p \ne 0$ for at least one $p \in 4$. Each $x_{p'}^{2}$ with $x_{p'} \ne 0$ is positive, and each $x_{p'}^{2}$ with $x_{p'} = 0$ equals $0$; a sum in which at least one summand is positive and the rest are positive or $0$ is positive, since $P$ is closed under addition and $u + 0 = u$. Hence $N(x) > 0$, and in particular $N(x) \ne 0$.

Real scalars pass through the product. For $\lambda \in \mathbb{R}$ the formula gives $\hat{\lambda}x = x\hat{\lambda} = (\lambda x_0, \lambda x_1, \lambda x_2, \lambda x_3)$, an element we abbreviate $\lambda x$; in particular $\hat{\lambda}\hat{\mu} = \widehat{\lambda\mu}$, $1_{\mathbb{H}}x = x1_{\mathbb{H}} = x$, and $(-1)x = -x$ by step 1.1.

The coefficients are the multiplication table: $\varepsilon_m(p,q) = [e_pe_q]_m$ for all $m,p,q \in 4$. Fix $p$ and $q$ and substitute $x = e_p$, $y = e_q$ into step 1.2: then $x_{p'} = 0$ for $p' \ne p$ and $y_{q'} = 0$ for $q' \ne q$, and a real product with a factor $0$ is $0$, so every one of the sixteen summands vanishes except the one indexed by $(p,q)$, which equals $\varepsilon_m(p,q) \cdot 1 \cdot 1 = \varepsilon_m(p,q)$.

Both distributive laws hold. By step 1.2 and distributivity in $\mathbb{R}$, $[(x+x')y]_m = \sum_{p,q}\varepsilon_m(p,q)(x_p + x'_p)y_q = \sum_{p,q}\bigl(\varepsilon_m(p,q)x_py_q + \varepsilon_m(p,q)x'_py_q\bigr) = [xy]_m + [x'y]_m$, the last equality being a regrouping of a finite sum of thirty-two real terms; the same computation in the second argument gives $x(y+y') = xy + xy'$.

The nine table checks that establish $e_{\gamma(p)}e_{\gamma(q)} = \Gamma(e_pe_q)$ for all $p,q \in \{1,2,3\}$, read off the table of The quaternions $\mathbb{H}$: real quadruples with componentwise addition and an explicit multiplication formula matching the table on $1, i, j, k$: $jj = -1 = \Gamma(-1) = \Gamma(ii)$; $jk = i = \Gamma(k) = \Gamma(ij)$; $ji = -k = \Gamma(-j) = \Gamma(ik)$; $kj = -i = \Gamma(-k) = \Gamma(ji)$; $kk = -1 = \Gamma(-1) = \Gamma(jj)$; $ki = j = \Gamma(i) = \Gamma(jk)$; $ij = k = \Gamma(j) = \Gamma(ki)$; $ik = -j = \Gamma(-i) = \Gamma(kj)$; $ii = -1 = \Gamma(-1) = \Gamma(kk)$. The cases with $p = 0$ or $q = 0$ are immediate, both sides being $e_{\gamma(q)}$ or $e_{\gamma(p)}$ respectively.

Claim 5: every field is a commutative ring by [L6], and the multiplication of $\mathbb{H}$ is not commutative by step 1.4; so $\mathbb{H}$ is not a field.

Real scalars pass through a triple product too: by step 1.2 and step 2.1, $[(\lambda x)y]_m = \sum_{p,q}\varepsilon_m(p,q)(\lambda x_p)y_q = \lambda[xy]_m$ and likewise $[x(\lambda y)]_m = \lambda [xy]_m$, so $(\lambda x)y = \lambda(xy) = x(\lambda y)$. Taking $\lambda = -1$ gives $(-x)y = -(xy) = x(-y)$.

Reduction of associativity to the sixty-four basis triples. Applying step 1.2 twice and rearranging, $[(xy)z]_m = \sum_{s,r}\varepsilon_m(s,r)[xy]_sz_r = \sum_{p,q,r}\Bigl(\sum_s \varepsilon_m(s,r)\varepsilon_s(p,q)\Bigr)x_py_qz_r$, and applying step 2.2 twice, $[(e_pe_q)e_r]_m = \sum_s\varepsilon_m(s,r)[e_pe_q]_s = \sum_s\varepsilon_m(s,r)\varepsilon_s(p,q)$; hence $[(xy)z]_m = \sum_{p,q,r<4}x_py_qz_r\,[(e_pe_q)e_r]_m$. The same computation with the other bracketing gives $[x(yz)]_m = \sum_{p,q,r<4}x_py_qz_r\,[e_p(e_qe_r)]_m$. Therefore, if $(e_pe_q)e_r = e_p(e_qe_r)$ holds for all $p,q,r \in 4$, then $(xy)z = x(yz)$ for all $x,y,z \in \mathbb{H}$.

Basis triples containing the index $0$. Since $e_0 = 1_{\mathbb{H}}$ is a two-sided identity by step 2.1, each of $(e_0e_q)e_r = e_qe_r = e_0(e_qe_r)$, $(e_pe_0)e_r = e_pe_r = e_p(e_0e_r)$ and $(e_pe_q)e_0 = e_pe_q = e_p(e_qe_0)$ holds. So only the twenty-seven triples with $p,q,r \in \{1,2,3\}$ remain.

The multiplication of $\mathbb{H}$ commutes with $\Gamma$: $\Gamma(x)\Gamma(y) = \Gamma(xy)$ for all $x,y$. By step 1.2 and $[\Gamma x]_p = x_{\gamma^{-1}(p)}$, $[\Gamma(x)\Gamma(y)]_{\gamma(m)} = \sum_{p,q}\varepsilon_{\gamma(m)}(\gamma(p),\gamma(q))x_py_q$, while $[\Gamma(xy)]_{\gamma(m)} = [xy]_m = \sum_{p,q}\varepsilon_m(p,q)x_py_q$; by step 2.2 the two families of coefficients are $[e_{\gamma(p)}e_{\gamma(q)}]_{\gamma(m)}$ and $[e_pe_q]_m = [\Gamma(e_pe_q)]_{\gamma(m)}$, which agree by step 2.4.

Reduction of the twenty-seven triples to nine. Suppose $(e_pe_q)e_r = e_p(e_qe_r)$ for a triple $(p,q,r)$. Applying step 3.4, $(e_{\gamma(p)}e_{\gamma(q)})e_{\gamma(r)} = \Gamma(e_pe_q)\Gamma(e_r) = \Gamma((e_pe_q)e_r)$ and $e_{\gamma(p)}(e_{\gamma(q)}e_{\gamma(r)}) = \Gamma(e_p)\Gamma(e_qe_r) = \Gamma(e_p(e_qe_r))$, so the identity holds for $(\gamma(p),\gamma(q),\gamma(r))$ as well. Since $\gamma$ restricted to $\{1,2,3\}$ is a cycle of length three, for each $p \in \{1,2,3\}$ there is exactly one $t \in \{0,1,2\}$ with $\gamma^{t}(p) = 1$; hence every triple in $\{1,2,3\}^{3}$ is obtained by iterating $\gamma$ from a triple whose first entry is $1$, and it suffices to check the nine triples $(1,q,r)$ with $q,r \in \{1,2,3\}$.

The nine remaining checks, using the table of The quaternions $\mathbb{H}$: real quadruples with componentwise addition and an explicit multiplication formula matching the table on $1, i, j, k$ and the sign rule of step 3.1: $(ii)i = (-1)i = -i$ and $i(ii) = i(-1) = -i$; $(ii)j = -j$ and $i(ij) = ik = -j$; $(ii)k = -k$ and $i(ik) = i(-j) = -(ij) = -k$; $(ij)i = ki = j$ and $i(ji) = i(-k) = -(ik) = j$; $(ij)j = kj = -i$ and $i(jj) = i(-1) = -i$; $(ij)k = kk = -1$ and $i(jk) = ii = -1$; $(ik)i = (-j)i = -(ji) = k$ and $i(ki) = ij = k$; $(ik)j = (-j)j = -(jj) = 1$ and $i(kj) = i(-i) = -(ii) = 1$; $(ik)k = (-j)k = -(jk) = -i$ and $i(kk) = i(-1) = -i$. All nine agree.

Multiplication on $\mathbb{H}$ is associative: by steps 3.3, 4.1 and 4.2 the identity $(e_pe_q)e_r = e_p(e_qe_r)$ holds for all sixty-four basis triples, and step 3.2 transfers it to all of $\mathbb{H}$.

Claim 1: by step 1.1 the additive structure is an abelian group; by step 5.1 and step 2.1 multiplication is associative with two-sided identity $1_{\mathbb{H}}$, so $(\mathbb{H},\cdot,1_{\mathbb{H}})$ is a monoid; and both distributive laws hold by step 2.3. So $\mathbb{H}$ is a ring.

Claim 4. First $1_{\mathbb{H}} \ne 0_{\mathbb{H}}$, because $1 \ne 0$ in $\mathbb{R}$. Let $x \ne 0_{\mathbb{H}}$ and put $\lambda := N(x)^{-1}$, which exists by step 1.6, and $y := \lambda\bar{x} = \widehat{N(x)^{-1}}\bar{x}$. Then $xy = x(\lambda\bar{x}) = \lambda(x\bar{x}) = \lambda\widehat{N(x)} = \widehat{\lambda N(x)} = \hat{1} = 1_{\mathbb{H}}$ by step 3.1, step 1.5 and step 2.1, and $yx = (\lambda\bar{x})x = \lambda(\bar{x}x) = 1_{\mathbb{H}}$ in the same way. So $x$ is a unit of the ring $\mathbb{H}$ with $x^{-1} = \widehat{N(x)^{-1}}\bar{x}$, and $\mathbb{H}$ is a division ring; by [L5] its units form a group, and by the description of a division ring that group is $\mathbb{H}\setminus\{0_{\mathbb{H}}\}$.

Claims 1 to 5 are established: claim 1 in step 6.1, claim 2 in step 1.4, claim 3 in step 1.5 together with step 1.6, claim 4 in step 7.1 and claim 5 in step 2.5.