So let's just look at it from a linear systems perspective:
$$
\begin{array}{ccl}
r(t) &=& Y_1(t) + Y_2(t)\\
&=& Z_1(t) \cos(2\pi f_0 t + \theta) + Z_2(t) \sin(2\pi f_0 t + \theta)\\
&=& [ 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) &=& [ X(t) * H(t) ] [\cos(2\pi f_0 t + \theta) + \sin(2\pi f_0 t + \theta)]\\
&=& \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) &=& \sqrt{2} [ a(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_b(t) &=& \sqrt{2} [ b(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_{a+b}(t) &=& \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) &=& \sqrt{2} [ X(t) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r_\delta(t) &=& \sqrt{2} [ X(t-\delta) * H(t) ] \cos(2\pi f_0 t + \theta - \pi/4)\\
r(t-\delta) &=& \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?