For two coordinate frames $O'$ and $O''$ both offset along the $z$-axis by $\pm R$ respectively, with corresponding offset spherical coordinates $r'$, $\theta'$, $r''$ and $\theta''$, and with prolate spheroidal coordinates $\xi=\frac{r''+r'}{2R}$ and $\eta=\frac{r''-r'}{2R}$

I have checked that the solid irregular prolate spheroidal harmonics $Q_n(\xi)P_n(\eta)$ can expressed as $$ Q_n(\xi)P_n(\eta)=\frac{1}{2}\sum_{k=0}^n\frac{(n+k)!}{k!^2(n-k)!}\left[ -\left(\frac{r'}{R}\right)^kQ_k(\cos\theta') + (-1)^{n+k}\left(\frac{r''}{R}\right)^kQ_k(\cos\theta'')\right] $$ where $Q_n$ and $P_n$ are the Legendre functions. The functions on the right hand side are individually solutions of Laplace's equation.

How can I go about proving the above expansion?

I have been banging my head on the table for hours... I tried using the integral form : $$ Q_n(\xi)P_n(\eta)= \frac{R}{2}\int_{-1}^1\frac{P_n(u)\mathrm{d} u}{\sqrt{\rho^2+(z-Ru)^2}} $$ ($\rho=r\sin\theta$) and tried applying differentiation/integration along $z$ to step up $n$ by one to hopefully use induction, but with no success.

I do have a proof that

$$ P_n(\xi)P_n(\eta)=\frac{1}{2}\sum_{k=0}^n\frac{(n+k)!}{k!^2(n-k)!}{k}\left[\left(\frac{r'}{R}\right)^kP_k(\cos\theta')+(-1)^{n+k}\left(\frac{r''}{R}\right)^kP_k(\cos\theta'')\right] $$ which looks suspiciously similar to the first equation.

Any help is appreciated.