# Computing residues of matrix valued functions symbolically

I have a somewhat specific application in mind, but I've run into this kind of problem a few times and never found an elegant solution.

Is there an efficient way to get Mathematica to recognize that some variables in a symbolic computation are non-commutative (matrices) while others are commutative (scalars) while still being able to do elementary calculus operations on the commutative variables?

Specifically, I have expansions of the form $$f_j(z) = \frac{\psi_{j,1}}{t \zeta_j(z)} + \frac{ \psi_{j,2} }{t^2 \zeta_j(z)^2} + \frac{ \psi_{j,3} }{t^3 \zeta_j(z)^3} +\dots, \qquad j=1,2,3,\dots \\$$ where the $\zeta_j(z)$ are explicit scalar functions with $\zeta_j(0) = 0$ and the $\psi_{j,k}$ are symbolic constant matrices.

In the computations I'm doing I need to compute things like $$G_{j,k}(w) =\operatorname{Res}\limits_{z=0} \frac{f_j(z) f_k(z)}{z-w} .$$ up to some order in $t$ as $t\to \infty$.

This is a toy example where the non-commutativity isn't a problem. In the full problem I then need to iterate the computation by computing similar residues on the $G_{j,k}$'s so the commutativity becomes an issue.

Is there a way to hand Mathematica the functions $\zeta_j(z)$ and have it compute $G_{j,k}(w)$ and $G_{k,j}(w)$ correctly, i.e., so that it doesn't think $G_{j,k} - G_{k,j} = 0$? I would like to do this purely symbolically, but as a minimal example one could take $\zeta_1(z) = z$ and $\zeta_2(z) = \sin z$.

My naïve attempt

Here's what I had hoped might work. I included a parameter $t$ in $f$ and an $O(t^{-5})$ term in $G$ to truncate the number of terms in the output.

$Assumptions = { Table[P[i] \[Element] Matrices[{2, 2}], {i, 1, 4}], (\[Zeta] | t) \[Element] Complexes }; \[Zeta][0] = 0; f = Sum[ \[Psi][j]/(t \[Zeta][w])^j, {j, 1, 4}]; G = Residue[f.f/(w - z), {w, 0}] + O[t, Infinity]^5  This fails though; Mathematica returns only the$O(t^{-5})\$ term (the residue is left unevaluated).

If I replace the last line with the scalar product

G = Residue[f*f/(w - z), {w, 0}] + O[t, Infinity]^5


then it outputs, but the order of matrix multiplication is lost.

-