© ESO 2020

1. Introduction

In the study of cosmic microwave background (CMB) anisotropies, the largest signal is the dipole. This is mainly due to our local motion with respect to the CMB rest frame and it has been previously measured in Kogut et al. (1993), Fixsen et al. (1996), and Hinshaw et al. (2009), and most recently in Planck Collaboration I (2020), Planck Collaboration II (2020), and Planck Collaboration III (2020). Taking the large dipole as being solely caused by our motion, the velocity is v = (369.82 ± 0.11) km s⁻¹ in the direction (l, b) = (264$\stackrel{°}{.}021\pm$ 0$\stackrel{°}{.}011,48\stackrel{°}{.}253\pm 0\stackrel{°}{.}005)$ (Planck Collaboration I 2020). A velocity boost has secondary effects, such as aberration and a frequency-dependent dipolar-modulation of the CMB anisotropies (Challinor & van Leeuwen 2002; Burles & Rappaport 2006). These two effects were first measured using Planck¹ data, as described in Planck Collaboration XXVII (2014). The frequency-dependent part of the dipolar-modulation signal, however, is agnostic to the source of the large CMB dipole. Therefore, its measurement is an independent determination of the CMB dipole. While it may be tempting to use this measure to detect an intrinsic dipole, it has been shown that an intrinsic dipole and a dipole induced by a velocity boost would have the same dipolar-modulation signature on the sky (Challinor & van Leeuwen 2002; Notari & Quartin 2015).

In Notari & Quartin (2015), it was pointed out that the frequency dependence of the dipolar modulation signal could be exploited to achieve a detection with stronger significance than that of Planck Collaboration XXVII (2014). The signal comes from a frequency derivative of the CMB anisotropies’ frequency function and, thus, has essentially the same frequency dependence as the thermal Sunyaev-Zeldovich (tSZ) effect. The dipole-induced quadropule, or “kinematic quadrupole”, would also have the same frequency dependence as the tSZ effect (Kamionkowski & Knox 2003); however, as the quadrupole is less well constrained, a significant detection was not made in this study. Therefore a map of the tSZ effect must contain a copy of the dipole-modulated CMB anisotropies, so an appropriate cross-correlation of the CMB anisotropies with the tSZ effect would be able to pull out the signal. In principle, it could also contribute a bias and source of noise in the bispectrum of the tSZ effect. This is potentially important because the tSZ effect is highly non-Gaussian and much of its information content lies in the bispectrum (Rubiño-Martín & Sunyaev 2003; Bhattacharya et al. 2012; Planck Collaboration XXI 2014; Planck Collaboration XXII 2016).

In this paper we further investigate the CMB under a boost, including tSZ effects (Chluba et al. 2005; Notari & Quartin 2015). We explicitly measure the dipole using a harmonic-space-based method, similar to that outlined in Notari & Quartin (2015), and using a new map-space-based analysis, with consistent results. We also estimate the contamination in the tSZ bispectrum, finding that it is a negligible source of noise.

The structure of the paper is as follows. We describe the nature of the signal we are looking for in Sect. 2. The data that we use, including the choice of CMB maps, tSZ maps, and masks, are described in Sect. 3. The analysis is presented in Sect. 4, and the results in Sect. 5, separately for the multipole-based and map-based methods. We briefly discuss some potential systematic effects in Sect. 6 and we conclude in Sect. 7. We discuss the issues related to the tSZ bispectrum in Appendix A, the results coming from the use of an alternative tSZ map in Appendix B, and how the results could become stronger if we used a less conservative multipole cut in Appendix C.

2. Signal

Here we derive the signal we are looking for. First let us introduce some useful definitions:

$$\begin{array}{cc}& x\equiv \frac{h\nu }{k_{\mathrm{B}}T};\hfill \end{array}$$(1)

$$\begin{array}{cc}\hfill & I\equiv \frac{2k_{\mathrm{B}}^3{T}^3}{h^2{c}^2}\frac{x^3}{e^x-1};\hfill \end{array}$$(2)

$$\begin{array}{cc}\hfill & f(x)\equiv \frac{x{e}^x}{e^x-1};\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill & Y(x)\equiv x\frac{e^x+1}{e^x-1}-4.\hfill \end{array}$$(4)

These are the dimensionless frequency, the Planck blackbody intensity function, the frequency dependence of the CMB anisotropies, and the relative frequency dependence of the tSZ effect, respectively. Here h is Planck’s constant, k_B is the Boltzmann constant, and c is the speed of light.

To first order, anisotropies of intensity take the form

$$\begin{array}{c}\hfill \delta I(\widehat{\mathit{n}})=If(x)[\frac{\delta T(\widehat{\mathit{n}})}{T_0}+y(\widehat{\mathit{n}})Y(x)],\end{array}$$(5)

where the first term represents the CMB anisotropies², and the second term is the tSZ contribution, entering with a different frequency dependence and parameterized by the Compton y-parameter,

$$\begin{array}{c}\hfill y=\int n_{\mathrm{e}}\frac{k_{\mathrm{B}}{\sigma }_{\mathrm{T}}{T}_{\mathrm{e}}}{m_{\mathrm{e}}{c}^2}\mathrm{d}s.\end{array}$$(6)

Here m_e the electron mass, σ_T the Thomson cross-section, ds the differential distance along the line of sight $\widehat{\mathit{n}}$, and n_e and T_e are the electron number density and temperature. Next we apply a Lorentz boost (β ≡ v/c) from the unprimed CMB frame into the primed observation or solar-system frame to obtain

$$\begin{array}{c}\hfill \delta I^{\prime }({\widehat{\mathit{n}}}^{\prime })=I^{\prime }f(x^{\prime })[\frac{\delta T^{\prime }({\widehat{\mathit{n}}}^{\prime })}{T_0^{\prime }}+y^{\prime }({\widehat{\mathit{n}}}^{\prime })Y(x^{\prime })].\end{array}$$(7)

Here $T_0^{\prime }$ is the new boosted blackbody temperature, and only differs from T₀ to lowest order by β²,

$$\begin{array}{c}\hfill T_0^{\prime }=T_0+\frac{{\beta }^2}{2}{T}_0;\end{array}$$(8)

thus to first order $T_0^{\prime }=T_0$. Taking each piece in turn, this transforms as (to first order)

