Analytic continuation of Riemann zeta $\zeta(s)$ from the complex $\mathbb{C}$ to quaternion $\mathbb{H}$? One way to define  Riemann zeta function is by the analytic continuation of
$$\zeta(s) = \sum_{n=1}^\infty \frac{1}{n^s} = \frac{1}{1^s} + \frac{1}{2^s} + \frac{1}{3^s} + \cdots$$
for the domain $Re(s)>1$ to the full complex plane in $\mathbb{C}$.
Thus, Riemann zeta function is defined for $s \in \mathbb{C}$ and $\zeta(s)  \in  \mathbb{C}$

My question is that do we gain anything new to do analytic continuation of Riemann zeta function such that a "modified Riemann zeta function" so
$s  \in  \mathbb{H}$ is in quaternion? and $\zeta(s)  \in  \mathbb{H}?$

Does this lead to any interesting result in the math literature?
Edit: more precisely, according to the comment, we seek for an analytic continuation of $\zeta(s)$ from the complex $\mathbb{C}/\{1\}$ to quaternion $\mathbb{H}/\{1\}$?
 A: You don't gain anything. 
To extend a holomorphic function $f(z)$ of a complex variable $z=x+yi$ to a function of a quaternion variable, if its series' coefficients are real then it's just defined by
$$ f(x+yi)=u+vi \implies  f(x+y\mathbf{t})=u+v\mathbf{t} \tag{$\circ$}$$
for unit vectors $\mathbf{t}$. (Every quaternion is expressible as $x+y\mathbf{t}$ for a unit vector $\mathbf{t}$, which is unique up to the signs of $y$ and choice of $\pm\mathbf{t}$). Equivalently, $f$ extends to quaternions by "rotating" the graph in $\mathbb{C}^2$ around to get a graph in $\mathbb{H}^2$. In other words, $f(pzp^{-1})=pf(z)p^{-1}$ for complex numbers $z$ and quaternions $p$ (note every quaternion is expressible as $pzp^{-1}$ for a complex number $z$ and quaternion $p$, but not uniquely).
The reason this happens is because the unit vectors (i.e. pure imaginary unit quaternions) $\mathbf{t}$ are precisely the square roots of $-1$ in $\mathbb{H}$, so algebraically they behave just like $i$ does in $\mathbb{C}$. If you look at the Dirichlet series definition of the zeta function $\zeta(s)$ for $\mathrm{Re}(s)>1$, they involve $1/n^s$ which is computed as $\exp(-\ln(n)s)$ Euler's formula $\exp(i\theta)=\cos\theta+\sin\theta\,i$ generalizes to quaternions since it follows entirely from $i$ being a square root of negative one. The same applies to the analytic continuation of $\zeta(s)$.
Same story for octonions.
In order to get something nontrivial, you would want to to start with a power series that has complex coefficients (so, isn't simply extended from a real variable function like $\zeta(s)$ is). There is extra freedom in how you define the monomials for a function of a quaternion variable, since each $a_nz^n$ may be replaced by
$$ \square z\square\cdots\square z\square $$
where the $\square$'s are complex numbers which multiply to $a_n$ and there are $n$ $z$s present.
However, doing this will not give you differentiable functions. In fact, the limit definition
$$ f'(p)=\lim_{h\to0}\frac{f(p+h)-f(p)}{h} $$
generalizes to quaternions in two ways: a "left" derivative and a "right" derivative, depending on which side of $\Delta f$ you put $h^{-1}$ (note $h\to0$ within $\mathbb{H}$ now). This turns out to be extremely restrictive: the only left or right differentiable quaternion functions are affine functions $f(q)=qa+b$ or $f(q)=aq+b$ respectively. It's a small miracle complex differentiable yields such a rich theory.
Moreover, say you start with a holomorphic function $f$, pick two complex numbers $\alpha$ and $\beta$ such that disk of convergence of the Taylor series around $\alpha$ includes $\beta$ and vice-versa. This gives you two different series (one in $z-\alpha$ and one in $z-\beta$), and (I'm pretty sure) these almost never give you the same function of a quaternion variable!
A: $\Bbb{H}$ is just a sub-algebra of $M_2(\Bbb{C})$.
For $A \in M_n(\Bbb{C})$ use the Jordan normal form to obtain $A = P J P^{-1} = P (D+N)P^{-1}$ where $D$ is diagonal and $DN=ND$ and $N^n = 0$. Let $f(s) = (s-1)\zeta(s)= \sum_{k=0}^\infty c_k s^k$ which is entire then $$P^{-1} f(A)P  =f(D+N)=\sum_{k=0}^\infty c_k (D+N)^k =\sum_{k=0}^\infty c_k \sum_{l=0}^{n-1} {k \choose l} D^{k-l}N^l= \sum_{k=0}^{n-1} \frac{N^k}{k!} f^{(k)}(D)$$
Note the obtained function $A \mapsto f(A)$ doesn't depend on the basepoint $s_0 = 0$ we chose to expand $f(s)$ in power series.
It is not hard to convince that something similar happens with a meromorphic function such as $\zeta(s)$ obtaining $$\zeta(A) = P \zeta(D+N)P^{-1}=   P \sum_{k=0}^{n-1} \frac{N^k}{k!} \zeta^{(k)}(D)P^{-1}$$
where $\zeta^{(k)}(D)$ is the matrix of $k$-th derivatives $$\zeta^{(k)}(D) = \pmatrix{\zeta^{(k)}(D_{11}) &  & \\  &  \zeta^{(k)}(D_{22}) & \\ & & \ldots}$$
If $q \in \Bbb{H}\subset M_2(\Bbb{C})$ then $q q^* = q^* q = N(q) I$ so that $$q = P DP^{-1}, \qquad \zeta(q)  =P \zeta(D) P^{-1}$$
