Meaning of power series expansion in perturbation theory I encountered a power series expansion of $x$ in $\epsilon$ when solving for the general solution to Mathieu's equation in this paper. 
$x(\xi, \eta) = x_0(\xi, \eta) + \epsilon x_1(\xi, \eta) + \epsilon^2x_2(\xi,\eta) + \cdots$
Is there an intuitive reason as to why we perform such power series expansion (are we assuming some perturbed $\tilde{x}(t)$ around the exact solution $x(t)$ using the perturbation parameter as $\epsilon$)?
If I didn't know what the right-hand side of the above expansion, how could I have figured it out? Is it a standard expansion for such problems?
I don't think I appreciate why we introduce a perturbation factor either. Is it simply because it makes an approximation in first-order convenient?
 A: Mathieu's equations is a particularly tricky example since the value of the parameter $\delta$ is very important. I recommend starting off with some simpler equations requiring the method of multiple scales.

Is there an intuitive reason as to why we perform such power series expansion (are we assuming some perturbed $\tilde x(t)$  around the exact solution $x(t)$  using the perturbation parameter as $\epsilon$)?

We are looking for good approximations to $x(t)$ in the limit $\epsilon\to0$. It turns out that there are a range of related techniques to find asymptotic approximations to $x(t)$. Sometimes these give power series expansions, but often they do not (they often do not converge, for example, and include dependence on $\epsilon$ other than just positive powers). The simplest case is a simple power series, which works (surprisingly) often.

If I didn't know what the right-hand side of the above expansion, how could I have figured it out? Is it a standard expansion for such problems?

Yes it is fairly standard and can be figured out. There are two parts to this expansion, the sequence of gauge functions in the expansion (1, $\epsilon$, $\epsilon^2$, etc.) which is the simplest expansion we use, and also the multiple time-scales, $\xi=t$ and $\eta=\epsilon t$. This is known as the method of multiple scales. For motivation, look at any perturbation theory textbook, there is usually a worked example showing why an expansion like $x_0(t)+\epsilon x_1(t)+\ldots$ is insufficient in some cases and that something like $x_0(\xi,\eta)+\epsilon x_1(\xi,\eta)+\ldots$ is necessary.

I don't think I appreciate why we introduce a perturbation factor either. Is it simply because it makes an approximation in first-order convenient?

We don't usually introduce the perturbation factor, it's part of the problem. For example, an equation describes the motion of some object might have a relatively small amount of drag, and the scaling of the drag becomes the perturbation factor, i.e.
$$ \frac{\mathrm d^2x}{\mathrm dt^2}=-1-\epsilon|x|^2$$ describes an object being thrown upwards which is affected by gravity and a small amount of air resistance (scaled by $\epsilon$). It's naturally a small effect and we expect a small change from the solution of $$ \frac{\mathrm d^2x}{\mathrm dt^2}=-1, $$ at least while $\epsilon$ is small.
