Are the basins of attraction of two static gravitational sources two half-planes?

The good news is that your conjecture is true. The bad news is that my proof offers little insight. I'm not sure a small proof by some symmetry argument is possible, since it's not simply the symmetry of the setup that counts, the result is specific to inverse-square forces. Any perturbation to the exponent of the force leads to the particle skipping from side to side, generally behaving quite chaotically. This is reminiscent of the Laplace-Runge-Lenz vector, I wouldn't be surprised if there is some correspondence.

Anyway, first off, I'm assuming $GM = m = 1$. By symmetry considerations, we can assume the motion happens in a plane, so we'll just use two coordinates $x$ and $y$. Now let's use elliptical coordinates, of course placing the static masses at the foci $x = -a$ and $x = +a$

$$ x = a \cosh \mu \cos \nu ,\\ y = a \sinh \mu \sin \nu .$$

Working through the algebra (and with the help of Mathematica) our Lagrangian turns out to be

\begin{align*} \mathcal{L} &= \frac{\dot{x}^2 + \dot{y}^2}{2} + \frac{1}{\sqrt{(x-a)^2+y^2}}+\frac{1}{\sqrt{(x+a)^2+y^2}} \\ & = \frac{a^2 (\cosh^2 \mu - \cos^2 \nu) (\dot{\mu}^2+\dot{\nu}^2)}{2}+\frac{1}{a} \frac{2 \cosh \mu}{\cosh^2 \mu - \cos^2 \nu} \\ &= a^2 u \frac{\dot{\mu}^2+\dot{\nu}^2}{2} + \frac{1}{a} \frac{2 \cosh \mu}{u}. \end{align*}

For convenience, I have defined $u(\mu, \nu) = \cosh^2 \mu - \cos^2 \nu$. To formulate the Hamiltonian, we need the conjugate momenta

$$ p_\mu = \frac{\partial \mathcal{L}}{\partial \dot\mu} = a^2 u \dot{\mu} \\ p_\nu = \frac{\partial \mathcal{L}}{\partial \dot\nu} = a^2 u \dot{\nu}.$$

The Legendre transform yields

$$ \mathcal{H} = p_\mu \dot{\mu} + p_\nu \dot{\nu} - \mathcal{L} = \frac{1}{u} \left( \frac{p_\mu^2 + p_\nu^2}{2a^2} + \frac{2 \cosh \mu}{a}\right).$$

As this Hamiltonian is time-independent, we can define a conserved energy $E$ (which will be negative in our case). At this point, I'm uncertain about how exactly the next step works, and how correct it is. I could be mistaken. Taking inspiration from this paper, we introduce a new fictitious time coordinate $\tau$ satisfying

$$ d\tau = u dt ,$$

and a new Hamiltonian

$$ \mathcal{H}' = u(\mathcal{H} - E).$$

Taking partial derivatives we find an equivalent of Hamilton's equations

$$ -\frac{\partial \mathcal{H}'}{\partial q} = -\frac{\partial u}{\partial q}(\mathcal{H}-E) - u \frac{\partial \mathcal{H}}{\partial q} = u \frac{dp}{dt} = \frac{dp}{d\tau} \\ \frac{\partial \mathcal{H}'}{\partial p} = u \frac{\partial \mathcal{H}}{\partial q} = u \frac{dq}{dt} = \frac{dq}{d\tau}, $$

where we used the fact that $\mathcal{H} = E$ on-shell in the second step of the first line. These equations allow us to solve the equations of motion in function of the fictitious time coordinate $\tau$. Given the implicit definition of $\tau$, based on the coordinates itself, this might not seem like it would help, but for this purpose qualitative arguments based on the general form of the Hamiltonian will suffice.

Substituting in everything, in our case the new Hamiltonian looks like

$$ \mathcal{H}' = \frac{p_\mu^2 + p_\nu^2}{2 a^2} - \frac{2 \cosh \mu}{a} - E \cosh^2 \mu + E \cos^2 \nu .$$

It turns out $\mu$ and $\nu$ completely decouple, and basically look like two harmonic oscillators (up to a coordinate transformation)! If the two oscillation periods are incommensurate, as they will be for generic initial conditions, the particle will ergodically fill the region it is energetically allowed inside, coming arbitrarily close to any point inside of it, including the massive bodies themselves. This agrees with my simulations, an example of which is pictured. Here the blue line shows the trajectory of the particle, and the green and orange lines are curves of respectively constant $\mu$ and $\nu$.

The part corresponding to $\nu$,

