How to derive equation (N.15) in Ashcroft and Mermin's Solid State Physics?

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.