# Singular matrix in derivation of stationary distribution of AR(1) process

I am trying to derive the stationary distribution of the AR(1) process $$\mathbf{x}_t = A \mathbf{x}_{t-1} + \varepsilon_t$$ with $$\varepsilon_t \sim \mathcal{N}(0, Q)$$ and $$A = \begin{pmatrix} 1 & c \\ 0 & 1 \\ \end{pmatrix}$$ where $$0 \le c \le 1$$.

As the noise term is Gaussian, the stationary distribution will be Gaussian as well, i.e. $$\lim_{t\rightarrow \infty} p(X_t) = \mathcal{N}(\mu,\Sigma)$$. In this case $$\mathbb{E}[x_t] = \mu = 0$$. However, I run into a problem when trying to determine the covariance matrix $$\Sigma$$ of the stationary distribution: \begin{align} \Sigma &= A^2 Var[\mathbf{x}_{t-1}] + Var[\varepsilon_t] \\ \Leftrightarrow \Sigma &= (I - A^2)^{-1} Q \end{align} $$I - A^2 = \begin{pmatrix} 0 & 2c \\ 0 & 0 \\ \end{pmatrix}$$ is not invertible. From a simulation I ran, I know that $$\Sigma$$ converges, so I assume the stationary distribution should exist. Is there a way I can derive it analytically?

I realized that the given AR(1) process is not wide-sense stationary as $$|A| = 1$$. As a result the variance of $$\mathbf{x}_t$$ depends on the time lag $$dt$$ and diverges as $$dt$$ goes to infinity.