The search for anisotropy in the gravitational-wave background with pulsar-timing arrays and astrometry

Pulsar-timing arrays (PTAs) are seeking gravitational waves from supermassive-black-hole binaries, and there are prospects to complement these searches with stellar-astrometry measurements. Theorists still disagree, however, as to whether the local gravitational-wave background will be isotropic, as arises if it is the summed contributions from many SMBH binaries, or whether it exhibits the type of anisotropy that arises if the local background is dominated by a handful (or even one) bright source. Here we derive, using bipolar spherical harmonics, the optimal PTA estimators for anisotropy in the GW background and simple estimates of the detectability of this anisotropy. We provide results on the smallest detectable amplitude of a dipole anisotropy (and several other low-order multipole moments) and also the smallest detectable amplitude of a"beam"of gravitational waves. Results are presented as a function of the signal-to-noise with which the GW signal is detected and as a function of the number of pulsars (assuming uniform distribution on the sky and equal sensitivity per pulsar). We provide results first for measurements with a single time-domain window function and then show how the results are augmented with the inclusion of time-domain information. The approach here is intended to be conceptually straightforward and to complement the results of more detailed (but correspondingly less intuitive) modeling of the actual measurements.

Pulsar-timing arrays (PTAs) are seeking gravitational waves from supermassive-black-hole binaries, and there are prospects to complement these searches with stellar-astrometry measurements. Theorists still disagree, however, as to whether the local gravitational-wave background will be isotropic, as arises if it is the summed contributions from many SMBH binaries, or whether it exhibits the type of anisotropy that arises if the local background is dominated by a handful (or even one) bright source. Here we derive, using bipolar spherical harmonics, the optimal PTA estimators for anisotropy in the GW background and simple estimates of the detectability of this anisotropy. We provide results on the smallest detectable amplitude of a dipole anisotropy (and several other low-order multipole moments) and also the smallest detectable amplitude of a "beam" of gravitational waves. Results are presented as a function of the signal-to-noise with which the GW signal is detected and as a function of the number of pulsars (assuming uniform distribution on the sky and equal sensitivity per pulsar). We provide results first for measurements with a single time-domain window function and then show how the results are augmented with the inclusion of time-domain information. The approach here is intended to be conceptually straightforward and to complement the results of more detailed (but correspondingly less intuitive) modeling of the actual measurements.
It is still not understood, though, whether the local GW signal due to SMBH mergers will be the type of stochastic background that arises as the sum of a large number of cosmological sources, or whether it will be dominated by just a handful-or even just one-source [22][23][24][25][26]. Roughly speaking, if there are ∼ N sources contributing to the signal, then the amplitude of anisotropy in the GW background should be ∼ N −1/2 . A first obvious step, after the initial detection of a gravitationalwave signal, will therefore be to seek the anisotropy in the background that may arise from a finite number of sources. Exotic sources might also lead to anisotropy [27].
Prior work [28][29][30] has developed tools to characterize and seek with PTAs anisotropy in the GW background that were then implemented in a null search [31].
This anisotropy was characterized (as it is here also) in terms of an uncorrelated and unpolarized background of gravitational waves with a direction-dependent intensity parametrized in terms of spherical-harmonic expansion of the intensity. Here we re-derive anisotropy-detection tools using mathematical objects (bipolar spherical harmonics; BiPoSHs [32][33][34]) developed for analogous problems in the study of the cosmic microwave background. The analysis here provides some simplifications and insights and also intuitive estimates for the smallest detectable signals. We provide numerical results for the smallest detectable dipole-anisotropy amplitude as a function of the signal-to-noise with which the isotropic signal is detected and as a function of the number of pulsars in the array. We restrict our attention to PTAs but describe how the detectability will be augmented with the inclusion of astrometry.
This paper is organized as follows: In Section II we describe the idealized observables that we model. In Section III we review the standard Hellings-Downs correlation function (and its harmonic-space equivalent, the timing-residual power spectrum) used to detect the GW background. Section IV introduces the bipolar-sphericalharmonic formalism and describes how to infer the Bi-PoSH amplitudes from the observables. Section V describes the model of an uncorrelated anisotropic background we consider here (and considered in Refs. [29,30]) and calculates the BiPoSH coefficients for the model in terms of the model's anisotropy parameters g LM . Section VI derives minimum-variance estimators for the spherical-harmonic coefficients g LM that parametrize the anisotropy and the variances (∆g LM ) 2 with which they can be measured. Section VII evaluates the smallest detectable anisotropy beginning with a dipole and then generalizing to anisotropies of higher-order multipole moment and then the anisotropy due to a beam of uncorrelated unpolarized gravitational waves from a specific direction. Section VIII describes how the previous results, obtained for a single timing-residual map, are generalized to incorporate the multiple maps that may be obtained from time-domain information. We discuss the extension to astrometry and make closing remarks in Section IX.

