Gravitational potential energy of a polytropic sphere using the Lane-Emden equation

I apologise if this is a naive question. I’ve been 2 days in a mental puzzle trying to catch possible mistakes solving for the gravitational potential energy of a star, assuming a polytropic model.

I’ve been advised to follow this solution but I still (before consulting the source) can’t tell what am I doing wrong. I’ll try to explain my approach in the next lines.

First of all the $vec{g}$ field

then

where also $vec{nabla}cdot vec{g}= -4pi G rho(vec{x})$, $vec{g}(vec{x}) = – vec{nabla}phi$ and $phi(vec{x}) = – Gint frac{rho(vec{x}’)}{|vec{x}-vec{x}’|}d^3x’$. We have the Poisson equation

$$vec{nabla}^2 phi(vec{x}) = 4pi G rho(vec{x})$$
in which $vec{nabla}^2 phi =0$ where the mass distribution is 0 (outside the sphere). In this context $W_{AB}$ would be $mint_A^B dphi = m(phi(B) – phi(A))$

In a discrete way

$$vec{g}(vec{x}) = – Gsumlimits_{i=1}^N m_i frac{vec{x} – vec{x_i}}{|vec{x} – vec{x}_i|^3}$$

then

$$W = frac{1}{2} sumlimits_{j=1}^N W_j = -frac{G}{2} sumlimits_{j=1}^N sumlimits_{i=1, i neq j}^N frac{m_i m_j}{|vec{x}_j -vec{x}_i|}$$

And assuming a polytropic model for $rho = rho(r)$, with spherical simmetry, the solution for $W$ follows

$$W = frac{1}{2} int rho(x) phi(x) d^3x = frac{1}{2} 4pi int_0^infty rho(r) phi(r) r^2 dr = frac{1}{2} int_0^infty phi(r) dM(r)$$

and integrating by parts(the whole procedure can be found at pg. 57)

$$phi(R) = – frac{GM}{R} qquad boxed{W = -G int_0^{infty} frac{M(r)}{r} dM(r)}$$

First question: is this valid?

The thing is that a polytropic equation for $rho$ (the density)

$$rho(r) = rho_c theta^n (r)$$
with $theta$ an adimensional quantity, $n$ a parameter from 0 to 5 and $rho_c$ the central density of the modelled star, satisfies

$$rho(r) = frac{dM}{dV} rightarrow dM(r) = rho(r) dV(r)$$
so $W$ becomes

where $-theta^n(r)$ can be integrated using Lane Emden’s equation:

$$frac{1}{xi^2}frac{d}{dr}left(r^2 frac{dtheta}{dxi}right) = – theta^n (r) qquad xi = frac{r}{r_n} qquad r_n^2 = frac{(n+1)P_c}{4pi G rho_c^2}$$

Integrating this and using Boundary conditions I end up with a nice closed expression

$$Phi_{mathrm{sist}}= -2pi G rho_c r_n^2 M = -frac{1}{2} frac{M P_c}{rho_c}(n+1)$$

Is this kind of approach familiar to anyone here? Probably there’s something wrong… anyways I’m interested in learning about these mistakes