I came across the following neat property:
For an operator $\hat{A}$ which is a linear combination of creation and annihilation operators, we have: $$ \langle e^{\hat{A}} \rangle = e^{\langle \hat{A^2} \rangle/2}.$$
Any help in approaching the proof of this identity would be appreciated.
I expect this to be a result of the Wick theorem. If you are considering an equilibrium situation with a quadratic Hamiltonian all odd moments will vanish. Thus you remain with even powers of your operator. If you now count the number of possible contractions you should get the correct result.
P.S.: I just tried and it works indeed. Note the helpful identity
$$ \frac{(2n-1)!!}{(2n)!} = \frac{1}{n! 2^n}$$
(SPOILER ALERT)
What $\left< \cdot \right>$ actually means is
$$ \left< B \right> \equiv \frac{1}{\mathcal{Z}} \text{Tr} \left\lbrace e^{- \beta H } B\right\rbrace$$
for some operator $B$, $\beta = 1/(k_B T)$, and $\mathcal{Z} = \text{Tr}\ e^{- \beta H }$. If your Hamiltonian now takes the form
$$ H = \sum_{nm} a^\dagger_n h_{nm} a_{m},$$
then you can prove the Wick theorem. Let $\alpha_m$ be either a bosonic creation or annihilation operator. I'll spare you the details but you will find that
$$\left< \alpha_1 \alpha_2 \cdots \alpha_{2n} \right> = \sum\limits_{\substack{\text{all possible} \\ \text{pairings}\ \pi}} \left< \alpha_{\pi(1)}\alpha_{\pi(2)} \right>\cdots \left< \alpha_{\pi(2n-1)}\alpha_{\pi(2n)} \right>$$
and $\left< \alpha_1 \alpha_2 \cdots \alpha_{2n+1} \right> = 0 $. This is all you need for the proof:
$$\left< e^{A} \right> \overset{(1)}{=} \sum_{n=0}^\infty \frac{1}{n!} \left< A^n\right> \overset{(2)}{=} \sum_{n=0}^\infty \frac{1}{(2n)!} \left< A^{2n}\right> \overset{(3)}{=} \sum_{n=0}^\infty \frac{(2n-1)!!}{(2n)!} \left< A^2\right>^n \overset{(4)}{=} \sum_{n=0}^\infty \frac{1}{n!2^n} \left< A^2\right>^n \overset{(1)}{=} e^{\left< A^2\right>/2} . $$
(1) Series expansion.
(2) All odd powers vanish.
(3) There are $(2n-1)(2n-3)(2n-5) \cdots (2n-2n) = (2n-1)!!$ ways of forming $n$ pairs of $\left<A^2 \right>$ by the means of the Wick theorem.
(4) The identity from above.
I didn't check if it holds for fermions as well. There might be a similar identity but the Wick decomposition will change.
Hint: OP's formula follows from a Wick-type theorem $$ \tag{1} T(f(\hat{A})) ~=~ \exp\left(\frac{1}{2}\hat{C}\frac{\partial}{\partial\hat{A}}\frac{\partial}{\partial\hat{A}} \right):f(\hat{A}): $$ between time-ordering $T$ and normal ordering $::$. Here $$\tag{2} \hat{C}~=~T(\hat{A}\hat{A})~-~:\hat{A}\hat{A}:$$ is a contraction. See e.g. this Phys.SE post and links therein.
