Time series decomposition results interpretation

I have a long multi-seasonal time series, and the stl() decomposition got me this:

The remainder is definitely not white noise. Then what should be the next step to decide the model?

Try the model with ARMA error term? but it seems that the stl() decomposition of the remainder term still get non-WN remainder, which confused me.