II. HARMONIC AND REAL-SPACE ANGULAR OBSERVABLES
PTA measurements are characterized by the temporal evolution of the timing residuals and the dependence of the observables as a function of position on the sky. Here we focus primarily on the angular structure. To simplify, we speak here of the "timing residuals" z(n) measured in a PTA as a function of positionn on the sky. These "timing residuals," in a more complete analysis, will be obtained from some convolution of the timing residuals (TRs) with a time-sequence window function (and there may well be a number of such timing residuals that are obtained from convolutions of the full timing-residual data with different time-sequence window functions-more on this in Section VIII). Strictly speaking, therefore, each appearance of a GW power spectrum P h (k) in the expressions below should be re- Any such timing residual z(n) can be expanded in terms of spherical harmonics Y m (n), which constitute a complete orthonormal basis for scalar functions on the two-sphere. The expansion coefficients are obtained from the inverse transform, The sum in Eq. (1) is only over ≥ 2, as the transversetraceless gravitational waves that propagate in general relativity give rise only to timing-residual patterns with ≥ 2. We assume that the timing residuals (convolved with a time-sequence window function) are real, and so z * m = (−1) m z ,−m . 1 Note that specification of z(n) is equivalent to specification of z m and vice versa-they are two different ways to describe the same observables.

III. POWER SPECTRUM AND CORRELATION FUNCTION
The timing residuals z(n,k) arising from a gravitational wave with polarization tensor h ab moving in directionk are given by z(n;k) = n a n b h ab 2(1 +k ·n) .
Strictly speaking, the timing residuals are observed as a function of time, but the angular pattern here is that after those time-domain data have been convolved with a time-domain window function so that the resulting map z(n) is then real.
As discussed in Refs. [21,35] (and below), the rotationally-invariant observed power spectrum C ∝ m |z m | 2 /(2 + 1) for this plane wave is Since Eq. (3) is a scalar and linear in h ab , the timing residuals from any collection of plane waves-i.e., any gravitational-wave signal-will have the power spectrum of Eq. (4). 2 If the timing residuals z(n) arise from a realization of a statistically isotropic gravitational-wave background, then the spherical-harmonic coefficients z m of the observed z(n) map will satisfy where the angle brackets denote and average over all realizations of the gravitational-wave background, and δ and δ mm are Kronecker deltas. Eq. (5) states that if the GW background is statistically isotropic then all of the z m are uncorrelated and that each z m is some number selected from a distribution of variance C . The resulting map, z(n), is then real after convolution with the appropriate time-domain window function. The timing-residual power spectrum is related to the two-point autocorrelation function, C(Θ) = z(n)z(m) n·m=cos Θ = 2 + 1 4π C P (cos Θ); (6) i.e., the product of the timing residuals in two different directions separated by an angle Θ, averaged over all such pairs of directions. The two-point autocorrelation function from a stochastic GW background is the classic Hellings-Downs curve, Again, the two-point autocorrelation function has this form regardless of whether the GW background is statistically isotropic or otherwise. Since the power spectrum C and two-point autocorrelation function C(Θ) do not depend on whether the background is isotropic or otherwise, the natural first step in any effort to detect a GW background is to establish from the data that these are nonzero. Formulas to derive C from (idealized) data are provided below.

