  We propose a novel methodology to select host galaxy candidates of future
pulsar timing array (PTA) detections of resolved gravitational waves (GWs) from
massive black hole binaries (MBHBs). The method exploits the physical
dependence of the GW amplitude on the MBHB chirp mass and distance to the
observer, together with empirical MBH mass-host galaxy correlations, to rank
potential host galaxies in the mass-redshift plane. This is coupled to a
null-stream based likelihood evaluation of the GW amplitude and sky position in
a Bayesian framework that assigns to each galaxy a probability of hosting the
MBHB generating the GW signal. We test our algorithm on a set of realistic
simulations coupling the likely properties of the first PTA resolved GW signal
to synthetic all-sky galaxy maps. For a foreseeable PTA sky-localization
precision of 100 squared degrees, we find that the GW source is hosted with 50%
(90%) probability within a restricted number of <50 (<500) potential hosts.
These figures are orders of magnitude smaller than the total number of galaxies
within the PTA sky error-box, enabling extensive electromagnetic follow-up
campaigns on a limited number of targets.
