What symmetry causes the Runge-Lenz vector to be conserved?

1) Hamiltonian Problem. The Kepler problem has Hamiltonian

$$ H~=~T+V, \qquad T~:=~ \frac{p^2}{2m}, \qquad V~:=~- \frac{k}{q}, \tag{1} $$

where $m$ is the 2-body reduced mass. The Laplace–Runge–Lenz vector is (up to an irrelevant normalization)

$$ A^j ~:=~a^j + km\frac{q^j}{q}, \qquad a^j~:=~({\bf L} \times {\bf p})^j~=~{\bf q}\cdot{\bf p}~p^j- p^2~q^j,\qquad {\bf L}~:=~ {\bf q} \times {\bf p}.\tag{2}$$

2) Action. The Hamiltonian Lagrangian is

$$ L_H~:=~ \dot{\bf q}\cdot{\bf p} - H,\tag{3} $$

and the action is

$$ S[{\bf q},{\bf p}]~=~ \int {\rm d}t~L_H .\tag{4}$$

The non-zero fundamental canonical Poisson brackets are

$$ \{ q^i , p^j\}~=~ \delta^{ij}. \tag{5}$$

3) Inverse Noether's Theorem. Quite generally in the Hamiltonian formulation, given a constant of motion $Q$, then the infinitesimal variation

$$\delta~=~ -\varepsilon \{Q,\cdot\}\tag{6}$$

is a global off-shell symmetry of the action $S$ (modulo boundary terms). Here $\varepsilon$ is an infinitesimal global parameter, and $X_Q=\{Q,\cdot\}$ is a Hamiltonian vector field with Hamiltonian generator $Q$. The full Noether charge is $Q$, see e.g. my answer to this question. (The words on-shell and off-shell refer to whether the equations of motion are satisfied or not. The minus is conventional.)

4) Variation. Let us check that the three Laplace–Runge–Lenz components $A^j$ are Hamiltonian generators of three continuous global off-shell symmetries of the action $S$. In detail, the infinitesimal variations $\delta= \varepsilon_j \{A^j,\cdot\}$ read

$$ \delta q^i ~=~ \varepsilon_j \{A^j,q^i\} , \qquad \{A^j,q^i\} ~=~ 2 p^i q^j - q^i p^j - {\bf q}\cdot{\bf p}~\delta^{ij}, $$ $$ \delta p^i ~=~ \varepsilon_j \{A^j,p^i\} , \qquad \{A^j,p^i\}~ =~ p^i p^j - p^2~\delta^{ij} +km\left(\frac{\delta^{ij}}{q}- \frac{q^i q^j}{q^3}\right), $$ $$ \delta t ~=~0,\tag{7}$$

where $\varepsilon_j$ are three infinitesimal parameters.

5) Notice for later that

$$ {\bf q}\cdot\delta {\bf q}~=~\varepsilon_j({\bf q}\cdot{\bf p}~q^j - q^2~p^j), \tag{8} $$

$$ {\bf p}\cdot\delta {\bf p} ~=~\varepsilon_j km(\frac{p^j}{q}-\frac{{\bf q}\cdot{\bf p}~q^j}{q^3})~=~- \frac{km}{q^3}{\bf q}\cdot\delta {\bf q}, \tag{9} $$

$$ {\bf q}\cdot\delta {\bf p}~=~\varepsilon_j({\bf q}\cdot{\bf p}~p^j - p^2~q^j )~=~\varepsilon_j a^j, \tag{10} $$

$$ {\bf p}\cdot\delta {\bf q}~=~2\varepsilon_j( p^2~q^j - {\bf q}\cdot{\bf p}~p^j)~=~-2\varepsilon_j a^j~. \tag{11} $$

6) The Hamiltonian is invariant

$$ \delta H ~=~ \frac{1}{m}{\bf p}\cdot\delta {\bf p} + \frac{k}{q^3}{\bf q}\cdot\delta {\bf q}~=~0, \tag{12}$$

showing that the Laplace–Runge–Lenz vector $A^j$ is classically a constant of motion

$$\frac{dA^j}{dt} ~\approx~ \{ A^j, H\}+\frac{\partial A^j}{\partial t} ~=~ 0.\tag{13}$$

(We will use the $\approx$ sign to stress that an equation is an on-shell equation.)

7) The variation of the Hamiltonian Lagrangian $L_H$ is a total time derivative

$$ \delta L_H~=~ \delta (\dot{\bf q}\cdot{\bf p})~=~ \dot{\bf q}\cdot\delta {\bf p} - \dot{\bf p}\cdot\delta {\bf q} + \frac{d({\bf p}\cdot\delta {\bf q})}{dt} $$ $$ =~ \varepsilon_j\left( \dot{\bf q}\cdot{\bf p}~p^j - p^2~\dot{q}^j + km\left( \frac{\dot{q}^j}{q} - \frac{{\bf q} \cdot \dot{\bf q}~q^j}{q^3}\right)\right) $$ $$- \varepsilon_j\left(2 \dot{\bf p}\cdot{\bf p}~q^j - \dot{\bf p}\cdot{\bf q}~p^j- {\bf p}\cdot{\bf q}~\dot{p}^j \right) - 2\varepsilon_j\frac{da^j}{dt}$$ $$ =~\varepsilon_j\frac{df^j}{dt}, \qquad f^j ~:=~ A^j-2a^j, \tag{14}$$

and hence the action $S$ is invariant off-shell up to boundary terms.

8) Noether charge. The bare Noether charge $Q_{(0)}^j$ is

$$Q_{(0)}^j~:=~ \frac{\partial L_H}{\partial \dot{q}^i} \{A^j,q^i\}+\frac{\partial L_H}{\partial \dot{p}^i} \{A^j,p^i\} ~=~ p^i\{A^j,q^i\}~=~ -2a^j. \tag{15}$$

