© ESO 2021

1. Introduction

The peculiar motion of an observer relative to an ideal reference frame at rest with respect to the cosmic background in a given frequency band, produces boosting effects in the anisotropy patterns at low multipoles with frequency spectral behaviours related to the spectrum of the isotropic monopole emission. The largest effect is on the dipole, that is, on the anisotropy at the ℓ = 1 multipole, which is mainly attributed to the solar system barycentre motion. The study of the dipole anisotropy spectrum is an alternative way to the achievement of absolute measurements for extracting information about the background monopole spectrum. This approach was originally proposed by Danese & De Zotti (1981) within the framework of cosmic microwave background (CMB) spectral distortions that possibly occurred in the cosmic plasma at different epochs. This method has been exploited by Balashev et al. (2015) in the context of future CMB anisotropy missions and, in particular, by De Zotti et al. (2016) with the aim of improving the characterization of the cosmic infrared background (CIB) spectrum. Numerical simulations have been performed to assess the impact of instrumental performance, potential residuals from imperfect foreground subtraction, and relative calibration uncertainties in the reconstruction of the types of signals described above (Burigana et al. 2018). This differential approach has been investigated to be applied to the analysis of the redshifted 21 cm line (Slosar 2017) and of its diurnal pattern in drift-scan observations (Deshpande 2018). Recent predictions of the cosmic dipole from four types of imprints that are expected from (or associated with) cosmological reionization – the diffuse free-free (FF) emission, the Comptonization distortion, the redshifted 21 cm line, and the radio extragalactic background, along with combinations of these types – have been presented by Trombetti & Burigana (2019).

In this work, we carry out an analysis of the effect of the observer peculiar motion on the frequency spectra of the monopole and of the anisotropy patterns at higher multipoles for the monopole component of various types of background radiation, ranging from the radio to the far-infrared (far-IR). For a blackbody spectrum, the amplitude of this effect decreases as β^ℓ at increasing ℓ, where β ≪ 1 is the module of the dimensionless peculiar velocity of the observer, which is defined by the vector β = v/c, c being the speed of light.

In general, the vector β is the sum of an almost constant component, β_C, that is due to the motion of the Solar System barycentre with respect to the cosmic background and of a time varying component, β_V = β_V(t), which is due to the motion of the observer relative to the Solar System barycentre reference frame. For ground-based or sub-orbital experiments, β_V is given by the combination of the motions of the Earth around the Solar System barycentre (β_ES ≃ 10⁻⁴) and of the experimental equipment on the Earth’s surface (β_EE ≃ 1.7 × 10⁻⁶ for an experiment located at the Earth equator). For past, planned, and proposed CMB space missions, β_V is given by the combination of the motion of the Earth or of the second Lagrangian point of the Sun-Earth system (L2) around the Solar System barycentre (β_L2 ≃ 1.01 β_ES) and of the motion of spacecraft around the Earth or around L2, according to the adopted trajectory. For example, the typical spacecraft velocity around L2 is β_L ≃ 1.7 × 10⁻⁶ for a Lissajous ‘orbit’ with a ‘radius’ of ≃3 × 10⁵ km described over a timescale of around six months. We note that the Solar System motion around the Galactic centre and the presence of local Universe gravitational fields imply time variations of β, but they can be neglected within the typical duration, τ_S, of a survey. For example, in the approximation of a uniform circular motion of the Solar System around the Galactic centre, with a rotation period P = 2π/ω ≃ 2.4 × 10⁸ yr and a velocity v_SG ≃ 220 km s⁻¹, that is, β_SG ≃ 7.3 × 10⁻⁴, the relative variation of the velocity at a timescale of τ_S ≃ 5 yr is Δv_SG/v_SG = Δβ_SG/β_SG ≃ ω² ⋅ (v_SG/ω) ⋅ τ_S/v_SG = ω τ_S ≃ 1.3 × 10⁻⁷.

For numerical estimates, we will typically assume that the CMB dipole is due to velocity effects only and we neglect the modulation of β introduced by the contribution of β_V. After the correction for the spacecraft motion around the Solar System barycentre, the nominal CMB dipole amplitude according to the Planck 2015 results (Planck Collaboration I 2016; Planck Collaboration V 2016; Planck Collaboration VIII 2016) is A_dip = (3.3645 ± 0.002) mK. When using the Low Frequency Instrument alone (Planck Collaboration II 2020), the most recent analysis of the Planck 2018 results gives an almost identical value of A_dip as in the 2015 release. On the other hand, when including, again, the High Frequency Instrument, the analysis of the Planck 2018 results gives A_dip = (3.36208 ± 0.00099) mK (Planck Collaboration III 2020; Planck Collaboration I 2020), which is a slightly lower value. These results are clearly compatible within the errors. Based on a joint analysis (Fixsen 2009) of the data from the Far Infrared Absolute Spectrophotometer (FIRAS) on board the Cosmic Background Explorer (COBE) and from the Wilkinson Microwave Anisotropy Probe (WMAP), we adopt T₀ = (2.72548 ± 0.00057) K for the current CMB effective temperature in the blackbody spectrum approximation, where $a{T}_0^4$ gives the current CMB energy density with a = 8πI₃k⁴/(hc)³, I₃ = π⁴/15, k and h the Boltzmann and Planck constants. We use the velocity v_C = (369.82 ± 0.11) km s⁻¹ given in Table 3 of Planck Collaboration I (2020), that is, β = β_C = v_C/c ≃ 1.2336 × 10⁻³ ≃ A_dip/T₀, to characterise the velocity of the Solar System barycentre with respect to the cosmic background.

The main contribution to the modulation of β coming from β_V(t), which produces an effect amounting to ≃8.1% of the global signal, is derived from the component of the observer motion due to the revolution of the Earth or of L2 around the Solar System barycentre (see above estimates).

We study the effect of peculiar motion in terms of spherical harmonic expansion up to a certain value of ℓ_max. In this way, we introduce a relative error in the prediction of the effect at a given ℓ that strongly decreases with ℓ_max. Neglecting the contributions from higher orders, the dipole anisotropy spectrum was estimated as the difference between the signal measured in the direction of motion and in its perpendicular direction (Danese & De Zotti 1981), namely, in terms of a very simple linear combination of the signals in two specific directions. In this work, we show how this concept can be generalized to derive the frequency behaviour of the anisotropy pattern up to higher multipoles. We provide both a recipe and explicit solutions that can be directly used for accurate and swift theoretical predictions of the individual multipole patterns and of the global pattern, allowing us to bypass the need for more computationally demanding approaches that are based on delicate numerical integrations or on map generation and inversion.

In Sect. 2, we introduce the adopted formalism and the set of equations we aimed to solve. We provide the explicit equations and solutions up to ℓ_max = 6 in Sect. 3 and ℓ_max = 4 in Sect. 5 in order to point out some general properties of the solutions. In Sects. 4 and 6, we work out these solutions for the particular case of a blackbody spectrum to make clear their simple link with the contributions coming from the various orders of β. In some cases, we show their equivalence with the corresponding explicit exact solution at the order of β corresponding to ℓ_max. In Sects. 7 and 8, we discuss the solutions up to low ℓ_max, that is, ℓ_max = 1 and 2, derived using just two or three colatitudes. Some remarks linking the general properties of the found solutions at various ℓ_max with the monopole spectrum integration and differentiation are given in Sect. 9. The main applications and results of the proposed method are presented in Sect. 10 for eight specific types of background: we give concise presentations of the monopole spectrum models adopted in this study, we describe the main features of the found solutions and, for two very different cases, we compare them with the results based on a numerical integration. Some applications to combinations of signals are discussed in Sect. 11. In Sect. 12, we briefly present a set of results related to all-sky maps and angular power spectra, also for the purpose of comparison with previous analyses based on map generation and inversion. In Sect. 13, we discuss how the developed method can be directly generalized from the case of maps obtained with a constant observer velocity to the case of maps derived from the average of data taken with a set of different observer velocities, as, for example, in the case of β modulated by β_V(t). In Sect. 14, we focus on the global pattern at microwave frequencies, discussing the superposition of the CMB intrinsic anisotropies and of the effects induced by the observer peculiar motion. The possibility of constraining the intrinsic dipole embedded in the kinematic dipole in the presence of CMB spectral distortions is then discussed in Sect. 15. Finally, in Sect. 16, we draw our main conclusions. Some technical aspects are provided in the three sections of the appendix.

2. Theoretical framework and formalism

The peculiar velocity effect on the frequency spectrum can be evaluated on the whole sky using the complete description of the Compton-Getting effect (Forman 1970). This is based on the Lorentz invariance of the photon distribution function. In this work, we are interested in the effects induced on the monopole (or global) signal that, by definition, is isotropic in an ideal reference frame at rest with respect to the CMB or, more generally, to the cosmic background under consideration. In principle, the CMB and the other cosmic backgrounds provide information on processes that possibly occurred at different epochs or that are differently weighted for a range of redshift shells. Thus, the above ideal reference frame should correctly refer to the corresponding cosmic phase.

At a given ν, the photon distribution function, η^BB/dist, for the considered type of spectrum needs to be computed with the frequency multiplied by the product $(1-\widehat{n}·\mathit{\beta })/{(1-{\beta }^2)}^{1/2}$. The notation ‘BB/dist’ stands for a blackbody spectrum or for any type of non-blackbody signal (or for combinations of signals). This accounts for all the possible sky directions, which are defined by the unit vector $\widehat{n}$, relative to the peculiar velocity of the observer, which is defined by the vector β in the reference frame at rest with respect to the considered cosmic background. This includes all the orders in β and the link with the geometrical properties induced at each multipole. We study the effect in terms of equivalent thermodynamic temperature, T_th(ν), defined as the temperature of the blackbody having the same η(ν) at the frequency ν,

$$\begin{array}{c}\hfill T_{\mathrm{th}}(\nu )=\frac{h\nu }{k}\frac{1}{ln(1+1/\eta (\nu ))}·\end{array}$$(1)

The observed signal map is then given by (Burigana et al. 2018)

$$\begin{array}{c}\hfill T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\nu ,\widehat{n},\mathit{\beta })=\frac{x{T}_0}{\ln (1+1/{(\eta (\nu ,\widehat{n},\mathit{\beta }))}^{\mathrm{BB}/\mathrm{dist}})},\end{array}$$(2)

where $\eta (\nu ,\widehat{n},\mathit{\beta })=\eta ({\nu }^{\prime })$ with ${\nu }^{\prime }=\nu (1-\widehat{n}·\mathit{\beta })/{(1-{\beta }^2)}^{1/2}$, x = hν/(kT_r) and T_r = T₀(1 + z) are the redshift invariant dimensionless frequency and the redshift-dependent effective temperature of the CMB.

The unit vector $\widehat{n}$ is associated to the polar coordinates θ (colatitude) and ϕ (longitude). The function $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\nu ,\widehat{n},\mathit{\beta })=T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\nu ,\theta ,\varphi ,\mathit{\beta })$ can be expanded in spherical harmonics. We adopt a reference system with the z axis parallel to the observer velocity and we can then simply replace β with β in the above dependencies. Thus,

$$\begin{array}{c}\hfill T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\nu ,\theta ,\varphi ,\beta )={\displaystyle \sum _{\ell =0}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{a}_{\ell ,m}(\nu ,\beta )Y_{\ell ,m}(\theta ,\varphi ),}\end{array}$$(3)

where Y_ℓ, m(θ, ϕ) are the spherical harmonics related to the associated Legendre polynomials, $P_{\mathcal{l}}^m(\text{cos}\theta )$, and the coefficients a_ℓ, m(ν, β) contain information on the background spectrum and the observer velocity.

In the adopted reference system, the isotropy of the background monopole, or, equivalently, of η, implies that $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ depends on θ but not on ϕ. Thus, in Eq. (3), we can take only the terms with m = 0 and in this case, $Y_{\ell ,m}(\theta ,\varphi )={\stackrel{\sim }{P}}_{\ell }^m(\cos \theta )$, where ${\stackrel{\sim }{P}}_{\ell }^m(\cos \theta )$ are the renormalized associated Legendre polynomials:

$$\begin{array}{c}\hfill {\stackrel{\sim }{P}}_{\ell }^m(\cos \theta )=\sqrt{\frac{2\ell +1}{4\pi }\frac{(\ell -m)!}{(\ell +m)!}}{P}_{\ell }^m(\cos \theta ).\end{array}$$(4)

In general, for a real function, the coefficients of the spherical harmonics expansion with m >  0 are related to the coefficients with m <  0 by the relation $a_{\mathcal{l},m}^{*}={(-1)}^m{a}_{\mathcal{l},-m}$, where the index denotes the complex conjugation. We note that for this problem and with the adopted reference system with the z axis parallel (or antiparallel) to the observer velocity, we are interested only in the non-vanishing coefficients with m = 0, but, in general, we can also see that the coefficients a_ℓ, m with m ≠ 0 do not vanish. The publicly available tools allow us to efficiently compute the a_ℓ, m passing from a reference system to another (see Górski et al. 2005).

Formally, the coefficients a_ℓ, m(ν, β) can be computed through an inversion of Eq. (3):

$$\begin{array}{cc}\hfill a_{\ell ,m}(\nu ,\beta )& ={\displaystyle {\int }_{\theta =0}^{\pi }{\int }_{\varphi =0}^{2\pi }{T}_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\nu ,\theta ,\varphi ,\beta ){\mathrm{e}}^{-\mathrm{i}m\varphi }{\stackrel{\sim }{P}}_{\ell }^m(\cos \theta )\sin \theta \mathrm{d}\theta \mathrm{d}\varphi }\hfill \\ \hfill & ={\displaystyle {\int }_{\theta =0}^{\pi }{\int }_{\varphi =0}^{2\pi }\left[T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\nu ,\theta ,\varphi ,\beta )-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist},\mathrm{rest}}(\nu )\right]}\hfill \\ \hfill & \times {\mathrm{e}}^{-\mathrm{i}m\varphi }{\stackrel{\sim }{P}}_{\ell }^m(\cos \theta )\sin \theta \mathrm{d}\theta \mathrm{d}\varphi +a_{\ell ,m}^{\mathrm{rest}}(\nu ),\hfill \end{array}$$(5)

where $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ is evaluated through Eq. (2) and m = 0. In the last part of Eq. (5), $a_{\ell ,m}^{\mathrm{rest}}(\nu )=\sqrt{4\pi }{T}_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist},\mathrm{rest}}(\nu )$ and $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist},\mathrm{rest}}(\nu )$ are evaluated in the background rest frame, that is, these are the intrinsic spherical harmonics expansion coefficients and the intrinsic (isotropic) background monopole expressed in equivalent thermodynamic temperature. For this problem, $a_{\ell ,m}^{\mathrm{rest}}(\nu )$ does not vanish only for ℓ = 0. This form in Eq. (5) is useful in the numerical computation (see also Sect. 10.1) because the integrand function becomes the difference between the equivalent thermodynamic temperatures in the reference frames in motion and at rest with respect to the background. For a general background spectrum, this approach requires a delicate and computationally demanding integration over θ. For a small β, it could be difficult to achieve the extreme precision needed to characterize the fine and small details of spectral features. We can instead consider Eq. (3) with m = 0 for a set of N directions, namely, of colatitudes θ_i with i = 0, N − 1, to construct a linear system of N equations in the N unknowns a_ℓ, 0(ν, β), with ℓ = 0, N − 1, that can be solved given the corresponding N values of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\nu ,{\theta }_{\mathrm{i}},\varphi ,\beta )$, provided that the determinant of the system coefficient matrix does not vanish. The solutions for the unknowns a_ℓ, 0(ν, β) can be then written as linear combinations of N signals, $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\nu ,{\theta }_{\mathrm{i}},\varphi ,\beta )$, that are evaluated for a given background monopole at N colatitudes.

With this simple scheme, we can fully characterize the observed signal map $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\nu ,\widehat{n},\mathit{\beta })$ up to the desired multipole component ℓ_max = N − 1. Let us assume, as a rule of thumb, that the amplitude of this effect decreases at increasing multipole as β^{ℓ ⋅ p}, with p ≈ 1. The value appropriate to the case of a blackbody spectrum is p = 1, as mentioned in Sect. 1, while, in general, the effective scaling with ℓ is frequency-dependent and related to the monopole spectrum shape, as discussed in next sections. Considering a spherical harmonic expansion up to ℓ_max leads to neglect the contributions from ℓ >  ℓ_max. Thus, given the above scaling, the relative error in the computation of the effect at a given ℓ ≤ ℓ_max is, at most, on the order of β^{(ℓ_max − ℓ + j) ⋅ p}. For a generic choice of the N colatitudes, we simply have j = 1. Since β is on the order of 10⁻³, adopting ℓ_max = 6, we expect to achieve an extremely high numerical accuracy that is sufficient for any application even in the very distant future; whereas, setting ℓ_max = 4 can be adequate for predicting the corresponding multipole patterns as part of the analysis in forthcoming and planned (or proposed) surveys since no relevant error is introduced by neglecting the contributions at higher multipoles. In general, an accuracy up to any desired order can be then achieved with this approach by just computing $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ in only a relatively small number of sky directions, N = ℓ_max + 1.

We note that $Y_{0,0}=\sqrt{1/(4\pi )}$ and that for ℓ = 0, and for even ℓ, the associated Legendre polynomials $P_{\mathcal{l}}^0(\text{cos}\theta )$ are symmetric with respect to θ = π/2, whereas for odd ℓ, they vanish at θ = π/2 and are antisymmetric with respect to θ = π/2. This suggests that the linear system of N equations using θ = π/2 and pairs of colatitudes symmetric with respect to θ = π/2 are expected satisfy the following properties: (i) for θ = π/2, all the coefficients multiplying a_ℓ, 0 are null for odd ℓ; (ii) for each pair of colatitudes – if ℓ is even, the coefficient multiplying a_ℓ, 0 is the same, whereas, it is the opposite result for odd ℓ.

As made evident in the next sections, these properties can be used to significantly simplify the explicit solution of the system because they allow us to combine the equations into two separate subsystems: (1) (N − 1)/2 + 1 equations for a_ℓ, 0, with ℓ = 0 and even ℓ and (2) (N − 1)/2 equations for a_ℓ, 0 with odd ℓ. For even ℓ_max, a choice of odd N = ℓ_max + 1 colatitudes θ_i that satisfies the above symmetry implies j = 2 (instead of 1) for even ℓ in the scaling, β^{(ℓ_max − ℓ + j) ⋅ p}, of the relative error of the method (see also the discussion at the end of Appendix A). For odd ℓ_max, the system can be built with N = ℓ_max + 1 colatitudes θ_i as above, but avoiding the inclusion of π/2. The system can be split into two separate subsystems of N/2 equations: (1) for ℓ = 0 and even ℓ and (2) for odd ℓ; and in this case, j = 2 for odd ℓ (see also Sect. 7 and the discussion at the end of Sect. 8). This property allows us to achieve a significant improvement in accuracy with respect to a generic choice of the N colatitudes θ_i.

3. Explicit solutions up to ℓ_max= 6

Explicitly expanding $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ in spherical harmonics up to ℓ_max = 6, we get

$$\begin{array}{cc}\hfill T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}=& a_{0,0}\sqrt{\frac{1}{4\pi }}\hfill \\ \hfill & +a_{1,0}\sqrt{\frac{3}{4\pi }}\cos \theta \hfill \\ \hfill & +a_{2,0}\sqrt{\frac{5}{4\pi }}\left(\frac{3}{2}{\cos }^2\theta -\frac{1}{2}\right)\hfill \\ \hfill & +a_{3,0}\sqrt{\frac{7}{4\pi }}\left(\frac{5}{2}{\cos }^3\theta -\frac{3}{2}\cos \theta \right)\hfill \\ \hfill & +a_{4,0}\sqrt{\frac{9}{4\pi }}\left(\frac{35}{8}{\cos }^4\theta -\frac{15}{4}{\cos }^2\theta +\frac{3}{8}\right)\hfill \\ \hfill & +a_{5,0}\sqrt{\frac{11}{4\pi }}\left(\frac{63}{8}{\cos }^5\theta -\frac{35}{4}{\cos }^3\theta +\frac{15}{8}\cos \theta \right)\hfill \\ \hfill & +a_{6,0}\sqrt{\frac{13}{4\pi }}\left(\frac{231}{16}{\cos }^6\theta -\frac{315}{16}{\cos }^4\theta +\frac{105}{16}{\cos }^2\theta -\frac{5}{16}\right),\hfill \end{array}$$(6)

