Invariance of action $\Rightarrow$ covariance of field equations? Is this statement true?
I have only seen examples of this, like the invariance of Electromagnetic action under Lorentz transformations.
How could we prove it?
The action is a scalar, $S$, so how I can't even transform it like $U^{-1}SU$...
In this answer we formally show that a (quasi)symmetry of an action implies a corresponding symmetry of its EOM$^{\dagger}$. The answer does not discuss form covariance of EOM. For further relations between symmetries of action, EOM, and solutions of EOM, see e.g. this Phys.SE post.
Let us first recall the definition of a quasi-symmetry of the action
$$ S_V[\phi]~:=~\int_V \! \mathbb{L}, \qquad \mathbb{L}~:=~{\cal L}~d^nx.\tag{1}$$
It means that the action (1) changes by a boundary integral
$$ S_{V^{\prime}}[\phi^{\prime}] +\int_{\partial V^{\prime}} \!d^{n-1}x~(\ldots) ~=~S_V[\phi]+ \int_{\partial V} \!d^{n-1}x~(\ldots) \tag{2}$$
under the transformation. In the following we will assume that the spacetime integration region $V$ is arbitrary.
Theorem. If a local action functional $S_V[\phi]$ has a quasi-symmetry transformation $$ \phi^{\alpha}(x)~~\longrightarrow~~ \phi^{\prime \alpha}(x^{\prime}), \qquad x^{\mu}~~\longrightarrow~~x^{\prime \mu}, \tag{3}$$ then the EOM $$ e_{\alpha}(\phi(x),\partial\phi(x),\ldots ; x)~:=~\frac{\delta S_V[\phi]}{\delta \phi^{\alpha}(x)}~\approx~0\tag{4}$$ must have a symmetry (wrt. the same transformation) $$ e_{\alpha}(\phi^{\prime}(x^{\prime}),\partial^{\prime}\phi^{\prime}(x^{\prime}),\ldots ; x^{\prime})~\approx~e_{\alpha}(\phi(x),\partial\phi(x),\ldots ; x).\tag{5}$$
I) Formal finite proof: This works both for a discrete and a continuous quasi-symmetry.
$$\begin{align} e_{\alpha}&(\phi(x),\partial\phi(x),\ldots ; x)\cr ~:=~&\frac{\delta S_V[\phi]}{\delta \phi^{\alpha}(x)}\cr ~\stackrel{(2)}{=}~&\frac{\delta S_{V^{\prime}}[\phi^{\prime}]}{\delta \phi^{\alpha}(x)}\cr ~\stackrel{{\ddagger}}{\sim}~&\int_{V^{\prime}}\!d^nx^{\prime}~\frac{\delta S_{V^{\prime}}[\phi^{\prime}]}{\delta \phi^{\prime\alpha}(x^{\prime})} \frac{\delta \phi^{\prime\alpha}(x^{\prime})}{\delta \phi^{\alpha}(x)}\cr ~=~&\int_{V^{\prime}}\!d^nx^{\prime}~e_{\alpha}(\phi^{\prime}(x^{\prime}),\partial^{\prime}\phi^{\prime}(x^{\prime}),\ldots ; x^{\prime}) \frac{\delta \phi^{\prime\alpha}(x^{\prime})}{\delta \phi^{\alpha}(x)}.\end{align} \tag{6}$$
II) Formal infinitesimal proof: This only works for a continuous quasi-symmetry. From the infinitesimal transformation (3)
$$\begin{align} \delta \phi^{\alpha}(x)~:=~&\phi^{\prime \alpha}(x^{\prime})-\phi^{\alpha}(x), \cr \delta x^{\mu}~:=~&x^{\prime \mu}-x^{\mu},\end{align} \tag{7}$$
we define a so-called vertical transformation
$$\begin{align} \delta_0 \phi^{\alpha}(x)~:=~&\phi^{\prime \alpha}(x)-\phi^{\alpha}(x)\cr ~=~&\delta \phi^{\alpha}(x)-\delta x^{\mu} ~d_{\mu}\phi^{\alpha}(x),\cr d_{\mu}~:=~&\frac{d}{dx^{\mu}},\end{align} \tag{8} $$
which transforms of the fields $\phi^{\alpha}(x)$ without transforming the spacetime points $x^{\mu}$. The quasi-symmetry implies that the Lagrangian $n$-form $\mathbb{L}$ transforms with a total spacetime derivative
$$\begin{align} \delta \mathbb{L}~=~&d_{\mu} f^{\mu}~d^nx, \cr \delta_0 \mathbb{L}~=~&d_{\mu}(f^{\mu}-{\cal L}~\delta x^{\mu})~d^nx.\end{align} \tag{9} $$
The EOM (4) are typically of second order, so let us assume this for simplicity. (This assumption is not necessary.) Then the infinitesimal transformation of EOM (4) reads
$$\begin{align} \delta e_{\alpha}(x)~=~&\delta_0 e_{\alpha}(x) +\delta x^{\mu} ~\underbrace{d_{\mu} e_{\alpha}(x)}_{\approx 0}\cr ~\approx~&\delta_0 e_{\alpha}(x) \cr ~=~&\frac{\partial e_{\alpha}(x)}{\partial\phi^{\beta}(x)}\delta_0\phi^{\beta}(x) +\sum_{\mu}\frac{\partial e_{\alpha}(x)}{\partial(\partial_{\mu}\phi^{\beta}(x))}d_{\mu}\delta_0\phi^{\beta}(x)\cr &+\sum_{\mu\leq \nu }\frac{\partial e_{\alpha}(x)}{\partial(\partial_{\mu}\partial_{\nu}\phi^{\beta}(x))}d_{\mu}d_{\nu}\delta_0\phi^{\beta}(x) \cr ~\stackrel{{\ddagger}}{\sim}~& \int_V\! d^ny~ \delta_0\phi^{\beta}(y)\frac{\delta e_{\alpha}(x)}{\delta \phi^{\beta}(y)}\cr ~=~&\int_V\! d^ny~ \delta_0\phi^{\beta}(y)\frac{\delta^2 S_V[\phi]}{\delta \phi^{\beta}(y)\delta \phi^{\alpha}(x)} \cr ~=~& \int_V\! d^ny~ \delta_0\phi^{\beta}(y)\frac{\delta^2 S_V[\phi]}{\delta \phi^{\alpha}(x)\delta\phi^{\beta}(y)} \cr ~=~& \frac{\delta}{\delta \phi^{\alpha}(x)} \int_V\! d^ny~ \delta_0\phi^{\beta}(y)\frac{\delta S_V[\phi]}{\delta \phi^{\beta}(y)} \cr &-\int_V\! d^ny~ \frac{\delta(\delta_0\phi^{\beta}(y))}{\delta \phi^{\alpha}(x)} \frac{\delta S[\phi]}{\delta \phi^{\beta}(y)} \cr ~\sim~& \frac{\delta(\delta_0 S_V[\phi]) }{\delta \phi^{\alpha}(x)} -\int_V\! d^ny~ \frac{\delta(\delta_0\phi^{\beta}(y))}{\delta \phi^{\alpha}(x)} e_{\beta}(y) \cr ~\approx~& \frac{\delta(\delta_0 S_V[\phi]) }{\delta \phi^{\alpha}(x)}\cr ~=~&0. \end{align} \tag{10}$$
In the very last step of eq. (10) we used that the infinitesimal variation
$$\begin{align} \delta_0 S_V[\phi]+\int_V\! d^nx~d_{\mu} \left({\cal L}~\delta x^{\mu} \right) ~=~&\delta S_V[\phi]\cr ~=~&\int_{\partial V} \!d^{n-1}x~(\ldots)\end{align} \tag{11}$$
of the action is a boundary integral by assumption (2), so that its functional derivative (10) must vanish (if it exists).
--
$^{\dagger}$ Terminology and Notation: Equations of motion (EOM) means Euler-Lagrange equations (1). The words on-shell and off-shell refer to whether EOM are satisfied or not. The $\approx$ symbol means equality modulo EOM.
$^{\ddagger}$ Warning: This step is not always justified. The $\sim$ symbol indicates that we have formally integrated by parts and ignored boundary contributions. Also we have assumed that the pertinent functional derivative is well-defined and exists. This caveat is the main shortcoming of the formal proof given here. The point is quite serious, e.g. in the case of a global (=$x$-independent) variation, which typically doesn't vanish on the boundary. So boundary contributions could in principle play a role.
However, instead of using functional derivatives and integrations, it is possible to prove eq. (10) $x$-locally
$$\begin{align} \delta_0 e_{\alpha}(x)~=~&\ldots \cr~=~&\underbrace{E_{\alpha(0)} d_{\mu}}_{=0}\left(f^{\mu}(x)-{\cal L}(x)~\delta x^{\mu}\right)\cr &- \sum_{k\geq 0} d^k\left( \underbrace{e_{\beta}(x)}_{\approx 0} \cdot P_{\alpha(k)}\delta_0\phi^{\beta}(x) \right)\cr ~\approx~& 0 \end{align}\tag{12}$$
using only higher partial field derivatives
$$ P_{\alpha(k)} ~:=~\frac{\partial }{\partial \phi^{\alpha(k)}}, \qquad k~\in~\mathbb{N}_0^n,\tag{13}$$
and higher Euler operators
$$ E_{\alpha(k)} ~:=~\sum_{m\geq k} \begin{pmatrix} m \cr k\end{pmatrix} (-d)^m P_{\alpha(m)}, \tag{14} $$
that all refer to the same spacetime point $x$. This $x$-local approach circumvent the problem of unaccounted boundary contributions.
