Jumping ahead on a wide orbit

A grazing encounter between a star and a black hole often leaves the star bound, on an orbit so eccentric and so wide that the next pericentre passage could be tens of thousands of dynamical times or more away. For all but the last little while of the orbit, the two bodies are a two-body problem with a star quietly ringing inside it.

StarSmasher can skip the part handled by a Kepler two-body orbit. The routine is jumpahead, in parallel_bleeding_edge/src/skipahead.f. It measures the orbit the two components are on, solves the two-body problem analytically, and puts the system back down at a smaller separation on the infalling branch of the same orbit, with every particle’s position and velocity relative to its own component untouched. The hydrodynamics resumes from there.

Where this has been used

The technique goes back to the binary-disruption simulations of Antonini, Lombardi & Merritt (2011), ApJ 731, 128, whose Section 3.3, “Timescale considerations and orbital advancement”, sets out the argument: the thermal timescale of a bound star is \(10^5\) to \(10^7\) yr, far longer than an orbital period, so the star’s structure barely changes over an orbit and nothing is lost by advancing it analytically. They wait at least eight days after periapsis before measuring the orbital elements, and check the approximation against runs that integrate the orbit in full.

Godet et al. (2014), ApJ 793, 105 use the same treatment for the repeated partial stripping of a donor by an intermediate-mass black hole in HLX-1. Their Section 6.1 states it compactly: once the donor “has retreated sufficiently far from the black hole to become stabilized (typically about 100 dynamical timescales after periapsis), we employ the analytic Kepler two-body result to advance the orbit to the same separation but now with the donor infalling toward the BH”, and “during this advancement of the orbit, we excise from the simulation any particles that have been stripped from the star”. Sections 5 and 6 of the same paper discuss how sensitive the orbital evolution is to the phase of the star’s oscillation at pericentre, which is worth reading before deciding what to do about artificial viscosity.

Kıroğlu et al. (2023), ApJ 948, 89 call it orbital regularization in their Section 2.2, and say most clearly why it is needed: for \(M_{\rm BH} > 100\,M_\odot\) at \(r_p = r_{\rm T}\) the remnant comes away with \(e > 0.999\) and \(a > 10^4\,R_\odot\), an orbital period of order ten years, or \(\sim 10^5\) dynamical timescales. They also record the choices that go with it: jump once the star has receded far enough that the orbital elements are well determined and the star is close to hydrostatic equilibrium; turn artificial viscosity off for the leg that follows, so that oscillations in the remnant can be followed cleanly; and accept that all debris bound to the black hole is treated as accreted, which is harmless when \(M_{\rm BH} \gg M_*\).

Turning it on

All four parameters below are set in sph.input.

tjumpahead: The time at which the jump happens. The default, 1d30, never fires. Any other value is a deliberate request and is honoured. The code stores it negated, and that negative sign is what tells changetf, which revises the run’s own schedule at every output, to leave the jump time alone. main.f compares against its absolute value, so you write it positive and never see the sign.
throwaway: Whether the debris is discarded. The default is .true., which is what the papers do. See Discarding the debris.
internal_energy_fraction: How much of a particle’s specific internal energy is allowed to help unbind it when the particles are sorted into components. The default is 0. See Which particle belongs to which body.
tf: Not required for a jump, but set it negative anyway. A negative tf lets the code revise its own stopping time, and makes it analyse the system at every output and write ecc.sph, which is how you follow the orbit.

That is the whole interface. Put tjumpahead in sph.input and run; the jump fires on the first iteration past it. It can go in from the start, or be added later and the run resumed from restartrad.sph, which is useful when you would rather look at the first passage before committing to a jump time.

Note

changetf can in principle schedule a jump on its own: when it sees a bound pair receding with an apocentre past 1000 code units it logs FUTURE CANDIDATE FOR JUMPING AHEAD. The line that would set a jump time there is commented out, so nothing follows from it, and the choice stays with you.

Choosing when to jump

Jump too early and the orbital elements are still changing and the star is not back in equilibrium. Jump too late and you have paid for the integration you were trying to avoid.

A separation that works

