There are two parts to your question. For the second part, assuming X(t) and Y(t') are independent, requires the cross spectral density to be zero. The cross spectral density is the Fourier transform of the cross correlation function. The cross correlation is the ensemble average of the time-shifted product of X(t) and Y(t'), and if these are independent zero-mean processes than the ensemble average is the product of the two means is zero, thus making the cross spectral density zero.
To find an expression for the cross spectral density start by showing the relationship of $X(t)$ to $S_X(f)$. A WSS random function $X(t)$ can be represented by a Fourier series over the interval $T$ $$X(t) = \sum_{m=-\infty}^{+\infty}A_me^{i2\pi mt/T + i\zeta _m}$$ where $A_m$ is the Fourier coefficient related to the spectrum $S_X(f_m)$, the frequency is given by $f_m = m/T$, and the random behavior of $X(t)$ is carried in the random phase $\zeta_m$. The phase $\zeta_m$ is a uniformly distributed random variable between $0$ and $2\pi$; a random value of $\zeta_m$ is assigned to each term in the Fourier series. An infinite number of ensembles of $X(t)$ can be created by choosing new values of the phases $\zeta_m$ from the distribution. The mean value of $X(t)$ is given by $A_0$ with $\zeta_0 = 0$, and for a zero-mean process $A_0 = 0$. (Note: most of the literature that I have seen erroneously treats $\eta_0$ as a random variable; this is wrong because this makes the mean value a member of a distribution instead of a constant as required by stationarity)
The autocorrelation $R(\tau)$ is given by the temporal average over $t$ or the ensemble averages over $\zeta_m$; here we use the ensemble averages $$R(\tau) = \langle X(t)X^*(t+\tau)\rangle = \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty}A_m A_n e^{i2\pi (m-n)t/T -i2\pi n\tau/T} \langle e^{i(\zeta _m-\zeta_n)}\rangle $$ The ensemble averages are zero everywhere except when $m = n$ so that the autocorrelation reduces to $$R(\tau) = \sum_{m=-\infty}^{+\infty}A_m^2 e^{-i2\pi m\tau/T} $$
An estimate of the $m^{th}$ component of the spectrum of $X(t)$ on the finite interval $T$ is given by taking the Fourier transform of the autocorrelation estimate$$S_m = \int_{-T/2}^{T/2}R(\tau) e^{i2\pi p\tau/T}d\tau = \sum_{m=-\infty}^{+\infty}A_m^2 \int_{-T/2}^{T/2} e^{-i2\pi (m-p)\tau/T}d\tau = A_m^2T $$ where the integral is zero except when $p = m$. If we think of $S_m$ as a sample of $S_X(f_m)$ then the Fourier coefficient $A_m$ can be related to the power spectral density$$A_m = \sqrt\frac{S_m}T = \sqrt\frac{S_X(f_m)}T$$ Substituting this value for $A_m$ into the autocorrelation function yields $$R(\tau) = \sum_{m=-\infty}^{+\infty}\frac{S_m}T e^{-i2\pi m\tau/T} $$ Define $2\pi /T = \Delta\omega$ and after substituting this into the expression for autocorrelation take the limit as $T$ goes to infinity $$R(\tau) = \lim_{T \to \infty}\sum_{m=-\infty}^{+\infty} S_m e^{-i\Delta\omega \tau} \Delta\omega = \frac{1}{2\pi}\int_{-\infty}^{+\infty}S(\omega) e^{-i\omega\tau}d\omega $$ where the Reimann sum has been converted to an integral. Assume all continuity and convergence requirements are met.
All of the foregoing was just to establish that $X(t)$ and $Y(t+\tau)$ can be written in terms of their spectra as $$X(t) = \sum_{m=-\infty}^{+\infty}\sqrt\frac{S_X(f_m)}T e^{i2\pi mt/T + i\zeta _m}$$ $$Y(t+\tau) = \sum_{n=-\infty}^{+\infty}\sqrt\frac{S_Y(f_n)}T e^{i2\pi n(t+\tau)/T + i\mu_n}$$ The cross correlation is the ensemble average of the product of $X(t)$ and the complex conjugate of $Y(t+\tau)$ $$R_{XY}(\tau) = \sum_{m=-\infty}^{+\infty}\sum_{n=-\infty}^{+\infty}\frac{\sqrt{S_X(f_m)S_Y(f_n)}}T e^{i2\pi (m-n)t/T -i2\pi n\tau/T} \langle e^{i(\zeta _m-\mu_n)}\rangle $$
To achieve a non-zero cross spectral density the ensemble average of the cross correlation function cannot be zero; this means that phases between the spectral components must be related. Otherwise, if the phases $\zeta_m$ and $\mu_n$ are independent the ensemble average will always be zero. So if the phases between components with the same frequency are assumed to be related by $\zeta_m - \mu_m = \delta_m$ and $\zeta_m - \mu_n = 0$ the cross correlation can be written as $$R_{XY}(\tau) = \sum_{m=-\infty}^{+\infty}\frac{\sqrt{S_X(f_m)S_Y(f_m)}}T e^{-i2\pi m\tau/T+i\delta_m} $$
Again, define $2\pi /T = \Delta\omega$ and after substituting this into the expression for autocorrelation take the limit as $T$ goes to infinity $$R_{XY}(\tau) = \lim_{T \to \infty}\sum_{m=-\infty}^{+\infty} \sqrt{S_X(f_m)S_Y(f_m)} e^{-i\Delta\omega \tau+i\delta(f_m)} \Delta\omega $$ $$ = \frac{1}{2\pi}\int_{-\infty}^{+\infty} \sqrt{S_X(\omega)S_Y(\omega)} e^{-i\omega\tau+i\delta(\omega)}d\omega $$ where the Reimann sum has been converted to an integral.
The cross power spectral density is just the Fourier transform of the cross correlation and it can be obtained by visual inspection of the above expression $$ S_{XY}(\omega) = \sqrt{S_X(\omega)S_Y(\omega)} e^{i\delta(\omega)} $$
A few observations on the above result. Although the question pertained to the spectra of the two processes, the more crucial parameter is the phase relationship $\delta(\omega)$. Also, consider relaxing the WSS requirement. One common mechanism that produces correlated processes is for one of them to be the time evolved version of the other; this only happens when WSS does not apply.
The solution of problem (taken from the book of Brockwell and Davis) will be given in the following.
The process $ Z_t, X_t$ and $Y_t$ are stationary processes and their spectral density
$ f_Z(\lambda), f_X(\lambda)$ and $f_Y(\lambda)$ are, respectively,
\begin{align*}
f_Z(\lambda)&=\frac{\sigma^2}{2\pi},\\
f_X(\lambda)&=|1-2\mathrm{e}^{-i\lambda}|^2\frac{\sigma^2}{2\pi}
=|1-\tfrac12\mathrm{e}^{-i\lambda}|^2\frac{2\sigma^2}{\pi},\\
f_Y(\lambda)&
=\Big|\sum_{j=0}^{k}a_j\mathrm{e}^{-ij\lambda} \Big|^2 f_X(\lambda)\\
&=\Big|\sum_{j=0}^{k}a_j\mathrm{e}^{-ij\lambda} \Big|^2 |1-\tfrac12\mathrm{e}^{-i\lambda}|^2\frac{2\sigma^2}{\pi} .
\end{align*}
Now if $ a_j=2^{-j}, j\ge1 $, then
\begin{align*}
&\bigg|1- \Big|\sum_{j=0}^{k}a_j\mathrm{e}^{-ij\lambda} \Big|^2 |1-\tfrac12\mathrm{e}^{-i\lambda}|^2\bigg| \\
&\qquad=\bigg|1- \Big|\sum_{j=0}^{k}(\tfrac12\mathrm{e}^{-i\lambda})^j \Big|^2
|1-\tfrac12\mathrm{e}^{-i\lambda}|^2 \bigg| \\
&\qquad =|1- |1-(2\mathrm{e}^{i\lambda})^{-(k+1)}|^2 |\\
&\qquad \le 2^{-2(k+1)}<\frac{\varepsilon}{2\sigma^2},\qquad
\text{if } k\ge \dfrac{-\log(\varepsilon/\sigma^2)}{\log4} .
\end{align*}
and
\begin{gather*}
\sup_{-\pi\le\lambda\le\pi}\Big|f_Y(\lambda)-\frac{2\sigma^2}{\pi}\Big| <\frac{\varepsilon}{2},\\
\Big|\frac{\mathsf{Var[Y_t]}}{2\pi}-\frac{2\sigma^2}{\pi}\Big|
=\Bigg|\frac{1}{2\pi}\int_{-\pi}^{\pi}\Big[f_Y(\lambda)-\frac{2\sigma^2}{\pi}\Big]\, \mathrm{d}\lambda \Bigg|\le \frac{\varepsilon}{2}.
\end{gather*}
Hence,
\begin{equation*}
\sup_{-\pi\le\lambda\le\pi}\Big|f_Y(\lambda)-\frac{\mathsf{Var}[Y_t]}{2\pi}\Big| <\varepsilon.
\end{equation*}
Best Answer
This result proves the conjecture is true. Consider a stationary complex random function $\zeta(t)$. Stationarity requires that the autocorrelation function depend on the difference $t –t’$ only and that the mean value must be a constant (not a member of a distribution function nor depend on $t$). Assume the power spectral density (hereafter called the spectrum) is absolutely continuous everywhere. The Fourier transform of $\zeta(t)$ itself does not exist in general, because stationary random functions are generally neither absolutely integrable nor square integrable, i.e. they do not satisfy Dirichlet’s condition. The correct representation of the random variable $\zeta(t)$ is either as a Fourier series or as a Fourier-Stieltjes integral [see, for example, Akira Ishimaru, Wave Propagation and Scattering in Random Media, (IEEE Press, 1997) Appendix A] Here, the Fourier-Stieltjes integral will be used $$\zeta(t) = \int_{-\infty}^{+\infty} e^{i\omega t}dA(\omega)$$ where $dA(\omega)$ is the random amplitude. Typically, the mean value and the autocorrelation function are found using ensemble averages, but this method provides little insight into the behavior of $dA(\omega)$. Hence, the average value and autocorrelation function for $\zeta(t)$ will be found using temporal averages. The ergodic hypothesis states that for a stationary process the temporal average and the ensemble average are equivalent.
$\mathbf Temporal Averages -$ The temporal average $\overline{\zeta(t)}$ is calculated over a finite interval $T$, and then the limit is taken as $T$ goes to infinity $$\overline {\zeta(t)} = \lim_{T\to \infty} \frac 1T \int_{-T/2}^{+T/2}\zeta(t) dt$$ Substitute the Fourier-Stieltjes integral into the equation and perform the $t$ integration to get $$\overline {\zeta(t)} = \lim_{T\to \infty} \frac 1T \int_{-\infty}^{+\infty} dA(\omega) \int_{-T/2}^{+T/2} e^{i\omega t}dt = \lim_{T\to \infty} \frac 1T \int_{-\infty}^{+\infty} dA(\omega) \frac {2sin(\omega T/2)} {\omega} $$ Use the following expressions to take the limit
$$ \lim_{T\to \infty} \frac {2\sin(\omega T/2)}{\omega} \rightarrow 2\pi\delta(\omega) \quad \mathrm{and} \quad \lim_{T\to \infty} \frac 1T = \lim_{T\to \infty} \frac {\Delta\omega} {2\pi}\rightarrow \frac {d\omega}{2\pi} $$ which gives for the temporal average $$\overline {\zeta(t)} = \int_{-\infty}^{+\infty} dA(\omega) \ \delta(\omega) \ d\omega = dA(0) $$ Finally, for a zero-mean process
$$\overline {\zeta(t)} = \langle \zeta \rangle = dA(0) = 0 $$ Thus, the mean value of $\zeta(t)$ obtained via temporal averaging is $dA(0)$ which is the random amplitude evaluated at the origin.
In the following analysis the random variable $\zeta(t)$ will be decomposed into the sum of zero-mean random variable $\zeta_0(t)$ and the mean $\langle\zeta\rangle$ so that the mean value is everywhere explicit. Now the temporal average for the autocorrelation function of $\overline{\zeta(t)\zeta^*(t’)}$ is given by $$ \overline{\zeta(t)\zeta^*(t+\tau)} = \overline{\zeta_0(t)\zeta_0^*(t+\tau)} + \langle \zeta \rangle^2 = \sigma^2C(\tau) + \langle \zeta \rangle^2 $$ $$ \lim_{T\to \infty} \frac 1T \iint_{-\infty}^{+\infty} e^{-i\omega\tau}dA(\omega)dA(\omega’) \int_{-T/2}^{+T/2} e^{i(\omega-\omega’)t}dt $$ Performing the $t$ integration yields $$ \sigma^2C(\tau) + \langle \zeta \rangle^2 = \iint_{-\infty}^{+\infty} e^{-i\omega\tau}dA(\omega)dA(\omega’) \lim_{T\to \infty} \frac 1T \frac {2sin[(\omega-\omega’) T/2]} {(\omega-\omega’)} $$ The following expressions are used in the limit process
$$ \lim_{T\to \infty} \frac {2\sin[(\omega-\omega’) T/2]}{(\omega-\omega’)} \rightarrow 2\pi\delta(\omega-\omega’) \quad \mathrm{and} \quad \lim_{T\to \infty} \frac 1T = \lim_{T\to \infty} \frac {\Delta\omega} {2\pi}\rightarrow \frac {d\omega}{2\pi} $$ Substituting these expressions into the autocorrelation gives $$ \sigma^2C(\tau) + \langle \zeta \rangle^2 = \iint_{-\infty}^{+\infty} e^{-i\omega\tau}dA(\omega)dA(\omega’) \delta(\omega-\omega’) \ d\omega $$ Finally, performing the $\omega$ integration yields $$ \sigma^2C(\tau) + \langle \zeta \rangle^2 = \int_{-\infty}^{+\infty} e^{-i\omega'\tau}|dA(\omega’)|^2 \qquad \mathrm {Equation \ (1)}$$ The quantity $dA(\omega’)$ can be related to the mean value $\langle\zeta\rangle$ and spectrum $S(\omega)$ by taking the inverse Fourier transform of the autocorrelation function
$$ \int_{-\infty}^{+\infty}\sigma^2C(\tau)e^{i\omega\tau} \ d\tau + \int_{-\infty}^{+\infty} \langle\zeta \rangle^2 e^{i\omega\tau} \ d\tau = \int_{-\infty}^{+\infty}|dA(\omega’)|^2\int_{-\infty}^{+\infty}e^{i(\omega-\omega’)\tau}\ d\tau$$ Performing the $\tau$ integration yields $$ S(\omega)+2\pi \langle \zeta \rangle^2 \delta(\omega)= 2\pi \int_{-\infty}^{+\infty} |dA(\omega’)|^2 \delta(\omega-\omega’) $$ Consider the case of $\omega = 0$ for a zero-mean process, i.e. $\langle\zeta\rangle = 0$; both sides of the expression reduce to $$ S(0)+0 = 2\pi \int_{-\infty}^{+\infty} |dA(\omega’)|^2 \delta(\omega’) $$ The Dirac delta function forces the terms on the right side to be evaluated at $\omega’ = 0$ to produce
$$ S(0) = 2\pi |dA(0)|^2\delta(\omega’) = 2\pi \langle \zeta \rangle ^2 \delta(\omega’) =0 $$ where the previous result $dA(0) = \langle\zeta\rangle$ has been used. Thus, it is shown that the spectrum for a zero-mean process is zero when the frequency is zero, i.e. $S(0) = 0$. It can be seen in Equation (1) for the general case that the required value of $|dA(\omega’)|^2$ must be given by $$|dA(\omega’)|^2 = \frac 1{2\pi}S(\omega’) \ d\omega’+|dA(0)|^2 \delta(\omega’) \ d\omega’ $$ $\mathbf Discussion -$ The vanishing of the spectrum at $\omega = 0$ establishes some new requirements on the autocorrelation function and spectrum. For a zero-mean random variable, when $\tau = 0$ the autocorrelation form of the Wiener-Khinchin theorem reduces to
$$ \sigma^2 = \frac 1{2\pi}\int_{-\infty}^{+\infty}S(\omega)\ d\omega $$ This expression has proven useful and has wide application. Now consider the spectrum form of the Wiener-Khinchin theorem when $\omega = 0$. This expression becomes
$$ S(0) = \int_{-\infty}^{+\infty}\sigma^2C(\tau)\ d\tau =0 $$ Besides the information on the behavior of the spectrum, this gives useful information on the autocorrelation function. For example, the autocovariance function must have both positive and negative excursions, and the area under these excursions must be equal otherwise the integral cannot be zero. This means that the simple Gaussian or exponential forms often used for the autocorrelation function are not valid for a zero mean process – they have no negative excursions anywhere.