where we omit, for simplicity, the dependence of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ on ν, θ and β and the dependencies of a_ℓ, 0 on ν and β.

To write the linear system of seven equations, we are able to choose among infinite possibilities, and the explicit form of the system (but not the solution up to the adopted maximum multipole) depends on the specific choice. Among the possible choices satisfying the symmetry properties described above, we selected a set of colatitudes θ_i such that the values of cos θ_i are rational numbers or just involve $\sqrt{2}$ in order to simplify the algebra: θ_i = 0, π/4, π/3, π/2, (2/3)π, (3/4)π, and π.

After a series of calculations, we derived the corresponding linear system. We obtain

$$\begin{array}{cc}& T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)=\sqrt{\frac{1}{4\pi }}{a}_{0,0}+\sqrt{\frac{3}{4\pi }}{a}_{1,0}\hfill \\ \hfill & +\sqrt{\frac{5}{4\pi }}{a}_{2,0}+\sqrt{\frac{7}{4\pi }}{a}_{3,0}+\sqrt{\frac{9}{4\pi }}{a}_{4,0}\hfill \\ \hfill & +\sqrt{\frac{11}{4\pi }}{a}_{5,0}+\sqrt{\frac{13}{4\pi }}{a}_{6,0},\hfill \end{array}$$(7)

$$\begin{array}{cc}& T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)=\sqrt{\frac{1}{4\pi }}{a}_{0,0}+\frac{\sqrt{2}}{2}\sqrt{\frac{3}{4\pi }}{a}_{1,0}\hfill \\ \hfill & +\frac{1}{4}\sqrt{\frac{5}{4\pi }}{a}_{2,0}-\frac{\sqrt{2}}{8}\sqrt{\frac{7}{4\pi }}{a}_{3,0}-\frac{13}{32}\sqrt{\frac{9}{4\pi }}{a}_{4,0}\hfill \\ \hfill & -\frac{17\sqrt{2}}{64}\sqrt{\frac{11}{4\pi }}{a}_{5,0}-\frac{19}{128}\sqrt{\frac{13}{4\pi }}{a}_{6,0},\hfill \end{array}$$(8)

$$\begin{array}{cc}& T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /3)=\sqrt{\frac{1}{4\pi }}{a}_{0,0}+\frac{1}{2}\sqrt{\frac{3}{4\pi }}{a}_{1,0}\hfill \\ \hfill & -\frac{1}{8}\sqrt{\frac{5}{4\pi }}{a}_{2,0}-\frac{7}{16}\sqrt{\frac{7}{4\pi }}{a}_{3,0}-\frac{37}{128}\sqrt{\frac{9}{4\pi }}{a}_{4,0}\hfill \\ \hfill & +\frac{23}{256}\sqrt{\frac{11}{4\pi }}{a}_{5,0}+\frac{331}{1024}\sqrt{\frac{13}{4\pi }}{a}_{6,0},\hfill \end{array}$$(9)

$$\begin{array}{cc}& T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)=\sqrt{\frac{1}{4\pi }}{a}_{0,0}+0·a_{1,0}\hfill \\ \hfill & -\frac{1}{2}\sqrt{\frac{5}{4\pi }}{a}_{2,0}+0·a_{3,0}+\frac{3}{8}\sqrt{\frac{9}{4\pi }}{a}_{4,0},\hfill \\ \hfill & +0·a_{5,0}-\frac{5}{16}\sqrt{\frac{13}{4\pi }}{a}_{6,0}\hfill \end{array}$$(10)

$$\begin{array}{cc}& T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(2/3)\pi )=\sqrt{\frac{1}{4\pi }}{a}_{0,0}-\frac{1}{2}\sqrt{\frac{3}{4\pi }}{a}_{1,0}\hfill \\ \hfill & -\frac{1}{8}\sqrt{\frac{5}{4\pi }}{a}_{2,0}+\frac{7}{16}\sqrt{\frac{7}{4\pi }}{a}_{3,0}-\frac{37}{128}\sqrt{\frac{9}{4\pi }}{a}_{4,0}\hfill \\ \hfill & -\frac{23}{256}\sqrt{\frac{11}{4\pi }}{a}_{5,0}+\frac{331}{1024}\sqrt{\frac{13}{4\pi }}{a}_{6,0}\hfill \end{array}$$(11)

$$\begin{array}{cc}& T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )=\sqrt{\frac{1}{4\pi }}{a}_{0,0}-\frac{\sqrt{2}}{2}\sqrt{\frac{3}{4\pi }}{a}_{1,0}\hfill \\ \hfill & +\frac{1}{4}\sqrt{\frac{5}{4\pi }}{a}_{2,0}+\frac{\sqrt{2}}{8}\sqrt{\frac{7}{4\pi }}{a}_{3,0}-\frac{13}{32}\sqrt{\frac{9}{4\pi }}{a}_{4,0}\hfill \\ \hfill & +\frac{17\sqrt{2}}{64}\sqrt{\frac{11}{4\pi }}{a}_{5,0}-\frac{19}{128}\sqrt{\frac{13}{4\pi }}{a}_{6,0},\hfill \end{array}$$(12)

$$\begin{array}{cc}& T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )=\sqrt{\frac{1}{4\pi }}{a}_{0,0}-\sqrt{\frac{3}{4\pi }}{a}_{1,0}\hfill \\ \hfill & +\sqrt{\frac{5}{4\pi }}{a}_{2,0}-\sqrt{\frac{7}{4\pi }}{a}_{3,0}+\sqrt{\frac{9}{4\pi }}{a}_{4,0}\hfill \\ \hfill & -\sqrt{\frac{11}{4\pi }}{a}_{5,0}+\sqrt{\frac{13}{4\pi }}{a}_{6,0}.\hfill \end{array}$$(13)

Equations (7)–(13) constitute the linear system that is to be solved; and the determinant of the coefficients of the system matrix, ≃ − 0.303, does not vanish. We can solve the system using the methods of elimination and substitution.

As anticipated, we can combine the above equations to split the system into two subsystems. By adding the left and right sides of Eqs. (7) and (13), of Eqs. (8) and (12), and of Eqs. (9) and (11), we get three equations which, complemented using Eq. (10), form a linear system involving only the four unknowns a_ℓ, 0 with ℓ = 0 and even ℓ. We can solve it by substitution. Equation (10) allows us to express $2\sqrt{1/(4\pi )}{a}_{0,0}$ as a combination of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)$, a_2, 0, a_4, 0, and a_6, 0 to be put in the other three equations. From the first one, we then express a_2, 0 as a combination of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)$, $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)$, $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )$, a_4, 0, and a_6, 0 to be put in the two other remaining equations. We then represent a_4, 0 as a combination of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)$, $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)$, $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)$, $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )$, $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )$, and a_6, 0, which allows us to derive first the solution for a_6, 0:

$$\begin{array}{cc}\hfill a_{6,0}=& \frac{64}{693}\sqrt{\frac{4\pi }{13}}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & -6\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & +8\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /3)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(2/3)\pi )\right)\hfill \\ \hfill & -6T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)].\hfill \end{array}$$(14)

As suggested in the introduction, a_6, 0 is written in terms of a linear combination of the set of values of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ computed for the seven adopted colatitudes. With substitution, we subsequently derive the solution for a_4, 0, a_2, 0 and, finally, for a_0, 0:

$$\begin{array}{cc}\hfill a_{4,0}& =\frac{8}{385}\sqrt{\frac{4\pi }{9}}[9\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & -10\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & -16\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /3)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(2/3)\pi )\right)\hfill \\ \hfill & +34T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)],\hfill \end{array}$$(15)

$$\begin{array}{cc}\hfill a_{2,0}& =\frac{1}{693}\sqrt{\frac{4\pi }{5}}[121\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & +396\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & -352\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /3)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(2/3)\pi )\right)\hfill \\ \hfill & -330T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)],\hfill \end{array}$$(16)

$$\begin{array}{cc}\hfill a_{0,0}& =\frac{\sqrt{4\pi }}{630}[29\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & +120\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & +64\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /3)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(2/3)\pi )\right)\hfill \\ \hfill & +204T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)].\hfill \end{array}$$(17)

Subtracting left and right sides of Eqs. (7) and (13), of Eqs. (8) and (12), and of Eqs. (9) and (11) we get three equations that form a linear system involving only the three unknowns a_ℓ, 0 with odd ℓ. From the difference between the first of these equations and the second equation multiplied by $\sqrt{2}$ and the difference between the first equation and the third equation multiplied by 2, we can write a system for a_3, 0 and a_5, 0. We first derive a_5, 0

$$\begin{array}{cc}\hfill a_{5,0}=& \frac{32}{189}\sqrt{\frac{4\pi }{11}}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & -3\sqrt{2}\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & +4\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /3)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(2/3)\pi )\right)]\hfill \end{array}$$(18)

and then, by substitution, a_3, 0 and a_1, 0

$$\begin{array}{cc}\hfill a_{3,0}& =\frac{2}{135}\sqrt{\frac{4\pi }{7}}[13\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & +15\sqrt{2}\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & -56\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /3)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(2/3)\pi )\right)],\hfill \end{array}$$(19)

$$\begin{array}{cc}\hfill a_{1,0}& =\frac{1}{210}\sqrt{\frac{4\pi }{3}}[29\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & +60\sqrt{2}\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & +32\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /3)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(2/3)\pi )\right)].\hfill \end{array}$$(20)

We note that the structure of the solutions for a_ℓ, 0 for ℓ = 0, and even ℓ, involving the sums of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ at pairs of colatitudes symmetric with respect to π/2 and $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ at π/2, as well as the structure of the solutions for odd ℓ, involving the differences of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ at pairs of colatitudes symmetric with respect to π/2, reflect the symmetry and antisymmetry properties that are discussed at the end of Sect. 2, together with the corresponding implications for the system solution accuracy.

The solutions expressed in Eqs. (14)–(20) can be compared with each of the Eqs. (7)–(13) for a given colatitude θ_i: as expected, the sum of products of the various coefficients that multiply a_ℓ, 0, for ℓ = 0, 6, in the equation for $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}({\theta }_{\mathrm{i}})$ with the coefficients in Eqs. (14)–(20) that multiply $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}({\theta }_{\mathrm{i}})$ gives exactly one. Remarkably, except for θ_i = π/2, where only ℓ = 0 and the even multipoles contribute to $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$, for all the other colatitudes θ_i the above sum is equally contributed for one half by ℓ = 0 and by the even multipoles and for one half by the odd multipoles. This is another property related to the symmetry with respect to π/2 of the set of colatitudes adopted.

4. Solutions for a blackbody up to ℓ_max= 6

Let us consider the specific case of the CMB, assumed to ideally exhibit a blackbody monopole spectrum with an effective temperature, T₀, in the CMB rest frame. In this case, the photon distribution function is

$$\begin{array}{c}\hfill {\eta }^{\mathrm{BB}}(\nu )={\eta }^{\mathrm{BB}}(x)=\frac{1}{e^x-1}·\end{array}$$(21)

Equation (2) then gives the well-known expression

$$\begin{array}{c}\hfill T_{\mathrm{th}}^{\mathrm{BB}}(\nu ,\theta ,\varphi ,\beta )=\frac{x{T}_0}{x^{\prime }}=\frac{T_0{(1-{\beta }^2)}^{1/2}}{1-\beta \cos \theta },\end{array}$$(22)

with x′ = hν′/(kT_r), highlighting that $T_{\mathrm{th}}^{\mathrm{BB}}(\nu ,\theta ,\varphi ,\beta )$ does not depend on ϕ nor on ν.

The observed CMB effective temperature averaged over the full sky, T_0, obs, is given by

$$\begin{array}{c}\hfill T_{0,\mathrm{obs}}=\frac{1}{4\pi }{\displaystyle {\int }_{\theta =0}^{\pi }{\int }_{\varphi =0}^{2\pi }\frac{T_0{(1-{\beta }^2)}^{1/2}}{1-\beta \cos \theta }\sin \theta \mathrm{d}\theta \mathrm{d}\varphi .}\end{array}$$(23)

For an observer at rest with respect to the CMB, β = 0 and then the substitution of the integration variable, θ, with a new variable, w = cos θ, obviously gives T_0, obs = T₀ and implies that in the expansion represented by Eq. (3), as specified by Eq. (6), the only non-vanishing contribution to T_0, obs comes from a term associated to the multipole coefficient a_0, 0. Since $Y_{0,0}=\sqrt{1/(4\pi )}$, $a_{0,0}=T_0\sqrt{4\pi }$.

For an observer in motion with respect to the CMB, β ≠ 0 and T_0, obs can be calculated by simply substituting the integration variable, θ, with w = 1 − βcos θ. We get

$$\begin{array}{c}\hfill T_{0,\mathrm{obs}}=\frac{1}{2}{(1-{\beta }^2)}^{1/2}{T}_0\frac{1}{\beta }\ln \frac{1+\beta }{1-\beta },\end{array}$$(24)

as already reported in Lucca et al. (2020) (see also Chluba 2011; Dai & Chluba 2014). By replacing ln[(1 + β)/(1 − β)] with its expansion in Taylor’s series up to β⁷, that is, with 2[β + (1/3)β³ + (1/5)β⁵ + (1/7)β⁷], we find:

$$\begin{array}{c}\hfill T_{0,\mathrm{obs}}=\frac{1}{2}{(1-{\beta }^2)}^{1/2}{T}_0·2·[1+(1/3){\beta }^2+(1/5){\beta }^4+(1/7){\beta }^6].\end{array}$$(25)

We now specify the coefficients, a_ℓ, 0, given by Eqs. (14)–(20) to the blackbody case using Eq. (22) to compute $T_{\mathrm{th}}^{\mathrm{BB}}$ at the seven considered colatitudes. After solving the algebra, we get

$$\begin{array}{cc}& a_{6,0}^{\mathrm{BB}}=\frac{128}{231}\sqrt{\frac{4\pi }{13}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{{\beta }^6}{(1-{\beta }^2)(2-{\beta }^2)(4-{\beta }^2)},\hfill \end{array}$$(26)

$$\begin{array}{cc}& a_{4,0}^{\mathrm{BB}}=\frac{16}{385}\sqrt{\frac{4\pi }{9}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{44{\beta }^4-17{\beta }^6}{(1-{\beta }^2)(2-{\beta }^2)(4-{\beta }^2)},\hfill \end{array}$$(27)

$$\begin{array}{cc}& a_{2,0}^{\mathrm{BB}}=\frac{2}{21}\sqrt{\frac{4\pi }{5}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{56{\beta }^2-50{\beta }^4+5{\beta }^6}{(1-{\beta }^2)(2-{\beta }^2)(4-{\beta }^2)},\hfill \end{array}$$(28)

$$\begin{array}{cc}& a_{0,0}^{\mathrm{BB}}=\frac{1}{2}\sqrt{4\pi }{(1-{\beta }^2)}^{1/2}{T}_0\hfill \\ \hfill & ·\frac{16-(68/3){\beta }^2+(118/15){\beta }^4-(204/315){\beta }^6}{(1-{\beta }^2)(2-{\beta }^2)(4-{\beta }^2)},\hfill \end{array}$$(29)

and

$$\begin{array}{cc}& a_{5,0}^{\mathrm{BB}}=\frac{64}{63}\sqrt{\frac{4\pi }{11}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{{\beta }^5}{(1-{\beta }^2)(2-{\beta }^2)(4-{\beta }^2)},\hfill \end{array}$$(30)

$$\begin{array}{cc}& a_{3,0}^{\mathrm{BB}}=\frac{4}{45}\sqrt{\frac{4\pi }{7}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{36{\beta }^3-23{\beta }^5}{(1-{\beta }^2)(2-{\beta }^2)(4-{\beta }^2)},\hfill \end{array}$$(31)

$$\begin{array}{cc}& a_{1,0}^{\mathrm{BB}}=\sqrt{\frac{4\pi }{3}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{8\beta -(46/5){\beta }^3+(71/35){\beta }^5}{(1-{\beta }^2)(2-{\beta }^2)(4-{\beta }^2)}·\hfill \end{array}$$(32)

In the above expressions, the factor (1 − β²)^1/2T₀ clearly comes from Eq. (22) while the three factors in the denominator come from the choice of the pairs of colatitudes θ symmetric to π/2, that have been set to 0 and π, π/4 and (3/4)π, π/3, and (2/3)π. We observe also that, because of the adopted ℓ_max = 6 and the separation of the system into two subsystems, the solutions for $a_{6,0}^{\mathrm{BB}}$ and $a_{5,0}^{\mathrm{BB}}$ (see Eqs. (26) and (30)) do not show at numerator additional terms coming from higher multipoles, while they appear in the solutions for $a_{\ell ,0}^{\mathrm{BB}}$ for ℓ ≤ 4 (see Eqs. (27)–(29), (31) and (32)).

It is evident that the coefficients, $a_{\ell ,0}^{\mathrm{BB}}$, do not depend on θ. Thus, the substitution of the integration variable θ with w = cos θ, in the expansion represented by Eq. (3) again gives T_0, obs and implies that the only non-vanishing contribution comes from a term associated to the multipole coefficient $a_{0,0}^{\mathrm{BB}}$. An expansion in Taylor’s series up to β⁶, gives 1/(1 − β²) = 1 + β² + β⁴ + β⁶, 1/(2 − β²) = (1 + β²/2 + β⁴/4 + β⁶/8)/2, and 1/(4 − β²) = (1 + β²/4 + β⁴/16 + β⁶/64)/4. It is then simple to verify that, at the same order in β, Eq. (29) gives exactly the result expressed by Eq. (25), as is required, in principle.

Equation (5) allows us to analytically derive the a_ℓ, 0(ν, β) for any ℓ for relatively simple dependencies of $T_{\mathrm{th}}^{\mathrm{BB}}(\nu ,\theta ,\varphi ,\beta )$, as in the case of the blackbody spectrum, namely, for Eq. (22). The form of the integrand in θ involves only the function sin θ/(1 − βcos θ) multiplied by polynomials in cos θ and when substituting the integration variable θ with w = 1 − βcos θ, the integrand consists only of functions as 1/w and powers of w. We omit the tedious calculation at ℓ ≥ 2. Instead, for ℓ = 1, we get:

$$\begin{array}{cc}\hfill a_{1,0}^{\mathrm{BB}}& =2\pi {(1-{\beta }^2)}^{1/2}{T}_0{\displaystyle {\int }_{\theta =0}^{\pi }\frac{\sin \theta }{1-\beta \cos \theta }\sqrt{\frac{3}{4\pi }}\cos \theta \mathrm{d}\theta }\hfill \\ \hfill & =2\pi \sqrt{\frac{3}{4\pi }}{(1-{\beta }^2)}^{1/2}{T}_0\frac{1}{{\beta }^2}\left[\ln \frac{1+\beta }{1-\beta }-2\beta \right],\hfill \end{array}$$(33)

where, by replacing ln[(1 + β)/(1 − β)] with its expansion in Taylor’s series up to β⁷, gives:

$$\begin{array}{c}\hfill a_{1,0}^{\mathrm{BB}}=\sqrt{4\pi }\sqrt{3}{(1-{\beta }^2)}^{1/2}{T}_0[\beta /3+{\beta }^3/5+{\beta }^5/7].\end{array}$$(34)

Performing a Taylor’s series expansion up to β⁶ for 1/(1 − β²), 1/(2 − β²) and 1/(4 − β²), it is simple to verify that at the same order in β, Eq. (32) for $a_{1,0}^{\mathrm{BB}}$ gives precisely the result expressed by Eq. (34).

Finally, we remember that in Sect. 3 we discussed deriving both a_0, 0 and a_1, 0 by substitution in the last step of the calculation to solve the corresponding linear subsystem. Thus, the consistencies discussed above for T_0, obs and $a_{1,0}^{\mathrm{BB}}$ at the adopted order also represent a further verification of the derived algebraic solutions.

5. Explicit solutions up to ℓ_max= 4

For many applications, a computation up to ℓ_max = 4 suffice to get the relevant information. We then provide here simpler solutions based on N = 5 equations using the set of colatitudes 0, π/4, π/2, (3/4)π and π. This also allows us to explicitly focus on some of the mentioned properties of the proposed method.

We construct a system formed only by Eqs. (7), (8), (10), (12) and (13) ignoring the terms associated to a_5, 0 and a_6, 0, and we again solve it with the methods of elimination and substitution in a way similar to that described in Sect. 3. The solutions are

$$\begin{array}{cc}\hfill a_{4,0}& =\frac{8}{35}\sqrt{\frac{4\pi }{9}}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & -2\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & +2T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)],\hfill \end{array}$$(35)

$$\begin{array}{cc}\hfill a_{2,0}& =\frac{1}{21}\sqrt{\frac{4\pi }{5}}[5\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & +4\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & -18T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)],\hfill \end{array}$$(36)

$$\begin{array}{cc}\hfill a_{0,0}& =\frac{\sqrt{4\pi }}{30}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & +8\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)\hfill \\ \hfill & +12T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)].\hfill \end{array}$$(37)

