Finding the radius of a neutron star that allows all points on its surface to be seen

Question

At first glance, this fascinating question may seem better placed on Physics Stack Exchange. But, since I am only questioning the mathematics of the solution I decided it was more appropriate to place it here.

How small should a neutron star be so that all points on its surface are visible?

Hints: use the form of the photon orbit equation. Consider a trajectory where $du/dφ = 0$ for $u = 1/R_n$ and $φ = 0$ where $R_n > R$ is the ‘radius’ of the neutron star. Determine numerically, what the ratio $R_n/R$ should be so that $u = 0$ for $φ = π$.

My interpretation of this question is to find the ratio of neutron star radius, $R_n$ to that of the Schwarzschild radius, $R$ given that all points on the neutron star are visible. Taking the most extreme point for photons emanating from a point directly behind the neutron star, since we can't see through the star this point would otherwise be invisible. However, if the curvature of spacetime is so large it will deflect light coming from the back of the star, then it is possible to see this point. If this point can be seen, then by that logic all points of the neutron stars' surface will be visible too.

Above is my attempt to draw the situation. The outer blue ring of the neutron star has a radius of $R_n$ and the inner (red) ring is characterized by the Schwarzschild radius, $R$.

Of course, the situation I drew above is highly contrived as we all know that the light rays emanating from the point at the back of the star will keep bending as there is no reason why they should end up in a perfect straight line.

---

Here is the author's solution which I just cannot understand:

$$\left(\frac{du}{d \varphi}\right)^2+u^2-Ru^3={R_n}^{-2}-R{R_n}^{-3}\tag{1}$$ so that $u=1/R_n$ when $du/d\varphi=0$. Separating the variables gives $$\int_{u=0}^{1/R_n}\frac{du}{\sqrt{{R_n}^{-2}-R{R_n}^{-3}-u^2+Ru^3}}=\int_{\varphi = 0}^{\pi}d\varphi=\pi\tag{2}$$ So that all the surface is visible it is assumed that light that grazes the surface at $φ=0$ escapes to $r=∞$ or $u=0$ at $φ=π$. Using the change of variable $p=Ru$ $$\int_{p=0}^{R/R_n}\frac{dp}{\sqrt{{(R/R_n)}^2-{(R/R_n)}^3 -p^2+p^3}}=\pi\tag{3}$$ To determine $x=R/R_n$ we need to solve $$\int_{p=0}^{R/R_n}\frac{dp}{\sqrt{{x}^2-x^3-p^2+p^3}}=\pi\tag{4}$$ Solve this numerically. The solution is $x ≈ 0.568$ giving the ratio $R_n/R = 1/x = 1.76$. That is if the radius of a neutron star is less than $1.76$ times its Schwarzschild radius the whole surface is visible due to the deflection of light. Current estimates on neutron star sizes give a ratio larger than $1.76$ suggesting that the entire surface of real neutron stars is not visible.

---

I have 2 questions regarding this solution:

1. Below equation $(1)$ it is written that "$u=1/R_n$ when $du/d\varphi=0$". But since $du/d\varphi=0$, a simple rearrangement of $(1)$ should mean that
   $$\frac{du}{d \varphi}=\pm\sqrt{{R_n}^{-2}-R{R_n}^{-3}-u^2+Ru^3}=0\tag{5}$$ If so then the denominator in the integrand of $(2)$ will diverge, so how can it ever be equal to $\pi$? Moreover, with $du/d\varphi =0$ then $\left(du/d\varphi\right)^2=0$ and by my logic there is no differential equation to solve. Put another way, since $u=1/R_n$ this implies that $x=p$, then substituting $x=p$ into equation $(4)$ leads to the same fate. Why was the negative root in $(5)$ discarded?
2. Just below equation $(4)$ it says "Solve this numerically". I really don't know how to do this and I'm not much of a programmer, but all I want to do is verify that this numerical integration results in $x\approx 0.568$. Could someone please show or explain to me how this is done?

---

Remark:

One reason why I am so sceptical about the authors' solution is while typesetting the solution I have already identified 2 typos which I show below underlined in red:

As can be seen above, a square root is missing in the denominator of the second equation and the ratio should strictly be written as '$1.76$'.

