# hard sphere collision condition with periodic boundary conditions

Hello Stack Exchange Physics,

I am curious if there is an analytical or efficient numerical solution for the collision of hard spheres in a rectangular unit cell with periodic boundary conditions.

Consider a system with two hard spheres. The first sphere is located at $\vec r_i$ with a velocity $\vec v_i$ and a radius $\sigma_i$. The second sphere is located at $\vec r_j$ with a velocity $\vec v_j$ and a radius $\sigma_j$. Note $\vec v_{ij} = v_i - v_j$, $\vec r_{ij} = r_i - r_j$, and $b_{ij} = r_{ij} \cdot v_{ij}$.

The hard sphere collision condition is well known. If $\vec r_{ij} \cdot \vec v_{ij} < 0$ and $b_{ij}^{2} - v_{ij}^{2}(r_{ij}^2-\sigma^2) > 0$, the spheres will collide. The time of collision is: $$\frac{-b_{ij}-\sqrt{b_{ij}^{2}-v_{ij}^{2}(r_{ij}^{2}-\sigma^{2})}}{v_{ij}^{2}}$$

How is this condition modified if the spheres are enclosed in a rectangular unit cell $L_x \times L_y \times L_z$ with periodic boundary conditions in the $x$, $y$, and $z$ directions? What is the most efficient way to determine if and when the spheres will collide?

Best, Eric

Only in very special cases will the balls not collide eventually. If the relative velocity cannot support a periodic orbit, then they will eventually collide.

So if the initial velocities are randomly distributed in the reals, since the density of the rationals is zero in the reals, then there is 0 probability that they don't collide.

Even if the relative velocity can support a periodic orbit, they could still collide. Conceptually, you can change coordinates so that one ball is sitting still at the center of the box. The other ball is travelling at the relative velocity.

Now just calculate line segments of its path, and check each line segment's distance from the origin. If it gets closer than the two radii before repeating, there is a collision. There may be a slightly faster way using number theory of the periods, but I'd have to think more on that.

Update:
Okay, putting this all together. It is possible to get an answer without simulating the entire path of the particle till it repeats. But it isn't all that simple either.

Let's change coordinate systems. Have one ball at rest at the origin in the center of the periodic cell ( $-L_x/2 < x \le L_x/2$, etc.), and the other ball is moving with velocity $\vec{v}$ and initial position $\vec{r}$. Also define the closest the two balls can come to each other $R$ (the sum of the two radii).

First the simple case. If $\vec{v}=0$, just check if $|r|<R$ to determine collision or not.

Otherwise if any component of $v$ is zero (for example $z$), check if that component of the relative position is out of range ($r_z > R$ or $r_z < -R$), in which case there is no collision.

Now we get to the more general cases.

In periodic space, a path can have the relative distance between the balls increase and decrease multiple times. So for discussion here let's consider the "impact parameter" just a function that takes in the current relative position and velocity, and gives the minimum distance of approach for the path in regular space. This way it is unambiguous and also makes it easy to calculate:

$b = |\vec{r}\times\vec{v}|/|v|$

In periodic space, we can imagine kind of hitting an imaginary boundary and then shifting the position to stay in our periodic cell. The velocity vector however is constant. So overall the impact parameter can change each time we hit a boundary. The locations where $b=|r|$ define a plane (where $\vec{r}\cdot\vec{v}=0$).

Now this plane in periodic space may hit the boundaries such that it looks broken/duplicated into parallel plane segments in the periodic cell, but if the closest approaching path to the origin goes through one of these planes it will necessarily also go through the plane segment on the origin. So minimizing on this plane segment is equivalent to find closest approach over all, allowing us to determine if the balls collide. Since for points on this plane $b=|r|$, I'll call this the impact plane.

There will always be at least two axes not parallel to v. Chose two of these coordinates to parameterize points on the plane. Without loss of generality let's call these x and y.

Now move the ball according to $v$, shifting position to remain in the cell every time it hits a boundary, until it hits the impact plane. Record its position as $A$.

