This is an answer to the question
Is this process ergodic? How could I demonstrate such thing?
regarding the random process $$\{Y(t) = I(t)\cos(2\pi f_0t)−Q(t)\sin(2\pi f_0t)\}$$
where $\{I(t)\}$ and $\{Q(t)\}$ are uncorrelated zero-mean stationary ergodic processes with identical autocorrelation function $R_I(z) = R_Q(z) = R(z)$.
I will take ergodicity as short-hand for mean ergodicity, that is, the requirement that for almost all sample paths $y(t)$ of the process
$$\lim_{T\to\infty} \frac{1}{2T} \int_{-T}^T y(t) \,\mathrm dt = E[Y(t)]=0.$$
Now, $y(t) = i(t)\cos(2\pi f_0t)−q(t)\sin(2\pi f_0t)$ where the sample paths $i(t)$ and $q(t)$ satisfy the ergodicity requirement, that is,
\begin{align}
\lim_{T\to\infty} \frac{1}{2T} \int_{-T}^T i(t) \,\mathrm dt &= E[I(t)]=0,\\
\lim_{T\to\infty} \frac{1}{2T} \int_{-T}^T q(t) \,\mathrm dt &= E[Q(t)]=0
\end{align}
but
$$\frac{1}{2T} \int_{-T}^T i(t) \,\mathrm dt \approx 0,\quad \frac{1}{2T} \int_{-T}^T q(t) \,\mathrm dt \approx 0$$
does not necessarily imply that
$$\frac{1}{2T} \int_{-T}^T i(t)\cos(2\pi f_0t) \,\mathrm dt \approx 0,\quad \frac{1}{2T} \int_{-T}^T q(t)\sin(2\pi f_0t) \,\mathrm dt \approx 0.$$
But with a further additional assumption about $\{I(t)\}$ and $\{Q(t)\}$ coupled with vigorous hand-waving (on my part; experts in ergodic theory might have solid reasons to see why what I am claiming is in fact completely true (or, God forbid, completely false!)), I can sort of make it work. The assumption is that $\{I(t)\}$ and $\{Q(t)\}$ are low-pass processes, that is, their common power spectral density $S(f)$ has value $0$ for all $f > f_L$ where $f_L \ll f_0$. What this says is that $i(t)$ and $q(t)$ are very slowly varying functions compared to the rapid oscillations of $\cos(2\pi f_0t)$ and $\sin(2\pi f_0t)$ and so over the short time interval $\left[\left(k-\frac{1}{2}\right)\frac{1}{f_0},\left(k+\frac{1}{2}\right)\frac{1}{f_0}\right]$corresponding to one period of the sinusoid of frequency $f_0$, $i(t)$ and $q(t)$ can be assumed to have approximately fixed value so that I can aver with a straight face that
\begin{align}
\int_{\left(k-\frac{1}{2}\right)\frac{1}{f_0}}^{\left(k+\frac{1}{2}\right)\frac{1}{f_0}} i(t)\cos(2\pi f_0t) \,\mathrm dt \approx i\left(\frac{k}{f_0}\right)\int_{\left(k-\frac{1}{2}\right)\frac{1}{f_0}}^{\left(k+\frac{1}{2}\right)\frac{1}{f_0}}\cos(2\pi f_0t) \,\mathrm dt = 0,\\
\int_{\left(k-\frac{1}{2}\right)\frac{1}{f_0}}^{\left(k+\frac{1}{2}\right)\frac{1}{f_0}} q(t)\cos(2\pi f_0t) \,\mathrm dt \approx q\left(\frac{k}{f_0}\right)\int_{\left(k-\frac{1}{2}\right)\frac{1}{f_0}}^{\left(k+\frac{1}{2}\right)\frac{1}{f_0}}\cos(2\pi f_0t) \,\mathrm dt = 0.
\end{align}
Thus, by dividing the long time interval $[-T,T]$ into short segments of length $\frac{1}{f_0}$ and breaking up the integral over $[-T,T]$ into the sum of integrals over these short segments, we can claim that, for large $T$,
$$\frac{1}{2T} \int_{-T}^T i(t)\cos(2\pi f_0t) \,\mathrm dt \approx 0,\quad \frac{1}{2T} \int_{-T}^T q(t)\sin(2\pi f_0t) \,\mathrm dt \approx 0,$$
hold, or at least are close enough for gummint purposes. In short,
$\{Y(t) = I(t)\cos(2\pi f_0t)−Q(t)\sin(2\pi f_0t)\}$ as defined above is an ergodic process if $\{I(t)\}$ and $\{Q(t)\}$ are assumed to be low-pass processes as compared to what we electrical engineers call the carrier frequency $f_0$.
I remark in conclusion that $S_Y(f)$, the power spectral density of $\{Y(t)\}$ is given by
$$S_Y(f) = \left.\left.\frac 12\right(S(f-f_0) + S(f+f_0)\right)$$
