I wrote a program; I like to know what I'm getting into before trying to prove things. I can already suggest that it appears all Hurwitz quaternions are expressible as $$ q = u (r^2 + s^2), $$ where $q,u,r,s$ are Hurwitz quaternions and $u$ is a unit (there are 24 of those). For any fixed norm, so far at least one out of six such quaternions are the sum of two squares, and multiplication on the left by a unit gives everything. The worst behavior is norm 21, only 136 out of 768 Hurwitz quaternions of norm 21 are the sum of two squares. About the sizes of things, in finding the sums of two squares (before multiplying by any unit): so far it has sufficed to take $r,s$ with norms less than double the norm of $q,$ so the norms of their squares are less that 4 times the square of the norm of $q.$ The extremes for this seem to occur when the norm of $q$ is already the square of a prime $p \equiv 3 \pmod 4.$ So, two examples are $q = r^2 + s^2$ with $$ q = \frac{ 5}{2 } - \frac{ 3 }{2 } i - \frac{ 1 }{2 } j - \frac{ 1 }{2 } k, \; \; r^2= \frac{ 19}{2 } + \frac{ 21 }{2 } i + \frac{ 7 }{2 } j + \frac{ 7 }{2 } k, \; \; s^2 = -7 - 12 i - 4 j - 4 k $$ with norms $9,225,225,$ and $225 / 81 \approx 2.777$ Next $q = r^2 + s^2$ with $$ q = \frac{ 13}{2 } - \frac{ 5 }{2 } i - \frac{ 1 }{2 } j - \frac{ 1 }{2 } k, \; $$ $$ r^2= -\frac{ 61}{2 } - \frac{ 165 }{2 } i - \frac{ 33 }{2 } j - \frac{ 33 }{2 } k, \; s^2 = 37 +80 i +16 j +16 k $$ with norms $49,8281,8281,$ and $8281 / 2401 \approx 3.449$ EDIT, Sunday: I ran norm 121 overnight out of curiousity, with bound $6 n^2,$ which may or may not have really been large enough to correctly count the two squares. An extreme was $q = r^2 + s^2$ with $$ q = \frac{ 7}{2 } + \frac{ 11 }{2 } i + \frac{ 1 }{2 } j - \frac{ 17 }{2 } k, \; $$ $$ r^2= \frac{ 407}{2 } + \frac{ 155 }{2 } i + \frac{ 341 }{2 } j + \frac{ 31 }{2 } k, \; \; \; s^2 = -200 -72 i -168 j -24 k $$ with norms $121,76729,73984,$ and $76729 / 14641 \approx 5.241$ If need be I can find out what the squares are squares of. Part of a big speed improvement was dropping that printout. Alright, sssssatistics. As I said, before multiplying on the left by the 24 units, the sums of two squares are not usually all items of that norm. For example, in norm $1,$ the six Hurwitz quaternions $\pm i, \pm j, \pm k$ are not the sum of two squares. Not my fault. norm two squares not total 1 18 6 24 2 6 18 24 3 68 28 96 4 24 0 24 5 84 60 144 6 24 72 96 7 144 48 192 8 18 6 24 9 162 150 312 10 42 102 144 11 180 108 288 12 88 8 96 13 228 108 336 14 48 144 192 15 432 144 576 16 24 0 24 17 180 252 432 18 84 228 312 19 392 88 480 20 120 24 144 21 136 632 768 22 120 168 288 23 432 144 576 24 96 0 96 25 390 354 744 26 90 246 336 27 724 236 960 28 184 8 192 29 276 444 720 30 192 384 576 31 600 168 768 32 24 0 24 33 564 588 1152 34 114 318 432 35 864 288 1152 36 240 72 312 37 564 348 912 38 168 312 480 39 848 496 1344 40 138 6 144 41 588 420 1008 42 192 576 768 43 792 264 1056 44 264 24 288 45 900 972 1872 46 216 360 576 47 936 216 1152 48 88 8 96 49 642 726 1368 ....................................................