I am considering a system of two interacting scalar fields: $\psi$, and $\phi$. The Lagrangian is given by:
\begin{equation}
\mathcal{L}[\psi]=\frac{1}{2}\partial_\mu\psi\partial^\mu\psi+\frac{1}{2}\partial_\mu\chi\partial^\mu\chi-\frac{{m_\psi}^2}{2}\psi^2-\frac{{m_\chi}^2}{2}\chi^2-\frac{\sigma^2}{4}\chi^2\psi^2
\end{equation}
$\psi$ is initially in a near vacuum state, while $\chi$ is initially supposed to be in a thermal state. For simplifying the calculations, the thermal state is taken as a state of definite occupation numbers, following the Bose-Einstein distribution. The interaction between $\chi$ and $\psi$ can be assumed to be small, and we can treat $\chi$ as an unchanging background for $\psi$. The exact definition of "small" is not given here, as it will be clear later that it is also part of my question.
I want to calculate, perturbatively, the excitation of $\psi$ as function of time due to its interaction with the thermal bath $\chi$. Since I need to know the dynamics of $\psi$ at an arbitrary time, the final state cannot be assumed to be an equilibrium or free state, so the Gell-Mann-Low theorem does not apply. The problem is I know neither a general method to calculate the effect of $\chi$ on $\psi$, nor an approximating method that still captures that effect. I am satisfied to know the least general method that gives a non-zero excitation, preferably using canonical approach. (I have a limited knowledge of QFT of non-equilibrium systems and finite temperature field theory, so I want to avoid them if possible.)
It would be most ideal for me to know what excitation effects we can capture considering different simplifying assumptions (for example, based on different levels of smallness of the coupling between the fields), and (for each simplifying case) the specific methods we can use to account for these effects. I know a complete answer is very unlikely, so a partial one would be very much appreciated still.
One canonical approach I tried is given below. Is there a different approach, more sophisticated perhaps, but still written in canonical/elementary language, that at least gives a non-zero $\psi$ particle production effect?
I understand that canonical approach must eventually fail. But for my purposes, accuracy is not very important. I am ready to decrease the complexity/generality of the problem, by adding more and more assumptions, perhaps until the canonical approach is adequate. Of course the more generality is preserved, the happier I am. So another way to put it is: What is the most complicated case that we can treat using canonical language?
(One approach that may produce result and can be treated using canonical language is of quasiparticle, but I don't know how to implement it. At some point the concept of particle is no longer valid, but if the coupling is weak enough, then $\chi$'s effect on $\psi$ might be sufficiently described by thinking of $\psi$ as a collection of dressed particles - quasiparticles.)
----------------------------------------------------------------------------------------------------------------------------------------
[The failed try] We derive the Euler-Lagrange equation out of the lagrangian above,
\begin{equation}
(\square+\lambda^2|\phi_0|^2)\psi=-\frac{\sigma^2}{2}\chi^2\psi
\end{equation}
Then, take the Fourier transform
\begin{equation}
\frac{d^2\psi_k(t)}{dt^2}+\omega^2(t)\psi_k(t)=J_k(t)
\end{equation}
Where $\psi_k(t)$ and $J_k(t)$ are the Fourier transforms of $\psi(x)$ and $-\frac{\sigma^2}{2}\chi^2\psi$ respectively.
and then I define the following "Hamiltonian"
\begin{equation}
\text{"}H_{\textbf{k}}\text{"}=\frac{1}{2}p_{\textbf{k}}^2+\frac{1}{2}\omega_{\textbf{k}}^2\psi_{\textbf{k}}^2-J_{\textbf{k}}(t)q_{\textbf{k}}
\end{equation}
with $q_{\textbf{k}}=\psi_{\textbf{k}}$, and $p_{\textbf{k}}=\dot{\psi}_{\textbf{k}}$. We assume that we can still define the annihilation operator as
\begin{equation}
a_{\textbf{k}}(t)=\sqrt{\frac{\omega_{\textbf{k}}}{2}}\left[q_{\textbf{k}}(t)+\frac{i}{\omega_{\textbf{k}}}p_{\textbf{k}}(t)\right]
\end{equation}
The evolution of this operator is given by the following Heisenberg equation
\begin{equation}
\frac{da_{\textbf{k}}}{dt}+i\omega_{\textbf{k}}a_{\textbf{k}}=\frac{i}{\sqrt{2\omega_{\textbf{k}}}}J_{\textbf{k}}(t)
\end{equation}
So
\begin{equation}
a_{\textbf{k}} (t)= a_{\textbf{k}}^0 e^{-i\omega_kt}+ \frac{i}{\sqrt{2 \omega_k}} \int^{t}_{0} dt' e^{i \omega_{\textbf{k}}(t'-t)} J_{\textbf{k}}(t')=a_{\textbf{k}}^0 e^{-i\omega_kt}+A_{\textbf{k}}(t)
\end{equation}
where $a_k^0$ is the evolving unperturbed annihilation operator, and we denoted the second term as $A_{\textbf{k}}(t)$.
The plan is take the expectation value of the time varying occupation number operator $N_k=a_k^\dagger a_k$. Since I do not know how to construct a state in the interacting theory, the best I can try is to sandwich the occupation number operator by a state $|s\big>$ obtained by taking the direct product of a free state of $\psi$ with definite occupation numbers and a free thermal state of $\chi$:
\begin{align}
\big<s|a_{\textbf{k}}^{\dagger}(t)a_{\textbf{k}}(t)|s\big>&\approx\big<s|(a_{\textbf{k}}^{0})^{\dagger}a_{\textbf{k}}^0|s\big>+\big<s|(a_{\textbf{k}}^0)^{\dagger}A_{\mathbf{k}}e^{i\omega_kt}+A^{\dagger}_{\textbf{k}}a_{\textbf{k}}^{0}e^{-i\omega_kt}|s\big>\\
&=\alpha_{\textbf{k}}+2\text{Re}\left[\big<s|(a_{\textbf{k}}^0)^{\dagger}A_{\mathbf{k}}e^{i\omega_kt}|s\big>\right]\\
&=\alpha_{\textbf{k}}+2\text{Re}\left[\frac{i}{\sqrt{2 \omega_k}} \int^{t}_{0} dt' e^{i \omega_{\textbf{k}}t'} .\frac{-\sigma^2}{2}\int d^3p\int d^3q\big<s|(a_{\textbf{k}}^0)^{\dagger}\chi_{\textbf{q}}\chi_{\textbf{p}-\mathbf{q}}\psi_{\textbf{k}-\textbf{p}}|s\big> \right]
\end{align}
Using $\psi_{\textbf{k}}(t)=\frac{1}{\sqrt{2}}\left[a_{\textbf{k}}v_{\textbf{k}}^*(t)+a^\dagger_{-\textbf{k}}v_{\textbf{k}}(t)\right]$ and $\chi_{\textbf{k}}(t)=\frac{1}{\sqrt{2}}\left[b_{\textbf{k}}V_{\textbf{k}}^*(t)+b^\dagger_{-\textbf{k}}V_{\textbf{k}}(t)\right]$, we have
\begin{equation}
\big<s|(a_{\textbf{k}}^0)^{\dagger}\chi_{\textbf{q}}\chi_{\textbf{p}-\mathbf{q}}\psi_{\textbf{k}-\textbf{p}}|s\big>e^{i\omega_kt}=\frac{1}{\sqrt{8}}\big<s|(a_{\textbf{k}}^0)^{\dagger}\left[a_{\textbf{k}-\textbf{p}}^0v_{\textbf{k}-\textbf{p}}^*+(a_{\textbf{p}-\textbf{k}}^0)^{\dagger}v_{\textbf{k}-\textbf{p}}\right]\left[b_{\textbf{q}}V_{\textbf{q}}^*+\\
b_{-\textbf{q}}^{\dagger}V_{\textbf{q}}\right]\left[b_{\textbf{p}-\textbf{q}}V_{\textbf{p}-\textbf{q}}^*+b_{\textbf{q}-\textbf{p}}^{\dagger}V_{\textbf{p}-\textbf{q}}\right]|s\big>
\end{equation}
\begin{equation}
=\frac{1}{\sqrt{8}}\big<s|\left[(a_{\textbf{k}}^0)^{\dagger}a_{\textbf{k}-\textbf{p}}^0v_{\textbf{k}-\textbf{p}}^*\right]\left[b_{\textbf{q}}b_{\textbf{q}-\textbf{p}}^{\dagger}V_{\textbf{p}-\textbf{q}}V_{\textbf{q}}^*+b_{-\textbf{q}}^{\dagger}b_{\textbf{p}-\textbf{q}}V_{\textbf{p}-\textbf{q}}^*V_{\textbf{q}}\right]|s\big>
\end{equation}
\begin{equation}
=\frac{1}{\sqrt{8}}\left[\alpha_{k}v_k^*\delta^{(3)}(\textbf{p})\right]\left[(\beta_q+1)V_q^*V_q\delta^{(3)}(\textbf{p})+\beta_qV_qV_q^*\delta^{(3)}(\textbf{p})\right]
\end{equation}
where $\alpha_{\textbf{k}}$ and $\beta_{\textbf{k}}$ are, respectively, the $\psi$ particle number and $\chi$ particle number number of $\textbf{k}$ mode. Therefore,
\begin{align}
\big<s|(a_{\textbf{k}}^0)^{\dagger}A_{\mathbf{k}}|s\big> e^{i \omega_{\textbf{k}} t}&=-\frac{i\sigma^2}{8\sqrt{\omega_k}} \int^{t}_{0} dt' e^{i \omega_{\textbf{k}}t'}\left[\frac{\alpha_k}{\sqrt{\omega_k}}e^{-i\omega_kt'}\right] \int d^3q \left[(2\beta_q+1)\frac{1}{\omega_q^{\chi}}\right]\\
&=imaginary
\end{align}
where we used $v_{\textbf{k}}=\frac{e^{i \omega_{\textbf{k}} t'}}{\sqrt{\omega_{\textbf{k}}}}$. Therefore $\big<N_{\textbf{k}}(t)\big>\approx\alpha_{\textbf{k}}$
We see that this approach produces a trivial result: no $\chi$ particle is produced. This is to be expected, since the $\psi_k$ operators only act on the $\psi$ part of the Hilbert space, and the $\chi_k$ operators only act on the $\chi$ parts. So, there is no mixing among the modes of $\psi$ and $\chi$.