The prescription behind the published intermediate-mass black hole runs is a separation,

\[r_{\rm jump} = 1.7 \left(\frac{M_{\rm BH}}{M_*}\right)^{1/3} \max(r_p,\, r_{\rm T}),\]

which for \(r_p \le r_{\rm T}\) is the same as \(1.7\,(M_{\rm BH}/M_*)^{2/3} R_*\). The coefficient comes from watching what the debris does. Measuring the star-black hole separation in \(r_p = r_{\rm T}\) runs with a 1 \(M_\odot\) star, and writing \(\sigma\) for \((M_{\rm BH}/M_*)^{1/3}\max(r_p,r_{\rm T})\):

Separation	What has happened by then
\(1.4\,\sigma\)	the first bound debris has fallen back to the black hole
\(1.7\,\sigma\)	the debris stream is just starting to self-intersect, and the disc radius has stopped growing
\(2.0\,\sigma\)	that material has been round once more

1.7 is the last moment before shocks from self-intersection begin to matter, which is what makes it the right place to jump in a run with the artificial viscosity turned off: there is nothing yet for the viscosity to do. It is the same criterion Kıroğlu et al. describe from the other end, noting that for \(r_p \gtrsim r_{\rm T}\) that separation is reached just as the stream starts to self-intersect.

For \(M_{\rm BH}/M_* = 5, 10, 100, 200, 500\) and 1000 the coefficient \(1.7(M_{\rm BH}/M_*)^{1/3}\) comes to 2.9, 3.7, 7.9, 9.9, 13.5 and 17. At the low end that is only a few tidal radii, close enough that the star may not have finished being disrupted, so look at a snapshot before trusting the number.

Turning a separation into a time

Barker’s equation is not required. tjumpahead is a time and nothing else, and the code knows nothing of \(r_{\rm jump}\) or of Barker’s equation. If you already know how long you want to wait (because you watched the first passage go by, or because a previous run of the same encounter told you) then write that time in and skip to the next section. What follows is only the way to turn a separation you have picked into a time, which you need when you are choosing the jump time before the run exists.

The first passage is near enough parabolic for Barker’s equation, the parabolic counterpart of Kepler’s equation (Pathan 2008, Math. Gaz. 92, 39, or Section 4.5 of Bate, Mueller & White, Fundamentals of Astrodynamics). Written in terms of \(x = r/r_p\) it gives the time since pericentre as

\[t = \frac{\sqrt{2}}{6\pi}\,(x+2)\sqrt{x-1}\; t_{\rm orb} = 0.07503\,(x+2)\sqrt{x-1}\; t_{\rm orb}, \qquad t_{\rm orb} = 2\pi\sqrt{\frac{r_p^3}{GM_{\rm BH}}}.\]

Putting \(x_{\rm jump} = 1.7(M_{\rm BH}/M_*)^{1/3}\) into it gives \(t_{\rm jump} = 1.9\,t_{\rm orb}\) for \(M_{\rm BH} = 100\,M_\odot\) and \(5.7\,t_{\rm orb}\) for \(1000\,M_\odot\), measured from pericentre. The same expression run with the starting separation gives the time from the start of the run to pericentre, so the two together give a tjumpahead before the run has been started.

A worked example

The run below takes under an hour on one GPU and uses only files that come with the repository. It sends the relaxed 8 \(M_\odot\) star from example_input/collision past a 100 \(M_\odot\) black hole, which is what hyp produces when sph.start2u is absent: the second body becomes a single point mass of mass mbh.

Working out the jump time

The star is \(R_* = 3.17\,R_\odot\), so \(r_{\rm T} = R_*(M_{\rm BH}/M_*)^{1/3} = 7.4\). Taking \(r_p = 12\), a grazing pass at about \(1.6\,r_{\rm T}\), the prescription gives \(r_{\rm jump} = 1.7 \times 2.32 \times 12 = 47\). With \(t_{\rm orb} = 2\pi\sqrt{12^3/100} = 26\) code units, Barker’s equation puts \(x = 47/12\) at 20 code units after pericentre, and the starting separation \(x = 60/12\) at 27 before it. So the jump wants to happen around \(t = 47\); the orbit here is bound rather than exactly parabolic, which brings pericentre in a little earlier, and 45 is a round number close enough.