$$\begin{array}{cc}& I^{\prime }f(x^{\prime })=If(x)[1+\beta \mu (Y(x)+3\beta \mu )],\hfill \end{array}$$(9)

$$\begin{array}{cc}\hfill & \frac{\delta T^{\prime }({\widehat{\mathit{n}}}^{\prime })}{T_0^{\prime }}=\frac{\delta T({\widehat{\mathit{n}}}^{\prime })}{T_0}+\beta \mu ,\hfill \end{array}$$(10)

$$\begin{array}{cc}\hfill & y^{\prime }({\widehat{\mathit{n}}}^{\prime })Y(x^{\prime })=y({\widehat{\mathit{n}}}^{\prime })[Y(x)-\beta \mu x\frac{\mathrm{d}Y(x)}{\mathrm{d}x}],\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill & {\widehat{\mathit{n}}}^{\prime }=\widehat{\mathit{n}}-\mathrm{\nabla }(\widehat{\mathit{n}}·\mathit{\beta }),\hfill \end{array}$$(12)

where μ = cos θ, and θ is defined as the angle from the direction β to the line of sight.

Equation (7) can thus be written as³

$$\begin{array}{cc}\hfill \delta I^{\prime }({\widehat{\mathit{n}}}^{\prime })=& If(x)(1+\beta \mu Y(x)+3\beta \mu )\hfill \\ \hfill & \times (\frac{\delta T({\widehat{\mathit{n}}}^{\prime })}{T_0}+\beta \mu +y({\widehat{\mathit{n}}}^{\prime })Y(x)-y({\widehat{\mathit{n}}}^{\prime })\beta \mu x\frac{\mathrm{d}Y(x)}{\mathrm{d}x})\hfill \end{array}$$(13)

or more explicitly, to first order in β,

$$\begin{array}{cc}\hfill \frac{\delta I^{\prime }({\widehat{\mathit{n}}}^{\prime })}{If(x)}=& \frac{\delta T({\widehat{\mathit{n}}}^{\prime })}{T_0}+\beta \mu [1+3\frac{\delta T({\widehat{\mathit{n}}}^{\prime })}{T_0}]\hfill \\ \hfill & +Y(x)[y({\widehat{\mathit{n}}}^{\prime })+\beta \mu \frac{\delta T({\widehat{\mathit{n}}}^{\prime })}{T_0}]\hfill \\ \hfill & +\beta \mu y({\widehat{\mathit{n}}}^{\prime })[3Y(x)+Y^2(x)-x\frac{\mathrm{d}Y(x)}{\mathrm{d}x}],\hfill \end{array}$$(14)

where we have split up each line on the right-hand side according to the frequency dependence. Assuming perfect component separation, and comparing with Eq. (5), the first line of Eq. (14) shows that the boost induces a pure dipole (βμ), an aberration effect ($\delta T({\widehat{\mathit{n}}}^{\prime })/T_0-\delta T(\widehat{\mathit{n}})/T_0$), and a dipolar modulation (3βμδT/T₀) of the CMB. The first effect is the classical CMB dipole, which has been measured many times most recently by Planck (Planck Collaboration I 2020), with the highest accuracy so far achieved. The effects of aberration and dipolar modulation (both frequency-independent and frequency-dependent parts) were measured in Planck Collaboration XXVII (2014) at a combined significance level of 5σ.

In the second line of Eq. (14) we see that the boost also induces a change in a map of the tSZ effect. The original y signal is aberrated ($y({\widehat{\mathit{n}}}^{\prime })-y(\widehat{\mathit{n}})$) and also gains a contribution from the dipolar modulated CMB (βμδT/T₀). This last effect is what we measure in this paper for the first time. Its expected signal can be seen in Fig. 1 (top right panel), along with the full y map obtained via the MILCA method (Fig. 1 top left Hurier et al. 2013). It is worth noting that although the contribution to the tSZ map is a dipolar modulation of the CMB anisotropies, this induces power in the y map that is modulated like a quadrupolar pattern (due to the lack of correlation between the CMB anisotropies and y signal). That is, a y map contains more power in the poles of the dipole, relative to the corresponding equator (see Fig. 1). It should be possible to pull out this signal compared with modulation patterns oriented in orthogonal directions (lower panels of Fig. 1).

Fig. 1. Map of the tSZ effect from the MILCA component-separation method in y-map units (top left) and the expected modulated CMB signal (top right) generated using the SMICA-NOSZ CMB map in units of T0. The bottom left and right figures are the CMB anisotropies modulated in orthogonal directions to the CMB dipole, as explained in the text after Eq. (17). The greyed out region shows the mask used for our analysis. NILC and 2D-ILC y maps, as well as 2D-ILC CMB modulated anisotropies, are not shown, since they look essentially the same as the maps presented here. Note that the map of the tSZ effect (top left) has a different scale bar when compared to the other three (i.e., the modulation signal is about 50 times weaker).

We note that the final line in Eq. (14) is simply the dipole modulation of the tSZ effect, with a peculiar frequency dependence Chluba et al. (2005). In principle one could generate a map of the anisotropies in this new frequency dependence and use the known CMB dipole to measure the y anisotropies again. Such a measurement would be correlated with the original y map, but would have independent noise properties and also have a very low amplitude. From a practical perspective it is unlikely that such a measurement would yield any significant information increase, since the signal is contaminated with relativistic tSZ and kSZ effects, and is suppressed by a factor of β (Chluba et al. 2005).

3. Data

We look for the signal by cross-correlating a template map derived from the CMB temperature data with a y map. Therefore it is important that the CMB map is free of y residuals and that the y map is free of CMB residuals, in order to avoid spurious correlations.