Answer 1

There is not many problems to compute $$I=\int\frac{dp}{\sqrt{{x}^2-x^3-p^2+p^3}}=\int \frac{dp}{\sqrt{(p-r_0)(p-r_1)(p-r_2)}}$$ $$I=-\frac{2}{\sqrt{r_1-r_0}}\,F\left(\sin ^{-1}\left(\frac{\sqrt{r_1-r_0}}{\sqrt{p-r_0}}\right)|\frac{r_0-r_2}{r_0-r_1}\right)$$ $$J=\int_0^x\frac{dp}{\sqrt{{x}^2-x^3-p^2+p^3}}=$$ $$-\frac{2}{\sqrt{r_1-r_0}}\Bigg[F\left(\sin ^{-1}\left(\frac{\sqrt{r_1-r_0}}{\sqrt{x-r_0}}\right)|\frac{r_0-r_2}{r_0-r_1}\right) -F\left(\sin ^{-1}\left(\frac{\sqrt{r_1-r_0}}{\sqrt{-r_0}}\right)|\frac{r_0-r_2}{r_0-r_1}\right) \Bigg]$$

The cubic equation $$p^3-p^2+(x^2-x^3)=0$$ has three real roots since its discriminant $$\Delta=(3 x-2)^2 (1-x) x^2 (3 x+1)$$ is always positive for $0\leq x \leq 1$. So, we have three real roots $$r_k=\frac 1 3 \left(1+2 \cos \left(\frac{2 \pi k}{3}-\frac{1}{3} \cos ^{-1}\left(\frac{27 x^3-27 x^2+2}{2} \right)\right)\right)$$ The $\cos ^{-1}\left(.\right)$ severely restrict the domain which is good for the numerical solver.

For any positive rhs larger than $\frac \pi 2$, the solution is $0\leq x \lt \frac 23$.

Pushing further the analysis,a good approximation is $$\color{red}{x\sim\frac{1}{27} \left(14 \sqrt{7} \cos \left(\frac{1}{3} \cos ^{-1}\left(\frac{271}{2401 \sqrt{7}}\right)\right)-17\right)}=0.5680821\cdots$$

Recomputing the integral, the lhs is $3.1415933$

Edit

After @Lutz Lehmann's comment $$p^3-p^2+(x^2-x^3)=(p-x) \Bigg[\left(p-\frac{1-x}{2}\right)^2-\frac{(1-x) (1+3 x)}{4} \Bigg]$$ simplifies the formal calculations using $$a=\frac{1-x}{2}\qquad \text{and} \qquad b=\frac 12\sqrt{(1-x) (1+3 x)}$$

Answer 2

Interesting. For a massless particle in the Schwarzchild metric we have$^\dagger$

$$\tag{1} \left( \frac{du}{d \phi} \right)^2=L^{-2}-u^2+R_su^3 $$

Where $u=1/r$, I have replaced the star's mass $M$ using the definition of the Schwarzchild radius $R_s=2GM/c^2$, and $L=r^2 \dot{\phi}$ is a constant along the trajectory. Evaluate (1) at the start of the trajectory where $du/d\phi=0$, with $u=1/r$ and $r=R_n$ is the initial radius. This yields the constant

$$\tag{2} L^{-2}=R_n^{-2}-R_s R_n^{-3} $$

Rewrite (1) so that

$$\tag{3} \frac{d\phi}{du}=\pm\left( L^{-2}-u^2+R_s u^3 \right)^{-1/2} $$

Integrate (3) from $\phi=\pi$, $u=1/R_n$ (this is the start of the trajectory) to $\phi=0$, $u=0$ (that is, $r \to \infty$).

$$\tag{4} 0-\pi=0-\pm \int\limits_{1/R_n}^0 \frac{du}{\sqrt{L^{-2}-u^2+R_su^3}} $$

The integrand is positive, so we must choose the $-$ sign (on reversing the integral bounds)

$$\tag{5} \pi=\int\limits_0^{1/R_n} \frac{du}{\sqrt{L^{-2}-u^2+R_su^3}}\stackrel{(2)}{=}\int\limits_0^{1/R_n} \frac{du}{\sqrt{R_n^{-2}-R_s R_n^{-3}-u^2+R_su^3}} $$