Setting up

Put the executable and the star in an empty directory:

$ mkdir jump && cd jump
$ cp ../parallel_bleeding_edge/parallel_bleeding_edge_gpu_sph .
$ cp ../example_input/collision/sph.start1u .

sph.init:

&INITT
INAME='hyp' ! "hyperbolic" collision (also works for parabolic and eccentric encounters)
&END

sph.input. Note the negative tf, and that sph.start2u is deliberately absent (so that a black hole is used instead):

&input
tf=-9999, ! a negative final time can help to make things more automatic
dtout=10, ! time in code units between output files
sep0=60, ! initial separation
rp=12.0d0, ! periapsis distance
e0=0.995d0, ! initial eccentricity
mbh=100.0d0, ! black hole mass in solar masses
tjumpahead=45.0, ! time in code units to do an orbital jump
&end

throwaway needs no line: its default is already .true.. Run it:

$ mpirun -np 4 ./parallel_bleeding_edge_gpu_sph

The first passage

Reading the out*.sph snapshots back gives the orbit tightening as it passes pericentre, which is the tidal energy going into the star:

out0000.sph  t=  0.06  r= 59.93  rdot=-1.69  a= 2515.6  e=0.99523
out0001.sph  t= 10.02  r= 42.12  rdot=-1.90  a= 2528.9  e=0.99526
out0002.sph  t= 20.02  r= 21.82  rdot=-2.11  a= 4166.8  e=0.99712
out0003.sph  t= 30.01  r= 15.67  rdot=+1.79  a= 1269.8  e=0.99054
out0004.sph  t= 40.00  r= 36.09  rdot=+1.98  a= 1199.7  e=0.99002

Pericentre is at \(t\approx24\), and by \(t=40\) the semimajor axis has dropped from 2516 to 1200. ecc.sph carries the same story one row per output. Once the star is receding, log0.sph adds this at every output; the numbers below are the ones from \(t=40\):

bound orbit with orbital period=   44277.381515463341
the stars will take a long time to orbit, we might give up: ecc=  0.99313867634119102

The jump

Everything that follows is from log0.sph. First the announcement, then compbest3 sorting the particles, which takes three passes to converge:

jumpping ahead at time t=   45.000532609753130
...
          nit        nchng           m1           m2           m3
            0        19702  7.98693      100.012      0.00000
            1           88  7.98637      100.013      0.00000
            2            0  7.98637      100.013      0.00000

The star is left with \(7.98637\,M_\odot\), and the black hole’s component has gained \(0.0130\,M_\odot\) of debris bound to it. Almost none of the stripped mass is unbound: the mejecta= reported at the previous output is \(3.9\times10^{-4}\).

Then the orbit, measured and solved:

eccentricity: r12=   45.7046568     am1=   7.98637468     am2=   100.012968
reduced mass mu=   7.3957953400447867
total orbital energy= -0.2275
total angular momentum=    376.1      -0.4639E-03   0.7433E-03    376.1
components of eccentricity vector=  -0.9932      -0.3172E-02  -0.1219E-05
apastron separation rmax=   3498.6677841672463
sep0=   22.852328392088165
1st simple check: -0.22751752517881471      -0.22751752517881466
semilatusrectum=   23.940548138731852   mu=   7.3957953400447867
  sep0=   22.852328392088165   ecc=  0.99315723880756224
cos=   4.7947739105261220E-002  sin= -0.99884984572992441
rdot(from semilatusrectum): -2.1069863876044512  rdot(from e): -2.1069863876044521

The separation has been halved, 45.70 to 22.85. \(\sin\theta\) is negative and so is \(\dot r\), which is what puts the star on the infalling branch, and the two independent ways of computing \(\dot r\) agree to sixteen digits. The orbital energy agrees with the value reconstructed from \(e\), \(\mu\) and \(L\) to the same precision.

