The sensitivity of ongoing searches for gravitational wave (GW) sources in
the ultra-low frequency regime (10−9~Hz to 10−7~Hz) using Pulsar
Timing Arrays (PTAs) will continue to increase in the future as more well-timed
pulsars are added to the arrays. It is expected that next-generation radio
telescopes, namely, the Five-hundred-meter Aperture Spherical radio Telescope
(FAST) and the Square Kilometer Array (SKA), will grow the number of well-timed
pulsars to O(103). The higher sensitivity will result in greater distance
reach for GW sources, uncovering multiple resolvable GW sources in addition to
an unresolved population. Data analysis techniques are, therefore, required
that can search for and resolve multiple signals present simultaneously in PTA
data. The multisource resolution problem in PTA data analysis poses a unique
set of challenges such as non-uniformly sampled data, a large number of
so-called pulsar phase parameters, and poor separation of signals in the
Fourier domain due to a small number of cycles in the observed waveforms. We
present a method that can address these challenges and demonstrate its
performance on simulated data from PTAs with 102 to 103 pulsars. The
method estimates and subtracts sources from the data iteratively using multiple
stages of refinement, followed by a cross-validation step to mitigate spurious
identified sources. The performance of the method compares favorably with the
global fit approaches that have been proposed so far.Comment: 15 pages, 8 figure