  An approximate analytical and non-linear solution of the Einstein field
equations is derived for a system of multiple non-rotating black holes. The
associated space-time has the same asymptotic structure as the Brill-Lindquist
initial data solution for multiple black holes. The system admits an
Arnowitt-Deser-Misner (ADM) Hamiltonian that can particularly evolve the
Brill-Lindquist solution over finite time intervals. The gravitational field of
this model may properly be referred to as a skeleton approximate solution of
the Einstein field equations. The approximation is based on a conformally flat
truncation, which excludes gravitational radiation, as well as a removal of
some additional gravitational field energy. After these two simplifications,
only source terms proportional to Dirac delta distributions remain in the
constraint equations. The skeleton Hamiltonian is exact in the test-body limit,
it leads to the Einsteinian dynamics up to the first post-Newtonian
approximation, and in the time-symmetric limit it gives the energy of the
Brill-Lindquist solution exactly. The skeleton model for binary systems may be
regarded as a kind of analytical counterpart to the numerical treatment of
orbiting Misner-Lindquist binary black holes proposed by Gourgoulhon,
Grandclement, and Bonazzola, even if they actually treat the corotating case.
Along circular orbits, the two-black-hole skeleton solution is quasi-stationary
and it fulfills the important property of equality of Komar and ADM masses.
Explicit calculations for the determination of the last stable circular orbit
of the binary system are performed up to the tenth post-Newtonian order within
the skeleton model.