Then the three energy analyses, before the jump, after it, and after the debris has been discarded. These are also the three lines written to jumpahead.sph, in the columns of energy*.sph:

           t           epot           ekin           eint           etot         ajtot
99822      -43.92045      17.29951       13.53356      -13.08738      376.6425
99822      -61.39812      34.77565       13.53356      -13.08891      376.6425
99822      -61.29525      34.73058       13.53313      -13.03154      376.2588

Line 1 to line 2 is the jump. The internal energy is identical to every digit printed, because nothing was done to the star, and so is the total angular momentum. The potential and kinetic energies both change, and should: the star has been moved from \(r=45.7\) to \(r=22.9\), so it is deeper in the black hole’s potential and moving faster. The total energy moves in the fifth digit.

Line 2 to line 3 is the mass removal (if throwaway is set to .true.):

now we throw away some particles:
mass thrownaway=   1.3416840884748770E-002
new ntot=       17729

etot falls from \(-13.0889\) to \(-13.0315\) and ajtot from 376.64 to 376.26, which is what the deleted particles carried off with them.

m1m2rp.sph records the result:

17728.000000000000        1.0000000000000000        12.000000000000000
7.9863746823827269        100.01296764968698        12.000000000000000

17728 SPH particles in the star, one point mass, and rp=12; then the component masses. The second component weighs \(100.0130\), the black hole plus the debris bound to it, and that is the mass the Kepler solve used. The point particle itself is still exactly \(100\), so the leg that follows is integrated with the original mbh.

What it bought

With \(a = 1755\) and \(e = 0.99316\) the orbital period is 44464 code units. Kepler’s equation gives the time of flight from \(r=45.7\) outbound, around apocentre at 3499, back to \(r=22.9\) infalling: 44439 code units, or 99.94% of a full period. The whole run up to the jump spans 45.

The snapshots after the jump show the star going straight back in and round again:

out0005.sph  t= 50.00  ntot=17729  r= 13.44  rdot=-1.30  a=1748.2  e=0.99313
out0006.sph  t= 60.00  ntot=17729  r= 24.75  rdot=+2.10  a= 881.6  e=0.98641

The second pericentre passage is over by \(t\approx53\), and the semimajor axis has halved again. Without the jump it would have arrived somewhere near \(t = 45000\).

How it works

Which particle belongs to which body

The Kepler solve needs two bodies, so before anything else the particles have to be divided between them. compbest3, in compbest3.f, determines the mass, the center of mass, and center of mass velocity of each of the two components.

The algorithm is the one set out in Section 2.2, “Analysis of Hydrodynamics”, of Kremer et al. (2022), ApJ 933, 203. The code does a two-body binding test applied particle by particle.

Each component \(j\) has a mass \(M_j\) and a centre. The centre of the star is its densest particle; the centre of the black hole and whatever is bound to it is the point particle itself. Particle \(i\) counts as bound to component \(j\) when

\[\frac{1}{2}v_{ij}^{2} + f_u\,u_i - \frac{GM_j}{d_{ij}} < 0,\]

where \(v_{ij}\) is the particle’s velocity relative to that component’s center-of-mass velocity, \(d_{ij}\) is its distance from that component’s center, \(u_i\) is its specific internal energy, and \(f_u\) is the sph.input parameter internal_energy_fraction. (The array holding \(f_u u_i\) is called enth in the source, but what it scales is specific internal energy and not enthalpy: whichever variable the run actually integrates, \(A\), \(\ln A\) or \(u\), is converted back to \(u_i\) first.)

Three rules settle the rest:

A particle bound to more than one component goes to the one whose center is closer, not to the one it is more tightly bound to.
A particle bound to neither body is ejecta, which the code calls component 4.
A point mass never changes component. Point masses are sorted once, on the first call, and stay put, so the black hole anchors its component for good.

The masses and centers are then recomputed from the new membership and the test repeated, until a pass moves nobody. That loop is the nit/nchng table in log0.sph.

A body left holding between one and four SPH particles is dissolved into the ejecta, on the grounds that it is not a star.

