Do gravitational sources move along ‘geodesics’? Assume we have a system of say two bodies which are orbiting each other. Now assume that we wish to find an equation of the orbits of the two gravitational sources. Do they follow a ‘geodesical’ path, if we assume that the sources may or may not be singularities, which in this case may require a puncture method or so I’ve heard. 
I have also read several articles which suggest that it does move ‘geodesically’ if one does not take into account the self-interaction of the space-time field. If we were to take into account the self-interaction, back-reaction and what nots, will it still move ‘geodesically’ in the metric of the entire manifold?
 A: I think this is largely a matter of terminology.
If we take for example a test particle moving in a Schwarzschild metric then we can calculate the geodesics in the usual way. However an real particle has a non-zero mass (or energy) and therefore it perturbs the metric. So the particle is not moving in a Schwarzschild metric but instead in a time dependent metric that is similar to the Schwarzschild metric but not identical. So the actual path of the particle will be different from the geodesic calculated assuming a pure Schwarzschild geometry.
Typically we assume the particle is small compared to the mass it orbits, so the perturbation to the background metric is small. Typically we would assume the perturbation is approximately linear so we get:
$$ g_{\alpha\beta} = s_{\alpha\beta} + h_{\alpha\beta} $$
where $s_{\alpha\beta}$ is the Schwarzschild metric and $h_{\alpha\beta}$ is the perturbation caused by the non-zero mass of the orbiting particle. In that case:


*

*the trajectory of the particle is not a geodesic of $s_{\alpha\beta}$

*but the trajectory of the particle is a geodesic of $g_{\alpha\beta}$
How significant the difference between the trajectories is depends on the relative masses. For spacecraft, or even planets, the difference is relatively small. For example the orbit of Mercury is correctly described to within experimental error without considering the perturbation to the metric due to Mercury's mass. However the spacetime geometric of merging black holes is totally different to a simple Schwarzschild calculation.
A: First note that in the case of block holes it is not even clear what "following a geodesic of the entire manifold" would mean as the interior of the black hole exists outside of the causal past of the rest of the space time. Even for extended bodies (e.g. stars) it is not immediately clear what this means as you would first have to identify a worldline for the body. (In Newtonian gravity this would its center of mass, but a priori it is not clear what to take in GR.)
That being said, what can be shown (see https://arxiv.org/abs/1405.5077) is that for extended bodies there always exists an "effective metric" and an "effective worldline" such that the world line is that of a test particle with the same multipole moments as the extended body through the effective spacetime. In particular, a spherically symmetric extended body will follow a geodesic of the effective spacetime.
For black holes the same statement has been proven (see https://arxiv.org/abs/1506.06245) perturbatively in the mass ratio of the objects (to any order). 
