Let's consider the Taylor expansion of $\sqrt{1-x}$, which is given by:
$$\sqrt{1-x} = 1-\frac{x}{2} – \frac{x^2}{8} + \mathcal{O}(x^3),$$
where $\mathcal{O(x^3)}$ means including all the higher order contributions from (and including) $x^3$ onwards.
However, let's say that we now have the expansion of $\sqrt{1-\frac{a}{x}}$ which we can obtain from the series above by the substitution $x \to \frac{a}{x}$. Note that technically the function is not analytical since it is not differentiable at $x=0$, but let's ignore this for now, as it is besides the point. We get:
$\sqrt{1-\frac{a}{x}} = 1 – \frac{\frac{a}{x}}{2}-\frac{(\frac{a}{x})^2}{8} +\mathcal{O(\frac{a}{x^3})},$
where we did an expansion in terms of powers of $\frac{a}{x}$.
Is this correct notation, in particular with regards to the $\mathcal{O}$?
The reason why I am asking this is because in the document that can be found here on page 9, they state $H = … + \mathcal{O}(S^0)$, which does not make sense to me. First of all because they are expanding in powers of $\frac{a}{S}$ and not ${S}$ ($x \to S$ in that context) and second of all because they did include the zeroth order term in the Taylor expansion (you can't exclude the zeroth order term in your expansion of course). The notation on page 10 also does not make sense to me.
Best Answer
In general when you want to expand something of the form: $(a+b)^\alpha$ with $a\gg b$, the usual method is to factor by $a$ and then reduce the problem of the Taylor expansion about $1$ (or equivalently, just use directly the Taylor formula with the n-th derivatives): $$ (a+b)^\alpha = a^\alpha(1+b/a)^\alpha = a^\alpha+\alpha a^{\alpha-1}b+O(a^{\alpha-2}b) $$ To use the Taylor expansion at $1$, you need to assume $b/a\to0$ and the subsequent terms are in powers of $b/a$. As always, there is no absolute size, only relative size. What you are perhaps more used to is to have $a$ fixed and let $b\to0$, but the previous method works just as well for $a\to\infty$ and $b$ fixed. It is the latter approach that is used in your reference on the Heisenberg model where they study the semiclassical limit with $S\to\infty$.
Expanding in terms of $a/S$ allows you just as in the previous example to reduce to a Taylor expansion with a small parameter. Since $a$ is fixed, it's still a power expansion in $S$. Secondly, the constant term of the square root Taylor expansion corresponds to the first term which is $O(S^2)$ due to the extra $S^2$ factor. Since you take the first term of the square root's Taylor expansion, you have precision $O(S^{-1/2})$ but due to the overall factors, you do get the constant error $O(S^0)$.
To make things more explicit, you start from equations $(10)-(12)$: $$ \begin{align} S_k^+ &= \sqrt{2S-n_k}a_k & S_k^- & = a_k^\dagger\sqrt{2S-n_k} & S_k^z &= S-n_k\\ &= \sqrt{2S}a_k + O(1/\sqrt S) & &= \sqrt{2S}a_k^\dagger + O(1/\sqrt S) & &= S-a_k^\dagger a_k \end{align} $$ So you get the expansion: $$ \begin{align} S_i\cdot S_j &= \frac12S_i^+S_j^-+\frac12S_i^+S_j^-+S_i^zS_j^z\\ &= \frac12(\sqrt{2S}a_i + O(1/\sqrt S))(\sqrt{2S}a_j^\dagger + O(1/\sqrt S))+\frac12(\sqrt{2S}a_i^\dagger + O(1/\sqrt S))(\sqrt{2S}a_j+ O(1/\sqrt S))+(S-a_i^\dagger a_i)(S-a_j^\dagger a_j) \\ &= S^2+S(a_ia_j^\dagger+a_i^\dagger a_j-a_i^\dagger a_i-a_j^\dagger a_j)+O(S^0) \end{align} $$ Multiplying by $-J$ and summing over all the edges of the graph, you get equation (26).
For the notation on page 10, you'll need to be more specific, but Physics SE might be a better place to ask. Basically, you want to prove that in the semiclassical limit, your Heisenberg model is a free theory and you interpret the independent modes as magnons.
Hope this helps.