Kremer et al. use \(f_u = 1\). The default here is 0, following Section 3.2 of Nandez, Ivanova & Lombardi (2014), ApJ 786, 39, who set the two side by side as the “conventional” and “abridged” definitions and adopt the abridged one. Their objection to counting \(u_i\) is that it declares material unbound on the strength of heat it never gets to spend: the particles it mislabels are shock-heated ones sitting beside the companion with a large \(u_i\) and almost no velocity, and they are still sitting there at the end of the run, having neither radiated that heat nor passed it to a neighbour. With \(f_u = 0\) the internal energy of the ejecta falls as adiabatic expansion says it should; with \(f_u = 1\) it stays high.

At a well-chosen jump time the setting should barely matter, and that is an argument for jumping late rather than for ignoring the parameter. What the test weighs is \(f_u u_i\) against \(GM_j/d_{ij}\), and by \(r_{\rm jump}\) those two are far apart for almost every particle: the star has had tens of dynamical times to return to hydrostatic equilibrium, so its material is bound by a wide margin, and the debris has expanded and cooled adiabatically, so \(u_i\) out there is small. Only material still hot from the pericentre shock lies near the boundary, and waiting is what removes it. If flipping \(f_u\) between 0 and 1 moves the masses in m1m2rp.sph appreciably, the jump was made too early.

The two-body solve

From the two component masses and their centre-of-mass positions and velocities the routine forms the reduced mass \(\mu\), the orbital energy \(E_{\rm orb}\), the angular momentum \(\mathbf{L}\), and the eccentricity vector

\[\mathbf{e} = \frac{\mathbf{v}\times\mathbf{L}}{k} - \hat{\mathbf{r}}, \qquad k = G m_1 m_2,\]

the Laplace-Runge-Lenz vector scaled to have magnitude \(e\). It is conserved in the Kepler problem and points at pericentre, so it fixes the orientation of the orbit as well as its shape. As a check, \(e^2\) computed from the energy and angular momentum is compared against \(|\mathbf{e}|^2\), and the run stops if they disagree by more than \(10^{-14}\).

The new separation is then chosen, and the true anomaly that goes with it follows from the orbit equation,

\[\cos\theta = \frac{1}{e}\left(\frac{\alpha}{r} - 1\right), \qquad \alpha = \frac{L^2}{\mu k},\]

with \(\theta\) taken negative. \(\theta = 0\) is pericentre, so a negative \(\theta\) is the pre-pericentre branch: the star is placed where it is falling in, not where it is climbing out. The radial velocity \(\dot r\) is taken negative for the same reason.

Note

The papers advance the orbit to the same separation the star had reached, and so did the version of the code the published runs were made with. The code as it stands halves it: sep0=0.5d0*r12, with the unhalved sep0=r12 commented out on the line below. Halving skips more of the orbit, at the cost of resuming the hydrodynamics closer in, and it breaks the property the published runs relied on, that the separation after a jump equals the separation before it, so each passage takes the same time to reach the next pericentre. That is the line to change to get it back.

Each particle is then translated and boosted by the difference between its component’s new and old centre-of-mass state. Every particle keeps its position and velocity relative to its own component, so the star’s internal structure, its spin, and whatever oscillation the encounter left it ringing with all carry across untouched. That is the approximation the whole method rests on, and it is what Antonini et al. justify with the thermal timescale argument.

Warning

The simulation clock does not advance. jumpahead changes t by about one timestep and no more. The orbital time that was skipped, which is most of a period, is simply not counted. Every time in log*.sph, energy*.sph, ecc.sph and the out*.sph headers after a jump is a time with the wide part of the orbit cut out of it. If you need real elapsed time, add the Kepler time of flight yourself, from the \(a\) and \(e\) printed in the log.

Discarding the debris

With throwaway=.true., the default, every particle that is neither a point mass nor a member of the surviving body’s component is deleted, the particle count is compacted, and the mass removed is reported as mass thrownaway=. That covers both the material thrown to infinity and the material left bound to the black hole. If fewer than nnopt particles survive, the run stops rather than continue with a star it cannot resolve. When there is no point mass anywhere, as in a star-star encounter, there is no accretor either, so both bodies are kept and only material bound to neither goes.