Next, we have

$$\begin{array}{cc}\hfill a_{3,0}& =\frac{2}{5}\sqrt{\frac{4\pi }{7}}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & -\sqrt{2}\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)],\hfill \end{array}$$(38)

$$\begin{array}{cc}\hfill a_{1,0}& =\frac{1}{10}\sqrt{\frac{4\pi }{3}}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & +4\sqrt{2}\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /4)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =(3/4)\pi )\right)].\hfill \end{array}$$(39)

While the structure of the solutions expressed by Eqs. (35)–(39) is analogous to the structure of Eqs. (15)–(17), (19) and (20), the different algebraic coefficients reflect the different choice adopted for the set of colatitudes.

6. Solutions for a blackbody up to ℓ_max= 4

Specifying the coefficients, a_ℓ, 0, given by Eqs. (35)–(39) to the blackbody case (see Eq. (22)), we compute $T_{\mathrm{th}}^{\mathrm{BB}}$ at the five considered colatitudes and we obtain

$$\begin{array}{cc}& a_{4,0}^{\mathrm{BB}}=\frac{16}{35}\sqrt{\frac{4\pi }{9}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{{\beta }^4}{(1-{\beta }^2)(2-{\beta }^2)},\hfill \end{array}$$(40)

$$\begin{array}{cc}& a_{2,0}^{\mathrm{BB}}=\frac{2}{21}\sqrt{\frac{4\pi }{5}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{14{\beta }^2-9{\beta }^4}{(1-{\beta }^2)(2-{\beta }^2)},\hfill \end{array}$$(41)

$$\begin{array}{cc}& a_{0,0}^{\mathrm{BB}}=\frac{\sqrt{4\pi }}{30}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{60-70{\beta }^2+12{\beta }^4}{(1-{\beta }^2)(2-{\beta }^2)},\hfill \end{array}$$(42)

and

$$\begin{array}{cc}& a_{3,0}^{\mathrm{BB}}=\frac{4}{5}\sqrt{\frac{4\pi }{7}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{{\beta }^3}{(1-{\beta }^2)(2-{\beta }^2)},\hfill \end{array}$$(43)

$$\begin{array}{cc}& a_{1,0}^{\mathrm{BB}}=\sqrt{\frac{4\pi }{3}}{(1-{\beta }^2)}^{1/2}{T}_0·\frac{2\beta -(9/5){\beta }^3}{(1-{\beta }^2)(2-{\beta }^2)}·\hfill \end{array}$$(44)

In this case (ℓ_max = 4), the solutions for $a_{4,0}^{\mathrm{BB}}$ and $a_{3,0}^{\mathrm{BB}}$ do not show, in terms of the numerator, additional higher multipoles terms, as they appear at ℓ ≤ 2. The algebraic coefficients appearing in Eqs. (40)–(44) and in Eqs. (27)–(29), (31) and (32) are different, but these sets of equations give exactly the same solutions when the ratios of their polynomials in β are computed up to the order of β⁴. Analogously, Eq. (42) gives for T_0, obs the same result of Eqs. (25) and (44) is equivalent to Eq. (34) when they are computed up to the same order of power in β.

7. Explicit solutions up to ℓ_max= 2 and 1

It is helpful to write the solutions for low values of ℓ_max.

From the N = 3 equations at the colatitudes 0, π/2 and π (Eqs. (7), (10), and (13)), ignoring the terms associated to a_ℓ, 0 for ℓ >  2, we get

$$\begin{array}{cc}& a_{2,0}=\frac{1}{3}\sqrt{\frac{4\pi }{5}}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & -2T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)],\hfill \end{array}$$(45)

$$\begin{array}{cc}& a_{1,0}=\frac{1}{2}\sqrt{\frac{4\pi }{3}}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)],\hfill \end{array}$$(46)

$$\begin{array}{cc}& a_{0,0}=\frac{\sqrt{4\pi }}{6}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)\hfill \\ \hfill & +4T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)].\hfill \end{array}$$(47)

Using only N = 2 equations at the colatitudes 0 and π (Eqs. (7) and (13)), neglecting the terms at ℓ >  1, we have

$$\begin{array}{cc}& a_{1,0}=\frac{1}{2}\sqrt{\frac{4\pi }{3}}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)],\hfill \end{array}$$(48)

$$\begin{array}{cc}& a_{0,0}=\frac{\sqrt{4\pi }}{2}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)+T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi )\right)].\hfill \end{array}$$(49)

With N = 2 equations but at the colatitudes 0 and π/2 (Eqs. (7) and (10)), which is clearly not symmetric with respect to π/2, instead we get:

$$\begin{array}{cc}& a_{1,0}=\sqrt{\frac{4\pi }{3}}[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =0)-T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)\right)],\hfill \end{array}$$(50)

$$\begin{array}{cc}& a_{0,0}=\sqrt{4\pi }[\left(T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}(\theta =\pi /2)\right)].\hfill \end{array}$$(51)

8. On dipole estimations based on two colatitudes

As proposed by Danese & De Zotti (1981), a suitable and observationally intuitive approximation for the dipole spectrum can be expressed in terms of the difference of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ in the direction of motion and in its perpendicular direction.

Equations (7) and (10) allow us to express this difference in terms of a combination of the coefficients a_ℓ, 0(ν, β) up to ℓ_max = 6:

$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_{0,\pi /2}{T}_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}=& \sqrt{\frac{3}{4\pi }}{a}_{1,0}+\frac{3}{2}\sqrt{\frac{5}{4\pi }}{a}_{2,0}+\sqrt{\frac{7}{4\pi }}{a}_{3,0}\hfill \\ \hfill & +\frac{5}{8}\sqrt{\frac{9}{4\pi }}{a}_{4,0}+\sqrt{\frac{11}{4\pi }}{a}_{5,0}+\frac{21}{16}\sqrt{\frac{13}{4\pi }}{a}_{6,0}.\hfill \end{array}$$(52)

Neglecting the contributions from ℓ >  1, Eqs. (50) and (52) are equivalent. We can also express the semi-difference in $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ measured in the direction of motion and in its opposite direction using Eqs. (7) and (13)

$$\begin{array}{c}\hfill \frac{1}{2}{\mathrm{\Delta }}_{0,\pi }{T}_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}=\sqrt{\frac{3}{4\pi }}{a}_{1,0}+\sqrt{\frac{7}{4\pi }}{a}_{3,0}+\sqrt{\frac{11}{4\pi }}{a}_{5,0}\end{array}$$(53)

and in leaving out the terms at ℓ >  1, Eq. (53) is equivalent to Eq. (48).

The estimation of a_1, 0(ν, β) through the simple difference of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ in only two directions can be performed using the two colatitudes θ = 0 and θ = π to automatically suppress the contributions from ℓ = 2 (and from higher even ℓ), as discussed at the end of Sect. 2. The same holds for any other pair of colatitudes symmetric with respect to π/2 (as can be derived, for example, by combining Eqs. (8) and (12) or Eqs. (9) and (11)). The solutions presented in Sect. 3 can instead be used to correct for the contributions from the odd terms at ℓ = 3 and 5 (and from the terms at even ℓ, when using Eq. (52)).

9. Solutions for a_ℓ, 0 and spectrum integration or differentiation

Equation (5) shows that the solution for a given a_ℓ, 0 is an integral of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ over θ. Since β is very small, when θ spans in the interval [0, π], around π/2, the values of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ in the integrand come from frequency values in a small interval around ν/(1 − β²)^1/2 (see Eq. (2) and the relation between ν and ν′), making the integral sensitive to the local variation of $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$. Formally, the solutions expressed by Eqs. (14)–(20) and (35)–(39) can be regarded as definitions of sets of weights assigned to a small number of values of function $T_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ in a given set of colatitudes, or corresponding frequencies, to compute the integrals that give the coefficients a_ℓ, m in Eq. (5).

To a first-order approximation, the dipole spectrum induced by the observer peculiar velocity is directly proportional to the first logarithmic derivative of the photon occupation number, η(ν), with respect to the frequency, ν (Danese & De Zotti 1981). This concept can be generalized to higher multipoles. Let us consider the partial derivative of Eq. (5) with respect to the frequency, ν. According to Leibniz’s rule, when performing the differentiation under the integral sign, a further multiplicative factor involving the product, β cos θ, enters in the integral over θ, other than factors depending on the form of η(ν). As is evident from Eq. (6), a further power of cos θ appears passing from ℓ to ℓ + 1 in the associated Legendre polynomials and, consequently, in the integrand function of a_ℓ, 0 (see Eq. (5)). Thus, the subsequent a_ℓ, 0(ν) at increasing ℓ is tightly related to the subsequent derivatives of η(ν) with respect to ν, or in other words, their frequency behaviours are particular sensitive to the local (in frequency space) monopole spectrum variation up to increasing derivative order.

It is interesting to note certain properties of the coefficients (or weights) in Eqs. (14)–(20), (35)–(39), (45)–(47), (48) and (49) that are related to the separation of odd and even multipoles in the system solution. As already discussed, this separation appears when we adopt sets of colatitudes θ symmetrically located around π/2. The central weight, applied to θ = π/2, is zero for odd ℓ but not for even ℓ. For angles θ symmetric with respect to π/2, the weights are opposite for odd ℓ and equal for even ℓ. The sum of the weights vanishes, except for ℓ = 0: in this case, the sum is exactly unit, when divided by the ‘normalization’ factor, $\sqrt{4\pi }$ (see also Eq. (6)). These properties are identical to those satisfied by the weights for the centred approximations at a grid point for the generation of finite difference formulas on arbitrarily spaced grids for any order of derivative (Fornberg 1988, 1998). Furthermore, we note that the relative weights in Eqs. (45), (46), and (48) are equivalent to the relative weights for the centred approximations at a grid point for the second and first order of derivative, the relative weights in Eq. (49) are equivalent to the relative weights for the centred approximations at the halfway point for the zero order of derivative, while the relative weights in Eqs. (50) and (51) are equivalent to the relative weights for the one-sided approximations at a grid point for the first and zero order of derivative. The different level of approximation in the estimate of a_1, 0 via Eqs. (52) and (53), neglecting terms at ℓ >  1, is clearly related to the different accuracies of the one-sided and centred scheme for numerical differentiation. Finally, the relative weights for a_0, 0 in Eq. (47) are not equivalent to relative weights for the zero order of derivative of the schemes mentioned above. This is of increasing evidence in the weights of the solutions at ℓ_max >  2. Remarkably, they do not satisfy the sign alternation appearing in the weights of the centred approximations at a grid point of finite difference formulas moving from the central node to the more external nodes. Indeed, they store the relations between the a_ℓ, 0 at different ℓ and the temperatures at the adopted set of colatitudes that originates from the system solution at the corresponding ℓ_max (this is analogous to the ‘mixing’ of derivatives discussed above).

10. Monopole spectrum models and single signal results

The method described can be applied to any type of signal and to combinations of signals, provided that they are summed in terms of additive quantities, such as the photon distribution function, η, or the antenna temperature,

$$\begin{array}{c}\hfill T_{\mathrm{ant}}(\nu )=\frac{h\nu }{k}\eta (\nu ).\end{array}$$(54)

In this work, we consider eight different types of monopole spectrum that can be represented in terms of analytical or semi-analytical functions.

We first focus on four types of signals characterized by a CMB-distorted photon distribution function, η^dist(ν) that is different from the blackbody, η^BB(ν), at the present temperature T₀. We then consider four types of extragalactic background superimposed onto the CMB blackbody spectrum. We give only a concise description of the various models, referring to the literature for further information. On the other hand, we report the equations relevant for a clear connection with Sect. 14.

We first consider the signals that are more relevant (or essentially relevant) at low frequencies (radio domain) and then those that are relevant over a very wide frequency range (up to the far-IR) or more important at increasing frequency. We compare the results based on the proposed method (the solutions in Sect. 3) with the computation based on direct numerical integration (see Eq. (5) and the discussion in Sect. 10.1). For simplicity, we perform the comparison (see also Appendix A) only for two representative cases, which were chosen because they are very different with regard to the spectrum features.

The results are presented in terms of the following quantities:

– The difference, ΔT_th, between the equivalent thermodynamic temperature of the intrinsic monopole spectrum for the adopted model and the CMB present temperature T₀.

– The ratio, $R=(a_{0,0}(\nu ,\beta )/\sqrt{4\pi })/T_{\mathrm{th}}(\nu )$, between the equivalent thermodynamic temperature of observed (see Eq. (6)) and intrinsic monopole, expressed in terms of the difference ΔR = R − R^BB, where R^BB ≃ (1−2.5362 × 10⁻⁷) is the same ratio but for the case of the blackbody, η^BB(ν), at the present temperature T₀ (see Eqs. (23)–(25) and the discussion in Sect. 4).

– The coefficients a_ℓ, 0(ν, β) from ℓ = 1 to ℓ_max = 6 (see Eqs. (14)–(16) and (18)–(20)) expressed in terms of their difference with the blackbody case (see Eqs. (26)–(28) and (30)–(32)).

10.1. Possible non-equilibrium imprint at low frequencies

An important extragalactic background signal that is much larger than the CMB background predicted for a blackbody spectrum at an equilibrium temperature in agreement with FIRAS results is observed at radio frequencies, in particular, below a few GHz (see e.g., Dowell & Taylor 2018). A signal excess could be also present at ≃3.3 GHz, as claimed by Singal et al. (2011) on the basis of the second generation of the Absolute Radiometer for Cosmology, Astrophysics, and Diffuse Emission (ARCADE 2) data. Models based on contributions by faint astrophysical sources, on interactions between dark matter (DM) and baryons, or on their combinations have been invoked to explain this background, possibly together with the pronounced absorption profile of the 21 cm redshifted line signal (see also Sect. 10.4), which has also been claimed by Bowman et al. (2018); see, for example, Seiffert et al. (2011), Barkana (2018), Muñoz & Loeb (2018), Ewall-Wice et al. (2018) and Mirabel (2019) (see also Subrahmanyan & Cowsik 2013; Hills et al. 2018; Sharma 2018).

Baiesi et al. (2020) proposed an alternative explanation of the signal excess in the low frequency background, involving a mechanism of stochastic frequency diffusion in the perspective of non-equilibrium statistical mechanics. The model implies a modification of the standard Kompaneets equation (Kompaneets 1957), explicitly considered by the authors in the limit that includes only the scattering, and a relaxation of the Einstein detailed balance relation. The resulting abundance of low frequency photons can be described by a stationary solution of the photon distribution function in the form:

$$\begin{array}{c}\hfill \eta (\nu )=\frac{1}{\exp \left[{\int }^{\nu }\mathrm{d}{\nu }^{\prime }\gamma ({\nu }^{\prime })/D({\nu }^{\prime })\right]-1}=\frac{1}{e^{\psi (\nu )}-1},\end{array}$$(55)

where γ and D are frequency dependent friction and diffusion terms and the function ψ(ν) can be approximated by

$$\begin{array}{c}\hfill \psi (\nu )=\frac{h\nu }{k{T}^{\star }}\frac{{(\nu /{\nu }_0)}^{\alpha }}{1+{(\nu /{\nu }_0)}^{\alpha }}·\end{array}$$(56)

Subtracting from the global extragalactic background signal the contribution by extragalactic radio sources, for instance assuming the model by Gervasi et al. (2008a) with an amplification factor of ≃1.3 in the resulting background (see also Sect. 10.5), and comparing the residual background with their almost complete collection of cosmic background absolute temperature data, they found: ν₀ ≃ 0.4 GHz for α = 3 (and T^⋆ = T₀ to fit high frequency data), ν₀ ≃ 0.35 GHz, and α ≃ 3.36 using both ν₀ and α as fit variables.

In Fig. 1, we show ΔT_th, ΔR and the coefficients a_ℓ, 0(ν, β) for ℓ = 1, 6, expressed in terms of Δa_ℓ, 0, derived for the two sets of best-fit parameters according to the solutions in Sect. 3 and in one case, also on the basis of Eq. (5). The computations were performed in quadruple precision. We first carried out some tests with a simple Gaussian quadrature scheme (Press 1992), using various accuracy parameter values and point numbers (e.g., with the accuracy parameter (EPS) set to 10⁻⁹ and 2048 points), and compare the result with the explicit analytical solution for a blackbody: we find unreliable results above ℓ = 3 (or above ℓ = 4 provided that Eq. (5) is written as in the last equality). We then performed the numerical integration using the very accurate and efficient routine D01AJF of the Numerical Algorithms Group (NAG) Numerical Library, available only in double precision, setting integration accuracy parameters to the smallest values and increasing the number of sub-intervals used by the routine and the related workspace allocation. However, we verified that splitting the integral in terms of sums of integrals over subsets of the integration intervals does not improve the accuracy at all.

There is very good agreement between the results found with the routine D01AJF and the solutions in Sect. 3 with ℓ_max = 6 (particularly at lower multipoles, where the lines are superimposed and indistinguishable). Their differences are compatible with a combination of higher order terms, that is, beyond ℓ = 6, and integration errors, that are missing in the solutions in Sect. 3 and present in the numerical results, respectively. The two types of differences clearly appear, respectively, at lower frequencies, where the signal is higher and the relative integration error is lower, and at higher frequencies, where the signal is lower and the relative integration error is higher. We report also the nominal integration error quoted by the routine D01AJF: the comparison with the above differences suggests that this error is likely very conservative. In Appendix A, we provide some results derived adopting a much larger value of β that implies much larger signals, relatively higher contributions from higher multipoles, as well as relatively lower numerical integration errors: the analysis clearly supports the above interpretation.

Figure 1 shows that the typical power law shape of the intrinsic monopole spectrum, subsequent to the subtraction of the blackbody at the present temperature T₀, is displayed also at higher multipoles, as already discussed in Trombetti & Burigana (2019) for the dipole. Remarkably, we find that the ratio, R, between observed and intrinsic monopole, ΔR in Fig. 1, is not frequency independent, as in the case of a blackbody, but exhibits a frequency dependence related to the assumed intrinsic monopole spectrum. At low frequencies, below ∼1 GHz, the values of ΔR are positive and with amplitudes comparable to |R^BB − 1| or even larger.

10.2. Comptonization distortion and free-free diffuse emission

Many types of sources of photon and energy injections in cosmic plasma generates Comptonization distortions (Zel’dovich et al. 1972) via electron heating, and in ionizing the matter, they also produce FF distortions. Although these signatures can be generated both before (see e.g., Chluba & Sunyaev 2012) and after the cosmological recombination epoch (see e.g., Stebbins & Silk 1986; Danese & Burigana 1994), the cosmological reionization associated with the early formation phases of bound structures is the most remarkable source of these distortions. Two key parameters quantify the amplitudes of these imprints that for a given model, are tightly coupled. They are the Comptonization parameter, u, proportional to the global fractional energy exchange between matter and radiation in the cosmic plasma (for small distortions u ≃ (1/4)Δε/ε_i), and the FF distortion parameter, y_B(x), defined by integrals over the relevant redshift interval. On the other hand, even in the context of the reionization process, a variety of astrophysical mechanisms can contribute to determine the final distortion levels (see e.g., De Zotti et al. 2016; Burigana et al. 2018, and references therein). The resulting distorted photon distribution function is well approximated by

$$\begin{array}{c}\hfill {\eta }^{\mathrm{FF}+\mathrm{C}}\simeq {\eta }_{\mathrm{i}}+u\frac{x/{\varphi }_{\mathrm{i}}{e}^{x/{\varphi }_{\mathrm{i}}}}{{(e^{x/{\varphi }_{\mathrm{i}}}-1)}^2}(\frac{x/{\varphi }_{\mathrm{i}}}{\tanh [x/(2{\varphi }_{\mathrm{i}})]}-4)+\frac{y_B(x)}{x^3},\end{array}$$(57)

