Understanding Big O notation in this PDE/Calculus of Variations proof I'm having trouble understanding the use of Big O notation in this proof that I've seen while learning about Calculus of Variations.  Here is the problem (the parts of the proof I am having trouble with are bolded):
The problem is in writing the Euler-Lagrange equation for the minimizer $u$ of the functional $$J(u) = \int \limits_{U} (1 + |\nabla u|^{2})^{\frac{1}{2}} \,dx.$$  We are minimizing over the subset of the domain of $J$ on which the functions equal some fixed function $\sigma(x)$ on the boundary of $U$.
If $\phi$ is any smooth test function with compact support in $U$, then on $\partial U$, $\phi(x) = 0$, so that for each $\epsilon \in \Bbb R$, $u + \epsilon \phi = \sigma(x)$ on $\partial U$.  Then if $u$ is the minimizer of $J$, we have $J(u) \leq J(u + \epsilon \phi)$.
Now, computing $J(u + \epsilon \phi)$ gives:
\begin{split} J(u + \epsilon \phi) &= \int \limits_{U} (1 + |\nabla (u + \epsilon \phi)|^{2})^{\frac{1}{2}} \,dx \\ &= \int \limits_{U} (1 + |\nabla u|^{2} + 2 \epsilon \nabla u \cdot \nabla \phi + \epsilon^{2} | \nabla \phi|^{2})^{\frac{1}{2}} \,dx \\ &= \int \limits_{U} (1 + |\nabla u|^{2})^{\frac{1}{2}}(1 + 2 \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + O(\epsilon^{2}))^{\frac{1}{2}} \,dx  \end{split}
Question 1:  What exactly is going on when the last term in the integrand is suddenly replaced by $O(\epsilon^{2})$ in the last step above??  How does this relate to the definition of Big O notation?
Now, we want to use the fact that $(1 + r)^{\frac{1}{2}} = 1 + \frac{1}{2}r + O(r^{2})$ as $r \to 0$ -- the use of Big O in this step/statement makes sense, as this is just by definition of the Taylor Expansion about $r = 0$; the approximation $1 + \frac{1}{2} r$ holds near $r = 0$ with remainder bounded by some constant times $r^{2}$.  So, we have:
$\int \limits_{U} (1 + |\nabla u|^{2})^{\frac{1}{2}}(1 + \underbrace{2 \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + O(\epsilon^{2})}_{r})^{\frac{1}{2}} \,dx$
$= \int \limits_{U} (1 + |\nabla u|^{2})^{\frac{1}{2}}(1 + \frac{1}{2} \left [ 2 \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + O(\epsilon^{2}) \right ] + O(\left [ 2 \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + O(\epsilon^{2}) \right ]^{2})) \,dx $
and somehow this last line just becomes
$=\int \limits_{U} (1 + |\nabla u|^{2})^{\frac{1}{2}}(1 +  \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + O(\epsilon^{2})) \,dx$
Question 2:  How did we go from the second to last line to the last line?  I have a feeling the Big O notation for $r$ means something different than the Big O notation for $\epsilon^{2}$...
Finally, distributing $(1 + |\nabla u|^{2})^{\frac{1}{2}}$, we find that:
$J(u) \leq J(u + \epsilon \phi) = J(u) + \int \limits_{U} \epsilon \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx + O(\epsilon^{2})$
Question 3:  Why can we pull out the ambiguous thing $O(\epsilon^{2})$ from the integral?
The last line gives that
$0 \leq \epsilon \int \limits_{U} \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx + O(\epsilon^{2})$
from which we deduce $\int \limits_{U} \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx = 0$ and then use integration by parts to get the desired Euler-Lagrange equation.
Question 4:  How can we deduce from the above inequality that $\int \limits_{U} \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx=0$?  What happened to the weird $O(\epsilon^{2})$ thing?  Note that $\epsilon$ can be positive or negative, and also for some reason the last step only considers $\epsilon$ small -- again, why?
 A: $O(x)$ just denotes a function which satisfies $|O(x)| < M|x|$ when $|x| < \delta$ for some constants $M$ and $\delta > 0$. However, unlike other symbols for which each occurence in the same context is expected to represent the same object, each occurrence of the $O$ notation may represent a different function from the others. It is used when we no longer have any interest in the exact form of the function, and merely need that bounding condition for the rest of the calculation.