IV. BIPOLAR SPHERICAL HARMONICS
There is, however, far more information in a map z(n) (or equivalently, its set of z m ) than that provided by the timing-residual power spectrum and Hellings-Downs correlation. The most general correlation between any two z m s can be written (see, e.g., Ref. [36,37]), where C is the (isotropic) power spectrum introduced above, m m |LM are Clebsch-Gordan coefficients, and the A LM are BiPoSH coefficients. Note that the power spectrum C can be identified as (−1) A 00 / √ 2 + 1. As Eq. (6) indicates, the Hellings-Downs curve C(Θ) considers information obtained only from the angular separation Θ between two directionsn andm, but disregards any information about the specific directionsn andm. This additional information is parametrized with BiPoSHs in terms of BiPoSH coefficients A LM that characterize departures from statistical isotropy. If there is a dipolar power anisotropy (higher flux of GWs from one direction than from the opposite direction), it is characterized by the L = 1 (dipolar) BiPoSHs, and the different M = 0, ±1 components provide the spherical-tensor representation of the dipole. A quadrupolar power asymmetry (e.g., as might arise if there were GWs coming from the ±ẑ direction) are characterized by the L = 2 BiPoSH coefficients, and so forth.

A. Measurement of BiPoSH coefficients
We suppose that the "data" come in the form of a collection of measured values z data m = z m + z noise m each of which has a contribution z m from the signal and another z noise m from measurement noise. We assume that the noise in each z noise m are uncorrelated and that each z m has a variance N zz (which we further assume to beindependent, a good approximation if the timing-residual noises in all pulsars are comparable).
Estimators for the BiPoSH coefficients are then obtained from This estimator has a variance, under the null hypothesis (for even L + + ) [36], where C data = C + N zz is the power spectrum of the map, which includes the signal and the noise. We will see below that we need consider only combinations with even + +L. If the map z(n) is real, then A LM = A LM (for even + + L), and the estimators A LM and A LM are the same. The covariance between any two other different A LM vanishes. By setting L = 0 and identifying C = (−1) A 00 / √ 2 + 1, we recover the power-spectrum estimator, which has a variance Under the null hypothesis, C data = N zz .

V. MODEL AND BIPOSH COEFFICIENTS
We now focus on understanding the , dependence of the BiPoSH coefficients A LM . To do so, we must understand the dependence of the observable z(n) on the gravitational-wave background.

A. Model of anisotropic background
In order to link measurements of the timing residuals to an underlying gravitational wave background, we need a model for the statistics of that background. Although there are an infinitude of ways the background can depart from statistical isotropy, we consider (as did Refs. [29,30]) here those that can be parametrized as where h s ( k) is the amplitude of the gravitational-wave mode with wavevector k and polarization s = +, ×.
With the Dirac delta function in this parametrization, we are still preserving the assumption that different Fourier modes are uncorrelated. We are also assuming that the frequency dependence of the GW background is the same in all directions 3 and that the + and × modes are still equally populated (i.e., that the background is unpolarized). The sum over spherical harmonics allows, however, for the most general angular dependence of the gravitational-wave flux, parametrized by sphericalharmonic coefficients g LM . Here, the gravitational-wave power spectrum is P h (k), and an isotropic background is recovered for g LM → 0 for all L > 0. In this model, the g LM are the spherical harmonic coefficients of the map of total gravitational-wave power.
Since the term in the brackets in Eq. (13) must be positive, the spherical-harmonic coefficients are restricted to be g L0 ≤ 4π/(2L + 1), and a roughly similar bound applies to √ 2 Reg LM and √ 2 Img LM for M = 0.
B. Resulting timing-residual BiPoSH coefficients (and angular power spectrum) We now calculate the BiPoSH amplitude that arises from a GW background of the form in Eq. (13), based on its imprint, Eq. (3). If the GW direction is taken to bê k =ẑ, then this becomes where h + and h × (both most generally complex) are the amplitudes of the + and × polarizations. This plane wave is described by spherical-harmonic coefficients, where we defined From this result, we can construct the spherical-harmonic coefficients for a plane wave in any other direction. To do so, we write where D ( ) mm (φ k , θ k , 0) are the Wigner rotation func-tions. 4 We thus infer, given Eq. (15), which restricts the m sum to ±2, that a gravitational wave moving in thek direction imprints a pulsar-timing-residual pattern described by spherical-harmonic coefficients, where h + and h × are the GW amplitudes for this wave, and the arguments of the rotation matrices are (φ k , θ k , 0). Given this result, we can now calculate the BiPoSH coefficients for a direction-dependent power spectrum of the form given in Eq. (13). We start by noting that a given gravitational-wave pattern is described by a set of amplitudes h + ( k) and h × ( k) for each possible wavevector k. The spherical-harmonic coefficients induced by this gravitational-wave pattern are The correlation between any two spherical-harmonic coefficients is therefore After performing the integral over directionsk we find an expression for z m z * m of the form Eq. (8) with and where in terms of Wigner-3j symbols, and we defined I ≡ [4π/(2π) 3 ] k 2 dk P h (k). The two terms in Eq. (20) cancel if + + L is odd, and so A LM is nonzero only for even + + L. There are two interesting features of Eq. (22). First, the dependence appears only in the factors z z H L ; as we will see below, this will allow us to write an optimal estimator for the anisotropy coefficients g LM . Second, the power spectrum and BiPoSH coefficients both depend in the same way on the power spectrum P h (k).

