How are multi-term products involving differential operators calculated?
Such as:
$$\bigg( \frac{1}{\epsilon^2} \frac{d^2}{d \xi^2} + \frac{2}{\epsilon} \frac{d^2}{d \xi d x_0} + \frac{d^2}{d x_0^2} \bigg)(y_0 + \epsilon y_1 + ...)$$
as given in:
http://www.naturalspublishing.com/files/published/m3183f1twe241h.pdf p. 1121 (or p. 3 in pdf).
Do they follow the same rules as in abstract algebra?
So the above would be like:
$$(a+b+c)(d+e+...)$$
$$=ad+ae+bd+be+cd+ce+...$$
However, since in the case of divergence:
Then I'm not sure if the algebra way applies here?

Simply distribute the terms, and make sure to note that the terms don’t commute, (ie. $AB \neq BA$).
For example:
$$(\dfrac{d}{dx} + \dfrac{d}{dy})(f(x) + g(y)) = (\dfrac{d}{dx} + \dfrac{d}{dy})f(x)+ (\dfrac{d}{dx} + \dfrac{d}{dy})g(y)$$
Or in the example you give: $$\bigg( \frac{1}{\epsilon^2} \frac{d^2}{d \xi^2} + \frac{2}{\epsilon} \frac{d^2}{d \xi d x_0} + \frac{d^2}{d x_0^2} \bigg)(y_0 + \epsilon y_1 + ...) = \bigg( \frac{1}{\epsilon^2} \frac{d^2}{d \xi^2} + \frac{2}{\epsilon} \frac{d^2}{d \xi d x_0} + \frac{d^2}{d x_0^2} \bigg)y_{0} + \bigg( \frac{1}{\epsilon^2} \frac{d^2}{d \xi^2} + \frac{2}{\epsilon} \frac{d^2}{d \xi d x_0} + \frac{d^2}{d x_0^2} \bigg)\epsilon y_{1} + ...$$
Differential Operators over a ring a generally non-commutative.
What you wrote on the bottom on the question is correct.