XMM-Newton observations of PSR B0656+14, PSR B1055-52 and Geminga have substantially increased the statistics available for these three isolated neutron stars, so apparently similar to deserve the nickname of "Three Musketeers" (Becker & Truemper, 1997). Here we shall take advantage of the EPIC statistics to perform phase resolved spectroscopy for all three objects. The phase-averaged spectrum of the three musketeers is best described by a three component model. This includes two blackbody components, a cooler one, possibly originating from the bulk of the star surface, and a hotter one, coming from a smaller portion of the star surface (a "hot spot"), plus a power law. The relative contributions of the three components are seen to vary as a function of phase, as the stars' rotation bring into view different emitting regions. The hot spots, which have very different apparent dimensions (in spite of the similarity of the three neutron stars polar cap radii) are responsible for the bulk of the phase variation. The amplitude of the observed phase modulation is also markedly different for the three sources. Another striking aspect of our phase-resolved phenomenology is the apparent lack of any common phase alignement between the observed modulation patterns for the two blackbody components. They are seen to vary in phase in the case of PSR B1055-52, but in anti-phase in the case of PSR B0656+14. These findings do not support standard and simplistic models of neutron star magnetic field configuration and surface temperature distribution.