To this end, we use the so-called 2D-ILC CMB temperature map (first used for kSZ detection⁴ in Planck Collaboration Int. XIII 2014), which was produced by the “Constrained ILC” component-separation method designed by Remazeilles et al. (2011) to explicitly null out the contribution from the y-type spectral distortions in the CMB map. We also use the SMICA-NOSZ temperature map, similarly produced with express intent of removing the y-type spectral distortions, and which was generated with the Planck 2018 data release (Planck Collaboration IV 2020). Likewise, we use the corresponding 2D-ILC y map, and the Planck MILCA y map, which explicitly null out the contributions from a (differential) blackbody spectral distribution in the y map (Hurier et al. 2013; Planck Collaboration XXII 2016). We also consider the Planck NILC y map, which does not explicitly null out the blackbody contribution. The extra constraint to remove CMB anisotropies in the 2D-ILC y map, or in the MILCA y map, is at the expense of leaving more contamination by diffuse foregrounds and noise. In Planck Collaboration XXII (2016), the Planck NILC y map was preferred over the 2D-ILC y map to measure the angular power spectrum of tSZ anisotropies, the CMB contaminant being negligible compared to diffuse foregrounds. Conversely, to measure the dipole modulation of the CMB anisotropies in the y map, the 2D-ILC and MILCA y maps are preferred over the Planck NILC y map because the latter does not fully null out the contribution from the CMB. This significantly contaminates the signal we are looking for, as can be seen in Appendix B. It is noted in Planck Collaboration XXII (2016) that for ℓ >  2000 the signal is dominated by correlated noise, and so we use the same cut as used in their analysis of ℓ_max = 1411, this is further justified in Sect. 4.2.

Figure 1 shows the mask used in our analysis. This is the union of the Planck 2018 data release common temperature confidence mask (Planck Collaboration IV 2020), and the corresponding y-map foreground masks (Planck Collaboration XXII 2016). This was then extended by 1° and apodized with a 200′ Gaussian beam. To account for any masked sections lost during the smoothing, the original mask was then added back. Tests were also done using the y-map point-source mask, with negligible changes seen in the results, and was thus omitted from the final analysis. This procedure aims to allow for the maximum signal while minimizing the foreground contamination. Various combinations of mask sizes and apodizations were also tested and final results were consistent, independent of the choice of mask.

4. Analysis

From Eq. (14) we see that a map of the tSZ effect (ℳ_SZ) contains the following terms (for each pixel, or direction $\widehat{\mathit{n}}$):

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{SZ}}=y+{\eta }^y+\beta \mu \frac{\delta T}{T_0},\end{array}$$(15)

where η^y is simply the noise in the y map, and we have neglected the aberration effect. Our goal is to isolate the final term in Eq. (15), which we do via a suitable cross-correlation with a CMB map. A map of the CMB (ℳ_CMB) contains the following terms:

$$\begin{array}{c}\hfill {\mathcal{M}}_{\mathrm{CMB}}=\frac{\delta T}{T_0}+{\eta }^T+3\beta \mu \frac{\delta T}{T_0},\end{array}$$(16)

where we have explicitly removed the full dipole term, and η^T is the noise in the CMB map. If we multiply our CMB map (Eq. (16)) with βμ and cross-correlate that with our tSZ map (Eq. (15)), then we can directly probe the dipole modulation. This of course neglects the noise and modulation terms in the CMB map, which we are justified in doing because the noise term is sub-dominant, except at very small scales (we make the restriction ℓ_max = 1411, so that the y map and CMB maps are still signal dominated Planck Collaboration XXII 2016), and because the modulation term becomes second order in β. Equivalently one could directly cross-correlate Eq. (15) with Eq. (16) and look for the signal in harmonic space from the coupling of ℓ and ℓ ± 1 modes.

In Planck Collaboration XXVII (2014) a quadratic estimator was used to determine the dipole aberration and modulation, in essence using the auto-correlation of the CMB fluctuation temperature maps, weighted appropriately to extract the dipole signal. The auto-correlation naturally introduces a correlated noise term, which must be well understood for this method to work. In this paper we take advantage of the fact that we know the true CMB fluctuations with excellent precision and therefore the signal that should be present in the y map. We can therefore exploit the full angular dependence of the modulation signal and remove much of the cosmic variance that would be present in the auto-correlation.

In order to implement this idea we define three templates, B_i (with i = 1, 2, 3) as

$$\begin{array}{c}\hfill B_i(\widehat{\mathit{n}})=\beta \widehat{\mathit{n}}·{\widehat{\mathit{m}}}_i\frac{\delta T}{T_0}(\widehat{\mathit{n}}),\end{array}$$(17)

where β = v/c is 1.23357 × 10⁻³ (Planck Collaboration I 2020) and ${\widehat{\mathit{m}}}_1,{\widehat{\mathit{m}}}_2,{\widehat{\mathit{m}}}_3$ are the CMB dipole direction, an orthogonal direction in the Galactic plane, and the third remaining orthogonal direction (see Fig. 1 and Planck Collaboration XXVII 2014, for a similar approach). Note that in Eq. (17), we simply use our CMB map in place of δT/T₀. We use two distinct methods to accomplish this, discussed in detail in Sects. 4.1 and 4.2.

In the region where the CMB is signal dominated we can regard δT/T₀ as fixed, and thus our templates B_i are fixed. Due to the presence of the CMB dipole, the signal B₁ should be present in the y map. We can therefore directly cross-correlate B₁ with our y map (Eq. (15)) to pull out the signal. Likewise, the cross-correlation of B₂ and B₃ with our y map should give results consistent with noise, although the coupling with the noise and mask leads to a bias that is recovered through simulations.

Our y simulations are generated by first computing the power spectra of our data y maps; specifically we apply the MASTER method using the NaMASTER routine (Alonso et al. 2019) to account for the applied mask (Hivon et al. 2002). Then we generate y maps using this power-spectrum with the HEALPix (Górski et al. 2005) routine synfast⁵. This is done separately for the 2D-ILC and MILCA maps because they have different noise properties (and thus different total power spectra). For our simulations that include the dipolar modulated CMB anisotropies we add the last term of Eq. (15). We finally apply a Gaussian smoothing of 5′ to model the telescope beam.

For each analysis method (to be described in the following subsections) we estimate the amplitude of the dipole (${\widehat{\beta }}_i$) in each of the three orthogonal directions⁶. We apply the same analysis on a suite of 1000 y simulations, generated with and without the dipolar modulation term in Eq. (15). We are then able to generate a covariance that appropriately contains the effects of the mask we use and are able to compute any bias that the mask induces. On the assumption (verified by our simulations) that the ${\widehat{\beta }}_i$ estimators are Gaussian, we are able to compute a value of χ² for the case of no CMB term and with the CMB term (see Table 1). We can then apply Bayes’ theorem along with our covariance to calculate the probability that each model is true (with or without the CMB term) and the posterior of our dipole parameters (β, l, b), summarized in Table 1 and Fig. 2.

