How to derive equation (N.15) in Ashcroft and Mermin's Solid State Physics? They state in their book on page 792 the following:

It can be proved, however, that if $A$ and $B$ are operators linear in the $u(R)$ and $P(R)$ of a harmonic crystal, then: $$\langle e^A e^B \rangle = \exp(\frac{1}{2}\langle A^2+2AB+B^2 \rangle)$$

They give as a reference the following paper by Mermin:
http://www.phys.ufl.edu/~maslov/phz6426/Mermin_JMathPhys_7_1038.pdf
But I don't see how to use this paper to prove the identity in the book.
Can anyone help with this?
P.S
$\langle \cdot \rangle$ denotes averaging.
 A: Answering to your question here since it's too long for a comment. I'll change a little the notation by the way. 
First recall that $\langle A \rangle$ is this context is the statistical average of the operator $A$ at thermal equilibrium, given by
\begin{equation} 
\langle A \rangle = \sum_{n=0}^\infty p_n \langle n|A|n\rangle \, ,
\end{equation}
where $p_n=e^{-\beta E_n }/Z$ with $\beta=1/k_\mathrm{B}T$ and $Z$ is the canonical partition function.
In our case (harmonic approximation) the energies are those of the quantum harmonic oscillator:
$$
E_n=\hbar \omega (n+1/2)\, ,
$$
and so it is easy to obtain:
$$
p_n= e^{-\beta \hbar \omega n}(1-e^{-\beta\hbar \omega }) = z^n (1-z) \, ,
$$
where we have defined $z\equiv e^{-\beta \hbar \omega}$.
Therefore the average becomes:
$$
\langle A \rangle = (1-z)\sum_{n=0}^\infty z^n\langle n|A|n\rangle \, .
$$
Now let us define a linear operator in the positions and momenta of the crystal (equivalently, linear in the creation and annihilation operators):
$$
A = c_1a+c_2a^\dagger \, .
$$
Now we can compute $\langle A^2 \rangle = \langle (c_1a+c_2 a^\dagger)^2\rangle$. According to the expression above:
$$
\langle A ^2\rangle = (1-z)\sum_{n=0}^\infty z ^n\langle n | (c_1a+c_2 a^\dagger)^2 | n\rangle =\\ =(1-z) \sum_{n=0}^\infty z ^n\langle n |(c_1^2a^2+c_2^2a^{\dagger 2}+c_1c_2aa^\dagger+c_1 c_2a^\dagger a)  | n\rangle \, .
$$
The non mixed terms (those with $a^2$ and $a^{\dagger 2}$) don't contribute to the sum so we end up with
$$
\langle A^2 \rangle = (1-z)c_1c_2\sum_{n=0}^\infty z ^n\langle n | \underbrace{(aa^\dagger+a^\dagger a)}_{[a,a^\dagger]-2a^\dagger a} | n\rangle=c_1c_2(1-z)\sum_{n=0}^\infty z^n(1+2n)\, .
$$
Thus, after a bit of algebra we obtain:
$$
\langle (c_1a+c_2 a^\dagger)^2\rangle = c_1c_2\left( 1+2\frac{z}{1-z}\right)=2c_1c_2 \left(\frac{1}{e^{\beta \hbar \omega}-1} +\frac{1}{2}\right)\, ,
$$
arriving at the result of the paper I linked in the comment. From there, the general case for several operators should not be too hard to prove as indicated in the paper. 
PS: hopefully there are no errors, it's been a while since I studied QM and I've forgotten many things!
Reference: Solid State Physics, G. Grosso, G. P. Parravicini (2nd. edition) pp. 430-435.
A: OP's formula$^1$ follows from Wick's theorem provided the underlying operators $(\hat{A}_i)_{i\in I}$ satisfy
$$ T(\hat{A}_i\hat{A}_i)~-~:\hat{A}_i\hat{A}_i:~=~\hat{C}_{ij} ~,$$
where the contractions $\hat{C}_{ij}$ are proportional to the identity operator, see e.g. this Phys.SE post and links therein. 
A generalization of OP's formula reads
$$  T( e^{\sum_{i\in I} \hat{A}_i }) ~=~e^{\frac{1}{2}\sum_{j,k\in I} \hat{C}_{jk} } :e^{\sum_{i\in I} \hat{A}_i }:~. $$
--
$^1$Here it is implicitly assumed that $\langle :\cdots:\rangle~=0$, and that we are studying correlators $\langle T(\cdots)\rangle$ with a $T$-ordering in them.
