Expansion of Green's Function in Spherical Coordinates in Jackson So I was reading about the expansion of the Green function in Spherical coordinates from Classical Electromagnetism by J.D. Jackson and I'm really confused about a subtle step that he makes to go from equations $(3.117)$ and $(3.118)$ to equation $(3.119)$ (the image is attached at the bottom of the question). Summarising the discussion, since we can expand any function of $(r,\theta, \varphi)$ in terms of the Spherical Harmonics $Y_{lm}(\theta,\varphi)$ and the radial function $U_{lm}(r)$ as -
$$F(r,\theta,\varphi)=\sum_{l=0}^\infty \sum_{m=-l}^{+l} (A_{lm}r^{l} +B_l r^{-(l+1)})Y_{lm}(\theta,\varphi) \tag{1}$$
then I can also write the Green function $G(\textbf{x},\textbf{x}')$ assuming that $\textbf{x}'$ is a fixed point and acts as some parameter. Then, I can expand the Green function too in terms of $Y_{lm}(\theta,\varphi)$ and $U_{lm}(r)$ as -
$$G(\textbf{x},\textbf{x}')=\sum_{l=0}^\infty \sum_{m=-l}^{+l} A_{lm} (r|r',\theta', \varphi') 
 Y_{lm}(\theta,\varphi) \tag{2} $$
where $A_{lm}$'s are the expansion coefficients with the $\textbf{x}'$ dependence wrapped into them  ($r',\theta', \varphi'$ act like parameters in this expansion). Now, Jackson says that using the equation ${\nabla_{r}}^2 G(\textbf{x},\textbf{x}')=-\delta(\textbf{x}-\textbf{x}'$), and expanding both sides in terms of $Y_{lm}$'s and comparing the coefficients,  we get that- $$A_{lm} (r|r',\theta', \varphi') =g_l (r,r'){Y_{lm}}^{*}(\theta ', \varphi ') \tag{3}$$
which I absolutely do not understand. There was no context at all about what $g_l$ is or why we could factorize $A_{lm}$ in this way. In fact, when I tried to do what Jackson says, i.e. expanding both sides of the Green Function equation and equating the coefficients, I instead get a differential equation as follows-
$${\nabla_r}^2 A_{lm}(r|r',\theta', \varphi') - \frac{l(l+1)}{r^2}A_{lm}(r|r',\theta', \varphi')= -\frac{1}{r^2} \delta{(r-r')} {Y_{lm}}^{*} (\theta ', \varphi ') \tag{4}$$
and I have no clue how $(3)$, the separation Jackson quotes comes from the above differential equation $(4)$. I'm really stuck on this part for days and I've looked at various online resources, and the usual thing that is done in most lecture notes is that it is often assumed that you can separate the $A_{lm}$ but it logically does not make very much sense to me. Any sort of help would be really appreciated. Thanks!
Edit: By secavara's suggestion, I expanded the coefficient $A_{lm}(r|r',\theta', \varphi')$ in terms of the conjugate of spherical harmonics as (expanding with respect to $r',\theta',\phi'$ assuming r to be some parameter)-
$$A_{lm}(r|r',\theta',\phi') = \sum_{\hat{l}=0}^\infty \sum_{\hat{m}=-\hat{l}}^\hat{l} H_{\hat{l}\hat{m}}(l,m,r,r') Y^*_{\hat{l}\hat{m}}(\theta',\phi') \tag{5}$$
Subbing this into the differential equation $(4)$ that I got, I get
$$\implies \sum_{\hat{l}=0}^{\infty} \sum_{\hat{m}=-\hat{l}}^{+\hat{l}} \bigg( {\nabla_r}^2 H_{\hat{l}\hat{m}} -\frac{l(l+1)}{r^2} H_{\hat{l}\hat{m}} \bigg) {Y_{\hat{l} \hat{m}}}^{*} (\theta' ,\phi') = -\frac{\delta(r-r')}{r^2} {Y_{lm}}^{*} (\theta', \phi') \tag{6}$$
Now if we compare the coefficients of LHS with RHS, we see that-
$$\bigg( {\nabla_r}^2 H_{\hat{l}\hat{m}} -\frac{l(l+1)}{r^2} H_{\hat{l}\hat{m}} \bigg) =0 $$
for all $\hat{l} \neq l$ or $\hat{m} \neq m $. But this merely suggests that coefficients other than $H_{ lm}$ satisfy another differential equation and not that they are $\textbf 0$. Any further suggestions would be appreciated!
From Jackson Page 120

 A: This is basically a compilation of the discussion I had with secavara and Anonjohn in the comment section so all credits go to them. Thanks, everyone!
So, starting from the differential equation $(4)$ which we get on substituting the expanded form of $G(\textbf{x},\textbf{x}')$ into the Green function equation ${\nabla_{r}}^2 G(\textbf{x},\textbf{x}')=-\delta(\textbf{x}-\textbf{x}')$ as
$${\nabla_r}^2 A_{lm}(r|r',\theta', \varphi') - \frac{l(l+1)}{r^2}A_{lm}(r|r',\theta', \varphi')= -\frac{1}{r^2} \delta{(r-r')} {Y_{lm}}^{*} (\theta ', \varphi ') \tag{1.1}$$
As suggested by secavara, just like we expanded $G(\textbf{x},\textbf{x}')$ in $(r,\theta,\phi)$ coordinates with the $\textbf{x}'$ dependence wrapped into the expansion coefficients $A_{lm}$'s , we can similarly expand $A_{lm}(r|r',\theta',\phi')$ by using the fact that as a function of $\theta′$ and $\phi′$, it can be expanded in terms of spherical harmonics (or their conjugates, for that matter) as follows-
$$A_{lm}(r|r',\theta',\phi')=\sum_{\tilde{l}=0}^\infty \sum_{\tilde{m}=-\tilde{l}}^{\tilde{l}}H_{\tilde{l}\tilde{m}} (l,m,r,r') {Y_{\tilde{l}\tilde{m}}}^*(\theta',\phi') \tag{1.2}$$
Substituting the form of $A_{lm}$ from equation $(1.2)$ into the differential equation in $(1.1)$, the differential equation becomes-
$$\bigg({\nabla_r}^2-\frac{l(l+1)}{r^2}\bigg)\bigg(\sum_{\tilde{l}=0}^\infty \sum_{\tilde{m}=-\tilde{l}}^{\tilde{l}}H_{\tilde{l}\tilde{m}} (l,m,r,r') {Y_{\tilde{l}\tilde{m}}}^*(\theta',\phi') \bigg)=-\frac{\delta(r-r')}{r^2} {Y_{lm}}^{*} (\theta', \phi')$$
$$\implies \sum_{\tilde{l}=0}^{\infty} \sum_{\tilde{m}=-\tilde{l}}^{+\hat{l}} \bigg( {\nabla_r}^2 H_{\tilde{l}\tilde{m}} -\frac{l(l+1)}{r^2} H_{\tilde{l}\tilde{m}} \bigg) {Y_{\tilde{l} \tilde{m}}}^{*} (\theta' ,\phi')=-\frac{\delta(r-r')}{r^2} {Y_{lm}}^{*} (\theta', \phi') \tag{1.3}$$
By comparing the coefficients in LHS and RHS of equation $(1.3)$, we see that for all $\tilde{l}\neq l$ or $\tilde{m} \neq m$, we get the following condition-
$$\bigg({\nabla_r}^2 - \frac{l(l+1)}{r^2}\bigg) H_{\tilde{l} \tilde{m}}=0 \qquad \forall \tilde{l}\neq l, \tilde{m} \neq m \tag{1.4}$$
Now the critical physical input that is required in this problem, as was pointed out by Anonjohn was that the Dirichlet boundary conditions for the Green function require that it vanishes at infinity which in turn implies that $H_{\tilde{l} \tilde{m}} \rightarrow 0$ as $|r|\rightarrow \infty$. Since the differential equation still holds when we analyze its asymptotic behavior at $|r|\rightarrow \infty$, we get that -
$${\nabla_r}^2 H_{\tilde{l} \tilde{m}}=0 \implies \frac{d^2}{dr^2}(H_{\tilde{l}\tilde{m}})=0 \implies H_{\tilde{l}\tilde{m}}=ar+b \qquad \text{(asymptotic behaviour)}\tag{1.5}$$
Thus, requiring that the green's function vanishes at infinity is tantamount to setting $a,b$ to $0$ which implies $H_{\tilde l \tilde m} =0$ for all $\tilde{l}\neq l$ or $\tilde{m} \neq m$. This trivially implies that the expansion of $A_{lm}$, as can be seen from equation $(1.2)$, would look like-
$$\boxed{A_{lm}(r|r',\theta',\phi')= H_{{l}{m}} (r,r') {Y_{\tilde{l}\tilde{m}}}^*(\theta',\phi') = g_l(r,r'){Y_{\tilde{l}\tilde{m}}}^*(\theta',\phi')} \tag{1.6}$$
since all the other $H_{\tilde{l}\tilde{m}}$'s are zero. Hence, in $(1.6)$, we get our desired factorization if we call $H_{lm}(r,r')\equiv g_l(r,r')$ (since it ultimately depends only on $l$ and not $m$ as can be seen from the differential equation). The differential equation satisfied by $g_l$ can be obtained on substituting the form of $(1.6)$ into the differential equation of the Green function to get-
$$\boxed{\frac{1}{r^2} \frac{d^2}{dr^2}(rg_l(r,r'))-\frac{l(l+1)}{r^2}g_l(r,r')=-\frac{\delta(r-r')}{r^2}} \tag{1.7}$$
Edit: I discussed this with some of my friends and a flaw in this derivation is the assumption that the Green function does vanish at infinity. In the general case, there is no mention of the outer boundary being at infinity. The region of interest can have an inner boundary of radius $R_1$ and an outer boundary of radius  $R_2$, and then the asymptotic analysis doesn't really make sense. So, I am still leaving the question as unanswered so that someone else could provide their suggestions on this too.
