Classical approximation of coupling to a quantum bath This question is a follow up / more focused variant of the question posed here. With some feedback from other users I realized that the question posed was too broad to answer rigorously, and I have thus decided to write a self contained question asking exactly what I am interested in.
Let us take a single two level system of transition frequency $\omega_{0}$. This system is linearly coupled to a large collection of $l$ harmonic oscillators, acting as an environmental bath. The Hamiltonian of the system is given by
\begin{equation}
H = -\frac{\omega_0}{2}\sigma_z + \sum_l(a^\dagger_l a_l + 1/2) + \sum_l \chi_l \sigma_z (a^\dagger_l + a_l)
\end{equation}
Such a bath can be described in terms of its spectral density $J(\omega) = \sum_l \vert\chi_l\vert^2 \delta(\omega-\omega_l)$ where we can assume the spectral density to be a continuous function due to the high number of modes. 
Now, as put forward in "Quantum simulator of an open quantum system using superconducting qubits: exciton transport in photosynthetic complexes" by Mostame et al. one can approach such a system in its classical limit as
\begin{equation}
H = -\frac{(\omega_0 + \delta\omega(t))}{2}\sigma_z 
\end{equation}
where the environmental bath is now fully contained in the time dependent frequency fluctuations $\delta\omega(t)$ that behaves according to certain statistics and has a power spectral density $S_X(f)$, both of which should be able to take on rather arbitrary forms like any classical noise signal. 
What I am interested in is how one derives 
\begin{equation}
\sum_l(a^\dagger_l a_l + 1/2) + \sum_l \chi_l \sigma_z (a^\dagger_l + a_l) \rightarrow -\frac{\delta\omega(t)}{2}\sigma_z 
\end{equation}
A few steps have been made in figuring out this derivation. As detailed in the answer to my previous question, user DanielSank poses that this type of coupling can essentially be seen as 
$$\sum_{l} \chi_l \sigma_z x_l $$
where $x_l$ is the position of phononic mode $l$. One can envision that for a hot bath, this position is essentially jiggling around, resulting in some time dependent term coupled to $\sigma_z$. I agree with this explanation, and I think it is a step in the right direction. Perhaps some high temperature limit comes into play.
As Daniel then further points out, it might be implicit in the equations is that the phononic modes are of much lower frequency than the two level system itself, $\omega_l \ll \omega_0$ because else we would probably not have simple $\sigma_z$ coupling. Moreover he pointed out that the Hamiltonian of the bath itself is most likely absorbed into the time dependent part by going into the interaction picture. His final point is that one might have to make some assumptions about the memory of the bath.
My question is how we can quantify this further, ideally with a derivation from the initial Hamiltonian to the final one with the required assumptions made explicit.
 A: $\def\ii{{\rm i}}