Fig. 2. Posteriors for the CMB dipole parameters using the map-space analysis (top) and harmonic-space analysis (bottom). The left panels use the SMICA-NOSZ CMB maps, whereas the right use the 2D-ILC CMB maps. MILCA y map results are shown in orange, and 2D-ILC y map results are shown in blue. Black dashed lines show the best-fit parameters from direct measurements of the CMB dipole. Dark and light contours enclose 68% and 95%, respectively. Titles for each panel give the best-fit results, along with the 68% uncertainties.

We estimate the covariance C_ij of the ${\widehat{\beta }}_i$ using the simulations⁷ and calculate the χ² as

$$\begin{array}{c}\hfill {\chi }_k^2={\displaystyle \sum _{\mathit{ij}}({\widehat{\beta }}_i-{⟨{\widehat{\beta }}_i⟩}_k)C_{\mathit{ij}}^{-1}({\widehat{\beta }}_j-{⟨{\widehat{\beta }}_j⟩}_k),}\end{array}$$(18)

where k denotes whether the expectation value in the sum is taken over the simulations that do or do not include the CMB term. For definiteness we define the null hypothesis H₀ (k = 0) to not include the CMB term, while hypothesis H₁ (k = 1) does include the CMB term. We can then directly calculate the probability that H_k is true given the data (${\widehat{\beta }}_i$) as

$$\begin{array}{cc}& P(H_k|{\widehat{\beta }}_i)=P({\widehat{\beta }}_i|H_k)P(H_k),\hfill \end{array}$$(19)

$$\begin{array}{cc}& P({\widehat{\beta }}_i|H_k)=\frac{1}{\sqrt{|2\pi C|}}{e}^{-{\chi }_k^2/2}.\hfill \end{array}$$(20)

We can calculate the odds ratio, O₁₀, on the assumption that the two hypotheses are equally likely,

$$\begin{array}{c}\hfill O_{10}\equiv \frac{P(H_1|{\widehat{\beta }}_i)}{P(H_0|{\widehat{\beta }}_i)}=\frac{e^{-{\chi }_1^2/2}}{e^{-{\chi }_0^2/2}}·\end{array}$$(21)

This quantity tells us to what degree H₁ should be trusted over H₀. Assuming that the two hypotheses are exhaustive, it is directly related to the probability that the individual hypotheses are true:

$$\begin{array}{cc}& P(H_0|{\widehat{\beta }}_i)=\frac{1}{1+O_{10}},\hfill \end{array}$$(22)

$$\begin{array}{cc}& P(H_1|{\widehat{\beta }}_i)=\frac{O_{10}}{1+O_{10}}·\hfill \end{array}$$(23)

These quantities and the ${\chi }_k^2$ values are given in Table 1.

We can also generate a likelihood for our parameters with the same covariance matrix:

$$\begin{array}{c}\hfill \mathcal{L}({\beta }_i)=\frac{1}{\sqrt{|2\pi C|}}{e}^{-{\chi }^2/2},\end{array}$$(24)

where we define the modified χ² as above. We can then apply Bayes’ theorem with uniform priors on the β_i, equating the posterior of β_i with Eq. (24). A simple conversion allows us to obtain the posterior of the parameters in spherical coordinates $(\stackrel{\sim }{\beta },\stackrel{\sim }{l},\stackrel{\sim }{b})$. We show this in Fig. 2 for our two analyses, using the 2D-ILC and MILCA maps.

In the following subsections we describe two methods of cross-correlation: the first we perform directly in map-space; and the second is performed in harmonic space. An advantage of using two independent methods is that their noise properties are different; for example, working in harmonic space introduces complications with masking, whereas in map space, although it may not be clear how to optimally weight the data, the estimator has less sensitivity to large-scale systematic effects. Thus, the advantage of using two approaches will become apparent when we try to assess the level of systematic error in our analysis.

4.1. Map-space method

First we apply our mask to the templates B_i and y map. Then we locate all peaks (i.e., local maxima or minima) of the template map B_i and select a patch of radius $2\stackrel{°}{.}0$ around each peak. Our specific implementation of the peak method follows earlier studies, for example Planck Collaboration VII (2020). The weighting scheme has not been shown to be optimal, but a similar approach was used for determining constraints on cosmic birefringence Contreras et al. (2017) and gave similar results to using the power spectra (and the issue of weighting is further discussed in Jow et al. 2019). Intuitively we would expect that sharper peaks have a higher signal, and hence that influences our choice for the weighting scheme described below. For every peak we obtain an estimate of ${\widehat{\beta }}_i$ by the simple operation

$$\begin{array}{c}\hfill {\widehat{\beta }}_{i,\mathrm{p}}=\beta \frac{{\sum }_{k\in D(p)}{B}_{i,k}{y}_k}{{\sum }_{k\in D(p)}{B}_{i,k}^2},\end{array}$$(25)

where D(p) is the collection of all unmasked pixels in a $2\stackrel{°}{.}0$ radius centred on pixel p, and p is the position of a peak. Equation (25) is simply a cross-correlation in map space and by itself offers a highly-noisy (and largely unbiased⁸) estimate.

We then combine all individual peak estimates with a set of weights (w_p) to give our full estimate:

$$\begin{array}{c}\hfill {\widehat{\beta }}_i=\frac{{\sum }_{\mathrm{p}}{w}_{i,\mathrm{p}}{\widehat{\beta }}_{i,\mathrm{p}}}{{\sum }_{\mathrm{p}}{w}_{i,\mathrm{p}}}·\end{array}$$(26)

The values of w_i, p depend solely on the templates B_i, and they can be chosen to obtain the smallest uncertainties. We choose w_p to be proportional to the square of the dipole, which ensures that peaks near the dipole direction (and anti-direction) are weighted more than those close to the corresponding equator. We further choose that the weights are proportional to the square of the Laplacian at the peak (Desjacques 2008); this favours sharply defined peaks over shallow ones. Finally we account for the scan strategy of the Planck mission by weighting by the 217 GHz hits map (denoted $H_{\mathrm{p}}^{217}$ Planck Collaboration VIII 2016), though this choice provides no appreciable difference to our results. The weights then are explicitly

$$\begin{array}{c}\hfill w_{i,\mathrm{p}}={|\widehat{\mathit{n}}·{\widehat{\mathit{m}}}_i|}_{\mathrm{p}}^2{\left({{\mathrm{\nabla }}^2(B_i)|}_{\mathrm{p}}\right)}^2{H}_{\mathrm{p}}^{217}.\end{array}$$(27)

