$E(r_t^2|r_{t-j},j=1,2,...)=E(\sigma^2_{t|t-1}\epsilon_t^2|r_{t-j},j=1,2,...)$. Now given $r_{t-j},j=1,2,...$, you know the values of $\sigma_{t|t-1}$ because $\sigma^2_{t|t-1}=\omega+\alpha r_{t-1}^2$. So you can take it out of expectation to have: $E(r_t^2|r_{t-j},j=1,2,...)=\sigma^2_{t|t-1}E(\epsilon_t^2|r_{t-j},j=1,2,...)$. Also $\epsilon_t$ is independent of $r_{t-j},j=1,2,...$, so (roughly speaking) any function of that is also independent of $r_{t-j},j=1,2,...$. In particular, $f(x)=x^2$. So, $\epsilon_t^2$ is independent of $r_{t-j},j=1,2,...$. Hence, $E(r_t^2|r_{t-j},j=1,2,...)=\sigma^2_{t|t-1}E(\epsilon_t^2)$. Finally note that $\epsilon_t\sim N(0,1)$. So $Var(\epsilon_t)=E(\epsilon_t^2)=1$. Therefore, $E(r_t^2|r_{t-j},j=1,2,...)=\sigma^2_{t|t-1}.$ Now we can change $t$ to $t+h$ to get $h$-step ahead forecasts, i.e. $E(r_{t+h}^2|r_{t-j},j=1,2,...)=\sigma^2_{t+h|t-1}.$ I think you can now answer your 2nd question. Hint: you condition on previous returns $r_t$ that makes your expectation easy to compute.
Why bother with GARCH(1,0)? The $q$ term is easier to estimate than the $p$ term (i.e. you can estimate ARCH($q$) with OLS) anyway.
Nevertheless, my understanding of the way MLE GARCH programs work is they will set the initial GARCH variance equal to either the sample variance or the expected value (that you derive for this case). Without any ARCH terms, the sample variance version would converge to the long-term one (depending on the size of $\delta$). I don't think there would be any change for the expected variance version. So, I'm not sure if you could say it is homoskedastic no matter what (it depends on how you choose the initial variance), but it likely would converge quickly to the expected value for common values of $\delta$.
Best Answer
Here are the equations defining an AR($p$)-ARCH($m$) model for a time series variable $x_t$:
$$ \begin{aligned} x_t &= \varphi_0 + \varphi_1 x_{t-1} + \dotsc + \varphi_p x_{t-p} + \varepsilon_t, \\ \varepsilon_t &= \sigma_t z_t, \\ z_t &\sim i.i.N(0,1), \\ \sigma_t^2 & = \alpha_0 + \alpha_1 \varepsilon_{t-1}^2 + \dotsc + \alpha_m \varepsilon_{t-m}^2. \end{aligned} $$
For the special case of $p=m$ just replace $p$ with $m$ in the top equation.