Let $\mathbb{H}$ denote the skew-field of the quaternions.
I am seeking for a proof that norm perserving $\mathbb{H}$-linear maps are precisely the symplectic linear maps, i.e. those preserving the symplectic form $$\langle h,h'\rangle = \sum_{i=1}^n h_i \bar{h}'_i,$$ meaning $A:\mathbb{H}^n\rightarrow \mathbb{H}^n$ satisfies $\langle Ah, Ah'\rangle = \langle h,h'\rangle$ for all $h,h'\in\mathbb{H}^n$ In the real and complex case this is achieved by the polarization identity. This would be solved by relating this symplectic form to the quaternionic norm. However, I cannot find any reference and I am lost trying to cancel terms: non-commutativity really hardens the calculation.
Note: by quaternionic norm I mean $\|h\| = \sqrt{\langle h, h\rangle}$, which agrees with the Euclidean norm under the identification $\mathbb{H}^n \simeq \mathbb{R}^{4n}$.
My preference is to apply matrices on the left of vectors, and hence apply scalars on the right, and hence my sesquilinear form $\langle u,v\rangle=\overline{u_1}v_1+\cdots+\overline{u_n}v_n$ is conjugate-linear in the first argument so that scalars slide according to the relation $\langle u\alpha,v\beta\rangle=\overline{\alpha}\langle u,v\rangle\beta$.
Expanding both sides of the identity $\|Au+Av\|^2=\|u+v\|^2$ yields
$$ \|Au\|^2+2\mathrm{Re}\langle Au,Av\rangle+\|Av\|^2=\|u\|^2+2\mathrm{Re}\langle u,v\rangle+\|v\|^2, $$
which implies $\mathrm{Re}\langle Au,Av\rangle=\mathrm{Re}\langle u,v\rangle$. IOW, $A$ preserves the real inner product $\mathrm{Re}\langle u,v\rangle$.
Define the difference $\Delta=\langle Au,Av\rangle-\langle u,v\rangle$. We have $\mathrm{Re}(\Delta)=0$. If we substitute $u\mapsto u\alpha$, $v\mapsto v\beta$ this becomes $\mathrm{Re}(\overline{\alpha}\Delta\beta)=0$. By picking $\alpha$ and $\beta$ appropriately, we can arrange for $\overline{\alpha}\Delta\beta=|\Delta|$, so we conclude $|\Delta|=0$, hence $\Delta=0$, or in other words $\langle Au,Av\rangle=\langle u,v\rangle$.
This trick is called "arbitrage."