Quaternionic and octonionic analogues of the Basel problem I asked this question in MSE around 3 months ago but I have received no answer yet, so following the suggestion in the comments I decided to post it here.
It is a well-known fact that
$$\sum_{0\neq n\in\mathbb{Z}} \frac{1}{n^k} = r_k (2\pi)^k$$
for any integer $k>1$, where $r_k$ are rational numbers which can be given explicitly in terms of Bernoulli numbers. For example, for $k=2$ the sum equals $\pi^2/3$ (this is essentially the Basel problem), and for $k=4$ it equals $\pi^4/45$. Note that for odd $k$ the sum vanishes.
The theory of elliptic curves with complex multiplication allows us to extend this result to systems of complex integers such as the Gaussian integers, or more generally the ring of integers in an imaginary quadratic number field of class number 1. Namely, for $k>2$ we have
$$\sum_{0\neq \lambda\in\mathbb{Z[\omega]}} \frac{1}{\lambda^k} = r_k \varpi^k,$$
where again $r_k$ are rational constants and $\varpi \in \mathbb{R}$ (the "complex $2\pi$") depends only on the ring $\mathcal{O}=\mathbb{Z[\omega]}$ and is an algebraic multiple of a so-called Chowla–Selberg period, given by a product of powers of certain gamma factors (note that the sum is always a real number since it is invariant under conjugation). For example, for the Eisenstein ($\omega = (1+\sqrt{3} i)/2$), Gaussian ($\omega = i$) and Kleinian ($\omega = (1+\sqrt{7} i)/2$) integers, we have respectively
$$\varpi_3 = 3^{-1/4} \sqrt{2\pi} \left(\frac{\Gamma(1/3)}{\Gamma(2/3)}\right)^{3/2}, \quad \varpi_4 = 4^{-1/4} \sqrt{2\pi} \left(\frac{\Gamma(1/4)}{\Gamma(3/4)}\right), \quad \varpi_7 = 7^{-1/4} \sqrt{2\pi} \left(\frac{\Gamma(1/7)\Gamma(2/7)\Gamma(4/7)}{\Gamma(3/7)\Gamma(5/7)\Gamma(6/7)}\right)^{1/2}.$$
For higher class numbers there is a similar formula, though in that case $r_k$ will in general not be rational but algebraic. A nice exposition of this result can be found in Section 6.3 of Zagier - Elliptic Modular Forms and Their Applications.

My question is whether this is still true for hypercomplex number systems, such as the Hurwitz integers or the octonionic integers.  Define $$S_k[\mathcal{O}] = \sum_{0\neq \lambda\in\mathcal{O}} \frac{1}{\lambda^k}$$ for $k>\operatorname{dim} \mathcal{O}$, where $\mathcal{O}$ is now an order in a totally definite rational quaternion/octonion algebra of class number 1. The restriction on $k$ is so that the sum converges absolutely.

Subquestion 1: Do we have $S_k[\mathcal{O}] = r_k \varpi^k$ for some rational sequence $r_k$ and some real number $\varpi$ depending only on $\mathcal{O}$ (a "quaternionic/octonionic $2\pi$")?

Obviously $\varpi$ will only be defined up to a nonzero rational factor. An equivalent question is whether $(S_m[\mathcal{O}])^n/(S_n[\mathcal{O}])^m$ is rational for any $m, n$ such that $S_n[\mathcal{O}]\neq 0$.

Subquestion 2: If so, can (some fixed choice of) $\varpi$ be expressed in terms of
known constants such as $\zeta'(-1)$ or $\zeta'(-3)$?

The reason I'm mentioning these particular constants is that in the previous cases (real and complex) the period $\varpi$ turns out to be equal to $e^{-\zeta'(\mathcal{O},0)/\zeta(\mathcal{O},0)}$ up to an algebraic factor, where the zeta function attached to the ring of integers $\mathcal{O}=\mathbb{Z}$ or $\mathbb{Z[\omega]}$ is defined as
$$\zeta(\mathcal{O},s) = \sum_{0\neq \lambda\in\mathcal{O}} |\lambda|^{-s}.$$
(This is in general not the same as the previous sums, note the absolute value.) In the case that $\mathcal{O}$ is instead a quaternionic or octonionic order, the logarithmic derivative of this zeta function at $s=0$ can be expressed in terms of $\zeta'(-1)$ or $\zeta'(-3)$ respectively, where $\zeta(s)$ is the ordinary Riemann zeta function.

