Let's consider the nondimensionalized Hamiltonian
$$\hat H=\frac{\hat{p}^2}2-\frac1r.\tag1$$
Its standard eigenfunctions diagonalize operators $\hat H$, $\hat L_z$ and $\hat L^2$. Laplace-Runge-Lenz operator can be defined as
$$\hat{\vec A}=\frac{\vec r}r-\frac12\left(\hat{\vec p}\times\hat{\vec L}-\hat{\vec L}\times\hat{\vec p}\right).\tag2$$
Its $z$-component $\hat A_z$ commutes with $\hat H$, $\hat L_z$, but doesn't commute with $\hat L^2$ nor with $\hat A^2$, although $\hat A^2$ does commute with $\hat H$, $\hat L_z$ and $\hat L^2$.
Since $\hat A^2$ commutes with the operators giving quantum numbers $n,l,m$ to the standard hydrogenic eigenstates, these standard eigenstates are also eigenstates of $\hat A^2$. So no new functions here. For these, we should look at $\hat A_z$.
2.1. Semi-numerical approach to finding exact eigenstates of $\hat A_z$
Since operators for $H$ and $L_z$ commute with $\hat A_z$, any eigenstate of $\hat A_z$ is a superposition of states with fixed $n$ and $m$ and different $l$ quantum numbers. This lets us find eigenstates of $\hat A_z$ for given $n,m$ without actually solving the eigenfunction PDE. We can e.g. numerically minimize variance of a sample of values of the expression for the eigenvalue
$$A_z=\frac{\hat A_z\sum\limits_l\alpha_l\psi_{n,l,m}}{\sum\limits_l\alpha_l\psi_{n,l,m}},\tag3$$
where $\psi_{n,l,m}$ are the standard simultaneous eigenfunctions of $\hat H$, $\hat L^2$ and $\hat L_z$, getting a set of weights $\alpha_l$. Then we can use the approximations we got to guess the exact values of these weights (and substitute into $(3)$, simplifying, to confirm the guess).
So, for $n=1$ we have the only eigenfunction, and it's obviously an eigenfunction of $\hat A_z$. For $n>1$ we have more basis functions, which lets us actually form more "interesting" eigenstates of $\hat A_z$. Here're some tables for $n=2,3,4$ (I calculated them using the above mentioned minimization procedure): $$n=2\\ \begin{array}{|c|c|c|} \hline \text{Eigenfunction of }\hat A_z&A_z&m\\ \hline \psi_{2,1,1} &0 &1 \\ \psi_{2,1,-1} &0 &-1 \\ \frac1{\sqrt2}(\psi_{2,0,0}+\psi_{2,1,0}) &-\frac12 & 0\\ \frac1{\sqrt2}(\psi_{2,0,0}-\psi_{2,1,0}) &\frac12 & 0\\ \hline \end{array}$$
$$n=3\\ \begin{array}{|c|c|c|} \hline \text{Eigenfunction of }\hat A_z&A_z&m\\ \hline \psi_{3,2,-2} & 0 & -2\\ \psi_{3,2,2} & 0 & 2\\ \frac1{\sqrt2}(\psi_{3,1,1}+\psi_{3,2,1}) & -\frac13 & 1\\ \frac1{\sqrt2}(\psi_{3,1,1}-\psi_{3,2,1}) & \frac13 & 1\\ \frac1{\sqrt2}(\psi_{3,1,-1}+\psi_{3,2,-1}) & -\frac13 & -1\\ \frac1{\sqrt2}(\psi_{3,1,-1}-\psi_{3,2,-1}) & \frac13 & -1\\ \sqrt{\frac13}\psi_{3,0,0}-\sqrt{\frac23}\psi_{3,2,0} & 0 & 0\\ \frac1{\sqrt3}\psi_{3,0,0}+\frac1{\sqrt2}\psi_{3,1,0}+\frac1{\sqrt6}\psi_{3,2,0} & -\frac23 & 0\\ \frac1{\sqrt3}\psi_{3,0,0}-\frac1{\sqrt2}\psi_{3,1,0}+\frac1{\sqrt6}\psi_{3,2,0} & \frac23 & 0\\ \hline \end{array}$$
$$n=4\\ \begin{array}{|c|c|c|} \hline \text{Eigenfunction of }\hat A_z&A_z&m\\ \hline \psi_{4,3,3} & 0 & 3\\ \psi_{4,3,-3} & 0 & -3\\ \frac1{\sqrt2}(\psi_{4,2,2}+\psi_{4,3,2}) & -\frac14 & 2\\ \frac1{\sqrt2}(\psi_{4,2,2}-\psi_{4,3,2}) & \frac14 & 2\\ \frac1{\sqrt2}(\psi_{4,2,-2}+\psi_{4,3,-2}) & -\frac14 & -2\\ \frac1{\sqrt2}(\psi_{4,2,-2}-\psi_{4,3,-2}) & \frac14 & -2\\ % \sqrt{\frac3{10}}\psi_{4,1,1}-\sqrt{\frac12}\psi_{4,2,1}+\sqrt{\frac15}\psi_{4,3,1} & \frac12 & 1\\ % \sqrt{\frac3{10}}\psi_{4,1,1}+\sqrt{\frac12}\psi_{4,2,1}+\sqrt{\frac15}\psi_{4,3,1} & -\frac12 & 1\\ % \sqrt{\frac3{10}}\psi_{4,1,-1}-\sqrt{\frac12}\psi_{4,2,-1}+\sqrt{\frac15}\psi_{4,3,-1} & \frac12 & -1\\ % \sqrt{\frac3{10}}\psi_{4,1,-1}+\sqrt{\frac12}\psi_{4,2,-1}+\sqrt{\frac15}\psi_{4,3,-1} & -\frac12 & -1\\ % \sqrt{\frac25}\psi_{4,1,1}-\sqrt{\frac35}\psi_{4,3,1} & 0 & 1\\ % \sqrt{\frac25}\psi_{4,1,-1}-\sqrt{\frac35}\psi_{4,3,-1} & 0 & -1\\ % \frac12\psi_{4,0,0}-\sqrt{\frac1{20}}\psi_{4,1,0}-\frac12\psi_{4,2,0}+\sqrt{\frac9{20}}\psi_{4,3,0} & \frac14 & 0\\ % \frac12\psi_{4,0,0}+\sqrt{\frac1{20}}\psi_{4,1,0}-\frac12\psi_{4,2,0}-\sqrt{\frac9{20}}\psi_{4,3,0} & -\frac14 & 0\\ % \frac12\psi_{4,0,0}+\sqrt{\frac9{20}}\psi_{4,1,0}+\frac12\psi_{4,2,0}+\sqrt{\frac1{20}}\psi_{4,3,0} & -\frac34 & 0\\ % \frac12\psi_{4,0,0}-\sqrt{\frac9{20}}\psi_{4,1,0}+\frac12\psi_{4,2,0}-\sqrt{\frac1{20}}\psi_{4,3,0} & \frac34 & 0\\ \hline \end{array}$$
The values of $m$ have been added to the table to make it more obvious that they, together with the eigenvalues $A_z$ and $n$, actually complete the CSCO. We can even plot configurations of the possible states using these quantum numbers, and see a pattern:
Looking at these plots, we can guess what the eigenvalues of $\hat A_z$ will be for the higher $n$. And then we won't even need the procedure of minimization of variance of the sample of $(3)$ to find the weights $\alpha_l$: we can simply solve this equation with a random sample of points in the domain, setting one of the weights to $1$, and those corresponding to the $m$ we're not interested in, to $0$ (number of points should be chosen to make number of equations equal to number of weights remaining unknown).
Now, $k$th value for $A_z$ appears to follow this formula:
$$A_z(n,m,k)=\frac{|m|-n-1+2k}n,\tag4$$ where
$$k=1,2,\dots,n-|m|.\tag5$$
2.2. Analytical approach to general solution
The general solution to eigenproblem of $\hat A_z$ can be found, if the Schrödinger's equation is expressed in the parabolic coordinates. Then the natural eigenfunctions there will be characterized by a different set of quantum numbers than usual: parabolic ones $n_1$ and $n_2$, and the usual magnetic quantum number $m$. The complete expression in terms of confluent hypergeometric function ${}_1F_1$ for the eigenfunctions and its derivation can be found e.g. in Landau and Lifshitz, "Quantum Mechanics. Non-relativistic theory" $\S37$ "Motion in a Coulomb field (parabolic coordinates)". The eigenvalues of $\hat A_z$ can be expressed in terms of the parabolic and magnetic quantum numbers as
$$A_z=\frac{n_1-n_2}n,\tag6$$
where
$$n=n_1+n_2+|m|+1\tag7$$
is the principal quantum number, and parabolic quantum numbers can have the values
$$n_{1,2}=0,1,...,n-|m|-1.\tag8$$
In terms of $n$, $m$ and $n_1$, $(6)$ can be rewritten as
$$A_z=\frac{|m|-n+2n_1+1}n,\tag9$$
which is consistent with $(4)$-$(5)$.
Eigenfunctions of $\hat A_z$ with high absolute values of eigenvalues look like bells in shape, with probability density symmetric along the $z$ axis. As $L_z$ is conserved, real and imaginary parts of the eigenfunctions oscillate when we go around the $z$ axis with $m\ne0$.
Here're plots of some of the eigenfunctions, with 3D density plot on the LHS and cross-sections in $y=0$ plane on the RHS:
$n=4,\,A_z=\frac34,\,m=0:$
$n=4,\,A_z=-\frac14,\,m=0:$
$n=4,\,A_z=-\frac14,\,m=2$, real part:
$n=3,\,A_z=\frac13,\,m=1$, real part:
4.1. Eigenfunctions of $\hat A_z$
Unfortunately, the eigenstates of $\hat A_z$ don't appear to be nicely related to classical eccentric orbits. In classical orbits, $\vec A$ is always in the plane of rotation. Classical orbital motion corresponds in quantum regime to the case of high expected values of $L_z$. But $\hat A_z$ commutes with $\hat L_z$, not with e.g. $\hat L_x$, so the "rotating" (in the sense of $e^{im\phi}$) eigenstates rotate roughly perpendicularly to direction of the orbital eccentricity. We can interpret these as eccentric standing waves in the $\theta$ direction, which is of course far from classical regime.
It's interesting to note that although with the high values of $A_z$ the $z$ component of $\vec A$ indeed dominates, giving the cross-section of the orbital somewhat elliptic shape, it's much less so for lower values. This is because of uncertainty in $A_x$ and $A_z$. See e.g. the following eigenfunctions. Ellipses show the supposed classical orbits with semi-major axis $a=n^2$ and LRL vectors $\vec A=A_z \vec e_z+s \sqrt{\langle A_x^2\rangle}\vec e_x$, where $s=-2,-1,0,1,2$.
$n=20,\, A_z=\frac{19}{20},\, m=0,\,\langle A_x^2\rangle=\frac{19}{800}\colon$
$n=20,\, A_z=\frac{9}{20},\, m=0,\,\langle A_x^2\rangle=\frac{159}{800}\colon$
$n=20,\, A_z=\frac{1}{20},\, m=0,\,\langle A_x^2\rangle=\frac{199}{800}\colon$
4.2. Eigenfunctions of $\hat A^2$
We may have better luck if we consider instead eigenstates of $\hat A^2$ (which, as mentioned above, doesn't commute with $\hat A_z$). These are also eigenstates of $\hat L^2$ and $\hat L_z$, so they are the familiar functions. Eigenvalues of $\hat A^2$ are
$$A_{n,l}^2=1-\frac{l(l+1)+1}{n^2}.$$
As we can see, consistently with classical intuition, eccentricity $e=\sqrt{A^2}$ decreases with increasing angular momentum. One might hope to see that at least these states will allow assignment of classical eccentricity. But, despite they do, there's a problem: the eigenstates of $\hat A^2$ are all almost symmetric with respect to rotations around $z$ axis — modulo the $\exp(im\phi)$ oscillations. So we never get anything resembling eccentric ellipses even in eigenstates of $\hat A^2$. Instead we get the following cross-sections in $xy$ plane of the real parts of wavefunctions (here $n=21$, $m=l$, $l$ changes from $0$ to $20$ with the step of $4$):
Since these states don't actually have any orientation of the LRL vector (aside from avoiding $z$ direction in case of high $m$ values), a better interpretation and drawing of classical eccentricity for them would be like this:
where the ellipses show some of the possible orbits one may get if e.g. one were to form a localized wave packet from similar states with this state being dominant.
To make the numerical procedure more understandable, I'll show an example of how one can get an eigenfunction of $\hat A_z$ and associated eigenvalue with $n=4$, $m=1$. The code in this section is in Wolfram Language. I did all the calculations in Mathematica 11.2, but the code is compatible with versions as old as Mathematica 9.
First, some definitions for the operators and functions we'll use here.
(* Components of momentum operator *)
px = -I D[#,x] &;
py = -I D[#,y] &;
pz = -I D[#,z] &;
(* z component of LRL vector operator *)
Az = Simplify[
z/Sqrt[x^2+y^2+z^2] # -
1/2 (z px@px@# + px[z px@#] - px[x pz@#] + z py@py@# +
py[z py@#] - py[y pz@#] - x pz@px@# - y pz@py@#)] &;
(* Hydrogenic wavefunction in spherical coordinates *)
ψ[n_,l_,m_,r_,θ_,ϕ_] = Sqrt[(n-l-1)!/(n+l)!] E^(-r/n) (2 r/n)^l 2/n^2 *
LaguerreL[n-l-1, 2 l+1, (2 r)/n] *
SphericalHarmonicY[l,m,θ,ϕ];
(* The same wavefunction converted to Cartesian coordinates *)
Ψ[n_,l_,m_,x_,y_,z_] = ψ[n, l, m, Sqrt[x^2+y^2+z^2],
ArcCos[z/Sqrt[x^2+y^2+z^2]],
ArcTan[x,y]];
Now the example test function for $n=4$, $m=1$.
(* Test function with parameters α, β and γ. Restricting arguments to
numeric to avoid attempts at symbolic evaluation, which can seriously slow
things down. Simplifying it to speedup calculations and reduce roundoff
errors. *)
test[x_Real,y_Real,z_Real,α_?NumericQ,β_?NumericQ,γ_?NumericQ] =
FullSimplify[Az@#/# &[α Ψ[4,1,1,x,y,z] +
β Ψ[4,2,1,x,y,z] +
γ Ψ[4,3,1,x,y,z]],
(x|y|z) ∈ Reals];
And finally minimization of its variance. Note that we don't need to actually calculate an integral as in variational methods: we only need a "close enough" approximation of the parameters, the rest can be left to Rationalize. So we use a coarse mesh of points to evaluate the function on. Note also that here we let Mathematica go to complex domain, although the parameters should be real-valued. This lets it avoid singularities in the function by simply going around them, and thus gives much faster convergence.
(* Table is generated not on integers to avoid problems like
division by zero on evaluation *)
With[{
var = Total[Abs[#-Mean@#]^2]& @ Flatten @
Table[test[x,y,z,α,βR + I βI,γR + I γI],
{x,-10.123,10,4},
{y,-10.541,10,5},
{z,-10.07,10,5}
]
},
{minVal,minim} = NMinimize[{var,Total[#^2] == 1 &[{α,βR,βI,γR,γI}] && α>0.1},
{α,βR,βI,γR,γI}]
]
{1.29281973036898*10^-9, {α -> 0.547724712837901, βR -> -0.707106153522197, βI -> 1.84829862406368*10^-7, γR -> 0.447211968838118, γI -> -1.8807768960726*10^-7}}
OK, so we see that indeed the imaginary parts are close to zero, so let's guess the form of actual parameters assuming that what we got are square roots of some rationals.
(* Tolerance of rationalization is chosen so at to
1) ignore numerical errors of minimization,
2) still give a good enough room to guess the correct number *)
Sqrt[Rationalize[#^2, 10^-4]]Sign[#]&[{α,βR,γR} /. minim]
{Sqrt[3/10], -(1/Sqrt[2]), 1/Sqrt[5]}
This is what we got. Let's check whether this is a correct guess.
FullSimplify[Az@#/# &[Sqrt[3/10] Ψ[4,1,1,x,y,z] -
Sqrt[1/2] Ψ[4,2,1,x,y,z] +
Sqrt[1/5] Ψ[4,3,1,x,y,z]],
(x|y|z) ∈ Reals]
1/2
Now we not only confirmed that our test function with the guessed values of parameters is an eigenfunction (since we got constant here), but also obtained the associated eigenvalue, $A_z=1/2$. This is entry #7 in the table above for $n=4$.
To find another set of parameters we can go the following way. First, we can guess that changing some signs in $\alpha$, $\beta$ and $\gamma$ might give us some more eigenfunctions. Indeed, it does, so using $+\sqrt{1/2}$ instead of $-\sqrt{1/2}$ for $\beta$ does result in an eigenfunction (entry #8 in the table, with eigenvalue $A_z=-1/2$).
Another approach at finding other eigenfunctions (useful when there are more parameters, e.g. for $n=4$, $m=0$ there are $4$) is using Orthogonalize to find a basis in the orthogonal subspace of parameters to the one we've already identified. Then we can use that basis to form our new set of parameters for NMinimize to work on. In our example case the situation is trivial, since the whole set of $n=4$, $m=1$ eigenfunctions consists of 3 elements, and we've already identified two of them, so no need in further NMinimize. So
Orthogonalize[{{Sqrt[3/10], -(1/Sqrt[2]), 1/Sqrt[5]},
{Sqrt[3/10], 1/ Sqrt[2], 1/Sqrt[5]},
{1, 1, 1}}] // FullSimplify
{{Sqrt[3/10], -(1/Sqrt[2]), 1/Sqrt[5]}, {Sqrt[3/10], 1/Sqrt[2], 1/Sqrt[5]}, {-Sqrt[(2/5)], 0, Sqrt[3/5]}}
The third element of the output list is the third eigenfunction (in the $|n,l,m\rangle$ basis). We can find that associated eigenvalue is $A_z=0$.
