Sum of Sine and Cosine with Random Phase as LTI System
So let's just look at it from a linear systems perspective:
$$
\begin{array}{ccl}
r(t) &amp;=&amp; Y_1(t) + Y_2(t)\\
&amp;=&amp; Z_1(t) \cos(2\pi f_0 t + \theta) + Z_2(t) \sin(2\pi f_0 t + \theta)\\
&amp;=&amp; [ X_1(t) * H_1(t) ] \cos(2\pi f_0 t + \theta) + [ X_2(t) * H_1(t)] \sin(2\pi f_0 t + \theta)
\end{array} 
$$
where $X_1$ and $X_2$ are the respective (unlabeled) inputs, $H_1$ and $H_2$ are the impulse responses of the filters, and $*$ is convolution.

You appear to be asking whether the system is LTI when $H_1 = H_2 = H$ and $X_1 = X_2 = X$.  In that case, we get
$$
\begin{array}{ccl}
r(t) &amp;=&amp; [ X(t) * H(t) ] [\cos(2\pi f_0 t + \theta) + \sin(2\pi f_0 t + \theta)]\\
&amp;=&amp; \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)
\end{array} 
$$

So it's linear: 
$$
\begin{array}{ccl}
r_a(t) &amp;=&amp; \sqrt{2} [ a(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_b(t) &amp;=&amp; \sqrt{2} [ b(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_{a+b}(t) &amp;=&amp; \sqrt{2} [ (a(t) + b(t)) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4) = r_a(t) + r_b(t)\\
\end{array} 
$$
But it's not time-invariant:
$$
\begin{array}{ccl}
r(t) &amp;=&amp; \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_\delta(t) &amp;=&amp; \sqrt{2} [ X(t-\delta) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r(t-\delta) &amp;=&amp; \sqrt{2} [ X(t-\delta) * H(t) ] \cos(2\pi f_0 (t-\delta) + \theta - \pi/4)\\
\end{array} 
$$
as $r_\delta(t) \not = r(t-\delta)$.  



EDIT

So, what happens if $X_1 = X_2 = X$ and $X\ \tilde{\ }\ N(0,\sigma^2)$ and $\theta\ \tilde{\ }\ U[0,2\pi)$?

I can see why $Y_1(t) = [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta)$ might be non-Gaussian.  

However, I cannot see why $Y_1 + Y_2 = \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)$ will be Gaussian?

Can you please update your question to explain?