Is there an approximate expression for the force between two black holes? Just curious: is there an approximate expression for the gravitational attraction between two Schwarzschild black holes of masses $M$ and $m$, held without relative speed at some center to center distance $d> 2G(M+m)/c^2$? Or is there even an exact expression?
In a book on relativity, a force expression is given as
$$F= \frac{GMm}{d^2} \frac{1}{\sqrt{1-\frac{2GM}{dc^2}}}.$$ But if I recall correctly, this is the force between a small test mass $m$ and a black hole $M$. Here, the question is about an expression for two black holes.  (A wild guess: Does the final expression just have two square roots, one with $M$ and one with $m$?)
In any case, the approximation should be more precise than the Newtonian expression $F= GMm/d^2$. The answer should take into account the two horizon radii. Newtonian gravity is not a good approximation at all when the distance $d$ is of the order of any of the two Schwarzschild radii. (Ideally, the force is given for two black holes held at constant distance, assuming that were possible.)
The question is not about dynamics, nor about the calculation of orbits. The question is only about the value of the force of attraction between black holes, e.g., while kept fixed. That way, no dynamics is involved. (If it helps to simplify the question, the black holes can also be of equal mass. They can also move in any desired way, such as towards each other. The black holes have no spin. The only quantity of interest in this question is the force.)
 A: 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.
