Here's a stab at what is going on. A note upfront: The approach I lay out here only gives the sames answers up to the $f$ orbitals as shown at Wikipedia and The Orbitron. It is just the best I could come up with. I would still greatly appreciate any other answers that can further resolve the issue.
Following the quantum mechanics convention on wikipedia, the complex spherical harmonics are given by
\begin{align}
Y_l^m(\theta, \phi) =& (-1)^m \sqrt{\frac{(2l + 1)}{4\pi}\frac{(l-m)!}{(l+m)!}}P_l^m(\cos\theta) e^{im\phi}\\
=& N_{lm}P_l^m(\cos(\theta))e^{im\phi}
\end{align}
Where $N_{lm}$ is a scalar normalization factor and $P_l^m$ is an Associated Legendgre Polynomial.
These polynomials have a closed form
\begin{align}
P_l^m(x) =& (-1)^m 2^l (1-x^2)^{m/2} \sum_{k=m}^l \frac{k!}{(k-m)!} \binom{l}{k}\binom{\frac{l+k-1}{2}}{l} x^{k-m}\\
=& M_{lm} (1-x^2)^{m/2} G_{l-m}(x)\\
=& M_{lm}K_{lm} (1-x^2)^{m/2} x^{l-m} H_{l-m-1}
\end{align}
Where $M_{lm}$ is a normalization constant, $G_{l-m}$ is a polynomial of order $l-m$ and
$$
G_{l-m} = K_{lm}(x^{l-m} + H_{l-m-1}(x))
$$
Where $K_{lm}$ is a constant and $H_{l-m-1}$ is a polynomial of order $l-m-1$.
If $x = \cos(\theta)$ then $(1-x^2)^{m/2} = (\sin(\theta))^m$. Let
$$
C_{lm} = N_{lm}M_{lm}K_{lm}
$$
We then have
$$
Y_l^m(\theta, \phi) = C_{lm} e^{im\phi}(\sin(\theta))^m\left((\cos(\theta))^{l-m} + H_{l-m-1}(\cos\theta)\right)
$$
We now re-express in cartesian components. For that we have
\begin{align}
e^{i\phi} =& \frac{x+i y}{\sqrt{x^2 + y^2}}\\
\sin(\theta) =& \frac{\sqrt{x^2 + y^2}}{\sqrt{x^2+y^2+z^2}} = \frac{\sqrt{x^2 + y^2}}{r}\\
\cos(\theta) =& \frac{z}{\sqrt{x^2 + y^2 + z^2}} = \frac{z}{r}
\end{align}
Importantly, we can see that $e^{i\phi}\sin(\theta) = (x+iy)/r$.
With this we get
\begin{align}
Y_l^m(x, y, z) =& C_{lm} \frac{(x+iy)^m}{r^m}\left(\frac{z}{r^{l-m}} + H_{l-m-1}\left(\frac{z}{r}\right)\right)\\
=& \frac{C_{lm}}{r^l}(x+iy)^m\left(z^{l-m} + r^{l-m}H_{l-m-1}\left(z/r\right)\right)
\end{align}
From this complex spherical harmonic we can extract two real spherical harmonics. These will be related to (the signs of the actual harmonics may differ from the expressions below because of various conventions)
\begin{align}
&\frac{C_{lm}}{r^l}\text{Re}\left\{(x+iy)^m\right\}\left(z^{l-m} + r^{l-m}H_{l-m-1}\left(z/r\right)\right)\\
&\frac{C_{lm}}{r^l}\text{Im}\left\{(x+iy)^m\right\}\left(z^{l-m} + r^{l-m}H_{l-m-1}\left(z/r\right)\right)\\
\end{align}
These are the expressions from which extract the subscript for the atomic orbital name.
Focusing on $(x+iy)^m$. For every term in this polynomial the exponents for $x$ and $y$ add up to $m$.
$\text{Re/Im}\left\{(x+iy)^m\right\}$ will in general include terms with both $x$ and $y$ but there is one term (or zero terms) that has the biggest exponent for $x$ and one term (or zero terms) that has the biggest exponent for $y$.
The procedure for coming up with the subscript is as follows.
- $z^{l-m}$ will appear in the subscript. This means we ignore all of $H_{l-m-1}(z/r)$.
- Take $\text{Re/Im}\left\{(x+iy)^m\right\}$ depending on which orbital you are considering. Find the term with the biggest $x$ exponent and the term with the biggest $y$ exponent and add these two terms together. You will have something like $Ax^ay^b + B x^cy^d$. If one of these terms is missing then one of $A$ or $B$ is zero.
- Now take $Ax^ay^b + B x^c y^d$ and factor out any possible scalars and factors of $x$ and $y$. You will now have something like $A x^a y^b(x^c + B y^c)$ or $A x^a y^b(Bx^c + y^c)$ (these are not the same $A, B, a, b, c$ as just before).
- Strip off $A$ from your final expression so you now have $x^a y^b (B x^c + y^c)$. This expression will appear in the subscript. There are two parts here, the factored out part $x^a y^b$ and the remaining prat $(Bx^c + y^c)$ or $(x^c + By^c)$
- We now put together the three parts.
- First put in the factored out part like $x^a y^b$ with $x$ preceding $y$.
- Then put the $z^{l-m}$ part.
- Finally put the remaining part $(B x^c + y^c)$ or $(x^c + B y^c)$.
- This should result in a subscript $x^a y^b z^{l-m} (B x^c + y^c)$ or $x^a y^b z^{l-m} (x^c + B y^c)$
A few notes:
First, and importantly, this scheme does not give the same answer as Wikipedia for two orbitals listed on that page.
- Wikipedia says $Y_{4, -3} = g_{zy^3} \propto (3x^2-y^2) yz$. I would say this should be $g_{yz(3x^2-y^2)}$.
- Similarly, Wikipedia says $Y_{4, 3} = g_{zx^3} \propto (x^2-3y^2)xz$. I would say this should be $g_{xz(x^2-3y^2)}$.
- My ordering also differs from that on Wikipeda. But honestly, i don't think they were following a pattern for the order on Wikipedia.
Wikipedia uses the same notation as The Orbitron, and I wonder if Wikipedia got their notation from this webpage. The Orbitron goes out to $i$ orbitals and I see similar discrepencies between the scheme used on that page and my approach for $g$ and $h$ orbitals. The author stops quoting abbreviated polynomials for $i$ orbitals for the relevant $m$ quantum numbers so I can't compare past that. I've reached out to the author of the Orbitron for further comment and I'm also trying to reach out to Wikipedia folks.
One advantage of my approach over whatever is going on on Wikipedia and the Orbitron is that the $x, y$ part of the coefficient is the same for orbitals with the same $m$ quantum number. They only differ in the factor of $z$.
Note that $a + b + c + l-m = l$, so all the exponents always add up to $l$ which is a nice spherical harmonic like result.
Finally, the issue remains of determining the $a$, $b$, $c$ and $B$ factors arising from $\text{Re/Im}\left\{(x+iy)^m\right\}$ alone. A computer could work these out without too much trouble but it would be nice to have closed form expressions for them. It should be possible to work it out using the binomial theorem.
Finally an alternative approach would be to simply take the terms with the biggest exponents in $\text{Re\Im}\{(x+iy)^m\}$ and add them if there max $x$ and $y$ exponents are the same, but this also doesn't agree with Wikipedia as it would replace something like $f_{y(3x^2-y^2)}$ with $f_{y^3}$, but it does explain wikipedia terms like $g_{zy^3}\propto (3x^2-y^2)yz$.
More notes: I've done a little more research and I really can't find any authoritative reference that gives names for orbitals higher than $f$. If anyone is familiar with one please let me know. I'm going to personally consider orbitals $g$ and higher to NOT have official cartesian polynomial abbreviation names and pursue and advertise the $l, m$ notation as being the proper way to refer to these higher orbitals. More motivation for this below.
For orbitals $f$ and below the $xy$ part of "abbreviated polynomial" subscript DOES allow you to identify the azimuthal nodal planes, for example, $p_x$ is zero when $x=0$ and $d_{xy}$ is zero when $x=0$ or $y=0$. This even holds for the intimidating $f_{x(x^2-3y^2)}$ which is zero if $x=0, \pm \sqrt{3}y$. However, when you move to orbitals with $m>4$ (for example $g$ orbitals with $(l, m) = (4, \pm)$), there are now 3 nodal planes, and any subscript that identifies all of them is going to require multiple terms and it is really not convenient to have expressions like that. It is likely for this reason authors don't typically give subscripts for orbitals at this order.
Frankly, for anyone reading this, my opinion is that the $l, m$ notation is far superior above basically the $p$ orbitals. The number $|m|$ tells you how many azimuthal nodal planes there are and the number $l - |m|$ tells you how many polar nodla planes there are. If $m$ is positive then the orbital has a set of lobes aligned along the $x$ axis. if $m$ is negative then the orbital is just like the one with $+m$ but rotated about the $z$ axis by $\pm 360/(4|m|)$ degrees ($\pm$ depending on if $m$ is positive or negative by the Condon-Shortley phase). No need to abbreviate a polynomial to convey this geometric information when it is already communicated by the $l, m$ numbers when properly understood.