\def\dd{{\rm d}}
\def\ee{{\rm e}}
\def\Tr{{\rm Tr}}
$As far as I know and would expect, the replacement of a quantum heat reservoir with a noisy classical field cannot be rigorously justified in general. However, for the simple problem of pure dephasing posed here, there is indeed a correspondence between the quantum and classical noise models. The quantum model has Hamiltonian $H = H_S + H_B + H_{SB}$, where the system Hamiltonian is $H_S = (\omega_0/2) \sigma^z$, the bath Hamiltonian is $$H_B = \sum_k \omega_k a_k^\dagger a_k,$$ and the interaction is $H_{SB} = \sigma^z X$, where $X$ is the collective bath coordinate
$$ X = \sum_k \chi_k(a_k+a_k^\dagger).$$
This is known as an independent boson model and is exactly solvable (see Emilio Pisanty's elegant general solution for the time-independent problem here), by close analogy with the solvable classical noise model discussed here.
Before presenting the solution, let us motivate the idea of replacing the quantum bath with a noisy classical field, as follows. Moving to an interaction picture generated by $H_S + H_B$, the Hamiltonian transforms to 
$$ H_{SB}(t) = \sigma^z X(t),$$
where (with $\hbar = 1$)
$$X(t) = \sum_k \chi_k ( \ee^{-\ii\omega_k t} a_k + \ee^{\ii\omega_k t} a_k^\dagger ).$$
Now, let's assume that the reservoir is initially in a thermal state $\rho_B = \ee^{-\beta H_B}/\mathcal{Z}$, with $\beta = 1/k_B T$ the inverse temperature and $\mathcal{Z} = \Tr[\ee^{-\beta H_B}]$ the partition function. Then, $X(t)$ has zero mean, $\langle X(t) \rangle = \Tr[X(t) \rho_B] = 0$. However, its fluctuations are non-vanishing, as quantified by the autocorrelation function
$$ \langle X(t) X(t') \rangle = \int_{-\infty}^\infty \dd \omega\; \ee^{-\ii\omega (t-t')} S(\omega).$$
Here the quantity $S(\omega)$ is defined by
$$
S(\omega) = \left \lbrace 
\begin{array}{ll}
J(\omega)[1 + n(\omega)] & (\omega > 0)\\
J(\lvert\omega\rvert) n(\lvert\omega\rvert) & (\omega < 0),
\end{array}\right.
$$
where $n(\omega) =(\ee^{\beta\omega} - 1)^{-1}$ is the Bose-Einstein distribution and the spectral density is 
$$ J(\omega) = \sum_k \lvert\chi_k\rvert^2 \delta(\omega- \omega_k).$$
(Note that really $S(\omega)$, not $J(\omega)$, plays the role analogous to the noise power spectral density, but unfortunately this terminology has been cemented by several decades of repetition and we're pretty much stuck with it.) In addition, the statistics of $X(t)$ are Gaussian, in the sense that all higher-order cumulants with respect to the initial state $\rho_B$ vanish. Thus, if we imagine replacing $X(t)$ in the above Hamiltonian with a classical zero-mean noisy field having Gaussian statistics and power spectral density $S(\omega)$, we might hope to reproduce the dynamics of the full quantum model. However, this idea neglects the fact that the quantum bath is a dynamical system which is perturbed by its interaction with the qubit, and so we cannot expect the noise statistics to remain independent of the qubit's state, nor even to remain Gaussian, as time progresses. Interestingly, for the independent boson model the statistics do in fact remain Gaussian, because the qubit simply acts to displace the equilibrium position of each bosonic mode, which in turn leads the qubit to decohere without affecting its mean energy.
So, on to the solution for the qubit's reduced density matrix $\rho_S(t) = \Tr_B[U(t) \rho(0) U^\dagger(t)]$, where $U(t)$ is the time evolution operator (in the interaction picture), assuming that the initial state is the (tensor) product $\rho(0) = \rho_S(0)\rho_B$. In order to compute the time evolution operator, the key observation is that $[X(t),X(t')]$ is a $c$-number which commutes with everything. Hence, using the Magnus expansion, one finds that
\begin{align} U(t) & = T\exp \left (-\ii \int_0^t\dd s\; H_{SB}(s)  \right) \\ &= \exp\left(-\ii \int_0^t\dd s\; H_{SB}(s) -\frac{1}{2}\int_0^t\dd s\int_0^s\dd s'\;[H_{SB}(s),H_{SB}(s')] \right),
\end{align}
i.e. on the second line we have reexpressed the intractable time-ordered exponential in terms of a normal exponential which can be readily evaluated. The second term in the exponent involving a commutator yields merely a negligible global phase. Dropping this phase factor, we obtain
\begin{align}
U(t) &= \exp\left\{\frac{1}{2}\sigma^z \sum_k \left[ \alpha_k(t) a_k^\dagger - \alpha^*_k(t) a_k \right]\right\},
\end{align}
where 
$$ \alpha_k(t) = \frac{2\chi_k(1 - \ee^{\ii\omega_k t})}{\omega_k}.$$
Thus, $U(t)$ describes a time-dependent displacement of each mode by an amount $\pm\alpha_k(t)/2$, conditioned on the state of the qubit. Explicitly, we can write
$$ U(t) = \lvert 1 \rangle \langle 1 \rvert \prod_k D(\alpha_k/2) +  \lvert 0 \rangle \langle 0 \rvert \prod_k D(-\alpha_k/2),$$
where $D(\alpha_k/2)$ denotes a standard displacement operator acting on the $k^{{\rm th}}$ mode.
Since $[H,\sigma^z] = 0$, the populations in the eigenbasis of $\sigma^z$ are constant, while the coherences decay as $\rho_{01}(t) = \langle 0 \rvert \rho_S(t) \rvert 1\rangle = \ee^{-\Gamma(t)}\rho_{01}(0)$, where we defined the decoherence function via
$$ \ee^{-\Gamma(t)} = \left\langle \prod_k D(\alpha_k)\right\rangle,$$
where the angle brackets denote an average with respect to the initial bath state $\rho_B$. Now, one can recognise the expectation value of the displacement operator as the Wigner characteristic function of the Gaussian bath state $\rho_B$. This characteristic function can be simply looked up, or otherwise derived by noting that $\Gamma(t)$ is essentially the cumulant generating function of the variable $\sum_k(\alpha_k a_k^\dagger - \alpha_k^* a_k)$, and for a Gaussian state only its second cumulant contributes, therefore
\begin{align} \Gamma(t) &= \sum_k \frac{4 \lvert\chi_k\rvert^2}{\omega_k^2} [1-\cos(\omega_k t)]\coth \left(\beta\omega_k/2\right)\\
& = \int_{-\infty}^\infty \dd\omega\; \frac{4[1-\cos(\omega t)]}{\omega^2} S(\omega). \end{align}
Alternatively, we can write these results in terms of a master equation in the Schroedinger picture:
$$ \dot{\rho} = -\frac{\ii}{2}\omega_0 [\sigma^z,\rho] + \frac{1}{2}\gamma(t)\left( \sigma^z \rho \sigma^z - \rho \right). $$
Here, dots denote time derivatives and the dephasing rate is $\gamma(t) = \dot{\Gamma}(t)$, or explicitly
$$ \gamma(t) = \int\dd\omega\; \frac{\sin(\omega t)}{\omega} S(\omega).$$
Comparison of these results with the answer here shows that the evolution is equivalent to that produced by a classical noisy field with stationary Gaussian statistics and power spectral density $S(\omega)$. However, actually generating such a classical noisy signal seems challenging, since $S(\omega)$ is not an even function of $\omega$ as it would be for a real-valued classical signal. This suggests that such a classical noise approximation makes sense only in the  high-temperature limit, where $n(\omega) \gg 1$, since $S(\omega)$ becomes approximately even for frequencies $\lvert \beta \omega\rvert \ll 1$.
