It is possible to get a general formula for stationary ARMA(p,q) autocovariance function. Suppose $X_t$ is a (zero mean) stationary solution of an ARMA(p,q) equation:
$$\phi(B)X_t=\theta(B)Z_t$$
Multiply this equation by $X_{t-h}$, $h>q$, take expectations and you will get
$$r(h)-\phi_1r(h-1)-...-\phi_pr(h-p)=0$$
This is a recursive equation, which has a general solution. If all the roots $\lambda_i$ of polynomial $\phi(z)=1-\phi_1z-...-\phi_pz^p$ are different,
$$r(h)=\sum_{i=1}^pC_i\lambda_i^{-h}$$
where $C_i$ are constants which can be derived from the initial conditions. Since $|\lambda_i|>1$ to ensure stationarity it is very clear why the autocorrelation function (which is autocovariance function scaled by a constant) is decaying rapidly (if $\lambda_i$ are not close to one).
I've covered the case of unique real roots of the polynomial $\phi(z)$, all other cases are covered in general theory, but formulas are a bit messier. Nevertheless the terms $\lambda^{-h}$ remain.
Answers to question 2 and 3 more or less follow from this formula. For $AR(1)$ process $r(h)=c\phi_1^h$ and when $\phi_1$ is close to one, i.e. close to non-stationarity, you get the behaviour you describle. The same goes for general formula, if the process is nearly unit-root one of the roots $\lambda_i$ is close to 1 and it dominates other terms, producing the slow decay.
Best Answer
I am not sure if it can be used for a time series, but there is a function
is.matrix_ergodic {popdemo}which tests ergodicity of a matrix.I also think that definition of white noise is stronger than that of ergodicity, because a process can be ergodic even if it is not independent, because it refers to asymptotic property. But independence is often formulated as third condition of whiteness, in addition to time independent first and second moments.