Commutator and Factorization of the Eigenfunctions I have stumbled upon the following statement:

Consider an Hamiltonian $H$ that is function of a multitude of operators: $H(\hat{O_1},\hat{O}_2,...,\hat{O}_n)$. If we can show that $H$ commutes with all of these operators
$$[H,\hat{O}_1]=[H,\hat{O}_2]=.....=[H,\hat{O}_n]=0$$
then we can write its eigenfunctions $\psi _H$ in a factorized form, as the product of the eigenfunctions of all the operators $\hat{O}$:
$$\psi_H=\psi_{O_1}\psi_{O_2}.....\psi_{O_n}$$

My question is: How can we prove this statement true?
(Strongly related) Bonus question: Is there a better way to state this theorem? Is there a more general form of it?

To make the situation as clear as possible let me give a concrete example of application of this theorem:

Consider the Hamiltonian:
$$H=\frac{p^2}{2m}+k\left(\vec{p}\cdot \vec{S}\right) \tag{1}$$
we want to find its spectrum and its eigenfunctions. This can seem challenging at first, but we can use the upper mentioned theorem: Notice in fact that
$$[H,p_i]=0 \ \ ; \ \ i=1,2,3$$
this means that we can write the eigenfunctions of $H$ as:
$$\psi=\psi_p\chi$$
where $\psi _p$ are the eigenfunctions of the operator $\vec{p}$ and $\chi$ are the eigenfunctions of the spin $\vec{S}$, so:
$$\psi=\psi_p\chi=\frac{1}{(2\pi\hbar)^{3/2}}\exp\left[\frac{i}{\hbar}\vec{p}\cdot \vec{x}\right]\chi \tag{2}$$
and from here the derivation of the eigenfunctions of $H$ is much simpler.


Edit in response to the comments: I am afraid I haven't expressed myself well; I will try to clarify even more:
I know that there is a strong link between commutation and factorizability of the eigenfunctions but I don't know what this link is precisely and I also of course don't know how to prove it. This is my problem.
 A: I think your theorem is a thoroughgoing misconception of the p.d.e. factorization, of, e.g.,  a spherically symmetric system. Nondimensionalizing all silly constants by absorbing them in the relevant units, you have something like $-\Delta +V(r)  -E=0$.
Its eigenvectors are not the product of the eigenvectors of all symmetry generators (operators commuting with the hamiltonian), here, among others, the three $\vec L$s. Instead, recall how this equation's variables separate in plain p.d.e. theory (cf. separation of variables ):
$$
0=-\Delta +V(r) -E\\ = -\frac{1}{r^2}\partial_r ~r^2\partial_r + \frac{1}{r^2} L^2(\theta,\phi) +V(r) -E ~~~\leadsto \\
  L^2 (\theta,\phi) = \partial_r ~r^2\partial_r  -r^2 V(r) +r^2 E .
$$
This is plain separation of variables: each side of the equation involves different variables, so its eigenvector structure is disjoint. The eigenvectors of the l.h.s. are the spherical harmonics, $Y_{lm}(\theta, \phi)$, with eigenvalues $l(l+1)$ and the eigenvectors of the r.h.side must be
functions of just r, but with the same eigenvalues, i.e.
$$
  -\frac{1}{r^2}\partial_r ~r^2\partial_r + \frac{l(l+1)  }{r^2}  +V(r) =E, 
$$
now written in a more familiar form, and with E to be determined, the new eigenvalue, for eigenfunctions $R_{nl}(r)$.
The radial eigenfunctions are inert under $L(\theta,\phi)^2$, but still contain its eigenvalues; and the spherical harmonics do not depend on r, but, of course, are not eigenfunctions of $L_x$, $L_y$, only of $L^2$ and $L_z$, which commute with the hamiltonian as well, so they are good symmetry charges for it.
Now, in other coordinate systems, and for special potentials, like the Hydrogen's you may be more efficient (cf. Pauli's original SO(4) solution of the problem; could do worse than studying this one.), but factorization of p.d.e.s is usually guided by symmetry, as you saw above. You should best consider selected eigenfunctions of the symmetry operators, and utilize those that entangle least with the rest of them.
Finally, in the trivial case where the symmetries commute among themselves, then, of course, the Hilbert space itself factorizes into a tensor product; whose tensor factors, and thus wave function, factor, and are exclusively operated upon by the corresponding eigen-operator, oblivious of the other tensor factors corresponding to the other operators. If this trivial case is the case your instructor discussed, it's hardly salutary to formalize  it so impossibly abstractly.


*

*Response to comment on suppositious example (1).

Let's write it in nondimensionalized units,
$$
0=-E+\vec p^2 +2k \vec p\cdot \vec S .
$$
Without loss of generality, for the purposes of illustrating the problem, take $\vec S= \vec \sigma /2$. Since the problem is manifestly spherically symmetric, we can always rotate the spin to the 3rd (z) direction, without affecting its eigenvalues!
$$
0=-E+\vec p^2 +k   p_z  \sigma_3~~.    
$$
The problem has separated to three decoupled pieces,
$$
0=-E+p_x^2+  p_y^2+(  p_z^2 +kp_z\sigma_3 ) .
$$
The first two pieces are scalar, but the third is a 2×2 matrix, so, acting on the space of 2-spinors.
The eigenvalues  $E=\vec p^2 \pm kp_z$ are the sum of the eigenvalues of each piece on the right, for the eigenfunctions constant(x), constant'(y) and constant times
$$
  \begin{bmatrix} 1\\0\end{bmatrix} ; \qquad  \begin{bmatrix} 0\\1\end{bmatrix}, 
