Trying to understand the relationship between Hydrogen atom, spherical harmonics and central field force in quantum mechanics I have a problem understanding three arguments in quantum mechanic:

*

*When we talk about a particle in a central field we have this kind of Hamiltonian:
$$H=\frac{p^2}{2m}+V(r)$$
if we use spherical coordinates we can write the Laplacian operator as:
$$\nabla^2=\frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r}+\frac{1}{r^2}(\frac{\partial^2}{\partial \theta^2}+\frac{cos(\theta)}{sin(\theta)}\frac{\partial}{\partial \theta}+\frac{1}{sin^2(\theta)}\frac{\partial^2}{\partial \phi^2})$$
and so $p^2$ in Schordinger's representation as:

$$p^2=-\hbar^2(\frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r})+\frac{L^2}{r^2}$$
so that the entire Schrodinger's equation becomes:
$$(\frac{\partial^2}{\partial r^2}+\frac{2}{r}\frac{\partial}{\partial r}-\frac{L^2}{\hbar^2r^2}-\frac{2mV(r)}{\hbar^2}+\frac{2mE}{\hbar^2})\psi(r,\theta,\phi)=0 \  (1)$$
Now the problem: my professore said that "

since a central Hamiltonian commutes with the operators $L^2$ and
$L_z$ (I know why), we write the solution of equation (1) in the factorised form:

$$\psi(r,\theta,\phi)=F(\theta,\phi)R(r)$$
First question: Why, from the fact that a central Hamiltonian commutes with the operators $L^2$ and $L_z$, we come to write the eigenfunction in that "separated" form? I think the answer is related to the fact that if 2 operators commute with each other they have a complete set of simultaneous eigenfunctions, but I can't understand why this lead to the "separate" form of the Hamiltonian eigenfuntion.


*The other topic that causes me the same confusion is the following one:

talking about the Hydrogen atom we have $$H=\frac{p_1^2}{2m_1}+\frac{p_2^2}{2m_2}+V(|q_1-q_2|)$$ which becomes, after a change of variables $$H=\frac{P^2}{2M}+\frac{p^2}{2\mu}+V(|q|) \ (2)$$
where $M=m_1+m_2$ and $\frac{1}{\mu}=\frac{1}{m_1}+\frac{1}{m_2}$
Now in my notes I've written that

the equation (2) is a separable variables equation, namely the sum of two terms $H_1$ which only depends from the center of mass variables and $H_2$ which only depends from relative variables, so  $H_1$ and $H_2$ commute each other. Since the total Hamiltonian is separable we can search the total eigenfunction in the form $\psi(X,x)=\phi(X)\psi(x)$  (3) where $X$ refers to center of mass variable and $x$ to relative variables.

I can now write equation $(2)$ in Schrodinger's representation:
$$(\frac{-\hbar^2\nabla^2_X}{2M}-\frac{-\hbar^2\nabla^2_x}{2\mu}+V(|X|)\psi(X,x)=E\psi(X,x)) \ (4)$$
then I replace equation $(3)$ into equation $(4)$ obtaining:
$$(-\frac{\hbar^2}{2M}\nabla^2_X\phi(X))\psi(x)+((\frac{-\hbar^2}{2\mu}+V(|x|))\psi(x))\phi(X)=E\phi(X)\psi(x)$$
I can now divide by $\phi(X)\psi(x)$:
$$-\frac{\hbar^2}{2M}\frac{\nabla^2_X\phi(X)}{\phi(X)}-\frac{\hbar^2}{2\mu}\frac{\nabla^2_x\psi(x)}{\psi(x)}+V(|x|)=E$$
which leads to two separated eigenvalues equations:
$$-\frac{\hbar^2}{2M}\nabla^2_X\phi(X)=E_1\phi(X)$$
$$-\frac{\hbar^2}{2\mu}\nabla^2_x\psi(x)+V(|x|)\psi(x)=E_2\psi(x)$$
which I can solve separately.
Here the question is:
Why the fact that $H=H_1+H_2$ and $H_1$ commutes with $H_2$ brings me to search the total eigenfunction in the form $\psi(X,x)=\phi(X)\psi(x)$?


