Exercise on bosonic vacuum 
Consider bosonic canonical transformation, generated by operator $S = e^{\lambda (a^{\dagger})^2}$. Show, that 
  \begin{equation}
b \equiv SaS^{-1} = a - 2\lambda a^{\dagger}. 
\end{equation}
  Calculate the norm of transformed vacuum $S|0>$ and show that the norm is finite only if $\lambda < 1/2$. 

since $a^\dagger$ commutes with itself $S^{-1} =  e^{-\lambda (a^{\dagger})^2}$ and the first question is easily answered. Let's show that $S|0>$ is indeed a vacuum: 
$$
b=SaS^{-1} \Rightarrow bS = Sa \Rightarrow bS|0> = Sa|0>. 
$$
Since $a|0>=0$ we've proved that S|0> is a vacuum state.
But after that I can't see where I'm wrong. Let's call $S|0> = |0_b>$:
$$
<0_b|0_b> = <0|S^{\dagger} S |0> = <0| e^{\lambda a^2} e^{\lambda (a^\dagger)^2} |0>
$$
I am assuming $\lambda$ is a real parameter. Let's rewrite the insides of the average over vacuum:
$$
e^{\lambda a^2} e^{\lambda (a^\dagger)^2} = e^{\lambda (a^\dagger)^2} e^{\lambda a^2} e^{\frac{\lambda^2}{2} [a^2,{a^\dagger}^2]}
$$
Since $<0|e^{\lambda (a^\dagger)^2} = <0|$ we can drop it. If my calculations are accurate (I hope so), the commutator is:
$$
[a^2,{a^\dagger}^2] = 2(2a^{\dagger} a + 1)
$$
and
$$
<0_b|0_b> = <0| e^{\lambda a^2}  e^{\lambda^2 (2a^\dagger a + 1) }|0> = e^{\lambda^2}
$$
Which is true for every $\lambda$, but this should have been an example of transformation that don't save the "finiteness" of vacuum norm, so I must have mistaken somewhere. 
 A: @J.G. solved the problem for you. As a pedagogical lark, and because it is a wonderful thing to know (I levitated, in the mid-70s, when I learned the method in Gilmore's standard text), I'll sketch how your method could have actually come close, if only you had not mangled the CBH formula. 
The degenerate form you wrote down only works (by dropping the plain wrong 1/2 in front of the commutator in the last exponent) if that commutator is central, i.e., commutes with everything else. While this is true for the Heisenberg algebra of oscillators, it  is not for the SU(1,1) of bilinears you have here and recognized in your question--think of it loosely  as a non-hermitian SU(2). The mutual commutators keep ever rotating each other. 
The proper CBH braiding identity to use, instead, is,
$$
G\equiv e^{\lambda a^2} e^{\lambda (a^\dagger)^2} = e^{\frac{\lambda}{1-4\lambda^2} (a^\dagger)^2} e^{ -\ln (1-4\lambda^2)\cdot  ~ (1/2+ a^\dagger a)} e^{\frac{\lambda}{1-4\lambda^2}  a^2} .
$$
(Read on, below, after the separating line, if you were interested in obtaining it.) 
For $\lambda <1/2$, all is fine; you can see that 
$$
\langle 0| G |0\rangle = \frac{1}{\sqrt{1-4\lambda^2}}~,
$$
by inspection, following your argument on trivial action on the vacuum.  
It diverges for $\lambda\to 1/2$. Also note the indeterminate $\infty \cdot 0$s in the exponents of the leftmost and rightmost exponentials you rightly ignored acting on the vacuum. Unless one were a strong geometer, beyond the call of duty here, one would not dare sensibly continue past the singularity we saw. 

So, now, we can indulge in the pedagogy of deriving the CBH formula, cf 
appendix of this answer.
You already noted the SU(1,1) (SU(2)), in your question, namely 
$$
L_-\equiv -a^2/2, \qquad L_+ \equiv a^{\dagger 2 }/2,  \qquad 2L_3\equiv   1/2 +a^\dagger a, 
$$
closing into the familiar Lie algebra
$$
[L_+, L_-]=2 L_3, \qquad [L_3, L_{\pm} ]=\pm L_{\pm}~~.
$$
But you already know the simplest faithful Pauli matrix irrep of this algebra, 
 $$
L_+=  \begin{pmatrix}
     0&1\\
     0&0
   \end{pmatrix} , \qquad     
L_-=  \begin{pmatrix}
     0&0\\
     1&0
   \end{pmatrix} , \qquad    
2L_3=  \begin{pmatrix}
     1&0\\
     0&-1
   \end{pmatrix} .     
$$


*

*The key point:A group element identity (product of exponentials of generators) holds for all representations; conversely, if it holds for a faithful irrep, such as this doublet (the Pauli matrices), it holds in general, for all reps, as the combinatorics of a putative CBH expansion would be identical, and the Lie algebra structure constants are in common. 
(This is nontrivial: it requires Poincaré's exponential theorem to the   effect  the CBH  series in the exponent is fully in the Lie algebra, and hence the representation is immaterial.)


So one need only derive the above group product braiding identity for the doublet irrep, and one is done!!
$$
e^{-2\lambda L_-}~ e^{2\lambda L_+} =   
 \begin{pmatrix}
     1&0\\
     -2\lambda&1
   \end{pmatrix}   \begin{pmatrix}
     1&2\lambda\\
     0&1
   \end{pmatrix}  =   \begin{pmatrix}
     1&2\lambda\\
     -2\lambda&1-4\lambda^2
   \end{pmatrix} \\ =  \begin{pmatrix}
     1&\frac{2\lambda}{1-4\lambda^2}\\
     0&1
   \end{pmatrix}    \begin{pmatrix}
     \frac{1}{1-4\lambda^2}&0\\
     0&1-4\lambda^2
   \end{pmatrix}   \begin{pmatrix}
     1&0\\
     \frac{-2\lambda}{1-4\lambda^2} &1
   \end{pmatrix}\\ = e^{\frac{2\lambda}{1-4\lambda^2}  L_+} ~e^{ -\ln (1-4\lambda^2)\cdot  ~  2L_3} ~e^{\frac{-2\lambda}{1-4\lambda^2}  aL_-} .
$$
Substituting the oscillator bilinear realization above yields the formula to be proven.
A: From $(a^\dagger)^n|0\rangle=\sqrt{n!}|n\rangle$ we have $e^{\lambda (a^\dagger)^2}|0\rangle =\sum_{n=0}^\infty\frac{\lambda^n\sqrt{(2n)!}}{n!}|2n\rangle$, which has norm $\sum_{n=0}^\infty\frac{\lambda^{2n}(2n)!}{n!^2}$. The Stirling approximation gives $\frac{(2n)!}{n!^2}\approx\frac{4^n}{\sqrt{n\pi}}$ for large $n$, so the convergence condition is $4\lambda^{2}< 1$. (The case $\lambda =\frac{1}{2}$ doesn't work because $\sum_n n^{-1/2}$ diverges.) In fact, the sum is $(1-4\lambda^2)^{-1/2}$.
