{'fraction_non_alphanumeric': 0.03590591958570093, 'fraction_numerical': 0.013838739840322233, 'mean_word_length': 4.388001860176717, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Insights into searches for anisotropies in the nanohertz gravitational-wave background

Within the next several years pulsar timing arrays (PTAs) are positioned to detect the stochastic gravitational-wave background (GWB) likely produced by the collection of inspiralling super-massive black holes binaries, and potentially constrain some exotic physics. So far most of the pulsar timing data analysis has focused on the monopole of the GWB, assuming it is perfectly isotropic. The natural next step is to search for anisotropies in the GWB. In this paper, we use the recently developed PTA Fisher matrix to gain insights into optimal search strategies for GWB anisotropies. For concreteness, we apply our results to EPTA data, using realistic noise characteristics of its pulsars. We project the detectability of a GWB whose angular dependence is assumed to be a linear combination of predetermined maps, such as spherical harmonics or coarse pixels. We find that the GWB monopole is always statistically correlated with these maps, implying a loss of sensitivity to the monopole when searching simultaneously for anisotropies. We then derive the angular distributions of the GWB intensity to which a PTA is most sensitive, and illustrate how one may use these"principal maps"to approximately reconstruct the angular dependence of the GWB. Since the principal maps are neither perfectly anisotropic nor uncorrelated with the monopole, we also develop a frequentist criterion to specifically search for anisotropies in the GWB without any prior knowledge about their angular distribution. Lastly, we show how to recover existing EPTA results with our Fisher formalism, and clarify their meaning. The tools presented here will be valuable in guiding and optimizing the computationally demanding analyses of pulsar timing data.


I. INTRODUCTION
Pulsar timing arrays (PTAs; see e.g. [1,2] for recent reviews) are poised to detect the nanohertz stochastic gravitational wave background (GWB) produced by the collection of inspiralling super-massive black hole binaries (SMBHB) within the next few years [3,4]. A detection of the GWB will provide a wealth of information about empirical scaling relations between galaxies and SMBH masses, and shed light on the dynamical processes that affect the evolution of SMBHBs after galaxy mergers (e.g., Refs. [5][6][7]). In addition, PTAs are sensitive to a range of low-frequency GWBs produced by more exotic but no less interesting processes, such as cosmic string networks, non-standard inflation, or phase transitions in the early Universe (see, e.g., Refs. [1,8,9]).
Three major PTAs have been operating for over a decade: the North American Nanohertz Observatory for Gravitational Waves (NANOGrav [10]), the European Pulsar Timing Array (EPTA [11]), and the Parkes Pulsar Timing Array (PPTA [12]). They moreover join effort within the International Pulsar Timing Array (IPTA [13]). These PTA collaborations have been setting increasingly stringent limits to the GWB amplitude, and refining their analysis methods as systematics are better understood. In the near future, the Square Kilometer Array (SKA) [14] is expected to detect hundreds of new low-noise millisecond pulsars [15], which will allow the construction a highly sensitive PTA [16]. With such large datasets, it is critical to have simple but robust and versatile forecasting tools, in order to understand the sensitivity of a PTA without having to run computationallydemanding simulations. In a previous paper [17] (hereafter, Paper I), we laid down the bases of a Fisher formalism to assess the sensitivity of a PTA to a GWB with generic frequency, angular, and polarization dependence. In the present work, we explore in detail the angular dependence aspect, using realistic PTA noise properties.
Most PTA data analyses assume that the GWB is isotropic [18][19][20][21][22]. Although this is expected to be a good approximation, some amount of anisotropy should arise from the finite number of loud, nearby SMBHBs, as well as the non-uniform distribution of cosmological sources tracing large-scale structure [7,23]. Other sources of nanohertz GWs such as cosmic string networks [24][25][26] or domain walls [27] may also produce an anisotropic GWB. A handful of studies have looked into characterizing the sensitivity of PTAs to GWB anisotropies. Refs. [23,28,29] computed the response of pulsar timing residual correlations to an anisotropic GWB, parametrized by its spherical-harmonic expansion. Refs. [30,31] developed a Bayesian analysis pipeline and tested it against simulated GWBs. Lastly, Ref. [32] (hereafter TMG15) conducted a search for GWB anisotropies by performing Markov Chain Monte Carlo (MCMC) analyses of EPTA data. In order to keep computational costs manageable, they limited their analysis to a subset of 6 EPTA pulsars determined to be the most sensitive arXiv:2010.13958v1 [gr-qc] 26 Oct 2020 for continuous gravitational-wave searches [33]. TMG15 found this dataset was not informative enough to constrain anisotropies, as their upper limits saturate the physical prior which imposes a positive GWB intensity.
The main thrust of this paper is to build on the Fisher formalism derived in Paper I, in order to assess the information content of a real PTA, and in particular its sensitivity to the angular distribution of the GWB intensity. We start by laying out the basic formalism in Section II, where we derive the weak-signal-limit Fisher matrix for a GWB with a fixed frequency dependence, but arbitrary angular dependence. To construct this Fisher matrix for a specific PTA requires an accurate characterization of the noise properties of each pulsar. To provide a specific application of our approach as well as make contact with the analysis in TMG15, we focus on the EPTA in this paper. In Section III we describe how we infer the characteristic noise strains of the EPTA pulsars using publicly available timing data, and the public code Hasasia [34,35]. As a first application of the Fisher formalism, we show how one can rank-order pulsar pairs -rather than individual pulsars -by their contribution to the signalto-noise ratio (SNR) of an isotropic GWB. In particular, in Fig. 3 we identify the best 44 EPTA pulsar pairs (out of 861 available) that contribute 90% of the SNR 2 .
Furthermore we show how to use the Fisher matrix to forecast the detectability of GWB anisotropies using standard bases in Section IV. Specifically, we first consider a toy model of a stochastic "hot spot" in a known direction on top of an isotropic background, then forecast the sensitivity of the EPTA to the GWB amplitude in a finite number of coarse pixels, as well as to its sphericalharmonic coefficients. Importantly, we show that a PTA with N pair pulsar pairs can constrain the projections of the GWB on at most N pair independent basis functions. For generic choices of bases (e.g. coarse pixels or spherical harmonics), these projections are, in general, statistically correlated with one another and with the amplitude of the GWB monopole. This implies that the sensitivity to each projection and to the monopole is systematically degraded as one expands the set of basis functions -e.g. by increasing the cutoff max in a spherical-harmonic expansion. In other words, the sensitivity to various GWB anisotropies strongly depends on the priors.
In Section V, we explicitly derive the N pair statistically independent functions of direction that a PTA is most sensitive to. These "principal maps" are simply the eigenmaps of the Fisher matrix, and provide a modelindependent and PTA-specific basis on which to decompose the GWB. They allow one to search under the lamppost provided by a given PTA. As an example application, we discuss how one may use principal maps to attempt to reconstruct the angular dependence of the GWB. We find, however, that a useful reconstruction likely requires a very large amplitude of the GWB, at least for the EPTA data we used. In Section VI, we propose an alternative approach to search for unknown anisotropies, in a modelindependent fashion, using a frequentist criterion.
The hurried reader may want to skip ahead to Section VII, where we apply our formalism to the subset of 6 EPTA pulsars used in TMG15. We derive sensitivity estimates for the low-order spherical-harmonic amplitudes in good agreement with their upper limits. We also clarify the meaning of their pixel-by-pixel upper-limit map, as well as their limits on high-order spherical harmonic amplitudes, for which the number of independent parameters is larger than N pair . We prove that the parameters they constrain are not the pixel-by-pixel or sphericalharmonic amplitudes of the true GWB, but rather of its projection on the N pair -dimensional space of sky maps observable by the PTA. As a consequence, the derived limits get worse with more pulsars, and cannot be interpreted as true upper limits. We sharply demonstrate our point by closely reproducing their results using the Fisher matrix and its principal maps. This comparison with existing results ought to provide a clear demonstration of the usefulness of our Fisher formalism to forecasting the sensitivity of real PTAs to anisotropies. We discuss our results and possible extensions in Section VIII. While we apply our formalism to the EPTA in this paper, we emphasize that our approach can be easily applied to any other PTA given the pulsar timing and parameter files. To this end, we will make a version of the Python script used in this paper public.

II. BASIC CONVENTIONS AND FORMALISM
A. Map notation and dot product in map space Throughout this paper, we denote byΩ a direction in the sky (i.e. a unit-length vector). We refer to functions M :Ω → M (Ω) as "maps" and denote them with bold symbols when referring to the function itself. We denote by 1 :Ω → 1 the monopole map. Lastly, we define the dot product of two maps M 1 , M 2 as follows: With this convention, the monopole map has unit norm, and the average of a map M over directions is just 1·M .

B. Stochastic gravitational-wave background
We denote by h ab (f,Ω) the GW strain at frequency f and with propagation directionΩ. We assume that it is a stationary Gaussian random field (as would be the case if it is generated by a large number of uncorrelated sources). It is therefore entirely determined by its rank-4 power spectrum, for which we give a complete geometric description in Paper I. Here we focus on the total intensity map I(f ) of the GWB, i.e. we assume that the strain power spectrum takes the form where the purely geometric (frequency-independent) rank-4 tensor I abcd (Ω) is given by We refer the reader to Paper I for the correspondence of our frame-independent expressions with the framedependent notation usually adopted in the literature, as well as the generalization of Eq.
(2) to a linearly and circularly polarized GWB. The GWB intensity is related to the characteristic strain h c (f ) through Equivalently, we may write where P (f ) is the (frequency-dependent) angular distribution of GWB intensity, with unit angle average.

C. Pulsar timing residuals correlations
Consider a pulsar p, in the directionp. The times of arrival (TOAs) of its pulses are a combination of a deterministic term due to intrinsic processes (such as the spin-down of the pulsar or variations in the position of the pulsar on the sky), a pulsar-specific intrinsic noise term N p , uncorrelated between different pulsars (due to measurement noise, variations of the dispersion measure, as well as any intrinsic stochastic processes, such as pulsar 'jitter'), and additional stochastic contributions correlated among pulsars. In this paper, we assume that the latter are exclusively generated by the GWB, and neglect additional sources of correlated noise, such as global clock and ephemeris errors [19,36,37]. After fitting for a deterministic timing model, the pulsar time residual R p is thus We define the (one-sided) cross-power spectrum R pq (f ) of the total time residual at pulsars p, q as follows: and similarly define R GW pq (f ) as the cross-power spectrum of the GW-induced timing residuals R GW p . The crosspower spectrum of the intrinsic noise is defined as In Paper I we derived the following simple compact form for the linear relation between R GW pq (f ) and I(f ): where the pairwise timing response function, γpq, is the map with values Using Eq. (6), we may rewrite Eq. (10) as where is the frequency-dependent "overlap reduction function". If the GWB is a pure monopole, the timing residual cross correlation is is proportional to the Hellings and Downs function [38], and only only depends on the dot productp ·q through the function Specifically, this function is given by, for µ ≡p ·q, Note that we do not impose any normalization on the Hellings and Downs function H(µ): it is simply defined as the response to an isotropic GWB intensity with unit amplitude, through Eq. (14). In particular, H(1) = 4/3.

D. Gaussian likelihood for the GWB intensity
In what follows we use capital indices to label unique pairs of distinct pulsars, e.g. I ≡ (p, q) = (q, p) for p = q. We then define R I ≡ R pq , γ I ≡ γpq, and H I ≡ H(p ·q).
Let us define R I (f ) to be unbiased estimators of the timing residual cross-power spectra at a given frequency f . Assuming that the only source of correlated timing residuals for distinct pulsars is the GWB, these estimators have mean value R I = R GW I . In practice, one cannot directly construct finite-variance estimators of the cross-power spectra, but only of finite bandpowers, from which one may then extract R I (f ) for sufficiently narrow bands, see Paper I for details. Mathematically, the covariance matrix of the estimators C IJ ≡ cov R I (f ), R J (f ) is inversely proportional to the bandwidth ∆f used to estimate the bandpowers: we denote it as where both C IJ and C IJ depend on frequency. Estimators of R I at frequencies separated by more than ∼ 1/T are uncorrelated, where T is the characteristic observation time. If the estimators are constructed from a sufficiently large number of uncorrelated data samples, we may approximate their probability distribution P as the product of uncorrelated multivariate Gaussians: where the second approximation holds provided the frequency bands are narrow, ∆f f . As is standard in Bayesian data analysis, one may interpret this probability distribution as the likelihood L of the signal I(f ) given the data (assuming flat priors): It is useful to write this probability explicitly as a Gaussian distribution for I. To do so, we first define the maps γ * I , as the unique linear combinations of the N pair pairwise timing response functions that satisfy γ * I · γ J = δ IJ . These are the dual maps (not to be confused with complex conjugates) of the pairwise timing response functions. We then define the following estimator of the intensity map: We may then rewrite the likelihood (18) in the form where G(f,Ω,Ω ) is a continuous inverse covariance "matrix", given by E. Weak-signal limit In Paper I we derived the weak-signal limit of C IJ : where T p (f ) is the transmission function, which accounts for the loss of information when fitting a deterministic timing model [39], and T pq ≡ min(T p , T q ) is the minimum observation time of the pulsar pair (p, q) (or more generally, the time of overlap between the two pulsar observations). This limit holds if the GWB-induced timing residual is subdominant to the intrinsic timing noise for each pulsar, and at every frequency, i.e. if In what follows we adopt the compact notation of Ref. [39] and define the characteristic noise strain for pulsar p as The factor of 3/4 arises from the Hellings and Downs normalization convention in Ref. [39]; we chose to use exactly the same definition for h c,p as in that paper to avoid any confusion. With this notation, we find the following expression for the weak-limit intensity inversecovariance: where, for I = (p, q), we have defined F. Application to a GWB with a power-law frequency dependence It is of course easier to search for signals with a specific frequency dependence. In the remainder of this paper, we will specialize to a GWB whose frequency and direction dependence factorize, and with a power-law frequency dependence, of the form where f yr ≡ (1 year) −1 . With this convention, and as can be seen from Eq. (5) the characteristic strain takes While we keep our expressions general, for numerical applications we will specialize to a power law α = 2/3. Indeed, this is what is expected from a population of SMBHBs on circular orbits, shrinking exclusively due to GW emission [40]. The true frequency dependence may differ from this simple form, especially at low-frequencies if the binaries are eccentric or if environmental effects such as stellar hardening and/or gas torques are driving the binaries to merge instead of GWs [1,5,41]. Note that our formalism can readily be applied to an arbitrary frequency dependence.
We now want to rewrite Eq. (20) as a likelihood for A. To do so, we define the dimensionless Fisher "matrix" as well as the inverse-variance-weighted estimator amplitude map The likelihood of the GWB amplitude A is then proportional to Equations (28) and (32) are fundamental to all the forthcoming calculations in this work. The total SNR 2 for an arbitrary GWB with amplitude A is given by We see that for each pulsar pair, the coefficient F I characterizes its overall sensitivity to a power-law GWB, and the pairwise-timing response function γ I characterizes its angular response. Before moving on, let us remark that, given N pair distinct pulsar pairs, there are at most N pair independent detectable components of the GWB. In particular, a PTA is completely blind to the projection of the GWB perpendicular to all the pairwise timing response functions. Specifically, given an amplitude A, we may always decompose it as a piece A || which is a linear combination of the γ I 's, and a piece A ⊥ orthogonal to all of them; the latter is completely undetectable by the PTA. The estimator A generalizes the "optimal statistic" A 2 GWB [19,42] to an anisotropic GWB. It is an unbiased estimator of the observable piece of the signal A || , but of course not of the full map A. We can therefore equally write A ≡ A || .

G. Contribution of anisotropies to timing residual autocorrelations
The GWB also contributes to each pulsar's timing residual autocorrelation (or power spectrum). These autocorrelations differ from the cross correlations in an essential aspect: one does expect a non-zero autocorrelation even in the absence of a GWB, due to intrinsic noise. In the analysis of real data, one must therefore search simultaneously for an appropriately parametrized intrinsic noise spectrum, for each pulsar, jointly with a GWB-induced noise common to all pulsars (a "common red noise"). Because of the unknown intrinsic noise, one therefore cannot hope to ever achieve a detection of the GWB with autocorrelations alone [43]. Still, autocorrelations do contain information about the GWB: at the very least, they can be used to set upper limits on its amplitude.
Here we expand on a point that we made in Paper I: a pulsar's autocorrelation term is sensitive to the GWB anisotropy. In other words, GWB anisotropies contribute a common red noise. While this fact can be readily inferred from the expressions for the overlap reduction functions given in Ref. [23], it is unclear whether this important point has been accounted for in previous searches for anisotropic GWB. Moreover, here we show that only the GWB monopole, dipole and quadrupole contribute to pulsars' autocorrelations.
From Eqs. (10) and (27), the GWB contribution to the single-pulsar timing residual power spectrum is proportional to γpp · A. As we saw in Paper I, the single-pulsar timing response function γpp(Ω) is a linear combination of a monopole, dipole projected alongp, and quadrupole twice projected ontop: Let us define the following geometric decomposition of the GWB dimensionless amplitude: where the tensor A ij 2 is symmetric and trace-free (thus has 5 independent components), and ∆A only contains contributions that are octupolar and higher order. The components A 0 , A i 1 and A ij 2 are uniquely defined through They are, respectively, the monopole, dipole and quadrupole components of the GWB. With this decomposition, we obtain Thus, as noted in Paper I, the auto-correlated power spectra constrain specific combinations of the GWB monopole, dipole and quadrupole. Importantly, when allowing for anisotropies, one must self-consistently propagate their impact on the single-pulsar noises, i.e. one must account for "common red processes" which arise from the dipole and quadrupole pieces of the GWB, in addition to its monopole. If autocorrelations are used to set upper limits on a common red process, their dependence on dipole and quadrupole GWB anisotropies systematically degrades the inferred upper limits on the monopole. One of our main objectives here is to explore the ability of PTAs to detect the GWB. While autocorrelations do contain information, they cannot be used to claim a detection [43]. As a consequence we focus on the correlations between different pulsars in the remainder of this paper.

III. DATASET USED IN THIS PAPER
Throughout most of this paper, to provide a concrete example, we will apply our formalism to the full EPTA dataset described in the last data release [44]. The EPTA is composed of the high-precision timing of 42 millisecond pulsars for up to 17 years with an overlap of 21 pulsars with the NANOGrav Nine-Year dataset [45]. We show the locations of the EPTA pulsars relative to the Galactic plane in Fig. 1.
Out of the 42 EPTA pulsars, six were determined to contribute 90% of the SNR 2 in simulated continuous GW searches [33], and were used for searches of anisotropies in the GWB in Ref. [32] (hereafter, TMG15). Specifically, these six pulsars are J0613-0200, J1012+5307, J1600-3053, J1713+0747, J1744-1134, and J1909-3744. Their locations in the sky are shown in Fig. 1, and their characteristic noise strains in Fig. 2. We will use this subset of the EPTA in Section VII to compare our sensitivity estimates to those of TMG15.
A. Evaluation of the characteristic noise strains As described in the previous section, the primary quantity characterizing pulsar noise properties is the characteristic noise strain. Here, we explain how we estimate the characteristic noise strains of all 42 EPTA pulsars.
We model the intrinsic white noise in the standard way: The subset of 6 highlighted pulsars is used when compare our results to those of Ref. [32] in Section VII. The background map is the Stockert and Villa-Elisa 1.4 GHz continuum map [46,47], which situates the positions of the pulsars relative to the Galactic plane.
where t i is the i th time-of-arrival of pulsar p, σ p,i is the uncertainty in the i th timing residual of pulsar p, and the white noise parameters, EFAC p,b and EQUAD p,b , are included for each observing system b (i.e. different telescope and/or backends). EFAC is dimensionless and EQUAD has dimensions of time (we use units of seconds for our analysis). Our priors on these parameters were a flatlinear prior on EFAC ∈ [0.01,10] and a flat-log prior on EQUAD ∈ [10 −8.5 , 10 −3 ] sec. The effects of variations in the dispersion measure (DM) and intrinsic red-noise (RN) were modeled with a power-law spectral density of the form where A RN is independent of radio frequency ν whereas A DM ∝ ν −2 . We use a flat-linear prior on the power-law indices, γ X ∈ [0, 7] and a flat-log prior on the dimensionless amplitudes A X ∈ [10 −20 , 10 −11 ], for both RN and DM. The prior on the power-law indices have been chosen so that they range from white noise (γ X = 0) to the steepest power-law for which the fit to the timing model removes any dependence on the functional form for these spectral densities at low frequencies (f < 1/T obs ) [48]. This range also covers the expected variation in the power-law index due to random walks in phase, period and period derivatives (which give γ RN = 2, 4, 6, respectively) [19]. Note that we do not include a "common red noise" process in our analysis. This is consistent with our weak-signal assumption, and allows us to fit for each pulsar noise properties independently. In order to extract values for these noise parameters we used the Parallel Tempering Markov-Chain Monte Carlo sampler PTMCMCSampler [49], the PTA software Enterprise [50], and the TOA and timing model parameters used in the EPTA data release [51] from the EPTA repository [52]. We set the noise parameters to equal their median values from MCMC chains that contain 10 6 Characteristic noise strains for 6 of the EPTA pulsars, produced using publicly available EPTA data [51] and Hasasia [34,35], see main text for details. The spikes at 1/yr and 2/yr arise from degeneracies with the Earth's orbital motion when fitting the pulsar position and distance. For J1713-0747, the spike at f ≈ 5.4/yr is at the 67.8-day period of this binary system. samples. We have confirmed that our noise parameters match well with the noise analysis described in [53].
After we have estimated the noise parameters for each pulsar we use Hasasia [34,35] to compute the total noise power spectrum σ 2 p (f ), transmission function T p (f ), and characteristic noise strain h c,p (f ), for each one of the EPTA pulsars. As an illustration, we show the characteristic noise strains of 6 of the EPTA pulsars in Fig. 2. We will use these specific 6 pulsars when comparing of our results with those of TMG15 in Section VII.

B. Best pulsar pairs for isotropic GWB searches
In what follows we determine which pulsar pairs contribute most to the total SNR 2 from a cross-correlation analysis. We will see that only a small number of pulsar pairs are needed to get most of the SNR 2 . A similar analysis was performed in Ref. [54] for the Parkes PTA; here we consider the EPTA and provide a much simpler, analytic calculation.
As can be seen from Eq. (33), for a general GWB, each pair I contributes F I (γ I · A) 2 to the SNR 2 . The relative ranking of pulsar pairs thus depends on the angular dependence of the GWB, and there is no universal ranking that holds for an arbitrary GWB. Since we expect the GWB to be predominantly isotropic, it is sensible to rank pulsar pairs under this assumption.
For an isotropic GWB A(Ω) = A 2 h , from Eq. (33) we see that the SNR 2 is given by Using the characteristic noise strains we computed for the EPTA dataset, we find A 95% h ≈ 2.5 × 10 −15 . If we restrict ourselves to the 6 EPTA pulsars used in Refs. [20,32], we find instead A 95% h ≈ 3.4 × 10 −15 . This lies between the 95% upper-limit values found in these references, of 3.0× 10 −15 [20] and 3.9 × 10 −15 [32]. This gives us confidence that our characteristic noise strains provide a realistic description of the data.
From Eq. (41), we see that for an isotropic GWB, each We have ranked the 861 pulsar pairs of the full EPTA dataset in terms of their contributions to the SNR 2 . In Fig. 3 we show the normalized cumulative contributions of the best 44 pulsar pairs, which contributed 90% of the SNR 2 . While this is three times as many pairs as what can be constructed with the 6 pulsars used in Refs. [20,32], these 44 pairs represent 5% of the total number of pairs in the full dataset, and should still constitute a manageable collection of data. Moreover, we find that the 6 pulsars used in Refs. [20,32] only contribute 26% of the total SNR 2 for an isotropic GWB, while the best 15 pairs would amount to 68% of the SNR 2 , as can be seen in Fig. 3.
This result is of significant importance to speed up future analyses of real data. It is well known that accounting for all the correlations between pulsar pairs in a full Bayesian analysis of timing residuals is computationally challenging. Here we see that it suffices to include a small, manageable number of pulsar pairs to recover most of the SNR 2 . Our simple Fisher formalism allows us to efficiently determine which pairs to use for any given dataset. This can be done not only for an isotropic GWB, but also for any assumed angular dependence, if desired.

IV. SENSITIVITY TO GWB ANISOTROPIES USING STANDARD BASES
A. General approach Suppose that we model the GWB amplitude to be a linear combination of known maps M n : where A n are scalar amplitudes. The maps M n could be, for instance, spherical harmonics, or Dirac functions centered at specific directions in the sky. They can be, in general, any set of linearly independent maps, and need not even be orthogonal nor normalized. We can always define the dual maps M * n (not to be confused with J1643-1224 and J1640+2224 complex conjugates), such that M n · M * m = δ mn . Given the estimator map A, we may then define A n ≡ M * n · A. The probability distribution L(A) given in Eq. (32) can then be interpreted as a joint probability distribution for the set of coefficients {A n }, and be recast in the form In other words, the coefficients A n are Gaussiandistributed, with mean A n , and covariance In particular, upon marginalizing over all other parameters, the probability distribution of any single parameter A n is a Gaussian with mean A n and variance (F −1 ) nn . This means that, to be able to claim a 95%-confidence detection of a non-zero coefficient A n , the underlying true parameter needs to satisfy This represents a threshold for detection, i.e. allows us to forecast the sensitivity of a given PTA. Note that (F −1 ) nn ≥ 1/F nn , unless the map M n is statistically uncorrelated with all other maps. Therefore the sensitivity to any given map amplitude A n is always degraded if one simultaneously searches for other amplitudes.
It is important to note that at most N pair map amplitudes can be simultaneously searched for. Indeed, the Fisher matrix F only has rank N pair , as can be seen by writing Eq. (28) formally as a sum of tensor products F = I F I γ I ⊗ γ I . Therefore, the matrix F mn is guaranteed to be singular if it has more than N pair rows and columns. In other words, at most N pair parameters can have finite variance. In practice, one can still set upper limits even with more than N pair parameters, if additional priors are included (e.g. physical priors, or a prior based on the auto-correlation functions of individual pulsars). But in order to be able to claim a detection, at most N pair parameters can be simultaneously searched for.

B. A toy problem: monopole and hot spot with a known direction
As a first application of this formalism, let us consider a simple toy problem, with only two basis maps. Suppose that external considerations lead us to expect a GWB "hot spot" in known directionΩ * , in a addition to a background monopole. Note that we are still considering a stochastic "hot spot", i.e. this cannot be generated by a single binary, but rather must originate from a concentration of many binaries radiating incoherently. Both the hot spot and monopole have an unknown amplitude.
Mathematically, we assume that the GWB amplitude takes the form Note that this is a different split between monopole and hot-spot than adopted in Ref. [29] and Paper I; with the convention (48), the amplitude of the hot spot is unbounded from above. One can interpret A 0 as the amplitude of a background monopole, on top of which the hot spot sits. The total monopole amplitude is then A 0 +A 1 .
In the notation of Section IV A, the maps M 0 , M 1 are The two maps M 0 , M 1 are not orthogonal, and more importantly, are in general not statistically independent, in the sense that M 0 · F · M 1 = 0. Applying the formalism of Section IV A, we must first compute the 2 by 2 inverse-covariance matrix F mn . For any givenΩ * , its components are The covariance matrix of the amplitudes (A 0 , A 1 ) is then the inverse of F . In Fig. 4, we show the minimum background monopole and hot-spot amplitudes detectable at the 95% confidence level, as a function of the direction of the hot spot. We see that the sensitivity to the background monopole amplitude A 95% 0 varies from 2.5 × 10 −15 to 3.3 ×10 −15 , depending on the assumed location of the hot spot. This is to be compared with the estimated EPTA sensitivity to a monopole-only GWB, A 95% h ≈ 2.5 × 10 −15 . We thus see that searching for a monopole simultaneously with another map, in this case a hot spot, systematically degrades the sensitivity to the monopole relative to the case where one assumes the GWB is a pure monopole. This is due to the statistical correlation (i.e. non-zero matrix element F 01 ) between the two maps.
It is important to emphasize that this analysis applies to a hot spot in a known direction. Without prior knowledge of the hot spot direction, it is not possible to set any constraint on its amplitude, as the number of free parameters would then be infinite, thus exceeding N pair . As mentioned, this is mostly a toy problem, as in practice it is highly unlikely that one would expect only one single spot at a time. Next we consider a more realistic setup, where we simultaneously constrain the amplitudes of the GWB in multiple coarse pixels.

C. Detectability of coarsely pixelized anisotropies
As explained earlier, we may simultaneously search for at most N pair independent amplitudes of the GWB. It is therefore not possible to simultaneously constrain the GWB in all pixels in the sky. We can, however, meaningfully analyze the amplitudes of GWB in N pix ≤ N pair coarse pixels. For this analysis, we will use the HEALPix pixelization [55] into N pix = 12N 2 side pixels of equal area 4π/N pix . We define the N pix coarse-pixel basis functions P n such that P n (Ω) = 1 ifΩ is within the n-th pixel, and 0 otherwise, and search for coarse-grained anisotropies of the form To compute the discrete Fisher matrix coefficients F nm , we need to first computed the dot products γpq · P n = 1 Npix γpq n , where M n is the average of map M in the n-th pixel. The latter is obtained by "degrading" the map to a resolution of N pix . The 861 distinct pulsar pairs afforded by the EPTA allow us to use N side as large as 8. We show the 95% sensitivity A 95% h,pix ≡ A 95% n in each coarse pixel in Fig. 5, for N side = 1, 2, 4, corresponding to N pix = 12, 48 and 192, respectively. For N side = 8, we found that the 768×768 coarse-pixel covariance matrix is numerically illconditioned, thus do not consider this case. We see that the pixel-by-pixel sensitivity degrades dramatically as N side is increased (especially in the less well-constrained parts of the sky). Note also that the lower N pix limits cannot be simply recovered by degrading the higher N pix limits, as the order of map-degrading and inverting the Fisher matrix cannot be interchanged. The average characteristic strain squared is given by which has variance translating to a 95% sensitivity A 95% h 2 = 2 var (A 2 h ). For the EPTA, we find A 95% h = (4.6, 5.0, 7.8) × 10 −15 for N side = 1, 2, 4, respectively. This is to be contrasted with A 95% h = 2.5 × 10 −15 if the GWB were known to be strictly isotropic. This highlights once again that the general form assumed for the GWB strongly affects the sensitivity to the monopole.
D. Detectability of spherical-harmonic amplitudes

Setup and physical prior
Following common practice [23,32], we may decompose GWB anisotropies on the basis of real spherical harmonics Y m (Ω), which we normalize such that Y m ·Y m = δ δ mm (i.e. with a factor √ 4π larger than the standard definition of spherical harmonics). Explicitly, we assume the GWB amplitude takes the form Note that this form assumes the rather stringent prior that harmonic coefficients for > max all vanish.
In order for A to represent a physical GWB, it needs to be positive everywhere. To our knowledge, there is no simple analytic expression for this physical prior. Nevertheless, we can derive an approximation as follows. We define the GWB anisotropy ∆A(Ω) ≡ A(Ω) − A 00 . It has the following variance over directions in the sky: where we have defined A simple approximate physical prior follows: we require that the monopole amplitude A 00 is larger than twice the rms anisotropy, implying This implies a conservative physical prior This simple criterion reproduces quite accurately the numerical results of TMG15 (see inset in their Fig. 1). Our conventions map to those of TMG15 through A m = . Note that the monopole coefficient is just the amplitude of the characteristic strain squared, A 00 = A 2 h , and that C 0 = A 2 00 = A 4 h .

Forecasted sensitivity of the EPTA
Let us start by pointing out that for a given maximum = max , the decomposition (56) includes ( max + 1) 2 independent coefficients. It thus follows that one may constrain the spherical-harmonic amplitudes from pulsar pair cross correlations only if max ≤ N pair − 1. For larger values of max , the A m are then only constrained by the physical prior. For the full EPTA, we may therefore use max ≤ 28. For the subset of 6 EPTA pulsars, one is restricted to max ≤ 2.
We compute the covariance of the amplitudes A m as described in Section IV A. Instead of showing the sensitivity to the individual A m coefficients, we estimate the detectability of the rotationally-invariant coefficients C defined in Eq. (58). Their mean and covariance matrix are given by The C are not Gaussian-distributed (especially for small ). Nevertheless, we may define an approximate 95% sensitivity estimate as follows: Note that even for = 0, this gives C 95% 0 = (1 + 2 √ 2)var(A 00 ) ≈ 3.8 var(A 00 ), which is very close to the correct 2-σ sensitivity C 95% 0 = 4 var(A 00 ). We show C 95% in Fig. 6, for several values of max . We see that for any given coefficient C , the sensitivity systematically degrades as max is increased. In particular, as anticipated, the sensitivity to the monopole significantly worsens as max is increased. This is because the monopole is correlated with other spherical harmonics, and as one enlarges the space of functions to be searched over, the uncertainty on the monopole amplitude increases.
Let us remark that the coefficients C are statistically correlated. We explicitly give their correlation coefficients in Tabs. I for max = 2 and 5. We see that these correlations coefficients are in general not small, and depend on the chosen max .
Finally, note that this analysis specifically estimates the minimum amplitudes necessary for a detection, through pulsar timing cross correlations. In addition, the monopole, dipole and quadrupole can be constrained by pulsar autocorrelations, as discussed in Section II G. In the limit that autocorrelations constrain the = 0, 1, 2 harmonic coefficients much more tightly than the cross correlations, the variance of the remaining coefficients should be obtained by inverting the Fisher matrix restricted to these coefficients, leading to a lower noise for the C with ≥ 3. We have checked explicitly that this lowers the noise by no more than ∼ 10 − 20%.
In conclusion, we have demonstrated that our Fisher formalism allows to forecast the sensitivity of a PTA to the spherical-harmonic amplitudes of the GWB, for a given cutoff max . We have applied this specifically for the EPTA, and shown that the minimum detectable amplitudes are systematically larger than allowed by the physical prior, given current upper limits on the monopole. This seems to indicate that spherical harmonics are a suboptimal choice of basis for anisotropy searches. We will specifically compare our results with those of TMG15 in Section VII. , as a function of the cutoff max, beyond which the coefficients are assumed to strictly vanish. The physical prior is computed for a monopole upper limit A h ≤ 4 × 10 −15 [32], and is systematically lower than the minimum detectable anisotropy amplitudes.

A. Motivation and formalism
As highlighted in the previous section, one of the challenges when using standard bases to decompose the GWB intensity is that the amplitudes of the basis maps are statistically correlated. Moreover, their covariance matrix depends on the number of maps considered. As a result, one cannot easily set model-independent limits on the map amplitudes.
If one is completely agnostic regarding the angular dependence of the GWB, it is best to search under the lamppost, i.e. look for the amplitudes of maps which are best-measured by a given PTA, and uncorrelated with one another. To do so, we construct the N pair unit-norm principal maps {M n }. They are defined to extremize SNR 2 = M n · F · M n , under the normalization constraint M n · M n = 1. Dimensionless correlation coefficients cov(C , C )/ var(C )var(C ), for max = 2 and max = 5, using the EPTA pulsars. The ellipses indicate that the tables are symmetric. This illustrates that the C 's are statistically correlated and that for a given pair , , the correlation coefficient depends on the assumed max.
optimization problem are simply the eigenmaps of the Fisher matrix, i.e. maps that satisfy where we denoted the corresponding eigenvalues by 1/Σ 2 n . Since F is symmetric, its eigenmaps are orthogonal to one another and, importantly, uncorrelated with one another: We may rewrite the Fisher matrix as a direct sum of tensor products of the principal maps: Any GWB amplitude A(Ω) can always be written as a linear combination of the N pair eigenmaps and a piece orthogonal to all of them: In that case, the total SNR 2 is given by From Eq. (69), we see that the minimum amplitude of the eigenmap M n detectable with SNR = 2 is A n = 2Σ n .
The eigenmaps for an orthonormal basis of the N pairdimensional space of observable maps for a given PTA. Maps A ⊥ which are orthogonal to this space are completely unobservable, i.e. have infinite noise.
To compute the eigenmaps in practice, we first write them in the form of linear combinations of the pairwise timing response functions: The eigenvalue problem (64) is then equivalent to the finite-dimensional eigenvalue-problem where the symmetric N pair × N pair matrix F IJ has elements We show the first 9 eigenmaps of the EPTA in Fig. 7. These maps look very different from any known analytic basis functions (such as spherical harmonics): indeed, they reflect the properties of a particular PTA. We show the eigenvalues (translated into a detectability threshold in characteristic strain) of the first 20 eigenmaps in Fig. 8. Note that beyond the first ∼ 20 eigenmaps, the noise eigenvalues Σ n start increasing exponentially with n.
One can use the principal maps as a basis to decompose the observable GWB, and follow the same formalism as in Section IV. The advantage of using this basis is that the map amplitudes are statistically uncorrelated. As a consequence, upper limits obtained on any given map amplitude do not depend on the cutoff. In particular, one could start by searching for the amplitude of the first principal map alone, then search simultaneously for the amplitudes of the first two maps, etc..., and not change any of the sensitivity estimates.
In principle, given N pair estimated timing residual cross-spectra R I (f ), one can build an estimate of the underlying GWB intensity from Eqs. (30)- (31). However, in the weak-signal limit, all the estimators A I are individually noise-dominated. Collectively, however, a PTA may still measure some independent pieces of information with high signal-to-noise ratio: these are precisely the amplitudes of the lowest-noise principal maps, A n ≡ M n · A. Concretely, using this expression with Eq. (31), we see that the A n are linear combinations of the timing residual cross-spectra R I (f ), appropriately integrated over frequencies.
We can then use the principal maps to attempt to "reconstruct" the GWB angular dependence. Provided some of the individual principal map amplitudes A n are measured with sufficiently high individual SNR n ≡ A n /Σ n (say, SNR n > 3), we can define the reconstructed map This procedure is analogous to the production of "dirty map" in radio interferometry, i.e. including the contributions of observed "visibilities", and setting the remaining (unobserved) piece to zero. Unlike interferometry, however, the ability to reconstruct a map is a strong function of its amplitude. Indeed, for a single interferometer, different visibilities typically have comparable noise, and as a consequence all contribute to the dirty map once any one of them is detected with sufficiently high SNR. In contrast, the noise of principal maps of a PTA increases steeply with the principal map number, see Fig. 8 for the principal maps of the EPTA. Note that this property is not a result of unequal pulsar noises: even for an array of equal pulsars, densely and isotropically distributed on the sky, the eigenvalues of the Fisher matrix are steep function of principal map index (see Fig. 4 in Paper I).
We illustrate the map reconstruction technique with the EPTA in Fig. 11. In the top row, we show the reconstructed maps obtained if the underlying GWB is a pure monopole, with amplitudes A h = 10 −14 , √ 3 × 10 −14 , 3 × 10 −14 , from left to right. The bottom row shows the reconstructed maps obtained if the underlying GWB is proportional to the map shown in Fig. 10, with the same monopole amplitude. We see that even with these very large amplitudes, the reconstructed maps have little resemblance with the underlying GWB intensity distribution. Note, also, that the overall amplitude of these maps (corresponding to total SNR of 30, 100 and 300, respectively, from left to right), is so large that they are most likely inconsistent with the weak-signal limit. The principal maps provide a useful basis to decompose the GWB if one is fully agnostic about its angular distribution. In practice, however, this is not the case: one does expect a significant monopole component in the GWB. In other words, physical considerations isolate a preferred map, which should always be included in GWB searches. This map is in general different from the principal maps, which are entirely based on the noise and geometric properties of the PTA, rather than external physical considerations. Fig. 9 shows the monopole components of the first 20 EPTA principal maps, as well as their correlation coef- ficient with the monopole. Even though the monopole has the largest projection on the first principal map, we see that it still has significant projections on and correlations with some of the higher-order eigenmaps. In other words, higher-order eigenmaps are (i) not anisotropic and (ii) statistically correlated with the monopole.
In order to alleviate issue (i), one could try to construct a set of "principal anisotropies", {A n }, which are orthogonal to the monopole, i.e. satisfy A n · 1 = 0. One would then find the unit-norm maps extremizing A n · F · A n under this additional constraint. This optimization problem admits N pair − 1 solutions, which, in addition to the monopole, form a basis of N pair maps which can be used to search for a monopole and anisotropies. However, the principal anisotropies derived in this fashion are in general statistically correlated with the monopole, A n · F · 1 = 0, thus would inflate the error bar on the monopole when included in a search.
To alleviate issue (ii), one could construct a set of "monopole-uncorrelated maps", {B n }, which are uncorrelated with the monopole, i.e. satisfy B n · F · 1 = 0. One would then find the unit-norm maps extremizing B n · F · B n under this additional constraint. This optimization problem also admits N pair − 1 solutions, which, in addition to the monopole, form a basis of N pair maps which can be used to search for a monopole and uncorrelated maps. However, the maps B n derived in this fashion in general have a non-zero projection on the monopole, B n · 1 = 0. Thus, detecting a non-zero amplitude for these maps would not imply that one has detected anisotropies in the data.
One could try and alleviate both issues simultaneously by searching for a set of normalized maps C n which are anisotropic and uncorrelated with the monopole, with extremal SNR 2 . As long as F · 1 is not colinear with 1 (which is the case if the monopole is not an eigenmap of the Fisher matrix), the two constraints are independent. The resulting optimization problem thus admits N pair − 2 solutions. To span the N pair -dimensional set of observable maps, one must supplement these monopoleuncorrelated principal anisotropies with the monopole and F · 1 − (F · 1)1 (properly normalized). Indeed, this additional map is orthogonal to the monopole and to all the C n maps (which stems from the condition that the C n are orthogonal to and uncorrelated with the monopole), thus linearly independent from all of them. However, this additional map is in general not statistically independent from the monopole nor from the maps C n . Therefore, this construction still does not alleviate the issue.
To conclude, unless the monopole is an eigenmap of the Fisher matrix, there is no good strategy to agnostically search for anisotropies within the standard Bayesian setup. In Section VI, we explore a frequentist approach to this problem.

VI. FREQUENTIST APPROACH TO AGNOSTIC ANISOTROPY SEARCHES
In the previous sections we have highlighted the difficulties in carrying a standard Bayesian search for anisotropies of pre-determined shapes alongside a monopole. Here we take a different approach, and derive a criterion to assess the presence of anisotropies in the data, regardless of their specific shape, without requiring a basis of maps on which to decompose the GWB.

A. Derivation of the frequentist criterion
Given the data translated into an estimator A for the GWB amplitude, one may seek the monopole amplitude that minimizes the χ 2 , given by The best-fit monopole amplitude is then simply The SNR of the best-fit monopole is then SNR bf = |A bf 0 | √ 1 · F · 1. This quantifies how well the best-fit monopole is detected relative to statistical noise. However, it does not quantify whether a pure monopole is a good fit to the data or not.
To quantify the goodness of fit, one can examine the χ 2 at this best-fit value. After simplification, one obtains where we have defined the projected Fisher matrix which satisfies 1 · F = F · 1 = 0. This implies that any monopole contribution to A drops out of χ 2 bf .
If the underlying GWB is purely isotropic, so that A is a monopole plus pure statistical noise, only the latter contributes to χ 2 bf , which has a chi-square distribution with N pair − 1 degrees of freedom. In particular, it has mean N pair − 1 and variance 2(N pair − 1).
In the limit that N pair 1, χ 2 bf is approximately Gaussian-distributed. A simple criterion for the detection of a anisotropy with 95% confidence (i.e. for the presence of a non-pure-statistical noise in A) is then Note that this criterion does not provide any information about the nature of the anisotropy (besides the fact that it must have some non-zero projections on one of the N pair − 1 eigenmaps of F ). It merely asserts that the data cannot be well described purely by a monopole. Consider a map A = A 1 M , where M = 1 is a dimensionless map. The minimum amplitude A 1 required for an anisotropy to be confidently detected in the data (without knowledge of its precise form) is thus This corresponds to a best-fit monopole A bf,min , hence to a minimum SNR for the best-fit monopole . (80)

B. Illustration with the EPTA
We illustrate this frequentist approach to the map M shown in Fig 10, which is normalized to have a unit monopole. First, we compute the best-fit monopole amplitude for this map, . We find A bf 0 = 0.96A 1 : the best-fit monopole amplitude turns out to be rather close to the actual value. This is because the projections of M and 1 on the first few principal maps happen to be similar in this case. This need not be always the case for other maps or PTAs. From Eq. (79), and using the EPTA Fisher matrix, we find A min 1 ≈ 1.2 × 10 −14 (this quantity has dimensions of characteristic strain). This corresponds to an SNR for the best-fit monopole amplitude 0.96 × A min 1 √ 1 · F · 1 ≈ 46. In other words, if the GWB has the angular dependence shown in Fig. 10, the full EPTA would only be able to establish the presence of anisotropies with 95% confidence (without being able to specify their precise form) after it detects the best-fit monopole with SNR 46.

VII. REPRODUCING TMG15'S RESULTS
To our knowledge, TMG15 is the only study of anisotropies in the GWB with real PTA data. In this section, we show how to use our Fisher formalism to reproduce their results. We moreover clarify the meaning of some of their upper limits. Indeed, as discussed earlier, a PTA is sensitive to N pair independent maps at most. This means that it is not possible to simultaneously constrain the GWB amplitude in more than N pair pixels in the sky, or constrain more than N pair of its spherical harmonic amplitudes. Yet, TMG15 present a 95% upperlimit map on the GWB in 12288 pixels, using only the 15 pairs afforded by 6 EPTA pulsars. They also present upper limits on spherical harmonic coefficients with max as large as 10, corresponding to 121 independent coefficients. Clearly, these cannot be true upper limits. In this section we shall understand precisely what these limits are, and reproduce them approximately using the Fisher matrix and its eigenmap decomposition.
A. Pixel-by-pixel upper limit maps Following TMG15, we define P ≡ A/A 2 h , which satisfies d 2Ω P (Ω) = 4π. With our notation, the N pair overlap reduction functions are given by Γ I ≡ 3 8 γ I · P . We now discretize the sky into N pix pixels centered atΩ i , with equal area ∆Ω. As in TMG15, we define c i ≡ P (Ω i ). The overlap reduction functions are then where the N pair × N pix matrix H has elements We can thus rewrite the relationship between the overlapreduction functions and the power distribution P (Ω) in matrix form, as Γ = H c, where Γ is the N pairdimensional vector with components Γ I and c is the the N pix -dimensional vector with components c i . Unless N pix = N pair , the matrix H is not square, and a fortiori not invertible. One can however define its Moore-Penrose pseudo-inverse H + . We may then define the N pix -dimensional vector Unless N pix = N pair , H + = H −1 and c = c. Instead, as we demonstrate in Appendix A, c i = P || (Ω i ) is the (discretized) power distribution projected on the N pair -dimensional space spanned by the γ I 's. In other words, c corresponds to the observable component of the GWB power, i.e. the only piece of P (Ω) that is measurable by a PTA. Importantly, it cannot be identified with c i = P (Ω i ) as done in TMG15. As a consequence, the map shown in TMG15 is only an upper-limit map for the projected (i.e. observable) GWB intensity. We refer to it as the Moore-Penrose pseudo-upper-limit map. We now explicitly construct an approximation of this pseudo-upper-limit map. Given estimators for the timing residual correlations, one may construct an estimator for the projected GWB map A || = A given by Eqs. (30)- (31). The pixel values of this estimator map are precisely In order to compute the variance of these values, it is best to first decompose it on the basis of principal maps M n : The N pair amplitudes A n are uncorrelated Gaussian variables with covariance We thus find From this we can estimate a 95%-sensitivity to A || (Ω i ), where the last equality defines the characteristic strain A 95% h (Ω i ).
We show A 95% h (Ω) in the top panel of Fig. 12. The qualitative and quantitative similarities with Fig. 2 of TMG15 are striking. The overall amplitude of our map is lower than that of TMG15, consistent with our lower 95% sensitivity estimate for the monopole amplitude A 95% h . Note, also, that we cannot hope to reproduce precisely the same map as TMG15, as their map is constructed from real data, i.e. a non-zero realization of the noise.
Let us point out that the pseudo-upper-limit map derived with this procedure gets worse if one includes more pulsars. Indeed, as can be seen from Eq. (86), the variance of A || (Ω i ) is dominated by the noisiest eigenmaps. As more pulsars are added, and the number of eigenmaps grows, the noisiest ones become increasingly noisy. Conversely, using only one pulsar pair results in a pseudo-upper-limit map with an overall amplitude lower by nearly one order of magnitude, as can be seen in the bottom panel of Fig. 12. In fact, in the right hemisphere, the map amplitude nearly vanishes, as does the pairwise-timing response function of the single pulsar pair. Clearly, adding more information should only tighten true upper limits, not make them weaker. This demonstrates once again that these Moore-Penrose pseudo-upper-limit maps cannot be interpreted as true pixel-by-pixel upper limits. 9 19 FIG. 12. Moore-Penrose pseudo-sensitivity maps produced with 6 EPTA pulsars (top) and one single EPTA pulsar pair (bottom). The qualitative and quantitative similarities of the top map with Fig. 2 of TMG15 are striking. The overall lower amplitude of our map relative to that of TMG15 is consistent with our lower 95% sensitivity to the monopole amplitude A h . These maps cannot be interpreted as actual upper limits on the GWB intensity in each pixel, as a PTA cannot constrain more than Npair independent components simultaneoulsy. This is further demonstrated by the fact that the pseudo-sensitivity map produced with only two pulsars (hence less data) has lower amplitude than the one produced with 6 pulsars.

B. Constraints to spherical harmonic coefficients
We now apply the formalism of Section IV D to the subset of 6 EPTA pulsars used in TMG15, and compare our results to theirs. With the 15 pairs afforded by 6 pulsars, cross-correlation data can simultaneously constrain spherical-harmonic coefficients only if max ≤ 2.
For these values of max , our 95% sensitivity estimates on C 1 and C 2 , shown in Fig. 13, are remarkably close to their 95% upper limits in the absence of physical prior (see their Fig. 1, right  we find that the 95% sensitivity to the monopole amplitude C 0 gets significantly degraded as max increases from 0 to 2, while the upper limit of TMG15 is virtually unaffected. We attribute this to their use of autocorrelations, which we do not include in our analysis. We note, however, that since the GWB dipole and quadrupole affect autocorrelations, increasing max from 0 to 1 and then 2 should result in poorer constraints on the monopole. It is unclear why this is not what is found in TMG15. Let us now understand TMG15's upper limits on the spherical-harmonic amplitudes for max > 2, for which there are more independent coefficients than independent data. In this case, TMG15 also define c m given the overlap reduction function Γ I through the Moore-Penrose pseudo inverse of the linear relation between Γ I and c m . Again, unless the number of harmonic coefficients is equal to N pair , c m = c m . For finite max , there is no simple intuitive interpretation of thec m , but for max → ∞, it is straightforward to show that the c m are the sphericalharmonic coefficients of the projected GWB. Specifically, using our normalization conventions, we define A m ≡ Using Eq. (84), we thus have Using Eq. (85), we thus obtain (89) We then define C ≡ m A 2 m /(2 + 1). The mean and variance of these coefficients is given by Eqs. (61) and (62), with the substitution cov(A m , A m ) → cov( A m , A m ), from which we can infer a 95% sensitivity C 95% as in Eq. (63). We show C 95% in Fig. 14. FIG. 14. 95% sensitivity to the Moore-Penrose pseudo C 's in the limit max → ∞, for 6 EPTA pulsars (upper curve) and 2 EPTA pulsars (lower curve). Compare the upper curve with the right panel of Fig. 1 of TMG15, for max = 7, 10. The C are obtained from the harmonic amplitudes of the projection of the GWB on the Npair-dimensional space of observable maps, and therefore are not equal to the true harmonic amplitudes of the GWB, see text for details.
This can be compared to the right panel of Fig. 1 of TMG15. The similarity of their max = 7, 10 results with ours (which holds in the limit max → ∞) is striking. Note again that we cannot expect to reproduce exactly TMG15's results, which are based on actual data, i.e. a particular realization of the noise. Here again, a sharp illustration of the fact that these Moore-Penrose pseudoupper-limits cannot be interpreted as upper limits on the C is the fact that C 95% decreases if one only includes 2 pulsars instead of 6, as can be seen in Fig. 14.

VIII. CONCLUSIONS
In this paper, we have built on the Fisher formalism we introduced in Paper I, and applied it to EPTA data, since this PTA is currently the only one with published limits on GWB anisotropy (TMG15). We showed through multiple examples how this framework can serve to understand the sensitivity of a PTA to a general, anisotropic GWB, and can be used to guide and optimize data analysis. We also showed how to recover existing EPTA results, and clarified their meaning, in the light of our improved understanding of the sensitivity of PTAs to GWB anisotropies. In particular, we pointed out that one cannot simultaneously detect -or constrain beyond physical priors -more than N pair components of the GWB, be they pixel values or spherical harmonic amplitudes.
The fundamental tool of our formalism is the Fisher "matrix" for the direction-dependent GWB intensity. This matrix condenses the essential noise characteristics of a PTA in a compact form. It can be computed given the positions of the pulsars and their characteristic noise strains. The latter requires analyzing the timing data of individual pulsars, but not conducting any crosscorrelation analysis with real data.
We started by using our formalism to find the pulsar pairs which are the most sensitive to an isotropic GWB with a power-law frequency dependence h c (f ) ∝ f −2/3 . We identified the 44 best EPTA pairs (out of the 861 available) that provide 90% of the SNR 2 . With the Fisher formalism, finding these best pairs for any PTA is a trivial exercise, and can be used to dramatically speed up real data analyses, at very little cost in sensitivity.
In the remainder of the paper we focused our attention to anisotropies in the GWB. We first considered the case where the GWB is assumed to be a linear combination of a finite number of known maps (no more than N pair ), such as spherical harmonics, coarse pixels, or hot spots in pre-determined directions. We showed that the sensitivity of a PTA to the component of the GWB on any of the basis maps systematically and quickly degrades as the number of basis functions is increased. This stems from the fact that the noise properties of standard basis maps are correlated; therefore, adding more functions to the search changes the overall covariance structure of the entire set of amplitudes. In particular, the ability to detect the monopole through pulsar cross correlations is systematically degraded as one searches simultaneously for anisotropies. We did not quantify the ability of a PTA to set upper limits through pulsar autocorrelations (or "common red noise" analyses), but pointed out that these are sensitive to the GWB monopole, dipole and quadrupole. Hence, upper limits on the monopole should also degrade as one accounts for anisotropies.
Having shown that standard basis maps are poorly adapted for agnostic GWB searches, we then constructed the "principal maps" of a PTA, which are the eigenmaps of the Fisher matrix. These maps form an orthonormal and uncorrelated basis of the N pair -dimensional space of observable maps. In other words, they represent the N pair statistically independent pieces of information accessible to a PTA, in the space of GWB maps. We showed how one can in principle "reconstruct" the observable part of the GWB as a linear combination of principal maps detected with sufficiently high significance, a procedure similar to the making of "dirty maps" in radio interferometry. However, in contrast with a radio interferometer, for which different "visibilities" have similar noise levels, principal maps have highly unequal noise properties. As a consequence, the ability to reconstruct the GWB is strongly dependent on its overall amplitude. For the EPTA, we found that even detecting a small number of principal maps would typically require a GWB amplitude largely exceeding current bounds.
Principal maps are a useful basis if one is fully agnostic about the angular dependence of the GWB: they allow for an optimal search "under the lamppost" of a given PTA. However, the GWB is expected to be predominantly isotropic; physical considerations therefore mandate including the monopole as part of the search, and being able to identify anisotropic components of the GWB. We showed that the principal maps are in general correlated with the monopole and have a non-vanishing isotropic component, and are thus not adapted for such searches. We argued that one cannot construct a complete set of anisotropies which are uncorrelated among themselves and with the monopole, because the latter is not an eigenmap of the Fisher matrix.
In summary, there are two options when searching for a general, direction-dependent GWB within a standard Bayesian framework. On one hand, one can include the monopole in GWB searches, alongside anisotropic basis functions, but then suffer from the loss of sensitivity resulting from correlations with the monopole. On the other hand, one may search under the lamppost of PTAs with principal maps, being fully agnostic about the angular dependence of the GWB, but then give up the special status of an isotropic GWB. To circumvent these limitations, we proposed a frequentist criterion to determine whether the GWB is consistent with isotropy, even if one cannot determine the specific type of anisotropy. This criterion simply relies on computing the chi-squared statistics of the best-fit monopole amplitude, and comparing it with its statistical noise. We illustrated this with a fictitious, highly anisotropic GWB map, and showed that the EPTA would have to first detect its best-fit monopole with a SNR of 46 before being able to confidently assert that the data is inconsistent with a purely isotropic GWB.
The prospects for current PTAs to detect fully unknown anisotropies in the GWB appear to be somewhat limited, but the situation may improve dramatically with the several hundreds (or even thousands) of millisecond pulsars that SKA is expected to detect [15]. Indeed, with such a large number of pulsars, the GWB monopole could be detected with a large SNR (while not exceeding the weak-signal limit), and one may then be able to search for anisotropies, even if their amplitude is significantly weaker than that of the monopole. In addition, there may be hope to detect anisotropy in the GWB through its statistical properties, for instance through its correlations with tracers of large-scale structure [69]. We will explore these avenues in future work.
In closing, let us note that our Fisher formalism has far broader applications than what we considered in this work. First, we assumed an overall known frequency dependence, identical for the monopole and anisotropies. It would be straightforward to generalize our work to arbitrary frequency dependence, possibly dependent on the angular structure -something which would be difficult to do with a full-on data analysis. One could also include other sources of correlated noise, such as global clock or ephemeris errors [36,37], and go beyond current studies which ascribe them the same frequency dependence as the GWB. In addition, we could generalize our Fisher matrix beyond the weak-signal limit, and consider, for instance, the weak-anisotropy limit. Lastly, one could easily extend the formalism to include polarized GWBs.
Overall, this formalism is a valuable tool that can be used to guide and cross-validate any full-blown data analysis. It will be especially important when SKA harvests hundreds of new millisecond pulsars, out of which PTAs with tens of thousands of independent pulsar pairs can be constructed. Once the GWB monopole is confidently detected by a PTA, the Fisher formalism will be extremely valuable in determining the next steps forward.

ACKNOWLEDGMENTS
We thank Jeffrey Hazboun, Stephen Taylor, and Joseph Romano for useful discussions.
This work used the Strelka Computing Cluster, which is run by Swarthmore College. YAH acknowledges support from the National Science Foundation through grant number 1820861. TLS acknowledges support from NASA (though grant number 80NSSC18K0728) and the Research Corporation. The Flatiron Institute is supported by the Simons Foundation.
Appendix A: Proof that ci = P || (Ωi) In this appendix, completing Section VII, we prove that the vector c ≡ H + Γ = (H + H) c is the GWB angular distribution projected on the N pair -dimensional space observable by a PTA. Thus, in general, c = c.
To see this, recall that the Moore-Penrose pseudoinverse matrix is defined as H + = H T (HH T ) −1 . As a consequence, H + H is clearly a projector, since This achieves the proof that H + H is the projector on the space spanned by the N pair maps γ I . Therefore, as announced, c i = P || (Ω i ) is the GWB angular power projected onto that N pair -dimensional space.