The post-Newtonian approximation is a small-velocity perturbation scheme, extensively used to model the evolution of slowly moving isolated systems in the weak field regime of general relativity. It is valid in the system near zone, a region that contains the source but is much smaller in size than its typical gravitational wavelength. Outside the matter support, the multipole expansion of the metric is constructed by solving iteratively the Einstein equations, relaxed to harmonic coordinates, with the help of the retarded Green function. The resulting expression for the gravitational field is bound to match the one derived from post-Newtonian calculations in their common domain of validity. This approach is appropriate to determine analytically both the phase and the amplitude of the waveform of an inspiraling binary of compact objects, including the spin and finite-size effects. To build unbiased waveforms suitable for the current ground-based detectors while preparing the data analysis of their successors, one must, most certainly, control the dynamics of those systems beyond the third-and-a-half post-Newtonian (3.5PN) order. In this context, we have achieved the computation of a 4PN Fokker Lagrangian for two point masses, treating the field divergences, notably those arising in the far zone, by means of dimensional regularisation. At this order, there also appears a non-local-in-time contribution to the action, which is due to the interaction of the source with the linear waves coming from its own radiation. This "wave tail" induces number of interesting subtleties.