Background
I am studying the paper Lee et al (2012), located at arXiv:1207.7174. In this paper, we study the spin-boson model under a polaron transformation. The Hamiltonian is $$H = \frac{\epsilon}{2}\sigma_z + \frac{\Delta}{2}\sigma_x + \sum_k\omega_k b_k^{\dagger}b_k + \sigma_z\sum_k \omega_k (b_k^{\dagger}-b_k)$$ and the polaron operator is defined by $$F = \sigma_z\sum_k\frac{g_k}{\omega_k}(b_k^{\dagger}-b_k).$$
Problem
I've derived the transformed Hamiltonian $H' = e^FHe^{-F}$ without difficulty. However, I am having trouble calculating Eq. 6 in the paper, which defined the expectation of the bath operators $$D_{\pm} = \exp\left[\pm\sum_k \frac{g_k}{\omega_k}(b_k^{\dagger}-b_k)\right]$$ as $R = \langle D_{\pm}^2\rangle_{H_B},$ where $H_B = \sum_k \omega_kb_k^{\dagger}b_k$. If we introduce the spectral density $$J(\omega) = \pi\sum_kg_k^2\delta(\omega-\omega_k),$$ the paper claims we can write $R$ as $$R = \exp\left[-2\int_0^{\infty}\frac{d\omega}{\pi}\frac{J(\omega)}{\omega^2}\coth(\beta\omega/2)\right].$$
Attempted Solution
For simplicity, I temporarily write $D_+ = e^C$, where $C = \sum_k\frac{g_k}{\omega_k}(b_k^{\dagger}-b_k)$. Then, I can write $R$ as $$R = \langle D_+^2\rangle_{H_B} = \sum_{n=0}^{\infty}\frac{1}{n!}\langle C^n\rangle_{H_B},$$ which reduced the problem to calculating the moments of C. However, at this point I am unsure how to proceed. Any help given would be immensely appreciated.