Under weak stationarity, it must be that the covariance function $C(Y_t,Y_{t+\tau})$ depends only on $\tau$, the mean function $E(Y_t)$ is constant for all $t$, and these moments exist (are finite!).
Assume exogeneity of $U_t$ with respect to $X_\tau$ for all $\tau$. Then you can write
$C(Y_t,Y_{t+\tau}) = b^2C(X_t,X_{t+\tau}) + C(U_t,U_{t+\tau})$. Assume $X_t$ is a stationary process, you get
$C(Y_t,Y_{t+\tau}) = k_X(\tau) + C(U_t,U_{t+\tau})$. Substitution gives
$C(Y_t,Y_{t+\tau}) = k_X(\tau) + C(U_t,\phi^{\tau}U_t + \sum^{\tau-1}_{l=0}\phi^{l-1}e_{t+\tau-l}) = k_X(\tau) + C(U_t,\phi^{\tau}U_t)$,
since $e_t$ is white noise. The last term is crucial for the weak stationarity of $Y_t$. You have
$C(U_t,\phi^{\tau}U_t) = \phi^{\tau}C(U_t,U_t) = \phi^{\tau}V(U_t)$.
How can you obtain $V(U_t)$? Transform $U_t = \phi L U_{t} +e_t$ into $U_t(1-\phi L) = e_t$. Hence,
$V(U_t) = V(\frac{e_t}{1-\phi L}) = V(\sum_{l=0}^{\infty}\phi^{l}e_{t-l})$, where the second equality holds only if $|\phi|\leq 1$.
Now, from the white noise assumption, you can simplify this last expression into
$\sum_{l=0}^{\infty} \phi^{2l} V(e_{t}) = \frac{1}{1-\phi^2}\sigma_e$. For $\phi=1$, you obtain that $V(U_t)$ is not finite, hence $C(Y_t,Y_{t+\tau})$ is not finite, and the process $Y_t$ is not stationary.
This shows that autocorrelated process may not be stationary.
However, it is also true that uncorrelated processes may not be stationary.
Consider the simple model: $Y_t = t + e_t$.
$E(Y_t) = t$, hence depends on $t$ which violates the stationarity assumption.
A linear process $X_t$ is defined to be causal if $X_t=\psi(B)w_t$ where $w_t$ are white noises and $\sum_{j=1}^{\infty}|\psi(j)|<\infty$.
$X_t$ is defined to be invertible if we can write $w_t=\pi(B) X_t$ where $\pi(B)=\pi_0 + \pi_1 B+\pi_2 B^2 + \cdots$ and $\sum_{j=0}^{\infty}|\pi(j)|<\infty$.
Apparently, an arbitrary $\text{AR}(p)$ model, $$X_t-\phi_1 X_{t-1}-\phi_2 X_{t-2}-\cdots - \phi_p X_{t-p}=w_t$$
automatically satisfies the requirement to be invertible since $\pi(B)=1-\phi_1 B - \cdots - \phi_p B^p$ and $\sum_{j=0}^{\infty}|\pi(j)|=1+\sum_{j=1}^{p}|\phi_j|<\infty$.
You can take a look at this note.
Best Answer
We first start with AR(1) with $|\phi|<1$. We write $$x_t=\phi^{-1}x_{t-1}+w_{t}=\phi(\phi x_{t-2}+w_{t-1})+w_{t}=\dots=\phi^{k}x_{t-k}+\sum_{j=1}^{k-1}\phi^{j}w_{t-j}$$ We assume $sum_t var(x_t)< \infty$, and we can write $$x_t=\sum \phi^j w_{t-j}$$ now$$E(x_t)=\sum \phi^j E(w_{t-j})=0 $$ autocovariance function is $$cov(x_{t+h},x_t)=E((\sum_{j}\phi^j w_{t+h-j})(\sum_k\phi^k w_{t-k}) )=E((w_{t+h}+\dots+\phi^h w_t+\phi^{h+1}w_{t-1}+\dots))\\(w_t+\phi w_{t-1}+\dots)=\sigma^2_w\sum_j \phi^{h+j}\phi^j=\sigma^2_w\phi^h\sum_j\phi^{2j}=\frac{\sigma_w^2\phi^k}{1-\phi^2} $$
The other way: $$x_t=\phi x_{t-1}+w_t\rightarrow(1-\phi B)x_t=w_t\rightarrow1-\phi B=0\rightarrow x=\frac{1}{\phi}\rightarrow |\frac{1}{\phi}|>1\rightarrow|\phi|<1$$ When $|\phi|>1$, the process is called explosive and one can prove that every explosive process has a stochastically equivalent causal process.
To generalize for AR(p):
$$x_t=\phi_1x_{t-1}+\phi_2x_{t-2}+\dots+\phi_px_{t-p}+w_t$$ We rewrite this as an AR(1) model $y_t=Fy_{t-1}+w_t$
[ \begin{bmatrix} x_{t} \\ x_{t-1} \\ \dots\\ x_{t-h+1} \end{bmatrix}=\begin{bmatrix} \phi_{1} & \phi_{2} & \phi_{3} & \dots & \phi_{p} \\ 1 & 0 & 0 & \dots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \dots & 1 \end{bmatrix}\begin{bmatrix} x_{t-1} \\ x_{t-2} \\ \dots \\ x_{t-p} \end{bmatrix}+\begin{bmatrix} w_{t} \\ 0 \\ \dots \\ 0 \end{bmatrix} ]
Lets call thus augmented AR(1) notation. We can write the AR model recursively $$y_t=F^ty_0+w_t+Fw_{t-1}+F^2w_{t-2}+\dots+F^{t-1}w_1+F^tw_0$$. Recall from eigenvalue decomposition: $F=Q\Lambda Q^{-1}$ and $F^j=Q\Lambda^jQ^{-1}$ where $\Lambda^j$ is \begin{bmatrix} \lambda_{1}^j & 0 & 0 & \dots & 0 \\ 0 & \lambda_{2}^j & 0 & \dots & 0 \\ \dots \\ 0 & 0 & 0 & \dots & \lambda_p^j \end{bmatrix} and is stable when $|\lambda_i|<1$ for all i.
The actual eigenvalues of F satisfy $$\lambda^p-\phi_1\lambda^{p-1}-\phi_2\lambda^{p-2}-\dots-\phi_{p-1}\lambda-\phi_p=0$$ So the eigenvalues are the reciprocal of the values that solve the characteristic polynomial of the AR(p), so the roots should be outside the unit circle