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$
EDIT, Sunday afternoon. Every time I raise the multiple of square bound, I get new sums of two squares, which makes the previous count incorrect. I am going to just post these as they print out: the norms are correct, each coefficient is doubled, so coefficients are either all odd or all even. Divide by two to get the actual coefficients.
13 : -1 1 -5 -5 729 : -45 21 15 15 784 : 44 -20 -20 -20
4.63905
norm two squares not total
13 252 84 336
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
....................................................