© L. Vacher et al. 2022

1. Introduction

Our present understanding of the primordial Universe relies on the paradigm of inflation (Brout et al. 1978; Starobinsky 1980; Guth 1981), introducing a phase of accelerated expansion in the first fractions of a second after the primordial singularity. Such a phenomenon is expected to leave a background of gravitational waves propagating in the primordial plasma during recombination, leaving a permanent mark imprinted in the polarization anisotropies of the cosmic microwave background (CMB): the primordial B-modes (Polnarev 1985; Kamionkowski et al. 1997; Seljak & Zaldarriaga 1997). The amplitude of the angular power spectrum of those primordial B-modes is characterized by the tensor-to-scalar ratio r, which is proportional to the energy scale at which inflation occurred (Lyth 1997). Hence, looking for this smoking gun of inflation allows us to test our best theories of fundamental physics in the primordial Universe at energy scales far beyond the reach of particle accelerators. In this scope, it is one of the biggest challenges of cosmology set out for the next decades. The best experimental upper limit on the r parameter so far is r < 0.032 (95% C.L., Tristram et al. 2021; Bicep/Keck Collaboration 2021; BICEP2/Keck & Planck Collaboration 2015).

The JAXA Lite (Light) satellite, used for the B-mode polarization and Inflation from cosmic background Radiation Detection (LiteBIRD) mission, is designed to observe the sky at large angular scales in order to constrain this parameter r down to δr = 10⁻³, including all sources of uncertainty (Hazumi 2018; LiteBIRD Collaboration 2020). Exploring this region of the parameter space is critical, because this order of magnitude for the tensor-to-scalar ratio is predicted by numerous physically motivated inflation models (for a review see e.g., Martin et al. 2014)

However, the success of this mission relies on our ability to treat polarized foreground signals. Indeed various diffuse astrophysical sources emit polarized B-mode signals above the primordial ones, the strongest being due to the diffuse polarized emission of our own Galaxy (Planck Collaboration IV 2020). Even in a diffuse region like the BICEP/Keck field, the Galactic B-modes are at least ten times stronger at 150 GHz than the r = 0.01 tensor B-modes targeted by the current CMB experiments (BICEP2 Collaboration & Keck Array Collaboration 2018).

The true complexity of polarized foreground emission that the next generation of CMB experiments will face is still mostly unknown today. Underestimation of this complexity can lead to the estimation of a spurious nonzero value of r (see e.g., Planck Collaboration Int. L 2017; Remazeilles et al. 2016).

At high frequencies (> 100 GHz), the thermal emission of interstellar dust grains is the main source of Galactic foreground contaminating the CMB (Krachmalnicoff et al. 2016; Planck Collaboration XI 2020). The canonical model of the spectral energy distribution (SED) of this thermal emission for intensity and polarization is given by the modified black body (MBB) law (Desert et al. 1990). This model provides a good fit to the dust polarization SED at the sensitivity of the Planck satellite (Planck Collaboration XI 2020) but it may not fully account for it at the sensitivity of future experiments (Hensley & Bull 2018). Furthermore, due to changes of physical conditions across the galaxy, spatial variations of the SEDs are present between and along the lines of sight. The former leads to what is known as frequency decorrelation in the CMB community (see e.g. Tassis & Pavlidou 2015; Planck Collaboration Int. L 2017; Pelgrims et al. 2021). Moreover, both effects lead to averaging MBBs when observing the sky (unavoidable line-of-sight or beam-integration effects). Because of the nonlinearity of the MBB law, those averaging effects will distort the SED, leading to deviations from this canonical model (Chluba et al. 2017).

Chluba et al. (2017) proposed a general framework called “moment expansion” of the SED to take into account those distortions, using a Taylor expansion around the MBB with respect to its spectral parameters (Taylor expansion of foreground SEDs was discussed in previous studies; see e.g., Stolyarov et al. 2005). This method is agnostic: it does not require any assumption on the real complexity of the polarized dust emission. The moment expansion approach thus provides a promising tool with which to model the unanticipatable complexity of the dust emission in real data.

Mangilli et al. (2021) generalized this formalism for the sake of CMB data analysis in harmonic space and for cross-angular power spectra and applied it successfully to complex simulations and Planck High-Frequency Instrument (HFI) intensity data. This latter work shows that the real complexity of Galactic foregrounds could be higher than expected, encouraging us to follow the path opened by the moment expansion formalism.

In the present work, we apply the moment expansion in harmonic space to characterize and treat the dust foreground polarized emission of LiteBIRD high-frequency simulations, using dust-emission models of increasing complexity. We discuss the ability of this method to recover an unbiased value for the r parameter, with enough accuracy to achievethe scientific objectives of the LiteBIRD mission.

In Sect. 2, we first review the formalism of moment expansion in map and harmonic domains. We then describe in Sect. 3 how we realize several sets of simulations of the sky as seen by the LiteBIRD instrument with varying dust complexity and how we estimate the angular power spectra. In Sect. 4, we describe how we estimate the moment parameters and the tensor-to-scalar ratio r in those simulations. The results are then presented in Sect. 5. Finally, we discuss those results and the future work that has to be done in the direction opened by moment expansion in Sect. 6.

2. Formalism

2.1. Characterizing the dust SED in real space

2.1.1. Modified black body model

The canonical way to characterize astrophysical dust-grain emission in every volume element of the Galaxy is given by the modified black body (MBB) function, consisting of multiplying a standard black body SED B_ν(T) at a given temperature T₀ by a power-law of the frequency ν with a spectral index β₀. The dust intensity map I_D(ν,n) observed at a frequency ν in every direction with respect to the unit vector n, can then be written as:

$$\begin{array}{c}\hfill I(\nu ,\mathit{n})={\left(\frac{\nu }{{\nu }_0}\right)}^{{\beta }_0}\frac{B_{\nu }(T_0)}{B_{{\nu }_0}(T_0)}A(\mathit{n})=\frac{I_{\nu }({\beta }_0,T_0)}{I_{{\nu }_0}({\beta }_0,T_0)}A(\mathit{n}),\end{array}$$(1)

where A(n) is the dust intensity template at a reference frequency ν₀¹. We know that the physical conditions (thermodynamic and dust grain properties) change through the interstellar medium across the Galaxy, depending, in an intricate fashion, on the gas velocity and density, the interstellar radiation field, the distance to the Galactic center (see e.g., Paradis et al. 2009; Ysard et al. 2015; Planck Collaboration XI 2014; Planck Collaboration IV 2020; Hutton et al. 2015; Fanciullo et al. 2015). This change of physical conditions leads to variations in β and T depending on the direction of observation n:

$$\begin{array}{c}\hfill I(\nu ,\mathit{n})=\frac{I_{\nu }(\beta (\mathit{n}),T(\mathit{n}))}{I_{{\nu }_0}(\beta (\mathit{n}),T(\mathit{n}))}A(\mathit{n}).\end{array}$$(2)

The SED amplitude and parameters (temperature and spectral index) are then different for every line of sight. It is therefore clear that, in order to provide a realistic model of the dust emission, the frequency and spatial dependencies may not be trivially separated.

2.1.2. Limits of the modified black body

The dust SED model given by the MBB has proven to be highly accurate (Planck Collaboration Int. XVII 2014; Planck Collaboration Int. XXII 2015). However, it must be kept in mind that this model is empirical and is therefore not expected to give a perfect description of the dust SED in the general case. Indeed, physically motivated dust grain emission models predict deviations from it (e.g., Draine & Hensley 2013). Surveys tend to show that the dust-emission properties vary across the observed 2D sky and the 3D Galaxy (Planck Collaboration XI 2020). Furthermore, in true experimental conditions, one can never directly access the pure SED of a single volume element with specific emission properties and unique spectral parameters. Averages are therefore made over different SEDs emitted from distinct regions with different physical emission properties, in a way that may not be avoided: along the line of sight; between different lines of sight, inside the beam of the instrument or; when doing a spherical harmonic decomposition to calculate the angular power spectra over large regions of the sky.

The MBB function is nonlinear, and therefore summing MBBs with different spectral parameters does not return another MBB function and produces SED distortions. For all these reasons, modeling the dust emission with a MBB is intrinsically limited, even when doing so with spatially varying spectral parameters. As a consequence, inaccuracies might appear when modeling the dust contribution to CMB data that will unavoidably impact the final estimation of the cosmological parameters.

2.1.3. Moment expansion in pixel space

A way to address the limitation of the MBB model in accurately describing the dust emission is given by the moment expansion formalism proposed by Chluba et al. (2017). This formalism is designed to take into account the SED distortions due to averaging effects by considering a multidimensional Taylor expansion of the distorted SED I(ν,p) around the mean values p₀ of its spectral parameters p = {p_i}. This is the so-called moment expansion of the SED, which can be written as

$$\begin{array}{cc}\hfill I(\nu ,\mathit{p})=I(\nu ,{\mathit{p}}_0)& +{\displaystyle \sum _i{\omega }_1^{p_i}{⟨{\partial }_{p_i}I(\nu ,\mathit{p})⟩}_{\mathit{p}={\mathit{p}}_0}}\hfill \\ \hfill & +\frac{1}{2}{\displaystyle \sum _{i,j}{\omega }_2^{p_i{p}_j}{⟨{\partial }_{p_i}{\partial }_{p_j}I(\nu ,\mathit{p})⟩}_{\mathit{p}={\mathit{p}}_0}}\hfill \\ \hfill & +\dots \hfill \\ \hfill & +\frac{1}{\alpha !}{\displaystyle \sum _{i,\dots ,k}{\omega }_{\alpha }^{p_i\dots p_k}{⟨{\partial }_{p_i}\dots {\partial }_{p_k}I(\nu ,\mathit{p})⟩}_{\mathit{p}={\mathit{p}}_0},}\hfill \end{array}$$(3)

where the first term on the right-hand side is the SED without distortion I(ν,p₀) evaluated at p = p₀, and the other terms are the so-called moments of order α, quantified by the moment parameters ${\omega }_{\alpha }^{p_i\dots p_k}$ for the expansion with respect to any parameter of p. Performing the expansion to increasing order adds increasing complexity to the SED I(ν,p₀).

For the MBB presented in Sect. 2.1.1, there are two parameters so that p = {β,T}. Thus the dust moment expansion reads

$$\begin{array}{cc}\hfill I(\nu ,\mathit{n})& =\frac{I_{\nu }({\beta }_0,T_0)}{I_{{\nu }_0}({\beta }_0,T_0)}\{A(\mathit{n})+{\omega }_1^{\beta }(\mathit{n})ln\left(\frac{\nu }{{\nu }_0}\right)+\frac{1}{2}{\omega }_2^{\beta }(\mathit{n}){ln}^2\left(\frac{\nu }{{\nu }_0}\right)\hfill \\ \hfill & +{\omega }_1^T(\mathit{n})(\mathrm{\Theta }(\nu ,T_0)-\mathrm{\Theta }({\nu }_0,T_0))+\dots \},\hfill \end{array}$$(4)

where the expansion has been written up to order two in β (with moment expansion parameters ${\omega }_1^{\beta }$ at order one and ${\omega }_2^{\beta }$ at order two) and to order one in T (with a moment expansion parameter ${\omega }_1^T$ at order one). The following expression has been introduced to simplify the black body derivative with respect to T:

$$\begin{array}{c}\hfill \mathrm{\Theta }(\nu ,T)=\frac{x}{T}\frac{e^x}{e^x-1},\mathrm{with}x=\frac{h\nu }{\mathit{kT}}·\end{array}$$(5)

The moment expansion in pixel space can be used for component separation and possibly crossed with other methods (see e.g., Remazeilles et al. 2021; Adak 2021). However, in the present work, we are interested in the modeling of the dust at the B-mode angular power spectrum level. Performing the moment expansion at the angular power spectrum level adds some complexity to the SEDs due to the additional averaging occurring when dealing with spherical harmonic coefficients. Indeed, these coefficients are estimated on potentially large fractions of the sky and probe regions with various physical conditions. On the other hand, the expansion at the power spectrum level possibly drastically reduces the parameter space with respect to performing the expansion in every sky pixel.

2.2. Characterizing the dust SED in harmonic space

2.2.1. Dust SED in spherical harmonic space

The expansion presented in Sect. 2.1.3 can be applied in spherical harmonic space using the same logic. The sky emission projection then reads

$$\begin{array}{c}\hfill I(\nu ,\mathit{n})={\displaystyle \sum _{\ell =0}^{\infty }\sum _{m=-\ell }^{\ell }{I}_{\ell m}^{\nu }{Y}_{\ell m}(\mathit{n}).}\end{array}$$(6)

Applying the moment expansion to the spherical harmonics coefficients, with respect to β and T, as in Eq. (4), leads to

$$\begin{array}{cc}\hfill I_{\ell m}^{\nu }=& \frac{I_{\nu }({\beta }_0(\ell ),T_0(\ell ))}{I_{{\nu }_0}({\beta }_0(\ell ),T_0(\ell ))}\{A_{\ell m}+{\omega }_{1,\ell m}^{\beta }ln\left(\frac{\nu }{{\nu }_0}\right)+\frac{1}{2}{\omega }_{2,\ell m}^{\beta }{ln}^2\left(\frac{\nu }{{\nu }_0}\right)\hfill \\ \hfill & +{\omega }_{1,\ell m}^T(\mathrm{\Theta }(\nu ,T_0(\ell ))-\mathrm{\Theta }({\nu }_0,T_0(\ell )))+\dots \},\hfill \end{array}$$(7)

where this time β₀(ℓ) and T₀(ℓ) are the averages of β and T at a given multipole ℓ over the sky fraction we are looking at. We note that the moment parameters ${\omega }_{\alpha ,\ell m}^{p_i}$ involved here are different from the ${\omega }_i^{p_i}(\mathit{n})$ appearing in Eq. (4) in the map space because they involve different averaging. In principle, the moment expansion in harmonic space can take into account the three kinds of spatial averages presented in Sect. 2.1.2.

As the dust spectral index and temperature are difficult to separate in the frequency range considered for CMB studies (i.e., Rayleigh-Jeans domain, see e.g. Juvela & Ysard 2012), the moment expansion in harmonic space has only been applied in the past with respect to β, with the temperature being fixed to a reference value T = T₀ (Mangilli et al. 2021; Azzoni et al. 2021). In the present paper, for the first time, the moment expansion in harmonic space is instead performed with respect to both β and T, as it was in real space in Remazeilles et al. (2021).

2.2.2. Cross-power spectra