We evaluate the Laplacian numerically in pixel space at pixel p. The weighting scheme closely resembles the bias factors that come about when relating peaks to temperature fluctuations (Bond & Efstathiou 1987), as used in Komatsu et al. (2011), Planck Collaboration Int. XLIX (2016) and Jow et al. (2019).

Combining Eqs. (26) and (27) gives us our estimates, ${\widehat{\beta }}_i$. We apply the method for each of our simulated y maps, in exactly the same way as for the data.

4.2. Harmonic-space method

The alternative approach is to directly cross-correlate Eq. (17) with the y map, and compare this to the auto-correlation of Eq. (17).

Our first step is identical to the previous method in that we mask the templates B_i and y maps. Under the assumption that the y map contains the template (B_i), the y multipoles are Gaussian random numbers with mean and variance given by

$$\begin{array}{cc}& s_{\ell m}^i={\displaystyle \int \mathrm{d}\mathrm{\Omega }\beta }{\widehat{\mathit{m}}}_i·\widehat{\mathit{n}}\frac{\delta T}{T_0}M(\mathrm{\Omega })Y_{\ell m}^{\ast },\hfill \end{array}$$(28)

$$\begin{array}{cc}\hfill & {\sigma }_{\ell }^2=C_{\ell }^y+N_{\ell }^y,\hfill \end{array}$$(29)

respectively, where M(Ω) is the mask over the sphere, Y_ℓm are the spherical harmonics, and the ${\widehat{\mathit{m}}}_i$ are as defined in Eq. (17). Thus we can obtain an estimate of β_i by taking the cross-correlation with inverse-variance weighting. We can demonstrate this simply by writing our y map as a sum of our expected signal plus everything else⁹,

$$\begin{array}{c}\hfill y_{\ell m}=\frac{{\beta }_i}{\beta }{s}_{\ell m}^i+{\eta }_{\ell m}^y.\end{array}$$(30)

Here our signal is of course given when β_i = βδ_1i and all sources of noise, such as tSZ, are given by ${\eta }_{\mathcal{l}m}^y$. We then cross-correlate with our template and sum over all multipoles with inverse-variance weighting. We explicitly consider noise in our template, that is our template (${\overline{s}}_{\ell m}^i$) is related to Eq. (29) via, ${\overline{s}}_{\ell m}^i=s_{\ell m}^i+{\eta }_{\ell m}^t$, where ${\eta }_{\mathcal{l}m}^t$ is the noise in our template. Then the cross-correlation looks like,

$$\begin{array}{cc}\hfill {\displaystyle \sum _{\ell m}}{\overline{s}}_{\ell m}^{i^{\prime }}{y}_{\ell m}^{\ast }/{\sigma }_{\ell }^2=& \frac{{\beta }_i}{\beta }{\displaystyle \sum _{\ell m}}{\overline{s}}_{\ell m}^{i^{\prime }}{(s_{\ell m}^i)}^{\ast }/{\sigma }_{\ell }^2\hfill \\ \hfill & +\sum _{\ell m}{\overline{s}}_{\ell m}^{i^{\prime }}{({\eta }_{\ell m}^y)}^{\ast }/{\sigma }_{\ell }^2,\hfill \end{array}$$(31)

and expanding it out, this becomes

$$\begin{array}{cc}\hfill {\displaystyle \sum _{\ell m}}{s}_{\ell m}^{i^{\prime }}{y}_{\ell m}^{\ast }/{\sigma }_{\ell }^2+{\displaystyle \sum _{\ell m}}{\eta }_{\ell m}^t{y}_{\ell m}^{\ast }/{\sigma }_{\ell }^2=& \frac{{\beta }_i}{\beta }{\displaystyle \sum _{\ell m}}{s}_{\ell m}^{i^{\prime }}{(s_{\ell m}^i)}^{\ast }/{\sigma }_{\ell }^2\hfill \\ \hfill & +\frac{{\beta }_i}{\beta }{\displaystyle \sum _{\ell m}}{\eta }_{\ell m}^t{(s_{\ell m}^i)}^{\ast }/{\sigma }_{\ell }^2\hfill \\ \hfill & +{\displaystyle \sum _{\ell m}}{s}_{\ell m}^{i^{\prime }}{({\eta }_{\ell m}^y)}^{\ast }/{\sigma }_{\ell }^2\hfill \\ \hfill & +{\displaystyle \sum _{\ell m}}{\eta }_{\ell m}^t{({\eta }_{\ell m}^y)}^{\ast }/{\sigma }_{\ell }^2.\hfill \end{array}$$(32)

The last term on the left and last three terms on the right are all statistically zero, since our template does not correlate with tSZ or any other types of noise and the noise in our template does not correlate with the template itself or with noise in the y map (by assumption). Hence we can solve for β_i neglecting those terms, to produce our estimator ${\widehat{\beta }}_i$:

$$\begin{array}{c}\hfill {\widehat{\beta }}_i=\beta {\displaystyle \sum _{i^{\prime }}{[\sum _{\ell m}^{{\ell }_{max}}{s}_{\ell m}^i{(s_{\ell m}^{i^{\prime }})}^{\ast }/{\sigma }_{\ell }^2]}^{-1}\sum _{\ell m}^{{\ell }_{max}}{s}_{\ell m}^{i^{\prime }}{(y_{\ell m})}^{\ast }/{\sigma }_{\ell }^2.}\end{array}$$(33)

It is important to note that in practice we do not have $s_{\mathcal{l}m}^i$, since we do not know the exact realization of noise in the CMB, so we instead use ${\overline{s}}_{\ell m}^i$. Using Weiner-filtered results would allow us to calculate $s_{\mathcal{l}m}^i$, but adds complexity in the masking process. We can compare what the Weiner-filtered results would be,

$$\begin{array}{c}\hfill {\displaystyle \sum _{\ell }\frac{{(C_{\ell }^{\mathit{TT}})}^2}{C_{\ell }^{\mathit{TT}}+N_{\ell }^{\mathit{TT}}}\frac{(2\ell +1)}{{\sigma }_{\ell }^2},}\end{array}$$(34)

to our results,