I calculated a few sums numerically for the ring of Hurwitz quaternions. The result is $$S_6[\mathcal{O}] \approx 10.76,\quad S_8[\mathcal{O}] \approx 1.196,\quad S_{12}[\mathcal{O}] \approx 23.9905.$$
Unfortunately the calculations take a lot of time, and the precision is not enough to determine whether e.g. $S_{12}[\mathcal{O}]/(S_6[\mathcal{O}])^2$ is rational to any degree of confidence.
I also found the recent paper Period Relations for Quaternionic Elliptic Functions by Z. Amir-Khosravi that refers to previous works by R. Fueter (Über die Quaternionenmultiplikation regulärer vierfachperiodischer Funktionen) and R. Krausshar (Generalized Analytic Automorphic Forms in Hypercomplex Spaces). A certain $3$-parameter family of quaternionic Eisenstein-like functions associated to a lattice in $\mathbb{R}^4$ is introduced, and shown to enjoy period-like relations resembling those in the complex case. Unfortunately, the form of these functions is restricted by quaternionic regularity to contain factors of the quaternionic norm (c.f. equations (2.5)-(2.7) in the paper), and as far as I can see they aren't directly related to the sums of pure powers I'm interested in.
 A: This isn't really a full answer, but it's too long for a comment, and perhaps it's informative all the same.
Your sum $S_k[\mathcal{O}]$ can be written as the value at $s = k$ of the sum
$$\sum_{0 \ne \lambda \in \mathcal{O}} \frac{\lambda^k}{Nm(\lambda)^s} = \sum_{n \ge 1} a^{(k)}_n n^{-s},$$
where $a^{(k)}_n := \sum_{N(\lambda) = n} \lambda^k$.
Now, I claim that $\sum_{n \ge 1} a^{(k)}_n q^n$ is the $q$-expansion of a modular form -- or something slightly more general, namely a quasi-modular form [*] -- of weight $k + 2$ and some level depending on $\mathcal{O}$; for the Hurwitz integers the level is $\Gamma_0(2)$. This should follow from thinking about Brandt matrices, which are a way of computing modular forms using quaternion algebras; see e.g. this article by Kimball Martin.
Anyway, once you know what to look for, it's now quite easy to recognise the sequences $(a^{(k)}_n)_{n \ge 1}$ for small $k$. For example, when $k = 6$, what you get is exactly the $q$-expansion of $12f_8$, where $f_8$ is the unique normalised modular cusp form of weight 8 and level 2. So $S_6[\mathcal{O}]$ is a value of the $L$-series of a modular form. In fact, we have $S_6[\mathcal{O}] = 12 L(f_8, 6) = 10.758540419274832757072...$, which agrees with your computations above. Similarly, unless I've slipped up in my computations, we have
$$S_8[\mathcal{O}] = 12 \big( L(f_{10}, 8) - L(f_8, 7) \big) = 1.18636076594110...$$
where $f_{10}$ is the cusp form of weight 10. Since the periods of $f_{10}$ and $f_{8}$ have essentially nothing to do with each other, this strongly suggests that there is no tidy algebraic relation between $S_6[\mathcal{O}]$ and $S_8[\mathcal{O}]$.
[*] Quasi-modular forms are not too scary: they're exactly the ring of functions you get by starting with genuine modular forms and throwing in the function $E_2 = 1 - 24\sum \sigma(n) q^n$.

EDIT. Further numerical experiments suggest the following explicit formula: if $\mathfrak{S}(m)$ denotes the set of normalised newforms of level 2 and weight $m$, then for every $k \ge 6$ we seem to have
$$S_k[\mathcal{O}] = 12\left( \sum_{f \in \mathfrak{S}(k+2)} L(f, k) - \sum_{f \in \mathfrak{S}(k)} L(f, k-1) \right).$$