VI. MINIMUM-VARIANCE ESTIMATORS OF ANISOTROPY
A. Isotropic signal-to-noise Before evaluating the smallest detectable anisotropy, we write the power spectrum in terms of the signal-tonoise ratio (SNR) with which the isotropic signal is detected; this will be useful below. To do so, we recall that the variance with which any given C can be measured is [2/(2 + 1)](N zz ) 2 . The signal-to-noise ratio SNR from a measurement that accesses multipole moments up to max is then given by and by using Eq. (21), we find I 2 [1.17(SNR)N zz ] 2 in the limit max → ∞, which turns out to be remarkably accurate for any finite max ≥ 3 given the very rapid decay of the summand with . We can thus fix the gravitationalwave amplitude I in terms of the signal-to-noise ratio with which the isotropic signal has been established, and the noise term, which will cancel in the estimator variance calculation below. 5

B. BiPoSH estimators and variance
The observables that we seek to obtain from the data are the anisotropy amplitudes g LM . Each estimator A LM provides an estimator, for g LM . The variance of each of these estimators is (for L + + even), 5 Eq. (24) also indicates that 93% of the total signal to noise in the detection of the GW signal comes from the quadrupole.
We then combine all the estimators (g LM ) with inversevariance weighting to obtain the minimum-variance estimator, Note that the sums here are only over pairs that have even + + L, | − | ≤ L ≤ + , and , ≤ max . Given the reality of z(n), we sum only over ≥ to avoid double counting the contributions from A LM and A LM . The variance (∆g LM ) 2 with which g LM can be measured is the inverse of the denominator; i.e., .
In this equation, the sum is now over all -pairs with | − | ≤ L ≤ + , + + L even, and , ≤ max , which we obtain by using 1 + δ = 2 for = and then including both > and < and dividing by 2 [37]. We can then use to obtain for the SNR → ∞ limit, The smallest g LM that can be distinguished at the ∼ 3σ level from the null hypothesis g LM = 0 is then g LM,min 3(∆g LM ).