Do that again, but start the ball at the origin. Move until it hits the impact plane. Record its position as $B$.

Define $c_x = B_x / L_x$, $c_y = B_y / L_y$.

If both $c_x$ and $c_y$ are irrational, the balls will collide.

If one is irrational (or is zero), then that coordinate will eventually reach anything (or is a known constant), so we only have one dimension left to optimize in. In that case find the $c_x = n/m$ such that gcd(n,m)=1 (denominator, numerator have no common factors). This means the step size of this one dimensional lattics is $s_x = L_x/m$. Then see if the point on the plane for $x=A_x - \text{floor}(A_x/s_x) s_x$ (and $y=0$ in the irrational case, or the intial $A_y$ in the case that $c_y$ is zero) is close enough to collide. Also check for $x=A_x - \text{ceil}(A_x/s_x) s_x$.

Okay, if we're here $c_x$ and $c_y$ are both rational and non-zero.

Write $c_x=n/m$ and $c_y=p/q$ such that gcd(n,m)=gcd(p,q)=1. Imagine if we took m of these steps. If gcd(m,q)=1, then m steps will wrap around to the same x value, but y will step $(mp \text{ mod } q)L_y/q$ where "mod q" converts an integer to the value in (-q/2,+q/2] which is equivalent mod q. Since both m and p were relatively prime to q, we can continue taking steps wrapping around till the net change was just $L_y/q$. We can do similarly for x as well. So we have a 2D lattice with primitive basis $s_{1x}=L_x/m,s_{1y}=0$ and $s_{2x}=0,s_{2y}=L_y/q$.

However if gcd(m,q)!=1, define $s_{1x} = L_x/m, s_{1y} = (\text{modinv}(n,m)p\text{ mod }q)L_y/q$ and similarly $s_{2x} = (\text{modinv}(p,q)n \text{ mod } m)L_x/m, s_{2y} = L_y/q$. Notice if $c_x=\pm c_y$ that means this again reduced to a one dimensional lattice. Anyway, now with the primitive lattice defined, find the points near to the origin on the plane and check if they are close enough to collide.

We parameterized the points on the plane segment with x,y so that we'd know how they wrapped. But now that we have them, in the two dimensional case it is probably best to write the lattice vectors not in terms of x,y on the plane, but in terms of some orthonormal basis in the plane. Since $b=|r|$ this will allow reasoning directly about the impact parameter instead of needing to find the z component that puts those x,y on the plane and then calculating the impact parameter each time. With that in hand, it should be possible to determine if there was a collision only checking maybe 6 or so points.

So, that was long. But it does look possible to check for a collision without simulating the entire path of the particle while waiting for a repeat, yet it remains more complicated than the non-periodic case.

• Thanks for the response, BuddyJohn. I agree, nearly all trajectories should result in a collision. The recursive method you suggest should work but may be inefficient. I am wondering if there is another way. Since the periodic boundaries only shift the position vector of the non-stationary sphere by the relevant unit cell length, I was wondering if some clever linear algebra could be used to determine the collision time with a single calculation. Jan 10, 2016 at 8:00
• Not with linear algebra. Before it could be done because after the change in coordinates, there is only one line, and we're just checking its impact parameter. With the periodic boundary conditions, it is equivalent to duplicating a ball an infinite number of times (a lattice), and checking if the other ball ever hits one. I'll try to work out the number theory method. But it will still involve a series of steps (a function in the programming sense, but not a nice algebraic expression). Jan 10, 2016 at 8:11
• Very interested in your number theory solution. Hope to hear back from you. Jan 10, 2016 at 8:16
• Wow, that got long. Only after writing it up did I notice you asked not only "if" the spheres collide, but "when". The when looks much harder to me, because then it is not enough to determine a collision will occur but (in terms of viewing the puzzle as infinite copies of the balls in 'normal space') we now need to determine precisely which one hit first. For the rational case, the above may still help. For the real valued case though, this starts to sound like an integer optimization problem. Jan 10, 2016 at 14:11