$$
respectively. Note the last eigenfunction hinges on the eigenspinors of $S_3$ only: you couldn't possibly diagonalize all three $S_x,S_y,S_z$ simultaneously. Actually, the first two have completely dropped out of the problem.
In momentum space, which you converted to a meretricious laconism by Fourier-transforming, there is hardly any insight into your suppositious theorem, and you properly invited me to not focus on it. Your hamiltonian is a 2×2 spin matrix and obviously its eigenvectors are 2-spinors.
A: I am going to try to answer my own question here, since I think I now understand what I was missing before.
We want to understand the link between commutation relation and factorizability of the eigenfunction. This link is all about the definition of separable Hamiltonians and about a proper interpretation of the meaning of commutation relations:
Think about an Hamiltonian of the form:
$$H=H_1+H_2 \tag{1}$$
we could have chosen to work with an Hamiltonian that is the sum of multiple terms, not just two, but the generalization from two terms to many terms is really easy so for now, to make our life easier, let's stick to an Hamiltonian of the form (1). This Hamiltonian is called separable if the following commutation relation holds true:
$$[H_1,H_2]=0$$
This is simply a definition. Now: the crucial fact is that for a separable Hamiltonian the following is true:

*

*Its eigenvalues $E$ are the sum of the eigenvalues of the two Hamiltonians $H_1,H_2$:
$$E=E_1+E_2$$

*Its eigenfunctions $\psi$ are the product of the eigenfunctions of $H_1,H_2$:
$$\psi=\psi _1 \psi _2$$
This is the relation between commutation and factorizability of the eigenfunctions! But we are of course not done, we have to prove that all this is true.
Let's start with the proof: the fact that $[H_1,H_2]=0$ means that we can find a common base of eigenvectors for $H_1,H_2$, this is a really famous theorem. So we have:
$$H_1|E_1,E_2\rangle=E_1|E_1,E_2\rangle$$
$$H_2|E_1,E_2\rangle=E_2|E_1,E_2\rangle$$
but this immediatly means that:
$$H|E_1,E_2\rangle=(H_1+H_2)|E_1,E_2\rangle=H_1|E_1,E_2\rangle+H_2|E_1,E_2\rangle=E_1|E_1,E_2\rangle+E_2|E_1,E_2\rangle$$
so:
$$H|E_1,E_2\rangle=(E_1+E_2)|E_1,E_2\rangle$$
And just like that we have proven point one! We now need to prove point two, the point about the factorization of the eigenfunction that started all this trouble. To prove this point we have to project our equations on a base, to switch from eigenvectors to eigenfunctions. The foundamental observation here is the following:

*

*If two distinct operators (for example the momentum $p_x$ and the spin $S_x$) commute, meaning that:
$$[p_x,S_x]=0$$
then we can think of this operators as acting on separate Hilbert spaces!
This is the key observation here, because it makes us understand that if we want to project we can't just write:
$$\langle x|H|E_1,E_2\rangle=(E_1+E_2)\langle x|E_1,E_2\rangle$$
this is wrong! We instead have to write:
$$\langle x_1,x_2|H|E_1,E_2\rangle=(E_1+E_2)\langle x_1,x_2|E_1,E_2\rangle$$
This is correct! But this compact notation $|a,b\rangle$ is not doing any favor to us in this context, instead is actively harming our understanding; better to switch to the proper notation:
$$|a,b\rangle = |a\rangle \otimes |b\rangle$$
using the fundamental fact (sometimes seen as a postulate) that the total Hilbert space is the tensor product of the Hilbert spaces. So we have:
$$(\langle x_1| \otimes \langle x_2|)H(|E_1\rangle \otimes |E_2\rangle)=(E_1+E_2)(\langle  x_1| \otimes \langle x_2|)(|E_1\rangle \otimes |E_2\rangle)$$
This can be rewritten as follows:
$$H'(\langle x_1| \otimes \langle x_2|)(|E_1\rangle \otimes |E_2\rangle)=(E_1+E_2)(\langle  x_1| \otimes \langle x_2|)(|E_1\rangle \otimes |E_2\rangle) \tag{2}$$
You can see that $(\langle x_1| \otimes \langle x_2|)$ jumped from the left to the right of the hamiltonian operator, this is due to the fact that we projected the hamiltonian operator on the space of $x_1,x_2$, meaning that $H'$ is now no longer an operator that acts on kets but insted an operator that acts on functions! (This procedure is really common in QM and $H'$ usually is just called $H$, even if it acts no longer on kets, this can be quite confusing if not explicitly stated) Ok, we have (2); now what? Now we remember the scalar product rule for products in different Hilber spaces! (you can find all about it here) Applying the definition of inner product (scalar product) to our equation (2) yields:
$$H\langle x_1|E_1\rangle \langle x_2|E_2\rangle=(E_1+E_2)\langle x_1|E_1\rangle \langle x_2|E_2\rangle \tag{3}$$
but since, of course, the eigenfunctions of $H_1,H_2$ are defined as (and this of course comes from the definition of wavefunction itself):
$$\psi _1 = \langle x_1| E_1 \rangle$$
$$\psi _2 = \langle x_2 | E_2 \rangle$$
follows that from (3) we can immediatly derive that (rewrite, really):
$$H\psi_1\psi_2=(E_1+E_2)\psi _1\psi_2$$
And this proves our precious point 2! We are done!
A: For the wave function to be a product, it is necesary and sufficient to deal with separated variables in the Schroedinger equation.
If the operators $\hat{O}_n$ commute with $\hat{H}$, then $\Psi_E$ may be an eigenstate for each operator in question.