where η_i is the photon occupation number at the dissipation process initial time denoted with the subscript i. Neglecting other processes, η_i can be assumed to have a Planckian distribution at the initial temperature defined by ϕ_i = ϕ(z_i) = (1 + Δε/ε_i)^−1/4 ≃ 1 − u, that is, η_i = 1/(e^x/ϕ_i − 1).

The global Comptonization distortion depends linearly on matter density, thus, assuming a uniform medium is not critical for computing u. Conversely, bremsstrahlung depends quadratically on matter density and, in the presence of a substantial intergalactic medium (IGM) matter density contrast, the FF distortion is amplified with respect to the case of a homogeneous medium (Cooray & Furlanetto 2004; Ponente et al. 2011; Trombetti & Burigana 2014) by a factor of ≃1 + σ²(z), that is, ${\mathrm{\Omega }}_{\mathrm{b}}^2(z)\to {\mathrm{\Omega }}_{\mathrm{b},\mathrm{homog}}^2(1+{\sigma }^2(z))$, with σ²(z) as the baryonic matter variance related to the thermal properties of the DM particles.

Following Trombetti & Burigana (2019), we consider two pairs of different FF and Componization distortion models to identify a plausible range of possible distortions. We first consider the ionization history of Gnedin (2000), resulting in a Thomson optical depth τ that is fully consistent with recent Planck results, along with a fixed cut-off value k_max = 100. We coupled it with two different levels of Comptonization distortion, characterized by u = 10⁻⁷, which is very close to that derived in Burigana et al. (2008) for the Gnedin (2000) model and corresponds to an almost minimal energy injection consistent with the current constraints on τ, and by u = 2 × 10⁻⁶, a value that accounts for possible additional energy injections by a broad set of astrophysical phenomena.

At long wavelengths, λ = c/ν ≳ λ₀ = 1.5 cm, y_B is well-described by a linear dependence on log λ, ${\mathit{y}}_B^{\mathrm{lin}}\simeq a^{\mathrm{lin}}\log \lambda +b^{\mathrm{lin}}$, while at λ = c/ν ≲ λ₀ a quadratic dependence, ${\mathit{y}}_B^{\mathrm{quad}}\simeq a^{\mathrm{quad}}{\log }^2\lambda +b^{\mathrm{quad}}\log \lambda +c^{\mathrm{quad}}$, works better. The coefficients a^lin and b^lin are given in Appendix C of Trombetti & Burigana (2014): for the adopted model a^lin ≃ 3.292 × 10⁻⁹, b^lin ≃ 2.070 × 10⁻⁹. To allow for continuous derivatives of y_B also at frequencies around the transition between the two regimes, thus avoiding to introduce spurious oscillations in the resulting a_ℓ, 0, we need to properly join the two representations. Combining them with exponential weights,

$$\begin{array}{c}\hfill y_B=y_B^{\mathrm{lin}}(1-e^{-{(\lambda /{\lambda }_0)}^d})+y_B^{\mathrm{quad}}{e}^{-{(\lambda /{\lambda }_0)}^d},\end{array}$$(58)

with d = 3/2, is suitable to this purpose. A best-fit (see Table C.1 of Trombetti & Burigana 2014) gives a^quad ≃ −4.657 × 10⁻¹⁰, b^quad ≃ 1.210 × 10⁻⁹, c^quad ≃ 2.841 × 10⁻⁹.

Larger FF distortions are expected from the integrated contribution of an ensemble of ionized halos at substantial redshifts, as in the model by Oh (1999) that predicts a value of y_B ∼ 1.5 × 10⁻⁶ at ν ∼ 2 GHz. We then consider a second pair of models rescaling the above FF representation to y_B(2 GHz) = 1.5 × 10⁻⁶, coupled with Comptonization distortions with u = 10⁻⁷ or 2 × 10⁻⁶.

In general, the frequency behaviour of y_B is significantly less model dependent than its overall amplitude. A power law representation of y_B with a single amplitude parameter is adopted, for simplicity, in Sect. 14. In this approximation

$$\begin{array}{c}\hfill y_B(\nu )\simeq A_{\mathrm{FF}}{(\nu /\mathrm{GHz})}^{-\zeta },\end{array}$$(59)

where, assuming ζ ≃ 0.15, suitable values of A_FF at low frequencies (several MHz ≲ ν ≲ some GHz) are respectively A_FF ≃ 7.012 × 10⁻⁹ and 1.664 × 10⁻⁶ (Trombetti & Burigana 2019). Assuming the same slope but values of A_FF multiplied by a proper factor, f, offers a reasonable approximation also at 30 GHz ≲ ν ≲ 100 GHz (on average we find f ∼ 0.578, while it ranges between ≃0.550 and ≃0.596). Of course, a better power law fit can be found jointly varying A_FF and ζ according to the considered frequency range. As examples, for the first model at 10 MHz ≲ ν ≲ 1 GHz (or at 30 GHz ≲ ν ≲ 100 GHz) we find A_FF ≃ 7.225 × 10⁻⁹ and ζ ≃ 0.143 (or A_FF ≃ 5.236 × 10⁻⁹ and ζ ≃ 0.214). The adopted intrinsic monopole models are shown in Fig. 2 (top-left panel) in terms of ΔT_th.

The results in Fig. 2 are derived multiplying y_B in Eq. (58) by a damping function, exp(−(λ/λ₁)^d)exp(−(λ/λ₂)^d), which is relevant at very short wavelengths (λ₁ = 0.09 cm and λ₂ = 0.05 cm, corresponding to ≃333 GHz and 600 GHz) to make the results at ν ≳ 400 GHz dependent essentially only on the Comptonization term. This does not appreciably affect the results shown in the various panels of Fig. 2 at ν ≲ 400 GHz. On the other hand, while a better theoretical characterization of the FF emission at very high frequencies is required for a proper estimate in this context, we note that at ν ≳ 400 GHz the signal associated to the CIB, discussed in Sect. 10.7, dominates over the other terms at any multipole.

The differences, Δa_ℓ, 0, derived for these models are positive at low frequencies, where the FF term dominates, and negative at high frequencies, where the Comptonization prevails. The transition from the FF to the Comptonization regime, which ranges from about 3 GHz to about 300 GHz, depends on the relative amplitude of the two parameters y_B and u. This generalizes the result already found by Trombetti & Burigana (2019) for the dipole: in particular, the transition frequency between the two regimes clearly increases with ℓ, increasing from a maximum value of ∼100 GHz for ℓ = 1 to a maximum value around ∼350 GHz for the highest values of ℓ. Again, the approximate power law shape of the intrinsic monopole spectrum at low frequencies is maintained at higher multipoles.

The ratio between observed and intrinsic monopole is frequency-dependent (see top-right panel of Fig. 2) and at low frequencies, that is, below ≈0.1 GHz, ΔR can be comparable in amplitude to |R^BB − 1| or even larger, mainly depending on the level of the FF diffuse emission.

10.3. Bose-Einstein-like distortion

Bose-Einstein-like distortions can be produced by a variety of early processes, including unconventional heating sources, which could occur before the end of the phase of kinetic equilibrium between radiation and matter. Under near-equilibrium conditions, the stationary solution of the standard Kompaneets equation including only Compton scattering is a Bose-Einstein (BE) photon distribution function (Sunyaev & Zel’dovich 1970):

$$\begin{array}{c}\hfill {\eta }_{\mathrm{BE}}=\frac{1}{e^{x_{\mathrm{e}}+\mu }-1},\end{array}$$(60)

with a frequency independent chemical potential, μ; here x_e = x/ϕ(z), ϕ(z) = T_e(z)/T_r = ϕ_BE(μ), with T_e(z) the electron temperature. For mechanisms intrinsically involving a negligible photon number density production or absorption, μ is related to the fractional energy exchanged in the plasma during the interaction, Δε/ε_i, where the subscript i denotes the process initial time. For small distortions, ϕ_BE ≃ (1−1.11 μ)^−1/4 and, for an almost instantaneous process, μ ≃ 1.4Δε/ε_i at the end of the dissipation phase. Photon production processes, such as bremsstrahlung and double (or radiative) Compton emission, are particularly efficient at low frequencies, making the chemical potential dependent on the frequency, μ = μ(x) (Sunyaev & Zel’dovich 1970; Illarionov & Sunyaev 1974; Danese & De Zotti 1980). In combination with photon diffusion by Compton scattering, they tend to decrease the value of μ.

At high frequencies, x_e ≳ 1, the relaxation to a BE stationary solution can be achieved for processes that have occurred at redshifts of z ≳ z_p, corresponding to a time Comptonization parameter of y_e ≳ y_p ≃ 1/4 (Danese & De Zotti 1980; Burigana et al. 1991a), where $y_{\text{e}}={\int }_1^{1+z}[t_{\text{exp}}/t_{\text{C}}][d(1+z^{\prime })/(1+z^{\prime })]$ and where t_exp is the cosmic expansion time and t_C = [1/(n_eσ_Tc)][m_ec²/(kT_e)] is the timescale for the achievement of the kinetic equilibrium, with n_e as the density of free electrons, m_e the electron mass, and σ_T the Thomson cross-section. At z ≲ z_p, if the dissipation mechanism is concluded, the evolution of the photon distribution function is mainly due to photon production processes that significantly affect the low-frequency spectral region up to the recombination epoch. Aside from this effect, a longer time is needed for the photon distribution function relaxation towards the final spectrum at low frequencies, x_e ≲ 1; it can be achieved for processes occurred at redshifts z ≳ z₁ corresponding to y_e ≳ y₁, with y₁ ≃ 5 for small distortions (Burigana et al. 1991a). Here, we define, using μ₀, the high frequency asymptotic value of μ at z₁, which substantially identifies the end of the kinetic equilibrium phase. For the above reasons, the observational constraints on the chemical potential are typically referred to μ₀, the constraints on μ = μ(z) at higher redshifts being derived theoretically through (semi-)analytical formulas or numerical methods (see e.g., Burigana et al. 1991b), according to the considered problem. The limits on μ(z) can be significantly relaxed at increasing redshifts, and the constraints on Δε/ε_i before the thermalization redshift (when even large distortions can be erased) are then set by cosmological nucleosynthesis. The current upper limit on μ₀ is mainly derived from FIRAS data at λ ≲ 1 cm, |μ₀|< 9 × 10⁻⁵ at the 95% confidence level (Fixsen et al. 1996), although jointly recovering early and late spectral distortion parameters and including measurements at longer wavelengths can marginally change this constraint (see e.g., Nordberg & Smoot 1998; Salvaterra & Burigana 2002; Gervasi et al. 2008b; Seiffert et al. 2011).

While a value of μ₀ ≃ few × 10⁻⁵ cannot be excluded by current data, the existence of much smaller BE-like distortions is predicted as a consequence of two unavoidable processes. The dissipation of primordial perturbations at small scales (Hu et al. 1994; Chluba et al. 2012a), generates a positive chemical potential with values of μ₀ between ∼10⁻⁹ and 10⁻⁷, mainly depending on the shape of spectrum of the primordial scalar perturbation, a wider range being achieved in some inflation models varying the amplitude of primordial perturbations at very small scales (Chluba et al. 2012b) that are not constrained by current CMB anisotropy data. The faster decrease of the matter temperature relative to the radiation temperature in an expanding Universe generates, instead, a negative chemical potential, because of the interaction of CMB photons with colder electrons, with an absolute value ≃3 × 10⁻⁹ (Chluba & Sunyaev 2012; Sunyaev & Khatri 2013).

According to the above discussion, we consider just three values of μ₀: 1.4 × 10⁻⁹, 1.4 × 10⁻⁵, and −2.8 × 10⁻⁹. We adopt here an updated implementation of the semi-analytical representation of BE-like distortions, suitable also at low frequencies, as proposed by Danese & De Zotti (1980) and described in detail in Burigana et al. (1995). We assume a cold DM plus cosmological constant (ΛCDM) model with the set of parameters based on the last Planck data release and derived in Planck Collaboration VI (2020) including CMB power spectra in combination with CMB lensing reconstruction (see their Table 2, column labelled ‘TT, TE, EE+lowE+lensing’). We adopt a universe with a Hubble constant, H₀ = 67.36 km s⁻¹ Mpc⁻¹, cosmological constant (or dark energy) and non-relativistic matter density parameters Ω_Λ = 0.6847, Ω_m = 0.3153, Ω_b [H₀/(100 km s⁻¹ Mpc⁻¹)]² = 0.02237 (implying a baryon density Ω_b = 0.0493017), and, according to the standard model, an effective number of relativistic neutrinos N_eff = 3.046. In principle, the fine accounting of the relativistic neutrinos contribution to the expansion rate in the presence of an energy injection should also require the specification of the heating redshift, z_h, particularly for z_h ≫ z₁ when μ could be significantly larger than μ₀ (Burigana et al. 1991a): for simplicity, we treat this aspect in numerical estimates as in the case equivalent to z_h ≃ z₁. We also assume ϕ(z_h) ≃ ϕ(z₁) = ϕ_BE(μ₀). The cosmic expansion time, t_exp, and the relevant rates depend on these parameters, that play the major role in determining the spectrum shape. We compute the bremsstrahlung term according to Karzas & Latter (1961), Rybicki & Lightman (2008) and Burigana et al. (1991a) but using, in its range of validity, the polynomial fitting formula for the non-relativistic exact Gaunt factor derived by Itoh et al. (2000) (see also Chluba et al. 2020a for recent improvements). We separately compute the contributions from ionized hydrogen (H⁺) and helium (He⁺⁺ and He⁺), counting accordingly the overall fraction of free electrons, given the helium mass fraction (f_He = 0.2454). We calculate the double Compton rate in the elastic limit according to Lightman (1981) and Thorne (1981), and using the cross-section by Gould (1984).

Considering the relevance of double Compton at high redshifts, we also include the correction factor, C_mr, for mildly relativistic thermal plasma in the soft photon limit, C_mr ≃ 1/[1 + 14.6kT_r/(m_ec²)], introduced by Chluba et al. (2007). In Appendix B, we provide a fitting formula that, in the limit of very small distortions, can be used to compute the double Compton Gaunt factor at a precision level better than ≃0.1−0.2% also at x ≳ 1. Finally, replacing the simple approximation of full ionization up to the hydrogen recombination with the introduction of the redshifts (∼6 × 10³ and ∼2 × 10³) at which He⁺⁺ and He⁺ disappear, resulting into a ∼two-steps helium recombination, introduces only a small correction (≲0.2%) in the final spectrum computation. The above details enter in the computation of the key redshifts z₁ and z_p (respectively, ≃5.38 × 10⁵ and ≃1.21 × 10⁵ with the adopted parameters), in the frequency-dependent optical depth of the universe for absorption, $y_{\text{abs}}={\int }_1^{1+z}[t_{\text{exp}}/t_{\text{abs}}][\text{d}(1+z^{\prime })/(1+z^{\prime })]$ (see Burigana et al. 1995), and in the characteristic dimensionless frequency, x_c, which quantifies the low-frequency damping of the chemical potential, μ(x) = μ₀ exp (−x_c/x_e), x_c being defined by t_abs(z₁) = t_C(z₁), where t_abs is the absorption timescale for photon production processes. For our purposes, a simple Gaussian quadrature scheme is accurate and efficient enough for computing the relevant integrals over z (we find advantageous to work with a logarithmic integration variable), while the NAG routine D01AJF can provide a better performance. The Brent’s method (Press 1992) is suitable to solve the equation for x_c (found to be ≃4.86 × 10⁻³ with the adopted parameters), given bracketing guesses based on the simple low frequency limit approximation.

The intrinsic monopole spectra so obtained for the three adopted values of μ₀ are shown in Fig. 3 (top-left panel) in terms of ΔT_th. As expected, the amplitude of |ΔT_th| is proportional to |μ₀| and, for μ₀ <  0, we find a spectrum shape, ΔT_th, opposite in sign with respect to the case μ₀ >  0. It is interesting to note the plateau at extremely low frequencies and the presence of two characteristic changes of sign in ΔT_th, corresponding to the well-known excess (or decrement) of signal of the BE-like spectrum with respect to the blackbody at temperature T₀ at low and high frequencies and the remarkable decrement (or excess) at intermediate frequencies for positive (or negative) values of μ₀. These sign changes also appear in the differences Δa_ℓ, 0 (see Fig. 3), but at two characteristic frequencies significantly increasing at increasing ℓ. We note that the sign change at the higher of the two characteristic frequencies occurs at ν >  1 THz for ℓ ≥ 4 and just for this reason it disappears in the corresponding plots. In addition, two further sign changes appear at each increase of an even multipole. They are located at frequencies between the smaller of the two above characteristic frequencies and the plateau at extremely low frequencies. Again, this pattern of sign changes is symmetric with respect to the sign of μ₀. These are remarkable features of the BE-like spectrum: they are almost independent of the value of μ₀, while their behaviour at low frequencies depend on the underlying cosmological parameters.

In Fig. 3, ΔR is displayed in the top-right panel. Over the whole frequency range, it is characterised by a module proportional to |μ₀| and much smaller than |R^BB − 1| even for values of |μ₀| not far from FIRAS limits. Again, ΔR depends on frequency. For BE-like distortions, the shape of ΔR is again symmetric with respect to the sign of μ₀, with a well-defined maximum (or minimum) located between two minima (or two maxima), clearly defined at lower frequencies and less pronounced at higher frequencies, for μ₀ >  0 (or μ₀ <  0).

10.4. 21 cm redshifted line

The 21 cm line corresponds to the spin-flip transition in the ground state of neutral hydrogen. This signal is described as the offset of the 21 cm brightness (i.e. antenna) temperature from the background temperature, T_back, along the observed line of sight at a frequency ν that, because of cosmic expansion, is related to the rest frame frequency, ν_21 cm = c/(21 cm), by ν = ν_21 cm/(1 + z). T_back is usually assumed equal to T_r but, in general, it could include potential distortions and other radiation backgrounds. It depends on the evolution of the gas spin temperature, T_S, which represents the excitation temperature of the 21 cm transition, of the fraction of neutral hydrogen, on the Hubble function, H(z), on Ω_m, on the matter density contrast, and on the comoving gradient of the line-of-sight component of the comoving velocity (Furlanetto et al. 2006). If T_S <  T_back (or T_S >  T_back), the gas is seen in absorption (or in emission). Since the signal detected at a given frequency corresponds to a specific redshift, the 21 cm line provides a tomographic view of the cosmic evolution. A rich set of 21 cm redshifted line models has been studied in Cohen et al. (2017), resulting in an wide envelope of predictions for $T_{\mathrm{ant}}^{21\mathrm{cm}}(\nu )$.

In this work, we consider only the pronounced absorption profile, with an almost symmetric U-shape centred at (78 ± 1) MHz, which was recently found by Bowman et al. (2018) based on an analysis of the data from the Experiment to Detect the Global EoR Signature (EDGES). The absorption feature was found to have an amplitude of $0.5_{-0.2}^{+0.5}$ K and a spread of the profile with a full width at half-maximum of ${19}_{-2}^{+4}$ MHz. These data support the presence of an ionizing background by 180 million years after the Big Bang and a phase of gas heating above the radiation temperature less than 100 million years later (Bowman et al. 2018). The explanation of the 21 cm redshifted line signal found by EDGES might require a substantial cooling of the IGM gas or an additional high-redshift extragalactic radio background, or a combination of them (see e.g., the references at the beginning of Sect. 10.1). The authors provide a suitable analytical representation of the EDGES absorption profile in terms of a flattened Gaussian characterized by a set of best-fit parameters, and we adopt here their expression of $T_{\mathrm{ant}}^{21\mathrm{cm}}(\nu )$.

The signals considered in Sects. 10.1–10.3, referring to intrinsic CMB spectral distortions, already include the contribution of the unperturbed CMB spectrum. In this section (as well as in Sects. 10.5–10.7), we are considering signals that are superimposed onto the CMB, assumed without spectral distortions, and we then add (in terms of η) the CMB blackbody at the current temperature T₀ to $T_{\mathrm{ant}}^{21\mathrm{cm}}(\nu )$ to construct the global signal to be studied as in previous sections.

We compute the coefficients a_ℓ, 0(ν, β) with the method described in Sect. 3 and, for comparison, also on the basis of Eq. (5), as in Sect. 10.1. The results are shown in Fig. 4.