It does not matter what takes the place of $x$ in this notation. Whether $r$ or $\epsilon^2$ or $1 + 2\epsilon$, or even just $1$. It means the same thing.
Question 1:
Note that the $O(\epsilon^2)$ has nothing to do with the square root. It merely represents a portion of the expression that happens to be inside the square root. The exact full expression inside the square root is $$1 + 2 \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + \epsilon^2\frac{|\nabla \phi|^{2}}{1 + |\nabla u|^{2}}$$
But since $\phi$ is smooth with compact support, $|\nabla \phi|^{2}/(1 + |\nabla u|^{2})$ is bounded by some $M$. So 
$$\left|\epsilon^2\frac{|\nabla \phi|^{2}}{1 + |\nabla u|^{2}}\right| < M\left|\epsilon^2\right|$$
I.e., it is $O(\epsilon^2)$. So the interior expression can be written as $$1 + 2 \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + O(\epsilon^2)$$
Question 2: 
Start with the expression inside the square in the second $O$: $$2 \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + O(\epsilon^{2})$$. Again by the compact support of $\phi$, there is some constant $N$ such that $$\left|2 \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}}\right| < N$$ and by definition there is some $M$ such that $\left|O(\epsilon^2)\right| < M\epsilon^2$ for small enough $\epsilon$. By restricting "small enough", we can obtain $\epsilon^2 < |\epsilon|$, so we have $O(\epsilon^2) < M|\epsilon|$. Therefore $$\left[2 \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + O(\epsilon^{2})\right]^2 < (N + M)^2\left|\epsilon^2\right|$$
Hence $$O(\left [ 2 \epsilon \frac{\nabla u \cdot \nabla \phi}{1 + |\nabla u|^{2}} + O(\epsilon^{2}) \right ]^{2}) = O(\epsilon^2)$$
And it should be easy to see that $O(\epsilon^2) + O(\epsilon^2) = O(\epsilon^2)$, which is how you get to the final form.
Question 3
You are correct that $O$-notation is too dodgy to simply say that $O(\epsilon^2)$ is constant with respect to $x$ and thus can be moved outside the integral. All we know about the $O$ function here is that bound on its behavior for small $\epsilon^2$. We don't know anything in general about it's behavior with respect to $x$. However when I called that $M$ constant in the definition, I meant constant with respect to all involved variables, not just $x$. And that is all we need.
Though you didn't explicitly state it, it is apparent that $U$ is a bounded region. In particular, $\int_U\,dx <\infty$. Therefore $$\left|\int_U O(\epsilon^2)\,dx\right| \le \int_U|O(\epsilon^2)|\,dx \le \int_U M\epsilon^2\, dx = \left(M\int_U\,dx\right)\epsilon^2$$
So $$\int_U O(\epsilon^2)\,dx = O(\epsilon^2)$$
Note that this doesn't work when $U$ has infinite measure.
Question 4
When $\epsilon > 0$, dividing by $\epsilon$ gives
$$0 \leq \int \limits_{U} \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx + O(\epsilon^{2}) / \epsilon < \int \limits_{U} \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx + M\epsilon$$
Let $\epsilon \to 0+$ to get $$0 \leq \int \limits_{U} \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx$$
When $\epsilon < 0$, then dividing by it gives
$$0 \geq \int \limits_{U} \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx + O(\epsilon^{2}) / \epsilon > \int \limits_{U} \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx - M|\epsilon|$$
Let $\epsilon \to 0-$ to get $$0 \geq \int \limits_{U} \frac{\nabla u \cdot \nabla \phi}{(1 + |\nabla u|^{2})^{\frac{1}{2}}} \,dx$$
Hence the result.