A. Results for dipole anisotropy
We now illustrate with the dipole L = 1. To do so, we note that the only nonvanishing -pairs are those with = ±1. We then choose to take = +1 and multiply by two and use (H L=1 , We take the sum on up to max −1 so that the maximum = + 1 corresponds to the largest multipole max that is measured. We then set max N on the sky; the sensitivity to higher-modes will be exponentially reduced. We then plot in Fig. 1 the smallest detectable (at the 3σ confidence level) dipole-anisotropy amplitude g 1M as a function of the signal-to-noise ratio SNR with which the isotropic signal is detected and for different numbers of pulsars. The results can be understood by noting that Eq. (31) becomes, in the SNR → ∞ limit and in the limit max 1, However, this asymptotic limit is reached only for very large SNR, given the very rapid decrease of C /N zz with . In more physical terms, the anisotropy is inferred through correlations between spherical-harmonic modes of different , and so individual modes of higher must be measured with high signal-to-noise. The steep dropoff of C with (each of the seven = 3 moments has a signalto-noise that is smaller by a factor of 5 than that for each quadrupole moment) requires that the isotropic signal (which is very heavily dominated by the quadrupole) be detected with very high significance. In the low-SNR limit, Eq. (31) is approximated, given the rapid decrease of the summand with in this low-SNR limit, this result is obtained for any max ≥ 3.
In practice, this SNR → 0 limit is somewhat academic (and optimistic), as the factor C =2 /N zz in the denominator of the summand in Eq. (31) is already 1.8 for SNR = 3. Thus, the numerical result is a bit larger, even for SNR = 3, than indicated by Eq. (33). The numerical results in Fig. 1 then indicate that the scaling with higher SNR is more like (SNR) −1/2 rather than (SNR) −1 at higher SNR, and that the SNR → ∞ limit (the "pulsarnumber-limited regime") is achieved for SNR 1000. This can be understood by noting that, for example, the C /N zz in the denominator of the summand in Eq. (31) does not reach unity until the signal-to-noise ratio grows, for = 4, to SNR 60, and for = 8, to SNR 1500. This then shows that the benefit of 16 ( 64) pulsars for this particular measurement is limited until the GW signal is detected at SNR 60 ( 1500).
If we wanted first to simply establish the existence of a dipole, without specifying its direction, then our observable would be the overall dipole amplitude, where we have included the factor of √ 4π so that d ≤ 1. The SNR with which this can be established is then √ 3 times that with which any individual g 1M can be measured, and so the smallest detectable (at 3σ) dipole has an amplitude again noting that the SNR → 0 limit is likely overly optimistic for SNR 3 and the SNR → ∞ limit is valid for max 1.

B. Higher L modes
The results for higher L of the smallest detectable g LM are easily obtained by numerical evaluation of Eq. (30) and shown for different max and (SNR) in Fig. 2. The qualitative dependence of the results are similar to those for g 1M,min , although the sensitivity to higher-L anisotropies is reduced a bit (e.g., by about 50% for L = 5) relative to the dipole sensitivity.

C. A gravitational-wave beam
Suppose that a gravitational-wave signal has been detected and that we wish to determine the fraction of the local gravitational-wave energy density coming from a specific direction. To be more precise suppose that we model the gravitational-wave signal as an isotropic uncorrelated background plus a flux of gravitational waves all coming from some specific direction (e.g., the direction of some specific SMBH binary candidate), which we take to be in theẑ direction, that makes up a fraction f of the local GW energy density. This situation is described by anisotropy coefficients g LM = √ 4πi L f δ M 0 .   Fig. 1), and the smallest detectable quadrupoleanisotropy amplitude, shown together with the smallest detectable (again, at 3σ) beam amplitude f obtained with measurements of gL0 up to L=8 is fmin 0.28 in the high SNR limit.
The minimum-variance estimator for the amplitude f is then obtained by summing the minimum-variance estimators for g L0 (scaled by √ 4πi L ), with inverse-variance weighting. In Fig. 3, we plot the smallest f using the results above for g LM,min for L ≤ 8, detectable with measurements of g LM up to L = 8, as a function of SNR from a single map and find it approaches f min 0.28 in the SNR 1000 regime.

