Comment (by Miguel Zilhão):

Thanks for looking into this. Indeed, I can confirm that setting SummationByParts::operator_type = "Minimal Bandwidth" with the 4th order accurate 2nd derivative operator seems to lead to a stable evolution. Should this maybe be the default choice?

With the optimal stencil, the instabilities seem to grow from the boundary. Could this somehow be due to a mismatch between the truncation error in the bulk and boundary regions?

‌

--
Ticket URL: https://bitbucket.org/einsteintoolkit/tickets/issues/2368/instabilities-with-summation-by-parts