Relying on the derivation made by Mangilli et al. (2021) and Eq. (7), we can explicitly write the cross-spectra between two maps M_νi and M_νj at frequencies ν_i and ν_j, using the moment expansion in β and T as follows:

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\ell }({\nu }_i\times {\nu }_j)& =\frac{I_{{\nu }_i}({\beta }_0(\ell ),T_0(\ell ))I_{{\nu }_j}({\beta }_0(\ell ),T_0(\ell ))}{I_{{\nu }_0}{({\beta }_0(\ell ),T_0(\ell ))}^2}·\{\hfill \\ \hfill 0^{\mathrm{th}}\mathrm{order}& \{\begin{array}{cc}& {\mathcal{D}}_{\ell }^{A\times A}\hfill \end{array}\hfill \\ \hfill 1^{\mathrm{st}}\mathrm{order}\beta & \{\begin{array}{cc}& +{\mathcal{D}}_{\ell }^{A\times {\omega }_1^{\beta }}[ln\left(\frac{{\nu }_i}{{\nu }_0}\right)+ln\left(\frac{{\nu }_j}{{\nu }_0}\right)]\hfill \\ \hfill & +{\mathcal{D}}_{\ell }^{{\omega }_1^{\beta }\times {\omega }_1^{\beta }}[ln\left(\frac{{\nu }_i}{{\nu }_0}\right)ln\left(\frac{{\nu }_j}{{\nu }_0}\right)]\hfill \end{array}\hfill \\ \hfill 1^{\mathrm{st}}\mathrm{order}T& \{\begin{array}{cc}& +{\mathcal{D}}_{\ell }^{A\times {\omega }_1^T}({\mathrm{\Theta }}_i+{\mathrm{\Theta }}_j-2{\mathrm{\Theta }}_0)\hfill \\ \hfill & +{\mathcal{D}}_{\ell }^{{\omega }_1^T\times {\omega }_1^T}({\mathrm{\Theta }}_i-{\mathrm{\Theta }}_0)({\mathrm{\Theta }}_j-{\mathrm{\Theta }}_0)\hfill \end{array}\hfill \\ \hfill 1^{\mathrm{st}}\mathrm{order}T\beta & \{\begin{array}{cc}& +{\mathcal{D}}_{\ell }^{{\omega }_1^{\beta }\times {\omega }_1^T}[ln\left(\frac{{\nu }_j}{{\nu }_0}\right)({\mathrm{\Theta }}_i-{\mathrm{\Theta }}_0)+ln\left(\frac{{\nu }_i}{{\nu }_0}\right)({\mathrm{\Theta }}_j-{\mathrm{\Theta }}_0)]\hfill \end{array}\hfill \\ \hfill 2^{\mathrm{nd}}\mathrm{order}\beta & \{\begin{array}{cc}& +\frac{1}{2}{\mathcal{D}}_{\ell }^{A\times {\omega }_2^{\beta }}[{ln}^2\left(\frac{{\nu }_i}{{\nu }_0}\right)+{ln}^2\left(\frac{{\nu }_j}{{\nu }_0}\right)]\hfill \\ \hfill & +\frac{1}{2}{\mathcal{D}}_{\ell }^{{\omega }_1^{\beta }\times {\omega }_2^{\beta }}[ln\left(\frac{{\nu }_i}{{\nu }_0}\right){ln}^2\left(\frac{{\nu }_j}{{\nu }_0}\right)+ln\left(\frac{{\nu }_j}{{\nu }_0}\right){ln}^2\left(\frac{{\nu }_i}{{\nu }_0}\right)]\hfill \\ \hfill & +\frac{1}{4}{\mathcal{D}}_{\ell }^{{\omega }_2^{\beta }\times {\omega }_2^{\beta }}\left[{ln}^2\left(\frac{{\nu }_i}{{\nu }_0}\right){ln}^2\left(\frac{{\nu }_j}{{\nu }_0}\right)\right]\hfill \end{array}\hfill \\ \hfill & +\dots \},\hfill \end{array}$$(8)

where we use the following abbreviation: Θ(ν_k,T₀(ℓ) ≡ Θk, so that Θ₀ = Θ(ν₀,T₀(ℓ)), and we defined the moment expansion cross-power spectra between two moments ℳ and 𝒩 as

$$\begin{array}{c}\hfill {\mathcal{C}}_{\ell }^{\mathcal{M}\times \mathcal{N}}={\displaystyle \sum _{m,m^{\prime }=-\ell }^{\ell }{\mathcal{M}}_{\ell m}{\mathcal{N}}_{\ell m^{\prime }},\mathrm{with}(\mathcal{M},\mathcal{N})\in \{A,{\omega }_1^{\beta },{\omega }_1^T,{\omega }_2^{\beta },\dots \}.}\end{array}$$(9)

In the remainder of this article, we use the 𝒟_ℓ quantity, which is a scaling of the angular power spectra, and is defined as

$$\begin{array}{c}\hfill {\mathcal{D}}_{\ell }\equiv \frac{\ell (\ell +1)}{2\pi }{\mathcal{C}}_{\ell }.\end{array}$$(10)

Equation (8) has been written using the expansion with respect to β at order two and T at order one, as in Eq. (7). Nevertheless, the terms involving power spectra between order two in β and order one in T have been neglected so as to match the needs of the implementation of our method in the following.

Hereafter, when we refer to “order k” at the angular power spectrum level, we are referring to moment expansion terms involving the pixel space moment up to order k. For example, $D_{\ell }^{A\times {\omega }_1^T}$ and $D_{\ell }^{{\omega }_1^{\beta }\times {\omega }_1^T}$ are order one, while $D_{\ell }^{A\times {\omega }_2^{\beta }}$, $D_{\ell }^{{\omega }_1^{\beta }\times {\omega }_2^{\beta }}$ and $D_{\ell }^{{\omega }_2^{\beta }\times {\omega }_2^{\beta }}$ are order two. At order zero, one retrieves the MBB description of the cross-angular power spectra SED 𝒟_ℓ(ν_i × ν_j) as a function of the frequencies ν_i and ν_j.

This formalism was originally introduced to analyze the complexity of intensity data in Mangilli et al. (2021). In the present work, we focus on B-mode polarization power spectra. This was put forward after analyzing the Planck and balloon-borne Large Aperture Submillimeter Telescope for Polarimetry (BLASTPol) data and finding that the polarization fraction appears to be constant in the far-infrared-to-millimetre wavelengths (Gandilo et al. 2016; Ashton et al. 2018). This allows us to assume that the same grain population is responsible for the total and polarized foreground emission (Guillet et al. 2018). As a result, intensity and polarization SED complexity may be similar. Nevertheless, Q and U can have a different SED because of the polarization angle frequency dependence (see e.g., Tassis & Pavlidou 2015; Ichiki et al. 2019) and so can E and B. This could be a limitation when analyzing the dust E and B with a single moment expansion, especially when SED variations occur along the line of sight. Even when trying to model a single polarization component –as we do in the present work, dealing only with B modes– it is not clear whether the distorted SED can be modeled in terms of β and T moments only. Further work needs to be done to assess this question. However, they should not impact the present study in which variations along the line of sight are not simulated.

Modeling the complexity of the foreground signals by means of the moment expansion of the B-mode angular power spectrum has already been successfully applied to Simons Observatory (The Simons Observatory collaboration 2019) simulated data (Azzoni et al. 2021). However, the approach taken by these latter authors is different from the one presented above. They apply a minimal moment expansion: assumptions are made to keep only the ${\mathcal{D}}_{\ell }^{{\omega }_1^{\beta }\times {\omega }_1^{\beta }}$ and ${\mathcal{D}}_{\ell }^{A\times {\omega }_2^{\beta }}$ parameters, which are modeled with a power-law scale dependence. These assumptions may not hold for experiments with higher sensitivity and observing wider sky patches. Furthermore, they assume a scale-invariant dust spectral index. In this work, on the other hand, we relax these assumptions in order to characterize the required spectral complexity of the dust emission for LiteBIRD.

3. Simulations and cross-spectra estimation

3.1. LiteBIRD

LiteBIRD is an international project proposed by the Japanese spatial agency (JAXA), which selected it in May 2019 as a strategic large class mission. The launch is planned for 2029 for a minimal mission duration of 3 years (Hazumi et al. 2020; LiteBIRD Collaboration 2021).

LiteBIRD is designed to realize a full sky survey of the CMB at large angular scales in order to look for the reionization bump of primordial B-modes and explore the scalar-to-tensor ratio (r) parameter space with a total uncertainty δ_r below 10⁻³, including foreground cleaning and systematic errors. LiteBIRD is composed of three telescopes observing in different frequency intervals: the Low-, Medium- and High-Frequency Telescopes (LFT, MFT and HFT). Each of the telescopes illuminates a focal plane composed of hundreds of polarimetric detectors. The whole instrument will be cooled down to 5 K (LiteBIRD Collaboration 2020) while the focal plane will be cooled down to 100 mK (Suzuki et al. 2018). In order to mitigate the instrumental systematic effects, the polarization is modulated by a continuously rotating half-wave plate. LiteBIRD will observe the sky in 15 frequency bands from 40 to 402 GHz. Table 1 gives the details of the frequency bands and their sensitivities in polarization (adapted from Hazumi et al. 2020, see Sect. 3.2.3).

3.2. Components of the simulations

We build several sets of LiteBIRD sky simulations. These multi-frequency sets of polarized sky maps are a mixture of CMB, dust, and instrumental noise. The simulations are made at the nine highest frequencies accessible by the instrument (≥ 100 GHz), where dust is the predominant source of foreground contamination. For every studied scenario, we built N_sim = 500 simulations, each composed of a set of N_freq = 9 pairs of sky maps (Q,U) built using the HEALPIX package, with N_side = 256 (Górski et al. 2005). All the signals will be expressed in μK_CMB units.

3.2.1. Cosmic microwave background signal

To generate the CMB signal, we use the Code for Anisotropies in the Microwave Background (CAMB, Lewis et al. 2000) to create a fiducial angular power spectrum from the best-fit values of cosmological parameters estimated by the recent Planck data analysis (Planck Collaboration I 2020).

For the B-modes, we consider the two different components of the spectrum: lensing-induced and primordial (tensor), so that ${\mathcal{D}}_{\ell }^{\mathit{BB}}={\mathcal{D}}_{\ell }^{\mathrm{lensing}}+r_{\mathrm{sim}}·{\mathcal{D}}_{\ell }^{\mathrm{tensor}}$, where ${\mathcal{D}}_{\ell }^{\mathrm{tensor}}$ refers to the tensor B-modes for r = 1 and r_sim labels the input values of the tensor-to-scalar ratio r contained in the simulation. We use two different values throughout this work: r_sim = 0, which is used in the present work as the reference simulations and r_sim = 10⁻² used for consistency checks when the CMB primordial signal is present.

For all simulations, we then generate the Stokes Q and U CMB polarization Gaussian realization maps $S_{\nu ,r_{\mathrm{sim}}}^{\mathrm{CMB}}$ from the angular power spectra using the synfast function of HEALPIX.

3.2.2. Foregrounds: dust

Our study focuses on high frequencies (≥ 100 GHz) only, where thermal dust emission is the main source of polarized foreground as mentioned in Sect. 1. We make use of two different scenarios of increasing complexity included in the PYSM (Thorne et al. 2017) and one of intermediate complexity not included in the PYSM:

• d0, included in the PYSM: the dust polarization Q and U maps are taken from $S_{\nu =353}^{\mathit{Planck}}$, the Planck 2015 data at 353 GHz (Planck Collaboration I 2016), extrapolated to a frequency ν using the MBB given in Eq. (1) with a temperature T₀ = T_d0 = 20 K and spectral index β₀ = β_d0 = 1.54 constant over the sky:

  $$\begin{array}{c}\hfill S_{\nu }^{\mathrm{dust}}=S_{\nu }^{\mathtt{d}}\mathtt{0}=\frac{I_{\nu }({\beta }_{\mathtt{d}}\mathtt{0},T_{\mathtt{d}}\mathtt{0})}{I_{{\nu }_0}({\beta }_{\mathtt{d}}\mathtt{0},T_{\mathtt{d}}\mathtt{0})}·S_{353}^{\mathit{Planck}},\end{array}$$(11)

• d1T, introduced here: the dust polarization Q and U maps are also taken from Planck Collaboration I (2016) but they are extrapolated to a frequency ν using the MBB given in Eq. (2), with spatially varying spectral index β(n), as in d1 and a fixed temperature T₀ = T_d1T = 21.9 K, obtained as the mean of the Planck COMMANDER dust temperature map (Planck Collaboration X 2016) on our f_sky = 0.7 sky mask:

  $$\begin{array}{c}\hfill S_{\nu }^{\mathrm{dust}}=S_{\nu }^{\mathtt{d}}\mathtt{1}\mathtt{T}=\frac{I_{\nu }(\beta (\mathit{n}),T_{\mathtt{d}}\mathtt{1}\mathtt{T})}{I_{{\nu }_0}(\beta (\mathit{n}),T_{\mathtt{d}}\mathtt{1}\mathtt{T})}·S_{353}^{\mathit{Planck}},\end{array}$$(12)

• d1, included in the PySM: similar to d1T with both a spatially varying temperature T(n) and spectral index β(n) obtained from the Planck data using the COMMANDER code (Planck Collaboration X 2016):

  $$\begin{array}{c}\hfill S_{\nu }^{\mathrm{dust}}=S_{\nu }^{\mathtt{d}}\mathtt{1}=\frac{I_{\nu }(\beta (\mathit{n}),T(\mathit{n}))}{I_{{\nu }_0}(\beta (\mathit{n}),T(\mathit{n}))}·S_{353}^{\mathit{Planck}}.\end{array}$$(13)

3.2.3. Instrumental noise

The band polarization sensitivities ${\sigma }_{Q,U}^{\mathrm{noise}}(\nu )$ are derived from the noise equivalent temperature (NET) values converted into μK arcmin for each telescope (LFT, MFT and HFT). As seen in Table 1, some frequency bands are overlapping between two telescopes. In this situation, we take the mean value of the two NETs, weighted by the beam full width at half maximum (FWHM) θ as:

$$\begin{array}{c}\hfill {\sigma }_{Q,U}^{\mathrm{noise}}({\nu }_{\mathrm{overlapping}})=\sqrt{\frac{1}{{\left(\frac{{\theta }_{\min }}{{\theta }_{\max }}{\sigma }_{Q,U}^{\mathrm{noise}}({\nu }_{{\theta }_{\min }})\right)}^{-2}+{\left({\sigma }_{Q,U}^{\mathrm{noise}}({\nu }_{{\theta }_{\max }})\right)}^{-2}}},\end{array}$$(14)

where θ_min is the smallest FWHM among the two and θ_max the largest. The band polarization sensitivities are displayed in Table 1. For every simulation, the noise component N_ν is generated in every pixel of the maps with a Gaussian distribution centered on zero, with standard deviation ${\sigma }_{Q,U}^{\mathrm{noise}}(\nu )$ weighted by the pixel size (and $\sqrt{2}·{\sigma }_{Q,U}^{\mathrm{noise}}(\nu )$ for the maps used to compute the auto-power spectra, see Sect. 3.4.2).

For simplicity, we choose to ignore beam effects in our simulations, assuming they can be taken into account perfectly. Simulations are thus produced at infinite (0 arcmin) resolution and no beam effect is corrected for when estimating the angular power spectrum. This is equivalent to convolving the maps by Gaussian beams of finite resolution and correcting the power spectra for the associated Gaussian beam window functions.

3.3. Combining signals and building the simulated maps

The simulated (Q,U) maps M_ν, for a given simulation, can be expressed as the sum:

$$\begin{array}{c}\hfill M_{\nu }=S_{\nu ,r_{\mathrm{sim}}}^{\mathrm{CMB}}+S_{\nu }^{\mathrm{dust}}+N_{\nu }.\end{array}$$(15)

Cosmic microwave background and noise are simulated stochastically: for each simulation, we generate a new realization of the CMB maps $S_{\nu ,r_{\mathrm{sim}}}^{\mathrm{CMB}}$ and the noise maps N_ν. The dust map $S_{\nu }^{\mathrm{dust}}$ is the same for each simulation, at a given frequency.

Hereafter, we use the notation d0, d1T, and d1 to refer to simulations containing only dust and LiteBIRD noise, d0c, d1Tc, and d1c for simulations including CMB, dust, and LiteBIRD noise and, finally, and c for the simulation containing only CMB and LiteBIRD noise. The different components present in these different simulation types are summarized in Table 2.

3.4. Angular power spectra of the simulations

3.4.1. Mask

A mask is applied on the simulated maps presented in Sect. 3.3 in order to exclude the Galactic plane from the power-spectrum estimation. The mask is created by setting a threshold on the polarized intensity ($P=\sqrt{Q^2+U^2}$) of the Planck 353 GHz map (Planck Collaboration I 2020)², smoothed with a 10° beam. In order to keep f_sky = 0.7, f_sky = 0.6, and f_sky = 0.5, the cut is applied at 121 μK, 80 μK, and 53 μK, respectively. We then realize a C2 apodization of the binary mask with a scale of 5° using NAMASTER (Alonso et al. 2019). The resulting Galactic masks are displayed in Fig. A.1. These masks are similar to those used in Planck Collaboration XI (2020).

3.4.2. Estimation of the angular power spectra

We use the NAMASTER³ software (Alonso et al. 2019) to compute the angular power spectra of each simulation. NAMASTER allows us to correct for the E to B leakage bias due to the incomplete sky coverage. Therein we use a purification process to suppress the effect of the E to B leakage in the variance. For every simulation, from the set of maps M_νi, we compute all the possible auto-frequency and cross-frequency spectra ${\mathcal{D}}_{\ell }({\nu }_i\times {\nu }_j)\equiv {\mathcal{D}}_{\ell }(M_{{\nu }_i}\times M_{{\nu }_j})$ with

$$\begin{array}{cc}\hfill {\nu }_i\times {\nu }_j\in \{& 100\times 100,100\times 119,100\times 140,\dots ,100\times 402,\hfill \\ \hfill & 119\times 140,\dots ,119\times 402,\hfill \\ \hfill & ⋮\hfill \\ \hfill & 337\times 337,337\times 402,\hfill \\ \hfill & 402\times 402\},\hfill \end{array}$$(16)

leading to N_cross = N_freq · (N_freq + 1)/2 = 45 cross-frequency spectra. These spectra are displayed in Fig. 1 for the case of the d1c simulation type.

Fig. 1. Mean value over the Nsim simulations of the B-mode angular power spectra 𝒟ℓ(νi × νj) for the d1c simulation type, with rsim = 0. The color bar spans all the Ncross spectra 𝒟ℓ(νi × νj), associated to their reduced cross-frequency ${\nu }_{\mathrm{red}.}=\sqrt{{\nu }_i{\nu }_j}$, from 100 GHz (dark red) to 402 GHz (dark blue). The input CMB lensing power spectrum is shown as a black dashed line.

In order to avoid noise auto-correlation in the auto-spectra (i.e., 𝒟_ℓ(ν_i × ν_j) when i = j), the latter are estimated in a way that differs slightly from what is presented in Sect. 3.2.3. We simulate two noise-independent data subsets at an observing frequency ν_i, with a noise amplitude $\sqrt{2}$ higher than that of the frequency band, and compute the cross-angular power spectrum between those. Thus, 𝒟_ℓ(ν_i × ν_i) is free from noise auto-correlation bias at the expense of multiplying the noise amplitude in the spectrum by a factor of two. This approach is similar to that commonly used by the Planck Collaboration (see e.g., Planck Collaboration Int. XXX 2016; Planck Collaboration XI 2020; Tristram et al. 2021).

The spectra are evaluated in the multipole interval ℓ ∈ [1,200] in order to be able to focus on the reionization and recombination bumps of the primordial B-modes spectra. The spectra are binned in N_ℓ = 20 bins of size Δℓ = 10 using NAMASTER. The same binning is applied throughout this article such that, in the following, the multipole ℓ denotes the multipole bin of size Δℓ = 10 centered on ℓ⁴.

From the sets of (Q,U) maps, NAMASTER computes the ${\mathcal{D}}_{\ell }^{\mathit{EE}}$, ${\mathcal{D}}_{\ell }^{\mathit{BB}}$, and ${\mathcal{D}}_{\ell }^{\mathit{EB}}$ angular power spectra; for the sake of the present analysis, we keep only ${\mathcal{D}}_{\ell }^{\mathit{BB}}$. Hence, when we discuss or analyze power spectra, we are referring to the B-mode power spectra ${\mathcal{D}}_{\ell }^{\mathit{BB}}$. All spectra are expressed in (μK_CMB)².

4. Best-fit implementation

In order to characterize the complexity of the dust SED that will be measured by LiteBIRD, we modeled the angular power spectra of our simulations described in Sect. 3 over the whole frequency and multipole ranges with the moment expansion formalism introduced in Sect. 2.

4.1. General implementation

For each multipole ℓ, we ordered the angular power spectra ${\mathcal{D}}_{\ell }^{\mathit{BB}}({\nu }_i\times {\nu }_j$) as in Eq. (16) in order to build a SED that is a function of both ν_i and ν_j. We fit this SED with models, as in Eq. (8) for example, using a Levenberg-Marquardt χ² minimization with mpfit (Markwardt 2009)⁵. All the fits performed with mpfit were also realized with more computationally heavy Monte Carlo Markov chains (MCMC) with emcee (Foreman-Mackey et al. 2013), giving compatible results, well within the error bars.

The reduced χ² minimization is given by

$$\begin{array}{c}\hfill {\chi }^2=\frac{1}{N_{\mathrm{d}.\mathrm{o}.\mathrm{f}.}}{\mathit{R}}^T{\mathbb{C}}^{-1}\mathit{R},\end{array}$$(17)

where N_d.o.f. is the number of degrees of freedom and ℂ is the covariance matrix of our N_sim simulations, represented in Fig. 2, of dimension (N_ℓ · N_cross)²:

$$\begin{array}{c}\hfill {\mathbb{C}}_{\ell ,{\ell }^{\prime }}^{i\times j,k\times l}=\mathrm{cov}({\mathcal{D}}_{\ell }^{\mathrm{sim}}({\nu }_i\times {\nu }_j),{\mathcal{D}}_{{\ell }^{\prime }}^{\mathrm{sim}}({\nu }_k\times {\nu }_l)).\end{array}$$(18)

Fig. 2. Correlation matrix (${\mathrm{Corr}}_{\ell {\ell }^{\prime }}\equiv {\mathbb{C}}_{\ell {\ell }^{\prime }}/\sqrt{{\mathbb{C}}_{\ell \ell }{\mathbb{C}}_{{\ell }^{\prime }{\ell }^{\prime }}}$) for the Nsim simulations in d1c. Every block represents a value of ℓ and contains the ordered Ncross = 45 cross-spectra. The red squares represent the truncation of the full covariance matrix applied in the analysis (kept entries in red, other entries set to zero).

The entire covariance matrix ℂ is, in general, not invertible. To avoid this, we kept only the ℓ = ℓ′ block-diagonal of ℂ with the strongest correlation values⁶, as well as the (ℓ = 6.5,ℓ′ = 16.5) off-diagonal blocks showing a significant anti-correlation, as illustrated in Fig. 2. It was then possible to invert the thus-defined truncated correlation matrix with the required precision most of the time.

In the case of the d1 simulation type, we experienced a fit convergence issue for ∼ 20% of the simulations, leading to a very large χ². In order to overcome this problem, two options lead to identical results: throwing away the outliers from the analysis or fitting using only the block-diagonal matrix (i.e., the ℓ = 6.5, ℓ′ = 16.5 block is set to zero). This last option solves the conversion issue while providing sufficient precision. The results presented in the following are using the block-diagonal matrix when the simulation type is d1.

Finally, in Eq. (17), R is the residual vector associated with every simulation of size N_ℓ × N_cross:

$$\begin{array}{c}\hfill \mathit{R}=\left(\begin{array}{c}{\mathcal{R}}_{\ell =6.5}(100\times 100)\\ {\mathcal{R}}_{\ell =6.5}(100\times 119)\\ ⋮\\ {\mathcal{R}}_{\ell =16.5}(100\times 100)\\ ⋮\\ {\mathcal{R}}_{\ell =196.5}(402\times 402)\\ \end{array}\right),\end{array}$$(19)

with ${\mathcal{R}}_{\ell }({\nu }_i\times {\nu }_j)={\mathcal{D}}_{\ell }^{\mathrm{sim}}({\nu }_i\times {\nu }_j)-{\mathcal{D}}_{\ell }^{\mathrm{model}}({\nu }_i\times {\nu }_j)$.

The expression used for the model to fit is given by:

$$\begin{array}{cc}\hfill {\mathcal{D}}_{\ell }^{\mathrm{model}}({\nu }_i\times {\nu }_j)=& {\mathcal{D}}_{\ell }^{\mathrm{dust}}({\beta }_0(\ell ),T_0(\ell ),{\mathcal{D}}_{\ell }^{\mathcal{M}\times \mathcal{N}}({\nu }_i\times {\nu }_j))\hfill \\ \hfill & +A_{\mathrm{lens}}·{\mathcal{D}}_{\ell }^{\mathrm{lensing}}+r·{\mathcal{D}}_{\ell }^{\mathrm{tensor}},\hfill \end{array}$$(20)

where A_lens is not a free parameter and will remain fixed to zero (when there is no CMB, simulation types d0, d1T, and d1) or one (when the CMB is included, simulation types d0c, d1Tc and d1c). We leave the question of the impact of dust modeling with moments on the lensing measurement for future work. In Eq. (20), the free parameters can thus be β₀(ℓ), T₀(ℓ), and ${\mathcal{D}}_{\ell }^{\mathcal{M}\times \mathcal{N}}({\nu }_i\times {\nu }_j)$ and the tensor-to-scalar ratio r. The estimated value of r is referred to as $\widehat{r}$

No priors on the parameters are used in order to explore the parameter space with minimal assumptions. Finally, a frequency-dependent conversion factor is included in ${\mathcal{D}}_{\ell }^{\mathrm{dust}}$ – from (MJy sr⁻¹)² to (μK_CMB)² – to express the dust spectra in (μK_CMB)² units. In those units, ${\mathcal{D}}_{\ell }^{\mathrm{lensing}}$ and ${\mathcal{D}}_{\ell }^{\mathrm{tensor}}$ are frequency-independent.

To mitigate the impact of outliers in our simulations, all the final values of the best-fit parameters and χ² distributions are represented by their median and median absolute deviations over N_sim values. For the tensor-to-scalar ratio $\widehat{r}$, we chose to represent all the best-fit values from the N_sim simulations in a histogram and we assume its distribution is normal. Fitting a Gaussian curve on this histogram and getting the mean and standard deviation gives us the final values of $\widehat{r}$ and ${\sigma }_{\widehat{r}}$ presented in the paper.

4.2. Implementation for the dust component

For the dust component, we consider four different fitting schemes, corresponding to four expressions for the dust model ${\mathcal{D}}_{\ell }^{\mathrm{dust}}$ in Eq. (20), which are referred to as “MBB”, “β-1”, “β-T”, and “β-2”. Each of them corresponds to a truncation of Eq. (8), keeping only some selected terms of the moment expansion: MBB stands for those of the modified black body, β-1 for those of the expansion in β at first order, β-2 for the expansion in β at second order, and β-T for the expansion in both β and T at first order. We chose the β-1 and β-2 truncations based on the studies of Mangilli et al. (2021) and Azzoni et al. (2021), where the dust SED moment expansion is performed only with respect to β. The β-T fitting scheme is instead the first-order truncation in both β and T, introduced here for the first time at the power spectrum level. The parameters fitted in each of these fitting schemes are summarized in Table 3. We note that the β-2 and β-T fitting schemes share the same number of free parameters. Finally, when we fit $\widehat{r}$ at the same time as the dust parameters, the fitting schemes will be referred to as rMBB, rβ-1, rβ-T, and rβ-2.

Different physical processes are expected to occur at different angular scales, leading to different SED properties. Thus, we estimate the dust-related parameters with one parameter per multipole bin. As an example, we estimate β₀ = β₀(ℓ) and T₀ = T₀(ℓ) to be able to take into account their scale dependence, at the cost of increasing the number of free parameters in our model. This is also true for the higher order moments. On the other hand, $\widehat{r}$ is not scale dependent and, when it is fitted, we add one single parameter over the whole multipole range.

In Mangilli et al. (2021), the first-order moment expansion parameter ${\mathcal{D}}_{\ell }^{A\times {\omega }_1^{\beta }}$ is considered to be the leading order correction to the MBB spectral index. We applied a similar approach in the present work, extending it to the dust temperature when it is fitted. In our pipeline, we proceed iteratively:

1. (i) we fit β₀(ℓ) and T₀(ℓ) at order zero (MBB), for each ℓ,

2. (ii) we fix β₀(ℓ) and T₀(ℓ) and fit the higher order parameters, as in Eq. (8), (iii) we update the β₀(ℓ) to β_corr(ℓ) (and T₀(ℓ) to T_corr(ℓ) in the case of β-T) as:

   $$\begin{array}{c}\hfill {\beta }_{\mathrm{corr}}(\ell )={\beta }_0(\ell )+\frac{{\mathcal{D}}_{\ell }^{A\times {\omega }_1^{\beta }}}{{\mathcal{D}}_{\ell }^{A\times A}},T_{\mathrm{corr}}(\ell )=T_0(\ell )+\frac{{\mathcal{D}}_{\ell }^{A\times {\omega }_1^T}}{{\mathcal{D}}_{\ell }^{A\times A}},\end{array}$$(21)

• (iv) and we iterate from (ii) fixing β₀(ℓ) = β_corr(ℓ), until ${\mathcal{D}}_{\ell }^{A\times {\omega }_1^{\beta }}$ converges to be compatible with zero (and T₀(ℓ) = T_corr(ℓ), until ${\mathcal{D}}_{\ell }^{A\times {\omega }_1^T}$ converges to zero in the case of β-T).

We used three such iterations, which we found to be sufficient to guarantee the convergence. As the moment expansion is a nonorthogonal and incomplete basis (Chluba et al. 2017), this iterative process is performed to ensure that the expansions up to different orders share the same β₀(ℓ) and T₀(ℓ) with ${\mathcal{D}}_{\ell }^{A\times {\omega }_1^{\beta }}=0$ and ${\mathcal{D}}_{\ell }^{A\times {\omega }_1^T}=0$.

5. Results

In this section, we present our evaluation of the best-fit parameters for the different fitting schemes presented in Sect. 4.1 on the B-mode cross-angular power spectra computed from the different simulation types presented in Sect. 3.3 and on the Galactic mask keeping f_sky = 0.7, which is defined in Sect. 3.4.1. We first tested the simulation types containing only dust and noise in order to calibrate the dust complexity of our data sets in Sect. 5.1. We then used CMB only plus noise simulations to assess the minimal error on $\widehat{r}$ in Sect. 5.2 and, finally, we explored the dust, CMB, and noise simulation types to assess the impact of the dust complexity on $\widehat{r}$ in Sect. 5.3.

5.1. Dust only

To evaluate the amplitude of the dust moment parameters contained in the dust simulations in the absence of CMB, we ran the fitting schemes presented in Sect. 4.1 in the three simulation types d0, d1T, and d1 presented in Sect. 3.3. In these cases, A_lens and r in Eq. (20) are both fixed to zero and the fitted parameters are given in Table 3 for every fitting scheme.

5.1.1. d0

The d0 dust maps presented in Sect. 3.2 extrapolate between frequency bands with a MBB SED with constant parameters over the sky: β_d0 = 1.54 and T_d0 = 20 K. We performed the fit with the four fitting schemes presented in Sect. 4.1.

In Fig. 3 the values of the reduced χ²(ℓ) for each fitting scheme are displayed. For every fitting scheme (MBB, β-1, β-T and β-2), the reduced χ² are close to 1 over the whole multipole range (slightly below 1 for the β-1, β-T and β-2 fitting scheme). This indicates that the MBB is a good fit to the cross-angular power spectra computed from the d0 maps with a spatially invariant MBB SED, as expected. Adding additional (higher order) parameters, such as with β-1, β-T and β-2, has no significant effect on the χ².

Fig. 3. Median of the reduced χ2 in every multipole bin ℓ, for all the Nsim simulations of d0 (top, orange), d1T (middle, green) and d1 (bottom, blue), on fsky = 0.7. The reduced χ2 values are reported for the four different fitting schemes: MBB (circles), β-1 (crosses), β-T (diamonds) and β-2 (triangles). The values for the four fitting schemes are shifted from each others by ℓ = 2, in order to distinguish them. The black dashed line represents ${\chi }_{\mathrm{red}}^2=1$.

In Fig. 4 we can see that the best-fit values of β₀(ℓ) and T₀(ℓ) are compatible with constant values β₀(ℓ) = β_d0 and T₀(ℓ) = T_d0, as expected for this simulated data set.

Fig. 4. Top: median of the best fit values of β0(ℓ) in d0 (orange), d1T (green), and d1 (blue) for the MBB (circles). βdo is marked by the dashed black line. Bottom: same as above but with T0(ℓ), the black dashed-lines being Td0 = 20 K and Td1T = 21.9 K.

The best-fit values of the dust amplitude and the moment-expansion parameters are presented in Figs. 5–8, respectively. The amplitude power spectrum is compatible with that of the dust template map used to build d0 and the moment-xpansion parameters are compatible with zero for every fitting scheme, as expected with no spatial variation of the SED. Therefore, the moment expansion method presented in Sect. 2 passes the null test in the absence of SED distortions, with the d0 simulated data set.

Fig. 5. Median of the best-fit values of ${\mathcal{D}}_{\ell }^{A\times A}$ for d0 (orange), d1T (green), and d1 (blue) using the MBB fitting scheme. The values for the three simulation types are shifted with respect to one another by ℓ = 2 in order to distinguish them. The black dashed line is the amplitude power spectrum of the dust template map used to build the three simulation sets d0, d1T, and d1.

Fig. 6. Best-fit values of the first-order moment ${\mathcal{D}}_{\ell }^{{\omega }_1^{\beta }\times {\omega }_1^{\beta }}$ for d0 (orange), d1T (green), and d1 (blue), fitting with β-1 (crosses), β-2 (triangles), and β-T (diamonds).

Fig. 7. Best-fit values of the second-order ${\mathcal{D}}_{\ell }^{A\times {\omega }_2^{\beta }}$, ${\mathcal{D}}_{\ell }^{{\omega }_1^{\beta }\times {\omega }_2^{\beta }}$ and ${\mathcal{D}}_{\ell }^{{\omega }_2^{\beta }\times {\omega }_2^{\beta }}$ moment parameters in d0 (orange), d1T (green), and d1 (blue) for β-2 (triangles).

Fig. 8. Best-fit values of the first-order ${\mathcal{D}}_{\ell }^{{\omega }_1^{\beta }\times {\omega }_1^T}$ and ${\mathcal{D}}_{\ell }^{{\omega }_1^T\times {\omega }_1^T}$ moment parameters in d0 (orange), d1T (green), and d1 (blue) for β-T (diamonds).

5.1.2. d1T

We now introduce, as a first layer of complexity, the spatial variations of the spectral index associated to a fixed temperature over the sky with the d1T simulation type. The dust temperature was fixed to T_d1T = 21.9 K while the spectral index β(n) was allowed to vary between lines of sight. The four different fitting schemes presented in Sect. 4.1 are fitted over the cross-spectra of our simulations as in Sect. 5.1.1.

The reduced χ²(ℓ) values for each fitting scheme can be found in Fig. 3. It can be seen that the MBB no longer provides a good fit for the dust SED, especially at low multipoles. Averaging effects of spatially varying SEDs are more important over large angular scales and thus SED distortions and moments are expected to be more significant at low multipoles. Indeed, the moments added to the fit in β-1 are enough to lower the reduced χ² such that it becomes compatible with 1 over almost all of the multipole range. The fitting schemes β-T and β-2, including more parameters than β-1, provide a fit of similar goodness, except in the multipole bin ℓ = 66.5 where they are closer to 1.

Figure 4 presents the best-fit values of β₀(ℓ) in the case of the MBB fit. For the sake of clarity, the values after iteration (see Sect. 4.1) for β-1, β-T, and β-2 are not shown, but they present comparable trends. We can see that the best-fit values of β₀(ℓ) for this d1T simulation type are no longer compatible with a constant. β₀(ℓ) fitted values show a significant increase at low ( < 100) multipoles, up to β₀ = (ℓ = 16.5) = 1.65. For ℓ > 100, β₀(ℓ) is close to a constant of value ∼ 1.53. This increase towards the low ℓ is correlated to the increase of the MBB χ² discussed in the previous paragraph. However, we note that in the lowest ℓ-bin, the β₀(ℓ) value is close to 1.53 and that the χ² of the MBB fit is close to unity.

The best-fit values of T₀(ℓ) are also presented in Fig. 4 in the case of the MBB fit. Here again, the values after iteration for the other fitting schemes are not presented, but are similar. The d1TT₀(ℓ) best-fit values oscillate around T_d1T = 21.9 K, without being strictly compatible with a constant value, as would be expected for this simulation type. This tends to indicate that the SED distortions due to the spectral index spatial variations are affecting the accuracy at which we can recover the correct angular dependence of the sky temperature.

The amplitude power spectrum is displayed in Fig. 5 for the MBB fitting scheme. The other fitting scheme results are not presented for clarity and would not be distinguishable from those of the MBB. The fitted ${\mathcal{D}}_{\ell }^{A\times A}$ is compatible with the one of the dust template map used to build the simulations.

All the parameters of the moment expansion with respect to β can be found in Figs. 6 and 7, and are now significantly detected, except for ${\mathcal{D}}_{\ell }^{{\omega }_2^{\beta }\times {\omega }_2^{\beta }}$. In Fig. 8, we can observe that the parameters of the moment expansion with respect to the temperature (only present in the β-T fit) remain undetected. The SED distortions due to the spatial variations of β are well detected, while no SED distortion linked to the temperature is seen, as expected for the d1T simulation type.

5.1.3. d1

We now discuss the d1 simulations, with the highest complexity in the polarized dust SED. In this more physically relevant simulation type, the dust emission is given by a MBB with variable index β(n) and temperature T(n) over the sky. We ran the four different fitting schemes on the d1 simulation type, as we did in Sects. 5.1.1 and 5.1.2.

The values of the reduced χ²(ℓ) are displayed in Fig. 3. For the MBB and β-1, the reduced χ² are not compatible with unity, especially at low multipole. This indicates that none of them are a good fit anymore for the spatially varying SED with β(n) and T(n). With β-2 and β-T, the χ²(ℓ) values become compatible with unity, except for the ℓ = 26.5 bin. We note that β-T provides a slightly better fit than β-2 in this bin.

Looking at the medians of the best-fit values of β₀(ℓ) for d1 in Fig. 4, we can see that the spectral index is changing with respect to ℓ, as discussed in Sect. 5.1.2, in a similar manner as for the d1T simulation type. The fitted temperature T₀(ℓ) values for d1 show an increasing trend from ∼ 17 to ∼ 20.5 K and from ℓ = 16.5 to ℓ ∼ 100. At higher multipoles, T₀(ℓ) is close to a constant temperature of 20.5 K. In d1, as for d1T, the angular scales at which we observe strong variations of β₀(ℓ) and T₀(ℓ) are the ones for which we observe a poor χ² for some fitting schemes. Also, as for d1T, the largest angular scale ℓ-bin, at ℓ = 6.5, shows β and T values close to the constant value at high ℓ, which are associated with χ² values closer to unity. The best-fit values of the amplitude ${\mathcal{D}}_{\ell }^{A\times A}$ are shown in Fig. 5. These are similar to those of the other simulation types.

The moment-expansion parameters fitted on d1 are shown in Figs. 6–8. For this simulation type, the moment parameters are all significantly detected with respect to both β and T. This was already the case with the Planck intensity simulations, produced in a similar way, as discussed in Mangilli et al. (2021). Their detections quantify the complexity of dust emission and SED distortions from the MBB present in the d1 simulation type, due to the spatial variations of β(n) and T(n).

5.2. CMB only

In order to calibrate the accuracy at which the r parameter can be constrained with the LiteBIRD simulated data sets presented in Sect. 3.3, we tested the simulation type with no dust component, $M_{\nu }^{\mathtt{c}}$, and with no tensor modes (r_sim = 0, only CMB lensing and noise). We fit the expression in Eq. (20) with ${\mathcal{D}}_{\ell }^{\mathrm{dust}}$ fixed to zero and A_lens fixed to one (i.e., r is the only parameter we fit in this case). Doing so over the N_sim simulations, we obtain $\widehat{r}=(0.7\pm 3.5)\times {10}^{-4}$. This sets the minimal value we can expect to retrieve for $\widehat{r}$ with our assumptions if the dust component is perfectly taken into account.

5.3. Dust and CMB

We now present our analysis of the simulations including dust, CMB (lensing), and noise (d0c, d1Tc and d1c) with no primordial tensor modes (r_sim = 0). As described above, we applied the four fitting schemes for the dust on the three simulation types, fitting $\widehat{r}$ and fixing A_lens to one (namely rMBB, rβ-1, rβ-T and rβ-2) simultaneously.

The best-fit values of β₀(ℓ), T₀(ℓ) and the moment expansion parameters ${\mathcal{D}}_{\ell }^{\mathcal{M}\times \mathcal{N}}$ derived with the simulation types d0c, d1Tc, and d1c are not discussed further when they are compatible with the ones obtained for the d0, d1T, and d1 simulation types and presented in Sect. 5.1.

5.3.1. d0c

For d0c, as for d0, we recover the input constant spectral index and temperature β_d0 and T_d0 at all angular scales for every fitting scheme. Furthermore, we do not detect any moment, when fitting rβ-1, rβ-T, and rβ-2. This simulation type therefore constitutes our null-test when $\widehat{r}$ and the dust parameters are fitted at the same time. The addition of the CMB lensing in the simulations and the addition of r to the fits thus does not lead to the detection of the moment parameters nor biases the recovery of the spectral index and the temperature.

The posterior distributions of the estimated tensor-to-scalar ratio $\widehat{r}$ are displayed in Fig. 9 and their mean and standard deviations are summarized in Table 4. We note that $\widehat{r}$ is compatible with the input value (r_sim = 0) for all the fitting schemes. For rMBB and rβ-1, the dispersion ${\sigma }_{\widehat{r}}$ is comparable with the CMB-only scenario discussed in Sect. 5.2. For rβ-T and rβ-2, the width of the distribution increases by a factor of ∼ 2 and ∼ 4, respectively.

Fig. 9. Top panel: posterior on $\widehat{r}$ in the d0c simulation type for the different fitting schemes: rMBB (blue, dotted line), rβ-1 (red, dashed line), rβ-T (green, solid line), and rβ-2 (yellow, dash-dotted line). The vertical black dashed line marks the value of rsim = 0. Central panel: same, but in the case of the d1Tc simulation type. Bottom panel: same, but in the case of the d1c simulation type.

5.3.2. d1Tc

The posterior distribution of $\widehat{r}$ in the case of the d1Tc simulation type is displayed in Fig. 9 for the four fitting schemes and the mean value and standard deviation of these distributions are summarized in Table 4. We can see that in the case of rMBB, we fit $\widehat{r}\pm {\sigma }_{\widehat{r}}=(99.7\pm 6.2)\times {10}^{-4}$. In that case, the input tensor-to-scalar ratio r_sim = 0 is not recovered and we obtain a bias on the central value of $\widehat{r}$ of $\sim 16{\sigma }_{\widehat{r}}$. As discussed in Sect. 2.1.2, this is expected because we know that the MBB is not a good dust model for a SED with spatially varying spectral index, as we also verify in Sect. 5.1.2 looking at the χ² values.

Using the rβ-1 fitting scheme allows us to recover $\widehat{r}=(-8.0\pm 6.4)\times {10}^{-4}$, where r_sim is recovered within $\sim 2{\sigma }_{\widehat{r}}$, while rβ-2 and rβ-T recover the input value within $1{\sigma }_{\widehat{r}}$ (with $\widehat{r}=(-3.6\pm 13.0)\times {10}^{-4}$ and $\widehat{r}=(0.7\pm 20.9)\times {10}^{-4}$, respectively). As in Sect. 5.3.1, the deviation remains similar between rMBB and rβ-1 and increases by a factor of ∼ 2 and ∼ 4 from rβ-1 to rβ-T and rβ-2, respectively.

5.3.3. d1c

In the case of the d1c simulation type, as in d0c and d1Tc, we fit $\widehat{r}$ in addition to the dust-related parameters. In that case, dust moment parameters are recovered as for d1 (see Sect. 5.1.3), except for the rβ-2 fitting scheme.

Figure A.2 compares the moment parameters between β-2 on the d1c simulations type, fitting only the dust-related parameters and rβ-2 on d1c when jointly fitting the dust parameters and $\widehat{r}$. We observe that ${\mathcal{D}}_{\ell }^{{\omega }_2^{\beta }\times {\omega }_2^{\beta }}$ is not consistently recovered when fitting $\widehat{r}$ in addition to the dust parameters.

A similar comparison can be found in Fig. A.3 for the moment parameters between β-T and rβ-T on the d1c simulation type. Using this fitting scheme, we can see that all the moments are correctly recovered when adding $\widehat{r}$ to the fit.

The $\widehat{r}$ posterior distributions in the case of d1c are displayed in Fig. 9 and summarized in Table 4. As discussed in Sect. 2.1.2 and observed in Sect. 5.3.3, the rMBB fit is highly biased, with $\widehat{r}=(125.1\pm 5.9)\times {10}^{-4}$ (by more than 21 ${\sigma }_{\widehat{r}}$). When fitting the rβ-1, this bias is significantly reduced ($\widehat{r}=(32.9\pm 6.5)\times {10}^{-4}$, $5{\sigma }_{\widehat{r}}$ away from r_sim = 0), illustrating the ability of the first-moment parameters to correctly capture part of the SED complexity. However, performing the expansion in both β and T with rβ-T allows us to recover r_sim without bias ($\widehat{r}=(-3.3\pm 11.7)\times {10}^{-4}$), highlighting the need for the description of the SED complexity in terms of dust temperature for this simulated data set where both β and T vary spatially. On the other hand, for rβ-2, a negative tension ($1.9{\sigma }_{\widehat{r}}$) can be observed: $\widehat{r}=(-37.4\pm 19.4)\times {10}^{-4}$. This tension is discussed in Sect. 6.5.

For d1c, the $\widehat{r}$ distribution widths roughly meet the foreground cleaning requirements of LiteBIRD presented in Sect. 3.1 for rMBB and rβ-1 but are higher for rβ-T and rβ-2. We also note that, with the same number of free parameters, all the standard deviations ${\sigma }_{\widehat{r}}$ slightly increase compared to the d0c simulation type. This is expected due to the increasing dust complexity.

6. Discussion

6.1. Lessons learnt

In Sect. 5, we apply the fitting pipeline introduced in Sect. 4 on LiteBIRD simulated data sets on f_sky = 0.7 and for r_sim = 0, including the various dust simulation types defined in Sect. 3.2. We fitted the estimated B-mode power-spectra with the four different fitting schemes summarized in Table 3. Our main results can be summarized as follows:

• The MBB fitting scheme provides a good fit for the dust component in the d0 and d0c simulation types. However, when the spectral index changes with the angular scale, such as in the d1T, d1Tc, d1, and d1c simulations, this approach no longer provides a good fit because of the complexity of the dust SED. As a consequence, in the rMBB case, r_sim cannot be recovered without a significant bias.

• The β-1 fitting scheme allows us to perform a good fit for the dust complexity using the d0 and d1T simulations but not for d1, while the rβ-1 fitting scheme yields estimates of $\widehat{r}$ close to r_sim within $1{\sigma }_{\widehat{r}}$ for d0c, and within $2{\sigma }_{\widehat{r}}$ for d1Tc, but presenting a bias of $\sim 6{\sigma }_{\widehat{r}}$ for d1c.

• The β-T fitting scheme provides a good fit for every dust model, while using the rβ-T fitting scheme allows us to recover $\widehat{r}$ values consistent with r_sim within $1{\sigma }_{\widehat{r}}$ for all the simulation types, but is associated with an increase of ${\sigma }_{\widehat{r}}$ by a factor ∼ 2 compared to the rβ-1 case.

• The β-2 fitting scheme also provides a good fit for each dust model, and the rβ-2 fitting scheme leads to values of $\widehat{r}$ compatible with r_sim within $1{\sigma }_{\widehat{r}}$ for all the simulation types but d1c. In this last case, there is a negative tension of $\sim 2{\sigma }_{\widehat{r}}$. For all the simulation types, there is an increase of ${\sigma }_{\widehat{r}}$ by a factor of ∼ 4 compared to the rβ-1 case.

The present analysis shows that the temperature could be a critical parameter for the moment expansion in the context of LiteBIRD.

Indeed, for simulations including a dust component with a spectral index and a temperature that both vary spatially, as in d1, the only fitting scheme allowing us to recover r_sim within $1{\sigma }_{\widehat{r}}$ is rβ-T, the expansion to first order in both β and T. This shows that expanding in β only, without treating T, is not satisfactory when looking at such large fractions of the sky. Indeed, when applying the β-2 fitting scheme, the ${\mathcal{D}}_{\ell }^{{\omega }_2^{\beta }\times {\omega }_2^{\beta }}$ parameter remains undetected for the d1T simulation type (Sect. 5.1.2), while it is significantly detected using the d1 simulation type (Sect. 5.1.3). Nevertheless, d1T and d1share the same template of β(n) (Sect. 3.2) and they only differ by the sky temperature (constant for d1T and varying for d1). This suggests that the observed ${\mathcal{D}}_{\ell }^{{\omega }_2^{\beta }\times {\omega }_2^{\beta }}$ with the d1 simulations originates from the temperature variations and not those in the spectral index. This observation shows that it is less convenient to use the β-2 fitting scheme than the β-T one in order to correctly recover the moment-expansion parameters and $\widehat{r}$ when temperature varies spatially.

Moreover, we saw that ${\sigma }_{\widehat{r}}$ is lower when using the fitting scheme rβ-T instead of rβ-2 for every simulation type, even if both have the same number of free parameters. This second observation additionally encourages an approach where the SED is expanded with respect to both β and T. Nevertheless, the uncertainty on $\widehat{r}$ we obtain in this case (${\sigma }_{\widehat{r}}=1.17\times {10}^{-3}$) is larger than the LiteBIRD requirements.

6.2. Increasing the accuracy on the tensor-to-scalar ratio

In Sect. 5.1.3 and Fig. 3, we see that the MBB and β-1 fitting schemes do not provide good fits for the d1 dust simulations, especially at low multipoles (ℓ ≲ 100). Conjointly, in Fig. 8, we can see that the β-T moment parameters are significantly detected for ℓ ≲ 100 and compatible with zero above that threshold, suggesting that their corrections to the SED are predominantly required at large angular scales.

This implies that we can improve the pipeline presented in Sect. 4 to keep only the required parameters in order to recover $\widehat{r}$ compatible with r_sim with a minimal ${\sigma }_{\widehat{r}}$. It can be achieved by applying the rβ-1 fitting sheme over the whole multipole range, while restricting the rβ-T-specific (${\mathcal{D}}_{\ell }^{{\omega }_1^{\beta }\times {\omega }_1^T}$ and ${\mathcal{D}}_{\ell }^{{\omega }_1^T\times {\omega }_1^T}$) moment-expansion parameters fit to the low multipoles range. We note that in order to correct the bias, it is still necessary to keep the rβ-1 moment parameters even at high multipoles, because the MBB does not provide a good fit even for ℓ ∈ [100,200], as we can see in Fig. 3. We define ℓ_cut as the multipole bin under which we keep all the rβ-T moment parameters and above which we use the rβ-1 scheme.

The best-fit values and standard deviations of $\widehat{r}$ for different values of ℓ_cut are displayed in Table 5. We can see that a trade-off has to be found: the smaller the ℓ_cut, the bigger the shift from r_sim, and the bigger the ℓ_cut, the higher the value of ${\sigma }_{\widehat{r}}$. The trade-off point seems to be found for ℓ_cut, allowing us to recover $\widehat{r}$ without tension, with ${\sigma }_{\widehat{r}}=8.8\times {10}^{-4}$. The error on r is thus reduced by more than ∼ 30% with respect to the nonoptimized fit and meets the LiteBIRD requirements.

6.3. Tests with smaller sky fractions

In all the results presented in Sect. 5, we were considering a sky fraction of f_sky = 0.7. This sky mask keeps a considerable fraction of the brightest Galactic dust emission. To quantify the impact of the sky fraction on our analysis, we ran the pipeline as in Sect. 5.3.3 with the different masks introduced in Sect. 3.4.1 (f_sky = 0.5 and f_sky = 0.6). This was done with the d1c simulation type.

The posteriors on $\widehat{r}$ for the different fitting schemes are displayed in Fig. 10 and Table 6. We can see that, while the rMBB fiting scheme always leads to biased estimates, the rβ-1 case allows us to recover $\widehat{r}$ at $1.25{\sigma }_{\widehat{r}}$ for f_sky = 0.5 and within $1{\sigma }_{\widehat{r}}$ for f_sky = 0.6. In the two situations, the results using the rβ-T and β-2 fitting schemes are both unbiased with estimates of $\widehat{r}$ compatible with r_sim within $1{\sigma }_{\widehat{r}}$. The ${\sigma }_{\widehat{r}}$ hierarchy between the rMBB, rβ-1, rβ-T, and rβ-2 fitting schemes is the same as for f_sky = 0.7 (see Sect. 5.3.3). Nevertheless, we observe that ${\sigma }_{\widehat{r}}$ increases as the sky fraction decreases, as does the statistical error (cosmic variance of the lensing and noise). The bias, on the other hand, decreases for all the fitting schemes with the sky fraction, which is expected because less dust emission contributes to the angular power spectra. The negative tension observed on the $\widehat{r}$ posterior in Sect. 5.3.3 for the rβ-2 case is not present when using smaller sky fractions. In Fig. A.5, the rβ-2 moment parameters are displayed. We can see that they are not significantly detected for the f_sky = 0.5 and 0.6, unlike for f_sky = 0.7. As we have seen that some of the moments in the β-2 fitting scheme failed to model SED distortions coming from temperature, we can suppose that, in our simulations, the temperature variations play a less significant role in the dust SED on the f_sky = 0.5 and 0.6 masks than they play in the f_sky = 0.7 one. As a consequence, they have a smaller impact on r when not properly taken into account.

Fig. 10. Top panel: posterior on $\widehat{r}$ in the d1c simulation type on fsky = 0.5 for the different fitting schemes: rMBB (blue, dotted line), rβ-1 (red, dashed line), rβ-T (green, solid line), and rβ-2 (yellow, dash-dotted line). The vertical black dashed line marks the value of rsim = 0. Bottom panel: same, in the case of the d1c simulation type on fsky = 0.6.

6.4. Tests with nonzero input tensor modes

We show in Sect. 5.3.3 that the rβ-T fitting scheme allows us to retrieve $\widehat{r}$ compatible with zero when r_sim = 0. We now want to assess the potential leakage of $\widehat{r}$ in the moment expansion parameters if r_sim ≠ 0. In this case, primordial tensor signals would be incorrectly interpreted as dust complexity. We run the pipeline as described in Sect. 5.3.3 with r_sim = 0.01, in the d1c simulation type. This value of r_sim = 0.01 is larger than the value targeted by LiteBIRD, but given the order of magnitude of the error on $\widehat{r}$ observed in the previous sections, a potential leakage could be left unnoticed using a smaller r_sim.

Looking at the final posterior on $\widehat{r}$ (Fig. 11 and Table 6), we can see that the results are comparable with the r_sim = 0 case, but centered on the new input value r_sim = 0.01. The rMBB fitting scheme gives a highly biased posterior of $\widehat{r}=(2.048\pm 0.077)\times {10}^{-2}$; the bias is reduced but still significant when using the rβ-1 scheme ($\widehat{r}=129.0\pm 8.3\times {10}^{-4}$); in the β-T case we get an estimate of $\widehat{r}=94.6\pm 15.1\times {10}^{-4}$ compatible with the input value of r_sim = 100 × 10⁻⁴; and finally, the β-2 fitting scheme leads to a negative $2{\sigma }_{\widehat{r}}$ tension ($\widehat{r}=62.5\pm 25.0\times {10}^{-4}$). This demonstrates the robustness of our method and its potential application to component separation. We note that the negative bias at second order is still present in the r_sim = 0.01 case, illustrating that setting a positive prior on $\widehat{r}$ would not have been a satisfying solution when r_sim = 0.

Fig. 11. Posterior on $\widehat{r}$ in the d1c simulation type with rsim = 0.01 and fsky = 0.7 for the different fitting schemes: rMBB (blue, dotted line), rβ-1 (red, dashed line), rβ-T (green, solid line), and rβ-2 (yellow, dash-dotted line). The vertical black dashed line marks the value of rsim.

6.5. Exploring the correlations between the parameters

We now examine the substantial increase in the dispersion on the $\widehat{r}$ posteriors between the rβ-1 fitting scheme on the one hand and the rβ-T and rβ-2 ones on the other. Indeed, in Sect. 5.3.3, we show that ${\sigma }_{\widehat{r}}$ is about two times greater when using the rβ-T scheme than the rβ-1 one, and about four times larger in the case of rβ-2, while the rβ-T and rβ-2 schemes share the same number of free parameters. Some other points to clarify are the shift on $\widehat{r}$ appearing for rβ-2 in the d1c scenario, discussed in Sect. 5.3.3, and the inability to correctly recover ${\mathcal{D}}_{\ell }^{{\omega }_2^{\beta }\times {\omega }_2^{\beta }}$ when $\widehat{r}$ is added to the fit illustrated in Fig. A.2.

The 2D-SED shapes of the parameters ${\mathcal{D}}_{\ell }^{\mathcal{N}\times \mathcal{M}}({\nu }_i\times {\nu }_j)$ in the (ν_i,ν_j) space⁷ are displayed in Fig. A.4. We used the nine frequencies of LiteBIRD presented in Sect. 3.2.3 and fixed β₀ = 1.54 and T₀ = 20 K. We also introduce the CMB 2D-SED shape with the black body function:

$$\begin{array}{c}\hfill B_{\mathrm{CMB}}({\nu }_i\times {\nu }_j)=\frac{B_{{\nu }_i}(T_{\mathrm{CMB}})B_{{\nu }_j}(T_{\mathrm{CMB}})}{B_{{\nu }_0}{(T_{\mathrm{CMB}})}^2},\end{array}$$(22)

where T_CMB = 2.726 K.

The 2D correlation coefficients between these 2D-SED shapes are displayed in Fig. 12. We present the correlations between the shapes of the parameters in the case of the rβ-T and rβ-2 fitting schemes. We can see that all the moment parameters in ${\omega }_2^{\beta }$ are strongly correlated with the CMB SED signal, while the ones in ${\omega }_1^T$ are not.

Fig. 12. Correlation matrices of the 2D-SED shapes of the CMB (BCMB(νi × νj) and dust moments ${\mathcal{D}}_{\ell }^{\mathcal{N}\times \mathcal{M}}({\nu }_i\times {\nu }_j)$ in the (νi,νj) space). Each element represents the Pearson correlation coefficient between any 2 of these 2D-SED shapes. The correlation matrices are displayed in the case of the β-2 fitting scheme (top panel) and the β-T one (bottom panel).

We showed that, when fitting β-2 on d1c, the SED distortions due to spatial variations of T are incorrectly detected by the second-order moment parameters with respect to the spectral index β. Due to the correlations highlighted above, those spurious moment parameters could then leak into $\widehat{r}$ when adding it to the fit in rβ-2. This explains both the negative shift on the $\widehat{r}$ posterior using β-2 in the d1c simulation type with f_sky = 0.7 presented in Sects. 5.3.3 and 6.4, and the inability to correctly recover the ${\omega }_2^{\beta }\times {\omega }_2^{\beta }$ dust moment parameter presented in Fig. A.2. In addition, it gives a natural reason for the surge of ${\sigma }_{\widehat{r}}$ when the second-order moments in β are added to the fit.

On the other hand, the moment parameters in ${\omega }_1^T$ are strongly correlated with the moments in ${\omega }_1^{\beta }$. This behavior is expected due to the strong correlation between β and T (see e.g., Juvela & Ysard 2012). However those moment parameters are less correlated with the CMB signal than the second-order parameters of β-2. This points out that the factor of ∼ 2 on ${\sigma }_{\widehat{r}}$ between β-T and β-2 is due to this correlation of the 2D-SED shapes. As the parameters in ${\omega }_1^T$ are highly correlated with one another, we expect them to be highly redundant in the fit. However, repeating the process described in Sect. 5.3.3 using only ${\mathcal{D}}_{\ell }^{A\times {\omega }_1^T}$ for β-T –which is equivalent to applying the β-1 fitting scheme with an iterative correction to the temperature T₀(ℓ)– gives a $\widehat{r}$ posterior similar to the one obtained for β-1 alone. Taking the other ${\omega }_1^T$ terms into account appears to be necessary in order to recover an unbiased distribution of $\widehat{r}$.

6.6. Adding synchrotron to the simulations

Thermal dust is not the only source of polarized foreground that must be considered for CMB studies. Although subdominant at high frequencies (≥ 100 GHz), the synchrotron emission due to accelerated electrons in the interstellar medium is still expected to represent a significant fraction of the total polarized signal.

In order to take one more step towards realistic forecasts for the LiteBIRD instrument, we add a synchrotron contribution to the d1c simulations presented in Sect. 3.3 using the s1 template included in the PYSM. In this scenario, the synchrotron SED for each line of sight is given by a power law of the form (in antenna temperature units)

$$\begin{array}{c}\hfill S_{\nu }^{\mathtt{s}}\mathtt{1}=A_{\mathtt{s}}\mathtt{1}(\mathit{n}){\left(\frac{\nu }{{\nu }^{\mathtt{s}}{\mathtt{1}}_0}\right)}^{{\beta }_{\mathtt{s}}\mathtt{1}(\mathit{n})},\end{array}$$(23)

where the amplitude A_s1(n) and the spectral index β_s1(n) maps are derived from the combination of the WMAP mission 23 GHz map Bennett et al. (2013) and Haslam 408 GHz map Haslam et al. (1981). ${\nu }^{\mathtt{s}}{\mathtt{1}}_0$ is defined as 23 GHz. The simulations containing synchrotron are referred to as d1s1c below.

If not treated in the fit, the presence of synchrotron is expected to induce a bias on the $\widehat{r}$ posterior distribution. Regarding the dust MBB discussed in Sect. 2.1.2, the synchrotron SED is expected to have distortions. However, as the synchrotron polarized emission is significantly lower than that of dust, in the frequency range considered in the present work, we expect the distortions to be small compared to the ones induced by dust and we leave their modeling to a further study.

In order to minimize the number of free parameters used for fitting the synchrotron emission, we model its power spectrum as a power law of the multipole ℓ (Krachmalnicoff et al. 2018). Therefore, combining with the synchrotron SED in Eq. (23), the synchrotron component of the cross-angular power spectra reads

$$\begin{array}{c}\hfill {\mathcal{D}}_{\ell }^{\mathrm{sync}}({\nu }_i\times {\nu }_j)=A_{\mathrm{s}}{\left(\frac{{\nu }_i{\nu }_j}{{\nu }_0}\right)}^{{\beta }_{\mathrm{s}}}{\ell }^{{\alpha }_{\mathrm{s}}},\end{array}$$(24)

where the amplitude coefficient A_s is treated as a free parameter while we fix β_s = −3 (median value of the s1 β_s map on our f_sky = 0.7 mask) and α_s = −1 (Krachmalnicoff et al. 2018).

When fitting the d1s1c simulations, we either use the rβ-T fitting scheme, neglecting the synchrotron component, or we add the synchrotron component in Eq. (24) to the model in Eq. (20). We refer to this latter case as the srβ-T fitting scheme. In Fig. 13, the $\widehat{r}$ posteriors derived from the d1s1c simulations are displayed with r_sim = 0 and f_sky = 0.7.

Fig. 13. Posterior on $\widehat{r}$ in the d1s1c simulation type with rsim = 0 and fsky = 0.7 for the different fitting schemes: rβ-T (green, solid line) and srβ-T (orange, dash-dotted line). The vertical black dashed line marks the value rsim = 0.

Using the rβ-T fitting scheme, we find $\widehat{r}=(143.1\pm 13.5)\times {10}^{-4}$. As expected, even at high frequencies, modeling the synchrotron component is critical and cannot be neglected in order to recover an unbiased value of $\widehat{r}$. On the other hand, using srβ-T fitting scheme, we recover $\widehat{r}=(-5.4\pm 13.2)\times {10}^{-4}$. This result is comparable with the one obtained for the d1c simulations in Sect. 5.3.3, with a minor increase in ${\sigma }_{\widehat{r}}$. We can therefore conclude that a model as simple as that of Eq. (24) is sufficient to take into account the s1 component at ν > 100 GHz and the corresponding SED distortions can be neglected in order to recover an unbiased value of $\widehat{r}$. In principle, as we know that the dust-synchrotron spatial correlation is significant at large scales (Planck Collaboration XI 2020), Eq. (20) should include a dust-synchrotron term (see e.g., Hensley et al. 2021). In our study, where we consider cross-spectra from 100 to 402 GHz, this dust-synchrotron term is subdominant, but it could be significant when considering cross-spectra between LiteBIRD’s extreme frequency bands (e.g., the 40 × 402 cross-spectrum). The moment expansion might be more complicated as well in this case, as we could expect some correlation between the dust and synchrotron moment-terms.

This result shows that a full polarized foreground content can be treated at high frequencies when using a power law SED for the synchrotron coupled with the moment expansion of the MBB up to first order in both β and T for the dust SED. A full study remains to be done in that direction using all the frequency bands of the LiteBIRD instrument. Eventually, Eq. (24) will also have to be expanded in moments with respect to its parameters. Doing so, one can expect to recover an unbiased value of $\widehat{r}$ associated with a decrease in ${\sigma }_{\widehat{r}}$ down to a value compatible with the full success criterion of the mission.

6.7. Limitations of this work and caveats

As discussed in Sect. 2.2.1, we neglected polarization effects through this work by treating the BB signal as an intensity signal. This is not problematic in the present work, because no variations along the lines of sight were present in the simulations. However, this point has to be addressed using complex simulations or real sky data.

The choice of reference frequency ν₀ used for the normalization of the MBB in Eq. (1), which is not discussed in this study, can potentially have a significant impact on the moment expansion and, in turn, on the measurement of $\widehat{r}$. Indeed, ν₀ is the pivot frequency of the moment expansion (moments are equal to zero at ν₀) and will determine the shape of the SED distortion around it. A poor choice for this reference frequency can have disastrous consequences for the moment fit: for example, if it is chosen far away from the observed bands, all the moments will become degenerated. In our case, the reference frequency (353 GHz) is within the observed frequency range (100 − 402 GHz), but we have not tried to optimize its position. In addition, the ν₀ pivot of our moment expansion coincides with the one used to extrapolate the dust template map in the PYSM and we have not quantified how much this impacts our results.

Finally, as pointed out several times in this work, the quantitative results depend strongly on the sky model of our simulations. Moreover, we lack dedicated sky models where we can control the complexity of the dust SED, either by directly including moments or by averaging the emission from the 3D structure of the Galaxy. However, both methods are beyond the scope of the present work.

7. Conclusion

Being able to precisely characterize the complexity of the Galactic thermal dust polarized SED has become critical for the measurement of the faint primordial B-mode signal of the CMB, especially at the sensitivity targeted by future CMB experiments such as the LiteBIRD satellite mission.

In this work, we applied the moment expansion formalism to the dust emission SED as a component-separation tool to recover the tensor-to-scalar ratio parameter r in LiteBIRD-simulated data. This formalism, proposed by Chluba et al. (2017) and implemented in harmonic space by Mangilli et al. (2021), allows us to deal with the spectral complexity of the Galactic dust signal by modeling its deviations from the canonical MBB model at the cross-angular power spectrum level. In the case of the data-driven realistic dust emission model –we explore (PYSM d1) here–, suitably taking into account the dust SED distortions prevents the spurious detection of the primordial B-mode signal.

We show that the dust spectral index β and dust temperature T spatial variations significantly distort the dust cross-power spectrum SED. The MBB is not a good model to describe the data in that case and the estimation of r is dramatically affected. In the case where no primordial signal is included in the simulated data sets, not taking into account the dust SED complexity leads to a highly significant spurious detection of r with LiteBIRD (from $\widehat{r}\simeq 5\times {10}^{-3}$ to 1.25 × 10⁻², with a 8.4 to 21.2σ significance, from 50 to 70% of the sky, respectively).

To overcome this obstacle, we applied the moment expansion formalism in order to model these SED distortions. We demonstrate that, at LiteBIRD sensitivity, the previously studied moment expansion with respect to the dust spectral index β (Mangilli et al. 2021; Azzoni et al. 2021) does not give satisfactory results. Indeed, expanding in β to first order (following the angular power spectrum definition of the order) leads to a significant bias on 70% of the sky ($\widehat{r}=(3.29\pm 0.65)\times {10}^{-3}$ when r_sim = 0 and $\widehat{r}=(1.29\pm 0.08)\times {10}^{-2}$, when r_sim = 10⁻²). At second order in β, we observe a ∼ 2 σ negative tension ($\widehat{r}=(-3.7\pm 1.9)\times {10}^{-3}$ when r_sim = 0 and $\widehat{r}=(6.25\pm 2.50)\times {10}^{-3}$, when r_sim = 10⁻²).

We introduce for the first time in this work the expansion of the dust angular cross-power spectra with respect to both β and T. We show that by using this expansion up to first order, we correctly model the dust SED distortions due to spatial variations of both β and T at the map level. This allows us to recover r parameter without bias, with $\widehat{r}=(-3.3\pm 11.7)\times {10}^{-4}$ if r_sim = 0 and $\widehat{r}=(0.95\pm 0.15)\times {10}^{-2}$ if r_sim = 10⁻². Thus, despite the known degeneracy between the dust spectral index and its temperature in the Rayleigh-Jeans domain, it is important to correctly model the latter in order to accurately retrieve the tensor-to-scalar ratio r at the unprecedented precision reached by experiments such as LiteBIRD.

Adding parameters to tackle the dust SED complexity means an increase in the error budget. Given the LiteBIRD bands and sensitivities we consider in this work (frequency bands above 100 GHz), the ideal sensitivity on r without delensing is ${\sigma }_{\widehat{r}}=3.4\times {10}^{-4}$. In the ideal case, where the dust β and T are constant over the sky (PYSM d0), separating the CMB from dust leads to ${\sigma }_{\widehat{r}}=3.9\times {10}^{-4}$ on 70% of the sky. Adding the expansion to first order in β does not significantly increase the error (${\sigma }_{\widehat{r}}=4.5\times {10}^{-4}$), but expanding to first order in both β and T multiplies it by a factor of ∼ 2 (${\sigma }_{\widehat{r}}=9.5\times {10}^{-4}$) and to second order in β by a factor of ∼ 4 (${\sigma }_{\widehat{r}}=16.4\times {10}^{-3}$). We show that the surge of ${\sigma }_{\widehat{r}}$ between the two latter cases, sharing the same number of free parameters, is due to strong correlations between the SED of the second-order moments in β and the CMB. This is an important point, as it could lead to some intrinsic limitation for component-separation algorithms based exclusively on the modeling of the SED. Furthermore, when dealing with real data, if the dust SED is complex enough to have significant second-order distortions with respect to β, CMB experiments might reach a dilemma: either include the second order in the modeling at the cost of losing sensitivity on r, or neglect it at the cost of a potential spurious detection. Coupling the SED-based separation with methods exploiting the diversity of spatial distribution between components (e.g., Regaldo-Saint Blancard et al. 2021) seems a natural way to overcome this issue.

Nevertheless, moment expansion at the cross-angular power spectrum level provides a powerful and agnostic tool, allowing us to analytically recover the actual dust complexity without making any further assumptions. We additionally show that this method is robust, in the sense that it can effectively distinguish the primordial tensor signal from dust when r_sim ≠ 0, as in the case of LiteBIRD simulations. The dust moments in β and T at first order are needed in order to retrieve a reliable measure of r; they are significantly detected for ℓ ≲ 100. We can therefore define a cut in ℓ above which we do not fit for the whole complexity of the dust (we fit only the expansion up to first order in β and not in β and T). Doing so, we can reduce the error on $\widehat{r}$ while keeping the bias negligible ($\widehat{r}=(-0.9\pm 8.8)\times {10}^{-4}$). We could imagine other ways to reduce the number of free parameters in our model (e.g., assuming a power-law of ℓ behavior for the moments, as in Azzoni et al. 2021) and hence reduce the error on r. However, this optimization really depends on the simulated sky complexity and has not been comprehensively explored in the present work.

The PYSM d1 sky simulations, being data-driven, are widely used by the CMB community as they contain some of the real sky complexity. Nevertheless, at high-Galactic latitudes, the dust spectral index and temperature templates from Planck are dominated by systematic errors (uncertainty on the assumed zero-level of the Planck intensity maps, residual cosmic infrared background (CIB), anisotropies, instrumental noise, etc.). Therefore, some of the complexity we observe far from the Galactic plane in this sky model is not real. On the other hand, the modeled SED of the dust is exactly a MBB in each pixel, and line-of-sight averages or more complex dust models are ignored. As a consequence, our method and CMB B-mode component-separation algorithms in general need to be confronted with more complex models in order to really assess their performances in a quantitative manner.

Finally, although we demonstrate that the synchrotron component can be tackled at frequencies above 100 GHz with a minimal model under our assumptions, a study over the full LiteBIRD frequency bands, including synchrotron and the potential moment expansion of its SED, will be considered as a natural next step for a further application.