We note the consistency between the results based on the integral given by Eq. (5) and our approach: the agreement is excellent up to ℓ_max = 6 (the green and red the lines are indistinguishable). Again, their differences are compatible with a combination of higher order terms, that is, beyond ℓ = 6, and integration errors, that are only present in the numerical results. The two types of differences clearly appear, respectively, where the signal is higher and the relative integration error is lower, and where the signal is lower and the relative integration error is higher. The latter point is also evident from the comparison with the nominal integration error quoted by the routine D01AJF. Remarkably, except for the numerical integration uncertainty, the spectral shape of the differences is very similar for all the odd, as well as the even multipoles because they mainly come from the contribution from ℓ = 7 and from ℓ = 8, respectively. In Appendix A, we repeat this analysis adopting a much larger value of β, that implies much larger signals, relatively higher contributions from higher multipoles and relatively lower numerical integration errors: the result clearly supports the above interpretation.

It is interesting to note the complexity, increasing with ℓ, of the features displayed in Fig. 4: they include the alternation of increasing and decreasing behaviours, the number of relative minima and maxima, and the changes of the sign of Δa_ℓ, 0.

Here, the result for Δa_1, 0 corrects the dipole spectrum published in Trombetti & Burigana (2019), expressed there in terms of ${\mathrm{\Delta }}_{0,\pi /2}{T}_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$ (see Eq. (52)), where, for EDGES, the profile $T_{\mathrm{ant}}^{21\mathrm{cm}}(\nu )$ was accounted in equivalent thermodynamic (not in antenna) temperature. This mere oversight significantly affected only the very small values of ${\mathrm{\Delta }}_{0,\pi /2}{T}_{\mathrm{th}}^{\mathrm{BB}/\mathrm{dist}}$, making the wings a bit steeper they are.

We note that the frequency, right above 70 MHz, corresponding to the change of sign of Δa_1, 0 is shifted at slightly larger values in the case of Δa_3, 0 and Δa_5, 0, while for the even multipoles it corresponds to a well-defined minimum that falls in the middle of a positive interval of the Δa_ℓ, 0 profile. The size of this frequency interval decreases as ℓ increases; the same holds, in the case of odd multipoles, for the interval identified by the two signs changes around the above minimum. Also, the increasing with ℓ of the number of sign changes in Δa_ℓ, 0 implies that, at frequencies outside the above interval, the size of each frequency range with unchanged sign of Δa_ℓ, 0 decreases with ℓ. These are remarkable features of the considered model.

Remarkably, ΔR (see top-right panel of Fig. 4), again exhibits a frequency dependence related to the assumed intrinsic monopole spectrum: it is positive in the inner frequency range, almost corresponding to the plateau of the absorption feature, and negative in the profile wings. In particular, |ΔR| is comparable to or, typically, greater than |R^BB − 1|, across most of the relevant frequency range.

10.5. Extragalactic radio background

An important radio background is produced by extragalactic sources. Differently from most of the signals discussed in previous sections, which are of intrinsically diffuse origin, this background, as well as those discussed in Sects. 10.6 and 10.7, results from the integrated contribution of discrete sources. With galactic surveys able to reach increasingly deeper flux density levels, it is possible to resolve a large fraction of this background. In spite of this, an observationally diffuse residual background comes from faint sources below the survey detection limits. Other than intrinsically interesting, this extragalactic background needs to be accurately known in order to understand the reionization imprints correctly. Remarkably, these classes of signals may also be tightly related. A notable extragalactic radio background is evident in, for example, the radio data by Dowell & Taylor (2018) and it was proposed by Seiffert et al. (2011) to explain the signal excess claimed by ARCADE 2.

We exploit several simple analytical representations of the extragalactic radio background. According to Trombetti & Burigana (2019), we assume the best-fit power-law model,

$$\begin{array}{c}\hfill T_{\mathrm{ant}}^{\mathrm{Back}}(\nu )\simeq 18.4\mathrm{K}{(\nu /0.31\mathrm{GHz})}^{-2.57},\end{array}$$(61)

by Seiffert et al. (2011).

A careful analysis and prediction of the extragalactic source radio background between 0.151 and 8.44 GHz, also including different source detection thresholds, was carried out by Gervasi et al. (2008a). We consider their best-fit single power-law model for the extragalactic source background signal, multiplied by a factor of 1.3 in order to approximately account for a larger contribution that is, ultimately, likely ascribed to the emerging of star-forming galaxies and radio-quiet active galactic nuclei at fainter flux densities. Indeed, the Lockman Hole Project and deep Low Frequency Array (LOFAR) imaging of the Boötes field support a certain flattening of differential number counts, N′(ν), at 1.4 GHz below ≈100 μJy (Prandoni et al. 2018) and at 0.15 GHz below ≈1 mJy (Retana-Montenegro et al. 2018). This may suggest an increase in N′(ν) of a factor of ∼2 with respect to the estimate of Gervasi et al. (2008a) at the faint flux densities and a ∼30% increase in the extragalactic radio background, which is proportional to ${\int }_{S_{\text{min}}}^{S_{\text{max}}}S{N}^{\prime }(\nu )\text{d}S$. We then adopt

$$\begin{array}{c}\hfill T_{\mathrm{ant}}^{\mathrm{Back}}(\nu )\simeq 1.3\times 0.88\mathrm{K}{(\nu /0.61\mathrm{GHz})}^{-2.707}.\end{array}$$(62)

Gervasi et al. (2008a) provided also an empirical analytical fit function of N′(ν) that can be used to estimate the remaining residual extragalactic radio background when a certain source detection threshold, S_max, is assumed. According to Trombetti & Burigana (2019), we exploit their differential number counts assuming S_max = 50 nJy, which almost corresponds to typical detection limits of the ultra deep reference continuum surveys planned for the Square Kilometre Array (SKA; Prandoni & Seymour 2015). We consider the above factor of ∼2 to be applied to N′(ν) by Gervasi et al. (2008a) at faint flux densities (we then label this case as a ‘High-radio background residual’) and fit the results found in the frequency range considered by the authors to find the corresponding estimate of the remaining residual extragalactic radio background,

$$\begin{array}{c}\hfill T_{\mathrm{ant}}^{\mathrm{Back}}(\nu )\simeq A{(\nu /\mathrm{GHz})}^{-2.65},\end{array}$$(63)

with A ≃ 4.7 mK. Different choices of S_max mainly reflects into the value of A.

As in previous section, we add (in terms of η) the CMB blackbody at the current temperature T₀ to these radio background models to construct the global signal. The intrinsic monopole spectra are shown in Fig. 5 in terms of ΔT_th (top-left panel) with the derived coefficients a_ℓ, 0(ν, β).

The typical power law shapes of the considered intrinsic monopole spectra, after the subtraction of the blackbody at the present temperature T₀, are also maintained at higher multipoles and the same holds for their relative amplitudes, as already noted in Trombetti & Burigana (2019) for the dipole.

Again, ΔR exhibits a frequency dependence that is related to the assumed intrinsic monopole spectrum. At low frequencies, below ∼1 GHz, ΔR can have an amplitude comparable to |R^BB − 1| or even larger.

10.6. Extragalactic millimetre background

The extragalactic radio source populations that mainly contribute to the radio and to the millimetre background are very different. Whereas steep-spectrum radio sources, particularly at high flux densities, are most important at radio frequencies, extragalactic compact sources with an almost flat, or possibly inverted, spectrum (primarily blazars, flat spectrum radio quasars, and BL Lacertae sources, where are not fully considered in Sect. 10.5) become increasingly relevant at wavelengths shorter than a few centimetres. They can be directly extracted by analyzing CMB maps. The products (Planck Collaboration VII 2013; Planck Collaboration XXVI 2016; Planck Collaboration XLV 2016; Planck Collaboration LIV 2018) from the Planck mission, complemented by the available ground-based data (see e.g., Mocanu et al. 2013) and, in the far-IR, via Herschel (see e.g., López-Caniego et al. 2013) observations, provide crucial information for the characterization of their number counts (De Zotti et al. 2015) (see also De Zotti et al. 2005; Tucci et al. 2011), while a substantial progress at fainter flux densities is expected from the next-generation space missions and deeper multi-frequency ground-based surveys (see e.g., De Zotti et al. 2018, and references therein).

We exploit their differential number counts, N′(ν), including both steep and flat spectrum sources, to estimate the corresponding millimetre background, $\propto {\int }_{S_{\text{min}}}^{S_{\text{max}}}S{N}^{\prime }(\nu )\text{d}S$, for which we find a power-law approximation, as in Sect. 10.5. We simply extrapolate the power law behaviours of the differential number counts at the lowest available flux density ranges to fainter flux densities. The (currently uncertain) characterization of N′(ν) at much fainter flux densities is relatively less relevant for the estimate of the global background than for the estimate of the residual millimetre background derived given a certain source detection threshold S_max. We consider two assumptions of S_max, namely: (1) 10 mJy (labeled as ‘Intermediate millimetre background residual’), a value that is similar to those typically considered for next-generation space missions; and (2) 100 μJy (labelled as ‘Low millimetre background residual’), representing an estimate of the potential improvement achievable with SKA while considering frequency extrapolation uncertainties from deep radio surveys to millimetre bands.

As expected due to the spectral shapes of the sources that are more relevant at millimetre wavelengths, we find a background spectrum slightly flatter in the millimetre than in the radio, with slopes (of the spectrum expressed in antenna temperature) between ≃ − 2.10 and −2.62, depending on the specific number counts model (with steeper spectra found for the number counts by Tucci et al. 2011) and also on the adopted source detection threshold.

For current numerical estimates we consider, for simplicity, a single slope based on number counts consistent with Planck results from 30 GHz to 857 GHz. We adopt:

$$\begin{array}{c}\hfill T_{\mathrm{ant}}^{\mathrm{mm}}(\nu )\simeq A{(\nu /\mathrm{GHz})}^{-2.19},\end{array}$$(64)

with A ≃ 32.4 mK for the millimetre background spectrum and ≃17.3 mK or ≃8.93 mK for residual millimetre background derived for S_max = 10 mJy or S_max = 100 nJy, respectively. We construct the global signal by adding (in terms of η) the CMB blackbody at the current temperature T₀.

Our results are shown in Fig. 6. The flattening of ΔT_th at ν ∼ 200 GHz and its increasing at higher frequencies are due to the representation in terms of equivalent thermodynamic temperature. We note the different behaviours of Δa_ℓ, 0 at odd and even multipoles, and the (minimal) increase with ℓ of the frequency where for odd ℓ, the change of sign of Δa_ℓ, 0 occurs, or for even ℓ, the minimum of Δa_ℓ, 0 is located.

Figure 6 reports the coefficients a_ℓ, 0(ν, β) for ℓ = 1, 6, expressed in terms of Δa_ℓ, 0, derived using the solutions given in Sect. 3. The typical power law shapes of the considered intrinsic monopole spectra, after the subtraction of the blackbody at the present temperature, T₀, that are evident in the figure at low frequencies, are kept also at higher multipoles and the same holds for their relative amplitudes, as already noted in Trombetti & Burigana (2019) for the dipole.

ΔR exhibits a dependence on the frequency as well as on the assumed intrinsic monopole spectrum: it is positive and its frequency shape and minimum location are in line with the behaviours of Δa_ℓ, 0 at even multipoles. On the other hand, the values of |ΔR| are less than |R^BB − 1| even at the highest frequencies.

Finally, we note that, despite the fact that the estimate of the extragalactic radio source background in the radio and at millimetre wavelengths, approximated respectively by Eqs. (62) and (64), are based on different models, the results found at ∼20 GHz, a frequency in the middle between the maximum and the minimum frequency of the approximations elaborated in the previous section and in this one, agree within a factor of two or better, as shown by the comparison between Figs. 5 and 6.

10.7. Cosmic infrared background

A precise analysis of CIB spectrum, which is still not well known, can provide a better understanding of the dust-obscured star-formation phase of galaxy evolution.

In spite of its absolute calibration precision of 0.57 mK, the FIRAS characterization of CIB amplitude and shape still presents a substantial uncertainty. According to Fixsen et al. (1998), a suitable analytic representation in terms of photon distribution function of the CIB spectrum at the present time, can be expressed by

$$\begin{array}{cc}\hfill {\eta }_{\mathrm{CIB}}& =\frac{c^2}{2h{\nu }^3}{I}_{\mathrm{CIB}}(\nu )=\frac{I_0{(\nu /{\nu }_0)}^{k_{\mathrm{F}}}}{e^{x_{\mathrm{CIB}}}-1}=\frac{I_0{(x_{\mathrm{CIB}}/x_{0,\mathrm{CIB}})}^{k_{\mathrm{F}}}}{e^{x_{\mathrm{CIB}}}-1},\hfill \end{array}$$(65)

where ν₀ = c/(100 μm) ≃ 3 THz and x_CIB = hν/(kT_CIB). The best-fit to FIRAS data gives k_F = 0.64 ± 0.12, T_CIB = (18.5 ± 1.2) K and I₀ = (1.3 ± 0.4) × 10⁻⁵ (Fixsen et al. 1998), I₀ setting the CIB spectrum amplitude. In the last equality of Eq. (65), x_0, CIB = hν₀/(kT_CIB) ≃ 7.78 (implying x_CIB ≃ 7.78 × ν/ν₀). The current uncertainty on CIB spectrum amplitude is quite high, with 1σ accuracy of about 30%. Indeed, the direct determination of the CIB spectrum is hard to obtain as it requires absolute intensity measurements and is limited by foreground signals.

We construct the global signal by adding η_CIB with the photon distribution function of the CMB blackbody at the current temperature T₀ and consider three simple cases corresponding to the above best-fit and 1σ limits of I₀. The results are shown in Fig. 7 in terms of ΔT_th, ΔR, and Δa_ℓ, 0. As shown by the comparison of Fig. 7 with Figs. 1–6 and anticipated in Sect. 10.2, the signal associated to the CIB, strongly increasing with frequency, at ν ≳ 400 GHz dominates over the other extragalactic contributions at any multipole.

As has already been found for extragalactic millimetre background, we note the different behaviours of Δa_ℓ, 0 at odd and even multipoles and, in the frequency range between ≃100 GHz and ≃200 GHz, the increase with ℓ of the frequency where, for odd ℓ, the change of sign of Δa_ℓ, 0 occurs or, for even ℓ, the shape of Δa_ℓ, 0 shows a clear steepening. At the highest frequencies, where the power law approximation of the intrinsic monopole spectrum breaks down and the spectrum changes its behaviour approaching its maximum, the shapes of Δa_ℓ, 0 show remarkable features. They are located at frequencies that decrease as ℓ increases: for example, for ℓ = 2 or 6 they occur at frequencies larger than ∼800 GHz or than ∼500 GHz.

Here, ΔR (top-right panel of Fig. 7) assumes positive values that are much smaller than |R^BB − 1| at ν ≲ 400 GHz, but significantly increase with ν at ν ≳ 400 GHz and become comparable to or larger than |R^BB − 1| at ν ≳ 700 GHz, achieving a maximum at a frequency ≃800−850 GHz that slightly increases as I₀ decreases.

11. Signal combination results

As already states previously in this paper, the method described in Sects. 2–9 and applied in Sect. 10 to specific emissions can be also used to predict the signatures expected from the desired combinations of signals, provided that they are summed in terms of additive quantities. A discussion of the imprints on the dipole spectrum left by combinations of backgrounds associated to cosmological reionization and relevant in the radio can be found in Trombetti & Burigana (2019). Of course, the number of the possible combinations of the models discussed in Sect. 10 is high. Here, we consider a couple of cases that are relevant at millimetre and sub-millimetre wavelengths.

We combine the model of Comptonization plus diffuse FF distortion with the highest values of u and y_B with an astrophysical background, namely, the adopted millimetre background model from extragalactic radio sources, an estimate of its residual signal given an assumption of source detection threshold (the ‘Low millimetre background residual’), and the best-fit CIB spectrum. The global signal is constructed by adding the photon distribution functions associated to the considered astrophysical background (Eqs. (64) or (65)) and to the Comptonization plus diffuse FF distortion (Eq. (57)), that already contains the initial unperturbed CMB spectrum, η_i. The results based on the solutions given in Sect. 3 are shown in Figs. 8 and 9 in terms of ΔT_th, ΔR, and Δa_ℓ, 0.

From these figures (see also Figs. 2, 6, and 7), it is possible to appreciate that at the frequencies where a certain component is much stronger than the other, that component also dominates in the combined signal, but this simplification does not hold where the components have comparable amplitudes.

The comparison of Fig. 8 with Figs. 2 and 6 shows that the millimetre background, as well as its residual for the assumed source detection threshold, dominates above a frequency of ∼500 GHz, with only a very little dependence on ℓ. Below a frequency ranging from few tens of GHz to ∼100 GHz the millimetre background and the considered model of Comptonization plus diffuse FF distortion give comparable effects: in particular, the millimetre background and the diffuse FF distortion have also similar power law behaviours (see Eqs. (59) and (64)). As a result, the difference of considering the millimetre background or its residual is more evident in the spectral shapes of their combinations, affecting, respectively, less or more the frequency where the spectra exhibit their sign change, that occurs at frequencies slightly higher than in the case of pure Comptonization plus diffuse FF distortion. The Comptonization distortion is, instead, clearly appreciable at intermediate frequencies, its characteristic plateau appearing in all the Δa_ℓ, 0.

The comparison of Fig. 9 with Figs. 2 and 7 shows analogous results. The CIB dominates over the Comptonization distortion above a frequency that ranges from ∼100 GHz to ∼200 GHz, slightly increasing with ℓ. Below that frequency, the Comptonization distortion emerges and the combined spectral shapes become flatter. They show also a remarkable change of sign at even ℓ from Figs. 2–6, where Δa_ℓ, 0 is positive for the CIB and negative for the Comptonization distortion. This does not occur for ΔT_th where the contributions from CIB is slightly larger than that from Comptonization, and for ΔR, that is positive in both the cases. At further decreasing frequencies, the FF diffuse emission emerges, and, consequently, the Δa_ℓ, 0 exhibit, at even and odd ℓ, the typical shapes already found in that case at the lowest frequencies as well as the typical change of sign corresponding to the transition from the range dominated by the FF term to that dominated by the Componization term. The frequency of this transition is only slightly larger than that found in the case of pure Comptonization plus diffuse FF distortion. This frequency shift is significantly smaller than that found above, which combines the Comptonization plus diffuse FF distortion with the millimetre background, because of the smaller additional contribution from the CIB at low frequencies.

Of course, the details in the above considerations depend also on the assumed model parameters. We discuss above some cases where the combined contribution from different signals may give not trivial effects. In general, refinements in modelling specific signals and combinations or couplings of them can be included in our method to improve the quality of the results.

12. Maps and angular power spectra

Having evaluated the coefficients a_ℓ, 0, it suffices then to compute the corresponding maps and angular power spectra. Here, we report just a few examples, also for comparison with the results presented in Burigana et al. (2018), which are based on map generation and inversion performed with the great set of publicly available tools implemented in the Hierarchical Equal Area and isoLatitude Pixelization of the sphere (HEALPix; Górski et al. 2005). Here, we adopt the same pixelization scheme to generate the maps using the a_ℓ, 0 coefficients: working in real space, based on Eq. (6), we immediately get the map for each multipole component in the reference system with the z axis parallel to the observer velocity. They can be also simply computed in any other reference system (e.g., in Galactic coordinates) given the angle of a pixel direction with respect to the observer velocity direction. Here, we generate the maps with n_side = 1024, corresponding to a pixel linear size of ≃3.44 arcmin.

We show the maps for ℓ = 3 and 6 in the case of a BE-like distortion with the maximum value of μ₀ (see Sect. 10.3), of a pure Comptonization distortion with the maximum value of u (see Sect. 10.2) and of the best-fit CIB spectrum (see Sect. 10.7). The three adopted frequencies, different for each type of signal, are selected to allow us to make an almost direct comparison with the maps displayed in Burigana et al. (2018) and to appreciate how the method presented here is suitable for a fast map computation, very precise up to the desired order, even where the relevant signal is weak and, in principle, more sensitive to numerical uncertainty. The maps (see Fig. 10) are indeed very clean and without numerical artifacts up to the highest ℓ thanks to the adopted analytical approach; each multipole pattern obviously reflects the corresponding scaling, see Eq. (6), and, in particular, for odd multipoles the maximum and minimum values of each map are equal in module.