With throwaway=.false. nothing is deleted, and the particles in neither body are each moved with whichever body they are more tightly bound to. The count is logged as throwaway is off, so N particles belonging to neither body are carried along. That is an approximation, since such a particle is not on either body’s orbit, but a better one than leaving the debris behind while both bodies are moved out from under it.

Because nothing is removed, nothing is renumbered either: with throwaway=.false. the particle order is not changed by a jump. Particle \(i\) in the snapshots after the jump is the same particle \(i\) as before it, and ntot is the same too, so a particle can be tracked straight through. With throwaway=.true. the survivors are compacted into a contiguous block, which keeps them in their existing relative order but shifts every index past the first deletion, so indices must not be compared across a jump.

The Kepler solve uses the component masses from compbest3, so debris bound to the black hole is counted in the black hole’s mass for the purpose of the orbit. The point particle’s own mass is not increased, so the leg after the jump is integrated with the original mbh. That is the approximation Kıroğlu et al. describe and defend on the grounds that \(M_{\rm BH}\gg M_*\); if your mass ratio is less extreme, this is the place to look first.

How much of the bound debris ought to count as accreted is a real question, and the answer is “not all of it”. Ayal, Livio & Piran (2000), ApJ 545, 772 find that around 75% of the returned debris becomes unbound again for supermassive black holes, and measurements on \(M_{\rm BH} = 500\,M_\odot\) runs here put the fraction that stays bound at 10% or less. Since the black hole is far more massive than anything being argued about, the hydrodynamics is essentially unaffected by the choice, and accretion rates can be rescaled afterwards by whatever factor a later reader prefers.

Warning

Discarding the debris is not a neutral bookkeeping choice for the orbit. Test runs that kept the material bound to the black hole ended up on markedly more tightly bound orbits, and in cases where discarding it led to the star being ejected after a few passages, keeping it postponed the ejection or prevented it. The question to ask is how the orbital period compares with the timescale on which the disc would really be cleared, by accretion and by feedback winds, neither of which is in the simulation. For intermediate-mass black holes the orbital period is long enough that discarding is the better approximation. For stellar-mass black holes it is much less clear, and the two settings are worth running against each other.

Checking that a jump was sane

The routine stops the run itself if the reconstructed orbit does not match the measured one, so the arithmetic needs no checking. What is worth checking is everything it does not guarantee:

The three lines in jumpahead.sph. eint and ajtot should be unchanged across the jump to every digit printed. epot and ekin will both change, since the separation changed. etot should barely move; if it moves by a percent, the two components are not well described as point masses, which usually means the jump was made too soon after pericentre.
The component masses in m1m2rp.sph. The excess of the accretor’s component over mbh is the debris treated as accreted, and is worth knowing: it is the mass the orbit was solved with but the point particle was not given.
That the star was really relaxed when you jumped. Nothing in the code checks this, and it is the assumption that matters most.
That enough particles are left. Each jump with throwaway=.true. removes some, and several passages in a row can leave a star without enough particles to find neighbours; an average neighbour count in single figures is not doing SPH any more, whatever the run says. The code stops with not enough particles to continue only when fewer than nnopt survive, which is later than you want to find out. ls -l out*.sph is the quick check: with throwaway on the snapshot size drops at every jump, so the file sizes show both where the jumps were and how much was lost.

There is also a self-test built into the routine. skipahead.f carries two commented-out lines, sintheta=-sintheta and rdot=-rdot, flagged in the source as a pair. Uncomment both and the Kepler solve puts the system back exactly where it already was, so the jump becomes the identity while everything around it still runs: the component split, the rotation, the particle move, the five checks. Anything that moves is then a bug. It is the first thing to reach for after touching the routine.

The same component finder runs in the post-processing. compbest3.f is shipped with the splot routines in splot_routines/ as well, where option 71 calls it to build massAndMore.out. Bound masses quoted from a jump and bound masses measured afterwards therefore come from the same algorithm. That copy has \(f_u = 0\) built in and reads no sph.input, so it agrees with a run left at the default and not with one that sets internal_energy_fraction=1.