Is there an explicit or closed form of the following matrix product series
$$ \sum_{n=0}^{\infty} \frac{A^n B^n}{n!} = ?$$
$A$ and $B$ are real valued square matrices that do not commute, i.e., $\left[A,B\right]\neq0$. $A$ is a diagonal matrix. $B$ is a symmetric matrix.
The power series representation $ \sum_{n=0}^{\infty} \frac{(A B)^n}{n!}$ of the matrix exponential $ \exp(A B) $ resembles to some extend the series above.
Note that there is an exact integral representation of the power series expression $\sum_{n=0}^{\infty} \frac{A^n B^n}{n!}$ that is \begin{align} A^{-1}\int_0^\infty dt e^{-t A^{-1}} \sum_{n=0}^\infty \frac{(2 \sqrt{t B})^{2n}}{4^n (n!)^2} = \sum_{n=0}^{\infty} \frac{A^n B^n}{n!} \end{align} Note that $ I_0 (x):=\sum_{n=0}^\infty x^{2n}/(4^n (n!)^2) $ is the zeroth order Bessel function of the first kind.
The underlying matrix ODE is \begin{align} \dot{Y}(t) = A Y(t) B \end{align} with the solution $Y(t) = \sum_{n=0}^{\infty} \frac{A^n B^n t^n}{n!}$.
However, I am looking for an expression that is explicit without any integral representation or ODE solution.
Also, it would be enlightening to consider the special case where A is the identity matrix plus a small and diagonal (perturbation) Matrix $\sigma$, i.e. $A\approx 1+\sigma$.
Any suggestions to solve this issue are welcome.
I don't really have any good ideas, but if it so happened that A and B anticommuted, like the Pauli matrices $\sigma_3$ and $\sigma_1$, then you may note the alternating signs in $$ (AB)^2= -A^2 B^2,\\ (AB)^3= -A^3 B^3,\\ (AB)^4= A^4 B^4, \\ (AB)^5= A^5 B^5,\\ (AB)^6= -A^6 B^6,\\ (AB)^7= -A^7 B^7,... $$ recursively.
You then have
$$ Y=I +ABt-\frac{(AB)^2t^2}{2!} -\frac{(AB)^3t^3}{3!} +\frac{(AB)^4t^4}{4!} +\frac{(AB)^5t^5}{5!}- \frac{(AB)^6t^6}{6!}- \frac{(AB)^7t^7}{7!} +...\\ =\cos (ABt) + \sin (ABt), $$ which you may refashion into exponentials, if you must...
It manifestly solves your differential equation.
Addendum to revised question on diagonal perturbation
For a diagonal perturbation, the perturbation series is straightforward, $$ A^n= I + n\sigma +\frac{n(n-1)}{2!}\sigma^2+... ~~~~\leadsto \\ Y= \sum_{n=0}^\infty (I + n\sigma +\frac{n(n-1)}{2!}\sigma^2+...)B^nt^n/n! \\ = \left (I +\sigma Bt+ \frac{\sigma^2 B^2t^2}{2!} +...\right ) e^Bt. $$ You may observe, however, that to sum the entire perturbation series in σ, you need, of course, to know the answer for A, since it has an identical commutator with B. Shifting A by the identity cannot change the structure of the problem.