The full Noether charge $Q^j$ (which takes the total time-derivative into account) becomes (minus) the Laplace–Runge–Lenz vector

$$ Q^j~:=~Q_{(0)}^j-f^j~=~ -2a^j-(A^j-2a^j)~=~ -A^j.\tag{16}$$

$Q^j$ is conserved on-shell

$$\frac{dQ^j}{dt} ~\approx~ 0,\tag{17}$$

due to Noether's first Theorem. Here $j$ is an index that labels the three symmetries.

9) Lagrangian Problem. The Kepler problem has Lagrangian

$$ L~=~T-V, \qquad T~:=~ \frac{m}{2}\dot{q}^2, \qquad V~:=~- \frac{k}{q}. \tag{18} $$

The Lagrangian momentum is

$$ {\bf p}~:=~\frac{\partial L}{\partial \dot{\bf q}}~=~m\dot{\bf q} \tag{19} . $$

Let us project the infinitesimal symmetry transformation (7) to the Lagrangian configuration space

$$ \delta q^i ~=~ \varepsilon_j m \left( 2 \dot{q}^i q^j - q^i \dot{q}^j - {\bf q}\cdot\dot{\bf q}~\delta^{ij}\right), \qquad\delta t ~=~0.\tag{20}$$

It would have been difficult to guess the infinitesimal symmetry transformation (20) without using the corresponding Hamiltonian formulation (7). But once we know it we can proceed within the Lagrangian formalism. The variation of the Lagrangian is a total time derivative

$$ \delta L~=~\varepsilon_j\frac{df^j}{dt}, \qquad f_j~:=~ m\left(m\dot{q}^2q^j- m{\bf q}\cdot\dot{\bf q}~\dot{q}^j +k \frac{q^j}{q}\right)~=~A^j-2 a^j . \tag{21}$$

The bare Noether charge $Q_{(0)}^j$ is again

$$Q_{(0)}^j~:=~2m^2\left(\dot{q}^2q^j- {\bf q}\cdot\dot{\bf q}~\dot{q}^j\right) ~=~-2a^j . \tag{22}$$

The full Noether charge $Q^j$ becomes (minus) the Laplace–Runge–Lenz vector

$$ Q^j~:=~Q_{(0)}^j-f^j~=~ -2a^j-(A^j-2a^j)~=~ -A^j,\tag{23}$$

similar to the Hamiltonian formulation (16).

While Kepler second law is simply a statement of the conservation of angular momentum (and as such it holds for all systems described by central forces), the first and the third laws are special and are linked with the unique form of the newtonian potential $-k/r$. In particular, Bertrand theorem assures that only the newtonian potential and the harmonic potential $kr^2$ give rise to closed orbits (no precession). It is natural to think that this must be due to some kind of symmetry of the problem. In fact, the particular symmetry of the newtonian potential is described exactly by the conservation of the RL vector (it can be shown that the RL vector is conserved iff the potential is central and newtonian). This, in turn, is due to a more general symmetry: if conservation of angular momentum is linked to the group of special orthogonal transformations in 3-dimensional space $SO(3)$, conservation of the RL vector must be linked to a 6-dimensional group of symmetries, since in this case there are apparently six conserved quantities (3 components of $L$ and 3 components of $\mathcal A$). In the case of bound orbits, this group is $SO(4)$, the group of rotations in 4-dimensional space.

Just to fix the notation, the RL vector is:

\begin{equation} \mathcal{A}=\textbf{p}\times\textbf{L}-\frac{km}{r}\textbf{x} \end{equation}

Calculate its total derivative:

\begin{equation}\frac{d\mathcal{A}}{dt}=-\nabla U\times(\textbf{x}\times\textbf{p})+\textbf{p}\times\frac{d\textbf{L}}{dt}-\frac{k\textbf{p}}{r}+\frac{k(\textbf{p}\cdot \textbf{x})}{r^3}\textbf{x} \end{equation}

Make use of Levi-Civita symbol to develop the cross terms:

\begin{equation}\epsilon_{sjk}\epsilon_{sil}=\delta_{ji}\delta_{kl}-\delta_{jl}\delta_{ki} \end{equation}

Finally:

\begin{equation} \frac{d\mathcal{A}}{dt}=\left(\textbf{x}\cdot\nabla U-\frac{k}{r}\right)\textbf{p}+\left[(\textbf{p}\cdot\textbf{x})\frac{k}{r^3}-2\textbf{p}\cdot\nabla U\right]\textbf{x}+(\textbf{p}\cdot\textbf{x})\nabla U \end{equation}

Now, if the potential $U=U(r)$ is central:

\begin{equation} (\nabla U)_j=\frac{\partial U}{\partial x_j}=\frac{dU}{dr}\frac{\partial r}{\partial x_j}=\frac{dU}{dr}\frac{x_j}{r} \end{equation}

so

\begin{equation} \nabla U=\frac{dU}{dr}\frac{\textbf{x}}{r}\end{equation}

Substituting back:

\begin{equation}\frac{d\mathcal A}{dt}=\frac{1}{r}\left(\frac{dU}{dr}-\frac{k}{r^2}\right)[r^2\textbf{p}-(\textbf{x}\cdot\textbf{p})\textbf{x}]\end{equation}

Now, you see that if $U$ has exactly the newtonian form then the first parenthesis is zero and so the RL vector is conserved.

Maybe there's some slicker way to see it (Poisson brackets?), but this works anyway.

The symmetry is an example of an open symmetry, i.e. a symmetry group which varies from group action orbit to orbit. For bound trajectories, it's SO(4). For parabolic ones, it's SE(3). For hyperbolic ones, it's SO(3,1). Such cases are better handled by groupoids.