Basically, in variational methods, we are interested in devising functions $f$ (usually "energy", i.e., squared field) of a given function $u(x)$ (usually a field), in such a way that "small" variations $\epsilon = \delta u$ of $u(x)$ give even smaller variations $\delta f$, i.e., $\delta f \sim {\cal O}(|\delta u|^2)=\epsilon^2$. Effectively, this leads to selection of functions that are stationary around an optimum, "stable" function $u_0$, provided the optimum is a minimum or maximum (but not, for example, a saddle point obtained from second-order derivatives). So this $f[u(x)]$ reaches a maximum or minimum at the targeted $u_0$.
The solution is found by solving Euler's equations, ${\cal L} u_0 = g$, where ${\cal L}$ is a linear operator which depends on the choice of $f$. The easiest way to find such an optimum is by considering quadratic functions of $u$, i.e., energy functions as the square of the field (parabolic functionals), because they reach a single (and hence global) maximum or minimum of $f$ for the optimum $u$.
So, just as with Taylor series expansions of functions, one can consider the expansion of functionals $f$ about the optimum function $u_0(x)$:
$$f(u) = f(u_0) + (u-u_0) \left ( \frac{d f(u)}{d u} \right )_{u=u_0} + \ldots
$$
Limiting to the first two terms when we search for the optimum, the expansion shows that the optimum as the solution of $df/du = 0$ at $u=u_0$ will depend on $u(x)$ (and hence on $x$) and on its first derivative $u^\prime$. This is our search space $f$ for $u_0$.
In principle, you could consider more complicated functions of $u$, that is $f(x,u,u^\prime, u^{\prime\prime}, \ldots)$ depending on higher-order partial derivatives. But this will not guarantee that the optimum is global. It also makes the search space of higher dimensions, which complicates the search. But it may lead to other better solutions.