Accounting for multiplicative seasonality by including external regressors in ARIMA-GARCH model using R?

Question

At the moment, I'm trying to model a time series in preparation for a multivariate analysis.

The time series comprises tweets per hour for a period of nine weeks. As you can see, and as is to be expected from this kind of data, a daily—seemingly multiplicative—seasonality is evident in the time series (there are less tweets during times where people sleep or typically work):

Consequently, the auto.arima()-function from the forecast::-package yields and ARIMA(2, 0, 0)(2, 1, 0)[24] as the best model for the series:

air_arima <- ts_merged_short %>%
  auto.arima(trace = TRUE, 
             approximation = F, 
             stepwise = F,  
             seasonal = T, 
             lambda = NULL)

Best model: ARIMA(2,0,0)(2,1,0)[24]

Series: . 
ARIMA(2,0,0)(2,1,0)[24] 

Coefficients:
  ar1      ar2     sar1     sar2
1.2015  -0.2855  -0.5578  -0.3049
s.e.  0.0250   0.0259   0.0265   0.0274

sigma^2 estimated as 25621:  log likelihood=-9674.21
AIC=19358.42   AICc=19358.46   BIC=19384.95

The model is, however, not adequate at all:

checkresiduals(air_arima, 24)

Ljung-Box test

data:  Residuals from ARIMA(2,0,0)(2,1,0)[24]
Q* = 118.96, df = 20, p-value = 4.441e-16

Model df: 4.   Total lags used: 24

The residuals aren't white noise. Instead, they display serial correlation as well as clusters of volatility:

My take from this is that ARIMA modeling is inadequate because it is unable to account for the time series' heteroscedasticity. GARCH-modeling on the other hand might be a sensible choice to accurately model the series.

To produce an adequate ARIMA-GARCH-model, I've looked into the rugarch::-package. Unfortunately, that package doesn't allow seasonal ARIMA-GARCH-models to be estimated, which is unfortunate since the data clearly shows multiplicative seasonality. But over at quantitative finance, someone pointed out that seasonality can be accounted for using the argument external.regressors in ugarchspec—however, they didn't provide details as to how to do it.

I'd therefore much appreciate step-by-step explanations or hint to resources and tutorials on how to fit SARIMA-GARCH models using R.

PS: I've tried Box-Cox-transforming the series—to no avail, except improved AIC and BIC.