Before doing anything numerically we should eliminate as many parameters as possible. The dimensionless variables $p=R_s u$ and $p_0=R_s/R_n$ are convenient. From (5) we have

$$\tag{6} \pi=\int\limits_0^{p_0} \frac{dp}{\sqrt{p_0^2-p_0^3-p^2+p^3}}:=I(p_0) $$

The integral (6) exists (for $p_0$ less than some maximum value that keeps the root positive) because the singularity of the integrand at the upper limit is integrable, the integrand behaves like $1/\sqrt{p_0-p}$ as $p\to p_0$. The simplest way to proceed is just evaluate $I(p_0)$ numerically for a range of $p_0$ values. Here is a plot using Mathematica's built in numerical integration:

We see the graphs intersect at about $p_0=0.57$. To proceed further, you are (numerically) finding the roots of the function $f(p_0)=\pi-I(p_0)$, where $I$ itself must be numerically evaluated. I find $p_0=0.568$ using Mathematica's built in FindRoot.

Aside: A back-of-the-envelope estimation to the solution of (6). The main contribution to $I$ will be near the upper limit, so expand

$$ p_0^2-p_0^3-p^2+p^3 \sim (2p_0-3p_0^2)(p_0-p) \qquad,\qquad p\to p_0 $$

Then the estimate of the integral is elementary, and we have

$$ \pi \stackrel{\text{approx.}}{=} \frac{1}{\sqrt{2p_0-3p_0^2}}\int\limits_0^{p_0} \frac{dp}{\sqrt{p_0-p}}=\frac{2\sqrt{p_0}}{\sqrt{2p_0-3p_0^2}} $$

This is an algebraic equation for $p_0$, which has the nonzero solution

$$ p_0=\frac{2}{3}\left(1-\frac{2}{\pi^2}\right)\approx 0.53 \\ \therefore p_0^{-1} \approx 1.88 $$

Which isn't too bad an estimate for the minimal amount of work required.

$\dagger$ See eg. Misner and Thorne Gravitation, eq 25.55 or B. Schutz A first course in GR eq 11.49, or probably any other GR book.

Answer 3

If you return to the second order equation $$ u''(φ)+u(φ)-\frac32Ru(φ)^2=0, $$ or with $p=Ru$ $$ p''(φ)+p(φ)-\frac32p(φ)^2=0, $$ the task amounts to solving the boundary value problem with $p'(0)=0$ and $p(\pi)=0$. The desired result is the value of $p(0)$.

As one wants to avoid the zero solution, add a constant component $T$ with boundary condition $T·p(0)=1$ (sometimes called the Rabinovich trick).

Now insert into the BVP solver of your choice, here python with scipy's solve_bvp

def sys(phi,p,T): return p[1], -p[0]+1.5*p[0]**2

def bc(p0,pp,T): return p0[1], pp[0], p0[0]*T-1

phi = np.linspace(0,np.pi,10)
p = [ 1-(phi/np.pi)**2, -2*phi/np.pi**2 ]

res = solve_bvp(sys, bc, x=phi, y=p, p=[1], tol=1e-10)
print(res.message, f"x={res.y[0,0]:15.12f}, 1/x={res.p[0]:15.12f}") 

which returns

The algorithm converged to the desired accuracy. x= 0.568082086989, 1/x= 1.760308981576

---

Alternatively in this approach, it is not necessary to first eliminate R and then introduce an auxilary T. One can also add $u(0)=1$ to the boundary conditions, thus also preventing the zero condition, and use R as the third constant state component. Then in the solution $x=R$.

def sys(phi,p,R): return p[1], -p[0]+1.5*R*p[0]**2

def bc(p0,pp,R): return p0[0]-1, p0[1], pp[0]

phi = np.linspace(0,np.pi,10)
p = [ 1-(phi/np.pi)**2, -2*phi/np.pi**2 ]

res = solve_bvp(sys, bc, x=phi, y=p, p=[1], tol=1e-8)
print(res.message, f"x={res.p[0]:15.12f}, 1/x={1/res.p[0]:15.12f}") 

which returns

The algorithm converged to the desired accuracy. x= 0.568082086988, 1/x= 1.760308981580