Gaussian integral of $\exp(A + Bx)$ with $A,B$ non-commuting matrices. Let $X_t$ be the Brownian motion on $\mathbb{R}$, $X_0 = 0$. Let $F_t = \exp(A + BX_t)$, where $A$ and $B$ are antihermitian, non-commuting, finite-dimensional matrices (the antihermitian property can be discarded, but I state it in case it might be useful). I am trying to find, hopefully, a closed-form of  $$\mathbb{E}[F_t] =  \frac{1}{\sqrt{2\pi t}}  \int_{\mathbb{R}}\exp(A+Bx)\exp\left(-\frac{x^2}{2t}\right) dx. $$
I tried by writing $\exp(A+Bx)$ as series, then using eq. (9) appearing in https://arxiv.org/pdf/1707.03861.pdf, $$ (A+B)^n = \sum_{k=0}^n \frac{n!}{k!(n-k)!}\{(A + d_B)^k 1 \}B^{n-k},$$ where $d_B(\bullet) = [B, \bullet]$ and $1$ is the identity matrix, and, finally, using  $$ \mathbb{E}[X_t^{2n}] = \frac{(2n)!}{2^nn!}t^n;  \quad n\in \mathbb{N},$$ with no success so far.
 A: You could use the Dyson expansion:
$$e^{A+xB} = e^A \{I+x\int_0^1 dt e^{-yA} B e^{yA}+x^2 \int_0^1 dy_1 \int_{y_1}^1 dy_2 (e^{-y_2 A} B e^{y_2 A}) (e^{-y_1 A}Be^{y_1 A})+\cdots\}  $$
Since $x$ is a scalar random variable it can be pulled out of the matrices, and the expectation can be calculated term by term:
$$e^{A}\{I + <x^2>\int_0^1 dt_1 \int_{y_1}^1 dy_2 (e^{-y_2 A} B e^{y_2 A}) (e^{-y_1 A}Be^{y_1 A})+\cdots\ \}$$ Note that the odd powers of $x$ drop out. The expressions can be written more concisely, e.g. defining $U(y)=e^{-yA}Be^{yA}$, using diagrams to keep track (always $y_n>y_{n-1}$), and the time ordering operator $\mathcal{T}$ to make it formally similar to a commutative expansion, but the bottom line is that you have to be able to integrate over terms like $e^{-yA}Be^{yA}$.
If you are lucky, you may be able to recognize the averaged series as an expansion of another matrix expression and resum. For example, note that the first non trivial term in the expansion is the same as the first term in ${1\over 2} (e^{A+\sqrt<x^2>B}+e^{A-\sqrt<x^2>B})$
