give

$$\begin{aligned}
\bar{\rho}(\Omega)= \int_{-\infty}^{\infty} dp'' d\zeta'' \exp[-i
\Omega \zeta''] f_2(\zeta'',p'')~, \cr && \label{ffft}\end{aligned}$$

where $\Omega = \omega/\omega_l$ is the conjugate variable of $\zeta$, whose meaning is that of normalized frequency.

The integrals in Eq. (\[ffft\]) cannot be easily performed, directly. As customary, one can transform the final variables $(\zeta'', p'')$ back to the initial variables $(\zeta,p)$ and perform the required integrations with respect to the old variables. This allows to use the fact that $f_2(\zeta'',p'') = f_0(\zeta,p)$. Since $d\zeta'' dp'' = d\zeta d p$ one obtains

$$\begin{aligned}
\bar{\rho}(\Omega) && = \int_{-\infty}^{\infty} dp d\zeta \exp[-i
\Omega \zeta''(\zeta,p)] f_0(\zeta,p) \cr && = \frac{N_0}{2\pi
\sigma_\zeta} \int_{-\infty}^{\infty} dp d\zeta \exp[-i \Omega
\zeta''(\zeta,p)] \exp\left[-\frac{p^2}{2}-\frac{\zeta^2}{2
\sigma_\zeta^2}\right]~, \cr && \label{ffft2}\end{aligned}$$

where $\zeta''(\zeta,p)$ can be obtained from Eq. (\[M1\]) and Eq. (\[M2\]), and reads

$$\begin{aligned}
\zeta''(\zeta,p) =&& \zeta + (B_1+B_2) p + A_1
(B_1+B_2)\sin(\zeta+\phi_1)\cr && +A_2 B_2 \sin\left(\zeta + B_1 p +
A_1 B_1 \sin(\zeta+\phi_1)+\phi_2\right) ~.\label{zeta2}\end{aligned}$$

Substituting Eq. (\[zeta2\]) into Eq. (\[ffft2\]) we find, explicitly:

$$\begin{aligned}
\bar{\rho}(\Omega) && = \frac{N_0}{2\pi \sigma_\zeta}
\int_{-\infty}^{\infty} dp d\zeta
\exp\left[-\frac{p^2}{2}-\frac{\zeta^2}{2 \sigma_\zeta^2}\right]~,
\cr && \times \exp\Bigg\{-i \Omega \Bigg[\zeta + (B_1+B_2) p + A_1
(B_1+B_2)\sin(\zeta+\phi_1)\cr && +A_2 B_2 \sin\left(\zeta + B_1 p +
A_1 B_1 \sin(\zeta+\phi_1)+\phi_2\right)\Bigg] \Bigg\}\label{ffft3}\end{aligned}$$

The following step consists in expanding the exponential 