{'fraction_non_alphanumeric': 0.03182823682498374, 'fraction_numerical': 0.005256994144437216, 'mean_word_length': 4.654208357857564, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Comparison of maximum likelihood mapping methods for gravitational-wave backgrounds

Detection of a stochastic background of gravitational waves is likely to occur in the next few years. Beyond searches for the isotropic component of SGWBs, there have been various mapping methods proposed to target anisotropic backgrounds. Some of these methods have been applied to data taken by the Laser Interferometer Gravitational-wave Observatories (LIGO) and Virgo. Specifically, these directional searches have focused on mapping the intensity of the signal on the sky via maximum likelihood solutions. We compare this intensity mapping approach to a previously proposed, but never employed, amplitude-phase mapping method to understand whether this latter approach may be employed in future searches. We build up our understanding of the differences between these two approaches by analysing simple toy models of time-stream data, and run mock-data mapping tests for the two methods. We find that the amplitude-phase method is only applicable to the case of a background which is phase-coherent on large scales or, at the very least, has an intrinsic coherence scale that is larger than that of the detector. Otherwise, the amplitude-phase mapping method leads to a loss of overall information, with respect to both phase and amplitude. Since we do not expect these phase-coherent properties to hold for any of the gravitational-wave background signals we hope to detect in the near future, we conclude that intensity mapping is the preferred method for such backgrounds.

Detection of a stochastic background of gravitational waves is likely to occur in the next few years. Beyond searches for the isotropic component of a stochastic gravitational-wave background, there have been various mapping methods proposed to target anisotropic backgrounds. Some of these methods have been applied to data taken by the Laser Interferometer Gravitational-wave Observatory (LIGO) and Virgo. Specifically, these directional searches have focused on mapping the intensity of the signal on the sky via maximum-likelihood solutions. We compare this intensity mapping approach to a previously proposed, but never employed, amplitude-phase mapping method to understand whether this latter approach may be employed in future searches. We build up our understanding of the differences between these two approaches by analyzing simple toy models of time-stream data, and run mock-data mapping tests for the two methods. We find that the amplitude-phase method is only applicable to the case of a background which is phase coherent on large scales or, at the very least, has an intrinsic coherence scale that is larger than the resolution of the detector. Otherwise, the amplitude-phase mapping method leads to an overall loss of information, with respect to both phase and amplitude. Since we do not expect these phase-coherent properties to hold for any of the gravitational-wave background signals we hope to detect in the near future, we conclude that intensity mapping is the preferred method for such backgrounds.

I. INTRODUCTION
Over the past two decades, several efforts have been made to understand the problem of gravitational wave (GW) map-making. The term "map-making" [1,2] is typically used when considering stochastic GW backgrounds (SGWBs), which are expected to be incoherent superpositions of GW signals arriving at a detector from all directions of the sky. These GWs are generated by a variety of underlying mechanisms, which include both compact or diffuse sources at both astrophysical or cosmological distances (see e.g., [3,4]).
The primary map-making method adopted in searches for anisotropy in stochastic signals is a maximum-likelihood approach that targets the overall background intensity [5], entirely discarding any time-domain phase information present in the data. This method, hereinafter referred to as intensity mapping, has been used on cross-correlated data from the Laser Interferometer Gravitational-wave Observatory (LIGO) detectors [6][7][8] and, more recently, also on data from the Virgo observatory [9]. We review the basics of this method in Sec. III. Intensity mapping is best suited for a background that is truly stochastic. In this case the signal is characterized by waves that have random time-domain phases but whose intensity varies as a function of frequency and, in principle, angular direction. This motivates the development of methods that compress the data by discarding the time-domain phase information of the signal. For noise-dominated detectors, such as LIGO-Virgo, this is achieved by integrating the equal-time cross-correlation of the time stream of pairs of detectors. Each detector pair has a characteristic correlated sky response to the GW signal, which is time-dependent following the Earth's rotation. The set of responses determine the mapping capabilities of the detector array.
Other methods have been proposed that solve for the phase information as well, effectively mapping the amplitude and phase-or real and imaginary parts-of the two GW polarization modes on the sky [10][11][12]. We refer to these as amplitude-phase mapping in what follows. A clear motivation for preserving phase information would be the case of a diffuse, anisotropic, but coherent GW background. It is often argued that the primordial background generated by an inflationary epoch satisfies this condition [13,14]. This is due to the squeezing of modes induced by any period of superhorizon evolution. However, it has been shown that, at accessible frequencies, any primordial coherence is wiped out by the effect of large-scale structure on the GWs as they propagate through the evolving universe [15,16]. Hence, reconstructing the phase information itself is effectively not necessary in stochastic analyses; however there remains a question of whether it is actually possible to determine the phase information for a stochastic source. This question is what we set out to answer with this paper.
In this paper we compare the application and validity of intensity and amplitude-phase mapping methods for SGWB analyses. We do this by analyzing mock data, focusing on the cross-correlation of measurements made by a network of detectors. It is important to note, at this stage, that a key element in any stochastic analysis is the accurate modeling and estimation of detector noise. We do not focus on this aspect here since our aim is to clarify the applicability of the mapping methods with respect to the signal properties. Our analysis assumes the simplest possible noise component and artificially high signal-to-noise ratios. Our conclusions, nonetheless, are independent of these choices. This paper is organized as follows: in Sec. II, we review the strain description of the GW signal and the key assumptions made throughout the paper. Sec. III details the two map-making methods. In Sec. IV we discuss different measurement scenarios to build an understanding of the issues encountered when reconstructing an anisotropic SGWB using intensity and amplitude-phase mapping methods. In Sec. V we apply the two methods to mock datasets and provide a comparison between the mapping fidelity. We conclude with Sec. VI, where we discuss our findings and future prospects.

II. GW SIGNAL
The stochastic GW metric perturbations at time t and position vector x may be written as an infinite superposition of plane waves having polarization A, frequency f , and arriving from directionn = (sin θ cos φ, sin θ sin φ, cos θ) [5]: where (θ, φ) are the standard angular coordinates on the 2-sphere and the spatial wave vector is written explicitly as k = −2π|f |n. We set the speed of light c = 1 here for simplicity. The metric perturbations h ij (t, x) are real; hence a reality condition is imposed on the complex, frequency-domain modes, h A (f,n) = h A (−f,n). In the above expression, we choose the linear polarization basis A = {+, ×}, where the orthogonal polarization basis tensors e A may be written as are the standard unit vectors tangent to the sphere. Since {h + , h × } are complex valued, we can write them in terms of either their real and imaginary components or their amplitude and temporal phase, all of which, in general, will depend on the frequency and direction of the waves. For a stochastic background, the metric perturbations h ij (t, x) and hence the Fourier components h A (f,n) are random fields, whose probability distributions define the statistical properties of the background. For the following discussion, we will assume that the background is (i) Gaussian, (ii) stationary, and (iii) unpolarized, which means that (i) the statistical properties of the random fields are completely characterized by their first-and second-order moments, (ii) there is no preferred origin of time (implying that random variables corresponding to different frequencies are statistically independent of one another), and (iii) the statistical properties of the background are invariant under rotations of the polarization tensors in the plane perpendicular ton (implying statistically independent and equivalent + and × polarization components). From these assumptions it follows that that {h R A , h I A } are statistically independent variants drawn from the same Gaussian distribution with zero mean and half the variance of h A , and that the phase ϕ A is uniformly distributed between 0 and 2π and is statistically independent of the amplitude A A , which is Rayleigh distributed being the square root of the sum of squares of two statistically independent Gaussian random variables each with zero mean and equal variance.
We will further assume (iv) that the background has no nontrivial phase coherence across the sky. This means that the phase of the GWB signal coming from two different directions on the sky are statistically independent of one another 1 : Although primordial backgrounds, such as those generated during an epoch of inflation, may be phase coherent as modes reenter the horizon and begin to propagate, this coherence is lost through propagation effects [16]. It is therefore reasonable to assume that phase incoherence is a generic feature of any stochastic background that may, or may not, have angular correlations in amplitude. In other words, the intrinsic angular scale ∆n over which the phases are correlated goes to 0, even though the correlation scale for the amplitude may be finite. Hence, any attempt to measure the GW phase using a detector with finite angular resolution will necessarily average the true phase over this angular resolution scale. This loss of information will degrade reconstruction of the amplitude if one tries to estimate it from the real and imaginary parts of the Fourier components. We will show this explicitly in Secs. IV and V. Putting together all of the above results, we can write where defines the intensity of the GWB as a function of frequency f and directionn. Here T is the total observation time.
The angle brackets . . . denote ensemble averaging over the random amplitudes and phases of the Fourier coefficients of the metric perturbations at a fixed spatial location x, under the assumption of time stationarity. In practice, this averaging is realized by averaging over all the available GW time-series data, assuming that the background is ergodic. The distribution of energy and matter in the universe, e.g., large-scale structure, which gives rise to the GW background is fixed with respect to this averaging process. The GW background may be isotropic, anisotropic, or statistically isotropic (i.e., invariant under arbitrary rotations of the sky) depending on the statistical properties of I(f,n) with respect to the sky directionn. For example, for both isotropic and anisotropic backgrounds, the intensity field is a deterministic quantity; it is independent of sky direction for a purely isotropic background-i.e., I(f,n) ≡ I(f ), and has preferred directions for an anisotropic background. For a statistically isotropic background, the intensity is a random field, assumed here to be approximately Gaussian, whose mean is independent of sky direction and whose quadratic expectation values depend only on the angular separation between two points on the sky as where P is the Legendre polynomial of order . In the above expressions, . . . Ω denotes averaging over different GW universes [17] (e.g., over different realizations of large-scale structure) which are drawn from a rotationally invariant probability distribution described by the angular power spectrum C (f ). We note that there is a subtle difference if one imposes statistical isotropy on the h A fields; for more details regarding this see [11]. In all cases, the intensity I(f,n) may be related to the fractional energy density parameter Ω GW (f,n) via [5] Ω which is the fundamental relation that allows one to connect GWB observations to the cosmological implications of the background. Integrating the above equation over direction on the sky yields which relates the monopole components of the fractional energy density parameter and the intensity of the background. Finally, we note that it is common in the literature to assume that the intensity factorizes as where E(f ) and I(n) encode the spectral and directional dependence of the background, respectively. We employ this assumption throughout as it considerably simplifies map-making.

III. MAP-MAKING METHODS WITH GW DETECTOR DATA
In this section we present two proposed maximum-likelihood map-making methods for GWBs in a common formalism. First we review the intensity mapping method which has been previously presented in several references, e.g., [5,18,19], and applied consistently to LIGO data [7][8][9][20][21][22]; then, we lay out the amplitude-phase mapping method which follows the presentation in [12].

A. Intensity mapping
Intensity mapping, such as that used for map-making with the LIGO-Virgo detectors, works with the crosscorrelation of time-coincident data directly to discard the autocorrelated noise terms which would otherwise dominate the calculation.
To start, we consider a set of detectors i = {1, . . . , N }. The data collected by detector i in a time segment τ may be considered as made up of separate signal and noise components, d τ i (t) = s τ i (t) + n τ i (t). We write the corresponding discrete Fourier transform asd where we have dropped the tilde and detector label in favor of the more concise boldface vector notation. As we are dealing with a real time-stream Fourier transformed into a discretized Fourier space, it is convenient to relabel frequencies as a discrete index, which picks out a single frequency mode in the transform. The signal component in the data is modeled as where R is the response function spanning detector space. Note that the h field has no τ dependence as it is assumed to be stationary. A pair of detectors i, j then observe d τ,f i , d τ,f j respectively; these form a baseline and we can consider a correlated data vector as spanning the space of different baselines directly, We can write down the likelihood for the residuals of the cross-correlated data as Note that we are using shorthand notation here, omitting τ and f everywhere in Eq. (18). As in [19], we take the signal model S for the cross correlation to be which may be derived directly from Eq. (9). Here Γ is the cross-correlated response vector to the GW intensity, obtained via the squared sum of the response terms above as The spectral dependence of I has been assumed factorizable here and absorbed into Γ for simplicity, assuming Eq. (15). The time dependence (encapsulated in τ ) and spatial dependence of the response function are fundamental to mapmaking, as they define the scan strategy of the set of baselines, which sets the resolution of the measurement. More details on this may be found in [1,[23][24][25][26].
The noise covariance matrix may be written in baseline space, assuming uncorrelated noise between detectors, where P i is the (two-sided) noise power spectrum in detector i. Maximizing the likelihood (18) above yields the mapping equation for GW intensity We refer to z as the projection of the dataset into pixel space, while M is the Fisher matrix of the mapping problem. These are constructed from the following quantities calculated at the individual f and τ , as each of these constitutes an independent measurement according to the likelihood in Eq. (18), which are then summed over all times, frequencies, and baselines ij in the set to obtain the full (and most informative) projection and Fisher matrix, Note that the integration over frequencies here requires an assumption for the frequency dependence of the GWB, which we have transferred to the response function. One could also choose to not integrate over frequencies, and solve for maps mode by mode, avoiding this frequency modeling step. At this stage, it is necessary to pick a working resolution on the sky in order to explicitly carry out the calculations. Using the healpix Python package to deal with pixelization, we can set a working resolution and translate sky directionn to pixel p. The data model then becomes The components of Eq. (23) can then be considered in the pixel domain, wheren → p.

B. Amplitude-phase mapping
Let us now take a step back and reconsider the data model in Eq. (16), which we can write compactly as d τ,f = s τ,f + n τ,f . Assuming zero-mean Gaussian noise in the strain, i.e., n τ,f = 0, Eq. (17) serves as our signal model. As {h + , h × } are complex valued, there are four independent (real) fields on the sky to estimate: We write down the likelihood in terms of the residuals in the strain as assuming each time segment τ and frequency f in the dataset are statistically independent of one another, where we are again employing shorthand notation, using boldface font to represent a vector that now spans just the space of the detectors (and not baselines). C n is the noise covariance C n = n ⊗ n . As the noise is modeled independently and is not part of the maximum likelihood estimation, it is possible to reduce the mapping problem to a closed-form χ 2 solution minimizing We solve for the four h fields separately; hence it is useful to decompose the signal model into real and imaginary components simply as where the spectral dependence of the stochastic field has been explicitly factored out of h A (n) and has been absorbed into the response term R so as to keep track of fewer dependencies. Minimizing χ 2 with respect to each field yields the maximum likelihood solution where Each component is, explicitly, For the sake of conciseness we do not include the frequency and time labels in each of the terms above; these mirror those in Eq. (23). Note that the projection z and Fisher matrix M here are not the same as for intensity mappinghowever, we have chosen to keep the same notation to draw the comparison between the two approaches. In the uncorrelated noise case the (noise) covariance matrix becomes diagonal in detector space, so calculations in Eq. (31) simplify considerably. Furthermore, as in the intensity mapping approach described above, the integration over frequencies in Eq. (30) requires an assumption for the spectral shape of the signal. However, in the amplitude-phase case one must take extra care: when broadband integrating here one must assume a certain phase-coherence across modes, or else the phase introduced in Eq. (6) will (correctly) average to zero and yield null sky maps. Finally, as in the previous section, we can discretize the sky and transform sky directionn to pixel p, to apply this method to data. The signal model becomes and similarly the terms in Eq. (31) may be translated into the pixel domain.

IV. INSTRUCTIVE CONSIDERATIONS
Before discussing the application of the map-making methods laid out above, let us start with some useful considerations about the nature of the measurement of broadband, stochastic GWs from all sky directions with a set of interferometers. In particular, let us focus on the notion that a detector naturally low-pass filters the signal, both in terms of temporal resolution and angular resolution; the subsequent digitization/sampling of the data needs to respect the maximum temporal and spatial frequencies present in the signal if aliasing of power is to be avoided.
It is useful to first break down the measurement into two components: the time or frequency-domain measurement, and the sky-dependent response. These two steps are inextricable in a real GW detector. However, to get an idea of the accessibility to each information component we can consider two simple scenarios separately: first, we can focus on the estimation of the amplitude and phase of a complex field in the frequency domain through time-domain measurements; then we can analyze the effect of a window-averaged response on the measurements of a complex field on the semicircle, in analogy with what happens with GW detectors observing the full sky.
When analyzing GW detector timestreams that contain a measurement of a stationary stochastic field, say h(t), we typically start by taking the Fourier transform (FT) of the data, as the frequency domain representation allows us to access information more efficiently. A way to see this is to note that the quadratic expectation value in the time and frequency domains are related by Thus, the Fourier transform of a stationary timeseries maps our measurement to the diagonal space for the autocorrelation of the field, where P (f ) is the (two-sided) power spectrum of the field. In other words, while in the time domain the correlation between the stationary field at different times depends (only) on the time difference t − t , in the conjugate space (frequency space) the correlation between different frequency components is a delta function, and thus the frequency dependence is totally "compactified" in the power spectrum. We can either write an estimator for a real field, for example, P where the phase ϕ has a certain functional dependence on direction. This is effectively a frequency-independent, unpolarized version of the GW strain in the Fourier domain given by Eq. (6). Note that, as in the relation t FT ←→ f , there is a useful conjugate space to directionn space, where a m are the spherical harmonic coefficients ofh(n), and C is the angular power spectrum ofh. The C s are related to C(n ·n ) via where P (x) are the Legendre polynomials. The spherical harmonic decomposition (SHD) introduced in Eq. (36) allows us to transform the field from directional spacen to degree and order m space, and, again, the correlation is diagonal in the latter. The measurement ofh, however, does not occur inn space directly, but in the time domain; the relation between these two spaces depends on the specific detector used to make the measurement. This is the crucial difference between measurements of a frequency-dependent GW observable and a direction-dependent one.
Let us now discuss an instructive example. Let us consider a case where we can fix the observation to a single value of φ on the sphere, such that we reduce the mapping problem to the estimation of a one-dimensional field on the semicircle.h is then parametrized solely by 0 ≤ θ ≤ π on the semicircle, such that effectively we are observing h(θ, φ = const.) ≡h(θ) = A(θ)e iϕ(θ) . We assume that our detector makes a time-dependent measurement r(t) of the field on the semicircle, filtered by the detector response R(t, θ), which scans the semicircle as time goes by. Note here the similarity with the model for the GW strain in Eq. (17), relevant to the measurement of anisotropic GWBs. The intrinsic coherence scale of the field, ∆θ, sets the input resolution of our example, inducing π/∆θ independent samples on the semicircle. We take the phase ϕ(θ) to be random at the intrinsic scale, such that the field has statistically independent phases at each sample. The amplitude A is assumed to vary much more slowly than the phase, such that the scale over which the amplitude varies appreciably is much larger than ∆θ. We model the instantaneous response of the detector at time t 0 as a top-hat function with width ∆θ d , and the time dependence is imposed by the scanning strategy, i.e., the angular function θ(t). In this simple example we take θ(t) = π t/τ where τ is the time period over the semicircle, such that consecutive measurements correspond to consecutive angular samples on the semicircle, and continuity between measurements at 0 and π is ensured. The measurement ofh is then mediated by the detector-induced window width ∆θ d , and the comparison between ∆θ and ∆θ d will determine to what degree the field is resolvable. ∆θ d is the angular resolution of the measurement, and as long as ∆θ d is sufficiently less than ∆θ, the field is overresolved and its amplitude and phase information may be measured exactly. However if ∆θ d > ∆θ there is loss of information as the window averages over the phases, which are randomly distributed along the semicircle. In the extreme case where ∆θ d ∆θ, the measurement becomes compatible with 0, and the amplitude of the field is entirely lost. Translating the measurement into the conjugate space, one may interpret this limit as an angular Nyquist frequency, Nyq ≡ π/(2∆θ d ), imposed by the detector. In this sense, the measurement only works when Nyq > max , where max ∼ π/∆θ is the maximum angular frequency of the signal. Note the analogy then between ∆θ d , and the sampling rate ∆t d described above. The useful estimator for this example is one that marginalises over the phases and aims for the field intensity C(θ) directly, which is obtained from Eq. (36) assuming total phase incoherence ofh(θ). Equal-time measurements r 1 (t), r 2 (t) made by two identical detectors 1, 2 then satisfy where the angle brackets introduce an expectation value similar to Eq. (9); i.e., they imply ensemble averaging over different data samples, and the Dirac delta in Eq. (40) has already been applied. Mapping then amounts to inverting Eq. (41) to estimate C(θ). Hence to preserve the amplitude information it is necessary to cross-correlate (i.e., square) signals before averaging over the response window. Note that the spatial averaging introduced in Eq. (38) is analogous to what happens in a GW detector: the directional information is not directly accessible in the measurement, but rather needs to be reconstructed based on observation features, such as the time dependence of R.
In the simple case above the only relevant coherence scale of the field is that of the phase; however, in general it is also necessary to worry about the coherence scale of the intensity, ∆θ C . Note this is directly related to the coherence scale of the amplitude, as in this case C = A 2 . When ∆θ C < ∆θ d , this will have an impact on the observation of C (or, equivalently, of A, taking the square root of the intensity) similar to what occurs for ϕ. However, note that as the intensity is positive definite, this field will never average to zero but would rather approach zero from above in the limit of ∆θ C ∆θ d .

V. MOCK-DATA MAPPING TESTS
Simplified mock-data mapping tests are presented here to illustrate the points made in the previous sections. We consider the two LIGO interferometers that have a well-known response function, illustrated, for example, in [18]. In these examples, the noise is taken to be white, Gaussian, stationary, and much weaker than the signal. This is to highlight the potential and shortcomings of the mapping algorithms based on the characteristics of the signal component only. As the noise is assumed to be known, it is not part of the estimation. In practice, all GW mapping attempts to date using ground-based detectors have relied on independent estimators for the noise, which is a valid approach when the signal is entirely subdominant and hence does not bias the noise estimation, at least on a segment-by-segment basis. This methodology needs to be revisited in the presence of competing signal and noise components-an example of this may be seen in [27] for the Laser Interferometer Space Antenna. The reconstructions presented here are obtained following the data-handling recipe of [19] which imitates the steps performed on real LIGO data [28,29]. In both intensity and amplitude-phase mapping, the mock data are generated in segments that represent Fourier transformed one-minute segments of a detector array time-stream. The one-minute time scale is chosen such that the response of the detectors may be considered constant on the sky throughout the segment. This not only sets a natural lower bound on the frequency range probed, but also imposes a pixelization scheme: in order for the reconstruction to be complete, the sky response should vary smoothly from segment to segment.
The most challenging part of this procedure is the safe inversion of the Fisher matrix. This arises since the spatial sampling of the sky is suboptimal, as explored, for example, in [19]. However, as long as the noise is subdominant it is possible to pseudo-invert M discarding a minimal number of modes and recovering the sky signal perfectly, as shown below. The situation capsizes in the presence of noise: in that case, it is necessary to choose a cutoff on the singular values of the eigenmodes of M , to differentiate between the signal and noise modes. This is a highly nontrivial problem, so we do not discuss it here-see the discussion in [9] for example.
In all mapping examples presented here, the signal is modelled as a diffuse point source on the sky, with varying phase as a function of direction. This is obtained in practice by associating a random phase to each pixel of the amplitude sky map, which is then scanned to prepare the data segments. This is explicitly performed via Eq. (33) in the case of amplitude-phase mapping, and Eq. (19) in the case of intensity mapping. The pixelization schemes used for data generation and map reconstruction need not match-in fact, one may not in general expect that the natural coherence scale of the field be comparable to that of the detector. The size of the unit pixel at injection, 4π/N in pix , is analogous to ∆θ in the example made in Sec. IV, while the size of the unit pixel at reconstruction, 4π/N out pix , is analogous to ∆θ d . For simplicity, we do not add a frequency-dependent phase term, although in reality it is present and may not be neglected in a real GW measurement; we discuss this further below. The top row shows reconstructed maps with N out side = 16, which corresponds to the case when the intrinsic coherence scale of the field is equal to that of the detector. The second and third rows show reconstructions with N out side = 4 and N out side = 2, respectively. The bottom row shows the corresponding recovered intensity maps for the three different N side cases. In the case when N out side = N in side , the recovered intensity matches exactly the one shown in Fig. 1; in the other cases, the intensity is degraded.
phases of the h A fields in each pixel are perfectly recovered. However, the second and third rows of Fig. 2 show two cases where the output resolution is lowered from the input N in side = 16 to N out side = 4 and N out side = 2 respectively; here, the phase components are averaged out within neighboring pixels, hence the phase information is lost. This leads to an overall loss of intensity information as well, by the same reasoning described for the example in Sec. IV. Note that style choices have been made to underline the effect: the values chosen for the color bars in the intensity I plots in Fig. 2 match the correct values one would obtain when coarse-graining the I in intensity map to the relevant N out side . The coarse-graining exercise is performed with the ud grade function of the healpy package [30].
This effect does not occur in the case of intensity mapping, since the phase information is not required to reconstruct the intensity on the sky. We have repeated the exercise above using the same input map as in Fig. 1 employing the intensity mapping algorithm, varying the values of N out side , and find that the I out maps agree with the ud graded input maps within less than 1%. To quantitatively compare intensity and amplitude-phase mapping, we show in Fig. 3 the difference between the recovered maps with intensity mapping, I out,I−map , and amplitude-phase mapping, I out,AP−map , in the cases where N out side < N in side . In both cases, it is clear that most of the overall power is lost when performing amplitude-phase mapping, and in particular in the case where N out side = 2, the injected point source is hardly discernible from the fluctuations around 0 in the pixels.
Hence, we confirm that in the case of intensity mapping, the overall GW intensity is conserved, and the structure  Fig. 2. On the right, the difference maps Iin − Iout are shown to quantitatively comparative the performance of the algorithms. Note that the input map has been ud graded to allow for this comparison; the actual input map is at N side = 16 as shown in Fig. 1. on the sky as well, as long as the coherence scale of the intensity is larger than that of the detector. In Fig. 4, we show an example where this is not the case: here, the injected signal is a very high resolution input map, N in side = 32, which is null everywhere apart from a tiny patch, mimicking a pointlike source, and is recovered at N out side = 8. Here indeed the signal from the (almost unique) incoming direction is averaged within the neighboring pixels, and hence the point value in that specific direction is much lower than the injection (by a factor of ∼ 10). Note, however, that the recovery is perfectly in line with a coarse-graining of the input map: the monopole is conserved in the operation, and the correct value is recovered in the single "hot" pixel. Hence, in this case there is as little loss of information as possible. But this may present an issue when the signal is competing with a high level of noise.
Adding frequency-dependent phase terms to the injected maps has no effect on the intensity mapping results, as the data are first cross-correlated in the frequency domain and the phase term cancels out. However, in the amplitudephase case this adds another level of complexity, as we need to keep track of the individual mode phases as well as those in each pixel. In fact, each frequency bin is independent, given the assumption of time stationarity, hence bandintegrating over frequency bins would then have a similar effect as attempting map-making with a detector resolution that is larger than the intrinsic coherence scale of the field: the signal is degraded as the recovered amplitude averages to zero with the number of bins.

VI. DISCUSSION
We have taken a pedagogical approach in discussing the advantages and shortcomings of two different maximumlikelihood map-making techniques-intensity mapping and amplitude-phase mapping. The intent was to demystify these techniques and clarify the potential of detection and mapping of GWBs in general. Intensity mapping targets the GWB intensity as a function of direction. It constrains the average GW power on the sky and any anisotropies about this average. Versions of the method presented here have been applied to LIGO-Virgo data and have been proposed for LISA [1,2,12,31]. This method is best suited to astrophysical stochastic backgrounds, which carry no phase information and hence solving for it would unnecessarily complicate the mapping procedure. On the other hand, amplitude-phase mapping allows one to do just that; however, we have shown that it is only possible to resolve the phase component if the intrinsic angular coherence scale of the signal is comparable to or larger than the angular resolution of the detector. When this is not the case, attempting to estimate the phase on the sky leads to loss of information, including the amplitude of the signal. Hence, this method is to be avoided in the case of stochastic backgrounds, where the intrinsic coherence scale is usually zero.
A possible exception to this rule may be the case of a stochastic background dominated by very few sources, such as that considered in [10]. This is a possibility for the stochastic background probed by pulsar timing arrays (PTAs), as the signal should be dominated by ∼ 10 2 − 10 3 sources on the sky, such that a single source should dominate a pixel of ∼ 40 deg 2 , which may be taken as the intrinsic coherence scale of the background. Even in this case, we are not quite there yet: achieving such an angular resolution with signal to noise ratio SNR∼3 would require a pulsar array with almost 10 4 pulsars [11], while current PTAs are monitoring ∼ 10 2 pulsars, setting the present resolution to ∼ 400 deg 2 . Hence, this remains a target for future observatories, such as the Square Kilometre Array.
We have provided a useful investigation into the relationship between the time-frequency and the pixel-angular frequency domains. Time samples t are conjugate to frequencies f , while pixels p which correspond to sky locationsn are conjugate to the angular "frequency" scale l. However, our sky observations are effectively carried out in the time domain. A detector naturally low-pass filters the data, both in terms of temporal resolution and angular resolution; however, while the temporal resolution is well-known and controlled by the experimenter, the angular resolution is determined by the geometry and motion of the detector array. Furthermore, the digitization and sampling of the data needs to respect the maximum temporal and spatial frequencies, if aliasing of power is to be avoided.
The first detection of a GWB is drawing near, whether it be with PTAs, possibly by the NANOGrav Collaboration [32], or with the LIGO-Virgo-Kagra ground-based interferometer network. Currently, searches are focusing on the isotropic background, as astrophysical backgrounds are not expected to be highly anisotropic [17,33]. However, understanding how to handle any anisotropy to confirm or rule out this hypothesis remains an essential task of our collaborations. Ideally, the isotropic GWB could be estimated as the monopole of GWB maps [as in Ref. 21, for example], and not under the strict assumption of zero anisotropy.