$$ \mathcal{H}_\nu' = \frac{p_\mu^2}{2a^2} - |E| \cos^2 \nu , $$

implies that $\cos^2 \nu > \cos^2 \nu_0$. Since the half-plane $x = 0$ corresponds to $\nu = \pm \frac{\pi}{2}$, where the cosine achieves a minimum, this also proves that a particle starting with zero initial velocity will be confined to its own side.

finally

In the following $G=M=d=1$.

The equipotentials look like

For a particle starting from rest, the contours also represent the level curves of the first constant of motion, total energy $E$. Each point on the contour represents a starting point for a trajectory. In addition, a trajectory can't leave the interior of the contour. Since the contours for $V\le-2$ (the "$\infty$" curve) don't cross over to the other plane, for trajectories in these regions with $|x|>0$, $E$ alone segregates the basins into different plains.

Now consider the second constant of motion,

$R=r_1^2r_2^2\dot{\theta_1}\dot{\theta_2}-2(\cos\theta_1+\cos\theta_2)$

where the symbols are as as in- As before we plot the level curves for this constant and see what they look like. Each curve indicates a starting point of the trajectory.Each pair$(E,R)$ uniquely identifies the entire trajectory.

We observe that,$R$ seems to divide the starting points into 3 classes:

$R\gtrless0$ corresponds to $x\lessgtr0$.

with $R=0$ implying $x=0 $ or $(y=0,0<|x|<1)$. So trajectories starting in the right half-plane never enter the left as this requires a constant of motion to change sign. The $R=0$ with $(y=0,0<|x|<1)$ is explicitly checked and found also to satisfy the above.

Since the equipotential bounds the trajectory on the right, the particle is doomed to meet its corresponding source of field.

---

The formula used for proving the correspondence was

$$ R(x_0,y_0,\dot{x}=0,\dot{y}=0)= -2\left(\frac{x_0+1}{(x_0+1)^2+y_0^2}+\frac{x_0-1}{(x_0-1)^2+y_0^2}\right) $$

---

Earlier attempts

---

1

An object($|x_0|>0$) starting from rest, will only roam within the equipotential it was on at the beginning. So for points within $\infty$ contour($V=-2$), the answer is in the affirmative: no crossing over to the other plane. This is evident from the direction of forces too.

Now consider the field plot for the direction of the horizontal component of force

Trivially, the answer is affirmative for $y_0=0$ as y component of force vanishes.

The red dots(numerically calculated) indicate the boundary (excluding $x=0$), say $B$, where the horizontal component switches sign. Within this region, the particle always moves away from $x=0$.

If we could show that the particle never crosses B, then all such trajectories must terminate. Since the only point of singularity for these would be the mass location so answer be in the affirmative for these points too.

Alas! this isn't true.

---

2

So instead,
If we can show that the particle never crosses $y=\pm m x$ where $m=y_0/x_0$ and $0<x<x_0$, the basins would indeed be half planes. This seems to be indicated by sims.

work in progress

To prove the this hypothesis, we''ll try the following program.

1. First we show that the force $F$ on $y=\pm mx$ always points into the region $mx>|y|$. As a result at least initially, the particle isn't going to cross.

2. Let energy $E(x,y,\dot{x},\dot{y})$ and something else, $R(x,y,\dot{x},\dot{y})$ be the 2 constants of motion.

3. We can eliminate $\dot{y}$ in $R$ using $E=E(x_0,y_0,0,0)=E_0$ to get $R(x,y,\dot{x})$.

4. For any trajectory, $R=R(x_0,y_0,0)=R_0$. Now comes the important part:

   If the trajectory intersects $y=m x$ at another point $x'$, then $R'=R(x',m x',\dot{x'})=R_0$. We then analyze the behaviour of $R'-R_0 \forall 0<x'<x_0 ,-\infty<\dot{x'}<\infty $ hoping to find no real root.

5. A similar procedure for $y=-m x$ would bound the particle on the left while the equipotential already bounds it on the right.

The $R$ being used is
$R=r_1^2r_2^2\dot{\theta_1}\dot{\theta_2}-2(\cos\theta_1+\cos\theta_2)$
where the parameters are as indicated in the main answer above.

Currently, this procedure seems to fail. e.g. for $m=1, x_0=y_0=3$, here's what the plot is:

(the blue plane is $R_0=R'$) So for the following $x,\dot{x}$ values

there seem to exist values $(x',\dot{x'})$ for which, the trajectory swivels back to the line.

However, we should note that the existence of a negative $\dot{x}$ can't alone imply crossing as some large value of ${\dot{y}}$ may still keep the trajectory inbound.