I will answer for the 1-D case, or particle mechanics, instead of field mechanics, but the idea is the same. The approach is similar to that of getting the Killing vector field of a metric, and this approach reduces to that when applied to purely kinetic Lagrangians. The objective is to get the so called Rund-Trautman identity
Theory and set up
Let us assume the following: the system is characterised by a configuration manifold $\mathcal Q$ of dimension $n\equiv dim \mathcal Q$ and a Lagrangian $L:T\mathcal Q \times\mathbb R\rightarrow \mathbb R$, possibly dependent on time.
Definition.- We say that a Lagrangian is (quasi-)invariant with respect to the $r-$parameter transformation
$$
\bar{x} = \phi(x,t) = x + \varepsilon^s\xi_s + o(\varepsilon) \simeq x^i + \varepsilon^s\xi_s^i + o(\varepsilon) ,
$$
$$
\bar{t} = \psi(x,t) = t + \varepsilon^s\tau_s + o(\varepsilon) \simeq t + \varepsilon^s\tau_s,
$$
with $s=1,\dots,r$ ; $\simeq$ means that we drop the higher orders on $\varepsilon$ and we use an index expression; and (quasi-) refers to if there is a divergence term $dG$ or not; if and only if
$$
S[x(t)]- S[\bar{x}(\bar t)] = \varepsilon^s G_s(x(t),t)|_a^b + o(\varepsilon).
$$
Lemma.- A Lagrangian L is invariant under a transformation $\iff$ the following $k$ equations hold:
$$
\frac{\partial L}{\partial t} \tau_s + L\frac{d \tau_s}{dt} + \frac{\partial L}{\partial x^i} \xi^i_s +\frac{\partial L}{\partial \dot{x}^i}(\frac{d \xi^i_s}{dt} - \dot x^i \frac{d\tau_s}{dt} ) = \frac{d}{dt}G_s.
$$
Sketch of the Proof
We can write the invariance condition as if we change the integration over $\bar t$ to $t$ in $S[\bar x]$ as
$$
L\left(\bar{x}(\bar t),\frac{d}{d\bar t}\bar{x}(\bar t), \bar{t}\right)\frac{d\bar{t}}{dt} - L(x,\dot x, t) \simeq \varepsilon^s\frac{d}{dt}G_s(x,t).
$$
where $G_s$ is a divergence term. After this differentiate this equation w.r.t. $\varepsilon^s$ at $\varepsilon^s=0$.
This lemma gives us a general relation between the transformation and the Lagrangian. It can be used in different ways:
- To check if a known transformation $(\phi,\psi)$ is a symmetry of a known Lagrangian $L$, and from here derive the Noether conserved quantities.
- If the Lagrangian is unknown, but the transformations are, you get a system of $r$ PDE's on the Lagrangian L, an can be used to impose symmetries,
- Lastly, we can find the symmetry transformations of a given Lagrangian L, considering that the $\dot x^i$ and their powers are independent, so the coefficients in the polynomial $P(x^i)$ are a system of PDE's to obtain $\xi$ and $\tau$.
The third point of view is the one that you want to find the symmetries of a Lagrangian.
Application
In the case of a natural Lagrangian:
$$
L \equiv \frac{1}{2} g_{ij}(x)\dot{x}^i\dot{x}^j - V(x)
$$
we compute the derivatives
$$
\partial_t L =0, \;\;\;\; \partial_{k}L = \frac{1}{2}\partial_{k}g_{ij}\dot{x}^i\dot{x}^j - \partial_{k}V, \;\;\;\; \partial_{,k}L = \frac{1}{2}g_{il}\left(\dot{x}^i\delta^l_k+\dot{x}^l\delta^i_k\right)
$$
Then the equations become:
$$
0\cdot \tau_s +
\left(
\frac{1}{2} g_{ij}\dot{x}^i\dot{x}^j - V
\right)
\frac{d \tau_s}{dt} +
\left(
\frac{1}{2}\partial_{k}g_{ij}\dot{x}^i\dot{x}^j - \partial_{k}V
\right) \xi^k_s +
\frac{1}{2}g_{il}\left(\dot{x}^i\delta^l_k+\dot{x}^l\delta^i_k\right)
\left(
\frac{d \xi^k_s}{dt} - \dot x^k \frac{d\tau_s}{dt}
\right)
= \frac{d}{dt}G_s.
$$
taking the powers of $\dot x$ as independent and demanding the equations to be satisfied always, we get the first order PDE's:
$$
1 \; : \;\; V\partial_t \tau_s -\partial_k V\xi^k_s = \partial_tG_s
$$
$$
x^i \; : \;\;
- V
\partial_i\tau_s +
\frac{1}{2}\left(g_{il}\partial_t\xi^l_s+g_{li}\partial_t\xi^l_s\right)
= \partial_iG_s
$$
$$
\dot x^i \dot x^j \; : \;\;
-\frac{1}{2} g_{ij}
\partial_t\tau_s +
\frac{1}{2}\partial_{k}g_{ij} \xi^k_s +
\frac{1}{2}\left(
g_{il}\partial_j\xi^l_s +
g_{lj}\partial_i\xi^l_s \right)
= 0
$$
$$
\dot{x}^i\dot{x}^j\dot x^k\; : \;\;
g_{ij}\partial_k\tau_s -
2g_{ik}\partial_j\tau_s= 0
$$
These are the Rund-Trautman identities, or generalised Killing equations to obtain the symmetries of the Lagrangian L. We note that for every $s=1\dots r$ has the same system of equations, and that the number of equations is highly dependant on the form of the Lagrangian.
In the next part I will illustrate one example that particularly interests me. More examples can be found in the references, particularly in [1,4].
Example
A particle in the Hyperbolic Poincaré Disc $\mathbb D$, that has a metric $\frac{2|dz|^2}{1-|z|^2}$, has a Lagrangian
$$
L = \frac{\dot x^2 + \dot y^2}{(1-(x^2+y^2))^2} = \gamma^2( \dot x^2 + \dot y^2)
$$
with $\gamma \equiv \frac{1}{1-(x^2+y^2)}$. Then $g_{ij}=\gamma^2\delta_{ij}$ and $V=0$. I only want time independent transformations, this means that I fix $\tau_s=0$ and $\partial_t=0$, then the R-T identities become
$$
1 \; : \;\; 0 - 0 = \partial_tG_s
$$
$$
x^i \; : \;\;
0 = \partial_iG_s
$$
$$
\dot x^i \dot x^j \; : \;\;
\frac{1}{2}\partial_{k}g_{ij} \xi^k_s +
\frac{1}{2}\left(
g_{il}\partial_j\xi^l_s +
g_{lj}\partial_i\xi^l_s \right)
= 0
$$
$$
\dot{x}^i\dot{x}^j\dot x^k\; : \;\; 0 = 0
$$
we only have one set of equations, the coefficients of the square terms. Taking into account that $\partial_kg_{ij} = 2\gamma (-\gamma^2)(-2x^k)= 4\gamma^3 x^k\delta_{ij} $ we have the system
$$
4\gamma x_k\delta_{ij} \xi^k_s +
\partial_j\xi^i_s +
\partial_i\xi^j_s
= 0
$$
And in components $(x,y)$ we get 3 equations:
$$
\partial_y\xi^y_s
= -2 \gamma \vec{x}\vec{\xi}_s
$$
$$
\partial_x\xi^x_s
= -2 \gamma \vec{x}\vec{\xi}_s
$$
$$
\partial_x\xi^y_s +
\partial_y\xi^x_s
= 0
$$
The second equation tells us directly that one family of solutions is given by the vector field $\vec{\xi}_s=s(-y,x)$, and it checks with the others. So rotations, that could be seen directly from the Lagrangian. If we sum the first two, we get
$$
\nabla\cdot \vec{\xi}_s = -4 \gamma \vec{x}\vec{\xi}_s
$$
Continuation
For the other possible symmetries you will have to wait, as I haven't computed them yet. Well, one and two can be writen as
$$
\frac{1}{\gamma}\partial_i(\gamma \xi^i)=0, \text{ no sum over } i
$$
and unfortunatley, to fulfil both at the same time we need
$$
\xi = \gamma^{-1}(f_1(y),f_2(x))
$$
and this doesn't fulfil whater the choice of $f_1,f_2$ se make. So it seams that we cant extract more variational simmetries this way. Because I know that there is a transformation that leaves the Lagrangian invariant:
$$
\Phi_{\alpha\in\mathbb R, a\in \mathbb C}(z) =\exp(i\alpha) \frac{z-a}{\bar{a}z-1}
$$
leaves the Lagrangian invariant, but has not appeared entirely, only the rotational parti, not the one dependent in $a$. This is a mistery.
Conclusion
There exists a method to find the symmetries of a Lagrangian, and it is cumbersome and involve huge PDE systems. For Lagrangian densities, check 1 and 3. Have fun and report if you find something interesting.
Update: And it might be that the equations are not integrable, as my previous case. So whe have an algorithm, but is cumbersome and might be that it still misses some. And this is strange, could it be the complex part.
Bibliography
Invariant Variational Problems, D.J. Logan, Elsevier
"Classical Noether’s theory with application to the
linearly damped particle", Raphaël Leone and Thierry Gourieux (LPM) arXiv:1412.7523v2 [math-ph]
Emmy Noether's Wonderful Theorem Dwight E. Neuenschwander, John Hopkings University Press
"Variational symmetries of Lagrangians", G.F. Torres del Castillo, C. Andrade Mirón, and R.I. Bravo Rojas, Rev. Mex. Fis. E 59(2) (2013) 140.