Doubt in the derivation of the field Euler-Lagrange equations I'm looking at a derivation of the Euler-Lagrange equations in a field setting, and one step in the proof is continually eluding me. Let $\phi(\vec x,t)$ be a field and $\mathscr L(\phi,\partial_\mu\phi)$ be a Lagrangian density where $\partial_\mu \equiv (c^{-1}\frac{\partial}{\partial t},\vec\nabla)$ is the indexed covariant spacetime differential operator. We define the action of $\mathscr L$ $$\mathcal S[\phi(\vec x,t)]\equiv \int dt\,d^d x\mathscr L.$$
We now pick a field $\phi'(\vec x,t) = \phi(\vec x, t) + \delta\phi(\vec x, t)$ and wish to extremize $\mathcal S$ - in particular, to solve the principle of least/extremal action $$\delta\mathcal S=0+\mathcal O(\delta\phi^2)$$ (here the higher order terms will be implicitly ignored).
$$\begin{align}0 = \delta\mathcal S \equiv \mathcal S[\phi']-\mathcal S[\phi]
\\ = \int dt\,d^d x\left[\mathscr L(\phi+\delta\phi,\partial_\mu\phi+\partial_\mu\delta\phi) - \mathscr L(\phi,\partial_\mu\phi)\right]
\\ = \int dt\, d^d x\left[\frac{\partial\mathscr L}{\partial\phi}\delta\phi+\frac{\partial \mathscr L}{\partial \partial_\mu\phi}\delta\partial_\mu\phi\right]
\\\delta\mathcal S = \int dt\, d^d x\left[\frac{\partial\mathscr L}{\partial\phi}\delta\phi+\frac{\partial \mathscr L}{\partial \partial_\mu\phi}\partial_\mu(\delta\phi)\right]\tag{1}
\end{align}$$
The majority of this is essentially trivial. However, I can't make heads nor tails of the next step - the only justification given for it is "we now integrate the second term by parts several times to split $(1)$ into two integrals". I know that in the regular Euler-Lagrange derivation you integrate by parts once w.r.t time and end up with a vanishing surface term, but this generalization for some reason makes me uneasy.
$$\delta\mathcal S = \int dt\,d^d x\left[\frac{\partial\mathscr L}{\partial\phi} - \partial_\mu\left(\frac{\partial\mathscr L}{\partial\partial_\mu\phi}\right)\right]\delta\phi\color{red}{+\int dt\,d^d x\left[\partial_\mu\left(\frac{\partial\mathscr L}{\partial\partial_\mu\phi}\delta\phi\right)\right]}$$
My question is simply: how did the red term come about, and why is it itself an integral? Did we somehow perform integration by parts a continuum of times?
 A: Let me polish things up a little bit. I assume you are doing something in relativity, so you are on a Lorentzian manifold $(M,g)$ where $g$ is the metric. You have not specified if your field is real valued, but I assume it is, so $\varphi:M\rightarrow \mathbb{R}$. Thus the configuration bundle of the theory is $M\times \mathbb{R}$, as usual in relativity take coordinates $x^{\mu}$ in $M$. The Lagrangian you are considering is first order, (i.e. it depends at most on the first derivative) so the natural geometrical setting to study the variation is the first jet bundle of $M\times \mathbb{R}$, that is $J^1(M\times \mathbb{R})$. If you find this confusing just ignore it and say that you are on a certain manifold whose local coordinated $(x^{\mu},\varphi,\varphi_{\mu})$, where $\varphi_{\mu}(x)=\partial_{\mu}\varphi(x)$ for some scalar field $\varphi$ at each $x \in M$. 
Now a Lagrangian is an object (a "kind" of $m$-form over $M$) locally of the form $\mathscr{L}(x,\varphi,\varphi_{\mu})d\sigma$ where $d\sigma=dx^1 \wedge dx^2 \wedge\ldots\wedge dx^m$ if we assume that dim$(M)=m$. In general we can integrate $m$-form over compact region of manifolds, so the action of $\mathscr{L}$ in $D\subset M$, compact region with regular boundary, along the field $\varphi$ is 
$$
\mathscr{A}_D(\mathscr{L},\varphi) = \int_D \mathscr{L}(x,\varphi(x),\varphi_{\mu}(x))d\sigma
$$
then you define what a deformation is: take a vertical vector field, (secretly $\delta\varphi$), over $M\times \mathbb{R}$ (i.e. some $X=\delta\varphi(x,\varphi)\frac{\partial}{\partial \varphi}$) compactly supported in the region $D$, to X one can always associate its one parameter group of diffeomorphisms $\Phi_s : M \times \mathbb{R}\rightarrow M \times \mathbb{R}$, which in this particular case is $id_M\times \phi_s$, with $\phi_s : \mathbb{R}\rightarrow \mathbb{R}$ having $\left. \frac{d}{ds}\right|_0\phi_s(x,\varphi)=\delta\varphi(x,\varphi)$. Finally the variation of the action along the vector field $X$ is 
$$
\delta\mathscr{A}_D(\mathscr{L},\varphi,X) =\left.\frac{d}{ds}\right|_0\ \int_D \mathscr{L}(x,\phi_s\circ\varphi(x),\partial_{\mu}(\phi_s\circ\varphi(x)))d\sigma
$$
Now remembering that $(x,\varphi,\varphi_{\mu})$ are coordinates we evaluate the above using smoothness for the commutation of derivatives
$$
\delta\mathscr{A}_D(\mathscr{L},\varphi,X) = \int_D \left( \frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi} \left.\frac{d}{ds}\right|_0(\phi_s\circ\varphi(x)) +\frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi_{\mu}} \partial_{\mu}\left.\frac{d}{ds}\right|_0(\phi_s\circ\varphi(x)) \right)d\sigma=\int_D \left( \frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi} \delta\varphi +\frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi_{\mu}} \partial_{\mu}\delta\varphi \right)d\sigma
$$
then you integrate by parts (here there is a subtle point since in general you are using local coordinates, and so you have to check that you can perform globally integration by parts, it turns out however that you can always do it by passing to covariant derifatives and then going back to usual ones without affect the integral apart for a boundary term which in this setup is zero) the second term and get 
$$
\delta\mathscr{A}_D(\mathscr{L},\varphi,X) =\int_D \left( \frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi}  -\partial_{\mu}\frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi_{\mu}}  \right)\delta\varphi d\sigma+\int_D \partial_{\mu}\left(\frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi_{\mu}} \delta\varphi\right)d\sigma
$$
the latter piece of the integral is an exact differential, so you use Stokes theorem and find that 
$$
\delta\mathscr{A}_D(\mathscr{L},\varphi,X) =\int_D \left( \frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi}  -\partial_{\mu}\frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi_{\mu}}  \right)\delta\varphi d\sigma+\int_{\partial D} \left(\frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi_{\mu}} \delta\varphi\right)d\sigma_{\mu}
$$
now recall that we assumed $\delta\varphi$ where compacly supported in $D$ so on its boundary these is zero and the boundary term vanish. 
Finally a field $\varphi $ is extremal for the action if for all compact regions $D$ of $M$ and for all $ X$ vertical vector fields (see above) compactly supported in $D$ we have 
$$
\delta\mathscr{A}_D(\mathscr{L},\varphi,X)=0
$$
that is 
$$
\frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi}  -\partial_{\mu}\frac{\partial\mathscr{L}(x,\varphi(x),\partial_{\mu}\varphi(x))}{\partial \varphi_{\mu}} =0.
$$
