ACF of differenced MA(p) process

I have an MA(4) process applied to the first order seasonal difference of $$Y_t$$ as follows:

$$(1-B^s) Y_t = (1+\theta_1B+\theta_2B^2+\theta_3B^3+\theta_4B^4) Z_t$$ where $$Z_t \sim N(0,\sigma^2)$$

This is equivalent to

$$Y_t = \frac{(1+\theta_1B+\theta_2B^2+\theta_3B^3+\theta_4B^4)}{(1-B^s)} Z_t$$ where $$Z_t \sim N(0,\sigma^2)$$

I understand how to refactor the denominator into an infinite order polynomial in the case when it has real roots - i.e. it can be decomposed into $$(1-\phi B)$$ terms which can then be inverted using the geometric sum identity $$a/(1-r) = \sum_k ar^k$$.

But in this case, because the denominator polynomial is a seasonal difference operator with S>2, the roots are complex (basically a unit circle on the complex plane).

So in short the question is how does one expand $$(1-B^S)^{-1}$$ into an infinite sequence? Or how do I determine the MA($$\infty$$) coefficients for the model above (numerically is fine, ie. happy to use numpy and truncate at some lag if it cannot be done analytically)

Assuming $$S$$ is an integer, your $$\{ Y_t\}$$ is a seasonal unit root process---it's SARIMA$$(0,0,4)\times(0,1,0)_S$$. The roots of the AR polynomial are $$S$$th-roots of unity. Such processes do not have MA representation---because such an expansion of $$(1−B^S)^{−1}$$ does not exist. Your second equation (heuristically suggestive and actually correct in the stationary case) does not hold in the unit root case.

(The question is really about existence of ACF for integrated, rather than differenced, MA process. If you difference an MA series, the ACF of the resulting over-differenced series can be computed in the standard way.)