*In a website (sorry, I can't find it anymore) I've read this when talking about spherical harmonics:


"if we want to determine simultaneous eigenfuntions of $L^2$ and $L_z$
we have to solve this system:
$$L^2|l,m\rangle=\hbar^2l(l+1)|l,m\rangle$$ $$L_z|l,m\rangle=\hbar m|l,m\rangle$$
in Schrodinger's representation:
$$-(\frac{1}{sin(\theta)}\frac{\partial}{\partial{\theta}}(sin(\theta)\frac{\partial}{\partial{\theta}})+\frac{1}{sin^2(\theta)}\frac{\partial^2}{\partial{\phi^2}})\psi_{lm}(r,\theta,\phi)=l(l+1)\psi_{lm}(r,\theta,\phi)$$
$$-i\frac{\partial}{\partial{\phi}}\psi_{lm}(r,\theta,\phi)=m\psi_{lm}(r,\theta,\phi)$$
we can note that operators $L^2$ and $L_z$ don't contain variable $r$
so the $\psi_{lm}(r,\theta,\phi)$ must be of the form
$\psi_{lm}(r,\theta,\phi)=f(r)Y_{lm}(\theta,\phi)$.

Third question: Why if operators $L^2$ and $L_z$ don't contain variable $r$ the $\psi_{lm}(r,\theta,\phi)$ must be of the form
$\psi_{lm}(r,\theta,\phi)=f(r)Y_{lm}(\theta,\phi)$?
The final question about all this is:
are the three problems I've written related in a certain way? because all three propose an eigenfunction written as a product of eigenfunctions but in the first case this is due to the fact that $H$, $L^2$, $L_z$ commute with each other. In the second case this is due to the fact that $H$ can be written as a sum of Hamiltonians that commute with each other and we speak of separation of variables. In the third case it is due to the fact that the variable $r$ does not appear in the operators $L^2$ and $L_z$. Do all three mean the same thing or are they three different things?
 A: Note: My approach lacks rigour and I'm rusty. I'm sure there are lots of holes, overly-complicated parts, and perhaps even just outright mistakes. However, I think the gist of the argument motivates why two operators which commute leads one to consider solutions of the form you mention.
Paraphrased: Why does does the fact that $\hat{A}$ and $\hat{B}$ commute in the Hamiltonian $\hat{H}=\hat{A} + \hat{B}$ mean that we can write the solution as factors of the solutions of $\hat{A}$ and $\hat{B}$?
Essentially, the fact that the Hamiltonian's parts, $\hat{A}$ and $\hat{B}$, commute means we can separate the general solution (of the form $\exp(-i\hat{H})$ into two exponential operators.
My approach lacks rigour, but I think it shows the idea:
For the equation
$\newcommand{\Ket}[1]{\left|#1\right>}$
$\newcommand{\Bra}[1]{\left<#1\right|}$
$\newcommand{\Braket}[2]{\left<{#1}|{#2}\right>}$
\begin{equation}
\hat{H}\Ket{\psi} = i\frac{\partial}{\partial t}\Ket{\psi} \tag{1}\label{SE}
\end{equation}
the general solution is known to be
\begin{align}
\Ket{\psi(t)} &= \exp\left(-it\hat{H} \right) \Ket{\psi(0)} \\
&= \exp\left(-it(\hat{A}+\hat{B}) \right) \Ket{\psi(0)}
\tag{2}
\end{align}
You can see this through substitution or look it up in any QM textbook or course notes.
The fact that $[A,B]=0$ means we can split the exponential operator, giving
\begin{align}
\Ket{\psi(t)} &= \exp\left(-it\hat{H} \right) \Ket{\psi(0)} \\
&= \exp\left(-it\hat{A} \right)\exp\left(-it\hat{B} \right) \Ket{\psi(0)}
\tag{3}
\end{align}
Let's try to write this in position representation.
First, let's recognise that we're dealing with more than one spatial coordinate and each operator is acting in mutually exclusive vector spaces, so that $\mathcal{H} = \mathcal{H}_A \otimes \mathcal{H}_B$. Alternatively, in position representation, we can see that $\hat{A}$ has an $x$-dependence and $\hat{B}$ has a $y$-dependence (and can be seen to commute by observation).
Also, the time-independent solutions for $\hat{A}$ and $\hat{B}$ separately in position representation are
\begin{align}
f(x) &= \Braket{x}{f} = \sum_i\lambda_i\Braket{x}{a_i} \\
g(y) &= \Braket{y}{g} = \sum_i\mu_i\Braket{x}{b_i} 
\tag{4}
\end{align}
where the eigenfunctions are
\begin{align}
\hat{A}\Ket{a_i} &= a_i \Ket{a_i} \\
\hat{B}\Ket{b_i} &= b_i \Ket{b_i}
\tag{5}
\end{align}
Then using $\Braket{x,y}{\psi} = \psi(x,y)$ with $\Bra{x,y} = \Bra{x}\otimes\Bra{y}$, we get
\begin{align} 
\Braket{x,y}{\psi} =& \left[\Bra{x}\exp\left(-it\hat{A} \right) \otimes \Bra{y}\exp\left(-it\hat{B} \right)\right] \Ket{\psi(0)}
\tag{6}
\end{align}
which we can write using the expansion of the identity as
\begin{align} 
\Braket{x,y}{\psi} =& \left[\sum_n\Bra{x}\exp\left(-it\hat{A} \right) \Ket{a_n}\Bra{a_n} \right. \\
&\otimes \left. \sum_m\Bra{y}\exp\left(-it\hat{B} \right)\Ket{b_m}\Bra{b_m}\right]\psi(0)
\tag{7}
\end{align}
where the operators can now act on their respective eigenstates, giving
\begin{align} 
\Braket{x,y}{\psi(t)} =& \sum_n\Braket{x}{a_n}\exp\left(-ita_n \right)\sum_m\Braket{y}{b_m}\exp\left(-itb_m \right) \\ &\int\int dx^\prime dy^\prime\Braket{a_n}{x^\prime}\Braket{b_m}{y^\prime} \Braket{x^\prime,y^\prime}{\psi(0)}
\tag{8}\label{FINALEQ}
\end{align}
where I've simplified the expression using $$\left(\Bra{x^\prime}\otimes\Bra{y^\prime}\right)\Ket{\psi(0)} = \Braket{x^\prime,y^\prime}{\psi(0)} \tag{9}$$
Note that $\Braket{x}{a_n}$ are the eigenfunctions of $\hat{A}$ and $\Braket{y}{a_m}$ are the eigenfunctions of $\hat{A}$ in position representation.
So, we can see that $\ref{FINALEQ}$ is simply a linear combination of factors of eigenstates of the individual operators $\hat{A}$ and $\hat{B}$:
\begin{align}
\Braket{x,y}{\psi(t)} =& \sum_{m,n}c_{m,n}f_n(x,t)g_m(y,t) \tag{10}
\end{align}
This is a general solution.
Ignoring the time-dependence, we can see that just one of those solutions (not constrained by any initial conditions) is $f_n(x,t)g_m(y,t) = \Braket{x}{a}\Braket{y}{b}$ i.e. the multiplication of the eigenstates of the operators separately. Generally, we can consider a linear combination of such and it's still a solution, and in this case, we've worked backwards from that.
I think this is also shows (maybe a bit strong given the lack of rigour... "suggests"?) that the independence of the variables leads us to consider states in different subspaces of $\mathcal{H}$, and hence operators in those different spaces commute, and therefore the general solution in position representation is formed of products of the separate solutions to those subspaces.
It doesn't matter that we're talking about $r$, $\phi$, and $\theta$, or something simpler like $x$ and $y$ - they're just specific representations of the fundamental states.  Therefore, we can propose solutions of the form $\psi(x,y) = f(x)g(y)$.
Now, knowing this, whenever we see a Hamiltonian where $H=H_1 + H_2$, and the parts commute, we can jump straight to the conclusion and state $\psi(x,y) = f(x)g(y)$, where $f(x)$ and $g(y)$ are the solutions, and then immediately go about calculating $f(x)$ and $g(y)$.