Given the map, we can use publicly available tools, such as the anafast facility of HEALPix, to compute the corresponding angular power spectrum and the a_ℓ, m coefficients in any reference system (see Goldstein 1984, and also Appendix C, for explicit formulas to transform the a_ℓ, m under rotation). Indeed, for the considered problem, the a_ℓ, m with m ≠ 0 do not vanish in reference systems with the z axis not parallel to the observer velocity direction. On the other hand, since the angular power spectrum, C_ℓ, is an invariant under rotation of the reference system, we can compute it simply working in the reference system with the z axis parallel to the observer velocity direction, where the a_ℓ, m with m ≠ 0 vanish, and the coefficients a_ℓ, 0 are evaluated as in the previous sections

$$\begin{array}{c}\hfill C_{\ell }(\nu ,\beta )={\displaystyle \sum _{m=-\ell }^{m=\ell }\frac{a_{\ell ,m}^2(\nu ,\beta )}{2\ell +1}=\frac{a_{\ell ,0}^2(\nu ,\beta )}{2\ell +1}·}\end{array}$$(66)

In Fig. 11, we only report the result derived for some of the considered models for ℓ from 1 to 4. We plot the quantity:

$$\begin{array}{c}\hfill \frac{\mathrm{\Delta }{a}_{\ell ,0}^2(\nu ,\beta )}{2\ell +1}=\frac{{[a_{\ell ,0}^{\mathrm{dist}}(\nu ,\beta )-a_{\ell ,0}^{\mathrm{BB}}(\beta )]}^2}{2\ell +1}·\end{array}$$(67)

Equation (67) gives the angular power spectrum of the map obtained from the difference between the map produced in a given model (i.e. for a CMB distorted spectrum or for an astrophysical background spectrum summed with the blackbody spectrum at the current temperature T₀) and the map obtained for the blackbody spectrum at the current temperature T₀.

There is a good agreement with the results reported by Burigana et al. (2018) – but in the frequency range between 60 and 600 GHz considered for the Cosmic Origins Explorer (COrE) in Fig. 14 (for the CIB) and in Fig. 12 (for BE and Comptonization distortions). As expected, the central panels of Fig. 11, that span from the radio to the sub-millimetre, show that the replacement of the pure BE spectrum with a BE-like spectrum and the inclusion of FF diffuse emission, which is not included in Burigana et al. (2018), is remarkable at lower frequencies (given the high FF model considered here, a little effect is already appreciable at the lowest frequencies of the COrE range, see also Fig. 2).

The signals more relevant at low and high frequencies are displayed in terms of $\mathrm{\Delta }{a}_{\mathcal{l},0}^2/(2\mathcal{l}+1)$ in the top and bottom panels of Fig. 11, respectively. All the well-defined minima of the various lines displayed in Fig. 11 correspond to the changes of sign of Δa_ℓ, 0 (see Figs. 2–4 and 6–7), but not in the case of the minimum at ≃80 MHz with respect to the EDGES profile at ℓ = 2 and 4 (appearing in the top panels of Fig. 11, which correspond to the positive minima at the same ℓ’s in Fig. 4).

Finally, we remark that the values of $\mathrm{\Delta }{a}_{\mathcal{l},0}^2/(2\mathcal{l}+1)$ at the minima corresponding to the sign changes of Δa_ℓ, 0 (as well as the values of Δa_ℓ, 0) should, in principle, go to zero at the corresponding frequencies; typically, this does not appear in the plots merely because of their frequency discretization.

13. Varying observer velocity

As discussed in Sect. 1, the observer velocity β = β_C + β_V(t) is modulated by the time dependent component β_V that, in real measurements, mainly comes from the revolution of the Earth or of L2 around the Solar System barycentre (β_V ≃ β_ES ≃ β_L2 ≃ 10⁻⁴ ≃ 0.081 β). The formalism and the general properties of the solutions presented in previous sections hold for any choice of β. Thus, we can consider a relatively small time interval, Δt_i, around a certain instant t_i, in which the time variation of β_V can be neglected and β can be considered as a constant. For example, in a time interval Δt_i ≃ 1 h, the relative change of β, which is given by Δβ/β ≃ (Δβ_V/β_V)(β_V/β) ≃ ω Δt_i (β_V/β) ≃ 0.0058%, is about one order of magnitude smaller than the current relative uncertainty on β_C, given the period of the main modulation P = 2π/ω = 1 yr.

In a time interval, Δt_i, only a certain fraction of the sky can be observed. The detailed pattern of the sequence of the observed sky positions is defined by the so-called scanning strategy, or observational strategy, which is specific of each considered experiment. Let us define, using $M_{\mathrm{i}}(\nu ,{\theta }_{\mathrm{i}}^{\prime },{\mathit{\beta }}_{\mathrm{i}})$, the all-sky map given by Eq. (6), where β is replaced by β_i = β(t_i) = β_C + β_V(t_i) and the colatitude θ is replaced by the angle ${\theta }_{\mathrm{i}}^{\prime }={\theta }_{\mathrm{i}}^{\prime }({\mathit{\beta }}_{\mathrm{i}},\theta ,\varphi )=\arccos [(\widehat{n}·{\mathit{\beta }}_{\mathrm{i}})/{\beta }_{\mathrm{i}}]$ between β_i and a sky direction corresponding to the unit vector $\widehat{n}$ defined by the colatitude θ and the longitude ϕ in any adopted reference system, S. The spherical harmonic coefficients to be used in Eq. (6) are given in Sect. 3 for ℓ_max = 6 (or, alternatively, in Sect. 5 for ℓ_max = 4 or in Sect. 7 for ℓ_max = 2 or 1).

The experiment scanning strategy defines the sequence of pointing maps, P_i(t_i, θ, ϕ), on a pixelized sky at the times t_i, with i = 1, n and n = τ_S/Δt_i where τ_S is the survey duration: if the direction defined by θ and ϕ is observed in the interval Δt_i then P_i = 1, otherwise P_i = 0. The averaged map corresponding to a survey based on a set of observer velocities β_i is then given by:

$$\begin{array}{c}\hfill M(\nu ,\theta ,\varphi )={\displaystyle \sum _{i=1}^n\frac{M_{\mathrm{i}}(\nu ,{\theta }_{\mathrm{i}}^{\prime },{\mathit{\beta }}_{\mathrm{i}})·P_{\mathrm{i}}(t_{\mathrm{i}},\theta ,\varphi )}{N(\theta ,\varphi )},}\end{array}$$(68)

where N(θ, ϕ) is the global number of observer velocities β_i for which the direction defined by θ and ϕ is observed.

As discussed in Sect. 2, we need to compute the signal in only ℓ_max + 1 directions to predict the map up to ℓ_max, and the number of signal evaluations required by Eq. (68) is then N_ev, ℓmax = (ℓ_max + 1) n = (ℓ_max + 1) τ_S/Δt_i. Of course, the map M(ν, θ, ϕ) can be directly computed through Eqs. (1) and (2) using β = β_C + β_V(t) at each time, without applying the discretization associated to the choice of the time interval Δt_i. In this case, since the receiver sampling time, τ_samp, is typically much smaller than Δt_i, the required number of signal evaluations, N_ev, d = τ_S/τ_samp, is much larger than N_ev, ℓmax. The ratio, r_ev = N_ev, ℓmax/N_ev, d, between the required number of signal evaluations in the two cases is r_ev = τ_samp(ℓ_max + 1)/Δt_i ≪ 1. In order to properly compare the computing times in the two approaches, we evaluate this ratio discretizing β = β_C + β_V(t) according to the choice of Δt_i also in the direct scheme based on Eqs. (1) and (2). Thus, N_ev, d reduces to N_ev, d = N_i, pixτ_S/Δt_i, where N_i, pix is the number of sky pixels observed in the time interval Δt_i, and we obtain r_ev = (ℓ_max + 1)/N_i, pix (again, r_ev ≪ 1).

The method based on Eq. (68) is particularly advantageous when applied to scanning strategies designed for future CMB missions that foresee the use of a huge number of receivers at the focal surface of a wide field of view telescope. The scanning strategies typically foreseen for these missions involve fast spacecraft spin axis precessions, in order to achieve in a short time a large sky coverage and to observe sky pixels with many orientations during the survey to improve the quality of map making results and of polarization analyses (see e.g., Delabrouille et al. 2018; Natoli et al. 2018).

Since ${\theta }_{\text{i}}^{\prime }$ depends on β_i, θ and ϕ, Eq. (68) implies that the time varying component β_V(t) will introduce in the map M(ν, θ, ϕ) a further modulation along θ and ϕ, which is superimposed onto the main pattern, even if we choose a reference frame S with the z axis parallel to the constant component β_C. Although the details of this modulation depend on the experiment scanning strategy, its amplitude amounts to ∼β_V/β ∼ 8.1% of the main signal due to the constant velocity component.

For a given model, fixing the other parameters and neglecting the terms in β^(ℓ+2) and beyond, the differences of the spherical harmonic coefficients Δa_ℓ, 0 scale proportionally to β^ℓ, that is, Δa_ℓ, 0(β) ≃ Δa_ℓ, 0(β_C) (β/β_C)^ℓ at the leading order. This property comes from the separation of the system into even and odd multipoles (see Sects. 2–7) and is explicit in the formulas found for the blackbody spectrum (see Sects. 4 and 6). The comparison of the results displayed in Figs. 1 and 4 with those reported in Figs. A.1 and A.2 where the velocity is multiplied by a factor 100 and 10, respectively, also suggests this scaling, providing that the signal is sufficiently smooth in frequency and except for smearing effects where the signal rapidly changes. We verified this property in our numerical results varying β within a ±10% because of the contribution of β_V, namely, for realistic deviations of β from β_C. This feature can further speed up the computation of $M_{\mathrm{i}}(\nu ,{\theta }_{\mathrm{i}}^{\prime },{\mathit{\beta }}_{\mathrm{i}})$ because it allows a fast estimation of the contribution of the time varying component of the velocity over two consecutive multipole patterns, for example, in the case of the dipole and quadrupole, as useful in many applications.

14. Global pattern

In general, the global sky pattern is a combination of intrinsic anisotropies and of the anisotropies induced by the peculiar observer motion. We focus here on the diffuse cosmic signals more relevant in the microwaves, where the intrinsic (mainly of primordial nature) anisotropies are better studied and the background frequency spectrum can be modelled in terms of small deviations from a blackbody.

Working within a reference frame at rest with respect to the background, in a given sky direction identified by θ and ϕ the frequency dependent equivalent thermodynamic temperature can be seen as a function of an effective temperature and of a set of P distortion parameters p_j with j = 0, P − 1

$$\begin{array}{c}\hfill T_{\mathrm{th}}(\nu ,\theta ,\varphi )=T_{\mathrm{th}}(\nu ,T_{\mathrm{BB}}(\theta ,\varphi ),p_j(\theta ,\varphi )).\end{array}$$(69)

Since both fluctuations and distortions are small, T_BB(θ, ϕ) ∼ T₀ and p_j(θ, ϕ) ∼ 0 and we can expand T_th in Taylor’s series around these values. At linear order

$$\begin{array}{cc}\hfill T_{\mathrm{th}}(\nu ,\theta ,\varphi )& ={(T_{\mathrm{th}})}_0+{\left(\frac{\partial T_{\mathrm{th}}}{\partial T_{\mathrm{BB}}}\right)}_0(T_{\mathrm{BB}}-T_0)+{\displaystyle \sum _{j=0}^{P-1}{\left(\frac{\partial T_{\mathrm{th}}}{\partial p_j}\right)}_0{p}_j}\hfill \\ \hfill & =T_0+{\left(\frac{\partial T_{\mathrm{th}}}{\partial \eta }\frac{\partial \eta }{\partial T_{\mathrm{BB}}}\right)}_0(T_{\mathrm{BB}}-T_0)\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{P-1}{\left(\frac{\partial T_{\mathrm{th}}}{\partial \eta }\frac{\partial \eta }{\partial p_j}\right)}_0{p}_j,}\hfill \end{array}$$(70)

where ()₀ denotes that the quantities are evaluated at T_BB(θ, ϕ) = T₀ and p_j(θ, ϕ) = 0, that is, for a blackbody with effective temperature T₀. Expanding the temperature fluctuation (T_BB(θ, ϕ) − T₀) and the distortion parameters p_j(θ, ϕ) in spherical harmonics with coefficients a_{BB, ℓ, m} and a_{j, ℓ, m}, we have

$$\begin{array}{cc}\hfill T_{\mathrm{th}}(\nu ,\theta ,\varphi )=& a_{\mathrm{BB},0,0}{Y}_{0,0}+{\displaystyle \sum _{j=0}^{P-1}{a}_{j,0,0}{Y}_{0,0}{\left(\frac{\partial T_{\mathrm{th}}}{\partial \eta }\frac{\partial \eta }{\partial p_j}\right)}_0}\hfill \\ \hfill & +{\displaystyle \sum _{\ell =1}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{a}_{\mathrm{BB},\ell ,m}{Y}_{\ell ,m}{\left(\frac{\partial T_{\mathrm{th}}}{\partial \eta }\frac{\partial \eta }{\partial T_{\mathrm{BB}}}\right)}_0}\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{P-1}\sum _{\ell =1}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{a}_{j,\ell ,m}{Y}_{\ell ,m}{\left(\frac{\partial T_{\mathrm{th}}}{\partial \eta }\frac{\partial \eta }{\partial p_j}\right)}_0,}\hfill \end{array}$$(71)

where a_BB, 0, 0Y_0, 0 = T₀ and $a_{j,0,0}{Y}_{0,0}={\overline{p}}_j$ are the average of T_BB(θ, ϕ) and p_j(θ, ϕ) over the full sky. Although not necessary, we adopt here for simplicity a reference system with the z axis parallel to the observer velocity in order to avoid rotations in the following considerations (see Eq. (84)). In the right-hand side of Eq. (71), the first line represents the (possibly distorted) monopole spectrum, the second line the intrinsic temperature fluctuations, the third line the intrinsic fluctuations of the distortion parameters and (see Eq. (1))

$$\begin{array}{c}\hfill \frac{\partial T_{\mathrm{th}}}{\partial \eta }=\frac{h\nu }{k}\frac{1}{\eta (1+\eta ){\ln }^2(1+1/\eta )},\end{array}$$(72)

implying that (see Eq. (21)):

$$\begin{array}{c}\hfill {\left(\frac{\partial T_{\mathrm{th}}}{\partial \eta }\right)}_0=T_0\frac{{(e^x-1)}^2}{x{e}^x},\end{array}$$(73)

and the functions (∂η/∂p_j)₀ (j = 0, P − 1) depend on the type of distortion. Equation (71) shows that the spectral shape of each distortion term in the monopole spectrum (added to T₀) and of the fluctuations of the corresponding distortion parameter is the same, but they are weighted differently depending on the coefficients a_{j, ℓ, m}.

Next, we consider, in the case of distorted spectra, the derivatives ∂η/∂p_j and ∂η/∂T_BB. In the case of a BE distortion (see Eq. (60)) with a frequency independent chemical potential μ₀ (i.e. neglecting for simplicity the spectrum modifications introduced by considering a BE-like distortion that are relevant at lower frequencies), and for small values of μ₀, with the approximation ϕ_BE ≃ (1−1.11μ₀)^−1/4 we have

$$\begin{array}{cc}& \frac{\partial \eta }{\partial {\mu }_0}\simeq \frac{e^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{BE}}+{\mu }_0}}{{(e^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{BE}}+{\mu }_0}-1)}^2}·(\frac{1.11}{4}\frac{h\nu }{k{T}_{\mathrm{BB}}{\varphi }_{\mathrm{BE}}}{\varphi }_{\mathrm{BE}}^4-1),\hfill \end{array}$$(74)

$$\begin{array}{cc}& \frac{\partial \eta }{\partial T_{\mathrm{BB}}}\simeq \frac{1}{T_{\mathrm{BB}}}\frac{e^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{BE}}+{\mu }_0}h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{BE}}}{{(e^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{BE}}+{\mu }_0}-1)}^2},\hfill \end{array}$$(75)

implying

$$\begin{array}{cc}& {\left(\frac{\partial \eta }{\partial {\mu }_0}\right)}_0\simeq \frac{e^x[(1.11/4)x-1]}{{(e^x-1)}^2},\hfill \end{array}$$(76)

$$\begin{array}{cc}& {\left(\frac{\partial \eta }{\partial T_{\mathrm{BB}}}\right)}_0\simeq \frac{1}{T_0}\frac{x{e}^x}{{(e^x-1)}^2}·\hfill \end{array}$$(77)

We derive now ∂η/∂p_j and ∂η/∂T_BB in the case of a Comptonization distortion with a small Comptonization parameter u and an initial Planckian spectrum, η_i with ϕ_i ≃ 1 − u (see Eq. (57)), combined with a FF distortion with y_B approximated by ${\mathit{y}}_B\simeq {\stackrel{\sim }{A}}_{\mathrm{FF}}{x}^{-\zeta }$ with ${\stackrel{\sim }{A}}_{\mathrm{FF}}=A_{\mathrm{FF}}{[(h/k)(\mathrm{GHz}/T_{\mathrm{BB}})]}^{\zeta }\simeq 0.5456A_{\mathrm{FF}}$ (see Eq. (59)). In principle, also the slope parameter ζ or, for the more general description in Eq. (58), an alternative set of five parameters could be included in the set of distortion parameters p_j. For simplicity, we include in the p_j only the most relevant emission amplitude parameter in the power law approximation. In the derivatives below, we report only the terms that will not vanish when specified at p_j = 0. The terms that are not multiplied by u in ∂η/∂u give

$$\begin{array}{cc}\hfill \frac{\partial \eta }{\partial u}=& -\frac{[h\nu /(k{T}_{\mathrm{BB}})]e^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{i}}}/{\varphi }_{\mathrm{i}}^2}{{(e^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{i}}}-1)}^2}\hfill \\ \hfill & +\frac{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{i}}{e}^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{i}}}}{{(e^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{i}}}-1)}^2}(\frac{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{i}}}{\tanh [h\nu /(k{T}_{\mathrm{BB}})/(2{\varphi }_{\mathrm{i}})]}-4),\hfill \end{array}$$(78)

and we get

$$\begin{array}{c}\hfill {\left(\frac{\partial \eta }{\partial u}\right)}_0=-\frac{x{e}^x}{{(e^x-1)}^2}+\frac{x{e}^x}{{(e^x-1)}^2}(\frac{x}{\tanh (x/2)}-4),\end{array}$$(79)

while

$$\begin{array}{c}\hfill \frac{\partial \eta }{\partial {\stackrel{\sim }{A}}_{\mathrm{FF}}}\simeq {\left(\frac{h\nu }{k{T}_{\mathrm{BB}}}\right)}^{-(3+\zeta )},\end{array}$$(80)

implying

$$\begin{array}{c}\hfill {\left(\frac{\partial \eta }{\partial {\stackrel{\sim }{A}}_{\mathrm{FF}}}\right)}_0\simeq x^{-(3+\zeta )}.\end{array}$$(81)

The terms that are not multiplied by u nor by ${\stackrel{\sim }{A}}_{\mathrm{FF}}$ in ∂η/∂T_BB give

$$\begin{array}{c}\hfill \frac{\partial \eta }{\partial T_{\mathrm{BB}}}\simeq \frac{1}{T_{\mathrm{BB}}}\frac{e^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{i}}}h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{i}}}{{(e^{h\nu /(k{T}_{\mathrm{BB}})/{\varphi }_{\mathrm{i}}}-1)}^2},\end{array}$$(82)

implying, again,

$$\begin{array}{c}\hfill {\left(\frac{\partial \eta }{\partial T_{\mathrm{BB}}}\right)}_0\simeq \frac{1}{T_0}\frac{x{e}^x}{{(e^x-1)}^2}·\end{array}$$(83)

Formally, in Eqs. (74), (75), (78), (80), and (82), T_BB, μ₀, ϕ_BE(μ₀) and ϕ_i(u) refer to a given sky direction identified by θ and ϕ. However, their relations with θ and ϕ vanish when they are specified for T_BB = T₀ and p_j = 0 as in Eqs. (76), (77), (79), (81), and (83), where indeed x = hν/(kT₀).

