Considering $X$ and $Y$ such that $[X,Y]=\lambda$, which is complex, and $\mu$ is another complex number, prove: $$e^{\mu(X+Y)}=e^{\mu X} e^{\mu Y} e^{-\mu^2\lambda/2}$$ My attempt (so far) is: Expand the exponent. $$\mu(X+Y)=\mu X+ \mu Y$$ and then split it. How can I introduce $\lambda$?
Taylor expansion: $$e^{\mu(X+Y)}=\sum\limits_{n=0}^\infty \frac{(\mu X+\mu Y)^n}{n!}=1+\mu X+\mu Y+\ldots$$
There's a trick to proving this result which you would certainly be forgiven for not spotting!
Consider the quantity $\exp(\mu X) \exp(\mu Y)$ which appears on the right hand side. Now differentiate this with respect to $\mu$. Say what!? Yeah, bear with me, just try it:
$$\frac{d}{d \mu} \exp(\mu X) \exp(\mu Y) = X\exp(\mu X) \exp(\mu Y) + \exp(\mu X) Y \exp(\mu Y)$$
This result is hopefully intuitive --- when dealing with matrices, the product rule still works just fine, but we need to be careful with the ordering of things, since matrices don't commute. Now we want to move that $Y$ in the middle of the second term out to the front, since this will allow us to group the $X$ and $Y$ together into $(X + Y)$, which is getting closer to what we want. To do this, we use the commutation relation $[X,Y] = \lambda$, which implies (as you should check) that $[\mu X, Y] = \mu \lambda$.
Of course, what we want to switch is not $\mu X$ with $Y$, but rather $\exp(\mu X)$ with $Y$. Let me just assert, for now, that the commutation relation we need is
$$ [\exp(\mu X), Y] = \mu\lambda \exp(\mu X) $$
Hence the RHS of my first equation becomes
$$X\exp(\mu X) \exp(\mu Y) + Y \exp(\mu X) \exp(\mu Y) + [\exp(\mu X), Y]\exp(\mu Y) $$ $$= (X + Y + \mu\lambda)\exp(\mu X) \exp(\mu Y) $$
What we have here is a differential equation, which doesn't look all that pleasant... But if I just substitute $\exp(\mu X) \exp(\mu Y)$ by, say, $F$, the structure of the equation becomes clearer:
$$ \frac{d F}{d \mu} = (X + Y + \mu \lambda) F $$
Admittedly, $F$ is a matrix here, but despite that, this sort of differential equation should be familiar to you --- it's a first-order linear ODE! Let's write it as
$$ \frac{d F}{d \mu} + p F = 0 $$
where $p$ is defined appropriately. Now the standard method is to multiply through by an integrating factor
$$\exp\left( \int^\mu p \, d \mu \right) $$
such that the differential equation becomes
$$ \frac{d}{d \mu} \left( F\exp\left( \int^\mu p \, d \mu \right) \right) = 0 $$
which is easily solved:
$$ F\exp\left( \int^\mu p \, d \mu \right) = c = \mathrm{constant\ of\ integration} $$
Actually computing the integral inside the exponential, and putting everything back in terms of $X$ and $Y$, gives us a result which looks very similar to the one we're after (note that I've pulled the exponential over onto the right hand side using the result $\exp(A) \exp(-A) = 1$, where $1$ here is interpreted as the identity matrix or the number $1$ as appropriate):
$$ \exp( \mu X) \exp( \mu Y) = c\exp\left(-\int^\mu p \, d\mu\right) = c \exp\left( \int^\mu (X + Y + \lambda \mu) \, d \mu \right) $$ $$ = c \exp( \mu X + \mu Y + \lambda \mu^2 / 2)$$
All that is left to do now is fix the integration constant, which is clearly seen to be $1$ just by setting $\mu = 0$, and take the term in $\mu^2$ over to the other side. It's legitimate to split the RHS into
$$\exp(\mu X + \mu Y) \exp( \lambda \mu^2 /2)$$
since $\lambda$, $\mu$ and $2$ are just numbers, and commute with everything. Note, however, that the essence of this result is that this same, seemingly standard splitting apart of an exponential doesn't work for matrices! And so we're done. Actually I've cheated a little here because I never proved the result $ [\exp(\mu X), Y] = \mu\lambda \exp(\mu X) $. Let me just give you a hint: first prove that $$[A^n, B] = [A, B]nA^{n-1}$$ for $[A, B] = x$, a number, and use the Taylor expansion of the exponential function.
This is a basic example of a BCH formula. There are many ways to prove it. For example, write the exponential as $$ \exp(\mu X + \mu Y) = \lim_{N\to \infty} \left(1 + \frac {\mu X+ \mu Y}N\right)^N = \dots $$ Because the deviations from $1$ scale like $1/N$, it is equal to $$ = \lim_{N\to \infty} \left[\left(1 + \frac {\mu X}N\right)\left(1 + \frac {\mu Y}N\right)\right]^N $$ Now, we need to move all the $X$ factors to the left and $Y$ factors to the right. Each factor $1+\mu X / N$ commutes with itself, and similarly for $1+\mu Y/N$, of course. However, sometimes the $Y$ factor appears on the left of the $X$ factor and we need to use $$ \left(1 + \frac {\mu Y}N\right)\left(1 + \frac {\mu X}N\right) = \left(1 + \frac {\mu X}N\right) \left(1 + \frac {\mu Y}N\right)\left ( 1 - \frac{\mu^2 (XY-YX)}{N^2} \right ) $$ plus terms of order $O(1/N^3)$ that will disappear in the limit. The only thing we need to count is the number of such permutations of the $(1+\mu Y/N)$ factors with the $(1+\mu X/N)$ factors.
It's not hard. An average $(1+\mu Y/N)$ factor stands on the left side from $N/2\pm O(1)$ of the $(1+\mu x/N)$ factors, and there are $N$ factors of the form $(1+\mu Y/N)$, so we produce $N^2/2\pm O(N)$ factors of the form $$ \left ( 1 - \frac{\mu^2[X,Y]}{N^2} \right) $$ which is a $c$-number that commutes with everything.
Now we just collect the factors inside the limit. On the left, we see $N$ factors $(1+\mu X/N)$ which combine to $\exp(\mu X)$, then we have on the right from them $N$ factors with the $Y$ that combine to $\exp(\mu Y)$, and then there are $N^2/2$ factors of the form $(1-\mu^2[X,Y]/N^2)$ which, in the $N\to \infty $ limit, combine to $\exp(-\mu^2[X,Y]/2)$.
You have obtained several interesting answers. Here is the one I prefer, which actually is a variation of gj255's answer.
Let us attack the problem by computing the function $$f(\mu) := e^{\mu X} Y e^{-\mu X}\:.$$ This function verifies $$f'(\mu) = e^{\mu X} X Y e^{-\mu X} - e^{\mu X} YX e^{-\mu X} = e^{\mu X} [X, Y] e^{-\mu X} = \lambda I\:.$$ Integrating the trivial differential equation we immediately have $$f(\mu) = f(0) + \lambda\mu I\:,$$ that is $$e^{\mu X} Y e^{-\mu X} = Y+\lambda \mu I\:. \tag{1}$$ Next we pass to the core of the problem with a similar reasoning. Define $$g(\mu) := e^{\mu(X+Y)}\:.$$ This function satisfies $$g'(\mu) = (X+Y)g(\mu)\:,\quad g(0)=I\:.\tag{2}$$ Finally define $$h(\mu) := e^{\mu X} e^{\mu Y} e^{-\mu^2\lambda/2}\:.$$ We immediately find that, exploiting (1) in the second identity, $$h'(\mu) = X h(\mu) + e^{\mu X} Y e^{\mu Y} e^{-\mu^2\lambda/2} - \lambda\mu h(\mu)= Xh(\mu) + Yh(\mu) + \lambda\mu h(\mu) - \lambda \mu h(\mu) = (X+Y) h(\mu)\:.$$ Summing up we have that $$h'(\mu) = (X+Y)h(\mu)\:,\quad h(0)=I\:.$$ This is exactly the same Cauchy problem as that in (2). From the uniqueness property of the solution we conclude that: $$h(\mu)= g(\mu)$$ which is our thesis.
If $X$ and $Y$ are matrices, all the presented proof is correct. If they are (unbounded) operators, some subtleties have to be adjusted concerning the topology used in computing the derivatives (this problem is in common with the other presented answers) and referring to the precise statement of uniqueness property of Cauchy problem. Roughly speaking everything goes right if dealing with the weak operator topology.
OP's identity generalizes to the truncated BCH formula
$$\tag{1} e^{\hat{A}}e^{\hat{B}}~=~e^{\hat{A}+\hat{B}+\frac{1}{2}\hat{C}} $$
where the commutator
$$\tag{2} \hat{C}~:=~[\hat{A},\hat{B}]$$
is assumed to commute with both $\hat{A}$ and $\hat{B}$,
$$\tag{3} [\hat{A},\hat{C}]~=~0\quad \text{and}\quad [\hat{B},\hat{C}]~=~0. $$
[In particular, the commutator $\hat{C}$ does not have to be proportional to the unit operator ${\bf 1}$ as OP assumes (v4).]
The identity (1) follows by setting $t=1$ in one of the following three identities
$$\tag{4a} e^{t\hat{A}+\hat{B}}e^{-\hat{B}} ~=~e^{t(\hat{A}-\frac{1}{2}\hat{C})}, $$
$$\tag{4b} e^{-\hat{A}}e^{\hat{A}+t\hat{B}} ~=~e^{t(\hat{B}-\frac{1}{2}\hat{C})}, $$
$$\tag{4c} e^{t(\hat{A}+\hat{B})} ~=~e^{t\hat{A}}e^{t\hat{B}}e^{-\frac{t^2}{2}\hat{C}}. $$
To prove equations (4a-4c), first notice that they are trivially true for $t=0$. Next differentiate their left-hand and right-hand sides wrt. $t$ in order to show that the left-hand and right-hand sides satisfy the same ODE
$$\tag{5a} \hat{f}^{\prime}(t)~=~(\hat{A}-\frac{1}{2}\hat{C})\hat{f}(t), $$
$$\tag{5b} \hat{g}^{\prime}(t)~=~\hat{g}(t)(\hat{B}-\frac{1}{2}\hat{C}), $$
$$\tag{5c} \hat{h}^{\prime}(t)~=~(\hat{A}+\hat{B}) \hat{h}(t), $$
respectively. See also answers by Valter Moretti and gj255. Useful identities are
$$\tag{6} e^\hat{X} \hat{Y} e^{-\hat{X}}~=~e^{[\hat{X},~\cdot~]}\hat{Y}, $$
and
$$ \tag{7} \frac{d}{dt}e^{\hat{X}} ~=~ \int_0^1\!ds~e^{(1-s)\hat{X}}\frac{d\hat{X}}{dt}e^{s\hat{X}}, $$
cf. e.g. my Phys.SE answer here.
Here are some hints.
use the definition of $\exp(x) = 1+x+\frac{1}{2}x^2+...$ on both sides.
Now compare coefficients of $X^mY^n$ on both sides. Start with the $XY$ term as a warm up.
Be careful to introduce commutators if you are swapping $X$s and $Y$s (Lie algebras are non-commutative in general).
You should be able to write down the proof now!
