Exponential operator expansion In my lectures, the professor discussed that for exponential linear operators it is
$$
\exp(\lambda A + \lambda B) \neq \exp(\lambda A)\exp(\lambda B)
$$
for $AB\neq BA$.
Now I know that the exponential operator has a similar definition as the exponential function. So there must be also the notion of a Taylor series for operators.
Let us assume that $\lambda \ll 1$ such that we can neglect $O(\lambda^3)$ but we have to keep to second order in the real number $\lambda$. If the exponent of the exponential operator consists of only one small operator with the form
$$
\exp(A+\lambda B)
$$
can we then expand to second order in $\lambda$ such that
$$
\exp(A+\lambda B)\approx \exp(A) + \exp(A)\lambda B + \frac{1}{2}\exp(A) \lambda^2 B^2 + .....?
$$
I believe that this is not the case because we also have
$$
\exp(A)\exp(\lambda B)\approx \exp(A) (I+\lambda B + \frac{1}{2} \lambda^2 B^2 + .....)
$$
which is the same as above and that contradicts $AB\neq BA$.
One can also apply the definition of the exponential operator directly by
$$
\exp(A+\lambda B) = \sum_{n=0}^\infty \frac{(A+\lambda B)^n}{n!}
$$
but I am not sure how to decompose the "$A$" part (assuming we can neglect $O(\lambda^3)$) and also how the Taylor expansion would work applied to such an operator without bumping into that contradiction.
Thanks in advance!
 A: The standard way to expand $e^{A+\lambda B}$ is the Dyson series https://en.wikipedia.org/wiki/Dyson_series.
I will give the derivation I know, though I am sure there must be a good online reference I just couldn't seem to track it down:
Consider  $e^{A+\lambda B}=e^{N((A+\lambda B)/N)}= \prod_{k=1}^N e^{(A+\lambda B)/N} $
Thus $${d\over d\lambda}e^{A+\lambda B}=\sum_k e^{(k-1)((A+\lambda B)/N)}{d\over d\lambda}[e^{(A+\lambda B)/N}]e^{(N-k)((A+\lambda B)/N)}$$
We are allowed to make $N$ as large as we wish, and in the limit we may expand the derivative to first order in $1\over N$ to get $${d\over d\lambda}[e^{(A+\lambda B)/N}]\rightarrow {d\over d\lambda}[1+{(A+\lambda B)/N}]= B/N$$
returning everything and going to the continuum limit ( sum to integral):
$${d\over d\lambda}e^{A+\lambda B} = \int_0^1 e^{\tau (A+\lambda B)}Be^{(1-\tau) (A+\lambda B)} d\tau $$ So to first order in $\lambda$ :
$e^{A+\lambda B}\sim e^A +\lambda \int_0^1  e^{\tau A}Be^{(1-\tau) A}d\tau+\cdots$
Higher orders are straightforward, but cumbersome.
A: The actual expansion is
$$\eqalign{\exp(A + \lambda B) &= \sum_{j=0}^\infty \frac{(A+\lambda B)^j}{j!}\cr
&= \sum_{j=0}^\infty \sum_{k=0}^j \frac{P_{j,k}(A,B)}{j!} \lambda^k}$$
where $P_{j,k}$ is the sum of the products of $j$ factors (in all orders) of which $k$ are $B$ and the other $j-k$ are $A$.
Thus the coefficient of $\lambda^0$ in the expansion of
$\exp(A+\lambda B)$ is $\exp(A)$, but the coefficient of $\lambda^1$ is
$$ \sum_{j=1}^\infty P_{j,1}(A,B)/j! = \sum_{j=1}^\infty \sum_{i=0}^{j-1} A^i B A^{j-i-1}/j!$$
It certainly isn't anything as simple as $\exp(A) B$.
EDIT: You can write
$$ \sum_{j=1}^\infty \sum_{i=0}^{j-1} \frac{A^i B A^{j-i-1}}{j!} = \sum_{i=0}^\infty A^i B f_i(A)$$
where
$$f_i(x) = \sum_{j=i+1}^\infty \frac{x^{j-i-1}}{j!} = \frac{\exp(x) (i! - \Gamma(i+1,x))}{x^{i+1} i!} $$
and $\Gamma(\cdot, \cdot)$ is the incomplete Gamma function.
EDIT: For the coefficient of $\lambda^2$, you want to consider all products with two $B$'s.  So
$$ \sum_{j=2}^\infty \sum_{i=0}^{j-2} \sum_{k=0}^{j-2-i}
  \frac{A^i B A^k B A^{j-2-i-k}}{j!} $$
