Ah, what the heck.
So I've been mulling this over for quite some time now: let [$]B[$] be the Black-Scholes generator, and [$]A[$] be SV generator terms, such as the vega, volga, vanna differential operators.
From semi-group theory it is justified to write the formal solution of an SV PDE as [$]C(\tau) = e^{\tau (A + B)} C(0)[$]. Kind of justified as [$] C(0) [$] is not smooth, but let's ignore that for the moment.
We know the Baker-Campbell-Hausdorff formula, but there is another (lesser known) formula called the Zassenhaus formula (see for instance https://arxiv.org/pdf/1702.04681.pdf) which expresses [$] e^{\tau(A+B)} [$] as a product of exponentiated differential operators, in particular one could write [$] e^{\tau(A+B)} = \cdots e^{\tau B} [$]. But [$] e^{\tau B} C(0) = C^{BS} (\tau)[$], which means that if the explicit form of " [$] \cdots [$] " is known then *maybe* a general perturbative solution could be found for any SV model.
Now I think the general form of "[$] \cdots [$]" can be found or is actually known (see link above to ArXiv paper), but I am not sure what the "best" way is to cut-off or Taylor expand the "[$] \cdots [$]" part, i.e. to what order in [$] \tau [$] or other parameters.
Thoughts? My thinking at the moment is that this approach could be too complex, not worth it, but maybe somebody else has some better insights.