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:

enter image description here

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.

  • $\begingroup$ The equation contains only the product b*v. Why investigate separately the effects of b and v? $\endgroup$ Commented Mar 18, 2019 at 21:55
  • $\begingroup$ I'm looking for the collision cross section as seen by a particle which is at speed $v$ at infinity, i.e. $\pi b_{\text{crit}}^2$, this only depends on $b$. $\endgroup$
    – tbfr416
    Commented Mar 18, 2019 at 22:10
  • 1
    $\begingroup$ What is the condition of capture? $\endgroup$ Commented Mar 18, 2019 at 23:21
  • 1
    $\begingroup$ @tbfr416 We need to check the equation. How is it received? $\endgroup$ Commented Mar 19, 2019 at 12:29
  • 1
    $\begingroup$ This equation 9.33 is for particles that start moving from a state of rest when r->Infinity. Therefore, their energy $E=m_0c^2$. In the scattering problem should be $E>m_0c^2$. $\endgroup$ Commented Mar 19, 2019 at 14:02

1 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

enter image description here

And now with bv=1.9999999

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

enter image description here

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.

  • $\begingroup$ That's odd, I'm quite sure the equation is correct, I've put a reference/link in the question comments to the section of the GR textbook which derives it. $\endgroup$
    – tbfr416
    Commented Mar 19, 2019 at 13:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.