Doing joint estimation is the preferred way. If you do estimation in two stages, a logical inconsistency arises. In the first stage of seasonal ARIMA estimation there is an effective assumption of conditional homoskedasticity. It is contradicted in the second stage when you explicitly model conditional heteroskedasticity using a GARCH model.

Functions ugarchspec and ugarchfit in package rugarch (see here for a vignette) allow specifying and estimating ARMA+GARCH models simultaneously for a variety of GARCH model classes. Unfortunately, seasonal ARMA models do not seem to be implemented there. Perhaps you could try seasonally adjusting your series before fitting an ARMA+GARCH model (although this would be suboptimal if the "true" model is seasonal ARIMA with conditionally heteroskedastic errors).

Doing joint estimation is the preferred way. If you do estimation in two stages, a logical inconsistency arises. In the first stage of seasonal ARIMA estimation there is an effective assumption of conditional homoskedasticity. It is contradicted in the second stage when you explicitly model conditional heteroskedasticity using a GARCH model.

If you have no MA terms in the ARIMA model, you will still get consistent parameter estimates even when neglecting GARCH errors, but these estimates will be inefficient. If you do have MA terms in the ARIMA model, the estimates of the parameters in the ARIMA model will not even be consistent.

Functions ugarchspec and ugarchfit in package rugarch (see here for a vignette) allow specifying and estimating ARMA+GARCH models simultaneously for a variety of GARCH model classes. Unfortunately, seasonal ARMA models do not seem to be implemented there. Perhaps you could try seasonally adjusting your series before fitting an ARMA+GARCH model (although this would be suboptimal if the "true" model is seasonal ARIMA with conditionally heteroskedastic errors).