It is really quite difficult. This gets you to the original paper by Cummings in Phys. Rev. Lett. 140 in 1965, page A1051: http://journals.aps.org/pr/abstract/10.1103/PhysRev.140.A1051 . You must notice first that if $P_e(t)$
is given as you wrote then the coefficients $\frac{e^{-\alpha^2}|\alpha|^{2n}}{{n!}}$ multiplying $\cos()^2$ terms are
localized (spike in values) around $n = n_{av} = |\alpha|^2$.
You can get it from the Stirling formula for large $n$ ($n>10$) by approximating
$n! \approx \sqrt{2 \pi n}\left(\frac{n}{e}\right)^n$
and than trying to find the maximum of the function
$\frac{e^{-n_{av}}|n_{av}|^{n}}{{n!}} \approx 1 /\sqrt{2 \pi n}\left(\frac{e}{n}\right)^n e^{-n_{av}}|n_{av}|^{n}$
with respect to $n$ as it was the continuous variable.
Calculating and putting the derivative to zero you get
$e^{-n_{av}} /(2 \sqrt{2 \pi})e^n n_{av}^n n^{-3/2-n}(1 + 2 n_{av} \log n - 2 n \log n) = 0 $
and the approximate solution, neglecting the small $1$ in the bracket, is
$n=n_{av}$.
Than you can Taylor-expand $g(n+1)^{1/2}$ up to the first order around $n_{av}$ i.e.
$g(n+1)^{1/2} \approx g (n_{av}+1)^{1/2} + g/[2 (n_{av}+1)^{1/2}](n-n_{av})$ to care only about the most contributing terms in the sum.
The energies get linear in $n$ now like for the harmonic oscillator.
You get both the revival time and the short time decay using this approximation
after some page of further math (which is not in the paper as it was the Letter).
Revivals are much simpler since now all cosines oscillate with the multiple
of the same frequency $g/[2 (n_{av}+1)^{1/2}]$ and therefore they are seen directly to rephase. The revival time is therefore
$2 \pi (n_{av}+1)^{1/2}/g$ while you get the decay time $1/g$ using small time expansion
of the sum terms for short times to sum them up to the real dropping exponent.
Paradoxically revivals were not noticed or ignored by Cummings in his original paper deriving $1/g$ collapse rate even if the approximate
formula he derived has them in it and they were discovered later by Eberly
in Phys. Rev. lett. 44, 1980, page 1323:
http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.44.1323 .
He derives a complicated compact formula for the arbitrary case of the nonzero detuning that is long time periodic behavior of the full revival after the time
$2 \pi (n_{av})^{1/2}/g \approx 2 \pi (n_{av}+1)^{1/2}/g$ that resembles slightly
the oscillatory motion of the weakly spreading free particle wave packet if one calculated the overlap of it with itself from the time $0$ (the autocorrelation function) and with the revival period while the fast oscillations are some phase.
To obtain the decay time: Now we have:
$P_e(t) \approx \sum_{n=0}^{\infty}{\frac{e^{-n_{av}}n_{av}^{n}}{n!}} \times$
$\cos^2{[g (n_{av}+1)^{1/2} + g/[2 (n_{av}+1)^{1/2}](n-n_{av})]}$
or
$P_e(t) \approx \sum_{n=0}^{\infty}{\frac{e^{-n_{av}}n_{av}^{n}}{n!}} \times$
$[\cos{[2 g (n_{av}+1)^{1/2} t + g/[ (n_{av}+1)^{1/2}](n-n_{av}) t]}/2+1/2]$
Now we use to formula for the $\cos (a+b)=\cos a \cos b - \sin a \sin b$ twice and drop all terms that
have $\sin (...t) \approx 0$ in it since they they are small and terms with $\cos (...t)$ except the most oscillatory with $2g$
not having the running $n$ we put to $1$ for the same reason.
The result is:
$P_e(t) \approx \sum_{n=0}^{\infty}{\frac{e^{-n_{av}}n_{av}^{n}}{n!}} \times$
$[\cos [2 g (n_{av}+1)^{1/2} t] \cos{[g/[ (n_{av}+1)^{1/2}]n t]}/2+1/2]$
Now using the formula for $\cos(x)= (e^{ix}+e^{-ix})/2$
the series can be readily summed up from the series exponent expansion
$e^x = \sum_{n=0}^{\infty} \frac{x^n}{n!}$
(note the
terms containing the exponent of the exponent):
$P_e(t) \approx (1/4)\cos [2 g (n_{av}+1)^{1/2} t] e^{-n_{av}} [e^{n_{av}e^{i g/[ (n_{av}+1)^{1/2}]t}}$
$+e^{n_{av}e^{-i g/[ (n_{av}+1)^{1/2}t]}}]$ +1/2$
The one more step is to expand again for small $t$ (we care only about the real
parts while the imaginary have the small $\sin (...t)$ again):
$e^{-i g/[ (n_{av}+1)^{1/2}t} \approx 1 - t^2 g^2 /2/ (n_{av}+1)$
and so
$e^{i g/[ (n_{av}+1)^{1/2}t} \approx 1 - t^2 g^2 /2/ (n_{av}+1)$
approximating
$n_{av}+1 \approx n_{av}$
we obtain
$P_e(t) \approx (1/2) \cos [2 g (n_{av}+1)^{1/2} t] e^{-g^2 t^2/2}+ 1/2$
so the decay time is
$1/g$
Best Answer
$$ \hat{a} \hat{a}^{\dagger}\left| n \right\rangle = (n+1) \left| n \right\rangle $$ let $$\hat{n}= \hat{a} \hat{a}^{\dagger}$$ so we have; $$B=\cos \left( \lambda t \sqrt{ \hat{a} \hat{a}^{\dagger}}\right)\left| n \right\rangle = \cos \left(\lambda t \hat{n}^{1/2}\right)\left| n \right\rangle =\sum_m (\lambda t)^{2m} \hat{n}^{m}/(2m!)\left| n \right\rangle =\sum_m (\lambda t)^{2m} (n+1)^{m}/(2m!)\\=\cos(\lambda t\sqrt{n+1})$$ you have B^2 in you case so QED.
We really do not need this fuss in general if you have a function of an operator and hit it to arguments eigenvector you just replace the operator in the argument with its eigenvalue. This can be shown with Taylor series for any nice function generally here we did for cos