Conditional evaluation: black hole cross section

Question

I'm trying to calculate the collision cross section of black holes, more precisely, I want the critical impact parameter $b_{\text{crit}}$ for a given initial velocity $v$ such that, at this speed, any smaller impact parameter would result in the particle crossing the event horizon. I have code which evaluates each step in the reasoning, but I'm having trouble getting it all to work together.

The idea is:

(1) from the impact parameter and initial speed, numerically integrate the equation of motion to get $r(\phi)$

(2) Get the distance of closest approach $R(b,v)$ as a function of the initial motion parameters.

(3) solve $R(b,v)=1$ varying $b$, the value of $b$ where this equation holds is $b_{\text{crit}}$.

The equation of motion is

$$ u''+u=\frac{1}{2(\gamma bv)^2}+\frac{3}{2}u^2 $$

We have here $u=1/r$. And I'm using units where $c=1$, $2GM=1$, $r_{\text{Schwarzschild}}=1$.

To do this I have the code

Func[b_, v_] := 
  NDSolve[{D[u[p], {p, 2}] + u[p] == 1/(2* (b*v) ^ 2) + 3/2* u[p]^ 2, 
    u[0] == 0, u'[0] == 0}, {u[p]}, {p, 0, 28}];


U[p_, b_, v_] := u[p] /. First[Func[b, v]];

So 1/U(b,v) gives me the motion for a given $b$ and $v$. Then to find the minimum I have

Rmin[bInit_, vInit_] := 
 MinValue[{(1/U[p, bInit, vInit]), 1 > p > 0}, p] 

Wherein lies the problem: In some cases, the particle is actually captured by the black hole, so $r$ will tend to zero and bad things will happen near the singularity. And in some cases the particle gets captured for very large $\phi$ (p in the above code). So I need some condition to set the range of p to look at.

One way to solve this would be to evaluate

NSolve[1/U[p, 100, 0.01] == 1, p]

Where I've chosen some values for $(b,v)$. If NSolve finds a solution to this, then the particle gets captured, so the answer will be a larger b, for fixed v, if NSolve finds no solutions to this, then the particle won't get captured and $b_{\text{crit}}$ is smaller. Maybe some sort of bisection would work but I'm not sure how to do this in Mathematica.

The ultimate objective is to get a function which returns $b_{\text{crit}}(v)$

=======

EDIT: Here is a diagram to make things a bit clearer:

If ever $r$ is less than one, then the particle gets captured. $\Delta\phi$ is the deflection from the straight line path that the particle would take if there was no central mass.

Answer 1

I think it is a fundamental property of your differential equation that there is a bifurcation in behavior at bv=2. Greater than that, you have no capture. Less than that, you have capture.

In playing with the relationship, it doesn't appear that there is a value for the product b v such that for some p you have Min[1/u[p]]==1. Rewriting your code in terms of a single term bv...

 Func[bv_] := NDSolve[{u''[p] + u[p] == 1/(2 (bv)^2) + 3/2 u[p]^2, u[0] == 0, 
                       u'[0] == 0}, {u[p]}, {p, 0, 28}];

Now look at the curves using Manipulate.

Manipulate[Plot[Evaluate[1/u[p] /. (Func[k] // First // First)], {p, 0, 28}, 
                PlotRange -> {{0, 28}, {0, 5}}], 
            {k, 1, 10}]

Here's a plot of 1/u[p] with bv=2

And now with bv=1.9999999

Plot[Evaluate[1/u[p] /. (Func[1.999999] // First // First)], {p, 0, 28}, 
      PlotRange -> {{0, 50}, {0, 5}}]

The value of bv=2 seems to be special, as when set, we can directly solve the differential equation.

DSolve[{u''[p] + u[p] == 1/8 + 3/2 u[p]^2, u[0] == 0, u'[0] == 0}, u[p], p]

 (* u[p] -> 1/2 Tanh[p/(2 Sqrt[2])]^2 *)

In playing around, no luck with other values.