$$\begin{array}{c}\hfill {\displaystyle \sum _{\ell }(C_{\ell }^{\mathit{TT}}+N_{\ell }^{\mathit{TT}})}\frac{(2\ell +1)}{{\sigma }_{\ell }^2}\end{array}$$(35)

and find the bias to be on the order of 2% for ℓ_max = 1411, justifying our use of the cut-off. Equation (33) is in fact a direct solution for β_i in the absence of noise, since it is the direct solution of Eq. (32) in the absence of noise.

Relating back to the map-space method, $s_{\mathcal{l}m}^i$ are the spherical harmonic coefficients of the templates denoted previously by B_i, and y_ℓm are the spherical harmonic coefficients of the y map. The values for $s_{\ell m}^i{(s_{\ell ,m}^{i^{\prime }})}^{\ast }$ and $s_{\mathcal{l}m}^i{(y_{\mathcal{l}m})}^{*}$ may be computed using the maps with the HEALPix (Górski et al. 2005) routine anafast. In the case of the term $s_{\ell m}^i{(s_{\ell ,m}^{i^{\prime }})}^{\ast }$ this results in a 3 × 3 matrix for each ℓ, with the cross-power spectrum for the three templates on the off diagonals.

In the absence of a mask M the signal |s_ℓm|² induces power in a cos²θ pattern. The presence of a mask (being largely quadrupolar in shape) induces power in a more complicated way, but has strong overlap with a cos²θ pattern as well. Therefore the application of a mask necessarily makes this method sub-optimal; however, since the template is masked in the same way, the method is unbiased.

We apply the method for each of our simulated y maps, in exactly the same way as for the data, in order to assess whether the dipole modulation is detected.

5. Results

The main results of this paper are presented in Tables 1 and 2 and Fig. 2. They show how consistent the data are with the presence (or non-presence) of the dipole term, and the recovered posteriors of the dipole parameters, respectively. In the following subsections we describe our results for each method in more detail.

5.1. Map-space results

First we compare the consistency of the data with our two sets of simulations (with and without the dipole term). This comparison shown in Fig. 3, with blue histograms being the simulations with the dipole term and orange histograms without. The data (black line) for 2D-ILC and MILCA can clearly be seen to be consistent with the simulations with the dipole term; this observation is made quantitative from examination of the χ² (see Tables 1 and 2). The map-space method is more susceptible to biases induced by the mask, particularly in the off-dipole directions, $\widehat{{\beta }_2}$ and $\widehat{{\beta }_3}$; this is due to subtle correlations between the mask and templates, but has only a small effect in those directions (at the level of a few tenths of σ), as can be seen in Fig. 3. Converted into the equivalent probabilities for Gaussian statistics, we can say that the dipole modulation is detected at the 5.0 to 5.9σ level.

Fig. 3. Histograms of ${\widehat{\beta }}_i/\beta$ values (with 1, 2, and 3 corresponding to the CMB dipole direction, Galactic plane, and a third orthogonal direction) using the map-space analysis for MILCA (top) and 2D-ILC (bottom) y maps, and for CMB template maps SMICA-NOSZ (left) and 2D-ILC (right). Blue histograms are simulations with the dipolar modulation term, and orange histograms are simulations without. Black vertical lines denote the values of the data, demonstrating that they are much more consistent with the existence of the dipolar modulation term than without it. Dashed lines show the 68% regions for a Gaussian fit to the histograms.

5.2. Harmonic-space results

Figure 4 is the equivalent of Fig. 3, but for the harmonic-space analysis. Similar to the previous subsection the data are much more consistent with the modulated simulations than the unmodulated simulations. Tables 1 and 2 contain the explicit χ² values and verify this quantitatively. The harmonic-space method is somewhat susceptible to biases induced by the mask, due to the complex coupling that occurs, mainly between the ℓ and ℓ ± 2 modes. This can be seen in the slight bias in the results for ${\widehat{\beta }}_2$ and ${\widehat{\beta }}_3$. Nevertheless, we can say that we confidently detect the dipole modulation at the 6.2 to 6.6σ level.

Fig. 4. As in Fig. 3, except now for the harmonic-space analysis.

6. Systematics

We have generated results using two distinct methods, namely the map-space method and harmonic-space method, with two distinct CMB maps and two distinct y maps, and have shown the results to be consistent with the presence of a dipole-modulation signal in the expected direction. Each test is subject to slightly different systematics, but since the results are consistent, we can conclude that there is likely no significant systematic interfering with the results. Further tests, relaxing the limits of ℓ_max = 1411 show that it is possible to achieve even higher levels of significance using smaller-scale data (see Appendix C). In that sense, the results in this paper are conservative; however, if it becomes possible to construct reliable y maps out to higher multipoles then it should be possible to achieve a detection of the dipole modulation at perhaps twice the number of σ as found here.

6.1. Residuals in the component separation

The NILC y maps are known to contain some remnant CMB contamination, unlike the MILCA and 2D-ILC y maps, which have been generated with the express purpose of eliminating the CMB contribution. This contaminates the signal we are looking for; the results from the NILC y maps may be seen in Appendix B. Any contamination remaining in the MILCA and 2D-ILC y maps is sufficiently low that is does not hide the dipole modulation signal.

6.2. Galactic foregrounds

It is known that the y maps are contaminated by Galactic foregrounds; however, as the results here are from a cross-correlation of the modulated CMB maps with the y maps such contamination does not have a large effect on the results. To further support this, a number of different mask sizes and combinations were tested, with the final mask selected because among those choices consistent with the more conservative masks, it gave the highest signal-to-noise ratio. Larger masks serve only to decrease the signal-to-noise of the data. This suggests that foregrounds have only a small effect on the detection of the dipole modulation. Foregrounds have been mentioned as a potential issue in previous results (Planck Collaboration XXII 2016).

7. Conclusions

Due to the existence of the CMB dipole, a tSZ map necessarily contains a contaminating signal that is simply the dipole modulation of the CMB anisotropies. This occurs because CMB experiments do not directly measure temperature anisotropies, but instead measure intensity variations that are conventionally converted to temperature variations. This contamination adds power to the tSZ map in a Y₂₀ pattern, with its axis parallel to the dipole direction. We have measured this effect and determined a statistically independent value of the CMB dipole, which is consistent with direct measurements of the dipole. Using a conservative multipole cut on the y map, the significance of the detection of the dipole modulation signal is around 5 or 6σ, depending on the precise choice of data set and analysis method. This is a significant improvement from the 2 to 3σ results in Planck Collaboration XXVII (2014). We also find that the contamination of the tSZ map contributes negligible noise to the bispectrum calculations (see Appendix A).

