Orbital angular momentum as sum of harmonic oscillators On section 7.3 of Ballentine's "Quantum Mechanics: A Modern Development" there is a really nice argument on why the eigenvalues of the total angular momentum operator must be integer, cf. e.g. this Phys.SE answer by NessunDorma. By defining the operators
$$q_1=\frac{Q_x+P_y}{\sqrt{2}},\quad q_2=\frac{Q_x-P_y}{\sqrt{2}}, \quad p_1=\frac{P_x-Q_y}{\sqrt{2}}, \quad p_2=\frac{P_x+Q_y}{\sqrt{2}},$$
the projection of the orbital angular momentum can be written as
$$L_z=Q_xP_y-Q_yP_x=\frac{1}{2}(p_1^2+q_1^2)-\frac{1}{2}(p_2^2+q_2^2),$$
i.e., a difference of harmonic oscillators. The reason we can say these are indeed independent harmonic oscillators is that the following commutation relations are satisfied:
$$[q_1,q_2]=[p_1,p_2]=0, \quad [q_a,p_b]=i\delta_{ab}.$$
Expressing this operator in terms of ladder operators, we can easily see that its eigenvalues would be of the form $(n_1+1/2)-(n_2+1/2)=n_1-n_2$, which is always an integer. This means the quantum number $m$ is an integer, which in turn implies $l$ is an integer.
If we make the change $p_2,q_2 \rightarrow iq_2,ip_2$, then we have
$$L_z=\frac{1}{2}(p_1^2+q_1^2)+\frac{1}{2}(p_2^2+q_2^2),$$
and the commutation relations are still satisfied (if we make instead $p_2,q_2 \rightarrow ip_2,iq_2$, then there is an unwanted minus sign in the expression for $[q_a,p_b]$). As a footnote, I think this is unphysical since $ip_2$ and $iq_2$ are not hermitian, but probably the following still makes mathematical sense:
We know that the eigenfunctions of $L_z$ are spherical harmonics $Y_l^m(\theta,\phi)$. On the other hand, the eigenfunctions of the harmonic oscillator are Hermite polynomials (with some factors).
My questions are:
If we write $L_z$ in terms of the operators $\{q_i,p_i\}$, can we say the eigenfunctions are products of Hermite polynomials (2D harmonic oscillator)? And if so, would this give a relation between the product of Hermite polynomials and spherical harmonics by means of a change of variables?
I hope there are no trivial mistakes in my reasoning.
 A: Ach, yes & no, but... this is the most ponderous summary of the classic Jordan map construction of SU(2) matrix generators there is.