In general, we have ((∂T_th/∂η)(∂η/∂T_BB))₀ = 1, as evident also from Eqs. (72), (77), and (83); the coefficients a_{BB, ℓ, m} are then the usual spherical harmonic expansion coefficients defining the CMB temperature anisotropies. Considering both the effects induced by the observer peculiar motion on the monopole and the intrinsic anisotropies, working in the observer reference system with the z axis parallel to the observer velocity, we can combine Eqs. (3) and (71) to derive a global anisotropy pattern

$$\begin{array}{cc}\hfill T_{\mathrm{th}}(\nu ,\theta ,\varphi ,\beta )=& {\displaystyle \sum _{\ell =1}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{a}_{\mathrm{BB},\ell ,m}{Y}_{\ell ,m}{\left(\frac{\partial T_{\mathrm{th}}}{\partial \eta }\frac{\partial \eta }{\partial T_{\mathrm{BB}}}\right)}_0}\hfill \\ \hfill & +{\displaystyle \sum _{j=0}^{P-1}\sum _{\ell =1}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{a}_{j,\ell ,m}{Y}_{\ell ,m}{\left(\frac{\partial T_{\mathrm{th}}}{\partial \eta }\frac{\partial \eta }{\partial p_j}\right)}_0}\hfill \\ \hfill & +{\displaystyle \sum _{\ell =0}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{a}_{\ell ,m}(\nu ,\beta )Y_{\ell ,m}(\theta ,\varphi ),}\hfill \end{array}$$(84)

where the coefficients a_ℓ, m(ν, β) are given in Sect. 3 (or in Sect. 5); we omitted here the first line of the right-hand side of Eq. (71) to obviously avoid a double counting of the monopole that is already included in the term with ℓ = 0 of the last line, where the effect of the observer motion with respect to a frame at rest with the CMB is also taken into account. The focus in this work on the analysis at low multipoles of the effects caused by the observer motion on the isotropic monopole component, neglects the Doppler and aberration effects on the anisotropies, that is, on the first and second lines of the right-hand side of Eq. (84). These effects couple multipoles ℓ to ℓ ± n, particularly in the correlation between ℓ and ℓ ± 1 (Challinor & van Leeuwen 2002; Burles & Rappaport 2006; Kosowsky & Kahniashvili 2011; Amendola et al. 2011; Chluba 2011; Dai & Chluba 2014), a property that has been used to independently constrain β. The effects are indeed more important at high multipoles: their main information comes from ℓ ≳ 100, where many modes can be exploited (Burigana et al. 2018) (see Planck Collaboration LVI 2020 for a recent application to Planck data based on the modulation of the thermal Sunyaev-Zel’dovich effect (Chluba et al. 2005; Notari & Quartin 2016)).

We can rewrite the coefficients a_ℓ, m(ν, β) as

$$\begin{array}{c}\hfill a_{\ell ,m}(\nu ,\beta )=\mathrm{\Delta }{a}_{\ell ,m}(\nu ,\beta )+a_{\ell ,m}^{\mathrm{BB}}(\beta ),\end{array}$$(85)

where $a_{\ell ,m}^{\mathrm{BB}}(\beta )$ refer to the case of a blackbody spectrum and $\mathrm{\Delta }{a}_{\ell ,m}(\nu ,\beta )=a_{\ell ,m}(\nu ,\beta )-a_{\ell ,m}^{\mathrm{BB}}(\beta )$ depend on the type of considered distortion.

15. Intrinsic dipole versus kinematic dipole

We focus here on the dipole anisotropy. Including both the effect induced by the observer peculiar motion on the monopole and the intrinsic anisotropies, and using a reference system with the z axis parallel to the observer velocity, the global dipole pattern is characterized by the coefficients:

$$\begin{array}{c}\hfill a_{1,m}^{\mathrm{glob}}(\nu ,\beta )=\mathrm{\Delta }{a}_{1,m}(\nu ,\beta )+a_{1,m}^{\mathrm{BB}}(\beta )+a_{\mathrm{BB},1,m}+{\displaystyle \sum _{j=0}^{P-1}{\stackrel{\sim }{a}}_{j,1,m}(\nu ),}\end{array}$$(86)

where

$$\begin{array}{c}\hfill {\stackrel{\sim }{a}}_{j,1,m}(\nu )=a_{j,1,m}{\left(\frac{\partial T_{\mathrm{th}}}{\partial \eta }\frac{\partial \eta }{\partial p_j}\right)}_0\end{array}$$(87)

and Δa_1, m(ν, β) = 0, $a_{1,m}^{\mathrm{BB}}(\beta )=0$ for m ≠ 0. In these equations, the superscript ‘BB’ refers to observer peculiar motion effects while the subscript ‘BB’ refers to intrinsic anisotropies. The coefficients Δa_1, 0(ν, β) and ${\stackrel{\sim }{a}}_{j,1,m}(\nu )$ (through ((∂T_th/∂η)(∂η/∂p_j))₀), that are related to the type of distortion, do not vanish only in the presence of deviations from a Planckian spectrum.

The typical amplitude of the coefficients Δa_1, 0(ν, β) is of the order of $a_{1,0}^{\mathrm{BB}}(\beta )$, that is, of ≈βT₀, multiplied by the amplitude of the monopole spectral distortion, ΔT_th (see the first two left panels from the top in Figs. 2 and 3). The coefficients a_j, 1, m characterize the fluctuations of the different types of distortion at ℓ = 1 and then depend significantly on the specific mechanism and not only on the corresponding average distortion parameter and spectral shape. Global spectral distortions are still unobserved, as well their fluctuations (obviously except for the Sunyaev-Zel’dovich effect on galaxy clusters), and it is then reasonable to assume ${\stackrel{\sim }{a}}_{j,1,m}(\nu )<a_{\mathrm{BB},1,m}$.

Assuming Gaussian random temperature fluctuations, the coefficients a_BB, 1, m are expected to have zero mean and variance given by the angular power spectrum C_ℓ at ℓ = 1. Currently, C₁ is unknown, but it is typically predicted to be of the order of the temperature anisotropy intrinsic quadrupole C₂. Constraining the intrinsic anisotropy power at ℓ = 1 is difficult, but very interesting in the context of future CMB surveys (see Yasini & Pierpaoli 2017 for a method based on the exploitation the leakage of the intrinsic dipole into the CMB monopole and quadrupole and Meerburg et al. 2017 for an analysis based on the observation of the small scale temperature fluctuations that result from gravitational lensing). This topic is related to the power at low multipoles, and in particular, to the low power of the quadrupole discovered by the COBE Differential Microwave Radiometer (DMR; Wright et al. 1996; Hinshaw et al. 1996) and then confirmed by WMAP (Bennett et al. 2013; Hinshaw et al. 2013) and Planck (Planck Collaboration V 2020). Reconstructing the intrinsic anisotropy power at very low multipoles is very important for inflationary models (see Planck Collaboration X 2020 for recent constraints) predicting power suppression at large scales (see e.g., Vilenkin & Ford 1982; Starobinskij 1992; Bridle et al. 2003; Contaldi et al. 2003; Efstathiou 2003; Sinha & Souradeep 2006) and for their connection with universe geometry and topology (see e.g., Linde 1995, 2003; Gratton et al. 2002; Ellis et al. 2002; Lasenby & Doran 2004; Levin 2002).

In general, the amplitudes of the coefficients a_BB, 1, m are significantly smaller than the amplitude of the coefficient $a_{1,0}^{\mathrm{BB}}(\beta )$, the observed dipole being dominated by the Doppler effect associated to our peculiar motion with respect to the CMB. In the absence of deviations from a Planckian spectrum, ${\stackrel{\sim }{a}}_{j,1,m}(\nu )=0$ and Δa_1, 0(ν, β) = 0, and any (relatively minor) frequency independent contribution from a_BB, 1, 0 is degenerate with $a_{1,0}^{\mathrm{BB}}(\beta )$. The same holds for a_BB, 1, m with m ≠ 0 because the sum of dipole terms is still a dipole and, for a Planckian spectrum, it is possible to find a rotation of the reference system that jointly drops the terms with m ≠ 0 at all frequencies (see Appendix C). Thus, the peculiar motion Doppler effect alone does not allow to distinguish the intrinsic dipole from the kinematic dipole, at least in the absence of a very accurate measure of β provided by other methods.

Now, we consider the presence of deviations from a Planckian spectrum. The only frequency dependent terms in Eq. (87) are Δa_1, 0(ν, β) and ${\stackrel{\sim }{a}}_{j,1,m}(\nu )$.

Let us assume, as a first simple case (A), that the very large angular scale fluctuations of the spectral distortion parameters are very small in amplitude, namely ${\stackrel{\sim }{a}}_{j,1,m}(\nu )\ll \mathrm{\Delta }{a}_{1,0}(\nu ,\beta )$. Thus, in the adopted reference system, the only relevant frequency dependence in the dipole pattern is for m = 0 and comes from Δa_1, 0(ν, β), while any different choice of the z axis, that is, not parallel (or not antiparallel) to the observer velocity, would imply that the same frequency dependence is polluted in the dipole at m ≠ 0.

In a more general case (B), when we relax the above assumption, that is, for a non-negligible frequency dependent contribution from ${\stackrel{\sim }{a}}_{j,1,m}(\nu )$, the situation is a bit more complex, but conceptually not so different (in particular for given prescriptions of the coefficients ${\stackrel{\sim }{a}}_{j,1,m}(\nu )$). In the adopted reference system, a frequency dependence related only to the fluctuations of the spectral distortions would appear at m ≠ 0, and the combination of the different frequency dependencies related to the peculiar motion effect and to the fluctuations of the spectral distortions, properly weighted, would appear at m = 0. For a reference system with the z axis not parallel (or not antiparallel) to the observer velocity these frequency dependencies are polluted in the dipole at any m, in a way related to the underlined frequency spectra.

Let us assume that the intrinsic dipole and the kinematic dipole are not aligned (the opposite is possible by chance, but it is very unlikely).

In case (A), one can search for a reference system that drops (or, more realistically, minimize in a statistical sense) the frequency dependence of the coefficients, $a_{1,m}^{\mathrm{glob}}$ for m ≠ 0, implying that its z axis is parallel to β. The components, $a_{1,m}^{\mathrm{glob}}$, with m ≠ 0 are then to be ascribed only to the intrinsic temperature fluctuation terms, a_BB, 1, m, and although the contribution from a_BB, 1, 0 remains hidden in the larger term, $a_{1,0}^{\mathrm{BB}}(\beta )$, the squares of the components, $a_{1,m}^{\mathrm{glob}}$, with m ≠ 0, allow us to provide an estimate of the intrinsic dipole angular power spectrum, C₁, although with a slightly larger cosmic variance because this estimate is based on only two, instead than three, coefficients. We note that (see Eq. (84)) similar considerations apply in the limit of that approximation, also for ℓ >  1, allowing us to exploit a larger number of modes m. On the other hand, because Δa_ℓ, m(ν, β) (see Eq. (85)) decreases as β^ℓ, the information from ℓ >  1 does not add relevant constraints in this scheme.

In case (B), relaxing the assumption ${\stackrel{\sim }{a}}_{j,1,m}(\nu )\ll \mathrm{\Delta }{a}_{1,0}(\nu ,\beta )$, it is possible to search for a reference system that drops, or minimize, the difference of the frequency dependence of the coefficients, $a_{1,m}^{\mathrm{glob}}$, with m ≠ 0, with the behaviour expected from the terms, ${\stackrel{\sim }{a}}_{j,1,m}(\nu )$, related to the fluctuations of the distortion parameters (see Eqs. (86) and (87)), added with the intrinsic temperature fluctuation terms, a_BB, 1, m. We also note that the frequency dependencies of ${\stackrel{\sim }{a}}_{j,1,m}(\nu )$ and Δa_1, 0(ν, β), although different, are physically connected, being related to the types of involved distortions, and this property can be exploited in the joint analysis of the modes with m = 0 and m ≠ 0, helping, at m ≠ 0, the discrimination between ${\stackrel{\sim }{a}}_{j,1,m}(\nu )$ and a_BB, 1, m.

Therefore, in the presence of spectral distortions, a very careful multifrequency analysis of the dipole pattern, namely, of the frequency behaviour of its spherical harmonic expansion coefficients, can be used, at least in principle, to set constraints on the intrinsic dipole embedded in the kinematic dipole.

We note that the inclusion of the time variation of β (see Sect. 13), although requiring specific implementations, does not conceptually modify the above considerations since it is possible to split the sky area observed in a survey into a proper set of sky areas, each one observed in a shorter time interval and, hence, with a negligible variation of β, without reducing the global statistical information contained in the survey.

16. Discussion and conclusion

The peculiar motion of an observer with respect to the cosmic background in a certain frequency band produces boosting effects in the background anisotropy pattern. In this work, we studied how the frequency spectrum of the background isotropic monopole emission is modified and transferred to the frequency spectra of the patterns at higher multipoles. We performed the analysis in terms of spherical harmonic expansion for various models of background radiation, ranging from the radio to the far-infrared.

Adopting a reference frame with the z axis parallel to the observer motion direction allows us to simplify the problem since it is thanks to this choice that only the spherical harmonic coefficients a_ℓ, m with m = 0 do not vanish. We derive the system of linear equations to obtain the a_ℓ, m up to a desired value of ℓ_max. For each observational frequency, the a_ℓ, m are written as linear combinations of the signals at the set of frequencies corresponding to the chosen N = ℓ_max + 1 colatitudes θ_i (Sect. 2). We explicitly write the system and provide the solutions up to ℓ_max = 6, as well as for other smaller values of ℓ_max. The symmetry property of the associated Legendre polynomials with respect to θ = π/2 is used to separate the system into two subsystems, one for ℓ = 0 and even multipoles and the other for odd multipoles, improving the solutions accuracy with respect to an arbitrary colatitudes choice (Sects. 3–8). We apply these solutions to the case of a blackbody and verify their agreement with the exact analytical solutions for ℓ = 0 and ℓ = 1 at the order in β related to ℓ_max (Sects. 4 and 6).

The structure of the solutions is discussed and compared with respect to the general properties of monopole spectrum integration and differentiation (Sect. 9). The coefficients of the solutions can be regarded as sets of weights assigned to a small number of function evaluations, according to adopted order of accuracy, to compute the integrals that define the a_ℓ, m coefficients. The coefficients of these solutions exhibit remarkable symmetry properties. Some of these properties are the same of the ones shown by the weights used in finite difference formulas to compute numerical derivatives. The implicit mixing of the spectrum derivatives in the solutions is reflected by the absence of that sign alternation property of the weights which appears in finite difference formulas. Indeed, the frequency behaviours of the a_ℓ, m coefficients are particular sensitive to the local monopole spectrum variation in a way characterized by derivative orders increasing with ℓ.

We applied the method to some models for different types of monopole spectra that can be represented in terms of analytical or semi-analytical functions (Sects. 10 and 11), namely: four types of CMB distorted photon distribution functions, that is, a non-equilibrium imprint at low frequencies, Comptonization, FF, and BE-like distortions; four types of extragalactic background signals, that is, 21 cm redshifted line, radio, and millimetre backgrounds from extragalactic sources, CIB, superimposed onto the CMB Planckian spectrum; some combinations of signal relevant at millimetre and sub-millimetre wavelengths, that is, Comptonization and FF distortions combined with millimetre background from extragalactic sources or with CIB.

For each model, we show the intrinsic monopole spectrum in terms of the difference, ΔT_th(ν), between the equivalent thermodynamic temperature for the model and the present CMB temperature T₀. Our results are presented in terms of: (i) the difference ΔR(ν, β) = R(ν, β) − R^BB(β), where $R(\nu ,\beta )=(a_{0,0}(\nu ,\beta )/\sqrt{4\pi })/T_{\mathrm{th}}(\nu )$ and R^BB(β) are the ratios between the equivalent thermodynamic temperature of observed and intrinsic monopole for the model and for the blackbody at the present temperature T₀; (ii) the differences Δa_ℓ, 0(ν, β) between the spherical harmonic coefficients computed for the model, a_ℓ, 0(ν, β), and for the blackbody, a_ℓ, 0(β); we also present (iii) all-sky maps and (iv) angular power spectra (Sect. 12) directly derived from the spherical harmonic coefficients, for some representative cases.

The results are in excellent agreement with those based on more computationally demanding numerical integrations or map generation and inversion (Sects. 10 and 12; see also Appendix A), and even more accurate. The method could be obviously implemented for any ℓ_max, however, since β is of the order of 10⁻³ and the coefficients a_ℓ, 0(ν, β) scale approximately as β^ℓ, the solutions presented for ℓ_max = 6 allows us to achieve an extremely high accuracy that is sufficient for any application even in the very distant future. We provide also explicit solutions for ℓ_max = 4 that are fully adequate for the analysis of forthcoming and proposed surveys. The only accuracy limitation of the proposed scheme derives from neglecting the contributions from higher multipoles, the largest relative errors appearing at ℓ_max from ℓ_max + 2 and at ℓ_max − 1 from ℓ_max + 1: the relative errors are very small even at ℓ ≲ ℓ_max, while, remarkably, they are fully negligible at the lowest multipoles.

The high number of possible model combinations may result in a variety of signatures in the global signal expected by cosmological plus astrophysical backgrounds. The simplicity and computational efficiency of the proposed method can significantly alleviate the computational effort needed for theoretical predictions and for the comparison with data from future projects for a plethora of cases of interest. These include, for example, signal combinations or possible differences between vectors β that could refer to specific backgrounds or that appear when the observer velocity variation is taken into account (Sect. 13).

We discuss the main features found for the considered background models at the various multipoles in wide frequency ranges. All the patterns at different multipoles are related to the observer peculiar motion, although with a signal amplitude that decreases with ℓ. Thus, they exhibit, both in the whole sky and in different sky regions, spatial correlations at different angular length scales and well-defined geometrical properties, that are related each other, and this would improve their joint analysis that can be optimized considering the observational method and specifications of a given project. As discussed in Trombetti & Burigana (2019) for the dipole in the radio domain, the analysis does not necessarily require the mapping of the entire sky or of a very large fraction of it. Future radio surveys, and in particular the excellent resolution and sensitivity offered by the SKA (Dewdney et al. 2016) for a variety of scientific themes (see e.g., Weltman et al. 2020) can be used to investigate the spectra of the multipoles patterns at low frequencies. In general, the ultra-precise comprehension and subtraction of the Galactic foreground emission is likely to act as the most difficult problem in cosmological analyses at large angular scales.

Finally, we move on to a discussion of the superposition of the CMB intrinsic anisotropies and of the effects induced by the observer peculiar motion focussing on their different frequency behaviours in the presence of CMB spectral distortions (Sect. 14). We find that they can be used, at least in principle, to set constraints on the intrinsic dipole embedded in the kinematic dipole, through a very careful multifrequency analysis of the corresponding spherical harmonic expansion coefficients (Sect. 15; see also Appendix C). Detailed studies, object of future works, can clarify the feasibility of this approach and, possibly, the required specifications.

In general, the above considerations do not strictly rely on the direct absolute determination of the monopole spectrum. Thus, although they come with challenges, these studies can, in principle, be pursued with an anisotropy project, such as the Lite (Light) satellite for the studies of B-mode polarization and Inflation from cosmic background Radiation Detection (LiteBIRD; Matsumura et al. 2014) or a mission like COrE (Delabrouille et al. 2018; de Bernardis et al. 2018) or the Probe of Inflation and Cosmic Origins (PICO; Hanany et al. 2019), provided that an extremely accurate relative and inter-frequency calibration and suppression of systematic effects are achieved (Natoli et al. 2018).

Clearly, the most advantageous observational chance is offered by the next generation of CMB space missions designed to perform ultra-accurate temperature measurements of the whole sky with a highly precise absolute calibration and a relatively high resolution mapping over a wide set of frequencies, as for example, in the Primordial Inflation Explorer (PIXIE; Kogut et al. 2011, 2016), in the Polarized Radiation Imaging and Spectroscopy Mission (PRISM; André et al. 2014), and, as in more recent proposals, to the National Aeronautics and Space Administration (NASA; Chluba et al. 2019a), to the Indian Space Research Organisation (ISRO; CMB Bharat¹), and to the European Space Agency (ESA; Delabrouille et al. 2019; Chluba et al. 2019b). These projects have the advantage of measuring the frequency spectrum of the relevant multipole patterns – starting from the monopole.

---

