How is modeling the time series error/variance, e.g. ARCH or GARCH models, different from modeling time varying forecast intervals? I'm having a hard time understanding the intuitive difference between modeling the volatility or variance of a time series as it is done in ARCH and GARCH models: 
$$Y_t = c+\epsilon_t+\phi_1Y_{t-1}+...+\phi_p Y_{t-p}+ \epsilon_t$$
and then we model $\epsilon_t$ as: 
$$\epsilon_t=\sigma_tz_t$$ with $z_t$ white noise, and $\sigma_t$ is modeled itself as its own time series process: 
$$\sigma_{t+h}=f(\sigma_t,\sigma_{t-1},\sigma_{t-2},...,\epsilon_t,\epsilon_{t-1},\epsilon_{t-2},...)$$
On one hand. 
On the other hand, using a sophisticated mean time series model where the forecast intervals are estimated "properly" as opposed to using some simplifying assumptions, so that there width is time dependent, e.g. when using FB Prophet and simulating the intervals using MCMC to get intervals that vary with the seasonality. 
Or better still, forecasting full densities, so that the entire shape of the distribution is time dependent, and therefore $\epsilon_t$ and $\sigma_t$ are time dependent as well.  
It seems to me that although the math and the terminology are different, whether we model the variance as its own explicit time series, or whether we just make sure that we model the forecast intervals as best as we can (using MCMC or by forecasting full densities, parametric or non parametric), the end result is the same: A time dependent uncertainty that is conditional on the previous uncertainty estimates. 
What am I missing? Why is volatility modeling considered a separate topic from simulating forecast intervals or modeling full forecast densities?  
 A: GARCH models the entire distribution of a time series with (potentially) time-varying mean, time-varying variance but constant distributional features otherwise (e.g. time-constant higher moments once time variation in the first two moments is accounted for). The distribution of $y_t$ modelled by GARCH is moving up and down due to the conditional mean equation (unless the equation specifies a constant conditional mean) and its spread increases and decreases due to the conditional variance equation; see the direct expressions of the density of $y_t$ under GARCH and ARMA in my answer to the question "What is the difference between GARCH and ARMA?". You can also include seasonal and other effects into the equation of $\sigma^2_t$ in the GARCH model, thereby allowing richer dynamics depending on external regressors.
What is specific to GARCH is that it offers a special kind of time-varying variance, the autoregressive one. Alternatively, if the conditional variance is not autoregressive but time varying in other ways, you can specify models other than GARCH for it. I have forgotten how the FB Prophet works, but perhaps it would be suitable for the task. On another note, you may wish to allow for time-variation in higher-order moments in addition to the first two.
Volatility modeling is considered a special topic as its focus is the time-variation in the second moment, something different than the usual focus on the first moment. If more moments were modeled explicitly in addition to the second moment, we would have more flexible modeling of the density, but when only the first two moments are time-varying with explicit equations for them and the subject-matter focus is more on the variance than the mean, it makes sense to call this volatility modeling.
