The thing that one should keep in mind, is that forces gravitate. So if we eliminate dynamics by considering a static system of black holes that are “held fixed” then the thing that prevents those black holes from falling toward each other would necessarily contribute to the curvature of spacetime and would distort the geometry of black holes.
Having said that, there is a class of static vacuum solutions of Einstein equations with axial symmetry, the Weyl metrics, that allows one to consider superposition of two static black holes held at a fixed distance apart. This is because the metrics are written in terms of a single potential, $\psi$, that satisfies a linear flat-space Laplace equation written in cylindrical coordinates. For a single Schwarzschild black hole the source of the potential is a rod of length $2M$ placed along the symmetry axis. By using as a source two rods along the same axis but at some distance apart we would obtain the potential for the “superposition” of two Schwarzschild black holes.
However, if we integrate all the components of the metric from the potential we would find that it must have cosmic string-like conical singularity(-ies) along the symmetry axis, either between black holes or stretching to infinity from one of (or from both) black holes. The physical meaning of such conical singularities is simple: this is a string (or a strut) or a system of strings/struts enabling the black holes to remain static without falling toward each other, conical singularity corresponding to $\delta$-like distributional stress-energy tensor. But this enables us to find the force that the string/strut exerts on the black holes from the conical singularity's angle deficit/excess. The result for the force between two black holes with masses $M$ and $μ$ is:
$$
F_z = \frac{M μ }{R^2 - (μ+M)^2},
$$
where $R$ is the coordinate separation between centers of the rods in Weyl coordinates and a unit system with $G=c=1$ is used. This force diverges when horizons almost touch each other, and in the limit of large separation it approaches the expected Newtonian expression. In the limit of one of the masses being much smaller than the other, the expression would approach the one provided by OP after changing from Weyl to Schwarzschild coordinates (see e.g. this paper for an example of such limit that also includes the $m^2$ corrections).
As a sanity check let us consider if two expressions agree with each other when applicable (when $μ\ll M$). The equations needed are in Section III of LaHaye & Poisson paper cited above. First, we note that the mass $μ$ in the above expression does not have the same meaning as $m$ in OP, because $μ$ is the contribution from the second black hole to the ADM mass of the entire system (which is thus $M+μ$), so for a light second black hole ($μ\ll M$) it is its Killing energy, whereas $m$ in the formula from the question is the norm of particle's 4-momentum:
$$
μ=E_\text{K}=m u_\alpha t^\alpha,
$$
where $u^\alpha$ is the 4-velocity and $t^\alpha$ is timelike Killing vector (in static coordinates $t^\alpha=\partial_t$). In Schwarzschild coordinates small mass located at $r=d$ would have:
$$
u^\alpha=\frac{t^\alpha}{\sqrt{f}},\quad E_\text{K}=m\sqrt{f},\quad a= \frac{M}{d^2\sqrt{f}},
$$
where $f=1-\frac{2M}{d}$, and $a$ is the norm of 4-acceleration, which has only radial component. To obtain the force needed to keep the particle still we just take $F=ma$, recovering expression from OP. That same particle when viewed in Weyl coordinates of background black hole (particle is assumed to be on the symmetry axis) would have $\rho=0$, $z=d-M=R$ ($R$ is the separation as defined above).
$$
\mu = m \sqrt{\frac{R-M}{R+M}},\quad F=\frac{Mm}{(R-M)^{1/2}(R+M)^{3/2}}=\frac{M\mu}{R^2-M^2},
$$
which coincides with linear in $\mu$ part of the above expression for the force $F_z$ as expected.
References
Discussion of Weyl metrics could be found in a book:
- J. B. Griffiths and J. Podolský, Exact space-times in Einstein’s general relativity (Cambridge University Press, Cambridge, 2009), Ch. 10.
Note, there are a lot of sources discussing Weyl metrics going back to original works of H. Weyl circa 1919 so here is a couple of papers that discuss the forces and has open access versions:
P.S. Letelier and S.R. Oliveira. Superposition of Weyl solutions: the equilibrium forces. Classical and Quantum Gravity 15, no. 2 (1998): 421, arXiv:gr-qc/9710122.
P. Krtouš and A. Zelnikov, Thermodynamics of two black holes, JHEP 02, 164 (2020) arXiv:1909.13467.