Acknowledgments

The Planck Collaboration acknowledges the support of: ESA; CNES and CNRS/INSU-IN2P3-INP (France); ASI, CNR, and INAF (Italy); NASA and DoE (USA); STFC and UKSA (UK); CSIC, MINECO, JA, and RES (Spain); Tekes, AoF, and CSC (Finland); DLR and MPG (Germany); CSA (Canada); DTU Space (Denmark); SER/SSO (Switzerland); RCN (Norway); SFI (Ireland); FCT/MCTES (Portugal); and ERC and PRACE (EU). A description of the Planck Collaboration and a list of its members, indicating which technical or scientific activities they have been involved in, can be found at http://www.cosmos.esa.int/web/planck/planck-collaboration. Some of the results in this paper have been derived using the HEALPix package and the NaMaster package and some plots were generated using the pygtc package.

References

Appendix A: The tSZ bispectrum

Fundamentally the modulation is a correlation between C_ℓ and C_ℓ ± 1. The signal considered here therefore shows up most prominently in the 4-point function (i.e., trispectrum) and thus we do not expect it to bias the measurements of the tSZ bispectrum; however, since the bispectrum is an important quantity for characterizing the tSZ signal, it is worth checking to ensure that the dipolar modulation does not add significant noise. In other words, we want to check if it is important to remove the dipole modulations before performing analysis of the tSZ bispectrum. Lacasa et al. (2012) and Bucher et al. (2010) describe in detail the calculation of the bispectrum and the binned bispectrum, and this is summarized below. The reduced bispectrum is given by

$$\begin{array}{cc}\hfill B_{{\ell }_1{\ell }_2{\ell }_3}=& {(N_{{\ell }_1{\ell }_2{\ell }_3})}^{-1/2}\hfill \\ \hfill & \times {\displaystyle \sum _{m_1{m}_2{m}_3}\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ m_1& m_2& m_3\end{array}\right)a_{{\ell }_1{m}_1}{a}_{{\ell }_2{m}_2}{a}_{{\ell }_3{m}_3},}\hfill \end{array}$$(A.1)

where

$$\begin{array}{c}\hfill \left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ m_1& m_2& m_3\end{array}\right)\end{array}$$

represent the Wigner-3j functions and

$$\begin{array}{c}\hfill N_{{\ell }_1{\ell }_2{\ell }_3}=\frac{(2{\ell }_1+1)(2{\ell }_2+1)(2{\ell }_3+1)}{4\pi }{\left(\begin{array}{ccc}{\ell }_1& {\ell }_2& {\ell }_3\\ 0& 0& 0\end{array}\right)}^2.\end{array}$$(A.2)

The normalized bispectrum is non-zero for terms where m₁ + m₂ + m₃ = 0, |ℓ₂ − ℓ₂|≥ℓ₃ ≥ ℓ₁ + ℓ₂, and ℓ₁ + ℓ₂ + ℓ₃ is even (this is due to the term before the sum with m₁, m₂, m₃ = 0). Typically, the binned bispectrum is analysed to reduce the number of terms calculated and saved, which constitutes only a small loss of information because the bispectrum is expected to vary slowly with ℓ (Lacasa et al. 2012). The data are binned by breaking down the interval from ℓ_min to ℓ_max into i bins, denoted by Δ_i. An average for the bispectrum of a particular bin can then be calculated using

$$\begin{array}{c}\hfill B_{i_1{i}_2{i}_3}=\frac{1}{{\mathrm{\Xi }}_{i_1{i}_2{i}_3}}{\displaystyle \sum _{{\ell }_1\in {\mathrm{\Delta }}_1}\sum _{{\ell }_2\in {\mathrm{\Delta }}_2}\sum _{{\ell }_3\in {\mathrm{\Delta }}_3}{B}_{{\ell }_1{\ell }_2{\ell }_3},}\end{array}$$(A.3)

where Ξ_i1i2i3 is the number of non-zero elements in the given bin. Both the bispectrum and the binned bispectrum may be calculated using an integral over the map space as well, rather than in harmonic space. This is achieved by first generating the binned scalemaps defined by

$$\begin{array}{c}\hfill y_{{\mathrm{\Delta }}_i}(\widehat{\mathit{n}})={\displaystyle \sum _{\ell \in {\mathrm{\Delta }}_i,m}{y}_{\ell m}{Y}_{\ell m}(\widehat{\mathit{n}}),}\end{array}$$(A.4)

where the sum goes from ℓ_min to ℓ_max in the bin Δ_i. We can then use

$$\begin{array}{c}\hfill B_{i_1{i}_2{i}_3}=\frac{1}{N_{i_1{i}_2{i}_3}}{\displaystyle \int {\mathrm{d}}^2}\widehat{\mathit{n}}{y}_{{\mathrm{\Delta }}_1}(\widehat{\mathit{n}})y_{{\mathrm{\Delta }}_2}(\widehat{\mathit{n}})y_{{\mathrm{\Delta }}_3}(\widehat{\mathit{n}})\end{array}$$(A.5)

which gives the weighted average of the bispectrum within the bins (Lacasa et al. 2012).

In Figs. A.1 and A.2 we show a subset of the binned normalized bispectra for the y maps, with and without the dipole modulation. For simplicity, since we are just comparing the results of two simulated maps, there are no non-Gaussianities and no mask applied. This analysis was performed using the MILCA y map and the SMICA-NOSZ CMB temperature map. Plots are constructed in the style suggested by Lacasa et al. (2012) for an ℓ_max of 500 (and an N_side of 512 to speed up computation).

Fig. A.1. Binned bispectrum, for simulated y maps, with ℓmax = 500 and bin sizes of 10, and “scalemaps” defined in Eq. (A.4), of Nside = 512. Left panels: bispectrum for a simulated y map with the dipole modulation; right panels: same with no dipole-modulation. The quantities P, F, and S are as defined in Eq. (A.13).

Fig. A.2. Absolute (left) and relative (right) difference between the bispectrum with and without the dipole-modulation term.

Useful definitions here are