Acknowledgments

We gratefully acknowledge financial support from the research program RITMARE SP3 WP3 AZ3 U02 and the research contract SMO at CNR/ISMAR, from the INAF PRIN SKA/CTA project FORmation and Evolution of Cosmic STructures (FORECaST) with Future Radio Surveys and from the ASI/Physics Department of the University of Roma–Tor Vergata agreement n. 2016-24-H.0 for study activities of the Italian cosmology community. We gratefully acknowledge the use of the NAG Numerical Library. Some of the results in this paper have been derived using the HEALPix (Górski et al. 2005) package. We also thank the anonymous referee for comments that helped improve the paper.

References

Appendix A: Numerical tests with amplified observer velocity

In Sects. 10.1 and 10.4, we compare, for two very different cases, the results based on the solutions described in Sect. 3 and on a direct numerical integration (see Eq. (5)) to find differences that are compatible with a combination of higher order terms, that is, beyond ℓ = 6, and integration errors, the latter becoming more relevant at increasing ℓ and decreasing signals, in connection with the low value of the observer velocity, β.

To better clarify this aspect, we perform the same type of comparison, but adopting a much larger value of β, in order to exploit much larger signals and relatively higher contributions from higher multipoles and to deal with relatively lower numerical integration errors. Here, we reconsider the cases of the non-equilibrium model and of the EDGES profile of redshifted 21 cm line, but assuming a value of β arbitrarily amplified respectively by a factor 100 and 10, this choice being motivated by their different signal amplitudes in their relevant frequency ranges.

We report the results for Δa_ℓ, 0 found using the equations in Sect. 3 and the Eq. (5) for ℓ = 1, 6, as well as the results based on Eq. (5) for ℓ = 7 and 8. They are shown in Figs. A.1 and A.2. As expected, and as evident from the figures, such higher values of β imply a strong reduction of the relative numerical integration error, making feasible an accurate computation of Δa_ℓ, 0 at larger ℓ. Thus, the differences at ℓ = 1, 6 between the results found with the two methods come essentially only from higher order terms, neglected in the equations of Sect. 3, that are obviously dominated by the contributions from ℓ = 7 at odd multipoles and from ℓ = 8 at even multipoles.

We note that the order of magnitude of Δa_7, 0 (or of Δa_8, 0) is equal to the order of magnitude of the differences between the results found with the two methods at odd (or even) multipoles. The spectral shapes of Δa_ℓ, 0 in the case of the non-equilibrium model are featureless, reflecting the original power law shape of the monopole spectrum when expressed in terms of ΔT_th. Conversely, the spectral shapes of Δa_ℓ, 0 in the case of the EDGES profile are rich in features, increasing in number at increasing ℓ. Remarkably, the spectral shape of the differences between the results found with the two methods reflects the spectral shape of Δa_7, 0 (or of Δa_8, 0) at odd (or even) multipole.

As discussed in Sect. 2, the above properties derive from the separation of the system in two subsystems, one at odd ℓ the other at even ℓ, associated to the adoption of a set of colatitudes symmetric with respect to π/2 (plus π/2). This suppresses in each a_ℓ, 0 at odd (or even) ℓ the contribution of the multipole immediately larger than maximum odd (or even) multipole considered in the system solution (see also, for comparison, the discussion in Sect. 7). This is particularly important in practice for the ‘robustness’ of the accuracy of the solutions presented in this work, suitable for the real (low) value of β, and especially in the case of spectra rich in features for which changes of sign of Δa_ℓ, 0 and extremely low values of Δa_ℓ, 0 may occur at different frequencies for different ℓ, possibly enhancing the relative contribution of the missing higher order terms.

Appendix B: Double Compton Gaunt factor

As anticipated in Sect. 10.3, here we provide an improved approximation for the double Compton Gaunt factor (see also Ravenni & Chluba 2020; Chluba et al. 2020b for recent studies) in the elastic limit, suitable for small BE-like distortions. Writing the double Compton term according to Burigana et al. (1995), and accounting separately in the rate for the correction factor, C_mr, for mildly relativistic thermal plasma in the soft photon limit (Chluba et al. 2007), the double Compton Gaunt factor can be written as (Burigana et al. 1991a)

$$\begin{array}{c}\hfill g_{\mathrm{DC}}(x_{\mathrm{e}})=\frac{{\int }_{2x_{\mathrm{e}}}^{\infty }{x}_{\mathrm{e}}^{\prime 4}[1+\eta (x_{\mathrm{e}}^{\prime }-x_{\mathrm{e}})]\eta (x_{\mathrm{e}}^{\prime })\left[w\frac{F(w)}{2}\right]\mathrm{d}{x}_{\mathrm{e}}^{\prime }}{{\int }_0^{\infty }[1+\eta (x_{\mathrm{e}}^{\prime })]\eta (x_{\mathrm{e}}^{\prime })x_{\mathrm{e}}^{\prime 4}\mathrm{d}{x}_{\mathrm{e}}^{\prime }},\end{array}$$(B.1)

where

$$\begin{array}{cc}\hfill w\frac{F(w)}{2}=& \frac{1}{2}(1-w)\hfill \\ \hfill & \times \left[1+{(1-w)}^2+\frac{w^2(1+w^2)}{{(1-w)}^2}+w^4+w^2{(1-w)}^2\right],\hfill \end{array}$$(B.2)

with w the ratio between the frequencies of the created and incident photons, 0 ≤ w ≤ 1/2 and wF(w)/2 → 1 as w → 1 (Gould 1984). For a (pure) BE spectrum with a frequency independent chemical potential μ and equilibrium temperature $T_{\text{BE}}=T_{\text{eq}}=[h/(4k)]{\int }_0^{\infty }(1+\eta )\eta {\nu }^4\text{d}\nu /{\int }_0^{\infty }\eta {\nu }^3\text{d}\nu$ (Peyraud 1968; Zel’dovich & Levich 1969), the integral at denominator in Eq. (B.1) is simply 4I₃f(μ), with $f(\mu )={\int }_0^{\infty }{\eta }_{\text{BE}}(x_{\text{e}})x_{\text{e}}^3\text{d}{x}_{\text{e}}/I_3$, I₃ = π⁴/15 and f(μ)→1 in the Planckian limit μ → 0. Thus, for a pure BE spectrum only the integral at numerator of Eq. (B.1) needs a numerical computation. We have performed this calculation using both a Gaussian quadrature scheme (Press 1992) and the NAG routine D01AJF, adopting a very large upper integration limit (set to 500, to have a good estimation also at very high frequencies) and working with a dimensionless frequency in logarithmic space as integration variable. We first consider the case of a blackbody spectrum and compare the results of two methods (see Fig. B.1): clearly, they are in excellent agreement. We then consider the case of a (pure) BE spectrum with μ₀ = 1.4 × 10⁻⁵, the largest value considered in this work. The relative differences between g_DC computed for this case and in the case of a blackbody spectrum are less than few ×10⁻²% in the whole frequency range, and obviously it decreases for decreasing values of μ₀. Although the very low frequencies give only a little relative contribution to the integral, a similar comparison for a BE-like spectrum is in principle a bit more difficult, because the spectrum at low frequencies depends also on g_DC. On the other hand, the above difference represents a good upper limit estimation also when applied to a BE-like spectrum with the same high frequency asymptotic value of μ₀, because a BE-like spectrum differs from the blackbody less than a pure BE spectrum.

An approximation for g_DC(x_e) which, as explicitly stated by the authors, works well only at low frequencies, where the rate is higher, was found by Burigana et al. (1991a)

$$\begin{array}{c}\hfill {\stackrel{\sim }{g}}_{\mathrm{DC}}(x_{\mathrm{e}})\simeq e^{-x_{\mathrm{e}}/2}.\end{array}$$(B.3)

Of course, for very small distortions x_e ≃ x. In the limit μ ≪ 1, Chluba & Sunyaev (2012) found the approximation

$$\begin{array}{cc}\hfill {\widehat{g}}_{\mathrm{DC}}(x)\simeq & e^{-2x}\hfill \\ \hfill & \times [1+(3/2)x+(29/24)x^2+(11/16)x^3+(5/12)x^4].\hfill \end{array}$$(B.4)

Figure B.1 shows the relative accuracies of these two formulas. The latter works significantly better than the former at x ≳ 0.3, providing a little improvement at smaller frequencies (the two approximations clearly agree to first order in Taylor’s series for small x). On the other hand, an error at some % level still remains when using Eq. (B.4) at x >  1.

We then search for a better description of the numerical result. Considering the very good accuracy of Eq. (B.4) at low frequencies, we consider the following expression:

$$\begin{array}{c}\hfill g_{\mathrm{DC}}(x)\simeq {\widehat{g}}_{\mathrm{DC}}(x){10}^{P(X)}{e}^{-{(\widehat{x}/x)}^m}+{\widehat{g}}_{\mathrm{DC}}(x)[1-e^{-{(\widehat{x}/x)}^m}],\end{array}$$(B.5)

where X = log₁₀x, P(X) is a polynomial of a certain degree d, $\widehat{x}$ and m are a dimensionless frequency and an exponent. P(X) modifies the approximation represented by ${\widehat{g}}_{\mathrm{DC}}(x)$ to better describe the numerical results at high frequencies, while $\widehat{x}$ and m define the exponential weights to assure a smooth transition from low to high frequencies, around $\widehat{x}$. We fit our numerical results (for instance, the ones derived with the D01AJF routine) for g_DC(x) in the Planckian limit to find the values of $\widehat{x}$, m, d, varying them on a simple three dimensional grid, and of the best-fit polynomial coefficients. We find m = 1.6, $\widehat{x}=0.23$, d = 17 and the following polynomial coefficients, from the power of order 0 to the power of order d:

$$\begin{array}{cc}& +1.137720\times {10}^{-3},-2.949735\times {10}^{-2},-5.578415\times {10}^{-2},\hfill \\ \hfill & +3.212930\times {10}^{-2},+8.794246\times {10}^{-2},-4.492392\times {10}^{-3},\hfill \\ \hfill & -5.746682\times {10}^{-2},-1.277239\times {10}^{-2},+1.573671\times {10}^{-2},\hfill \\ \hfill & +7.068248\times {10}^{-3},-9.391122\times {10}^{-4},-1.073426\times {10}^{-3},\hfill \\ \hfill & -1.716636\times {10}^{-4},+8.114468\times {10}^{-6},-4.003431\times {10}^{-7},\hfill \\ \hfill & -1.886616\times {10}^{-6},-3.587674\times {10}^{-7},-2.055208\times {10}^{-8}.\hfill \end{array}$$

The accuracy of this approximation is also displayed in Fig. B.1: it is always better than ≃1% and better than ≃0.1−0.2% in the whole relevant frequency range.

Appendix C: Reference system rotation and dipole terms with m≠ 0

Let us consider a real field, T(θ, ϕ), defined by the spherical harmonics coefficients A_ℓ, m in a given reference system, S, with polar coordinates θ and ϕ, that is:

$$\begin{array}{c}\hfill T(\theta ,\varphi )={\displaystyle \sum _{\ell =0}^{{\ell }_{\max }}\sum _{m=-\ell }^{\ell }{A}_{\ell ,m}{Y}_{\ell ,m}(\theta ,\varphi ),}\end{array}$$(C.1)

and another reference system S′, with coordinates θ′ and ϕ′, that is defined, with respect to S, by the Euler angles α_E, β_E, γ_E. The range of α_E and γ_E is defined modulo 2π radians: we adopt the range [0, 2π] (but another widely adopted choice is, for example, [ − π, π]). The range of β_E covers π radians: we adopt the range [0, π] (but it could be, for example, [ − π/2, π/2]). In the reference system S′

$$\begin{array}{c}\hfill T({\theta }^{\prime },{\varphi }^{\prime })={\displaystyle \sum _{\ell =0}^{{\ell }_{\max }}\sum _{j=-\ell }^{\ell }{D}_{\ell ,j}{Y}_{\ell ,j}({\theta }^{\prime },{\varphi }^{\prime }),}\end{array}$$(C.2)

where D_ℓ, j are the corresponding (rotated) spherical harmonics coefficients. According to Goldstein (1984)

$$\begin{array}{c}\hfill D_{\ell ,j}={\displaystyle \sum _{m=-\ell }^{\ell }{Q}_{\ell ,m,j}{A}_{\ell ,m},}\end{array}$$(C.3)

where

$$\begin{array}{cc}\hfill Q_{\ell ,m,j}=& {(-1)}^{\ell -m}{\left(\begin{array}{c}2\ell \\ \ell +m\end{array}\right)}^{1/2}{\left(\begin{array}{c}2\ell \\ \ell +j\end{array}\right)}^{-1/2}{i}^{j-m}{e}^{i(m{\alpha }_{\mathrm{E}}+j{\gamma }_{\mathrm{E}})}\hfill \\ \hfill & \times {\cos }^{m+j}({\beta }_{\mathrm{E}}/2){\sin }^{j-m}({\beta }_{\mathrm{E}}/2)J_{l-j}^{m+j,j-m}(-\cos {\beta }_{\mathrm{E}});\hfill \end{array}$$(C.4)

here i² = −1 and $J_r^{s,t}(z)$ are the Jacobi polynomials

$$\begin{array}{c}\hfill J_r^{s,t}(z)=\frac{1}{2^r}{\displaystyle \sum _{k=0}^r\left(\begin{array}{c}r+s\\ k\end{array}\right)\left(\begin{array}{c}r+t\\ r-k\end{array}\right){(z-1)}^{r-k}{(z+1)}^k}\end{array}$$(C.5)

expressed in explicit polynomial form. For r = 0, we simply have $J_0^{s,t}(z)=1$.

Defining A_1, 0 = a_1, 0, A_1, 1 = b_1, 1 + i c_1, 1 (and, for symmetry, A_1, −1 = −b_1, 1 + i c_1, 1), after some algebra one gets

$$\begin{array}{cc}& \mathrm{Re}(D_{1,1})=\frac{\sqrt{2}}{2}{a}_{1,0}\sin {\beta }_{\mathrm{E}}\sin {\gamma }_{\mathrm{E}}\hfill \\ \hfill & +\frac{1+\cos {\beta }_{\mathrm{E}}}{2}[b_{1,1}\cos ({\alpha }_{\mathrm{E}}+{\gamma }_{\mathrm{E}})-c_{1,1}\sin ({\alpha }_{\mathrm{E}}+{\gamma }_{\mathrm{E}})]\hfill \\ \hfill & +\frac{1-\cos {\beta }_{\mathrm{E}}}{2}[b_{1,1}\cos ({\alpha }_{\mathrm{E}}-{\gamma }_{\mathrm{E}})-c_{1,1}\sin ({\alpha }_{\mathrm{E}}-{\gamma }_{\mathrm{E}})],\hfill \end{array}$$(C.6)

$$\begin{array}{cc}& \mathrm{Im}(D_{1,1})=-\frac{\sqrt{2}}{2}{a}_{1,0}\sin {\beta }_{\mathrm{E}}\cos {\gamma }_{\mathrm{E}}\hfill \\ \hfill & +\frac{1+\cos {\beta }_{\mathrm{E}}}{2}[b_{1,1}\sin ({\alpha }_{\mathrm{E}}+{\gamma }_{\mathrm{E}})+c_{1,1}\cos ({\alpha }_{\mathrm{E}}+{\gamma }_{\mathrm{E}})]\hfill \\ \hfill & -\frac{1-\cos {\beta }_{\mathrm{E}}}{2}[b_{1,1}\sin ({\alpha }_{\mathrm{E}}-{\gamma }_{\mathrm{E}})+c_{1,1}\cos ({\alpha }_{\mathrm{E}}-{\gamma }_{\mathrm{E}})].\hfill \end{array}$$(C.7)

A rotation of the reference system S → S′, that is, a set of Euler angles α_E, β_E, γ_E, that makes D_1, 1 = 0 = D_1, −1 or, equivalently, for which the whole dipole signal is along the z′ axis of the S′ reference system, does not depend on the choice of the angle γ_E that specifies the (last) rotation around the z′ axis. Thus, setting for simplicity γ_E = 0, the condition D_1, 1 = 0 is satisfied for

$$\begin{array}{c}\hfill b_{1,1}\cos {\alpha }_{\mathrm{E}}-c_{1,1}\sin {\alpha }_{\mathrm{E}}=0\end{array}$$(C.8)

and

$$\begin{array}{c}\hfill -\frac{\sqrt{2}}{2}{a}_{1,0}\sin {\beta }_{\mathrm{E}}+\cos {\beta }_{\mathrm{E}}[c_{1,1}\cos {\alpha }_{\mathrm{E}}+b_{1,1}\sin {\alpha }_{\mathrm{E}}]=0.\end{array}$$(C.9)

Equations (C.8) and (C.9) imply

$$\begin{array}{cc}& \sin {\alpha }_{\mathrm{E}}=\frac{\sqrt{2}}{2}\tan {\beta }_{\mathrm{E}}\frac{a_{1,0}{b}_{1,1}}{b_{1,1}^2+c_{1,1}^2}=\pm \frac{b_{1,1}}{\sqrt{b_{1,1}^2+c_{1,1}^2}},\hfill \end{array}$$(C.10)

$$\begin{array}{cc}& \cos {\alpha }_{\mathrm{E}}=\frac{\sqrt{2}}{2}\tan {\beta }_{\mathrm{E}}\frac{a_{1,0}{c}_{1,1}}{b_{1,1}^2+c_{1,1}^2}=\pm \frac{c_{1,1}}{\sqrt{b_{1,1}^2+c_{1,1}^2}},\hfill \end{array}$$(C.11)

where the last equalities of Eqs. (C.10) and (C.11) come from the obvious condition sin² α_E + cos² α_E = 1 that gives

$$\begin{array}{c}\hfill \tan {\beta }_{\mathrm{E}}=\pm \frac{\sqrt{2(b_{1,1}^2+c_{1,1}^2)}}{a_{1,0}};\end{array}$$(C.12)

the sign is related to the orientation of the (starting) reference system, S, with respect to the dipole pattern.

Finally, simply considering that the angular power spectrum, C_ℓ, is invariant under rotation, for ℓ = 1, we have

$$\begin{array}{cc}\hfill \sqrt{C_1}& =\frac{\sqrt{D_{1,0}^2+D_{1,1}^2+D_{1,-1}^2}}{\sqrt{3}}=\frac{D_{1,0}}{\sqrt{3}}=\frac{\mathrm{Re}(D_{1,0})}{\sqrt{3}}\hfill \\ \hfill & =\frac{\sqrt{A_{1,0}^2+A_{1,1}^2+A_{1,-1}^2}}{\sqrt{3}}=\frac{\sqrt{a_{1,0}^2+2(b_{1,1}^2+c_{1,1}^2)}}{\sqrt{3}}·\hfill \end{array}$$(C.13)

Let us consider now, see Eq. (86), $A_{1,m}=a_{1,m}^{\mathrm{BB}}(\beta )+a_{\mathrm{BB},1,m}$, as in the case of frequency independent $a_{1,m}^{\mathrm{glob}}$ coefficients evaluated for a Planckian spectrum, but in a reference system S with the z axis only roughly aligned with the observer peculiar velocity β, that is easily identified from the sky area where the observed (large scale pattern) temperature is not far from its maximum (this implies A_1, 0 = a_1, 0 >  0 and, with the adopted ranges for the Euler angles, this also defines a unique, positive choice of the sign in Eq. (C.12) as well as in Eqs. (C.10) and (C.11), since β_E should be relatively small). Thus, the above expressions allow to simply find a reference system S′ where D_1, 1 = 0 = D_1, −1 or, equivalently, the global dipole signal is along the z′ axis for all the frequencies (and, in Eq. (C.13), the positive sign of D_1, 0 means that the z′ axis points towards the direction of the maximum of the dipole pattern signal).

Considering, instead, $A_{1,m}=a_{1,m}^{\mathrm{glob}}(\nu ,\beta )$ with $a_{1,m}^{\mathrm{glob}}(\nu ,\beta )$ given by Eq. (86) but in the presence of deviations from a Planckian spectrum. The coefficients a_1, 0, b_1, 1 and c_1, 1 in the above expressions are then frequency dependent, it is no longer possible to find a reference system S′ satisfying D_1, 1 = 0 = D_1, −1 for all the frequencies, and, in principle, this property can be used to constrain the intrinsic dipole, as discussed in Sect. 15.

All Figures