Iterative time-domain method for resolving multiple gravitational wave sources in Pulsar Timing Array data

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(10^3)$. 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 $10^2$ to $10^3$ 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.