$$\begin{array}{cc}& {\sigma }_1={\ell }_1+{\ell }_2+{\ell }_3,\hfill \end{array}$$(A.6)

$$\begin{array}{cc}\hfill & {\sigma }_2={\ell }_1{\ell }_2+{\ell }_1{\ell }_3+{\ell }_2{\ell }_3,\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill & {\sigma }_3={\ell }_1{\ell }_2{\ell }_3,\hfill \end{array}$$(A.8)

$$\begin{array}{cc}\hfill & {\stackrel{\sim }{\sigma }}_2=12{\sigma }_2/{\sigma }_1^2-3,\in [0,1],\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill & \stackrel{\sim }{{\sigma }_3}=27{\sigma }_3/{\sigma }_1^3,\in [0,1],\hfill \end{array}$$(A.10)

$$\begin{array}{cc}\hfill & F=32({\stackrel{\sim }{\sigma }}_2-{\stackrel{\sim }{\sigma }}_3)/3+1,\hfill \end{array}$$(A.11)

$$\begin{array}{cc}\hfill & S=\stackrel{\sim }{{\sigma }_3},\hfill \end{array}$$(A.12)

$$\begin{array}{cc}\hfill & P={\sigma }_1,\hfill \end{array}$$(A.13)

where P is the perimeter, each plot represents the results of a particular perimeter size, F is plotted along the y-axis of the panels and S is plotted along the x-axis of the panels.

Our main goal is to determine whether the dipole modulation contamination of the y maps is significant, and to what degree it is significant for current and future analysis as data improves. For this purpose a subset of the tested perimeter values are plotted, for data with the dipole modulation and without, and the absolute value of the differences. It does not appear that the dipole modulation has a noticeable effect on the bispectrum results.

Appendix B: NILC y-map results

In creating the Plancky maps using NILC, the choices were optimized for removal of the contamination by CMB, foregrounds, and noise. With the MILCA y maps there was an additional constraint added to fully eliminate the CMB, at the expense of adding more foregrounds and noise contamination. For this reason the CMB contamination in the NILC y maps is too high for us to robustly detect the dipole modulation. The 2D-ILC y map was also produced with the express intent of removing all CMB contamination, and both it and the MILCA maps clearly show that the dipole modulation is present. For completeness, here we present the effect of the contamination in the NILC y maps in Fig. B.1. The dipole modulation signal is seen to be completely hidden by the CMB contamination.

Fig. B.1. As in Figs. 3 (top) and 4 (bottom), except using the NILC y maps. Top panels: results from the map-space method, bottom panels: those from the harmonic-space method.

Appendix C: Increased ℓ_max results

In our analysis for the harmonic-space method the results were truncated at ℓ_max = 1411, since this is the recommendation from Planck Collaboration XXII (2016) to avoid the correlated noise and foreground contaminations present in higher ℓ. If we were to assume that the simulations model the data properly up to a higher ℓ_max, and that we also trust the data up to this higher ℓ_max, then we would be able to achieve a greater significance than reported in the conclusions. This can be seen in the simulation results using the MILCA y map and the SMICA-NOSZ CMB templates. These particular results are from 500 simulations for N_side = 1024, and ℓ_max = 2750. The significance appears to be at the > 12σ level. To do the analysis fully at this ℓ_max the Weiner filter would also need to be applied to the CMB maps, as without it the bias would be much larger than the 2% found in our analysis.

Fig. C.1. As in Fig. 4, but with ℓmax = 2750 compared to ℓmax = 1411 used in the paper. If we were to trust the y map out to these multipoles, then these results would have a significance of > 12σ.

All Tables

All Figures

Fig. 1. Map of the tSZ effect from the MILCA component-separation method in y-map units (top left) and the expected modulated CMB signal (top right) generated using the SMICA-NOSZ CMB map in units of T0. The bottom left and right figures are the CMB anisotropies modulated in orthogonal directions to the CMB dipole, as explained in the text after Eq. (17). The greyed out region shows the mask used for our analysis. NILC and 2D-ILC y maps, as well as 2D-ILC CMB modulated anisotropies, are not shown, since they look essentially the same as the maps presented here. Note that the map of the tSZ effect (top left) has a different scale bar when compared to the other three (i.e., the modulation signal is about 50 times weaker).

In the text

Fig. 2. Posteriors for the CMB dipole parameters using the map-space analysis (top) and harmonic-space analysis (bottom). The left panels use the SMICA-NOSZ CMB maps, whereas the right use the 2D-ILC CMB maps. MILCA y map results are shown in orange, and 2D-ILC y map results are shown in blue. Black dashed lines show the best-fit parameters from direct measurements of the CMB dipole. Dark and light contours enclose 68% and 95%, respectively. Titles for each panel give the best-fit results, along with the 68% uncertainties.

In the text

Fig. 3. Histograms of ${\widehat{\beta }}_i/\beta$ values (with 1, 2, and 3 corresponding to the CMB dipole direction, Galactic plane, and a third orthogonal direction) using the map-space analysis for MILCA (top) and 2D-ILC (bottom) y maps, and for CMB template maps SMICA-NOSZ (left) and 2D-ILC (right). Blue histograms are simulations with the dipolar modulation term, and orange histograms are simulations without. Black vertical lines denote the values of the data, demonstrating that they are much more consistent with the existence of the dipolar modulation term than without it. Dashed lines show the 68% regions for a Gaussian fit to the histograms.

In the text

	Fig. 4. As in Fig. 3, except now for the harmonic-space analysis.
In the text

	Fig. A.1. Binned bispectrum, for simulated y maps, with ℓmax = 500 and bin sizes of 10, and “scalemaps” defined in Eq. (A.4), of Nside = 512. Left panels: bispectrum for a simulated y map with the dipole modulation; right panels: same with no dipole-modulation. The quantities P, F, and S are as defined in Eq. (A.13).
In the text

	Fig. A.2. Absolute (left) and relative (right) difference between the bispectrum with and without the dipole-modulation term.
In the text

	Fig. B.1. As in Figs. 3 (top) and 4 (bottom), except using the NILC y maps. Top panels: results from the map-space method, bottom panels: those from the harmonic-space method.
In the text

	Fig. C.1. As in Fig. 4, but with ℓmax = 2750 compared to ℓmax = 1411 used in the paper. If we were to trust the y map out to these multipoles, then these results would have a significance of > 12σ.
In the text