VIII. MULTIPLE MAPS
So far, we have assumed that there is a single timingresidual map z(n) obtained by convolving the timedomain data with a single window function. Suppose, however, that the time-domain data are convolved with n w different time-domain window functions that have negligible overlap in frequency space (or in phase). For example, if we were to have measurements performed, every two weeks for ∼ 10 years, yielding ∼ 250 measurements for each pulsar, the time-domain window functions could be taken to be the ∼ 250 different time-domain Fourier modes. In this case, we will have n w ∼ 250 statistically independent timing-residual maps z i (n), with i = 1, 2, . . . , n w . If the Hellings-Downs power spectrum is detected with signal-to-noise ratio (SNR) i in each individual map i, then the signal-to-noise ratio (squared) for the entire experiment, after co-adding all the information, will be (SNR) 2 = i (SNR) 2 i . The optimal estimator for any given g LM is then obtained by adding (with inverse-variance weighting) the estimators g i LM from each map i; i.e., we augment Eq. (28) with an additional sum over i and replace the SNR, the power spectrum C , and noise power spectrum N zz by those-(SNR) i , C i , and N zz i -associated with the ith map: .
(36) The SNR, power spectrum, and noise power spectrum for each map are related by In the SNR → 0 limit, these replacements result (given i (SNR) 2 i = (SNR) 2 ) in the same anisotropy sensitivity as inferred for a single map in this limit in Eq. (33). If, however, the signal-to-noise ratio SNR i in some number n high of maps is high enough (e.g., (SNR) i 60 for N p = 16 or (SNR) i 1500 for N p = 64) that pulsar-numberlimited regime is reached in each individual map, then the sensitivity to anisotropy can be improved by a factor √ n high relative to that, Eq. (32), as shown in Fig. 4. It must be kept in mind that the improvement shown in Fig. 4 possible with additional maps is achieved only if the total SNR is split evenly among all of these many maps.
The remaining question, then, is how the total signalto-noise is distributed among the maps. In the bestcase scenario, it will be distributed equally among the n w maps. If so, then sensitivity to anisotropy could be improved, in principle, by the factor √ n w over that in Eq. (32), as shown in Fig. 4. This improvement would require, however, that the total signal-to-noise be ∼ √ n w larger than that (SNR 60 for N p = 16 and SNR 1500 for N p 64) for a single map. Given the likely (given the most promising astrophysical scenarios) decrease of the signal with GW frequency, however, the signal-to-noise will probably be dominated by a small subset of the maps (those at the lowest frequencies). If so, then n high may be far smaller than n w , and the sensitivity, from multiple maps, to anisotropy will be only marginally improved over the single-map estimate in Eq. (32).

