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.
