Gravitational potential inside infinite grid of masses

Question

Assume an infinite two-dimensional grid with nodes at $i,j\in \mathbb{Z}$, and each node $(i,j)$ hosting a mass $M_{ij}$. Can we find the total gravitational force $F(x,y)$ in closed form a small mass $m_{xy}$ at point $(x,y)$ feels on the plane, using Newton's formula for the force between two masses: $F=G\frac{m_1\cdot m_2}{r^2}$?

Preliminary notes:

Temporarily reducing the problem to the unit square and assuming unit masses $M_{ij}$, it is fairly obvious that the lines shown below (red) serve as axes of symmetry, so for each grid square, $F(x,y)=0$, iff $(x,y)$ falls on any of these lines.

Are there more lines that exhibit symmetry, different from the 8 red ones for each square?

Answer 1

The force $F_{ij}$ from mass $M_{ij}$ at position $r_{ij} = (i, j)^t$ applied to mass $m$ at position $r = (x,y)^t$ should be \begin{align} F_{ij} &= G \, m \, M_{ij} \frac{r_{ij}-r}{\lVert r_{ij}-r \rVert^3} \\ &= G \, m \, M_{ij} \frac{(i-x,j-y)}{\left(\sqrt{(i-x)^2+(j-y)^2}\right)^3} \end{align} The total force then is \begin{align} F &= \sum_{i,j\in\mathbb{Z}} F_{ij} \\ &= G \, m \, \sum_{i,j\in\mathbb{Z}} M_{ij} \frac{r_{ij}-r}{\lVert r_{ij}-r \rVert^3} \\ &= G \, m \, \sum_{i,j\in\mathbb{Z}} M_{ij} \frac{(i-x,j-y)}{\left(\sqrt{(i-x)^2+(j-y)^2}\right)^3} \end{align}

That sum looks ugly, so I asked my girl Ruby what she thinks about the resulting force for the simplified case $G = M_{ij} = m = 1$. I told her my view of the problem and she replied with this:

n = 98: r = Vector[0.5, 0.5], F = Vector[-0.014357405607155703, -0.01435740560715777]
n = 99: r = Vector[0.5, 0.5], F = Vector[-0.014213111904535183, -0.014213111904536712]
n = 100: r = Vector[0.5, 0.5], F = Vector[-0.014071689664539606, -0.01407168966453928]

n = 98: r = Vector[0.25, 0.25], F = Vector[-4.6817100348218865, -4.6817100348219425]
n = 99: r = Vector[0.25, 0.25], F = Vector[-4.681637888660634, -4.681637888660691]
n = 100: r = Vector[0.25, 0.25], F = Vector[-4.681567178203503, -4.681567178203558]

n = 98: r = Vector[0.25, 0.0], F = Vector[-14.708761955833403, -1.0966027348333074e-16]
n = 99: r = Vector[0.25, 0.0], F = Vector[-14.7086898091545, 1.0235546134620965e-16]
n = 100: r = Vector[0.25, 0.0], F = Vector[-14.708619098200186, 9.540979117872439e-17]

n = 98: r = Vector[0.0, 0.25], F = Vector[2.872627537318234e-15, -14.708761955833696]
n = 99: r = Vector[0.0, 0.25], F = Vector[1.2367358656270588e-16, -14.708689809154844]
n = 100: r = Vector[0.0, 0.25], F = Vector[-4.756347496848862e-15, -14.708619098200511]

So the resulting force at $(0.5,0.5)$ seems to vanish etc. This would make a nice 2D plot but I was too lazy to elaborate further.