IX. DISCUSSION
We have discussed the search for anisotropy in a PTA-detected gravitational-wave signal in terms of bipolar spherical harmonics for idealized measurements parametrized in terms of the number of pulsars (assumed to be uniformly distributed on the sky) and the signalto-noise ratio (SNR) with which the isotropic signal is established. We focussed our attention first on the case of a single timing-residual map z(n) (obtained from the convolution of the data with a single time-domain window function) and then discussed the generalization to multiple maps (which take into account more of the timedomain information).
We considered a search for anisotropy in an uncorrelated and unpolarized GW background in which the anisotropy is independent of GW frequency. In this case, the anisotropy is parametrized entirely in terms of spherical-harmonic coefficients g LM . We derived the optimal estimators for these g LM for idealized measurements in which N p pulsars are distributed roughly uniformly on the sky and the same timing-residual noise in each pulsar. We then obtain the variance with which each g LM can be determined; this variance is expressed in terms of the signal-to-noise with which the isotropic signal is detected and in terms of the number of pulsars.
The main qualitative upshot of the analysis is that the isotropic signal will have to be established very well before there is any possibility to detect anisotropy. The reason stems from the the fact that the anisotropy is obtained (for odd L) through cross-correlation of sphericalharmonic modes z m of the timing-residual map of different and from the fact (inferred from Eq. (24)) that 94% of the (SNR) 2 for the isotropic signal comes from = 2, with only 6% coming from higher modes. Our numerical results in Fig. 1 show that with a single map it will require the isotropic signal to be established with SNR 1000 before even the maximal dipole anisotropy can be distinguished, at the 3σ level, from a statistically isotropic background. This would, moreover, require 60 pulsars spread uniformly over the sky. The sensitivity to a dipole signal can be improved with more pulsars and/or (as Fig. 4 shows) with multiple maps, constructed with different statistically-independent timedomain window functions. This latter improvement can be achieved requires, however, only if the signal-to-noise is spread evenly among these different maps. Fig. 2 indicate the additional challenge facing a search for higherorder anisotropy.
When discussing the prospects to detect anisotropy, we must be careful to state clearly the question we are trying to answer. Here we have focused on the sensitivity to departures from statistical isotropy parametrized in terms of spherical-harmonic coefficients g LM , under the null hypothesis of a statistically-isotropic background. This sensitivity is limited not only by measurement noise, but also by cosmic variance. In our null hypothesis of a statistically isotropic signal, the spherical-harmonic coefficients z m for the map are selected, in the limit of no-noise measurements, from a distribution with variance C . Departures from statistical isotropy show up, roughly speaking, in terms of disparities between the amplitudes of the different m modes for a given . The conclusion of our analysis is that this is difficult to establish given the variance C under the null hypothesis.
A measurement that is consistent with statistical isotropy may still well exhibit some evidence that the local GW background is a realization that exhibits anisotropy. Suppose, for example, that we had precise measurement of the five timing-residual quadrupole moments z 2m and found that the z 22 and z 2,−2 components were significantly larger than the other three quadrupole moments. Our calculation [obtained by evaluating Eq. (30) with only the = = 2 term] indicates that it would be impossible to infer from this measurement any departure from statistical isotropy. Still, if such a result were observed, it would indicate that the local GW signal is coming primarily from the ±ẑ direction. If there was indeed a strong candidate GW source (e.g., a SMBH-binary candidate) in theẑ direction, then this observation would provide some evidence that the GW signal was coming predominantly from that source.
Our initial calculations explored the detectability of anisotropy from a single timing-residual map obtained by convolving the data with a single time-domain window function. If, however, multiple maps that explore different GW frequency ranges can be obtained, then there are prospects to co-add the anisotropy estimators from those maps to improve upon the pulsar-number limit that arises from a single map. Significant improvement in this way requires, however, that the z m are measured with high SNR in multiple maps.
We also considered the prospects to measure the fraction f of the local GW intensity that comes from a given direction. We conclude also that measurement of f will be similarly challenging: For example, we found a value f min 0.28 for the smallest detectable fraction for a survey with 64 pulsars with SNR 1000 for the isotropic signal. This calculation leaves out ingredients (e.g., timing, polarization, and source evolution) that, if included in the analysis, might improve the ability to localize a point source. Still, the rough conclusions and scalings of the sensitivities with SNR and pulsar number should translate to those for a more complete pointsource search.
We also note, for possible comparison with previous work in configuration space [28][29][30], that a sky described by BiPoSHs A LM has a two-point correlation function,  (39) are the bipolar spherical harmonics (BipoSHs). These Bi-PoSHs constitute a complete orthonormal basis for functions ofn andm in terms of total-angular-momentum states labeled by quantum numbers L and M composed of angular momentum states with l and l ; they are an alternative to the outer product of the {l, m} and {l , m } bases.
The analysis presented here should be straightforwardly generalized to astrometric GW searches. As shown in Ref. [21], the E-mode map from an astrometry survey provides the same information as a timing-residual map. Therefore, everything said here about a timingresidual map can be applied equally to the E-mode map. The higher density of stellar astrometric sources on the sky may ultimately allow higher max but the advantage of this higher max for anisotropy searches can be capitalized upon only with a sufficiently high SNR. The B modes in the astrometry map can provide additional statistically independent information and, when combined with the E modes and/or timing residuals, can conceiv-ably improve the sensitivity to anisotropy by a factor of √ 2.
The numerical results we find for the sensitivity to anisotropy are most likely optimistic, given the idealizations assumed here. Uneven distribution of pulsars on the sky and/or pulsar-to-pulsar variations in the timingresidual noises will degrade the sensitivity. There is another, more subtle, caveat: The estimator in Eq. (27), and the expression, Eq. (28), for its variance, are derived under the assmption that the different A LM are statistically independent. Although the covariance between any two different A LM vanishes (except for those with ↔ ), they are not statistically independent. The variance will thus most generally be a bit larger, and the sensitivity to anisotropy a bit degraded. We anticipate, though, that for the low-L values considered here that this will be a relatively small (perhaps ∼ 10%) effect, although this should be evaluated further with Monte Carlo simulations of isotropic GW signals.
We hope that the approach developed here provides a conceptually straightforward way to understand the search for anisotropy in the GW background and aids in the development of observational/analysis strategies for the PTA search for gravitational waves. It will be interesting in future work to compare the results here to those obtained from detailed simulations of the PTA analysis pipeline, as well as with those inferred from a fully Bayesian approach (see similar applications for the cosmic microwave background [38, 39], for example). It will also be interesing to extend the analysis here to seek anisotropies in the polarization of the GW background, as parametrized, for example, by GW Stokes parameters [40], or anisotropies in the frequency dependence of the GW background.
During the preparation of this work, we learned of related work [41] in preparation that addresses prospects to detect anisotropy with a focus on developing a formalism to produce maps of the gravitational wave background from pulsar timing array measurements. We plan to follow up with detailed comparison of that work and the formalism described here.