Most of the difficulty and confusion here lies in the change of bases, which runaway abstraction of notation fosters. The notional spaces of the two oscillators are not the spacetime indexed by the generators of rotations, so our everyday $\theta, \phi$ space angles! They are oscillating in auxiliary spaces similar to those of QFT oscillators, likewise unconnected to our spacetime!
I will thus illustrate all this with a specific example: let's look at the $\ell =2$, spin two representation consisting of the three $2\ell +1=5$ -dimensional matrices, 5×5, satisfying $[L^j,L^k]= i \epsilon^{jkr}L^r$, where I've nondimensionalized $\hbar$ for sanity.
In Fock space, life is simple: spin two is but the symmetric spin addition of four spin-1/2s, two $a_1$ and two $a_2$ oscillators, so $n_1=2, ~ n_2=2$.
The building block of the 5×5 matrices is the image of the  three 2×2 Pauli matrices  in this map,
$${\vec  L} \equiv    {\mathbf  a}^\dagger \cdot\frac{ {\vec  \sigma } } {2} \cdot    {\mathbf a}^{\,}   ~,$$
for two-vectors ${\mathbf a},{ \mathbf a}^\dagger$,
the starting point of Schwinger’s 1952 treatment of the theory of quantum angular momentum, predicated on the action of these operators on  Fock states built of arbitrary higher powers of such operators.
Here the two oscillators are tensor-squared to
$$L^2\equiv {\vec  L} \cdot {\vec  L} = \frac{n_1+n_2}{2} \left ( \frac{n_1+n_2}{2}+1\right )I_5~, 
$$
the eigenvalue here being 6, as expected for spin 2, of course.
For instance, acting on an (unnormalized) Fock eigenstate,  recalling   $L_+=a_1^\dagger a_2$ and  $L_-=a_2^\dagger a_1$, observe
$$L^2~   a^{\dagger k}_1  a^{\dagger n}_2  |0\rangle=   \frac{k+n}{2} \left ( \frac{k+n}{2}+1\right ) ~  a^{\dagger k}_1  a^{\dagger n}_2 |0\rangle ~,$$
while
$$L_z ~   a^{\dagger k}_1  a^{\dagger n}_2  |0\rangle=   \frac{1}{2} \left ( k-n\right ) ~  a^{\dagger k}_1  a^{\dagger n}_2 |0\rangle ~,$$
so that, for $l=(k+n)/2$  ,   $m=  (k−n)/2$, this is proportional to the eigenstate $ |l,m\rangle$,
$$|l,m\rangle= \frac{a_1^{\dagger ~(l+m)} a_2^{\dagger ~ (l-m)} }{\sqrt{(l+m)!~(l-m)!}}|0\rangle~  . $$
In our example, $l=2$, so
$$|2,m\rangle= \frac{a_1^{\dagger ~2+m} a_2^{\dagger ~ 2-m} }{\sqrt{(2+m)!~(2-m)!}}|0\rangle~  . $$
Now recall the standard bases connection to now two disjoint one dimensional spaces (!)
$$
\langle x_j|a_i^{\dagger ~~n}|0\rangle = \delta_{ij}\psi_n^{(i)}(x_i), 
$$
where the $\psi_n (x)$ are the Hermite functions related to the Hermite polynomials in the customary way, n indicates the energy excitation index, and (i) labels the oscillator 1,2, they came from, but now they are the same functions.   So, e.g., take m=2,
$$
(\langle x_1| \otimes \langle x_2 |) |2,2\rangle= \psi_4 (x_1)\psi_0(x_2),
$$
etc.
Likewise, see here
$$
\langle \theta,\phi| l,m\rangle= Y_l^m(\theta, \phi),
$$
with this explicit form, for the rest
Consider $Y_2^2(\theta,\phi)= \frac{1}{4}  \sqrt{\frac{15}{2\pi} } \sin^2\theta ~ e^{2i\phi}$. This is a function with arguments on the celestial sphere, so to speak, and having little to do with the two notional spaces in which the Hermite functions of the two Schwinger oscillators live.

*

*Indeed, the two very different realizations describe the same states and the same 5×5 matrices (in our illustration), through the employment of motions on two very different manifolds; but these hardly induce a natural connection between associated Legendre polynomials and Hermite polynomials (and hence functions, above).

NB. Still, there is an important connection between Laguerre and Hermite polynomials, but that's yet another construction... phase space, where the Laguerres of Wigner functions are bilinears of the Hermites of the wave functions. Totally outranges your focus here.

Note in response to comment by @ytlu
The above review/deconstruction was meant to highlight the tenuousness of the connection between Hermite polynomials and spherical harmonics. The five states I mapped above, ($y\equiv \cos\theta$),
$$
Y^2_2 \sim e^{i2\phi} P_2^2 (y) \leftrightarrow \psi_4(x_1)\psi_0(x_2), \\
Y^1_2 \sim e^{i\phi} P_2^1 (y) \leftrightarrow \psi_3(x_1)\psi_1(x_2), \\
Y^0_2 \sim  P_2^0 (y) \leftrightarrow \psi_2(x_1)\psi_2(x_2), \\
...,
$$
indeed provide equivalent realizations, in two variables, of the very same matrix of derivatives.
I provide the nontrivial quantities with flat measure,
$\int dx ~\psi_n(x) \psi_k(x)=\delta_{nk}$ (not polynomials!), and $\int dy ~P_l^m(y)~ P_k^m(y)\propto \delta_{lk}$, to compare apples with apples: so they provide an obvious correspondence dictionary, with evident correspondence rules. But I don't think that was the change of variables in the sky the OP was asking for.
